# How to access FEM shape functions?

A couple of days ago I asked here about surface meshes and plotting on surfaces.

Now I have another question: How can I access the surface or boundary element shape functions?

I would like to approximate the scalar stream function s on my surface by finite elements, where I specify values of s on the nodes and approximate it with shape functions. That would allow my to calculate the current density as the curl of the normal vector multiplied with s, and from there I can calculate magnetic fields, inductance, eddy currents on other surfaces and so on.

Names["NDSolveFEM*Shape*"]
(*
{"ElementShapeFunction","ElementShapeFunctionDerivative",
"FEMShapeFunctionTest","FindShapeFunction",
"GetIntegratedShapeFunction","GetIntegratedShapeFunctionDerivative",
"IntegratedShapeFunction","IntegratedShapeFunctionDerivative"}
*)


seems to indicate that shape functions exist in the FEM universe, but couldn't find any documentation.

I know that I am intending to use the FEM package in an unconventional way, but if someone could help me here, that would be great!

• Have you tried asking Wolfram directly? These functions are not documented, and perhaps somebody here has spent considerable time testing them out via trial and error (which is not impossible!), but you are nevertheless more likely to obtain that kind of information directly from the horse's mouth. Apr 18, 2016 at 20:25
• Would there be interest in a "FEM theory" tutorial? Apr 19, 2016 at 1:12
• @user21 Absolutely! As I understand it, conceptually most FEM packages work from bottom to top, i.e. from elements to the solution, whereas Mathematica works top to bottom, where you start with your PDE and drill down to the finite elements. Apr 19, 2016 at 6:06
• That's a nice way of putting it. For the tutorial, it quite a bit of work to write on but I'll keep it in mind. Apr 19, 2016 at 9:31
• @user21, I understand that it might take an entire book to explain everything nicely, but I also want to assure you that such an effort will be appreciated by a lot of people dealing with FEM. Apr 21, 2016 at 5:21

There are no surface element shape functions. There are, however, the normal shape functions.

Needs["NDSolveFEM"]


This gives you the shape functions for implemented elements (see documentation)

elementOrder = 1;
ElementShapeFunction[TriangleElement, elementOrder][r, s]
{1 - r - s, r, s}

ElementShapeFunction[TriangleElement, 2][r, s]
{1 + 2 r^2 - 3 s + 2 s^2 + r (-3 + 4 s), r (-1 + 2 r),
s (-1 + 2 s), -4 r (-1 + r + s), 4 r s, -4 s (-1 + r + s)}


This gives you the derivative of the shape function:

ElementShapeFunctionDerivative[TriangleElement, elementOrder][r, s]
{{-1, 1, 0}, {-1, 0, 1}}


This gives you the integrated shape function:

integrationOrder = 2;
IntegratedShapeFunction[TriangleElement, elementOrder, \
integrationOrder]
{{{0.6666666666666667, 0.16666666666666666,
0.16666666666666666}}, {{0.1666666666666667, 0.6666666666666666,
0.16666666666666666}}, {{0.16666666666666674,
0.16666666666666666, 0.6666666666666666}}}


These are the integration points and weights:

ElementIntegrationPoints[TriangleElement, integrationOrder]
{{0.16666666666666666, 0.16666666666666666}, {0.6666666666666666,
0.16666666666666666}, {0.16666666666666666, 0.6666666666666666}}

ElementIntegrationWeights[TriangleElement, integrationOrder]
{0.16666666666666666, 0.16666666666666666, 0.16666666666666666}


New shape functions can be found with FindShapeFunction. For that we need the base polynomial they should use and the base coordinates of the element.

MeshElementBasePolynomial[TriangleElement, elementOrder, {r, s}]
{1, r, s}
MeshElementBasePolynomial[TriangleElement, 2, {r, s}]
{1, r, s, r^2, r s, s^2}
MeshElementBaseCoordinates[TriangleElement, elementOrder]
{{0, 0}, {1, 0}, {0, 1}}

FindShapeFunction[
MeshElementBasePolynomial[TriangleElement, elementOrder, {r, s}],
MeshElementBaseCoordinates[TriangleElement, elementOrder], {r, s}]
{1 - r - s, r, s}


So if we want to find the shape function of a nine node quad element we use the 2nd order quad element add a coordinate and a term to the polynomial:

qp = Join[MeshElementBasePolynomial[QuadElement, 2, {r, s}], {r^2*s^2}]
{1, r, s, r s, r^2, r^2 s, r s^2, s^2, r^2 s^2}
qc = Join[MeshElementBaseCoordinates[QuadElement, 2], {{0, 0}}]
{{-1, -1}, {1, -1}, {1, 1}, {-1, 1}, {0, -1}, {1, 0}, {0, 1}, {-1,
0}, {0, 0}}
sf = FindShapeFunction[qp, qc, {r, s}]
{1/4 (-1 + r) r (-1 + s) s, 1/4 r (1 + r) (-1 + s) s,
1/4 r (1 + r) s (1 + s),
1/4 (-1 + r) r s (1 + s), -(1/2) (-1 + r^2) (-1 + s) s, -(1/2)
r (1 + r) (-1 + s^2), -(1/2) (-1 + r^2) s (1 + s), -(1/
2) (-1 + r) r (-1 + s^2), (-1 + r^2) (-1 + s^2)}


This new shape function evaluates to 1 at the nodes and the sum is 1:

Function[{r, s}, Evaluate[sf]] @@@ qc
{{1, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1,
0, 0, 0, 0, 0, 0}, {0, 0, 0, 1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0,
0, 0, 0}, {0, 0, 0, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 1, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 1}}
Total[sf] // Simplify
1


Update:

Let's look at how the shape functions are used to map the integration points into each element in global space.

Make a mesh:

order = 1;
mesh = ToElementMesh[Disk[], "MeshOrder" -> order];
mesh["Wireframe"]


Get the integrated shape function:

intOrder = 4;
isf = IntegratedShapeFunction[TriangleElement, order, intOrder];


As a side note, the integrated shape function is the same as the shape function evaluated at the integration points of the mother element:

sf = ElementShapeFunction[TriangleElement, order];
sf[r, s]
ip = ElementIntegrationPoints[TriangleElement, intOrder];
isf[[All, 1]] === (sf @@@ ip)


Get the element coordinates:

coords = mesh["Coordinates"];
eleCoords =
GetElementCoordinates[coords,
ElementIncidents[mesh["MeshElements"]][[1]]];


Map the integrated shape functions into the elements:

mappedCoords = (isf[[All, 1]] . #) & /@ eleCoords;


These are the coordinates at which the PDE coefficients are evaluated, for a shape function of order and an integration order of intOrder.

Visualize the mapped integration points:

Show[Graphics[{Red, Point /@ mappedCoords}], mesh["Wireframe"]]


• But how we can transform shape functions in base coordinate to real coordinates? My experiments with using Jacobi matrix didn't give me right results. Nov 19, 2018 at 0:02
• It seems there is a missing definition coords = mesh["Coordinates"]. (Added.) Dec 28, 2022 at 2:24
• @AntonAntonov, thanks. Dec 28, 2022 at 6:55

For a given mesh you get the nodal shape functions as follows:

Needs["NDSolveFEM"]
mesh = ToElementMesh[Disk[] ];

pi=mesh["Coordinates"];


shape functions of single nodes:

nodefun = Map[ElementMeshInterpolation[mesh, #] &,IdentityMatrix[Length[pi]]];


examplary for node #75

Plot3D[nodefun[[75]][x, y], Element[{x, y}, mesh],Mesh->False]


Through[nodefun[x, y]] gives a list of all nodal shapefunctions.

Hope it's what you're looking for?