(* Helper package for minimal surface notebooks http://www.indiana.edu/~minimal Matthias Weber Public Version 0.1 13/24/2004 This package must be distributed without charge. Copyright is by the author. Changes must be notified in this header. This package is compatibel with MMA 4.0 upwards. *) BeginPackage[ "Own`Mesh`", { "Utilities`FilterOptions`"} ]; StereographicProjection::usage = "StereographicProjection[z] returns the image of z under stereographic \ \ \ projection."; Domain::usage = "Domain is a 2-dimensional graphics format for triangle meshes."; Mesh3D::usage = "Mesh3D is a 3-dimensional graphics format for rectangle meshes."; Vertices::usage = "DomainVertices[list] is a list of coordinates."; RectangularDomain::usage = "RectangularDomain[x,y] creates a Domain object representing a \ \ rectangular \ \ grid with vertical and horizontal lines given by the lists \ x \ and y."; PolarDomain::usage = "PolarDomain[x,y] creates a Domain object representing a polar \ \ grid \ \ with radii and angles given by the lists x and y."; NRange::usage="NRange"; MeshApply::usage= "MeshApply[f,mesh] evaluates the complex function f on a 2-dimensional mesh and returns the image mesh."; TestMesh::usage="TestMesh[mesh] returns basic information about a mesh."; MeshBoundingBox::usage="MeshBoundingBox[mesh] returns the bounding box \ \ extremal vertices of the mesh."; DomainToGraphics::usage = "DomainToGraphics[mesh] converts a Domain object in a graphics object."; MeshPlot3D::usage = "MeshPlot3D[{fn},m] evaluates the function fn on the mesh vertices of m \ \ and \ \ returns a Mesh3D object."; MeshJoin::usage = "MeshJoin[mesh1,mesh2] joins two meshes,thereby mergin the points \ \ in \ \ merge into pairs."; Mesh3DToGraphics3D::usage = "Mesh3DToGraphics3D[m] converts a Mesh3D object into a Graphics3D \ \ \ \ object."; MeshRotate::usage = "MeshRotate[mesh,edge] rotates a mesh about an edge"; MeshRotateZ::usage="MeshRotateZ[mesh,phi,trans] applies scre motion by angle \ \ phi and tanslation trans along z-axis"; StraightLine::usage = "StraightLine[v1,v2] represents a straight line through \ \ v1 and v2."; Plane::usage="Plane[n-d] represents a plane with normal vector n at distance \ \ d from the origin."; MeshReflect::usage = "MeshReflect[mesh,plane] reflects a mesh at a plane"; MeshTranslate::usage="MeshTranslate[mesh,v] translates mesh by vector."; MeshScale::usage="MeshScale[mesh,factor] scales the mesh by factor."; MeshFlip::usage="MeshFlip[mesh] changes the orientation of the normals."; AlmostEqual::usage="Returns true if the two arguments are almost equal"; Bjoerling::usage="Bjoerling[c,n,z] solves the Bjoerling problem for c and n \ wrt to z"; MeshNormalize::usage=""; Begin["`Private`"]; eTolerance = 10.^-7; norm[v_]:=v.v; tonorm[v_]:=v/Sqrt[v.v]; Bjoerling[c_,n_,z_]:= Module[{w},c[z]-I Integrate[Cross[n[w],D[c[w],w]],{w,0,z}]] StereographicProjection[z_] := {(2*Re[z])/(1 + Abs[z]^2),(2*Im[z])/(1 + \ \ Abs[z]^2), 1 - 2/(1 + Abs[z]^2)}; Vertices[dm_/;(Head[dm]==Domain)] := Flatten[List@@dm,1]; Vertices[msh_/;(Head[msh]==Mesh3D)] :=Flatten[Transpose[List@@msh][[1]],2]; (* mesh generating functions *) RectangularDomain[x_, y_] := Domain@@ Outer[List,x,y]; PolarDomain[x_, y_] := Domain@@ Outer[#1{ Cos[#2],Sin[#2]}&,x,y]; reim[z_] := {Re[z], Im[z]}; ComplexForm[fn_][{x_, y_}] := reim[fn[x + I y]]; MeshApply[fn_, v_] := Domain@@Map[ComplexForm[fn],List@@v,{2}]; AlmostEqual[x_,y_]:=Abs[x-y]<0.00001; (* no normals: *) nrml[x_,y_,z_]:=tonorm[Cross[x-z,y-z]]; eps=0.001; MeshPlot3D[{fn_}, {x_, y_}, v_] := Module[{f1,n1}, f1 = Function[{x, y}, fn]; n1=Function[{x,y},nrml[f1[x+eps,y],f1[x,y+eps],f1[x,y]]]; Mesh3D[ {Apply[f1, List@@v, {2}],Apply[n1, List@@v, {2}]} ] ]; (* yes normals: *) MeshPlot3D[{fn_, nn_}, {x_, y_}, v_ ]:= Module[{f1, n1}, f1 = Function[{x, y}, fn]; n1 = Function[{x, y}, nn]; Mesh3D[ {Apply[f1, List@@v, {2}],Apply[n1, List@@v, {2}]} ] ]; SetAttributes[MeshPlot3D, HoldFirst]; NRange[a_,b_,n_]:=Range[N[a],N[b],N[b-a]/N[n-1]]; (* Mesh manipulation *) (*MeshJoin[m1_, m2_] := Mesh3D@@ Join[List@@m1,List@@m2];*) MeshJoin[m1__] := Mesh3D@@ Join[m1]; (* rotate p about line through v1 and v2 *) rotate[p_, {v1_, v2_}] := -p + 2v1 + 2(p - v1).(v1 - v2)/((v1 - v2).(v1 - v2)) (v1 - v2); (* differential of former map, use for normals *) nrotate[p_, {v1_, v2_}] := -p + 2p.(v1 - v2)/((v1 - v2).(v1 - v2)) (v1 - v2); meshrotate0[{v_,n_}, {v1_, v2_}] := {Map[rotate[#, {v1, v2}] &, v,{2}], -Map[nrotate[#, {v1, v2}] &, n,{2}] }; MeshRotate[l_, StraightLine[v1_, v2_]] := Mesh3D@@Map[meshrotate0[#,{v1,v2}]&,List@@l]; zrotate[{v_, n_}, phi_, trans_] := {Map[zrotate0[#, phi, trans] &, v, {2}], Map[zrotate0[#, phi, 0] &, n, {2}]}; zrotate0[{x_, y_, z_}, phi_, trans_] := {Cos[phi]x + Sin[phi]y, -Sin[phi]x + Cos[phi]y, z + trans}; MeshRotateZ[l_, phi_, trans_] := Mesh3D @@ Map[zrotate[#, phi, trans] &, l]; meshscale0[{v_,n_},lm_]:={lm v, n}; MeshScale[ l_, lm_] := Mesh3D@@ Map[meshscale0[#,lm]&,l]; meshflip0[{v_,n_}]:={v, -n}; MeshFlip[ m_] := Mesh3D@@ Map[meshflip0,m]; meshtranslate0[{v_,n_},lm_]:={Map[lm+#&,v,{2}], n}; MeshTranslate[ l_, lm_] := Mesh3D@@ Map[meshtranslate0[#,lm]&,l]; minmax[l_]:={Min[l],Max[l]}; MeshBoundingBox[l_]:=Transpose[minmax/@Transpose[Flatten[Transpose[List@@l][[\ \ 1]],2]]]; reflect[p_,{pn_, d_}]:=p-2(p.pn-d)pn; nreflect[p_,pn_]:=p-2(p.pn)pn; meshreflect0[{v_,n_},pn_, d_]:={Map[reflect[#, {pn, d}] &, v,{2}], \ \ Map[nreflect[#, pn] &, n,{2}]}; MeshReflect[ l_, Plane[pn_, d_]] := Mesh3D@@ Map[meshreflect0[#,pn,d]&,List@@l]; MeshReflectP[ Plane[pn1_, d1_],Plane[pn_, d_]] := Plane[ nreflect[pn1, pn], d1-2d(pn1.pn) ]; MeshFlip[l_, sign_] := Mesh3D @@ Map[meshflip0[#, sign] &, l]; meshflip0[{v_, n_}, sign_] := { v, sign*n}; equalPBC[Plane[n1_,d1_],Plane[n2_,d2_],eps_:0.001]:= norm[Cross[n1,n2]]+(d1-d2)^2