Conservation laws in the finite element method for domains with curved boundaries

Question

Consider, for concreteness, the finite element (FE) method applied to stationary heat conduction in a domain $\Omega \subset \mathbb{R}^3$. Let the heat flux (thermal energy per unit area and time) be denoted by $\mathbf{q}$, and let the source term (applied thermal energy per unit volume and time) be denoted by $Q$. The strong form of the equation is simply given by the balance law $$1) \nabla \cdot \mathbf{q}=Q.$$ Suppose that we prescribe Neumann conditions on the whole boundary, i.e. $\mathbf{q}\cdot \mathbf{n} = g$ on $\partial \Omega$, where $\mathbf{n}$ denotes the outward unit normal, and $g$ denotes some function defined on $\Omega$. Integrating the conservation law 1) over $\Omega$ and using Gauss' theorem along with the boundary conditions, shows that $g$ and $Q$ must be related according to $$2)\int_\Omega QdV =\int_{\partial \Omega}gdA.$$ Consider now the FE formulation of 1). Let's assume a constitutive relation $\mathbf{q} = -\mathbf{D}\nabla T$, where $\mathbf{D}$ denotes the constitutive matrix and $T$ the temperature field. Let $\mathbf{N} = [N_1 \ ... \ N_n]$ denote the global shapefunction vector. The FE formulation then takes the form $$3)\mathbf{K}\mathbf{a}=\mathbf{f}_b+\mathbf{f}_l=\mathbf{f},$$ where $\mathbf{K}$ is the stiffness matrix given by $$4)\mathbf{K} = \int_{\Omega}(\nabla \mathbf{N})^T\mathbf{D}(\nabla \mathbf{N})dV$$ and we have the boundary force vector $\mathbf{f}_b$ and the load force vector $\mathbf{f}_l$, given by $$5)\mathbf{f}_b = -\int_{\partial \Omega}\mathbf{N}^TgdA $$ and $$6)\mathbf{f}_l = \int_{\Omega}\mathbf{N}^TQdV, $$ respectively.

Question 1

My first question, which was not really addressed in the book I'm using, concerns how the integrals in 4), 5) and 6) are actually defined. The shape functions are of course only defined on the elements in the FE mesh, whose domain we can call $M$. If our elements have straight boundaries but $\Omega$ does not have a straight boundary, then $M \neq \Omega$ and $\partial M \neq \partial \Omega$. Therefore, the integrals in 4), 5) and 6) should really be written as $$4')\mathbf{K} = \int_{M}(\nabla \mathbf{N})^T\mathbf{D}(\nabla \mathbf{N})dV,$$ $$5')\mathbf{f}_b = -\int_{\partial M}\mathbf{N}^TgdA $$ and $$6')\mathbf{f}_l = \int_{M}\mathbf{N}^TQdV, $$ respectively. However, then another issue appears. While the integrals in 4') and 6') are well-defined, since $M \subset \Omega$ and $Q$ is defined on $\Omega$, the integral in 5') is not well-defined a-priori, since $\partial M \not\subset \partial \Omega$ while $g$ is defined on $\partial \Omega$. Say for example that $\Omega$ is the quarter ball $\Omega = \{(x,y,z)|x^2+y^2+z^2 \leq 1, x,y,z, \geq 0\}$ and that we use a single element in the form of a tetrahedron, so that $M$ is the tetrahedron with vertices $(0,0,0),(1,0,0),(0,1,0)$ and $(0,0,1)$. Then the "spherical part" $(\partial\Omega)_s =\{(x,y,z)|x^2+y^2+z^2=1\}\cap \partial \Omega$ of $\partial\Omega$ does not coincide with the corresponding face $(\partial M)_s = \{(x,y,z)|x+y+z=1\}\cap \partial M$ of $\partial M$. If $g$ is constant, it probably makes sense to define $\mathbf{f}_b$ as $$5'')\mathbf{f}_b = -g\int_{\partial M}\mathbf{N}^TdA \ (\mathrm{for} \ g \ \mathrm{constant}).$$ However, say that we have defined $g$ on $(\partial \Omega)_s$ by $g|_{(\partial \Omega)_s} = \sqrt{x^2+y^2+z^2-0.9}$, which is well defined. However, we cannot just use the same "formula" $g(x,y,z) = \sqrt{x^2+y^2+z^2-0.9}$ when integrating over $(\partial M)_s$ since e.g. at the point $(1/3,1/3,1/3)$ the argument of the square root becomes negative. In this case, how should one then actually define the integral $-\int_{(\partial M)_s}\mathbf{N}^TgdA$?

Question 2

My second question concerns conservation laws in the FE method, and connects to the issues raised in the first question. In the book that I am using, the authors claim that the FE method exactly reproduces the conservation law in 1). They note that the shape functions satisfy the identity $$7) \sum_{i=1}^n N_i = 1.$$ If we define a total force vector $\mathbf{f}$ by $\mathbf{f}=\mathbf{f}_b+\mathbf{f}_l$ and use equation 7), then from 2), 5) and 6) we obtain $$8)\sum_{i=1}^n f_i = 0. $$

However, as we discussed above, we cannot actually use 5) and 6) to calculate $\mathbf{f}_b$ and $\mathbf{f}_l$, but we must rather use 5') and 6'). Doing so, the relation 8) will of course no longer hold. I argued this point in a course I took based on the aforementioned book, and the reply I got (unless I misunderstood it) was that the FE method does indeed exactly implement the conservation law 1). However, based on the arguments just presented, I don't see how that can be the case.

Thus, my question is whether the FE method does indeed correctly implement the conservation law 1), and if so, how?

Answer 1

The short answer is yes, your intuition is correct, I would also have argued if someone had said the same thing. The approximation of a non polytopal domain by a straight mesh/ triangulation is called a "variational crime", due to the fact that consistency/ Galerkin orthogonality properties are violated.

Unfortunately your experience is very common in finite element literature and my colleague referred to this as the "smooth polygon" assumption, where an author assumes a domain is smooth, so that nice properties hold for the PDE, but then they consider example PDEs on a polytopal domain, so that it can be meshed exactly. For example IF the domain is polytopal, then the deduction is correct, and the conservation law can be replicated.

Otherwise, one must use modified finite element spaces, either by using a curvilear mesh, or a straight edge Triangulation that uses nonaffine maps to define the finite element space. Also, in the special case where $Q$ and $g$ are constant, then it can be conserved, as the source functions are then geometry independent.

Something else that also seems to have been overlooked is that the pure Neumann problem is singular, since any constant fuction is in the kernel of the operator. To overcome this, either a lagrange multiplier should be used (and therefore the variational formulation changes, and so the conservation law won't be replicated), or one has to modify the basis functions so that they integrate to zero (so that the FE solution integrates to zero, and we remove the singularity). If this is the case, then we lose the property $\sum N_i=1$, and the main deduction cannot be made.

---

Update: It seems that majoritively, your question is about what I have referred to as a "variational crime" in the inexact approximation of the computational domain, rather than necessarily the violation of a conservation law. I will try to summarise a few points:

1. Variational crime: - generally corresponds to any case in which the variational form is inexact. Polytopal approximation of a curved domain is one such example, but also so is e.g. numerical integration which ubiquitous in FEM and many other applications as it allows for the efficient and accurate approximation of integrals. See e.g. these lecture notes: https://people.maths.ox.ac.uk/suli/fem.pdf

2. Accuracy: In the case of the polytopal approximation of a curved domain, e.g., the unit disk, one will notice a cap on convergence rates, I have observed in many cases that if one approximates the unit disk with a simplicial mesh, even when the solution is smooth, the convergence rates are capped at $O(h^2)$

3. Literature: a good bit of the original literature is from the 70s & 80s:

Capturing the domain exactly:

Christine Bernardi: Optimal Finite-Element Interpolation on Curved Domains, https://www.jstor.org/stable/2157925

Scott, L. Ridgway: Finite element techniques for curved boundaries, https://dspace.mit.edu/handle/1721.1/12182

Isoparametric (polynomial) domain approximation: Optimal Isoparametric Finite Elements and Error Estimates for Domains Involving Curved Boundaries, M. Lenoir, https://www.jstor.org/stable/2157524.

I can also shamelessly plug some of my own work here:

• https://arxiv.org/abs/1903.08735 - Curved DGFEM
• https://arxiv.org/pdf/1807.03535.pdf - Optimal transport Monge Ampere equations (experiments demonstrate the $O(h^2)$ cap.)

4. Implementation: This is a good question. One approach would be to extend the known boundary data into the full domain, and approximate that on the boundary instead. In the case of piecewise linear finite elements, all of the boundary degrees of freedom would lie on the true boundary, so in that case you would be able to interpolate the data.

Another point to consider here: if your true domain is not convex, the polytopal approximation could lie both inside and outside of the true domain.

Answer 2

There are at least three different disciplines involved with the question:

• CA. Classical Analysis. Think of limits, continuity, derivatives, infinitesimals
• NA. Numerical Analysis. Think of Finite Elements, Finite Differences, Finite Volumes
• PR. Physical Reality.

Let's start with the last one. Stationary heat conduction inside a solid body is described by a (set of coupled) partial differential equation(s). However, it is clear from the outset that the infinitesimal quantities in such equations must be considered as a mere illusion, an idealization to be precise. Especially the differential volumes are not really infinitely small. Far from that! They should contain a lot of molecules, if we want the approximations, yes, the approximations of Classical Analysis to be considered as valid. There is a beautiful text about this (thanks to Jean Baptiste Perrin) in Benoit Mandelbrot's 'The Fractal Geometry of Nature', first few pages (6,7). Unfortunately, I don't have the book in my shelves and I can't exactly quote. Anyway, while thinking about it, two things should become clear, if PR is to be adopted as our truth frame of reference:

1. The infinitesimal volumes of Classical Analysis are too small when compared to PR
2. Likewise the finite elements, differences, volumes of Numerical Analysis are too large

With a pinch of salt, because boundaries between the disciplines are not always crisp. But I hope the message is clear: Classical Analysis is not "better" than Numerical Analysis per se. Neither is NA superior to CA (but it seems that the majority of people does not think that anyway).
If I understand it well, the OP's questions boil down a great deal to the following: trying to comprehend what the relationship is between numerically and analytically. I think the best approach is that CA and NA are quite different worlds, two entirely different ways of doing mathematics, though (hopefully) corresponding with the same physical reality in the end.
Perhaps the most important difference between the analytical world and the numerical world is that their function spaces are completely different. In the analytical world, functions most of the time have continuous n-th order (partial) derivatives. In the numerical world, functions are often piecewise linear. This means that continuity already stops at low orders of the derivatives. With all this in mind, it's probably less difficult to accept that equations 4) 5) 6) belong to CA and dfferent equations 4') 5') 6') belong to NA.

Question 1.

But let's end the prose for a while and calculate. As a generalization of the OP's tetrahedron, consider one with vertex $(0,0,0)$ at the origin and three vertices at the surface of the unit sphere; let the latter form an equilateral triangle with edge lengths $=a$. Then with some elementary geometry and algebra we find for the distance $\,h\,$ between the origin and the midpoint of any such equilateral triangle: $$ h = \sqrt{1-a^2/3} $$ For the sample tetrahedron $(0,0,0),(1,0,0),(0,1,0),(0,0,1)$ we calculate $\,h=1/\sqrt{3}\,$.
Now there are two ways to look at the OP's dilemma concerning $\,g(x,y,z) = \sqrt{x^2+y^2+z^2-0.9}\,$:

1. Stay within the world of NA and accept that the triangles employed at the surface of the unit sphere form a piecewise linear interpolation of the values at the vertices and nothing else. Then the solution is that $g(\mbox{midpoint})$ is the mean ($1/3)$ of the $3$ function values $\sqrt{0.1}$ at the vertices, which is simply $\sqrt{0.1}$.
2. Staying within the world of CA and demanding that the argument of the square root remains positive. The common trick to accomplish this is mesh refinement, in our case as follows: $$ h^2-0.9 \ge 0 \quad \Longrightarrow \quad 1-a^2/3-0.9 \ge 0 \quad \Longrightarrow \quad a\le\sqrt{0.3} $$

Question 2.

Concerning the second question, it is indeed the case, or rather: it should be the case that conservation laws hold in the analytical world as well as in the numerical world. There is one major difference, though. The conservation laws in the analytical world hold for an infinitesimal volume as well, but the conservation laws in the numerical world only hold for finite volumes. I do not consider these isues as trivial and have worked out some in answers to the following question:

• What is the difference between Finite Difference Methods, Finite Element Methods and Finite Volume Methods for solving PDEs?

Perhaps there are other ways, but I have done it as follows. Conservation laws hold by definition if finite volume methods (FVM) are employed. So if we are able to to establish a mapping between finite elements and finite volumes, then it is for sure that conservation laws also hold with FEM.