16
$\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$
  • $\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 '16 at 20:25
  • 3
    $\begingroup$ Would there be interest in a "FEM theory" tutorial? $\endgroup$ – user21 Apr 19 '16 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 '16 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 '16 at 9:31
  • 2
    $\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 discontentment♦ Apr 21 '16 at 5:21
16
$\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
| improve this answer | |
$\endgroup$
  • $\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 '18 at 0:02

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.