BeginPackage["BVP2`"];

ReadMesh::usage="{nodes,elements,segments}=ReadMesh[filename]
Read all information for a mesh generated by EasyMesh
filename is a string with the name of the mesh files without extension, 
e.g. EasyMesh/test2";

FEMSolve::usage="points=FEMSolve[a, b, f, gDirichlet, gNeumann,nodes,elements,segments];
Solve a second order boundary value problem";

FEMGradient::usage="gradient=FEMGradient[points,elements];
compute the values of the gradient at the nodal points";

MeshGraphics::usage="gmesh=MeshGraphics[nodes,elements] generates a graphics object to visualize the mesh by Show[gmesh]";
FEMSurface::usage="g=FEMSurface[points,elements] generates a 3D graphics object to be displayed with Show[g]";

FEMSurfaceNoMesh::usage="g=FEMSurfaceNoMesh[points,elements] generates a 3D graphics object to be displayed with Show[g]\n the mesh is not shown";

LevelGraphics2D::usage="glevel=LevelGraphics2D[points,elements,levels] generates a graphics object to visualize the mesh by Show[glevel]";

LevelGraphics3D::usage="glevel=LevelGraphics3D[points,elements,levels] generates a graphics object to visualize the mesh by Show[glevel]";


(*
uses the globale variables  nodes, elements, segments and Node2Degree
*)

Begin["`Private`"];

Aglobal={}
bglobal={};


ReadMesh[filename_String]:=Module[{instream,n,nodes,elements,segements},
		(* read all the nodes *)
		instream=OpenRead[StringJoin[filename,".n"]];
n=Read[instream,Number];
nodes=Table[{},{n}];
For[k=1,k<=n,k+=1,
	     tmp=Read[instream,{Number,Character,Real,Real,Number}];
      nodes[[k]]=Drop[tmp,{2}];
   	nodes[[k,1]]+=1; (* numbering starts at 1 *)
      ];
Close[instream];
		
		(* read the elements *)
instream=OpenRead[StringJoin[filename,".e"]];
n=Read[instream,Number];
elements=Table[{},{n}];
For[k=1,k<=n,k+=1,
	     tmp=Read[instream,{Number,Character,
				Number,Number,Number,
				Number,Number,Number,
				Number,Number,Number,
				Real,Real,Number}];
      elements[[k]]={tmp[[1]]+1,tmp[[3]]+1,tmp[[4]]+1,tmp[[5]]+1,Last[tmp]}
			    (* numbering starts at 1 *)
];
Close[instream];

(* read all segments *)

instream=OpenRead[StringJoin[filename,".s"]];
n=Read[instream,Number];
segments=Table[{},{n}];
For[k=1,k<=n,k+=1,
	     tmp=Read[instream,{Number,Character,
				Number,Number,Number,
				Number,Number}];
      segments[[k]]=Drop[Drop[tmp,2],{3,4}];
	    segments[[k,1]]+=1; 
	 segments[[k,2]]+=1; (* numbering starts at 1 *)
];
Close[instream];
segments=Select[segments,Last[#]>0&];
(* return all results *)
{nodes,elements,segments}]


(*
Numbering of the degrees of freedom: FindDOF[]
*)
FindDOF[nodes_]:=Module[{counter,Node2Degree},
counter=1;
Node2Degree=Table[0,{Length[nodes]}];
For[k=1,k<=Length[nodes],k+=1,
   If[nodes[[k,4]]==0, Node2Degree[[k]]=counter;counter+=1;];
   If[nodes[[k,4]]==2, Node2Degree[[k]]=counter;counter+=1;]
];
{counter-1,Node2Degree}];


(*
Integration over one triangle: ElementContribution[]
*)
ElementContribution[corners_List,aCoeff_Function,bCoeff_Function,fCoeff_Function]:=
Module[{area,a,b,c,M},
a=Flatten[{corners[[2]]-corners[[1]],0}];
b=Flatten[{corners[[3]]-corners[[1]],0}];
c=Cross[a,b];
area=Abs[c[[3]]]/2;
M={{corners[[3,2]]-corners[[2,2]],corners[[1,2]]-corners[[3,2]],corners[[2,2]]-corners[[1,2]]},
   {corners[[2,1]]-corners[[3,1]],corners[[3,1]]-corners[[1,1]],corners[[1,1]]-corners[[2,1]]}};
{(aCoeff[corners[[1,1]],corners[[1,2]]]+aCoeff[corners[[2,1]],corners[[2,2]]]+
  aCoeff[corners[[3,1]],corners[[3,2]]])/(12*area)*
   Transpose[M].M+
 area/3* {{bCoeff[corners[[1,1]],corners[[1,2]]],0,0},
          {0,bCoeff[corners[[2,1]],corners[[2,2]]],0},
          {0,0,bCoeff[corners[[3,1]],corners[[3,2]]]}},
 area/3* {fCoeff[corners[[1,1]],corners[[1,2]]],
          fCoeff[corners[[2,1]],corners[[2,2]]],
          fCoeff[corners[[3,1]],corners[[3,2]]]}
}
]

(*
Integration over one edge: ElementConributionEdge[]
*)
(*
ElementContributionEdge[corners_List,g_Function]:=
Module[{},
Sqrt[(corners[[1,1]]-corners[[2,1]])^2+(corners[[1,2]]-corners[[2,2]])^2]/2*
    {g[corners[[1,1]],corners[[1,2]]] , g[corners[[2,1]],corners[[2,2]]] }
]
*)

ElementContributionEdge[corners_List,g_Function]:=
Module[{p1,p2,m,w},
  m=(corners[[1]]+corners[[2]])/2;
  p1=m+(corners[[1]]-corners[[2]])/Sqrt[3]/2;
  p2=m-(corners[[1]]-corners[[2]])/Sqrt[3]/2;
  w=N[(1-1/Sqrt[3])/2];
Sqrt[(corners[[1,1]]-corners[[2,1]])^2+	(corners[[1,2]]-corners[[2,2]])^2]/2*
    {g[p1[[1]],p1[[2]]]*(1-w)+g[p2[[1]],p2[[2]]]*w,
     g[p1[[1]],p1[[2]]]*w    +g[p2[[1]],p2[[2]] ]*(1-w)}
]



(*
Assembling the equations and solving
*)
insertElement[element_List,aCoeff_Function,bCoeff_Function,fCoeff_Function,
    gCoeff_Function,nodes_List,Node2Degree_List] :=
Block[{elMat,elVec,corners,dofs,tMat},
corners={{nodes[[element[[1]],2]],nodes[[element[[1]],3]]},
         {nodes[[element[[2]],2]],nodes[[element[[2]],3]]},
         {nodes[[element[[3]],2]],nodes[[element[[3]],3]]}};
{elMat,elVec}=ElementContribution[corners,aCoeff,bCoeff,fCoeff];
dofs={Node2Degree[[element[[1]]]],Node2Degree[[element[[2]]]],Node2Degree[[element[[3]]]]};
For[k1=1,k1<=3,k1+=1,
   If[dofs[[k1]]>0,
     For[k2=1,k2<=3,k2+=1,
        If[dofs[[k2]]>0,Aglobal[[dofs[[k1]],dofs[[k2]]]]+=elMat[[k1,k2]];,
        (*else*)    						
           bglobal[[dofs[[k1]]]]+= 
           (elMat[[k1,k2]]*gCoeff[corners[[k2,1]],corners[[k2,2]] ]);
        ](*if k2*)
     ](*for k2*)
;
     bglobal[[dofs[[k1]]]]+= elVec[[k1]];
     ](*if k1*);
](*for k1*);
](*Block*)


insertElementEdge[edge_List,g_Function,nodes_List,Node2Degree_List] :=
Block[{elVec,corners,dofs,k},
corners={{nodes[[edge[[1]],2]],nodes[[edge[[1]],3]]},
         {nodes[[edge[[2]],2]],nodes[[edge[[2]],3]]}};
elVec=ElementContributionEdge[corners,g];
dofs={Node2Degree[[edge[[1]]]],	Node2Degree[[edge[[2]]]]};
For[k=1,k<=2,k+=1,
   If[dofs[[k]]>0,bglobal[[dofs[[k]]]] -= elVec[[k]] ];
](*for k*);
](*Block*)


FEMSolve[a_Function,b_Function,f_Function,g1_Function,g2_Function,nodes_List,elements_List,segments_List]:=
   Block[{uglobal,dof,Node2Degree},
   {dof , Node2Degree} = FindDOF[nodes];
Aglobal=Table[0,{dof},{dof}];
bglobal=Table[0,{dof}];
For[k=1,k<=Length[elements],k+=1,
   insertElement[Take[elements[[k]],{2,4}],a,b,f,g1,nodes,Node2Degree]] ;
For[k=1,k<=Length[segments],k+=1,
   If[segments[[k,3]]==2,insertElementEdge[Take[segments[[k]],2],g2,nodes,Node2Degree]]];
uglobal=LinearSolve[Aglobal,-bglobal];
Table[{nodes[[k,2]],nodes[[k,3]],
      Which[Node2Degree[[k]]>0,uglobal[[Node2Degree[[k]]]],
           True,g1[nodes[[k,2]],nodes[[k,3]]]]},
     {k,Length[nodes]}]
]

FEMGradient[points_,elements_]:=
  Module[{area,gradpoints,a,b,c,corners,grad,angle},
    gradpoints=Table[Flatten[{points[[k,{1,2}]],0,0,0}],{k,Length[points]}];
    For[k=1,k<=Length[elements],k++,
	corners={points[[elements[[k,2]],{1,2}]],
		 points[[elements[[k,3]],{1,2}]],
		 points[[elements[[k,4]],{1,2}]]};
      a=Flatten[{corners[[2]]-corners[[1]],0}];
      b=Flatten[{corners[[3]]-corners[[1]],0}];
      c=Cross[a,b];
      area=Abs[c[[3]]]/2;
      M={{corners[[3,2]]-corners[[2,2]],corners[[1,2]]-corners[[3,2]],
            corners[[2,2]]-corners[[1,2]]},
         {corners[[2,1]]-corners[[3,1]],corners[[3,1]]-corners[[1,1]],
            corners[[1,1]]-corners[[2,1]]}};
      grad=-1/(2*area)*M.
	{points[[elements[[k,2]],3]],points[[elements[[k,3]],3]],
	 points[[elements[[k,4]],3]]};
      angle={0,0,0};
      angle[[1]]=ArcCos[(a.b)/
            (Sqrt[a[[1]]^2+a[[2]]^2]*Sqrt[b[[1]]^2+b[[2]]^2])];
      c=b-a;
      angle[[2]]=ArcCos[-(a.c)/
            (Sqrt[a[[1]]^2+a[[2]]^2]*Sqrt[c[[1]]^2+c[[2]]^2])];
      angle[[3]]=N[Pi-angle[[1]]-angle[[2]]];
      gradpoints[[elements[[k,2]],{3,4}]]+=grad*angle[[1]];
      gradpoints[[elements[[k,2]],5]]+=angle[[1]];
      gradpoints[[elements[[k,3]],{3,4}]]+=grad*angle[[2]];
      gradpoints[[elements[[k,3]],5]]+=angle[[2]];
      gradpoints[[elements[[k,4]],{3,4}]]+=grad*angle[[3]];
      gradpoints[[elements[[k,4]],5]]+=angle[[3]];		
   ];
   For[k=1,k<=Length[points],k++, 
      gradpoints[[k]]= gradpoints[[k,{3,4}]]/gradpoints[[k,5]]];
   gradpoints
]


MeshGraphics[nodes_,elements_]:=
	Module[{},
	Graphics[Table[
	Line[{{nodes[[elements[[k,2]] ]][[2]],nodes[[elements[[k,2]] ]][[3]]},
    {nodes[[elements[[k,3]] ]][[2]],nodes[[elements[[k,3]] ]][[3]]},
    {nodes[[elements[[k,4]] ]][[2]],nodes[[elements[[k,4]] ]][[3]]},
    {nodes[[elements[[k,2]] ]][[2]],nodes[[elements[[k,2]] ]][[3]]}}],
	{k,1,Length[elements]}]] 
]

FEMSurface[points_,elements_]:=
Module[{},
Table[Graphics3D[
      Polygon[{points[[elements[[k,2]]]],points[[elements[[k,3]]]],
          points[[elements[[k,4]]]]}]],{k,Length[elements]}]
]


FEMSurfaceNoMesh[points_,elements_]:=
Module[{},
Table[Graphics3D[{EdgeForm[],
      Polygon[{points[[elements[[k,2]]]],points[[elements[[k,3]]]],
              points[[elements[[k,4]]]]}    ]}],
      {k,Length[elements]}]
]


LevelGraphics2D[points_,elements_,levels_]:=
	Module[{v1,v2,v3,lines,inside,coord},
coord={1,2}; (* allows for identical code for 2D and 3D graphics *)
lines={};
inside=False;
Do[ (* j loop over all levels*)
Do[ (* i loop over all elements*)
  v1=points[[elements[[i,2]],3]];
  v2=points[[elements[[i,3]],3]];
  v3=points[[elements[[i,4]],3]];
  If[levels[[j]]<Min[{v1,v2,v3}] || levels[[j]]>Max[{v1,v2,v3}],Null,
    (* there is a level curve *)
    Which[
      v1==v2&&v2==v3,Null,
      v1==v2&&v1==levels[[j]],
             AppendTo[lines,{points[[elements[[i,2]],coord]],
                             points[[elements[[i,3]],coord]]} ],
      v1==v3&&v1==levels[[j]],
             AppendTo[lines,{points[[elements[[i,2]],coord]],
                             points[[elements[[i,4]],coord]]} ],
      v2==v3&&v2==levels[[j]],
             AppendTo[lines,{points[[elements[[i,3]],coord]],
                             points[[elements[[i,4]],coord]]} ],
      True, inside=True
      ];
    If[inside,
      If[v1!=v2,s1=(levels[[j]]-v2)/(v1-v2),s1=2];    
      If[v2!=v3,s2=(levels[[j]]-v3)/(v2-v3),s2=2];    
      If[v3!=v1,s3=(levels[[j]]-v1)/(v3-v1),s3=2];
      Which[
        0<=s1&&s1<=1,
           Which[
             0<=s2&&s2<=1, (* segments 12 and 23 *)
               AppendTo[lines,
                 {s1*points[[elements[[i,2]],coord]]+
                   (1-s1)*points[[elements[[i,3]],coord]],
                  s2*points[[elements[[i,3]],coord]]+
                  (1-s2)*points[[elements[[i,4]],coord]]}],
             True,         (* segements 12 and 31 *)
               AppendTo[lines,
                 {s1*points[[elements[[i,2]],coord]]+
                   (1-s1)*points[[elements[[i,3]],coord]],
                  s3*points[[elements[[i,4]],coord]]+
                  (1-s3)*points[[elements[[i,2]],coord]]}]
             ],
        True, (* segement 23 and 31 *)
          AppendTo[lines,
            {s2*points[[elements[[i,3]],coord]]+
             (1-s2)*points[[elements[[i,4]],coord]],
            s3*points[[elements[[i,4]],coord]]+
             (1-s3)*points[[elements[[i,2]],coord]]}]
        ]
      ];
  ]
,{i,Length[elements]}]
,{j,Length[levels]}];

lines=Map[Line,lines];
Map[Graphics,lines]
]

LevelGraphics3D[points_,elements_,levels_]:=
	Module[{v1,v2,v3,lines,inside,coord},
coord={1,2,3}; (* allows for identical code for 2D and 3D graphics *)
lines={};
inside=False;
Do[ (* j loop over all levels*)
Do[ (* i loop over all elements*)
  v1=points[[elements[[i,2]],3]];
  v2=points[[elements[[i,3]],3]];
  v3=points[[elements[[i,4]],3]];
  If[levels[[j]]<Min[{v1,v2,v3}] || levels[[j]]>Max[{v1,v2,v3}],Null,
    (* there is a level curve *)
    Which[
      v1==v2&&v2==v3,Null,
      v1==v2&&v1==levels[[j]],
             AppendTo[lines,{points[[elements[[i,2]],coord]],
                             points[[elements[[i,3]],coord]]} ],
      v1==v3&&v1==levels[[j]],
             AppendTo[lines,{points[[elements[[i,2]],coord]],
                             points[[elements[[i,4]],coord]]} ],
      v2==v3&&v2==levels[[j]],
             AppendTo[lines,{points[[elements[[i,3]],coord]],
                             points[[elements[[i,4]],coord]]} ],
      True, inside=True
      ];
    If[inside,
      If[v1!=v2,s1=(levels[[j]]-v2)/(v1-v2),s1=2];    
      If[v2!=v3,s2=(levels[[j]]-v3)/(v2-v3),s2=2];    
      If[v3!=v1,s3=(levels[[j]]-v1)/(v3-v1),s3=2];
      Which[
        0<=s1&&s1<=1,
           Which[
             0<=s2&&s2<=1, (* segments 12 and 23 *)
               AppendTo[lines,
                 {s1*points[[elements[[i,2]],coord]]+
                   (1-s1)*points[[elements[[i,3]],coord]],
                  s2*points[[elements[[i,3]],coord]]+
                  (1-s2)*points[[elements[[i,4]],coord]]}],
             True,         (* segements 12 and 31 *)
               AppendTo[lines,
                 {s1*points[[elements[[i,2]],coord]]+
                   (1-s1)*points[[elements[[i,3]],coord]],
                  s3*points[[elements[[i,4]],coord]]+
                  (1-s3)*points[[elements[[i,2]],coord]]}]
             ],
        True, (* segement 23 and 31 *)
          AppendTo[lines,
            {s2*points[[elements[[i,3]],coord]]+
             (1-s2)*points[[elements[[i,4]],coord]],
            s3*points[[elements[[i,4]],coord]]+
             (1-s3)*points[[elements[[i,2]],coord]]}]
        ]
      ];
  ]
,{i,Length[elements]}]
,{j,Length[levels]}];

lines=Map[Line,lines];
Map[Graphics3D,lines]
]

End[];
EndPackage[];
