Tent Functions

TriangularSurfacePlot

First a picture of a linear interpolation of a function on a triangulated domain using some Mathematica functions.

Clear[a,b,c,d,f,nInt,nEdge,Δx,Δy]
{{a,b},{c,d}}={{7.0,18.0},{1.0,3.0}};
nInt=40;
nEdge=3;
{Δx,Δy}={b-a,d-c}/nEdge;

InteriorData={};
While[
Length[InteriorData]<nInt,
x=RandomReal[{a,b}];
y=RandomReal[{c,d}];
AppendTo[InteriorData,{x,y}]
]
EdgeData=Join[
Table[ {a,y},{y,c+Δy,d-Δy,Δy}],
Table[ {b,y},{y,c+Δy,d-Δy,Δy}],
Table[ {x,c},{x,a+Δx,b-Δx,Δx}],
Table[ {x,d},{x,a+Δx,b-Δx,Δx}]
];
CornerData={{a,c},{a,d},{b,c},{b,d}};
Data=Join[InteriorData,EdgeData,CornerData];
ListPlot[
{InteriorData,EdgeData,CornerData},
PlotStyle->PointSize[0.02]
]

What does a liner interpolation of  function look like on this rectangular domain?

Needs["ComputationalGeometry`"]
Clear[x,y,f]
f[{x_,y_}]=x Sin[x + y];
HeightData=Table[
{x,y}=Data[[i]];
{x,y,f[{x,y}]},
{i,Length[Data]}
];
TriangularSurfacePlot[HeightData,
BoxRatios->{1,1,1}
]

Triangulation Structure

What do the underlying triangles look like?

PlanarGraphPlot[Data]

Triangulation Description

How do you descibe them!

DelaunayTriangulation[Data]

Reference Triangle

Our target is to linearly map a domain triangle {{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])},{(Subscript[x, 3]),(Subscript[y, 3])}} onto the reference triangle{{0,0},{0,1},{1,0}} by mapping each of the vertices {(Subscript[x, i]),(Subscript[y, i])} to the matching vertex of the reference triangle. The following takes the reference triangle onto a general domain triangle.

FromRef[ {p1_,p2_,p3_},{X_,Y_}]:= Module[
{A},
A= Transpose[ {p2-p1,p3-p1}];
A.{X,Y}
]

Here is a verification!

{
ToRef[{{x1,y1},{x2,y2},{x3,y3}},{0,0}],
ToRef[{{x1,y1},{x2,y2},{x3,y3}},{1,0}],
ToRef[{{x1,y1},{x2,y2},{x3,y3}},{0,1}]
}

We will need the inverse transformation as well.  Since it is written in the matrix form all we need to do is invert the matrix A.

ParametricPlot3D[
s{1,0,t},
{s,0,1},
{t,0,1}
]

Clear[f,x,y,ϵ]
ϵ=0.00000000000;
f[{x_,y_}]= (Sqrt[9-(x^2)-(y^2)])
ParametricPlot3D[
{x,y}=FromRef[{{0,0}, {1,1},{2,0}},{s,s*t}];
{x,y,f[{x,y}]},
{s,ϵ,1-ϵ},
{t,ϵ,1-ϵ},
MaxRecursion->0,
BoxRatios->{1,1,1},
PlotRange->All,
MeshShading->{{Red},{Green}}
]

?ParametricPlot3D

?MeshShading

GridP

x2 y1-x3 y1-x1 y2+x3 y2+x1 y3-x2 y3

{X,Y}/.Solve[
ToRef[{{x1,y1},{x2,y2},{x3,y3}},{X,Y}]=={x,y},
{X,Y}
][[1]]//TableForm

Initializations

The following utilities are to provide automatic means to build stiffness matrices etc. from triangulations. The points are specified in a list of points in (ℝ^n) and then the connectivity is specified in a seperate structure.  Mathematica  likes to use adjacency lists and incidence matrices.  It will be simpler for us to use a specific list of triangles.  For example, {{1,2,3},{3,4,5}} would specify a triangulation where the only triangles are {1,2,3} and {3,4,5}.

AdjacencyListToTriangles

The idea is simple to make a list of all the triangles from an adjacency list:
A) For each vertex get the list of adjacent vertices. 
B) Look for triangles using the adjacency lists to those vertices. 
C)Accumulate the triangles. 
D) Delete the vertex you examined
E) Repeat

Needs["ComputationalGeometry`"]
AdjacentTriangles[ i_,AdjListIn_,IgnoreIn_]:=Module[
{j,k,p,AdjList,Triangles},
{p,AdjList}=AdjListIn[[i]];
AdjList=Complement[AdjList,IgnoreIn]; 
Triangles=AdjListIn[[AdjList,2]];
Triangles=Map[Intersection[AdjList,#]&,Triangles,{1}];
Union[
Flatten[
Table[
Table[
Sort[{p,AdjList[[j]],Triangles[[j]][[k]]}],
{k,Length[Triangles[[j]]]}],
{j, Length[Triangles]}
],
1]
]
]

Needs["ComputationalGeometry`"]
AdjacencyListToTriangles[ AdjListIn_]:=Module[
{i,CompletedVertices,Temp},
CompletedVertices={};
Temp=Table[
CompletedVertices=Append[CompletedVertices,i];
AdjacentTriangles[i,AdjListIn,CompletedVertices],
{i,Length[AdjListIn]}
];
Apply[Union,Temp]
]

CheckTrianglePlot

I need to be able to check that I have done my triangles correctly.  The following plots the Delaunay triangulation using planar graph plot and then superimposes my triangulation on it.

Clear[CheckTrianglePlot,InPoints]
CheckTrianglePlot[InPoints_]:=Module[
{MyTriang,MyPic,
p1,p2,p3,
i},
MyTriang=AdjacencyListToTriangles[DelaunayTriangulation[InPoints]];
MyPic=Table[
{p1,p2,p3}=Part[InPoints,MyTriang[[i]]];
{Hue[(i/(Length[MyTriang]))],Polygon[{p1,p2,p3}]},
{i, Length[MyTriang]}
];
Show[
Graphics[MyPic],
PlanarGraphPlot[InPoints,
Frame->True],
Frame->True
]
]

Testing the bits so far.

Points=Table[
{RandomReal[], RandomReal[]},
{6}];
DelAdjList=DelaunayTriangulation[Points];
MyTriang=AdjacencyListToTriangles[DelAdjList]
CheckTrianglePlot[Points]

Triangle Reference Map

The mapping from the reference triangle onto a triangle specified by {{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])},{(Subscript[x, 3]),(Subscript[y, 3])}}={(Subscript[p, 1]),(Subscript[p, 2]),(Subscript[p, 3])} is fairly straight forward.

Clear[FromRef]
FromRef[P_,{p1_,p2_,p3_}]:= Module[
{A,p},
p=p1;
A=Transpose[{p3-p,p2-p}];
A.P+p
]

Testing this is easy.

FromRef[{0,0},{{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])},{(Subscript[x, 3]),(Subscript[y, 3])}}]
FromRef[{0,1},{{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])},{(Subscript[x, 3]),(Subscript[y, 3])}}]
FromRef[{1,0},{{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])},{(Subscript[x, 3]),(Subscript[y, 3])}}]

The area tranformation is the Jacobian of the reference map or in other words the determinant of the matrix A is

Clear[Jac]
Jac[{p1_,p2_,p3_}]:= Module[
{A},
A=Transpose[{p3-p1,p2-p1}];
Det[A]
]

Testing the code

Jac[{{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])},{(Subscript[x, 3]),(Subscript[y, 3])}}]//Simplify

Creating Slab Fns

The slab function has a simple representation on the reference triangle of 1-X-Y.  All we have to do is move it to the current space using the Reference Map

Clear[RefSlab,TriangulationSlab]
RefSlab[{X_,Y_}]:=If[0<X<1∧0<Y<1-X,(1-X-Y)]
TriangulationSlab[P_,{p1_,p2_,p3_}]:=Append[FromRef[P,{p1,p2,p3}],RefSlab[P]]

Checking our computations. First the reference slab!

Plot3D[RefSlab[{X,Y}],{X,-0.1,1.1},{Y,-0.1,1.1},
PlotRange->All,
PlotPoints->40,
AxesLabel->{"X","Y","Z"}
]

Now the current version plotted from the Reference configuration.

ParametricPlot3D[
TriangulationSlab[{X,Y},{{1,3},{2,4},{-1,3}}],
{X,-0.1,1.1},{Y,-0.1,1.1},
PlotRange->All,
PlotPoints->40,
AxesLabel->{"x","x","x"},
BoxRatios->{1,1,1}
]

Rectangular Parametrization of the Reference Triangle

There is an alternative parmeterization of the reference triangle in terms of S=X+Y and θ=ArcTan[Y/X].  The inverse mapping is X=S/(1+Tan[θ]) and Y=S Tan[θ]/(1+Tan[θ]) which is useful.  The parametr range is 0<S<1and 0<θ<π/2.  You can get away with using 2 data points in each direction because the slab functions are flat.

XR[S_,θ_]:= S/(1+Tan[θ]);
YR[S_,θ_]:= S Tan[θ]/(1+Tan[θ])

Checking the computation.

ParametricPlot3D[{XR[S,θ],YR[S,θ],1-XR[S,θ]-YR[S,θ]},
{S,0,1},{θ,0,π/2},
PlotPoints->2]

Plotting Tent Functions

Given a Triangulation the TentFunction for a node is just a bunch of slab functions.

Points=Table[
{RandomReal[], RandomReal[]},
{8}];
DTriang=DelaunayTriangulation[Points];
MyTriangles=AdjacencyListToTriangles[DTriang]
CheckTrianglePlot[Points]

Now the slab function for a triangle and node on that triangle is fairly easy to create.  To build the tent function one needs to know which triangles meet at that vertex and then simply add up the appropriate slab functions.

i=5;
MatchingTriangles =AdjacentTriangles[i,DTriang,{}]

Sorting them so that i is first

MatchingTriangles=Table[
Sort[MatchingTriangles[[i]],Abs[#1-5]<Abs[#2-5]&],
{i,Length[MatchingTriangles]}
]

Here is the tent function using the X, Y parameterization of the reference triangle.

ParametricPlot3D[
Table[
TriangulationSlab[{X,Y},Points[[MatchingTriangles[[i]]]]],
{i, Length[MatchingTriangles]}],
{X,-0.1,1.1},{Y,-0.1,1.1},
PlotRange->All,
AxesLabel->{"x","y","z"},
BoxRatios->{1,1,1}
]

It took a long time to plot and does not quite match cause it used the rectangular parameterization and had to work out what to do with all the Nulls.  here it is with the S, θ parametrization.

π/2.

ParametricPlot3D[
Table[
TriangulationSlab[{(S/(1+Tan[θ])),S ((Tan[θ])/(1+Tan[θ]))},Points[[MatchingTriangles[[i]]]]],
{i, Length[MatchingTriangles]}],
{θ,0.01,1.57},{S,0.01,0.99},
PlotRange->All,
AxesLabel->{"x","y","z"},
BoxRatios->{1,1,1}
]

I cheated in a tiny bit on the rane to avoid the edges and a few Nulls.

Building the Stiffness Matrix A

The stiffness matrix is defines by
(Subscript[A, i,j])=(∫_Ω)∇(Subscript[ϕ, i]).∇(Subscript[ϕ, j])
where (Subscript[ϕ, i]) is the tent function at the ith node. It should be noted that ∇ϕ is a constant vector on each of the faces and that provide i≠j the integral reduces to an integral over two of the base triangles. Thus this does not require a real integration.

Since we have a representation of the slab function for each triangle in terms of the reference slab function 1-X-Y we can compute the required gradients explicitly using the chain rule.

Gradϕ[{p1_,p2_,p3_}]:=Module[
{A},
AInv=Inverse[{p3-p1,p2-p1}];
-AInv.{1,1}
]

Table[{(Subscript[x, i]),(Subscript[y, i])},{i,3}]

(Subscript[gϕ, 1])=Gradϕ[{{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])},{(Subscript[x, 3]),(Subscript[y, 3])}}]
(Subscript[gϕ, 2])=Gradϕ[{{(Subscript[x, 2]),(Subscript[y, 2])},{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 3]),(Subscript[y, 3])}}]
Simplify[(Subscript[gϕ, 1]).(Subscript[gϕ, 2])]

The area transformation is given by the Jacobian which is the same as the thing sitting on the bottom of the fraction.

Building the Force Fector F

The stiffness matrix is defines by
(Subscript[F, i])=(∫_Ω)f (Subscript[ϕ, i])
where (Subscript[ϕ, i]) is the tent function at the ith node.  This one does require a real integration that we do in the reference triangles.

The integral is the sum of integrals of slab functions time f over the individual triangles.  So if we can do a single slab function we can compute the F

Dec 18th

Plot3D[
1-{1,1}.{u1,u2},
{u1,-1,1},
{u2,-1,1},
RegionFunction->Function[{u1,u2,z},And[0<u1<1, 0<u2<1-u1]]
]

{p1,p2,p3}={{0,0},{1, 1},{3, 2}}
A= Transpose[{
p2-p1,
p3-p1
}
];

A.{0,0}+p1
A.{1,0}+p1
A.{0,1}+p1

Inverse[A].(p1-p1)
Inverse[A].(p2-p1)
Inverse[A].(p3-p1)

Show[
Plot3D[
{
1-{1,1}.Inverse[A].({x1,x2}-p1),
{0,1}.Inverse[A].({x1,x2}-p1)
},
{x1,-6,6},
{x2,-6,6},
RegionFunction->Function[{x1,x2,z},
{u1,u2}=Inverse[A].({x1,x2}-p1);
And[0<u1<1, 0<u2<1-u1]]
],
Graphics3D[{
PointSize[0.02],
Point[Append[p1,1]],
Red,
Point[Append[p2,0]],
Point[Append[p3,0]]
}
],
PlotRange->All
]