Surfacic PDEs on curved elements¶

In this document, we show how to solve a surfacic partial differential equation (PDE) on curved geometrical element. The treated example will be a linear transport equation on a parabola : $$ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial s} = 0,$$ where $s$ is the curvilinear abscissa on the parabola.

Note that dealing with curved elements when integrating volumic PDEs is less hard than integrating surfacic PDEs. This is due to the fact that a surface in a 3D space is a hyperplane. The same applies for a line in a 2D space (our example here).

The maths¶

Temporal discretization¶

The choice for the temporal discretization has really no direct influence on the problem we are treating. I mean that whatever the temporal discretization you select, it does not influence how you treat the curvature of your geometric domain.

We choose an explicit backward scheme. The discrete equation reads: $$\frac{u(s,t^{n+1}) – u(s,t^n)}{\Delta t} + c \frac{\partial u(s,t^n)}{\partial s} = 0.$$

From now, $u(s,t^n)$ will be noted $u^n(s)$.

Spatial discretization¶

The spatial discretization will be accomplished with the Finite Element Method (FEM). We won’t discuss the details of the procedure here. I’ll be really quick.

Let’s compute the weak form of our problem, computing the scalar product of our equation by a test function $v$ on the geometrical (curved!) domain $\Omega$: $$\int_\Omega u^{n+1}(s) v(s) ds = \int_\Omega u^{n}(s) v(s) ds – c \Delta t \int_\Omega \frac{\partial u^{n}(s)}{\partial s} v(s) ds.$$

Now, interpolating $u$ in each element our domain by known shape functions we have : $ u^n(s) = \sum_j u_j^n \lambda(s)$, where $u_j^n$ is the value of $u^n$ on node $j$. Using $v = \lambda_i$, we obtain: $$\sum_j u_j^{n+1} \int_\Omega \lambda_i(s) \lambda_j(s) ds = \sum_j u_j^n \int_\Omega \lambda_i(s) \lambda_j(s) ds – c \Delta t \sum_j u_j^n \int_\Omega \lambda_i(s) \frac{\partial \lambda_j}{\partial s}(s) ds$$

Calling $M$ the matrix whose elements are $\lambda_i(s) \lambda_j(s) ds$ and $D$ the matrix whose elements are $\int_s \lambda_i(s) \frac{\partial \lambda_j}{\partial s}(s) ds$, we have: $$M \mathbf{u}^{n+1} = M \mathbf{u}^n – c \Delta t D \mathbf{u}^n ~~~\rightarrow~~~ \mathbf{u}^{n+1} = (I – c \Delta t M^{-1} D) \mathbf{u}^n $$

If we set $A = I – c \Delta t M^{-1} D$, and add boundary conditions in it, we finally have: $$ \mathbf{u}^n = A^n \mathbf{u}^0$$

So all we have to do to solve our surfacic PDE is to be able to compute the integrals in $M$ and $D$ on the curved domain $\Omega$.

Integration formula¶

Since $\lambda_i$ are often polynomials, computing the integrals analytically could be possible. However, as soon as more complicated function multiplies these polynomials, this is harded. Imagine that instead of a constant transport velocity $c$ we had chosen a sinusoidal velocity $ c(s) = \sin(s) $. Then we would have to integrate $\int \sin(s) \lambda_i \partial_s \lambda_j$, wich is way harder!

For this reason, we use numerical integration applying quadrature rules: $$\int_\hat{\Omega} f(\xi_i) 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, for instance our curved 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(s) ds = \int_\hat{\Omega} |J_{rl}(\xi)| \left(f \circ F_{rl}\right)(\xi) d\xi,$$ where $F_{rl}$ is the reference to local mapping, i.e $F_{rl}(\xi) = s$, and $J_{rl}$ is the determinant of the jacobian matrix of this mapping.

Reference to local mapping for a parabola¶

Here we describe how to determine the reference element to local element mapping in the case of a parabola. The parabola will be defined by three nodes : its two boundary nodes, $p_1$ and $p_2$, and an additionnal node $p_3$ lying on the parabola. The first step is to map any parabola, whose three nodes are defined by (x,y) in the global frame, to a parabola where the two extremities lie on a local « $\xi$-axis ». So we take a parabola, we build a local frame whose $\xi$-axis is defined by the two parabola extremities and whose $\eta$-axis define a plane (with local $\xi$-axis) containing the third parabola point. More over we scale the axis to obtain $\xi(p_1) = -1$ and $\xi(p_2) = 1$.

We could call this operation the global to local mapping.

Once we have a mapped parabola defined by $\eta = a (\xi-1)(\xi+1)$, ($a$ depends on $p_3$) we can more easily build the $F_{rl}(\xi)=s$ mapping function. Indeed, for a curve defined by $\eta(\xi)$ in an orthonormal frame, the curvilinear is defined by: $$s(\xi) = \int_{s(p_1)}^{\xi} \sqrt{1 + \eta’^2(t)} dt$$

For our parabola, $\eta’^2(\xi) = 4 a^2 \xi^2 $. So if we decide that $s(p_1) = 0$, we have: $$s(\xi) = \int_0^{\xi} \sqrt{1 + 4 a^2 t^2} dt$$

I don’t detail the maths but we obtain: $$s(\xi) = \frac{1}{2a} \left( S(2a\xi) – S(-2a) \right), $$ where $S(\xi) = \xi \sqrt{1 + \xi^2} + \sinh^{-1}(\xi)$

To conclude, don’t forget that we scaled the local $\xi$-axis by $2 / l$ where $l$ is the distance between $p_1$ and $p_2$. So the final result is: $$ F_{rl}(\xi) = \frac{l}{4a}\left( S(2a\xi) – S(-2a) \right)$$.

Open discussion : interpolating in the reference element or in the local element?¶

This part is unclear, skip it if you are not sure to understand what I mean.

We said above that we interpolate $u$ with shape functions : $u(s) = \sum_i u_i \lambda_i(s)$. This being said, we can either choose to define $\lambda_i(s)$ explicitely in the local element (with a Lagrange polynomial for instance), or we can choose to define $\hat{\lambda}_i(\xi)$ explicitely, where $\hat{\lambda}$ is the shape function in the reference element. The relation between the two is simply : $\hat{\lambda} = \lambda \circ F_{rl}$, and equivalently $\lambda = \hat{\lambda} \circ F_{lr}$.

But note that if $\hat{\lambda}$ is a polynomial, $\lambda$ is not necessarily of the same order (if $F_{lr}(s) = \xi^2$ for instance), and it is not even necessarily a polynomial at all. Hence saying « I use a n-order spatial discretization » is ambiguous : it can mean a n-order in the reference element or in the local element.

In the next part, you will see that I had to choose to use Lagrange polynomials in the local element because I wasn’t able to determine $F_{lr}$. And if you don’t know $F_{lr}$, you can’t compute $D$ because $\partial_s \lambda = \partial_s F_{lr} [(\partial_s \hat{\lambda}) \circ F_{lr}]$ so you need at least to know $\partial_s F_{lr}$.

The code¶

Now, let’s build a demo example of what we described, solving our linear transport equation on a parabola. We will be using the Julia language. If you skipped the first part, I will remind some definitions here so it might be a bit redundant.

Let’s start by some usefull import:

We also need to define the sum and the product of two functions since we will use this later:

The curved element will be a parabola, defined by three nodes : its two boundary nodes, $p_1$ and $p_2$, and an additionnal node $p_3$ lying on the parabola. Note that the parabola is here aligned with the x-axis but it doesn’t necessarily has to.

Now the first thing to build is the reference element to local element mapping. The definition has been given if the first part of the document. We start by defining the « primitive » function, called $S$ in the previous section:

Now we can define the mapping.

Note that the function requires that the three nodes are not aligned. But it could be easily modified to handle this configuration (the mapping is then way easier to define).

Now that we have the mapping, we can compute it’s derivative (hence the jacobian since the element is mono-dimensionnal). But we will show this later.

Before that, let’s implement method to compute an integral, i.e applying a quadrature rule:

In fact the for loop of this function could be simply replaced by dot(weights, f(nodes)) but to keep this document as simple as possible, all the functions defined here are scalar only.

Rembember the integration formula? Yes, this one: $$\int_\Omega g(s) ds = \int_\hat{\Omega} |J_{rl}(\xi)| \left(g \circ F_{rl}\right)(\xi) d\xi.$$

Let’s define an integration function that applies this formula, knowing the mapping F and the absolute value of the determinant ot its jacobian matrix J:

It is this simple !

Let’s have a break and play with this things. We can for instance compute the length of our parabola by integrating the function $s \rightarrow 1$ (you can compare the result using a trapz function on the parabola with its polynomial definition):

You are not convinced? Let’s integrate the function $s \rightarrow s$ along the parabola, i.e $\int_\Omega s ds$ and compare with the analytical result:

Note that the precision of the result depend on the quadrature order.

We have now the capacity to integrate any function on our parabola. So let’s define shape functions in the local element to implement the FEM. We will be using Lagrange polynomials of order 1:

Now we can write a function that returns the $A$ matrix (without boundary condition) of the transport equation we want to solve. We defined in the first part of this document:

And now the final demonstration. Consider the linear transport equation $\partial_t u + \partial_s u = 0$. If you input $ u = u_{in}$ on a boundary of your domain, this information will be transported along your element at $1m/s$ (since $c=1$). So if your element measures $3$ meters, then after $3$ seconds the value of $u$ will be $u_{in}$ everywhere on you element. And this is true for curved elements.

The first example is for a « flat » element, a line, of length $1$ meter : the segment $[0, 1]$. Let’s input $666$ for $u$ and check that after one second, the information has reached $x=1$:

Ok but what about the parabola we’ve draw at the beginning? Let’s try: