Element Spaces

Currently the only finite element spaces we have are piecewise linear.  We want to have some more choices.  The key to thinking about the spaces is the connection between a basis for the space (which we can consider in a reference if it is convenient) and the nodal (hopefully physical) values we want to use.

The trick is to make sure that there is a one to one mapping.  This means that we need to think about invertibility of mappings to do anything more complicated than piecewise linear.

We are going to revisit the piecewise linear example from this new perspective

2D: Triangulation

First we need to make the forwrad and backwards maps from a general triangle.

RefToReal[X_,{p1_,p2_,p3_}]:= p1 + Transpose[{p2-p1,p3-p1}].X
RealToRef[x_,{p1_,p2_,p3_}]:=Module[
{A= Transpose[{p2-p1,p3-p1}],a,b,c,d},
{{a,b},{c,d}}=A;
(1/(a d - b c)){{d,-b},{-c,a}}.(x-p1)
]

Piecewise linears

The peicewise linear functions are defined (as before) in terms of the reference slab function.

RefSlab[{X1_,X2_}]:=1 - X1 - X2

Here is a picture of one of them.

Triangle= RandomReal[{-1,1},{3,2}];Chop[Simplify[RefToReal[RealToRef[{x1,x2},Triangle],Triangle]]];
Chop[Simplify[RealToRef[RefToReal[{X1,X2},Triangle],Triangle]]];
ParametricPlot3D[{
Append[RefToReal[{X1,X2},Triangle],RefSlab[{X1,X2}]],
Append[RefToReal[{X1,X2},Triangle],0]
},
{X1,0,1},{X2,0,1},
RegionFunction->Function[{x1,x2,u,X1,X2},(X1+X2<1)],
BoxRatios->{1,1,0.5},
BoundaryStyle->Directive[Thick,Red]
]

Here is a picture of the three of them

Triangle= RandomReal[{-1,1},{3,2}];
ParametricPlot3D[{
Append[RefToReal[{X1,X2},RotateLeft[Triangle,0]],RefSlab[{X1,X2}]],
Append[RefToReal[{X1,X2},RotateLeft[Triangle,1]],RefSlab[{X1,X2}]],
Append[RefToReal[{X1,X2},RotateLeft[Triangle,2]],RefSlab[{X1,X2}]],
Append[RefToReal[{X1,X2},Triangle],0]
},
{X1,0,1},{X2,0,1},
RegionFunction->Function[{x1,x2,u,X1,X2},(X1+X2<1)],
BoxRatios->{1,1,0.5},
BoundaryStyle->Directive[Thick,Green]
]

One of the three slab functions is more complicated than the others. We picked 1-X1-X2 but X1 and X2 are easier!

EasyRefSlab[{X1_,X2_}]:= X1
Triangle= RandomReal[{-1,1},{3,2}];
ParametricPlot3D[{
Append[RefToReal[{X1,X2},RotateLeft[Triangle,0]],EasyRefSlab[{X1,X2}]],
Append[RefToReal[{X1,X2},RotateLeft[Triangle,1]],EasyRefSlab[{X1,X2}]],
Append[RefToReal[{X1,X2},RotateLeft[Triangle,2]],EasyRefSlab[{X1,X2}]],
Append[RefToReal[{X1,X2},Triangle],0]
},
{X1,0,1},{X2,0,1},
RegionFunction->Function[{x1,x2,u,X1,X2},(X1+X2<1)],
BoxRatios->{1,1,0.5},
BoundaryStyle->Directive[Thick,Green]
]

Piecewise Quadratics

From now on we are going to use the easy slab function.  In terms of the reference triangle there are two flavors of quadratics: A) The quadratics that are 1 at one midpoint, zero at the other midpoints and zero at the corners. B) The quadratics that are 1 at one corner, zero at the other corners and zero at the midpoints.

RefSlab[{X1_,X2_}]:=X1
RefQuad[{X1_,X2_},Midpt]:=4X1(1-X1-X2)
RefQuad[{X1_,X2_},Corner]:= (-1+X1+X2) (-1+2 X1+2 X2)

Here is a picture of these three on the reference triangle.

Plot3D[{
RefSlab[{X1,X2}],
RefQuad[{X1,X2},Corner],
RefQuad[{X1,X2},Midpt]
},
{X1,0,1},
{X2,0,1},
RegionFunction->Function[{X1,X2,z},X1+X2<1],
PlotRange->All,
BoundaryStyle->Green
]

Here is a picture of all six quadratics on a general triangle.

Triangle= RandomReal[{-1,1},{3,2}];
ParametricPlot3D[{
Append[RefToReal[{X1,X2},RotateLeft[Triangle,0]],RefQuad[{X1,X2},Midpt]],
Append[RefToReal[{X1,X2},RotateLeft[Triangle,1]],RefQuad[{X1,X2},Midpt]],
Append[RefToReal[{X1,X2},RotateLeft[Triangle,2]],RefQuad[{X1,X2},Midpt]],
Append[RefToReal[{X1,X2},RotateLeft[Triangle,0]],RefQuad[{X1,X2},Corner]],
Append[RefToReal[{X1,X2},RotateLeft[Triangle,1]],RefQuad[{X1,X2},Corner]],
Append[RefToReal[{X1,X2},RotateLeft[Triangle,2]],RefQuad[{X1,X2},Corner]],
Append[RefToReal[{X1,X2},Triangle],0]
},
{X1,0,1},{X2,0,1},
RegionFunction->Function[{x1,x2,u,X1,X2},(X1+X2<1)],
BoxRatios->{1,1,0.5},
BoundaryStyle->Directive[Thick,Green]
]

Piecewise Cubics

Looking at the thrid entry of Table 3.7 in terms of the reference triangle there are several flavors of cubics.

P3Dots=Flatten[
Table[
Table[{i/3,j/3},{i,0,3-j}],
{j,0,3}],
1
]
Show[
Graphics[{
{LightBlue,Polygon[{{0,0},{0,1},{1,0}}]},
{Red, PointSize[0.02], Table[Point[P3Dots[[i]]],{i, Length[P3Dots]}]},
{Table[Text[i,P3Dots[[i]],{-1.5,-1.5}],{i, Length[P3Dots]}]}
}
],
Frame->True
]

I think I should automate some of the process of generating the reference triangle cubics.  Here is a basis (dimension 10) for the space of cubics alonmg with the values of those cubics on the dots defined up above.

fs[{X1_,X2_}] :=  {1, X1, X2, (X1^2),(X2^2),X1 X2, (X1^3),(X1^2)X2, X1 (X2^2), (X2^3)}
A=Map[fs,P3Dots];
MatrixForm[A]

A general cubic is a linear combination of the basis I wrote down.  The matrix A is just the values of my basis at the points I care about.  What I want is the linear combinations that give simple values at the dots.  This is easy if I think about it abstractly.  Here are three of the special functions.

f1[{X1_,X2_}]=LinearSolve[A, {1,0,0,0,0,0,0,0,0,0}].fs[{X1,X2}]
f2[{X1_,X2_}]=LinearSolve[A, {0,1,0,0,0,0,0,0,0,0}].fs[{X1,X2}]
f3[{X1_,X2_}]=LinearSolve[A, {0,0,1,0,0,0,0,0,0,0}].fs[{X1,X2}]
Plot3D[{
f1[{X1,X2}],
f2[{X1,X2}],
f3[{X1,X2}]
},
{X1,0,1},
{X2,0,1},
RegionFunction->Function[{X1,X2,z},X1+X2<1],
PlotRange->All,
BoxRatios->{1,1,0.5},
BoundaryStyle->Directive[Thick,Green]
]

Here are all of them.

fCubicNodes[{X1_,X2_}]=Transpose[Inverse[A]].fs[{X1,X2}];
Plot3D[
fCubicNodes[{X1,X2}],
{X1,0,1},
{X2,0,1},
RegionFunction->Function[{X1,X2,z},X1+X2<1],
PlotRange->All,
BoxRatios->{1,1,0.5},
BoundaryStyle->Directive[Thick,Green]
]

Here is a picture of all 10 of them on a general triangle.

Triangle= RandomReal[{-1,1},{3,2}];
Temp3DPlot[{X1_,X2_}]= Table[
Append[RefToReal[{X1,X2},Triangle],fCubicNodes[{X1,X2}][[i]]],
{i, 1, 10}];
ParametricPlot3D[Temp3DPlot[{X1,X2}],
{X1,0,1},{X2,0,1},
RegionFunction->Function[{x1,x2,u,X1,X2},(X1+X2<1)],
BoxRatios->{1,1,0.5},
BoundaryStyle->Directive[Thick,Green]
]

Every so often you get a skinny triangle lined up with the diagonal.  It looks cute :)

Of course, a bunch of these are symmetry related!

More Piecewise Cubics

Looking at the fourth entry of Table 3.7 in terms of the reference triangle there are several flavors of cubics.  I now need second derivatives along the sides (this excludes the mixed derivative at the origin etc.) and function values.  The comp

Clear[fs,X1,X2,D2fs,P3]
P3={{0,0},{1,0},{0,1}};
fs[{X1_,X2_}] :=  {1, X1, X2, (X1^2),(X2^2),X1 X2, (X1^3),(X1^2)X2, X1 (X2^2), (X2^3)};
D2fs[{X1_,X2_}]= D[fs[{X1,X2}],{{X1,X2},2}];
A=Join[
Map[fs,P3],
{fs[{1/3,1/3}],
D2fs[{0,0}].{1,0}.{1,0},
D2fs[{0,0}].{0,1}.{0,1},
D2fs[{1,0}].{1,0}.{1,0},
D2fs[{1,0}].{(1/(Sqrt[2])),-(1/(Sqrt[2]))}.{(1/(Sqrt[2])),-(1/(Sqrt[2]))},
D2fs[{0,1}].{0,1}.{0,1},
D2fs[{0,1}].{(1/(Sqrt[2])),-(1/(Sqrt[2]))}.{(1/(Sqrt[2])),-(1/(Sqrt[2]))}
}];
MatrixForm[A]

I think I should automate some of the process of generating the reference triangle cubics.  Here are the basis functions

fCubicNodes[{X1_,X2_}]=Transpose[Inverse[A]].fs[{X1,X2}];
Plot3D[
fCubicNodes[{X1,X2}],
{X1,0,1},
{X2,0,1},
RegionFunction->Function[{X1,X2,z},X1+X2<1],
PlotRange->All,
BoxRatios->{1,1,0.5},
BoundaryStyle->Directive[Thick,Green]
]

Here is a picture of all 10 of them on a general triangle.

Triangle= RandomReal[{-1,1},{3,2}];
Temp3DPlot[{X1_,X2_}]= Table[
Append[RefToReal[{X1,X2},Triangle],fCubicNodes[{X1,X2}][[i]]],
{i, 1, 10}];
ParametricPlot3D[Temp3DPlot[{X1,X2}],
{X1,0,1},{X2,0,1},
RegionFunction->Function[{x1,x2,u,X1,X2},(X1+X2<1)],
PlotRange->All,
BoxRatios->{1,1,0.5},
BoundaryStyle->Directive[Thick,Green]
]

When you look at this you see 4 big functions and know there must be 6 little ones.  Ths big ones are the function values ones and are the first 4 because of the order we wrote the equations down in.  Here is a plot of the small six all on their own.

Triangle= RandomReal[{-1,1},{3,2}];
Temp3DPlot[{X1_,X2_}]= Table[
Append[RefToReal[{X1,X2},Triangle],fCubicNodes[{X1,X2}][[i]]],
{i, 5, 10}];
ParametricPlot3D[Temp3DPlot[{X1,X2}],
{X1,0,1},{X2,0,1},
RegionFunction->Function[{x1,x2,u,X1,X2},(X1+X2<1)],
PlotRange->All,
BoxRatios->{1,1,0.5},
BoundaryStyle->Directive[Thick,Green]
]

2D: Long Skinny Triangles

The idea is to use a coarser approximation on one side (presumably the short one) of a long skinny triangle.  I need to work out what functions I can use. I suspect things need to be degenerate it is clear that I need to be careful.

On my ref triangle the following set 
	{1, X1, X2, X1 (1 - X1 - X2),X2(1-X1-X2)}
is a basis for a five dimensional subset of the quadratics.

PStarDots= {{0,0}, {1,0},{0,1},{1/2,0},{0,1/2}};
Show[
Graphics[{
{LightBlue,Polygon[{{0,0},{0,1},{1,0}}]},
{Red, PointSize[0.02], Table[Point[PStarDots[[i]]],{i, Length[PStarDots]}]},
{Table[Text[i,PStarDots[[i]],{-1.5,-1.5}],{i, Length[PStarDots]}]}
}
],
Frame->True
]

fs[{X1_,X2_}] :=  {1, X1, X2, X1 (1- X1 -X2), X2 (1- X1 -X2)}
A=Map[fs,PStarDots];
MatrixForm[A]
Det[A]

fStarQuadraticNodes[{X1_,X2_}]=Transpose[Inverse[A]].fs[{X1,X2}];
Plot3D[
fStarQuadraticNodes[{X1,X2}],
{X1,0,1},
{X2,0,1},
RegionFunction->Function[{X1,X2,z},X1+X2<1],
PlotRange->All,
BoxRatios->{1,1,0.5},
BoundaryStyle->Directive[Thick,Green]
]

2D: More Long Skinny Triangles

The idea is to use a longer approximation on one side (presumably the longets one) of a long skinny triangle.  I need to work out what functions I can use. I suspect things need to be degenerate it is clear that I need to be careful.

On my ref triangle the following set 
	{1, X1, X2,  X1 X2}
is a basis for a four dimensional subspace of the quadratics.  It works with the obvious points

PStarDots= {{0,0}, {1,0},{0,1},{1/2,1/2}};
Show[
Graphics[{
{LightBlue,Polygon[{{0,0},{0,1},{1,0}}]},
{Red, PointSize[0.02], Table[Point[PStarDots[[i]]],{i, Length[PStarDots]}]},
{Table[Text[i,PStarDots[[i]],{-1.5,-1.5}],{i, Length[PStarDots]}]}
}
],
Frame->True
]

fs[{X1_,X2_}] :=  {1, X1, X2, X1 X2}
A=Map[fs,PStarDots];
MatrixForm[A]
Det[A]

fStarQuadraticNodes[{X1_,X2_}]=Transpose[Inverse[A]].fs[{X1,X2}];
Plot3D[
fStarQuadraticNodes[{X1,X2}],
{X1,0,1},
{X2,0,1},
RegionFunction->Function[{X1,X2,z},X1+X2<1],
PlotRange->All,
BoxRatios->{1,1,0.5},
BoundaryStyle->Directive[Thick,Green]
]