
(* 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<eps;

projection[a_,b_]:=Outer[Times,b,a];

    

(* Conversion to display graphics formats *)

insertValues[f_, v_] := Map[v[[#]] &, f];

conv2D3D[l_]:=Map[Append[#,0]&,l];



Domain /: MakeBoxes[Domain[___], StandardForm] := 
    RowBox[{"\[SkeletonIndicator]", "Domain", "\[SkeletonIndicator]"}];

Mesh3D /: MakeBoxes[Mesh3D[___], StandardForm] := 
    RowBox[{"\[SkeletonIndicator]", "Mesh3D", "\[SkeletonIndicator]"}];
    
(*SphereCarpet /: MakeBoxes[SphereCarpet[___], StandardForm] := 
    RowBox[{"\[SkeletonIndicator]", "SphereCarpet", \
\
"\[SkeletonIndicator]"}];*)

LineArray /:MakeBoxes[LineArray[___], StandardForm] := 
    RowBox[{"\[SkeletonIndicator]", "LineArray", "\[SkeletonIndicator]"}];


Edges /: MakeBoxes[Edges[___], StandardForm] := 
    RowBox[{"\[SkeletonIndicator]", "Edges", "\[SkeletonIndicator]"}];


DomainToGraphics[v_] := 
    Graphics[Join[Map[Line,List@@ v],Map[Line, Transpose[List@@v]]]];
    
    
trp[{{a_,b_},{c_,d_}}]:=Polygon[{c,d,b,a}];

Mesh3DToGraphics3D[msh_] := 
    Graphics3D[Flatten[Map[
    Map[trp,Flatten[
          Map[Partition[#,2,1]&,Transpose[Map[Partition[#,2,1]&,#[[1]]]]],1]]&\
\
,List@@msh],1]]

      
End[];

EndPackage[];
