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

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?

• I ‘m planning another answer, yes. Please be patient. – Han de Bruijn Jun 13 at 11:25
• @HandeBruijn Ok, sounds good. Thanks! – WaltherArgyris Jun 13 at 11:58
• And here is a List of FEM related publications by this author, if you are interested. – Han de Bruijn Jun 15 at 9:21
• Please take care while reading Zienkiewicz. You know, my personal experience is in numerical analysis with (internal) fluid flow and heat transfer. Zienkiewicz may be good in structural mechanics, but when trying to apply his FEM on my sort of problems, I found that the content of the book can be often highly misleading. I have an old version at my disposal; don't know if later editions of the book show significant improvement concerning the issues encountered at that time. – Han de Bruijn Jun 15 at 12:46
• @HandeBruijn Dear Han, sorry for my late feedback. I will have a look at your reply. Thanks for your comments on Zienkiewicz' books! – WaltherArgyris Jun 15 at 19:53

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:

1. 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.

• Thank you for your answer (by the way, it would interesting to hear from the downvoter if there is anything wrong about the answer)! I will hold off a little on accepting your answer and awarding the bounty, in order to see if there is further input. I hope you don't mind. – WaltherArgyris Jun 12 at 8:30
• I have a few questions that I would like to ask to your answer: – WaltherArgyris Jun 12 at 8:31
• 1) Do you have any good references to suggest where this is further discussed? Is it properly discussed in standard books like Zienkiewicz and Taylor? I'm very new to FEM, so I don't have much experience with the standard literature. 2)Say that my domain is not polytypal and that I don't have access to finite elements with curved domain boundaries. How do I actually make a correct implementation then? How would f_b actually be defined for non-constant g? 3)You're right that my example is singular. However, I could have added some Dirichlet conditions without changing the nature of my question. – WaltherArgyris Jun 12 at 8:41
• Ah, the panacea of curvilinear meshes. One should be extremely cautious with them, as can easily be demonstrated for the 1-D case: Understand 1D FEM solution using quadratics elements – Han de Bruijn Jun 14 at 19:22
• @WaltherArgyris I've added some more to answer which hopefully answers some of your questions and points in the right direction :) – Ellya Jun 14 at 21:02

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:

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.

• I think it is important to note that many PDE, with (hyperbolic) conservation laws as a key example, do not admit smooth or $C^k$ solutions. In fact hyperbolic PDEs have driven the development of the discontinuous galerkin FEM, due to the irregularity of solutions. – Ellya Jun 14 at 21:05
• @Ellya. You're quite right. I am not satisfied with my answer anyway and I'll withdraw it soon. – Han de Bruijn Jun 15 at 9:10
• @HandeBruijn Thank you for your reply, and sorry for my late feedback. It is indeed an interesting observation that CA is also an idealization. However, the nature of my question was rather to address the issues I raised in my post under the assumption of CA being a good description. Still, I appreciate your input, especially regarding how to practically do the implementation. It would be interesting to know how it's done in commercial software, but maybe it's a trade secret. – WaltherArgyris Jun 15 at 20:18
• @WaltherArgyris: okay, then I leave the answer, for what it’s worth. – Han de Bruijn Jun 16 at 8:16
• @HandeBruijn Ok, thanks! – WaltherArgyris Jun 16 at 9:34