What condition ensures the global continuity of the solution in the FEM?

Question

I know this is trivial but I don't seem to understand it. In which step of the FE formulation do we enforce the global continuity of the solution? Or in other words, how the construction of the local basis functions on each element allows to glue the solution with continuity across neighbouring elements?

To explain what I mean I'll use a toy example. Considering the following abstract variational problem $$ \text{find}\;u\in H^1(\Omega)\;\text{s.t.}\\ a(u,v)=f(v)\quad \forall v\in H^1(\Omega) $$ the Galerkin method allows me to formulate it in a finite-dimensional subspace $V_h\subset H^1(\Omega)$; in the FEM's pipeline I then go ahead and discretise the domain into a partition $\mathcal{T}_h$ of geometrical elements $\mathcal{K}$ (e.g. triangles) and specifically choose the subspace to be $V_h=\{u_h\in\mathcal{C}(\Omega)\;\text{s.t.}\;u_h\vert_{\mathcal{K}}\in\mathbb{P}_1\;\forall\mathcal{K}\in\mathcal{T}_h\}$ where the solution is locally interpolated by a linear function. Skipping a few steps I have now reduced the variational problem into a discrete algebraic one by testing the weak form against each basis function of $V_h$ (test functions) and by expressing the solution in such basis (trial functions) $$ \sum_{j=1}^nu_j\,a(\varphi_j,\varphi_k)=f(\varphi_k),\;\forall k=1,\dots,n $$ The above is a linear system of $n$ equations in $n$ unknowns which coincide precisely with the value of the numerical solution on the interpolation nodes since the property $\varphi_j(\mathbf{x_k})=\delta_{jk}$ (one unique global basis function with compact support on each node of the partition) entails that $$ u_h(\mathbf{x}_k)=\sum_{j=1}^nu_j\varphi_j(\mathbf{x}_k)=\sum_{j=1}^nu_j\delta_{jk}=u_k $$ The number of degrees of freedom $n$ coincides with the number of nodes of interpolation of either the solution itself or some of its derivatives; this depends on the type of FE element I choose but here, for the sake of simplicity I'll just refer to Lagrangian elements.

Now the crucial step (the one I get lost): to assemble the linear system I must compute the bilinear form $a(\cdot,\cdot)$ and linear form $f(\cdot)$ which let's pretend are simple volume integrals. To do that I enforce the local interpolation choice that I made when selecting $V_h$, i.e. I use linear basis functions. Then I map every physical triangle $\mathcal{K}$ into a reference triangle $\hat{\mathcal{K}}$ via a bijective affine transformation (this will then require the computation of the determinant of the Jacobian etc...). For the current example of Lagrangian $\mathbb{P}_1$ element the basis functions are: $$ \hat{\varphi}_1=1 - \hat{x} - \hat{y},\quad\hat{\varphi}_2=\hat{x},\quad\hat{\varphi}_3=\hat{y} $$ All this is absolutely clear to me, including the steps that come next. What is not clear to me is how this choice of local linear interpolation reflects into the global continuity of the solution. I read a lot of both textbooks and papers but, although they do discuss the properties of FEs of class $\mathcal{C}^0$ (like the Lagrangian ones) and $\mathcal{C}^1$ (like the Hermitian ones), they never explicilty address the link that there is between the construction of ad-hoc local basis functions and the glueing process that comes next.

I am sure they don't because it's a very obvious result however in my brain I see at the locality of the choice of the basis function and the globality of the property of continuity to be totally disconnected and arbitrary. If anything in my reasoning above is wrong or ambiguous please do point out and I'll amend it.

Answer 1

For the linear elements you consider, you are mapping the shape functions from the reference cell to each of the cells of your mesh. The important properties you are using here are:

• The shape functions on the reference cells are linear functions.
• As a consequence, the shape functions are also linear on each edge of the reference cell.
• The mapping is linear.
• As a consequence, the shape functions on actual cells are also linear and, in particular, linear on each edge of each edge.

If you now have shape functions on neighboring cells that agree in their value at the nodes of the common edge between the two cells, then the fact that they are linear along the edge means that they agree along the entire edge. In other words, the function is continuous.

Answer 2

After much work and frustration, although eased by the answers and comments above that guided my researches to more detailed queries, I finally found the solution to my problem.

It consists of a mix of the answers gave by both Wolfgang and Aage and also the comment of knl under the original question. The latest contribution of tietäjä, in the form of a literature reference, also came useful to double-check the exactness of the solution.

DISCLAIMER

The full answer turns out to be trivial (as expected) and I apologise to the community if my question wasn't clear enough to highlight the core of my doubt. I am reformulating the problem here in a more general and concise framework in the hope of clarifing what I was looking for in the post:

Given a FE (according to its axiomatic definition provided by [Ciarlet, 1978]) in the form of the triple $\{\mathcal{K}, V_{\mathcal{K}}, \mathcal{L}_{\mathcal{K}}\}$, where $V_{\mathcal{K}}$ is a functional space of (finite) dimension $N_{\mathcal{K}}$ and $\mathcal{L}_{\mathcal{K}}$ is a set of $N_{\mathcal{K}}$ (unisolvent) linear functionals $\ell_j(v)$ for $V_{\mathcal{K}}$ (degrees of freedom), then what conditions on the local basis functions built on $\mathcal{K}$ ensures that the derived global basis functions, defined over $\mathcal{T}=\cup\,\mathcal{K}$, will be s.t. the solution is continuous?

This problem is closely related to the well-known issue of building FEs of class $\mathcal{C}^0(\Omega)$ (s.a. the Courant or quadratic Lagrangian ones) or $\mathcal{C}^1(\Omega)$ (s.a. the HCT or Argyris ones); they are so related (if not identical) in fact that the answers provided by the aforementioned members all somewhat referred to the latter. My problem however was concerning the conditions that the local basis functions must satisfy "on their own" and not when compared to those of contiguous elements. The proofs that I read in many references were all addressing the continuity of the edge restriction of the basis functions by approaching the problem of glueing the values of the degrees of freedom along the edge shared between two continguos elements.

To reiterate, this was not the problem I was looking for: after all, when we specify the shape of the local basis functions (in practice and code implementations) we do so on a reference/unitary element without (apparently) imposing the edge-wise continuity conditions that are used in the proofs mentioned above. This is not true obsiously since those basis functions (e.g. the Courant ones) are indeed edge-wise continuous, however the reason why I posted the question is that I was not able to find a clear and intuitive reason for it.

In the following I will provide a short, brief answer followed by a full one complemented by a numerical example. I am posting the answer as a community wiki for two reasons:

• there is a strong possibility that I'll make some mistakes both conceptual and practical;
• it seems to me more appropriate not getting any of the potential upvotes from it since it was indeed an exploitation of the collective pieces of knowledge I gathered from the previous answers provided by other members

Short answer

The following definition specifies the local conditions that the set of linear functionals (DOFs) of a given FE must satisfy to ensure global continuity independently from contiguous elements.

The FE of class $\mathcal{C}^0(\Omega_h)$ is a triple $\{\mathcal{K}, V_{\mathcal{K}}, \mathcal{L}_{\mathcal{K}}\}$ s.t. for each edge $L$ of $\mathcal{K}$ there is a subset of degrees of freedom $\mathcal{L}_L:=\mathcal{L}_{\mathcal{K}}\vert_L\subset\mathcal{L}_{\mathcal{K}}$ restricted to such edge for which

• $\forall\ell\in\mathcal{L}_L\;\:\ell(v)=\lambda(v\vert_L)$ i.e. the restriction of the DOFs on the edge $L$ only depends on the values taken by their argument $v\in V_{\mathcal{K}}$ on said edge;
• $\forall\ell\in\mathcal{L}_L\;\:\ell(v)=0\iff v\vert_L\equiv0$ i.e. the restriction of the DOFs is unisolvent for $V_{\mathcal{K}}\vert_L$

In practicality, when we define a FE, given a geometrical element $\mathcal{K}$ (e.g. a triangle), a local approximation space $V_{\mathcal{K}}$ (e.g. $\mathbb{P}_1$) and a set of unisolvent DOFs (e.g. $\ell_j(v)=v(\mathbf{x}_j$) we must ensure that

1. The local basis functions $\hat{\phi_k}$ constructed on the reference/unitary element $\hat{\mathcal{K}}$ are s.t. $$\ell_j(\hat{\phi_k})=\delta_{jk},\;\;\forall j,k=1,\dots,N_{\mathcal{K}}$$
2. The 1-to-1 correspondance $j_g\leftrightarrow\{j,\mathcal{K}\}$ between the global numbering $j_g=1,\dots,N$ of the DOFs (e.g. nodes) and the local numbering $j=1,\dots,N_{\mathcal{K}}$ on the reference element is consistent for each geometrical element of the partition. This is done by specifying the same order by which the global nodes of an element are numbered (e.g. counterclockwise) on the local element regardless on which one is the first.

Full answer

Let's start with the traditional proof of FEs of class $\mathcal{C}^0(\Omega)$ that is commonly found on FEs textbooks and literature. If we consider two continuos elements $\mathcal{K}_1,\mathcal{K}_2$ s.t. $L:=\partial_{\mathcal{K}_1}\cap\partial_{\mathcal{K}_2}$ and the local functions $v_1\in V_{\mathcal{K}_1},v_2\in V_{\mathcal{K}_2}$ defined on those respectively then, in order for a global function $v$ defined on both elements (s.t. $v\equiv v_1$ on $\mathcal{K}_1$ and $v\equiv v_2$ on $\mathcal{K}_2$) to be continuous, the restrictions on $L$ of $v_1$ and $v_2$ must exactly coincide. This is intuitive enough even with a simple geometrical mind-experiment: suppose $\mathcal{K}_1, \mathcal{K}_2\subset\mathbb{R}^2$ and take $v_1,v_2\in\mathbb{P}_1$, i.e. they are portions of flat planes hanging over the two elements (suppose they are triangles); to obtain a solution $v$ that is continuous over the union of the two elements we must glue together the two planes along the common edge. If $v_1,v_2$ are generic functions this would require pointwise matching of the values taken along the edge; however, in this example, we are (locally) approximating the FE solution using linear polynomials and since the restriction of a bivariate linear ploynomial on an line is itself a linear polynomial we conclude that only need to ensure that $v_1\equiv v_2$ on two separate points of the edge. Here's the connection between the local degrees of freedom and the global continuity property of the basis functions of the functional space: in fact if we choose Lagrangian FEs, then the (unisolvent set of) DOFs defined on $\mathbb{P}_1$ coincide with the values taken by $v_1,v_2$ on the nodes of the element. Therefore we conclude that, if we set $\mathbf{x}_A,\mathbf{x}_B$ to be the endpoints of $L$, we have

$$v\in\mathcal{C}^0(\mathcal{K}_1\cup\mathcal{K}_2)\iff\ell^{(\mathcal{K}_1)}_{j}(v_1)=v_1(\mathbf{x}_j)\equiv v_2(\mathbf{x}_j)=\ell^{(\mathcal{K}_2)}_{j}(v_2),\;\;j\in\{A,B\}$$

This is what Wolfgang was referring to in its answer; Aage went further generalising on the choice of the DOFs to be, for example, the values of the polynomials at the edges midpoints of the triangles (rather than their vertices). This FE, as he mentioned, is clearly not class $\mathcal{C}^0$ since along each edge we cannot specify enough constraints on the DOFs to enforce continuity.

This is easily generalisable to higher-order Lagrangian approximations (i.e. quadratic, cubic, quintic polynomial spaces) since we will only need to increase the number of edge's nodes on which to enforce the equivalence of the DOFs. Furthermore we can also build $\mathcal{C}^1$ FEs by imposing very similar conditions on some of the derivatives of the polynomials (which is the case for Hermitian elements). What's wrong with this? Nothing except the fact that we are imposing an equivalence condition for the DOFs when considering a couple of contiguos elements and not when building the FE on its own. Take the basis functions $\hat{\varphi}(\hat{x},\hat{y})$; they are defined on their own and there is no apparent and starightforward way to verify that those do indeed satisfy the above continuity condition.

The subtle issue (perhaps not su subtle but at least for me it was not evident from the beginning) was then partially addressed by knl comment and finally complemented by Wolfgang comment under his answer. The definition given in the short answer above connects the global continuity with a local condition on the DOFs and thus on the local basis functions.

Let's address the point 1. first. The global basis functions (which we indicate with $\varphi_j$) are defined over a (compact) support that corresponds to the union of contiguous elements of the partition that each share $\mathbf{x}_j$ as one of their node (vertex, midpoint, ...). This function is only piecewise polynomial, in particular it's a polynomial of degree e.g. 1 when we consider its restriction to a single element of its support. If we say that such support is made of $N_j$ elements then we may also state that the global basis function itself is the union (which is an abuse of notation for the sake of clarity) of the local basis functions $\phi_{s}^{(\mathcal{K})}$ defined on each triangle $\mathcal{K}$ for the same node $\mathbf{x}_j$. For $\varphi_j$ to be continuous we need each local basis function to glue with continuity across each shared edge. If we thus specify, for the reference element, that the local basis functions must satisfy $\ell_j(\phi_k)=\delta_{jk}$ then, of the infinitely many linear polynomials that we can construct on the reference element for each node we are interested only in those that have unitary value on the prescribed node and null values on the remaining twos.

Now onto point 2. to complete the discussion. It is clear that 1. alone cannot automatically glue together the local basis functions along the shared edges of the support of $\varphi_j$. To do this we must choose a specific local order of numbering for the global nodes of an element and be consistent for every other element. This allows us to map coherently each physical element to the reference element in such a way that different local basis functions, defined on different contiguous triangles but for the same global node, agree on their shape along the shared edge.

This condition is easily verified with a numerical example. Take $\Omega_h=[-1,1]\times[-1,1]$ s.t. its partition is made of $8$ identical triangles as in the figure below.

We choose Courant FE i.e. $V_\mathcal{K}=\mathbb{P}_1$ and $\mathcal{L}_{\mathcal{K}}=\{\ell_j(v)=v(\mathbf{x}_j)\}$ for which we know that the basis functions on the reference element are the ones mentioned in the original answer. The global basis function associated with node $\mathbf{x}_{1}=(0,0)$ (i.e. global index $j_g=1$) is the hat/tent function displayed in the figure below.

We can think of it as the intersection between 4 planes:

$$ f_1(x,y)=1-x\\ f_2(x,y)=1+x\\ f_3(x,y)=1-y\\ f_4(x,y)=1+y $$ Our objective is to build the global basis function $\varphi_1$ Now, we must verify that the set of reference basis functions $\hat{\varphi}_j,\;j=1,2,3$ (which it's easy to see that they do satisfy the condition ${\varphi}_j(\hat{\mathbf{x}_k})=\delta_{jk}$), coupled with a consistent counterclockwise local numbering of the global elemental nodes to the unitary element does indeed $\varphi_1$.

We first specify the affine transformation $(x,y)^T=\mathbf{F}[\hat{x},\hat{y}]=\mathbf{B}\hat{\mathbf{x}}+\mathbf{x}_1$ that bijectively maps any physical element $\mathcal{K}$ into the reference $\hat{\mathcal{K}}$ since $\mathbf{B}=\begin{pmatrix}x_2-x_1 & x_3-x_1 \\ y_2-y_1 & y_3-y_1\end{pmatrix}$. The nodes $\mathbf{x}_j,\;j=1,2,3$ correspond to the local numbering of the global indices of the mesh in the previous figure; for each triangle we can choose any global node that we want to be $\mathbf{x}_1$ as long as we keep the same counterclockwise numberinh order for the remaining two nodes (i.e. be consistent). For the sake of simplicity we set $j_g=1\leftrightarrow\{j=1,\mathcal{K}\}$ i.e. for each triangle we set the global origin $\mathbf{x}_1=(0,0)^T$ to be the starting node for the local numbering. This choice in fact reduces the affine map to a mere rotation and scaling of the triangle, disregarding the translation, i.e. $\mathbf{F}=\mathbf{B}\hat{\mathbf{x}}$.

Let's compute $\phi_1^{(V)}$ and $\phi_1^{(IV)}$ and show that they in fact correspond to the restriction of plane $f_1(x,y)$ to the contiguous triangles $\mathcal{K}_V$ and $\mathcal{K}_{IV}$. The same operations can be done for all couples of triangles to check the glueing of the local basis functions.

Triangle V

$$ j_g=1\to j=1,\quad j_g=6\to j=2,\quad j_g=7\to j=3\\ \mathbf{B}=\begin{pmatrix}1 & 1 \\ 0 & 1\end{pmatrix}\;\Rightarrow\;\mathbf{B}^{-1}=\begin{pmatrix}1 & -1 \\ 0 & 1\end{pmatrix}\;\Rightarrow\;\mathbf{F}^{-1}=\mathbf{B}^{-1}\mathbf{x}=\begin{pmatrix}x-y\\y\end{pmatrix}\\ \phi_1^{(V)}=\hat{\varphi}_1\circ\mathbf{F}^{-1}=(1-\hat{x}-\hat{y})\circ\begin{pmatrix}x-y\\y\end{pmatrix}=1-(x-y)-(y)=1-x\equiv f_1(x,y)\vert_{\mathcal{K}_{V}} $$

Triangle IV

$$ j_g=1\to j=1,\quad j_g=5\to j=2,\quad j_g=7\to j=6\\ \mathbf{B}=\begin{pmatrix}1 & 1 \\ -1 & 0\end{pmatrix}\;\Rightarrow\;\mathbf{B}^{-1}=\begin{pmatrix}0 & -1 \\ 1 & 1\end{pmatrix}\;\Rightarrow\;\mathbf{F}^{-1}=\mathbf{B}^{-1}\mathbf{x}=\begin{pmatrix}-y\\x+y\end{pmatrix}\\ \phi_1^{(IV)}=\hat{\varphi}_1\circ\mathbf{F}^{-1}=(1-\hat{x}-\hat{y})\circ\begin{pmatrix}-y\\x+y\end{pmatrix}=1-(-y)-(x+y)=1-x\equiv f_1(x,y)\vert_{\mathcal{K}_{IV}} $$

We can easily demonstrate that, as long as we keep consistent the ordering of the nodes it does not matter which node we do number first; for example for Triangle V let's change the first node not to be the global origin but the upper right corner. Now the global origin will coincide with the second locally numbered node according to the convention we specified. Therefore we'll have

Triangle V

$$ j_g=7\to j=1,\quad j_g=1\to j=2,\quad j_g=6\to j=3\\ \mathbf{B}=\begin{pmatrix}-1 & 0 \\ -1 & -1\end{pmatrix}\;\Rightarrow\;\mathbf{B}^{-1}=\begin{pmatrix}-1 & 0 \\ 1 & -1\end{pmatrix}\;\Rightarrow\;\mathbf{F}^{-1}=\mathbf{B}^{-1}(\mathbf{x}-\mathbf{x}_1)=\begin{pmatrix}-(x-1)\\(x-1)-(y-1)\end{pmatrix}\\ \phi_1^{(V)}=\hat{\varphi}_1\circ\mathbf{F}^{-1}=(\hat{x})\circ\begin{pmatrix}1-x\\x-y\end{pmatrix}=1-x\equiv f_1(x,y)\vert_{\mathcal{K}_{V}} $$

which reresults in the identical global basis function.

Answer 3

For the P1-element, global continuity is enforced by choosing function values at the vertices as degrees of freedom. If you instead choose for example function values at the midpoint of the edges you don't get global continuity (only continuity at the midpoint of the edges). Global continuity thus depends on the choice of local basis functions and the choice of degrees of freedom.