getting started with dofs, quadrature points, and etc....

[gatling-nrl]

(I've been following this FEM tutorial and this extremely thorough pdf to get myself up to speed on both FEM in general and how it is implemented in skfem. They are quite handy but the nomenclature does sometimes differ from sfkem and the entire FEniCS code base in currently in shambles while they refactor with an intense focus on parallelization. If you're new here, in comparison to FEniCS and SfePy, skfem is extremely lightweight and the source code is (relatively) approachable.)

I feel like skfem is lacking in some of the fundamental documentation, especially for newbies. So I'm working on that and @kinnala and @gdmcbain please correct anything I've got wrong below!

First up is DoFs vs quadrature points... and clearly they are different.

import gmsh
import matplotlib.pyplot as plt
import meshio
import numpy as np

import skfem
from skfem.helpers import dot, grad
import skfem.io

fig, ax = plt.subplots(1,2,figsize=(12,6))
ax[0].plot(mesh.p[0], mesh.p[1], 'ok', markersize=15)

for t in mesh.t.T:
    t = [t[0], t[1], t[2], t[0]]
    ax[0].plot(mesh.p[0][t], mesh.p[1][t], 'k')

basis = skfem.CellBasis(mesh, skfem.ElementTriP3(), intorder=6)

for d in basis.doflocs.T[basis.nodal_dofs.reshape(-1)]:
    ax[0].plot(d[0], d[1], 'or')
for d in basis.doflocs.T[basis.facet_dofs.reshape(-1)]:
    ax[0].plot(d[0], d[1], 'og')
for d in basis.doflocs.T[basis.interior_dofs.reshape(-1)]:
    ax[0].plot(d[0], d[1], 'ob')

for x in basis.X.T:
    plt.plot([0,0,1,0],[0,1,0,0], 'k')  # draw a reference triangle to make it easier to understand
    ax[1].plot(x[0], x[1], 'or')

So lots of things to digest here... first, the quadrature points are in local coordinates for one "reference" cell, as seen in the plot on the right. Next up is the various categories of DoF in skfem as drawn on the left, which are clearly in 'global' or 'physical' coordinates:

• basis.nodal_dofs in red
• basis.facet_dofs in green
• basis.interior_dofs in blue

And then there's

• basis.element_dofs which is all of the above

1d and 3d meshes use slightly different labeling as shown below. Any DoF not fitting one of the labelled locations is an interior_dof.

Now to recap so far, there is an element order to consider, the * in skfem.ElementTriP*(), which changes the DoFs per element. P1 is special because the doflocs are right on the mesh vertexes. All of the others (including P0) move the DoF locations around. There is also an integration order to consider (the intorder argument to skfem.CellBasis) that controls the number of quadrature points. (As best I can tell, skfem exclusively uses equispaced Lagrange quadrature.)

We can use what we have to far to compute fun integrals...

a_domain = 2  # domain
a_tool = 1  # tool
gmsh.initialize()
try:
    domain_tag = gmsh.model.occ.add_rectangle(x=-a_domain, y=-a_domain, z=0, dx=2*a_domain, dy=2*a_domain)
    tool_tag = gmsh.model.occ.add_rectangle(x=-a_tool, y=-a_tool, z=0, dx=2*a_tool, dy=2*a_tool)
    gmsh.model.occ.cut([(2,domain_tag)], [(2,tool_tag)], removeTool=True)
    gmsh.model.occ.synchronize()
    gmsh.model.mesh.generate(2)
    gmsh.write('tmp.msh')
finally:
    gmsh.clear()
    gmsh.finalize()
mesh = skfem.io.from_meshio(meshio.read('tmp.msh')).refined(3)
basis = skfem.CellBasis(mesh, skfem.ElementTriP4(), intorder=19)
@skfem.Functional
def test1(w):
    x, y = w['x'].value
    return x**2 + y**2
print(skfem.asm(test1, basis))
sigma = basis.zeros()
sigma[:] = .5
@skfem.Functional
def test2(w):
    sigma = w['sigma']
    x, y = w['x'].value
    return sigma*x**2 + y**2
print(skfem.asm(test2, basis, sigma=sigma))
mesh

But here is where I start to get really lost.

(1) What if I have a functional form of sigma, e.g. sigma = np.sin(x) + np.exp(-(y**2))?
(2) I used sigma=basis.zeros() making sigma shaped like DoFS, but then i didn't need .interpolate() to get it interpolated to the quadrature points... why?
(3) Are gradients of sigma readily available in the Functional?
...
(4) I also have similar questions about FacetBasis but I'll have to expand this tomorrow.
...
(5) And I'm also trying to verify that the 'nodes' (in terms of the dual basis) used by skfem are always point evaluation nodes and that's not something that can be changed for e.g. ElementTriP1.
...
(6) Some other random confusion comes from using x to hold the solution in the examples, even though u is the variable used in the forms. Then u seems to get used again as a name of various dofs, but this doesn't make any sense to me at all.

[kinnala]

Thanks for this. The docs (scikit-fem.readthedocs.io) definitely need more information on these things under Advanced topics. The idea in the docs is that if there is a specific "goal oriented" tutorial, it will go to How-to guides. More general discussion around a topic will go to Advanced topics (we could change the name also to reflect better its contents).

(1) What if I have a functional form of sigma, e.g. sigma = np.sin(x) + np.exp(-(y**2))?

There are some options:

1. You could project this function onto a finite element basis. In finite element literature you often see that the "data" (as they call these prescribed fields) is projected onto the finite element space using L2 projection and then some "data oscillation" terms appear (which is the difference of the original function and its projection) in the error estimates. There is skfem.utils.projection for performing the L2 projection.
2. You could just substitute the global quadrature points into x and y inside the form. This is available via the last argument, e.g., w.x[0] and w.x[1].
3. It is also possible to evaluate this function at the vertices of the mesh and then consider the result as a piecewise linear finite element function. I'm sometimes lazy and do it although it may have some caveats (e.g., if the "data" is infinite at the vertices, the above approaches work while this approach doesn't). Even though not suggested, it's easy to do since you get x and y from mesh.p[0] and mesh.p[1].

[gatling-nrl]

Ah! So in

@skfem.Functional
def test1(w):
    x, y = w['x'].value
    return x**2 + y**2

x and y are quadrature points (mapped to gobal aka physical coordinates) and not doflocs?

[kinnala]

Yes!

[kinnala]

(2) I used sigma=basis.zeros() making sigma shaped like DoFS, but then i didn't need .interpolate() to get it interpolated to the quadrature points... why?

form.assemble(basis, sigma=sigma) will automatically call basis.interpolate(sigma) if sigma is ndarray. It makes sense only if you use same Basis for assembling and interpolating.

[kinnala]

Inside the form all inputs must use the same quadrature points (or you cannot add or multiply them). The point of basis0 = basis1.with_element(ElementTriP0()) is to use the same quadrature points for basis0 as are used for basis1. You could also do basis0 = Basis(mesh, ElementTriP0(), quadrature=basis1.quadrature) which in most cases is the same thing.

[kinnala]

What I mean by "in most cases" is that with_element also shares Mapping and information about the subset of elements/facets if basis1 does not integrate over the whole mesh.

[gatling-nrl]

Do the cell shape of the mesh and the integration order completely specify the quadrature points?

i.e. in basis(any_mesh, any_ElementTri, intorder=N) the ElementTri* and the N control the quadrature points?

So (except for the issue of mapping you mention) in the following basis0a and basis0b are guaranteed to be identical?

d = 10
basis1 = skfem.CellBasis(mesh, skfem.ElementTriP1(), intorder=d)
basis0a = basis1.with_element(ElementTriP0())
basis0b = skfem.CellBasis(mesh, ElementTriP0(), intorder=d)

[kinnala]

Yes I think so as long as elements or mapping kwargs are not specified.

[gdmcbain]

Do the cell shape of the mesh and the integration order completely specify the quadrature points?

No, or rather yes in scikit-fem if intorder is specified by no in the finite element method in general. One can specify any quadrature scheme that one likes. And there are a lot of possible quadrature schemes! See https://github.com/nschloe/quadpy for a mind-boggling compendium. There are in general multiple schemes of the same order. Here a quadrature scheme can be specified via the

scikit-fem/skfem/assembly/basis/abstract_basis.py

Line 45 in 68bc2aa

quadrature: Optional[Tuple[ndarray, ndarray]] = None,

attribute of an AbstractBasis (or a subclass) which is interpreted as

scikit-fem/skfem/assembly/basis/abstract_basis.py

Lines 75 to 81 in 68bc2aa

	if quadrature is not None:
	self.X, self.W = quadrature
	else:
	self.X, self.W = get_quadrature(
	refdom,
	intorder if intorder is not None else 2 * self.elem.maxdeg
	)

(i.e. intorder is ignored).

Over the years in the finite element method there have been many uses of special quadrature rules for various motivations. One selection at random from a vast literature:

• Hansbo, P. (1998). [A new approach to quadrature for finite elements incorporating hourglass control as a special case](doi: 10.1016/S0045-7825(97)00257-0). Computer Methods in Applied Mechanics and Engineering, 158:301–309

At the moment I'm working on one myself. One variant of ‘lumping the mass matrix’ (for transient or modal analysis) works by using a quadrature rule for which the abscissæ are at the nodes of the mesh; on a MeshTri, this looks like

basis["p_lumped"] = skfem.CellBasis(
    mesh, basis["p"].elem, quadrature=(basis["p"].elem.doflocs.T, np.full(3, 1 / 6))
)
``

The motivation for this is that although it is less accurate than Gaussian quadrature, the resulting mass matrix is diagonal so that Cholesky factorization is very cheap.  (Particularly important because `scipy.sparse.linalg` lacks a routine for Cholesky factorization. #86)

[kinnala]

(3) Are gradients of sigma readily available in the Functional?

They should be if you use basis.interpolate(sigma) (or if this gets automatically called inside Form.assemble due to sigma being ndarray).

[kinnala]

(5) And I'm also trying to verify that the 'nodes' (in terms of the dual basis) used by skfem are always point evaluation nodes and that's not something that can be changed for e.g. ElementTriP1.

Inside ElementTriP1 the basis is defined explicitly: https://github.com/kinnala/scikit-fem/blob/master/skfem/element/element_tri/element_tri_p1.py#L18

Here the dual basis is only conceptual and not used in practice. However, some elements are defined through the DOF functionals, e.g., ElementTriArgyris: https://github.com/kinnala/scikit-fem/blob/master/skfem/element/element_tri/element_tri_argyris.py#L37

Not all elements are defined as subclasses of ElementGlobal (like ElementTriArgyris) because those basis functions are created directly on the global element and it is not as numerically stable as using a reference element and most likely a bit slower.

It would technically be possible to build the local basis dynamically in Python using the DOF functionals and still have it computationally efficient. Trying this out is on my loooong TODO list.

[kinnala]

Should you somehow see sigma changing rapidly in the solution?

[gatling-nrl]

sigma isn't meant to be an unknown. Just a variable coefficient, and in this problem, just two horizontal blocks of 1 and 100.

[gatling-nrl]

Hmm, I made those changes and got the same result as you show. But sigma wasn't having the influence I expect... so I tweaked the definition of sigma_func a bit and added a plot of sigma_projected.

import matplotlib.pyplot as plt
import meshio
import numpy as np

import skfem
from skfem.helpers import dot, grad
import skfem.io
import skfem.models
import skfem.visuals.matplotlib

import logging
logging.basicConfig(format='%(levelname)s %(asctime)s %(name)s %(message)s')
logging.getLogger('skfem').setLevel(logging.INFO)

mesh = skfem.io.meshio.from_file('ring.msh', None).refined(3)
basis = skfem.CellBasis(mesh, skfem.ElementTriP2(), intorder=6)

float_equality = 1e-8
def is_on_inner_left_edge(x):
    return np.logical_and(
        abs(x[1]) <= 1.,
        abs(x[0] + 1.) <= float_equality
    )
fbasis = skfem.FacetBasis(mesh, basis.elem, intorder=6, facets=mesh.facets_satisfying(is_on_inner_left_edge))

sigma_background = 1
sigma_conductor = 10
def sigma_func(x):
    condcutor_region = np.logical_and(
        abs(x[1]) <= 1.,
        np.logical_and(
            x[0] <= -1,
            x[0] > -1.5
        ),
    )
    return np.where(condcutor_region, sigma_conductor, sigma_background)
sigma_projected = skfem.utils.projection(sigma_func, basis)
sigma_interpolated = basis.interpolate(sigma_projected)

def bc_val(x):
    return is_on_inner_left_edge(x).astype(float)
inner_left_dofs = fbasis.get_dofs(is_on_inner_left_edge)
bc_projected = basis.zeros()
bc_projected[inner_left_dofs] = skfem.utils.projection(bc_val, fbasis, I=inner_left_dofs)


outer_dofs = fbasis.get_dofs(
    lambda x: np.logical_or(
        abs(x[0]) >= 2-float_equality,
        abs(x[1]) >= 2-float_equality)
)
bc_dofs = np.hstack([outer_dofs, inner_left_dofs])

@skfem.BilinearForm
def my_laplace(u, v, w):
    return dot(w['sigma']*grad(u), grad(v))
A_ = skfem.asm(my_laplace, basis, sigma=sigma_interpolated)
b_ = basis.zeros()

A, b = skfem.enforce(A_, b_, x=bc_projected, D=bc_dofs)
u_projected = skfem.solve(A, b)
u_interpolated = basis.interpolate(u_projected)

fig, ax = plt.subplots(1, 2, figsize=(14,7))
skfem.visuals.matplotlib.plot(basis, sigma_projected, ax=ax[0])
skfem.visuals.matplotlib.plot(basis, u_projected, ax=ax[1], nrefs=2, colorbar=True)

This looks okay-ish, but when I use sigma_conductor = 100 it is clearly wrong

sigma_conductor = 1e6 would be more typical and looks different, but equally bad as 100.

edit: added plot titles.

[gdmcbain]

With such a strong ratio of coefficients, I'd really want to use predefined subdomains and boundaries rather than those inferred from predicates on coordinates, as recommended recently in #776 and as in

scikit-fem/docs/examples/ex17.py

Lines 58 to 60 in 68bc2aa

	conductivity = basis0.zeros()
	for s in mesh.subdomains:
	conductivity[basis0.get_dofs(elements=s)] = thermal_conductivity[s]

Is there a reason for not predefining the subdomains and boundaries?

However, this isn't the issue here. The issue here is that unlike in ex17 you are defining σ via sigma_projected which is ElementTriP2 and so supported on the vertices and midpoints of triangles and therefore on the interface between subdomains. In ex17, the thermal conductivity is only defined in the interior of elements (where all the quadrature points are) and so in the interior of subdomains. What value should σ have on the interface between subdomains? The arithmetic mean? The harmonic mean? One value or the other? Any of these choices will lead to bleed, i.e. will lead to some of the high σ spilling over into the low-σ subdomain; do not construct sigma_projected.

Another thing that I would be thinking about with such a strong ratio of coefficients is domain decomposition: solve the subproblems on the subdomains separately as in

scikit-fem/docs/examples/ex26.py

Line 37 in 68bc2aa

core_basis = Basis(mesh, basis.elem, elements=mesh.subdomains['core'])

and then couple the subproblems back together again. The subproblems should be really easy and, given the ratio of coefficients, an iterative scheme for the global coupling should converge in a couple of steps. (I intended ex26 to be the first step in such an approach, but the next step has been a while coming, sorry.)

[gatling-nrl]

With such a strong ratio of coefficients,

This is all wrong/red herring. (I strongly suspect my not yet understanding FacetBasis in whatever went wrong above.) @gdmcbain I normally agree that having some numbers be 1 and some be 1e6 is a recipe for computational disaster. It's a lot like our float == mantra. But it didn't seem right in this particular case.

In my previous experience with Laplace equation in electrostatics, it was quite okay to have huge contrast in the conductivity (or permittivity) coefficient. It didn't bother my pervious finite difference work and it didn't bother this exact same simulation in COMSOL.

It turns out it doesn't bother skfem either.

mesh = skfem.io.meshio.from_file('ring.msh', None)
basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
basis_p0 = basis_p1.with_element(skfem.ElementTriP0())

def is_on_inner_left_edge(x):
    return np.logical_and(
        abs(x[1]) <= 1.,
        np.isclose(abs(x[0] + 1.), 0),
    )
def is_on_outer_boundary(x):
    return np.logical_or(
        np.isclose(abs(x[0]), 2),
        np.isclose(abs(x[1]), 2),
    )
    
source_dofs = basis_p1.get_dofs(facets=is_on_inner_left_edge)
boundary_dofs = basis_p1.get_dofs(facets=is_on_outer_boundary)
bc_vals = basis_p1.zeros()
bc_vals[source_dofs] = 1
bc_dofs = np.hstack([boundary_dofs, source_dofs])

def is_condcuctor_region(x):
    return np.logical_and(
        np.logical_and(
            x[0] >= -1.5,
            x[0] <= 1.0,
        ),
        abs(x[1]) <= 1,
    )

sigma = basis_p0.zeros()
sigma[:] = 1
sigma[basis_p0.get_dofs(elements=is_condcuctor_region)] = 1e6

@skfem.BilinearForm
def laplace(u, v, w):
    return dot(w['sigma']*grad(u), grad(v))

A = skfem.asm(laplace, basis_p1, sigma=basis_p0.interpolate(sigma))
A, b = skfem.enforce(A, D=bc_dofs, x=bc_vals)
u = skfem.solve(A, b)

plt.subplots(figsize=(12,12))
skfem.visuals.matplotlib.plot(basis_p1, u, ax=plt.gca(), shading='gouraud')
skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())

I hope soon to write up more extended documentation about FacetBasis, after I understand how to use it!

[kinnala]

(6) Some other random confusion comes from using x to hold the solution in the examples, even though u is the variable used in the forms. Then u seems to get used again as a name of various dofs, but this doesn't make any sense to me at all.

When defining the forms you could use any names, e.g.,

@BilinearForm
def mass(p, q, _):
    return p * q

At the linear algebra level x comes from the usual notation of writing down the matrix system as Ax = b. I'm somehow very used to that because my background is in applied mathematics.

However, using strings like u, u^1, etc. inside Element.dofnames and also when querying for DOFs in DofsView is a suboptimal solution and made pretty hastily when we needed more advanced DOF queries. Ideally we would have some kind of well defined language (object hierarchy? enums?) for describing the DOFs and then use that instead. Any ideas are welcome.

[gatling-nrl]

I will think on this and write back as I learn from the outside in. I understand Ax=b but then you have del u = f for Poisson equation and now everything is a mess. And we only have so many letters. As an outsider, skfem uses confusing, terse, and sometimes inconsistent notation in the code. (I have to say I really like skfem and mean no offense... I only want to share plainly what has made it hard to use for a newcomer.) We need way more comments and/or long variable names.

[gdmcbain]

Yes I tend to avoid using x for unknowns because it does often have to represent a Cartesian coordinate in space. It's true that writing u for unknown whether it's the unknown function in the full infinite-dimensional Hilbert (or whatever—Sobolev?) space or again the vector in ℝⁿ representing that in some system of coordinates is an abuse of notation but it is done fairly often in the mathematical literature too; similarly with nondimensionalization or Fourier or other integral transforms. One starts by marking the new version of the variable by putting an h subscript on it to show that it's in the fintie element projection or a star to show that it's normalized or a hat to show that it's transformed to frequency-space but sooner or later one drops this to prevent being overwhelmed.

In the Python, I don't think we ever need to have a variable for the unknown in the full Sobolev space (unless there were to be some kind of symbolic preprocessing. #721), or for its abstract projection into a finite element subspace, only the ndarray in a particular basis, so the above doesn't apply too much to the Python code and u can be used without collision for the ndarray, and there is a sense in which this object represents the unknown (temperature, velocity, whatever) of the problem.

I think how I would arrange this in a library or ecosystem of libraries is that the linear algebraic libraries like scipy.sparse.linalg or krylov or PyAMG #273 would use x for the unknown (and they do) and functions like skfem.utils.solve should follow suit, as it does:

scikit-fem/skfem/utils.py

Lines 202 to 207 in 68bc2aa

	def solve(A: spmatrix,
	b: Union[ndarray, spmatrix],
	x: Optional[ndarray] = None,
	I: Optional[ndarray] = None,
	solver: Optional[Union[LinearSolver, EigenSolver]] = None,
	**kwargs) -> Solution:

but then at the more application-level user scripts, I'd call the unknown u and pass that to the function as x=u as in

scikit-fem/docs/examples/ex13.py

Line 45 in 68bc2aa

u = solve(*condense(A, x=u, D=basis.get_dofs({"positive", "ground"})))

The formal parameter and argument don't have to use the same symbol and perhaps shouldn't when they're in different conceptual universes (linear algebra versus partial differential equations); what's u in the calling PDE scope is x in the called LA scope. I don't know; that's my rationalization anyway.

[kinnala]

Related question: what details should we add to https://scikit-fem.readthedocs.io/en/latest/advanced.html#dofindexing to make it more useful?

[gatling-nrl]

In a functional analytic setting, DOFs are functionals that take in a function and return a value

So the first entry in x is (conceptually) the result of some g(f). If DOF[0] is associated with point (1,2) then g(f)=f(1,2) and so x[0]==f(1,2)?

But what is f? It's the polynomial covering the mesh cell associated with DOF[0]? And, if we're in P1, that just ends up working out to be equal to u at that location, when I'm solving del^2 u = y? And so it is possible to take shortcuts and interact with x in ways you cannot do when using P2 or P0?

[kinnala]

So in the examples, it sounds like you can't just change TriP1 to TriP2 or TriP0 or etc. and expect things to work, for numerous reasons. I think being able to just blindly change out the element like that isn't even a design goal of skfem?

Yeah I think it is not possible to do in general. Maybe for a subset of elements but not for everything we have implemented.

So the first entry in x is (conceptually) the result of some g(f). If DOF[0] is associated with point (1,2) then g(f)=f(1,2) and so x[0]==f(1,2)?

There are different types of DOFs. What you are describing is a point evaluation which is the simplest DOF. It could also be df/dx(1, 2) or some kind of integral integral_over_edge(f).

But what is f?

It is the finite element function, what you are drawing with skfem.visuals.matplotlib.plot.

We can also consider f and the DOFs separately in each element. However, in x we combine those DOFs that are shared between elements to get continuity.

[gdmcbain]

So in the examples, it sounds like you can't just change TriP1 to TriP2 or TriP0 or etc. and expect things to work, for numerous reasons. I think being able to just blindly change out the element like that isn't even a design goal of skfem?

Well, I wouldn't say that. There are good reasons why it's not always going to work straight away, as @kinnala says, but looking at how the docs/examples have evolved numerous commits will be seen which have made it more possible to switch the elements, e.g. from P1 to P2 or Tri to Quad, mostly because this is something that I do quite often. And I would say that facilitating that kind of switching is very much a design goal as far as I'm concerned.

One such feature is using skfem.utils.projection instead of just constructing the subset of DOFs using mesh.p, as mentioned above; projection is quite general, I believe.

• 024d5bd ‘…also doflocs->L2_projection’ (2020-04-28)

Another, enabling going from P1 to P2, is outputting .nodal_dofs instead of just assuming that the DOFs correspond to mesh.p; that does still only work for lagrangian elements.

• 9f98415 ‘generalize ex19 so that it works for Q1|Q2’ (2021-07-21)

An obstacle to going further in that direction is that not many postprocessing tools (I'm thinking of things like ParaView) understand finite elements so it's not possible to export a Raviart–Thomas function, for example. (I'd love to be corrected on this.) We have to project down to something commonly understood like ElementTriP1.

[gatling-nrl]

I've been trying to write a tutorial about the very basics of FEM

Can you take a quick look at section 1.7 of https://finite-element.github.io/L1_introduction.html#practical-implementation and let me know if the terminology is consistent with skfem? I see they use Ku=f where I think you use Ax=b. Are there other important differences?

[kinnala]

Yes, that seems pretty close to what I wrote in branch fem-intro-tutorial. I'll check it more carefully soon.