Dealing with volumic curved elements¶

The maths¶

To be improved :

The mapping¶

Consider a reference volumic element $\hat{\Omega}$ or order $n$ and the $n$ associated shape functions $\hat{\lambda}_i$. We restrict this explanation to elements of order $n$ with $n$ nodes $\hat{x}_i$ (excluding Hermite elements for instance).

Now considering a local volumic element of order $n$ and its $n$ nodes $x_i$. The reference to local element mapping read : $$x = F_{rl}(\hat{x}) = \sum_{i=1}^n \hat{\lambda}_i(\hat{x}) x_i.$$

So the shape functions in the reference element and the nodes of the local element give the mapping.

Integration formula¶

Using Finite Element Method, one has to compute integrals in the local element : $\int_\Omega f(x) dx$. To compute an integral on an element, we often use numerical integration by applying quadrature rules: $$\int_\Omega f(x) dx = \sum_{i \in N} \omega_i f(\xi_i),$$ where $\omega_i$ and $\xi_i$ are the $N$ weights and nodes of the selected quadrature rule.

Now, the quadrature rules are not known for every type of geometrical elements, they are actually often determined for simple geometrical elements called reference element : the $[-1,1]$ segment, a right triangle, a square…

To compute an integral on a arbitrary element, we have to apply a transformation to this local element to obtain the reference element and apply the quadrature rule. To do so, we need this very important formula: $$\int_\Omega f(x) dx = \int_\hat{\Omega} |J_{rl}(\hat{x})| \left(f \circ F_{rl}\right)(\hat{x}) d\hat{x},$$ where $F_{rl}$ is the reference to local mapping, i.e $F_{rl}(\hat{x}) = s$, and $J_{rl}$ is the determinant of the jacobian matrix of this mapping.

Summary¶

To use a curved volumic element, we only need to know a reference element and the associated shape functions. Then the mapping is directly obtained and the integration formula can be applied.

The code¶

Import necessary modules: LinearAlgebra, ForwardDiff for automatic differentiation, Cubature to apply quadrature on ‘multi-dimensionnal’ elements (otherwise in 1D there is FastGaussQuadrature).

We overload the sum and product operator of two functions (in order for ForwardDiff to apply without harm):

We create an alias for the sum of functions, just for fun:

Now we define two structures for nodes and elements. These structures could be parametrized by the size of the dimensionnal space ($\mathbb{R}^2$, $\mathbb{R}^3$…), but for the sake of clarity everything will be in $\mathbb{R}^2$ here. A Node is simply an array containing the coordinates. An Element is an ordered list of nodes. The chosen order is arbitrary, here we will choose gmsh convention (more details below).

A convenient function (specially to plot) is immediatly declared to convert an Element into a matrices containing all its nodes coordinates:

We can now specify two particular entities : the quadrilateral of order 2, defined by its nine nodes; and the reference quadrilateral (=square) of order 2. The order of the nodes is conform to gmsh convention : the four nodes of the square, the four mid-nodes and the center. The choice of $[-1,1] \times [-1,1]$ is arbitrary, the only thing that matter is that the shape-functions (see below) are built on these nodes.

Let’s plot the nodes and edges of our reference element:

Now the nine shape functions associated to the reference element or order 2. These are Lagrange $x,y$ polynomials who cancel in every nodes but one, where their value is 1. Important detail here : because we will apply automatic differentiation on a combination of these functions, we use a vector variable x instead of a couple x,y.

With these shape functions, the reference element to local element mapping can be easily defined. Indeed this mapping is: $$ \begin{pmatrix} x \\ y \end{pmatrix} = F_{rl}(\hat{x}, \hat{y}) = \sum_{i=1}^n \hat{\lambda}_i(\hat{x}_i, \hat{y}_i) \begin{pmatrix} x_i\\ y_i \end{pmatrix}$$ where $n$ is the number of nodes.

The building the mapping only consists in gathering the shape functions and computing the sum weighted by the coordinates of the local element:

Let’s illustrate our mapping. First we define 9 nodes of our curved element. This element represent a portion of annulus defined by $1 \leq r \leq 2$ and $0 \leq \theta \leq \pi / 4$. This can be obtained with gmsh : define your geometry, then ask for a second order mesh and read the output nodes.

Then we can build a quadrilateral of order 2 and retrieve the mapping:

Now plot the local element using the mapping and the nodes. We draw the edges of the local element by applying the mapping on the edge of the reference element. We also « translate » the local element to improve the plot.

Now that we have the mapping, we juste need a function to compute an integral in the local element : first we apply the local to reference element transformation (see formula in the first part of this document) and then we compute the integral in the reference element with a quadrature (using the Cubature package):

Here we chose to provide the integral function with the mapping and the jacobian, but it could have been retrieved from the element. This jacobian is simply computed from the mapping $F_{rl}$ using automatic differentiation with the ForwardDiff package (see below).

Let’s play with this function. First, what is the area of our reference element?

What is the area of the annulus portion? Analytically, the calculus is $\mathrm{A} = \int_{\theta = 0}^{\pi/4} \int_{r=1}^2 r dr d\theta = 3 \pi / 8 $. In cartesian coordinates, it is just the integral of the $x \rightarrow 1$ function.

Note that the two results are close but not exactly the same : this is likely due to the approximation of the annulus section by a second order quadrilateral.

And finally, let’s integrate a more complex function, such as $(x,y) \rightarrow \sqrt{x^2 + y^2}$. Analytically, we have : $\int_\Omega \sqrt{x^2 + y^2} dxdy = \int_{\theta = 0}^{\pi/4} \int_{r=1}^2 r^2 dr d\theta = 7\pi / 4$.