Utilities

MeshReader

ReadMesh[name_]:= Module[
{NodeData,ElData},
NodeData=Take[Import[StringJoin[name,".node"],"Table"],{2,-2}][[All,{2,3}]];
ElData =Take[Import[StringJoin[name,".ele"],"Table"],{2,-2}][[All,{2,3,4}]];
{NodeData, ElData}
]

Back to A

Here is a messh with almost 400 vertices.

Directory[];
SetDirectory["/home/campus03/struther/Desktop/Link to Courses/4620/ClassSampleMathematica/Sep21"];
{NodeData,ElData}=ReadMesh["A400"];
Show[
ListPlot[NodeData],
Graphics[
Table[
{Hue[RandomReal[]],Polygon[NodeData[[ElData[[i]]]] ]},
{i, Length[ElData]}
]
],
Frame->True,
AspectRatio->Automatic
]

Here is a plot of a function over the same mesh

f[{x_,y_}]:=Cos[10 x y];
fNodeData =Table[
{x,y}=NodeData[[i]];
{x,y,f[{x,y}]},
{i, Length[NodeData]}
];
Show[
Graphics3D[
Table[
Polygon[fNodeData[[ElData[[i]]]]],
{i, Length[ElData]}
]
],
Axes->True,
Boxed->True,
BoxRatios->{1,1,0.5},
AxesLabel->{"x","y","z"}
]

We need to be able to to build the stiffness matrix A and load vector b for this mesh. The key is a little bit (in this case one triangle) at a time.

2D Reference Map

So the first thing we need to do is build affine maps to go backwards and forwards between our reference triangle 
	{{0,0},{1,0},{0,1}}
and a general 2D triangle
	{p1,p2,p3}.

As we showed on Friday this is fairly easy and the formula is 
	X->p1+[({{p2-p1, |, p3-p1}})].X

RefToReal2[X_,{p1_,p2_,p3_}]:= p1 + Transpose[{p2-p1,p3-p1}].X
RealToRef2[x_,{p1_,p2_,p3_}]:=Inverse[ Transpose[{p2-p1,p3-p1}]].(x-p1)
{p1,p2,p3}={ {0.2, 0.4},{1.2, 1.4},{-0.2, -0.5}};
pp1=RefToReal2[{0,0},{p1,p2,p3}];
pp2=RefToReal2[{1,0},{p1,p2,p3}];
pp3=RefToReal2[{0,1},{p1,p2,p3}];
{
Show[
Graphics[{
{Pink, Polygon[{p1,p2,p3}]},
Text[1,p1],
Text[2,p2],
Text[3,p3]
}
],
Frame->True,
AspectRatio->Automatic
],
Show[
Graphics[{
{Cyan, Polygon[{pp1,pp2,pp3}]},
Text[1,pp1],
Text[2,pp2],
Text[3,pp3]
}
],
Frame->True,
AspectRatio->Automatic
]
}

The inverse of this map is
	x->(([({{p2-p1, |, p3-p1}})])^(-1)).(x-p1) 
It is easy to check this works.

RealToRef2[x_,{p1_,p2_,p3_}]:=Inverse[ Transpose[{p2-p1,p3-p1}]].(x-p1)
{p1,p2,p3}={ {0.2, 0.4},{1.2, 1.4},{-0.2, -0.5}};
{P1,P2,P3}=Map[RealToRef2[#,{p1,p2,p3}]&,{p1,p2,p3}]
Show[
Graphics[{
{Yellow, Polygon[{P1,P2,P3}]},
Text[1,P1],
Text[2,P2],
Text[3,P3]
}
],
Frame->True,
AspectRatio->Automatic
]

Tent functions are easy on the reference triangle and hard on the real triangle.

Tent[X_]:= 1-{1,1}.X
{
Plot3D[ 
Tent[{X1,X2}],
{X1,0,1},{X2,0,1},
PlotRange->{0,1},
BoxRatios->Automatic],
ParametricPlot3D[ 
{u(1+ v),u(1-  v),Tent[ {u(1+ v),u(1-  v)}]},
{u,0,1/2},{v,-1,1},
BoxRatios->Automatic
]
}

The solution is to pull everything back to the reference triangle and do the integrations there.  For instance, to compute
	(∫_(Subscript[▲, i]))f[x](Subscript[ϕ, i])[x]dx
with the (i^th) triangle described by the vertices {p1,p2,p3} we can compute and then adjust the area by the constant Jacobian of the mapping.

f[{x1_,x2_}]:=12+Sin[ x1 x2]
{p1,p2,p3}={{0.1, 0.3},{0.4,0.6},{0.3, 0.7}}
Abs[Det[Transpose[{p2-p1,p3-p1}]]]NIntegrate[
Evaluate[f[RealToRef2[{X1,X2},{p1,p2,p3}]] Tent[{X1,X2}]],
{X1,0,1},{X2, 0, 1-X1}]

In some sense, the gradients are easier because they are constants.  Make sure you are checking that things are the correct size and plausible.  It is very easy to get the Jacobians upside down.  Is the red correct in the formula above.

In Progress: How do you even tell if you are inside a given triangle?

A={{0,1},{-1,0}};
A.{x,y}

{p1,p2,p3}

Clear[Test]
Test[x_,{p1_,p2_,p3_}]:= Module[
{n12,n13,n23, A={{0,1},{-1,0}}},
n12=A.(p1-p2);
n13=A.(p1-p3);
n23=A.(p2-p3);
And[
Sign[x.n12]==Sign[p3.n12],
Sign[x.n13]==Sign[p2.n13],
Sign[x.n23]==Sign[p1.n23]
]
]
Plot3D[
Boole[Test[{x1,x2},{p1,p2,p3}]],
{x1,0.1,1},
{x2,0.1,1}
]

Test[ (p1+p2+p3)/3,{p1,p2,p3}]
Test[ p1,{p1,p2,p3}]

p1
p2
p3

ListPlot[{p1,p2,p3},
PlotRange->All,
PlotStyle->Red]

p1
p2
p3

3D Reference Map

So the first thing we need to do is build affine maps to go backwards and forwards between our reference tetrahedron 
	{{0,0,0},{1,0,0},{0,1,0},{0,0,1}}
and a general 3D tetrahedron
	{p1,p2,p3,p4}.

The formula is 
	X->p1+[({{p2-p1, |, p3-p1, |, p4-p1}})].X

RefToReal[X_,{p1_,p2_,p3_,p4_}]:= p1 + Transpose[{p2-p1,p3-p1,p4-p1}].X
RealToRef[x_,{p1_,p2_,p3_,p4_}]:=Inverse[ Transpose[{p2-p1,p3-p1,p4-p1}]].(x-p1)
Tet={p1,p2,p3,p4}={ {0.2, 0.4,0.3},{1.2, 1.4,0.1},{-0.2, -0.5,0.2},{0.1,0.2,0.3}};
RefTet={{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
Map[RefToReal[#,{p1,p2,p3,p4}]&, RefTet]==Tet
Map[RealToRef[#,{p1,p2,p3,p4}]&, Tet]==RefTet

MyTet[{p1_,p2_,p3_,p4_}]:= Module[
{Faces},
Faces=Subsets[{p1,p2,p3,p4},{3}];
Map[Polygon, Faces]
]
Show[
Graphics3D[
MyTet[RefTet]
],
Axes->True,
AxesLabel->{"x1","x2","x3"}
]

Subsets

==p1
pp2=RefToReal[{1,0,0},{p1,p2,p3,p4}]==p2
pp3=RefToReal[{0,1,0},{p1,p2,p3,p4}]==

{
Show[
Graphics[{
{Pink, Polygon[{p1,p2,p3}]},
Text[1,p1],
Text[2,p2],
Text[3,p3]
}
],
Frame->True,
AspectRatio->Automatic
],
Show[
Graphics[{
{Cyan, Polygon[{pp1,pp2,pp3}]},
Text[1,pp1],
Text[2,pp2],
Text[3,pp3]
}
],
Frame->True,
AspectRatio->Automatic
]
}

The inverse of this map is
	x->(([({{p2-p1, |, p3-p1}})])^(-1)).(x-p1) 
It is easy to check this works.

RealToRef[x_,{p1_,p2_,p3_}]:=Inverse[ Transpose[{p2-p1,p3-p1}]].(x-p1)
{p1,p2,p3}={ {0.2, 0.4},{1.2, 1.4},{-0.2, -0.5}};
{P1,P2,P3}=Map[RealToRef[#,{p1,p2,p3}]&,{p1,p2,p3}]
Show[
Graphics[{
{Yellow, Polygon[{P1,P2,P3}]},
Text[1,P1],
Text[2,P2],
Text[3,P3]
}
],
Frame->True,
AspectRatio->Automatic
]