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]).

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: Quads 
Reference set 0≤(Subscript[X, 1])≤1 and 0≤(Subscript[X, 2])≤1

Our basis functions are products of cubics in X1 and X2

Clear[X1,X2,Basis]
Basis[{X1_,X2_}]=Flatten[
Table[
(X1^i)(((1-X1))^(3-i)) (X2^j)(((1-X2))^(3-j)),
{i,0,3},
{j,0,3}]
];
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,1},{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,2]]
];
Det[A]

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

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}][[{5,9}]];
ShapeFunction[{X1_,X2_},T]=ShapeFunctions[{X1,X2}][[13]];
Plot3D[
ShapeFunction[{X1,X2},V],
{X1,0,1},
{X2,0,1},
AxesLabel->{"X1","X2","Z"}
]
Plot3D[
ShapeFunction[{X1,X2},G],
{X1,0,1},
{X2,0,1},
AxesLabel->{"X1","X2","Z"}]
Plot3D[
ShapeFunction[{X1,X2},T],
{X1,0,1},
{X2,0,1},
AxesLabel->{"X1","X2","Z"}
]