Fundamentals of elements.

[gatling-nrl]

I'm trying to write down some "elements for beginners" type notes. Are these statements correct?

---

Note: Tri can be changed to Quad when using quadratic meshes instead of triangular meshes.

ElementTriP0 Piecewise constant over the mesh cell. Discontinuous (having two values) on the boundary between cells.

ElementTriP1 Piecewise linear over the mesh cell. Continuity enforced at boundary. Derivative is discontinuous at boundary.
* For beginners, start here. This element can solve a huge space of problems, especially with refined meshes.

ElementTriPx Piecewise quadratic, etc. over the mesh cell. Continuity enforced at boundary. Derivative is also continuous.
* Only use when solution is quite smooth, because of the enforcement of continuous derivate(s).

ElementTriRT0 Lowest order Raviart-Thomas element. Divergence of a function in ElementTriRT0 will be in ElementTriP0.

ElementDG(element) will remove the continuity constraints at the boundary. Note that ElementDG(ElementTriP0()) == ElementTriP0()

---

Edit: should I expand the Element class docstrings with this info so it ends up in https://scikit-fem.readthedocs.io/en/latest/api.html#module-skfem.element?

Edit 2: I assume .. autosummary:: will only include the first line of the docstring... perhaps the element names should be links to respective element pages with notes, caveats, and an example usage of that element?

[gatling-nrl]

What is meant by "cutting connectivity" in the description of ElementDG?

scikit-fem/skfem/element/element_dg.py

Lines 4 to 5 in c25c9a8

	class ElementDG(Element):
	"""Turn a finite element discontinuous by cutting connectivity."""

[kinnala]

Try plotting:

from skfem import *
from skfem.visuals.matplotlib import *
import numpy as np
m = MeshTri().refined(2)
basis1 = Basis(m, ElementTriP1())
basis2 = Basis(m, ElementDG(ElementTriP1()))
plot(basis1, np.random.rand(basis1.N), shading='gouraud')
plot(basis2, np.random.rand(basis2.N), shading='gouraud')

[kinnala]

"cutting connectivity" between neighboring elements.

[kinnala]

Not sure if it is the best terminology, maybe "by turning all DOFs into interior DOFs" would be better.

[gatling-nrl]

"by turning all DOFs into interior DOFs"

By interior DOFs, you mean interior to the element, and are saying turn them into interior DOFs because the interior DOFs already had no continuity constraints on them, correct?

So... the "connectivity" is the notion that triangle A is "connected" to triangle B by facet 1 and that A and B should have the same value on either side of facet 1? This notion enters in through assembly, where it controls a row in the stiffness matrix, effectively restricting the space of possible solutions? And by "cutting connectivity" you remove this step during assembly?

[kinnala]

ElementTriPx Piecewise quadratic, etc. over the mesh cell. Continuity enforced at boundary. Derivative is also continuous. * Only use when solution is quite smooth, because of the enforcement of continuous derivate(s).

Derivative is not continuous. These are all "H^1-conforming elements" meaning that only the value is continuous. Elements that have continuous derivatives are called "H^2-conforming elements", e.g., ElementTriArgyris.

[gatling-nrl]

Why are higher order elements only advisable when the solution is smooth?

[kinnala]

Smoothness roughly means that are there are no kinks in the solution. High-order polynomials do not describe well such kinks and start to oscillate when using them to approximate kinks.

Consider this function:

I try fitting 50 degree polynomial:

Looks pretty OK, but FEM actually minimizes error in the derivative so let's compare the derivatives:

Those oscillations happen because of the kink.

This was done using numpy polyfit which has some regularization trying to avoid the oscillations. However, typically in FEM you do not have such things and it just won't work too well. Actually, the error is bounded by something which depends on the size of the element to the power of the polynomial order minus some kind of smoothness index. Specifying this exactly would need some Sobolev space theory.

[gatling-nrl]

The oscillations remind me of a Fourier series trying to represent such step function. I see what you mean, about the high order polynomials having a lot of extra freedom they don't need.

This was done using numpy polyfit which has some regularization

I think I saw this play out when I reproduced your plots and numpy complained with a warning

RankWarning: Polyfit may be poorly conditioned

Low order polynomials do a horrible job too though... the blue line is a quadratic.

Is your point that in FEM this could (should?) be solved by refining the mesh so that it can use different quadratics (or even linear functions) in the different regions, rather than increasing the order of the element?

[kinnala]

Yes. There is a theory about "hp-FEM" which tries automating where you should split the element and where you should increase the polynomial order.

[kinnala]

Edit: should I expand the Element class docstrings with this info so it ends up in https://scikit-fem.readthedocs.io/en/latest/api.html#module-skfem.element?

Sure! Makes sense. I think that is from skfem.element.__init__ docstring.

Edit 2: I assume .. autosummary:: will only include the first line of the docstring... perhaps the element names should be links to respective element pages with notes, caveats, and an example usage of that element?

Possibly yes but maintaining such a list can be tedious. Any ideas welcome.

[gdmcbain]

Note: Tri can be changed to Quad when using quadratic meshes instead of triangular meshes.

Almost, except that the P is dropped from the Lagrange elements, e.g. ElementTriP1 → ElementQuad1. (I have made this mistake in both directions many times using automatic find-and-replace…)

[gatling-nrl]

Is there a reason? I though the P meant "polynomial"... Is ElementQuadP1=ElementQuad1 worth adding?

Edit: So the search/replace would be TriP <-> Quad ?

[gdmcbain]

Edit: So the search/replace would be TriP <-> Quad ?

Except for ElementTriRT0 and ElementQuadRT0. The other non-P elements aren't shared by Tri and Quad.

Yes, I think ℙ is for polynomial, but often the polynomial quadrilateral (or more generally ‘tensor product’, so hexahedra in three dimensions) elements are written with a ℚ; e.g. in Ern & Guermond (2004, §1.2.4) and Arnold's Periodic Table. Also the ℙ is only used for simplicial elements (line segment, triangle, tetrahedron). There is a technical reason for not identifying ℚ with ℙ: whereas ℙ contains all the polynomials of a certain degree, ℚ is bigger because of the tensor product; e.g. ℚ₁ contains xy but not x² or y².