21
$\begingroup$

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["NDSolve`FEM`*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!

Improve this question
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – MarcoB
    Apr 18, 2016 at 20:25
  • 5
    $\begingroup$ Would there be interest in a "FEM theory" tutorial? $\endgroup$
    – user21
    Apr 19, 2016 at 1:12
  • $\begingroup$ @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. $\endgroup$
    – AxelF
    Apr 19, 2016 at 6:06
  • $\begingroup$ 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. $\endgroup$
    – user21
    Apr 19, 2016 at 9:31
  • 5
    $\begingroup$ @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. $\endgroup$
    – J. M.'s persistent exhaustion ♦
    Apr 21, 2016 at 5:21

2 Answers 2

Reset to default
23
$\begingroup$

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

Load the package:

Needs["NDSolve`FEM`"]

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"]

enter image description here

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"]]

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – Андрей Кротких
    Nov 19, 2018 at 0:02
  • 2
    $\begingroup$ It seems there is a missing definition coords = mesh["Coordinates"]. (Added.) $\endgroup$
    – Anton Antonov
    Dec 28, 2022 at 2:24
  • $\begingroup$ @AntonAntonov, thanks. $\endgroup$
    – user21
    Dec 28, 2022 at 6:55
1
$\begingroup$

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

Needs["NDSolve`FEM`"]
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]

enter image description here

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

Hope it's what you're looking for?

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.