The major difference between deriving variational formulations in 2D and 3D compared to 1D is the rule for integrating by parts. A typical second-order term in a PDE may be written in dimension-independent notation as

$$ \nabla^2 u \quad\hbox{or}\quad \nabla\cdot\left( a(\x)\nabla u\right) \tp $$ The explicit forms in a 2D problem become $$ \nabla^2 u = \nabla\cdot\nabla u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}, $$ and $$ \nabla\cdot\left( a(\x)\nabla u\right) = \frac{\partial}{\partial x}\left( a(x,y)\frac{\partial u}{\partial x}\right) + \frac{\partial}{\partial y}\left( a(x,y)\frac{\partial u}{\partial y}\right) \tp $$ We shall continue with the latter operator as the form arises from just setting $ a=1 $.

The general rule for integrating by parts is often referred to as Green's first identity:

$$ \begin{equation} -\int_{\Omega} \nabla\cdot (a(\x)\nabla u) v\dx = \int_{\Omega} a(\x)\nabla u\cdot\nabla v \dx - \int_{\partial\Omega} a\frac{\partial u}{\partial n} v \ds, \tag{95} \end{equation} $$ where $ \partial\Omega $ is the boundary of $ \Omega $ and $ \partial u/\partial n = \normalvec\cdot\nabla u $ is the derivative of $ u $ in the outward normal direction, $ \normalvec $ being an outward unit normal to $ \partial\Omega $. The integrals $ \int_\Omega ()\dx $ are area integrals in 2D and volume integrals in 3D, while $ \int_{\partial\Omega} ()\ds $ is a line integral in 2D and a surface integral in 3D.

$$ \begin{align} \v\cdot\nabla u + \alpha u &= \nabla\cdot\left( a\nabla u\right) + f, \quad\x\in\Omega,\\ u &= u_0,\quad\x\in\partial\Omega_D,\\ -a\frac{\partial u}{\partial n} &= g,\quad\x\in\partial\Omega_N \tp \end{align} $$ The vector field $ \v $ and the scalar functions $ a $, $ \alpha $, $ f $, $ u_0 $, and $ g $ may vary with the spatial coordinate $ \x $ and must be known.

Such a second-order PDE needs exactly one boundary condition at each point of the boundary, so $ \partial\Omega_N\cup\partial\Omega_D $ must be the complete boundary $ \partial\Omega $.

Assume that the boundary function $ u_0(\x) $ is defined for all $ \x\in\Omega $. The unknown function can then be expanded as

$$ u = B + \sum_{j\in\If} c_j\baspsi_j,\quad B = u_0 \tp $$ The variational formula is obtained from Galerkin's method, which technically implies multiplying the PDE by a test function $ v $ and integrating over $ \Omega $:

$$ \int_{\Omega} (\v\cdot\nabla u + \alpha u)v\dx = \int_{\Omega} \nabla\cdot\left( a\nabla u\right)\dx + \int_{\Omega}fv \dx \tp $$ The second-order term is integrated by parts, according to

$$ \int_{\Omega} \nabla\cdot\left( a\nabla u\right)v \dx = -\int_{\Omega} a\nabla u\cdot\nabla v\dx + \int_{\partial\Omega} a\frac{\partial u}{\partial n} v\ds \tp $$ The variational form now reads

$$ \int_{\Omega} (\v\cdot\nabla u + \alpha u)v\dx = -\int_{\Omega} a\nabla u\cdot\nabla v\dx + \int_{\partial\Omega} a\frac{\partial u}{\partial n} v\ds + \int_{\Omega} fv \dx \tp $$ The boundary term can be developed further by noticing that $ v\neq 0 $ only on $ \partial\Omega_N $,

$$ \int_{\partial\Omega} a\frac{\partial u}{\partial n} v\ds = \int_{\partial\Omega_N} a\frac{\partial u}{\partial n} v\ds, $$ and that on $ \partial\Omega_N $, we have the condition $ a\frac{\partial u}{\partial n}=-g $, so the term becomes

$$ \int_{\Omega} (\v\cdot\nabla u + \alpha u)v\dx = -\int_{\Omega} a\nabla u\cdot\nabla v \dx - \int_{\partial\Omega_N} g v\ds + \int_{\Omega} fv \dx \tp $$

$$ (\v\cdot\nabla u, v) + (\alpha u,v) = - (a\nabla u,\nabla v) - (g,v)_{N} + (f,v) \tp $$ The subscript $ \,{}_N $ in $ (g,v)_{N} $ is a notation for a line or surface integral over $ \partial\Omega_N $.

$$ \begin{align*} \sum_{j\in\If} ((\v\cdot\nabla \baspsi_j, \baspsi_i) &+ (\alpha \baspsi_j ,\baspsi_i) + (a\nabla \baspsi_j,\nabla \baspsi_i))c_j = \\ & (g,\baspsi_i)_{N} + (f,\baspsi_i) - (\v\cdot\nabla u_0, \baspsi_i) + (\alpha u_0 ,\baspsi_i) + (a\nabla u_0,\nabla \baspsi_i) \tp \end{align*} $$ This is a linear system with matrix entries

$$ A_{i,j} = (\v\cdot\nabla \baspsi_j, \baspsi_i) + (\alpha \baspsi_j ,\baspsi_i) + (a\nabla \baspsi_j,\nabla \baspsi_i) $$ and right-hand side entries

$$ b_i = (g,\baspsi_i)_{N} + (f,\baspsi_i) - (\v\cdot\nabla u_0, \baspsi_i) + (\alpha u_0 ,\baspsi_i) + (a\nabla u_0,\nabla \baspsi_i), $$ for $ i,j\in\If $.

In the finite element method, we usually express $ u_0 $ in terms of basis functions and restrict $ i $ and $ j $ to run over the degrees of freedom that are not prescribed as Dirichlet conditions. However, we can also keep all the $ c_j $, $ j\in\If $, as unknowns drop the $ u_0 $ in the expansion for $ u $, and incorporate all the known $ c_j $ values in the linear system. This has been explained in detail in the 1D case.

Transformation to a reference cell in 2D and 3D

$$ \begin{equation} \int_{{\Omega}^{(e)}} a(\x)\nabla\basphi_i\cdot\nabla\basphi_j\dx, \end{equation} $$ where the $ \basphi_i $ functions are finite element basis functions in 2D or 3D, defined in the physical domain. Suppose we want to calculate this integral over a reference cell, denoted by $ \tilde\Omega^r $, in a coordinate system with coordinates $ \X = (X_0, X_1) $ (2D) or $ \X = (X_0, X_1, X_2) $ (3D). The mapping between a point $ \X $ in the reference coordinate system and the corresponding point $ \x $ in the physical coordinate system is given by a vector relation $ \x(\X) $. The corresponding Jacobian, $ J $, of this mapping has entries $$ J_{i,j}=\frac{\partial x_j}{\partial X_i}\tp $$

The change of variables requires $ \dx $ to be replaced by $ \det J\dX $. The derivatives in the $ \nabla $ operator in the variational form are with respect to $ \x $, which we may denote by $ \nabla_{\x} $. The $ \basphi_i(\x) $ functions in the integral are replaced by local basis functions $ \refphi_r(\X) $ so the integral features $ \nabla_{\x}\refphi_r(\X) $. We readily have $ \nabla_{\X}\refphi_r(\X) $ from formulas for the basis functions in the reference cell, but the desired quantity $ \nabla_{\x}\refphi_r(\X) $ requires some efforts to compute. All the details are provided below.

$$ \frac{\partial \refphi_r}{\partial X} = \frac{\partial \basphi_i}{\partial X} = \frac{\partial \basphi_i}{\partial x}\frac{\partial x}{\partial X} + \frac{\partial \basphi_i}{\partial y}\frac{\partial y}{\partial X}, $$ and $$ \frac{\partial \refphi_r}{\partial Y} = \frac{\partial \basphi_i}{\partial Y} = \frac{\partial \basphi_i}{\partial x}\frac{\partial x}{\partial Y} + \frac{\partial \basphi_i}{\partial y}\frac{\partial y}{\partial Y} \tp $$ We can write these two equations as a vector equation $$ \left[\begin{array}{c} \frac{\partial \refphi_r}{\partial X}\\ \frac{\partial \refphi_r}{\partial Y} \end{array}\right] = \left[\begin{array}{cc} \frac{\partial x}{\partial X} & \frac{\partial y}{\partial X}\\ \frac{\partial x}{\partial Y} & \frac{\partial y}{\partial Y} \end{array}\right] \left[\begin{array}{c} \frac{\partial \basphi_i}{\partial x}\\ \frac{\partial \basphi_i}{\partial y} \end{array}\right] $$ Identifying $$ \nabla_{\X}\refphi_r = \left[\begin{array}{c} \frac{\partial \refphi_r}{\partial X}\\ \frac{\partial \refphi_r}{\partial Y} \end{array}\right], \quad J = \left[\begin{array}{cc} \frac{\partial x}{\partial X} & \frac{\partial y}{\partial X}\\ \frac{\partial x}{\partial Y} & \frac{\partial y}{\partial Y} \end{array}\right], \quad \nabla_{\x}\basphi_r = \left[\begin{array}{c} \frac{\partial \basphi_i}{\partial x}\\ \frac{\partial \basphi_i}{\partial y} \end{array}\right], $$ we have the relation

$$ \nabla_{\X}\refphi_r = J\cdot\nabla_{\x}\basphi_i,$$ which we can solve with respect to $ \nabla_{\x}\basphi_i $:

$$ \begin{equation} \nabla_{\x}\basphi_i = J^{-1}\cdot\nabla_{\X}\refphi_r\tp \end{equation} $$ On the reference cell, $ \basphi_i(\x) = \refphi_r(\X) $, so

$$ \begin{equation} \nabla_{\x}\refphi_r(\X) = J^{-1}(\X)\cdot\nabla_{\X}\refphi_r(\X)\tp \end{equation} $$

This means that we have the following transformation of the integral in the physical domain to its counterpart over the reference cell:

$$ \begin{equation} \int_{\Omega}^{(e)} a(\x)\nabla_{\x}\basphi_i\cdot\nabla_{\x}\basphi_j\dx \int_{\tilde\Omega^r} a(\x(\X))(J^{-1}\cdot\nabla_{\X}\refphi_r)\cdot (J^{-1}\cdot\nabla\refphi_s)\det J\dX \end{equation} $$

Numerical integration

Integrals are normally computed by numerical integration rules. For multi-dimensional cells, various families of rules exist. All of them are similar to what is shown in 1D: $ \int f \dx\approx \sum_jw_if(\x_j) $, where $ w_j $ are weights and $ \x_j $ are corresponding points.

The file numint.py contains the functions quadrature_for_triangles(n) and quadrature_for_tetrahedra(n), which returns lists of points and weights corresponding to integration rules with n points over the reference triangle with vertices $ (0,0) $, $ (1,0) $, $ (0,1) $, and the reference tetrahedron with vertices $ (0,0,0) $, $ (1,0,0) $, $ (0,1,0) $, $ (0,0,1) $, respectively. For example, the first two rules for integration over a triangle have 1 and 3 points:

>>> import numint
>>> x, w = numint.quadrature_for_triangles(num_points=1)
>>> x
[(0.3333333333333333, 0.3333333333333333)]
>>> w
[0.5]
>>> x, w = numint.quadrature_for_triangles(num_points=3)
>>> x
[(0.16666666666666666, 0.16666666666666666),
 (0.66666666666666666, 0.16666666666666666),
 (0.16666666666666666, 0.66666666666666666)]
>>> w
[0.16666666666666666, 0.16666666666666666, 0.16666666666666666]

Rules with 1, 3, 4, and 7 points over the triangle will exactly integrate polynomials of degree 1, 2, 3, and 4, respectively. In 3D, rules with 1, 4, 5, and 11 points over the tetrahedron will exactly integrate polynomials of degree 1, 2, 3, and 4, respectively.

Convenient formulas for P1 elements in 2D

We shall now provide some formulas for piecewise linear $ \basphi_i $ functions and their integrals in the physical coordinate system. These formulas make it convenient to compute with P1 elements without the need to work in the reference coordinate system and deal with mappings and Jacobians. A lot of computational and algorithmic details are hidden by this approach.

Let $ \Omega^{(e)} $ be cell number $ e $, and let the three vertices have global vertex numbers $ I $, $ J $, and $ K $. The corresponding coordinates are $ (\xno{I},\yno{I}) $, $ (\xno{J},\yno{J}) $, and $ (\xno{K},\yno{K}) $. The basis function $ \basphi_I $ over $ \Omega^{(e)} $ have the explicit formula

$$ \begin{equation} \basphi_I (x,y) = \half\Delta \left( \alpha_I + \beta_Ix + \gamma_Iy\right), \tag{96} \end{equation} $$ where

$$ \begin{align} \alpha_I &= \xno{J}\yno{K} - \xno{K}\yno{J}, \tag{97}\\ \beta_I &= \yno{J} - \yno{K}, \tag{98}\\ \gamma_I &= \xno{K} - \xno{J}, \tag{99}\\ 2\Delta &= \det\left(\begin{array}{rrr} 1 & \xno{I} & \yno{I} \\ 1 & \xno{J} & \yno{J} \\ 1 & \xno{K} & \yno{K} \end{array}\right) \tp \tag{100} \end{align} $$ The quantity $ \Delta $ is the area of the cell.

The following formula is often convenient when computing element matrices and vectors:

$$ \begin{equation} \int_{\Omega^{(e)}} \basphi_I^{p}\basphi_J^{q}\basphi_K^{r} dx dy = {p!q!r!\over (p+q+r+2)!}2\Delta \tag{101} \tp \end{equation} $$ (Note that the $ q $ in this formula is not to be mixed with the $ q(e,r) $ mapping of degrees of freedom.)

As an example, the element matrix entry $ \int_{\Omega^{(e)}} \basphi_I\basphi_J\dx $ can be computed by setting $ p=q=1 $ and $ r=0 $, when $ I\neq J $, yielding $ \Delta/12 $, and $ p=2 $ and $ q=r=0 $, when $ I=J $, resulting in $ \Delta/6 $. We collect these numbers in a local element matrix:

$$ \frac{\Delta}{12} \left[\begin{array}{ccc} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2 \end{array}\right] $$

The common element matrix entry $ \int_{\Omega^{(e)}} \nabla\basphi_I\cdot\nabla\basphi_J\dx $, arising from a Laplace term $ \nabla^2u $, can also easily be computed by the formulas above. We have

$$ \nabla\basphi_I\cdot\nabla\basphi_J = \frac{\Delta^2}{4}(\beta_I\beta_J + \gamma_I\gamma_J) = \hbox{const},$$ so that the element matrix entry becomes $ \frac{1}{4}\Delta^3(\beta_I\beta_J + \gamma_I\gamma_J) $.

From an implementational point of view, one will work with local vertex numbers $ r=0,1,2 $, parameterize the coefficients in the basis functions by $ r $, and look up vertex coordinates through $ q(e,r) $.

Variational formulations in 2D and 3D

Transformation to a reference cell in 2D and 3D

Numerical integration

Convenient formulas for P1 elements in 2D