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.

C1 Element: 2D Triangle 
General Quintics
Reference set 0≤(Subscript[X, 1])≤1 and 0≤(Subscript[X, 2])≤1-(Subscript[X, 1]).

Basic Scalar Elements

For convenience we will use the Bernstein-Bezier representation.  For the piecewise linear this means 
	{1,X1,X2}⟷{X1,X2,1-X1-X2}:={(Subscript[λ, 1]),(Subscript[λ, 2]),(Subscript[λ, 3])}
For Quadratics this means
	{1,X1,X2,(X1^2),X1 X2,(X2^2)}⟷(Subscript[{((Subscript[λ, 1])^i)((Subscript[λ, 2])^j)((Subscript[λ, 3])^(2-i-j))}, i,j=0… 2])
For Quintics this means
	(Subscript[{((Subscript[λ, 1])^i)((Subscript[λ, 2])^j)((Subscript[λ, 3])^(5-i-j))}, i=0… 5,j=0…5-i])
The advantages of this are that the formulas are that we will get symmetric formulas for the three corners.

Clear[X1,X2,Basis]
Basis[{X1_,X2_}]= Flatten[
Table[
(X1^i)(X2^j)(((1-X1-X2))^(5-i-j)),
{i,0,5},
{j,0,5-i}]
]
DBasis[{X1_,X2_}]=D[Basis[{X1,X2}],{{X1,X2}}];
DDBasis[{X1_,X2_}]=D[DBasis[{X1,X2}],{{X1,X2}}];

The computation is basically the same as before.  the only difference is that now at the vertices we need derivatives and second derivatives while at the midpoints we need normal derivatives.

Pts={{0,0},{0,1},{1,0}};
MidPts={{1/2,0},{1/2,1/2},{0,1/2}};
Ns =  {{0,-1},{1/(Sqrt[2]),1/(Sqrt[2])},{-1,0}};
A=Join[
Map[ Basis, Pts],
Map[DBasis,Pts][[All,All,1]],
Map[DBasis,Pts][[All,All,2]],
Map[DDBasis,Pts][[All,All,1,1]],
Map[DDBasis,Pts][[All,All,2,2]],
Map[DDBasis,Pts][[All,All,1,2]],
MapThread[DBasis[#1].#2&,{MidPts, Ns}]
];
Det[A]

Because of the order of the equations: 
A) The basis functions 1,2,3 are going to be the "value" functions at the corners.  
B) The basis functions 4,5,6 are the X1 slope functions at the corners.
C) The basis functions 7,8,9 are the X2 slope functions at the corners.
D) The basis functions 10,11,12 are the X1-X1 second derivative functions at the corners. 
E) The basis functions 13,14,15 are the X2-X2 second derivative functions at the corners.
D) The basis functions 16,17,18 are the X1-X2 second derivative functions at the corners.
F) The basis functions 19,20,21 are the normal derivative functions at the MidPts.

ShapeFunctions[{X1_,X2_}]=FullSimplify[Transpose[Inverse[A]].Basis[{X1,X2}]];
TableForm[ShapeFunctions[{X1,X2}],TableHeadings->{Automatic}]

As always we only need one representative of each flavor.  I am going to tag the functions as V, G, H, and M.

ShapeFunction[{X1_,X2_},"V"]=ShapeFunctions[{X1,X2}][[1]];
ShapeFunction[{X1_,X2_},"G"]=ShapeFunctions[{X1,X2}][[{4,7}]];
ShapeFunction[{X1_,X2_},"H"]=ShapeFunctions[{X1,X2}][[{10,13,16}]];
ShapeFunction[{X1_,X2_},"M"]=ShapeFunctions[{X1,X2}][[19]];
Plot3D[
ShapeFunction[{X1,X2},"V"],
{X1,0,1},
{X2,0,1},
RegionFunction->Function[{X1,X2,Z}, X1+X2<1],
AxesLabel->{"X1","X2","Z"}
]
Plot3D[
ShapeFunction[{X1,X2},"G"],
{X1,0,1},
{X2,0,1},
RegionFunction->Function[{X1,X2,Z}, X1+X2<1],
AxesLabel->{"X1","X2","Z"}]
Plot3D[
ShapeFunction[{X1,X2},"H"],
{X1,0,1},
{X2,0,1},
RegionFunction->Function[{X1,X2,Z}, X1+X2<1],
AxesLabel->{"X1","X2","Z"}
]
Plot3D[
ShapeFunction[{X1,X2},"M"],
{X1,0,1},
{X2,0,1},
RegionFunction->Function[{X1,X2,Z}, X1+X2<1],
AxesLabel->{"X1","X2","Z"}
]

C1 Element: 2D Triangle 
General Quintics
General Triangle

Now we need to move these to a general triangle from the reference map. Our mapping is
	(OverVector[x])= (OverVector[b])+H (OverVector[X])  or equivalently (OverVector[X])= (H^(-1)) ((OverVector[x])-(OverVector[b]))
and the corners of our general triangle are 
	(Subscript[(OverVector[p]), 1])=(OverVector[b])+H (Subscript[(OverVector[e]), 1])  and (Subscript[(OverVector[p]), 2])=(OverVector[b])+H (Subscript[(OverVector[e]), 2]) and (Subscript[(OverVector[p]), 3])=(OverVector[b])+{0,0}.
This means that H=[(Subscript[(OverVector[p]), 1])⋮ (Subscript[(OverVector[p]), 2])] just like before.

Value Functions

Our translation for the value functions is straightforward. All we need to do is stuff (H^(-1))(OverVector[x])  into our functions and we are in business.  In other words 
	ϕ[(OverVector[x]),V]=Φ[(H^(-1))(OverVector[x]),V]

p1={1,0.9};
p2={0,1.2};
H= Transpose[{p1,p2}];
HInv=Inverse[H];
ParametricPlot3D[{
Join[H.{X1,X2},  {0}],
Join[H.{X1,X2},  {ShapeFunction[{X1,X2},"V"]}]
},
{X1,0,1},
{X2, 0, 1-X1 },
AspectRatio->1
]

You can see that the derivatives are going to have to scale.

Gradient Function

Our translation for the gradient functions is less-straightforward because of the chain rule.  It is clear that our gradients in the current configuration are a linear combination of the gradients in the current configuration mapped into the current configuration.  In other words, 
	ϕ[(OverVector[x]),G]=M Φ[(H^(-1))(OverVector[x]),G]
where M is a 2×2 matrix that we need to determine to make the gradient functions nodal. Differentiating and using the chain rule  gives
	(Subscript[∇, (OverVector[x])])ϕ[(OverVector[x]),G]=M(Subscript[∇, (OverVector[X])])Φ[(H^(-1))(OverVector[x]),G](H^(-1))
now plugging in our special point (which for our explicit formulas is {0,0} in the reference and (OverVector[b]) in the real) gives the identity for both the gradients:  the reference version because that is the way we made it and the current version because that is what we want.
	I=M I (H^(-1))
We clearly need M=H.

Note that all the other nodal values are zero automatically.

p1={0,1};
p2={1,2.3};
H= Transpose[{p1,p2}];
HInv=Inverse[H];
ParametricPlot3D[{
Join[H.{X1,X2},  {0}],
Join[H.{X1,X2},  {(ShapeFunction[{X1,X2},"G"].H)[[1]]}],
Join[H.{X1,X2},  {(ShapeFunction[{X1,X2},"G"].H)[[2]]}]
},
{X1,0,1},
{X2, 0, 1-X1 },
AspectRatio->1
]

Hessian Functions

Our translation for the Hessian functions is similar.  It is clear that our Hessian entries (three of them) in the current configuration are a linear combination of the functions in the reference configuration mapped into the current configuration.  In other words, 
	ϕ[(OverVector[x]),H]=M ∇Φ[(H^(-1))(OverVector[x]),H]
where M is a 3×3 matrix that we need to determine to make the Hessian functions nodal.  Differentiating and using the chain rule  gives
	(Subscript[∇, (OverVector[x]) (OverVector[x])])ϕ[(OverVector[x]),H]=M(Subscript[∇, (OverVector[X]) (OverVector[X])])Φ[(H^(-1))(OverVector[x]),H](H^(-1))(H^(-1))=M (Subscript[∇, (OverVector[X]) (OverVector[X])])Φ[(H^(-1))(OverVector[x]),H](H^(-2))
where (Subscript[∇, (OverVector[x]) (OverVector[x])])ϕ[(OverVector[x]),H] is 3×2×2 now plugging in our special point (which for our explicit formulas is {0,0} in the reference and (OverVector[b]) in the real) gives the 3×3 tensor 
	(({{(({{1, 0}, {0, 0}}))}, {(({{0, 0}, {0, 1}}))}, {(({{0, 1}, {1, 0}}))}}))
for both the Hessians:  the reference version because that is the way we made it and the current version because that is what we want.  It is clear that we need

Hmm=D[ ShapeFunction[{X1,X2},"H"],{{X1,X2},2}]/.{X1->0,X2->0}

p1={1,0.9};
p2={0,1.2};
H= Transpose[{p1,p2}];
HInv=Inverse[H];
ParametricPlot3D[{
Join[H.{X1,X2},  {0}],
Join[H.{X1,X2},  {(H.H.ShapeFunction[{X1,X2},"H"])[[1]]}],
Join[H.{X1,X2},  {(ShapeFunction[{X1,X2},"H"].H.H)[[2]]}]
},
{X1,0,1},
{X2, 0, 1-X1 },
AspectRatio->1
]