Skip to content

Finite Element

As stated in tutorials, FEM approximates all functions w as with finite element basis functions \phi_k(x,y) and numbers w_k (k=0,\cdots,M-1). The functions \phi_k(x,y) are constructed from the triangle T_{i_k}, and called shape functions.

In FreeFem++ the finite element space is easily created by :

1
fespace IDspace(IDmesh,<IDFE>);

or with \ell pairs of periodic boundary conditions in 2D :

1
2
3
4
fespace IDspace(IDmesh,<IDFE>,
    periodic=[[la1, sa1], [lb1, sb1],
              ...
              [lak, sak], [lbk, sbl]]);

and in 3D :

1
2
3
4
fespace IDspace(IDmesh,<IDFE>,
    periodic=[[la1, sa1, ta1], [lb1, sb1, tb1],
              ...
              [lak, sak, tak], [lbk, sbl, tbl]]);

where IDspace is the name of the space (e.g. Vh), IDmesh is the name of the associated mesh and <IDFE> is an identifier of finite element type.

In 2D we have a pair of periodic boundary conditions, if [la_i, sa_i], [lb_i, sb_i] is a pair of int, and the 2 labels la_i and lb_i refer to 2 pieces of boundary to be in equivalence. If [la_i, sa_i], [lb_i, sb_i] is a pair of real, then sa_i and sb_i give two common abscissa on the two boundary curves, and two points are identified as one if the two abscissa are equal.

In 2D, we have a pair of periodic boundary conditions, if [la_i, sa_i, ta_i], [lb_i, sb_i, tb_i] is a pair of int, the 2 labels la_i and lb_i define the 2 pieces of boundary to be in equivalence. If [la_i, sa_i, ta_i], [lb_i, sb_i, tb_i] is a pair of real, then sa_i, ta_i and sb_i, tb_i give two common parameters on the two boundary surfaces, and two points are identified as one if the two parameters are equal.

Note

The 2D mesh of the two identified borders must be the same, so to be sure, use the parameter fixedborder=true in buildmesh command (see fixedborder).

As of today, the known types of finite elements are:

  • [P0, P03d] piecewise constant discontinuous finite element (2d, 3d), the degrees of freedom are the barycenter element value.

  • [P1, P13d] piecewise linear continuous finite element (2d, 3d), the degrees of freedom are the vertices values.

  • [P1dc] piecewise linear discontinuous finite element

    Warning

    Due to an interpolation problem, the degree of freedom is not the vertices but three vertices which move inside T(X)= G + .99 (X-G) where G is the barycenter.

  • [P1b, P1b3d] piecewise linear continuous finite element plus bubble (2d, 3d)

    The 2D Case:

    The 3D Case: where \lambda^{K}_{i}, i=0,..,d are the d+1 barycentric coordinate functions of the element K (triangle or tetrahedron).

  • P1bl,P1bl3d piecewise linear continuous finite element plus linear bubble (2d, 3d). The bubble is built by splitting the K, a barycenter in d+1 sub element. (need load "Element_P1bl")

  • [P2, P23d] piecewise P_{2} continuous finite element (2d, 3d) where P_{2} is the set of polynomials of \R^{2} of degrees \le 2.

  • [P2b] piecewise P_{2} continuous finite element plus bubble

  • [P2dc] piecewise P_{2} discontinuous finite element

    Warning

    Due to an interpolation problem, the degree of freedom is not the six P2 nodes but six nodes which move inside T(X)= G + .99 (X-G) where G is the barycenter.

  • [P2h] quadratic homogeneous continuous (without P1).

  • [P3] piecewise P_{3} continuous finite element (2d) (needs load "Element_P3") where P_{3} is the set of polynomials of \R^{2} of degrees \le 3.

  • [P3dc] piecewise P_{3} discontinuous finite element (2d) (needs load "Element_P3dc") where P_{3} is the set of polynomials of \R^{2} of degrees \le 3.

  • [P4] piecewise P_{4} continuous finite element (2d) (needs load "Element_P4") where P_{4} is the set of polynomials of \R^{2} of degrees \le 4.

  • [P4dc] piecewise P_{4} discontinuous finite element (2d) (needs load "Element_P4dc") where P_{4} is the set of polynomials of \R^{2} of degrees \le 3.

  • [P0Edge] piecewise P_{0} discontinuous finite element (2d) contained on each edge of the mesh.

  • [P1Edge] piecewise P_{1} discontinuous finite element (2d) (needs load "Element_PkEdge") P_1 on each edge of the mesh.

  • [P2Edge] piecewise P_{2} discontinuous finite element (2d) (needs load "Element_PkEdge") P_2 on each edge of the mesh.

  • [P3Edge] piecewise P_{3} discontinuous finite element (2d) (needs load "Element_PkEdge") P_3 on each edge of the mesh.

  • [P4Edge] piecewise P_{4} discontinuous finite element (2d) (needs load "Element_PkEdge") P_4 on each edge of the mesh.

  • [P5Edge] piecewise P_{5} discontinuous finite element (2d) (needs load "Element_PkEdge") P_5 on each edge of the mesh.

  • [P2Morley] piecewise P_{2} non conform finite element (2d) (needs load "Morley") where P_{2} is the set of polynomials of \R^{2} of degrees \le 2.

    Warning

    To build the interplant of a function u (scalar) for this finite element, we need the function and 2 partial derivatives (u,u_x, u_y), creating this vectorial finite element with 3 components (u,u_x,u_y).

    See our example for solving the BiLaplacien problem:

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    load "Morley"
    
    // Parameters
    int nn = 10;
    real h = 0.01;
    
    real f = 1;
    
    // Mesh
    mesh Th = square(nn, nn);
    Th = adaptmesh(Th, h, IsMetric=1);
    
    // Fespace
    fespace Vh(Th, P2Morley); //The Morley finite element space
    Vh [u, ux, uy], [v, vx, vy];
    
    // Macro
    macro bilaplacien(u, v) (dxx(u)*dxx(v) + dyy(u)*dyy(v) + 2.*dxy(u)*dxy(v)) //
    
    // Problem
    solve bilap ([u, ux, uy], [v, vx, vy])
        = int2d(Th)(
              bilaplacien(u, v)
        )
        - int2d(Th)(
              f*v
        )
        + on(1, 2, 3, 4, u=0, ux=0, uy=0)
        ;
    
    // Plot
    plot(u, cmm="u");
    
  • [HCT] P_3 C^1 conforms finite element (2d) (needs load "Element_HCT") one 3 sub triangles.

    Lets call \mathcal{T}^\triangle_{h} the sub mesh of \mathcal{T}_{h} where all triangles are split in 3 at the barycenter. where P_{3} is the set of polynomials of \R^{2} of degrees \le 3. The degrees of freedom are the values of the normal derivative at the mid-point of each edge BERNADOU1980.

    Warning

    To build the interplant of a function u (scalar) for this finite element, we need the function and 2 partial derivatives (u,u_x, u_y), creating this vectorial finite element with 3 components (u,u_x,u_y) like in previous Finite Element.

  • [P2BR] (needs load "BernadiRaugel") the Bernadi Raugel Finite Element is a Vectorial element (2d) with 2 components, see BERNARDI1985. It is a 2D coupled Finite Element, where the Polynomial space is P_1^2 with 3 normal bubble edge functions (P_2). There are 9 degrees of freedom:

    • 2 components at each of the 3 vertices and
    • the 3 flux on the 3 edges.
  • [RT0, RT03d] Raviart-Thomas finite element of degree 0.

    The 2D Case: The 3D Case: where by writing \textrm{div }\mathbf{w}=\sum_{i=1}^d\p w_i/\p x_i with \mathbf{w}=(w_i)_{i=1}^d: and where \alpha^1_{K}, \alpha^2_{K}, \alpha^3_{K}, \beta_{K} are real numbers.

  • [RT0Ortho] Raviart-Thomas Orthogonal, or Nedelec finite element type I of degree 0 in dimension 2

  • [Edge03d] 3d Nedelec finite element or Edge Element of degree 0.

    The 3D Case: where by writing \textrm{curl}\mathbf{w}=\vectthree{\p w_2/\p x_3-\p w_3/\p x_2}{\p w_3/\p x_1-\p w_1/\p x_3}{\p w_1/\p x_2-\p w_2/\p x_1} with \mathbf{w}=(w_i)_{i=1}^d: and \alpha^1_{K},\alpha^2_{K},\alpha^3_{K},\beta^1_{K},\beta^2_{K},\beta^3_{K} are real numbers.

  • [Edge13d] (needs load "Element_Mixte3d") 3d Nedelec finite element or Edge Element of degree 1.

  • [Edge23d] (needs load "Element_Mixte3d") 3d Nedelec finite element or Edge Element of degree 2.

  • [P1nc] piecewise linear element continuous at the mid-point of the edge only in 2D (Crouzeix-Raviart Finite Element 2D).

  • [P2pnc] piecewise quadratic plus a P3 bubble element with the continuity of the 2 moments on each edge (needs load "Element_P2pnc")

  • [RT1] (needs load "Element_Mixte")

  • [RT1Ortho] (needs load "Element_Mixte")

  • [RT2] (needs load "Element_Mixte")

  • [RT2Ortho] (needs load "Element_Mixte")

  • [BDM1] (needs load "Element_Mixte") the Brezzi-Douglas-Marini finite element

  • [BDM1Ortho] (needs load "Element_Mixte") the Brezzi-Douglas-Marini Orthogonal also call Nedelec of type II , finite element

  • [FEQF] (needs load "Element_QF") the finite element to store functions at default quadrature points (so the quadrature is qf5pT in 2D and is qfV5 in 3d). For over quadrature you have the following corresponding finite element's quadrature formula.

    • FEQF1 \mapsto qf1pT,
    • FEQF2 \mapsto qf2pT,
    • FEQF5 \mapsto qf5pT,
    • FEQF7 \mapsto qf7pT,
    • FEQF9 \mapsto qf9pT,
    • FEQF13d \mapsto qfV1,
    • FEQF23d \mapsto qfV2,
    • FEQF53d \mapsto qfV5

You can use this element to optimize the storage and reuse of functions with a long formula inside an integral for non linear processes.

Use of fespace in 2D#

With the 2D finite element spaces

X_{h} = \left\{ v \in H^{1}(]0,1[^2) |\; \forall K \in \mathcal{T}_{h}\quad v_{|K} \in P_{1} \right\}
X_{ph} = \left\{ v \in X_{h} |\; v\left(\vecttwo{0}{.}\right) = v\left(\vecttwo{1}{.}\right) , v\left(\vecttwo{.}{0}\right) = v\left(\vecttwo{.}{1}\right) \right\}
M_{h} = \left\{ v \in H^{1}(]0,1[^2) |\; \forall K \in \mathcal{T}_{h}\quad v_{|K} \in P_{2} \right\}
R_{h} = \left\{ \mathbf{v} \in H^{1}(]0,1[^2)^{2} |\; \forall K \in \mathcal{T}_{h}\quad \mathbf{v}_{|K}(x,y) = \vecttwo{\alpha_{K}}{\beta_{K}} + \gamma_{K}\vecttwo{x}{y} \right\}

when \mathcal{T}_h is a mesh 10\times 10 of the unit square ]0,1[^2, we only write in FreeFem++ :

1
2
3
4
5
6
mesh Th = square(10, 10);
fespace Xh(Th, P1); //scalar FE
fespace Xph(Th,P1,
            periodic=[[2, y], [4, y], [1, x], [3, x]]); //bi-periodic FE
fespace Mh(Th, P2); //scalar FE
fespace Rh(Th, RT0); //vectorial FE

where Xh, Mh, Rh expresses finite element spaces (called FE spaces) X_h,\, M_h,\, R_h, respectively.

To use FE-functions u_{h},v_{h} \in X_{h}, p_{h},q_{h} \in M_{h} and U_{h},V_{h} \in R_{h}, we write :

1
2
3
4
5
6
7
8
Xh uh, vh;
Xph uph, vph;
Mh ph, qh;
Rh [Uxh, Uyh], [Vxh, Vyh];
Xh[int] Uh(10);         //array of 10 functions in Xh
Rh[int] [Wxh, Wyh](10); //array of 10 functions in Rh
Wxh[5](0.5,0.5);        //the 6th function at point (0.5, 0.5)
Wxh[5][];               //the array of the degree of freedom of the 6th function

The functions U_{h}, V_{h} have two components so we have

Use of fespace in 3D#

With the 3D finite element spaces

X_{ph} = \left\{ v \in X_{h} |\; v\left(\vecttwo{0}{.}\right) = v\left(\vecttwo{1}{.}\right) , v\left(\vecttwo{.}{0}\right) = v\left(\vecttwo{.}{1}\right) \right\}
M_{h} = \{ v \in H^{1}(]0,1[^2) |\; \forall K \in \mathcal{T}_{h}\quad v_{|K} \in P_{2} \}
R_{h} = \left\{ \mathbf{v} \in H^{1}(]0,1[^2)^{2} |\; \forall K \in \mathcal{T}_{h}\quad \mathbf{v}_{|K}(x,y) = \vecttwo{\alpha_{K}}{\beta_{K}} + \gamma_{K}\vecttwo{x}{y} \right\}

when \mathcal{T}_h is a mesh 10\times 10\times 10 of the unit cubic ]0,1[^2, we write in FreeFem++ :

1
2
3
4
5
6
7
8
9
mesh3 Th = buildlayers(square(10, 10),10, zbound=[0,1]);
//label: 0 up, 1 down, 2 front, 3 left, 4 back, 5 right
fespace Xh(Th, P1); //scalar FE
fespace Xph(Th, P1,
        periodic=[[0, x, y], [1, x, y],
                  [2, x, z], [4, x, z],
                  [3, y, z], [5, y, z]]); //three-periodic FE
fespace Mh(Th, P2); //scalar FE
fespace Rh(Th, RT03d); //vectorial FE

where Xh, Mh, Rh expresses finite element spaces (called FE spaces) X_h,\, M_h,\, R_h, respectively.

To define and use FE-functions u_{h},v_{h} \in X_{h}, p_{h},q_{h} \in M_{h} and U_{h},V_{h} \in R_{h}, we write:

1
2
3
4
5
6
7
8
Xh uh, vh;
Xph uph, vph;
Mh ph, qh;
Rh [Uxh, Uyh, Uyzh], [Vxh, Vyh, Vyzh];
Xh[int] Uh(10);             //array of 10 functions in Xh
Rh[int] [Wxh,Wyh,Wzh](10);  // array of 10 functions in Rh
Wxh[5](0.5,0.5,0.5);        //the 6th function at point (0.5, 0.5, 0.5)
Wxh[5][];                   //the array of the degree of freedom of the 6th function

The functions U_{h}, V_{h} have three components, so we have

Note

One challenge of the periodic boundary condition is that the mesh must have equivalent faces. The ::freefem buildlayers mesh generator splits each quadrilateral face with the diagonal passing through the vertex with maximum number, so to be sure to have the same mesh one both face periodic the 2D numbering in corresponding edges must be compatible (for example the same variation). \codered

By Default, the numbering of square vertex is correct.

To change the mesh numbering you can use the change function like:

1
2
3
4
5
6
7
8
{
    int[int] old2new(0:Th.nv-1); //array set on 0, 1, .., nv-1
    fespace Vh2(Th, P1);
    Vh2 sorder = x+y; //choose an order increasing on 4 square borders with x or y
    sort(sorder[], old2new); //build the inverse permutation
    int[int] new2old = old2new^-1; //inverse the permutation
    Th = change(Th, renumv=new2old);
}

The full example is in Examples.

Lagrangian Finite Elements#

P0-element#

For each triangle (d=2) or tetrahedron (d=3) T_k, the basis function \phi_k in Vh(Th, P0) is given by

\phi_k(\mathbf{x})= \left\{ \begin{array}{cl} 1 & \textrm{ if }(\mathbf{x})\in T_k\\ 0 & \textrm{ if }(\mathbf{x})\not\in T_k \end{array} \right.

If we write

1
2
Vh(Th, P0);
Vh fh = f(x,y);

then for vertices q^{k_i},\, i=1,2,.. d+1 in Fig. 1(a), f_h is built as

fh= \displaystyle f_h(x,y)=\sum_k f(\frac{\sum_i q^{k_i}}{d+1}) \phi_k

See Fig. 3 for the projection of f(x,y)=\sin(\pi x)\cos(\pi y) on Vh(Th, P0) when the mesh Th is a 4\times 4-grid of [-1,1]^2 as in Fig. 2.

P1-element#

Fig. 1: P_1 and P_2 degrees of freedom on triangle T_k
P1P2

For each vertex q^i, the basis function \phi_i in Vh(Th, P1) is given by

\begin{eqnarray*} &&\phi_i(x,y)=a^k_i+b^k_ix+c^k_iy \textrm{ for }(x,y)\in T_k,\\ &&\phi_i(q^i)=1,\quad \phi_i(q^j)=0 \textrm{ if }i\neq j \end{eqnarray*}

The basis function \phi_{k_1}(x,y) with the vertex q^{k_1} in Fig. 1(a) at point p=(x,y) in triangle T_k simply coincide with the barycentric coordinates \lambda^k_1 (area coordinates) :

If we write

1
2
Vh(Th, P1);
Vh fh = g(x.y);

then

fh = \displaystyle f_h(x,y)=\sum_{i=1}^{n_v}f(q^i)\phi_i(x,y)

See Fig. 4 for the projection of f(x,y)=\sin(\pi x)\cos(\pi y) into Vh(Th, P1).

Fig. 2: Test mesh Th for projection Fig. 3: Projection to Vh(Th, P0)
P0P1P2nc peojP0

P2-element#

For each vertex or mid-point q^i. The basis function \phi_i in Vh(Th, P2) is given by :

\begin{eqnarray*} &&\phi_i(x,y)=a^k_i+b^k_ix+c^k_iy+d^k_ix^2+e^k_ixy+f^f_jy^2\textrm{ for }(x,y)\in T_k,\\ &&\phi_i(q^i)=1,\quad \phi_i(q^j)=0\textrm{ if }i\neq j \end{eqnarray*}

The basis function \phi_{k_1}(x,y) with the vertex q^{k_1} in Fig. 1(b) is defined by the barycentric coordinates:

\phi_{k_1}(x,y) = \lambda^k_{1}(x,y)(2\lambda^k_1(x,y)-1)

and for the mid-point q^{k_2}

\phi_{k_2}(x,y) = 4\lambda^k_1(x,y)\lambda^k_4(x,y)

If we write :

1
2
Vh(Th, P2);
Vh fh = f(x.y);

then :

fh = \displaystyle f_h(x,y)=\sum_{i=1}^{M}f(q^i)\phi_i(x,y)\quad (\textrm{summation over all vertex or mid-point})

See Fig. 5 for the projection of f(x,y)=\sin(\pi x)\cos(\pi y) into Vh(Th, P2).

Fig. 4: projection to Vh(Th, P1) Fig. 5: projection to Vh(Th, P2)
projP1 projP2

P1 Nonconforming Element#

Refer to THOMASSET2012 for details; briefly, we now consider non-continuous approximations so we will lose the property

w_h\in V_h\subset H^1(\Omega)

If we write

1
2
Vh(Th, P1nc);
Vh fh = f(x.y);

then

fh = \displaystyle f_h(x,y)=\sum_{i=1}^{n_v}f(m^i)\phi_i(x,y)\quad (\textrm{summation over all midpoint})

Here the basis function \phi_i associated with the mid-point m^i=(q^{k_i}+q^{k_{i+1}})/2 where q^{k_i} is the i-th point in T_k, and we assume that j+1=0 if j=3:

\begin{eqnarray*} &&\phi_i(x,y)=a^k_i+b^k_ix+c^k_iy~\textrm{for }(x,y)\in T_k,\\ &&\phi_i(m^i)=1,\quad \phi_i(m^j)=0\textrm{ if }i\neq j \end{eqnarray*}

Strictly speaking \p \phi_i/\p x,\, \p \phi_i/\p y contain Dirac distribution \rho \delta_{\p T_k}.

The numerical calculations will automatically ignore them. In THOMASSET2012, there is a proof of the estimation

The basis functions \phi_k have the following properties.

  1. For the bilinear form a defined in Fig. 6 satisfy

  2. f\ge 0 \Rightarrow u_h\ge 0

  3. If i\neq j, the basis function \phi_i and \phi_j are L^2-orthogonal: which is false for P_1-element.

See Fig. 6 for the projection of f(x,y)=\sin(\pi x)\cos(\pi y) into Vh(Th, P1nc).

Fig. 6: Projection to Vh(Th, P1nc) Fig. 7: Projection to Vh(Th, P1b)
p1nc projP1b

Other FE-space#

For each triangle T_k\in \mathcal{T}_h, let \lambda_{k_1}(x,y),\, \lambda_{k_2}(x,y),\, \lambda_{k_3}(x,y) be the area cordinate of the triangle (see Fig. 1), and put

\begin{equation} \beta_k(x,y)=27\lambda_{k_1}(x,y)\lambda_{k_2}(x,y)\lambda_{k_3}(x,y) \end{equation}

called bubble function on T_k. The bubble function has the feature: 1. \beta_k(x,y)=0\quad \textrm{if }(x,y)\in \p T_k.

  1. \beta_k(q^{k_b})=1 where q^{k_b} is the barycenter \frac{q^{k_1}+q^{k_2}+q^{k_3}}{3}.

If we write :

1
2
Vh(Th, P1b);
Vh fh = f(x.y);

then

fh = \displaystyle f_h(x,y)=\sum_{i=1}^{n_v}f(q^i)\phi_i(x,y)+\sum_{k=1}^{n_t}f(q^{k_b})\beta_k(x,y)

See Fig. 7 for the projection of f(x,y)=\sin(\pi x)\cos(\pi y) into Vh(Th, P1b).

Vector Valued FE-function#

Functions from \R^{2} to \R^{N} with N=1 are called scalar functions and called vector valued when N>1. When N=2

1
fespace Vh(Th, [P0, P1]) ;

makes the space

V_h=\{\mathbf{w}=(w_1,w_2)|\; w_1\in V_h(\mathcal{T}_h,P_0),\, w_2\in V_h(\mathcal{T}_h,P_1)\}

Raviart-Thomas Element#

In the Raviart-Thomas finite element RT0_{h}, the degrees of freedom are the fluxes across edges e of the mesh, where the flux of the function \mathbf{f} : \R^2 \longrightarrow \R^2 is \int_{e} \mathbf{f}.n_{e}, n_{e} is the unit normal of edge e.

This implies an orientation of all the edges of the mesh, for example we can use the global numbering of the edge vertices and we just go from small to large numbers.

To compute the flux, we use a quadrature with one Gauss point, the mid-point of the edge.

Consider a triangle T_k with three vertices (\mathbf{a},\mathbf{b},\mathbf{c}). Lets denote the vertices numbers by i_{a},i_{b},i_{c}, and define the three edge vectors \mathbf{e}^{1},\mathbf{e}^{2},\mathbf{e}^{3} by sgn(i_{b}-i_{c})(\mathbf{b}-\mathbf{c}), sgn(i_{c}-i_{a})(\mathbf{c}-\mathbf{a}), sgn(i_{a}-i_{b})(\mathbf{a}-\mathbf{b}).

We get three basis functions :

\begin{equation} \boldsymbol{\phi}^{k}_{1}= \frac{sgn(i_{b}-i_{c})}{2|T_k|}(\mathbf{x}-\mathbf{a}),\quad \boldsymbol{\phi}^{k}_{2}= \frac{sgn(i_{c}-i_{a})}{2|T_k|}(\mathbf{x}-\mathbf{b}),\quad \boldsymbol{\phi}^{k}_{3}= \frac{sgn(i_{a}-i_{b})}{2|T_k|}(\mathbf{x}-\mathbf{c}), \end{equation}

where |T_k| is the area of the triangle T_k. If we write

1
2
Vh(Th, RT0);
Vh [f1h, f2h] = [f1(x, y), f2(x, y)];

then

fh = \displaystyle \mathbf{f}_h(x,y)=\sum_{k=1}^{n_t}\sum_{l=1}^6 n_{i_lj_l}|\mathbf{e^{i_l}}|f_{j_l}(m^{i_l})\phi_{i_lj_l}

where n_{i_lj_l} is the j_l-th component of the normal vector \mathbf{n}_{i_l},

\{m_1,m_2,m_3\} = \left\{\frac{\mathbf{b}+\mathbf{c}}{2}, \frac{\mathbf{a}+\mathbf{c}}{2}, \frac{\mathbf{b}+\mathbf{a}}{2} \right\}

and i_l=\{1,1,2,2,3,3\},\, j_l=\{1,2,1,2,1,2\} with the order of l.

Fig. 8: Normal vectors of each edge
RT0

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
// Mesh
mesh Th = square(2, 2);

// Fespace
fespace Xh(Th, P1);
Xh uh = x^2 + y^2, vh;

fespace Vh(Th, RT0);
Vh [Uxh, Uyh] = [sin(x), cos(y)]; //vectorial FE function

// Change the mesh
Th = square(5,5);

//Xh is unchanged
//Uxh = x; //error: impossible to set only 1 component
           //of a vector FE function
vh = Uxh;//ok
//and now vh use the 5x5 mesh
//but the fespace of vh is always the 2x2 mesh

// Plot
plot(uh);
uh = uh; //do a interpolation of uh (old) of 5x5 mesh
         //to get the new uh on 10x10 mesh
plot(uh);

vh([x-1/2, y]) = x^2 + y^2; //interpolate vh = ((x-1/2)^2 + y^2)
Fig. 9: vh Iso on mesh 2\times 2 Fig. 10: vh Iso on mesh 5\times 5
onoldmesh onnewmesh

To get the value at a point x=1,y=2 of the FE function uh, or [Uxh, Uyh], one writes :

1
2
3
4
5
6
7
8
real value;
value = uh(2,4); //get value = uh(2, 4)
value = Uxh(2, 4); //get value = Uxh(2, 4)
//OR
x = 1; y = 2;
value = uh; //get value = uh(1, 2)
value = Uxh; //get value = Uxh(1, 2)
value = Uyh; //get value = Uyh(1, 2)

To get the value of the array associated to the FE function uh, one writes

1
2
3
4
real value = uh[][0]; //get the value of degree of freedom 0
real maxdf = uh[].max; //maximum value of degree of freedom
int size = uh.n; //the number of degree of freedom
real[int] array(uh.n) = uh[]; //copy the array of the function uh

Warning

For a non-scalar finite element function [Uxh, Uyh] the two arrays Uxh[] and Uyh[] are the same array, because the degree of freedom can touch more than one component.

A Fast Finite Element Interpolator#

In practice, one may discretize the variational equations by the Finite Element method. Then there will be one mesh for \Omega_1 and another one for \Omega_2. The computation of integrals of products of functions defined on different meshes is difficult.

Quadrature formula and interpolations from one mesh to another at quadrature points are needed. We present below the interpolation operator which we have used and which is new, to the best of our knowledge.

Let {\cal T}_{h}^0=\cup_k T^0_k,{\cal T}_{h}^1=\cup_k T^1_k be two triangulations of a domain \Omega. Let

V({\hbox{${\cal T}$}_{h}^i}) =\{ C^0(\Omega_h^i)~:~f|_{T^i_k}\in P_0\},~~~i=0,1

be the spaces of continuous piecewise affine functions on each triangulation.

Let f\in V({\cal T}_{h}^0). The problem is to find g\in V({\cal T}_{h}^1) such that

g(q) = f(q) \quad \forall q\hbox{~vertex of ~} {\cal T}_{h}^1

Although this is a seemingly simple problem, it is difficult to find an efficient algorithm in practice.

We propose an algorithm which is of complexity N^1\log N^0, where N^i is the number of vertices of \cal T_{h}^i, and which is very fast for most practical 2D applications.

Algorithm

The method has 5 steps. First a quadtree is built containing all the vertices of the mesh {\cal T}_{h}^0 such that in each terminal cell there are at least one, and at most 4, vertices of {\cal T}_{h}^0. For each q^1, vertex of {\cal T}_{h}^1 do:

  1. Find the terminal cell of the quadtree containing q^1.
  2. Find the the nearest vertex q^0_j to q^1 in that cell.
  3. Choose one triangle T_k^0\in{\cal T}_{h}^0 which has q^0_j for vertex.
  4. Compute the barycentric coordinates \{\lambda_j\}_{j=1,2,3} of q^1 in T^0_k.

    • if all barycentric coordinates are positive, go to Step 5
    • otherwise, if one barycentric coordinate \lambda_i is negative, replace T^0_k by the adjacent triangle opposite q^0_i and go to Step 4.
    • otherwise, if two barycentric coordinates are negative, take one of the two randomly and replace T^0_k by the adjacent triangle as above.
  5. Calculate g(q^1) on T^0_k by linear interpolation of f:

Fig. 11: To interpolate a function at q^0, the knowledge of the triangle which contains q^0 is needed. The algorithm may start at q^1\in T_k^0 and stall on the boundary (thick line) because the line q^0q^1 is not inside \Omega. But if the holes are triangulated too (doted line) then the problem does not arise.
fastInterpolate

Two problems need to be solved:

  • What if q^1 is not in \Omega^0_h ? Then Step 5 will stop with a boundary triangle. So we add a step which tests the distance of q^1 with the two adjacent boundary edges and selects the nearest, and so on till the distance grows.

  • What if \Omega^0_h is not convex and the marching process of Step 4 locks on a boundary? By construction Delaunay-Voronoï's mesh generators always triangulate the convex hull of the vertices of the domain. Therefore, we make sure that this information is not lost when {\cal T}_{h}^0,{\cal T}_{h}^1 are constructed and we keep the triangles which are outside the domain on a special list. That way, in step 5 we can use that list to step over holes if needed.

Note

Sometimes, in rare cases, the interpolation process misses some points, we can change the search algorithm through a global variable searchMethod

1
2
3
4
searchMethod = 0; // default value for fast search algorithm
searchMethod = 1; // safe search algorithm, uses brute force in case of missing point
// (warning: can be very expensive in cases where a lot of points are outside of the domain)
searchMethod = 2; // always uses brute force. It is very computationally expensive.

Note

Step 3 requires an array of pointers such that each vertex points to one triangle of the triangulation.

Note

The operator = is the interpolation operator of FreeFem++, the continuous finite functions are extended by continuity to the outside of the domain.

Try the following example :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
// Mesh
mesh Ths = square(10, 10);
mesh Thg = square(30, 30, [x*3-1, y*3-1]);
plot(Ths, Thg, wait=true);

// Fespace
fespace Ch(Ths, P2);
Ch us = (x-0.5)*(y-0.5);

fespace Dh(Ths, P2dc);
Dh vs = (x-0.5)*(y-0.5);

fespace Fh(Thg, P2dc);
Fh ug=us, vg=vs;

// Plot
plot(us, ug, wait=true);
plot(vs, vg, wait=true);

Fig. 12: Extension of a continuous FE-function Fig. 13: Extension of discontinuous FE-function
us-ug vs-vg

Keywords: Problem and Solve#

For FreeFem++, a problem must be given in variational form, so we need a bilinear form a(u,v), a linear form \ell(f,v), and possibly a boundary condition form must be added.

1
2
3
4
problem P (u, v)
    = a(u,v) - l(f,v)
    + (boundary condition)
    ;

Note

When you want to formulate the problem and solve it in the same time, you can use the keyword solve.

Weak Form and Boundary Condition#

To present the principles of Variational Formulations, also called weak form, for the Partial Differential Equations, let's take a model problem: a Poisson equation with Dirichlet and Robin Boundary condition.

The problem: Find u a real function defined on a domain \Omega of \R^d (d=2,3) such that

\begin{equation} \begin{array}{rcll} -\nabla\cdot(\kappa \nabla u) &=& f & \mbox{ in }\Omega\\ a u + \kappa \frac{\p u}{\p n} &=& b & \mbox{ on }\Gamma_r\\ u &=& g & \mbox{ on }\Gamma_d \end{array} \end{equation}

where

  • if d=2 then \nabla.(\kappa \nabla u) = \p_x(\kappa \p_x u ) + \p_y(\kappa \p_y u ) with \p_x u = \frac{\p u}{\p x} and \p_y u = \frac{\p u}{\p y}

  • if d=3 then \nabla.(\kappa \nabla u) = \p_x(\kappa \p_x u) + \p_y(\kappa \p_y u) + \p_z(\kappa \p_z u) with \p_x u = \frac{\p u}{\p x}, \p_y u = \frac{\p u}{\p y} and , \p_z u = \frac{\p u}{\p z}

  • The border \Gamma=\p \Omega is split in \Gamma_d and \Gamma_n such that \Gamma_d \cap \Gamma_n = \emptyset and \Gamma_d \cup \Gamma_n = \p \Omega,

  • \kappa is a given positive function, such that \exists \kappa_0 \in \R ,\quad 0 < \kappa_0 \leq \kappa.

  • a a given non negative function,

  • b a given function.

Note

This is the well known Neumann boundary condition if a=0, and if \Gamma_d is empty. In this case the function appears in the problem just by its derivatives, so it is defined only up to a constant (if u is a solution then u+c is also a solution).

Let {v}, a regular test function, null on \Gamma_d, by integration by parts we get :

where if d=2 the \nabla{ v} . \nabla u = (\frac{\p u}{\p x}\frac{\p { v}}{\p x}+\frac{\p u}{\p y}\frac{\p { v}}{\p y}),

where if d=3 the \nabla{ v} . \nabla u = (\frac{\p u}{\p x}\frac{\p { v}}{\p x}+\frac{\p u}{\p y}\frac{\p { v}}{\p y} + \frac{\p u}{\p z}\frac{\p { v}}{\p z}),

and where \mathbf{n} is the unitary outer-pointing normal of the \Gamma.

Now we note that \kappa \frac{ \p u}{\p n} = - a u + b on \Gamma_r and v=0 on \Gamma_d and \Gamma = \Gamma_d \cup \Gamma_n thus

- \int_{\Gamma} {v} \kappa \frac{ \p u}{\p n} = \int_{\Gamma_r} a u v - \int_{\Gamma_r} b v

The problem becomes:

Find u \in V_g = \{w \in H^1(\Omega) / w = g \mbox{ on } \Gamma_d \} such that

where V_0 = \{v \in H^1(\Omega) / v = 0 \mbox{ on } \Gamma_d \}

Except in the case of Neumann conditions everywhere, the problem \eqref{eqn::v-poisson} is well posed when \kappa\geq \kappa_0>0.

Note

If we have only the Neumann boundary condition, linear algebra tells us that the right hand side must be orthogonal to the kernel of the operator for the solution to exist. One way of writing the compatibility condition is:

\int_{\Omega} f \,d\omega + \int_{\Gamma} b \,d\gamma=0

and a way to fix the constant is to solve for u \in H^1(\Omega) such that: where \varepsilon is a small parameter ($ \sim \kappa\; 10^{-10} |\Omega|^{\frac2d} $).

Remark that if the solution is of order \frac{1}{\varepsilon} then the compatibility condition is unsatisfied, otherwise we get the solution such that $\int_\Omega u = 0 $, you can also add a Lagrange multiplier to solve the real mathematical problem like in the Lagrange multipliers example.

In FreeFem++, the bidimensional problem \eqref{eqn::v-poisson} becomes :

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
problem Pw (u, v)
    = int2d(Th)( //int_{Omega} kappa nabla v . nabla u
          kappa*(dx(u)*dx(v) + dy(u)*dy(v))
    )
    + int1d(Th, gn)( //int_{Gamma_r} a u v
          a * u*v
    )
    - int2d(Th)( //int_{Omega} f v
          f*v
    )
    - int1d(Th, gn)( //int_{Gamma_r} b v
          b * v
    )
    + on(gd, u=g) //u = g on Gamma_d
    ;

where Th is a mesh of the bi-dimensional domain \Omega, and gd and gn are respectively the boundary labels of boundary \Gamma_d and \Gamma_n.

And the three dimensional problem \eqref{eqn::v-poisson} becomes

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
macro Grad(u) [dx(u), dy(u), dz(u) ]//
problem Pw (u, v)
    = int3d(Th)( //int_{Omega} kappa nabla v . nabla u
          kappa*(Grad(u)'*Grad(v))
    )
    + int2d(Th, gn)( //int_{Gamma_r} a u v
          a * u*v
    )
    - int3d(Th)( //int_{Omega} f v
          f*v
    )
    - int2d(Th, gn)( //int_{Gamma_r} b v
          b * v
    )
    + on(gd, u=g) //u = g on Gamma_d
    ;

where Th is a mesh of the three dimensional domain \Omega, and gd and gn are respectively the boundary labels of boundary \Gamma_d and \Gamma_n.

Parameters affecting solve and problem#

The parameters are FE functions real or complex, the number n of parameters is even (n=2*k), the k first function parameters are unknown, and the k last are test functions.

Note

If the functions are a part of vectorial FE then you must give all the functions of the vectorial FE in the same order (see Poisson problem with mixed finite element for example).

Note

Don't mix complex and real parameters FE function.

Warning

Bug : The mixing of multiple fespace with different periodic boundary conditions are not implemented. So all the finite element spaces used for tests or unknown functions in a problem, must have the same type of periodic boundary conditions or no periodic boundary conditions. No clean message is given and the result is unpredictable.

The parameters are:

  • solver= LU, CG, Crout,Cholesky,GMRES,sparsesolver, UMFPACK ...

    The default solver is sparsesolver (it is equal to UMFPACK if no other sparse solver is defined) or is set to LU if no direct sparse solver is available.

    The storage mode of the matrix of the underlying linear system depends on the type of solver chosen; for LU the matrix is sky-line non symmetric, for Crout the matrix is sky-line symmetric, for Cholesky the matrix is sky-line symmetric positive definite, for CG the matrix is sparse symmetric positive, and for GMRES, sparsesolver or UMFPACK the matrix is just sparse.

  • eps= a real expression. \varepsilon sets the stopping test for the iterative methods like CG. Note that if \varepsilon is negative then the stopping test is: if it is positive, then the stopping test is

  • init= boolean expression, if it is false or 0 the matrix is reconstructed. Note that if the mesh changes the matrix is reconstructed too.

  • precon= name of a function (for example P) to set the preconditioner. The prototype for the function P must be

    1
    func real[int] P(real[int] & xx);
    
  • tgv= Huge value (10^{30}) used to implement Dirichlet boundary conditions.

  • tolpivot= sets the tolerance of the pivot in UMFPACK (10^{-1}) and, LU, Crout, Cholesky factorisation (10^{-20}).

  • tolpivotsym= sets the tolerance of the pivot sym in UMFPACK

  • strategy= sets the integer UMFPACK strategy (0 by default).

Problem definition#

Below v is the unknown function and w is the test function. After the "=" sign, one may find sums of:

  • Identifier(s); this is the name given earlier to the variational form(s) (type varf ) for possible reuse.

    Remark, that the name in the varf of the unknown test function is forgotten, we use the order in the argument list to recall names as in a C++ function,

  • The terms of the bilinear form itself: if K is a given function,

  • Bilinear part for 3D meshes Th

    • int3d(Th)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{T } K\,v\,w

    • int3d(Th, 1)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th},T\subset \Omega_{1}}\int_{T} K\,v\,w

    • int3d(Th, levelset=phi)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{T,\phi<0} K\,v\,w

    • int3d(Th, l, levelset=phi)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th},T\subset \Omega_{l}}\int_{T,\phi<0} K\,v\,w

    • int2d(Th, 2, 5)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{(\p T\cup\Gamma) \cap ( \Gamma_2 \cup \Gamma_{5})} K\,v\,w

    • int2d(Th, 1)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th},T\subset \Omega_{1}}\int_{T} K\,v\,w

    • int2d(Th, 2, 5)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{(\p T\cup\Gamma) \cap (\Gamma_2 \cup \Gamma_{5})} K\,v\,w

    • int2d(Th, levelset=phi)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{T,\phi=0} K\,v\,w

    • int2d(Th, l, levelset=phi)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th},T\subset \Omega_{l}}\int_{T,\phi=0} K\,v\,w

    • intallfaces(Th)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{\p T } K\,v\,w

    • intallfaces(Th, 1)(K*v*w) = \displaystyle\sum_{{T\in\mathtt{Th},T\subset \Omega_{1}}}\int_{\p T } K\,v\,w

    • They contribute to the sparse matrix of type matrix which, whether declared explicitly or not, is constructed by FreeFem++.

  • Bilinear part for 2D meshes Th

    • int2d(Th)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{T } K\,v\,w

    • int2d(Th, 1)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th},T\subset \Omega_{1}}\int_{T} K\,v\,w

    • int2d(Th, levelset=phi)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{T,\phi<0} K\,v\,w

    • int2d(Th, l, levelset=phi)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th},T\subset \Omega_{l}}\int_{T,\phi<0} K\,v\,w

    • int1d(Th, 2, 5)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{(\p T\cup\Gamma) \cap ( \Gamma_2 \cup \Gamma_{5})} K\,v\,w

    • int1d(Th, 1)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th},T\subset \Omega_{1}}\int_{T} K\,v\,w

    • int1d(Th, 2, 5)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{(\p T\cup\Gamma) \cap ( \Gamma_2 \cup \Gamma_{5})} K\,v\,w

    • int1d(Th, levelset=phi)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{T,\phi=0} K\,v\,w

    • int1d(Th, l, levelset=phi)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th},T\subset \Omega_{l}}\int_{T,\phi=0} K\,v\,w

    • intalledges(Th)(K*v*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{\p T } K\,v\,w

    • intalledges(Th, 1)(K*v*w) = \displaystyle\sum_{{T\in\mathtt{Th},T\subset \Omega_{1}}}\int_{\p T } K\,v\,w

    • They contribute to the sparse matrix of type matrix which, whether declared explicitly or not, is constructed by FreeFem++.

  • The right hand-side of the Partial Differential Equation in 3D, the terms of the linear form: for given functions K,\, f:

    • int3d(Th)(K*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{T} K\,w

    • int3d(Th, l)(K*w) = \displaystyle\sum_{T\in\mathtt{Th},T\in\Omega_l}\int_{T} K\,w

    • int3d(Th, levelset=phi)(K*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{T,\phi<0} K\,w

    • int3d(Th, l, levelset=phi)(K*w) = \displaystyle\sum_{T\in\mathtt{Th},T\subset\Omega_{l}}\int_{T,\phi<0} K\,w

    • int2d(Th, 2, 5)(K*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{(\p T\cup\Gamma) \cap ( \Gamma_2 \cup \Gamma_{5}) } K \,w

    • int2d(Th, levelset=phi)(K*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{T,\phi=0} K\,w

    • int2d(Th, l, levelset=phi)(K*w) = \displaystyle\sum_{T\in\mathtt{Th},T\subset \Omega_{l}}\int_{T,\phi=0} K\,w

    • intallfaces(Th)(f*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{\p T } f\,w

    • A vector of type real[int]

  • The right hand-side of the Partial Differential Equation in 2D, the terms of the linear form: for given functions K,\, f:

    • int2d(Th)(K*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{T} K\,w

    • int2d(Th, l)(K*w) = \displaystyle\sum_{T\in\mathtt{Th},T\in\Omega_l}\int_{T} K\,w

    • int2d(Th, levelset=phi)(K*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{T,\phi<0} K\,w

    • int2d(Th, l, levelset=phi)(K*w) = \displaystyle\sum_{T\in\mathtt{Th},T\subset\Omega_{l}}\int_{T,\phi<0} K\,w

    • int1d(Th, 2, 5)(K*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{(\p T\cup\Gamma) \cap ( \Gamma_2 \cup \Gamma_{5}) } K \,w

    • int1d(Th, levelset=phi)(K*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{T,\phi=0} K\,w

    • int1d(Th, l, levelset=phi)(K*w) = \displaystyle\sum_{T\in\mathtt{Th},T\subset\Omega_{l}}\int_{T,\phi=0} K\,w

    • intalledges(Th)(f*w) = \displaystyle\sum_{T\in\mathtt{Th}}\int_{\p T } f\,w

    • a vector of type real[int]

  • The boundary condition terms:

    • An "on" scalar form (for Dirichlet) : on(1, u=g)

      Used for all degrees of freedom i of the boundary referred by "1", the diagonal term of the matrix a_{ii}= tgv with the terrible giant value tgv (=10^{30} by default), and the right hand side b[i] = "(\Pi_h g)[i]" \times tgv, where the "(\Pi_h g)g[i]" is the boundary node value given by the interpolation of g. \codered

      Note

      if \mathrm{tgv} < 0 then we put to 0 all term of the line i in the matrix, except diagonal term a_{ii}=1, and b[i] = "(\Pi_h g)[i]".

    • An "on" vectorial form (for Dirichlet) : on(1, u1=g1, u2=g2)

      If you have vectorial finite element like RT0, the 2 components are coupled, and so you have : b[i] = "(\Pi_h (g1,g2))[i]" \times tgv, where \Pi_h is the vectorial finite element interpolant.

    • A linear form on \Gamma (for Neumann in 2d) -int1d(Th)(f*w) or -int1d(Th, 3)(f*w)

    • A bilinear form on \Gamma or \Gamma_{2} (for Robin in 2d) int1d(Th)(K*v*w) or int1d(Th,2)(K*v*w)

    • A linear form on \Gamma (for Neumann in 3d) -int2d(Th)(f*w) or -int2d(Th, 3)(f*w)

    • A bilinear form on \Gamma or \Gamma_{2} (for Robin in 3d) int2d(Th)(K*v*w) or int2d(Th,2)(K*v*w)

Note

  • If needed, the different kind of terms in the sum can appear more than once.
  • The integral mesh and the mesh associated to test functions or unknown functions can be different in the case of linear form.
  • N.x, N.y and N.z are the normal's components.

Warning

It is not possible to write in the same integral the linear part and the bilinear part such as in int1d(Th)(K*v*w - f*w).

Numerical Integration#

Let D be a N-dimensional bounded domain.

For an arbitrary polynomial f of degree r, if we can find particular (quadrature) points \mathbf{\xi}_j,\, j=1,\cdots,J in D and (quadrature) constants \omega_j such that

\begin{equation} \int_{D}f(\mathbf{x}) = \sum_{\ell =1}^L c_\ell f(\mathbf{\xi}_\ell) \end{equation}

then we have an error estimate (see CROUZEIX1984), and then there exists a constant C>0 such that

\begin{equation} \left|\int_{D}f(\mathbf{x}) - \sum_{\ell =1}^L \omega_\ell f(\mathbf{\xi}_\ell )\right| \le C|D|h^{r+1} \end{equation}

for any function r + 1 times continuously differentiable f in D, where h is the diameter of D and |D| its measure (a point in the segment [q^iq^j] is given as

\{(x,y)|\; x=(1-t)q^i_x+tq^j_x,\, y=(1-t)q^i_y+tq^j_y,\, 0\le t\le 1\} ).

For a domain \Omega_h=\sum_{k=1}^{n_t}T_k,\, \mathcal{T}_h=\{T_k\}, we can calculate the integral over \Gamma_h=\p\Omega_h by

\int_{\Gamma_h}f(\mathbf{x})ds =int1d(Th)(f)

=int1d(Th, qfe=*)(f)

=int1d(Th, qforder=*)(f)

where * stands for the name of the quadrature formula or the precision (order) of the Gauss formula.

Quadature formula on an edge
L (qfe=) qforder= Point in [q^iq^j](=t) \omega_\ell Exact on P_k, k=
1 qf1pE 2 1/2 \|q^iq^j\| 1
2 qf2pE 3 (1\pm\sqrt{1/3})/2 \|q^iq^j\|/2 3
3 qf3pE 6 (1\pm\sqrt{3/5})/2
1/2
(5/18)\|q^iq^j\|
(8/18)\|q^iq^j\|
5
4 qf4pE 8 (1\pm\frac{\sqrt{525+70\sqrt{30}}}{35})/2
(1\pm\frac{\sqrt{525-70\sqrt{30}}}{35})/2
\frac{18-\sqrt{30}}{72}\|q^iq^j\|
\frac{18+\sqrt{30}}{72}\|q^iq^j\|
7
5 qf5pE 10 (1\pm\frac{\sqrt{245+14\sqrt{70}}}{21})/2
1/2
(1\pm\frac{\sqrt{245-14\sqrt{70}}}{21})/2
\frac{322-13\sqrt{70}}{1800}\|q^iq^j\|
\frac{64}{225}\|q^iq^j\|
\frac{322+13\sqrt{70}}{1800}\|q^iq^j\|
9
2 qf1pElump 2 0
+1
\|q^iq^j\|/2
\|q^iq^j\|/2
1

where |q^iq^j| is the length of segment \overline{q^iq^j}.

For a part \Gamma_1 of \Gamma_h with the label "1", we can calculate the integral over \Gamma_1 by

\int_{\Gamma_1}f(x,y)ds =int1d(Th, 1)(f)

=int1d(Th, 1, qfe=qf2pE)(f)

The integrals over \Gamma_1,\, \Gamma_3 are given by

\int_{\Gamma_1\cup \Gamma_3}f(x,y)ds =int1d(Th, 1, 3)(f)

For each triangle T_k=[q^{k_1}q^{k_2}q^{k_3}], the point P(x,y) in T_k is expressed by the area coordinate as P(\xi,\eta):

\begin{eqnarray*} &&|T_k|=\frac12 \left| \begin{array}{ccc} 1&q^{k_1}_x&q^{k_1}_y\\ 1&q^{k_2}_x&q^{k_2}_y\\ 1&q^{k_3}_x&q^{k_3}_y \end{array} \right|\quad D_1=\left| \begin{array}{ccc} 1&x&y\\ 1&q^{k_2}_x&q^{k_2}_y\\ 1&q^{k_3}_x&q^{k_3}_y \end{array} \right| \quad D_2=\left| \begin{array}{ccc} 1&q^{k_1}_x&q^{k_1}_y\\ 1&x&y\\ 1&q^{k_3}_x&q^{k_3}_y \end{array} \right| \quad D_3=\left| \begin{array}{ccc} 1&q^{k_1}_x&q^{k_1}_y\\ 1&q^{k_2}_x&q^{k_2}_y\\ 1&x&y \end{array} \right|\\ &&\xi=\frac12 D_1/|T_k|\qquad \eta=\frac12 D_2/|T_k|\qquad \textrm{then } 1-\xi-\eta=\frac12 D_3/|T_k| \end{eqnarray*}

For a two dimensional domain or a border of three dimensional domain \Omega_h=\sum_{k=1}^{n_t}T_k,\, \mathcal{T}_h=\{T_k\}, we can calculate the integral over \Omega_h by

\int_{\Omega_h}f(x,y) =int2d(Th)(f)

=int2d(Th, qft=*)(f)

=int2d(Th, qforder=*)(f)

where * stands for the name of quadrature formula or the order of the Gauss formula.

Quadature formula on a triangle
L qft= qforder= Point in T_k \omega_\ell Exact on P_k, k=
1 qf1pT 2 \left(\frac{1}{3},\frac{1}{3}\right) \|T_k\| 1
3 qf2pT 3 \left(\frac{1}{2},\frac{1}{2}\right)
\left(\frac{1}{2},0\right)
\left(0,\frac{1}{2}\right)
\|T_k\|/3
\|T_k\|/3
\|T_k\|/3
2
7 qf5pT 6 \left(\frac{1}{3},\frac{1}{3}\right)
\left(\frac{6-\sqrt{15}}{21},\frac{6-\sqrt{15}}{21}\right)
\left(\frac{6-\sqrt{15}}{21},\frac{9+2\sqrt{15}}{21}\right)
\left(\frac{9+2\sqrt{15}}{21},\frac{6-\sqrt{15}}{21}\right)
\left(\frac{6+\sqrt{15}}{21},\frac{6+\sqrt{15}}{21}\right)
\left(\frac{6+\sqrt{15}}{21},\frac{9-2\sqrt{15}}{21}\right)
\left(\frac{9-2\sqrt{15}}{21},\frac{6+\sqrt{15}}{21}\right)
0.225\|T_k\|
\frac{(155-\sqrt{15})\|T_k\|}{1200}
\frac{(155-\sqrt{15})\|T_k\|}{1200}
\frac{(155-\sqrt{15})\|T_k\|}{1200}
\frac{(155+\sqrt{15})\|T_k\|}{1200}
\frac{(155+\sqrt{15})\|T_k\|}{1200}
\frac{(155+\sqrt{15})\|T_k\|}{1200}
5
3 qf1pTlump \left(0,0\right)
\left(1,0\right)
\left(0,1\right)
\|T_k\|/3
\|T_k\|/3
\|T_k\|/3
1
9 qf2pT4P1 \left(\frac{1}{4},\frac{3}{4}\right)
\left(\frac{3}{4},\frac{1}{4}\right)
\left(0,\frac{1}{4}\right)
\left(0,\frac{3}{4}\right)
\left(\frac{1}{4},0\right)
\left(\frac{3}{4},0\right)
\left(\frac{1}{4},\frac{1}{4}\right)
\left(\frac{1}{4},\frac{1}{2}\right)
\left(\frac{1}{2},\frac{1}{4}\right)
\|T_k\|/12
\|T_k\|/12
\|T_k\|/12
\|T_k\|/12
\|T_k\|/12
\|T_k\|/12
\|T_k\|/6
\|T_k\|/6
\|T_k\|/6
1
15 qf7pT 8 See TAYLOR2005 for detail 7
21 qf9pT 10 See TAYLOR2005 for detail 9

For a three dimensional domain \Omega_h=\sum_{k=1}^{n_t}T_k,\, \mathcal{T}_h=\{T_k\}, we can calculate the integral over \Omega_h by

\int_{\Omega_h}f(x,y) =int3d(Th)(f)

=int3d(Th,qfV=*)(f)

=int3D(Th,qforder=*)(f)

where * stands for the name of quadrature formula or the order of the Gauss formula.

Quadature formula on a tetrahedron
L qfV= qforder= Point in T_k\in \R^3 \omega_\ell Exact on P_k, k=
1 qfV1 2 \left(\frac{1}{4},\frac{1}{4},\frac{1}{4}\right) \|T_k\| 1
4 qfV2 3 G4(0.58\ldots,0.13\ldots,0.13\ldots) \|T_k\|/4 2
14 qfV5 6 G4(0.72\ldots,0.092\ldots,0.092\ldots)
G4(0.067\ldots,0.31\ldots,0.31\ldots)
G6(0.45\ldots,0.045\ldots,0.45\ldots)
0.073\ldots\|T_k\|
0.11\ldots\|T_k\|
0.042\ldots\|T_k\|
5
4 qfV1lump G4(1,0,0) \|T_k\|/4 1

Where G4(a,b,b) such that a+3b=1 is the set of the four point in barycentric coordinate and where G6(a,b,b) such that 2a+2b=1 is the set of the six points in barycentric coordinate

Note

These tetrahedral quadrature formulae come from http://www.cs.kuleuven.be/~nines/research/ecf/mtables.html

Note

By default, we use the formula which is exact for polynomials of degree 5 on triangles or edges (in bold in three tables).

This possible to add an own quadrature formulae with using plugin qf11to25 on segment, triangle or Tetrahedron.

The quadrature formulae in D dimension is a bidimentional array of size N_q\times (D+1) such that the D+1 value of on row i=0,...,N_p-1 are w^i,\hat{x}^i_1,...,\hat{x}^i_D where w^i is the weight of the quadrature point, and 1-\sum_{k=1}^D \hat{x}^i_k ,\hat{x}^i_1,...,\hat{x}^i_D is the barycentric coordinate the quadrature point.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
load "qf11to25"

// Quadrature on segment
real[int, int] qq1 = [
    [0.5, 0],
    [0.5, 1]
];

QF1 qf1(1, qq1); //def of quadrature formulae qf1 on segment
//remark:
//1 is the order of the quadrature exact for polynome of degree < 1

//Quadrature on triangle
real[int, int] qq2 = [
    [1./3., 0, 0],
    [1./3., 1, 0],
    [1./3., 0, 1]
];

QF2 qf2(1, qq2); //def of quadrature formulae qf2 on triangle
//remark:
//1 is the order of the quadrature exact for polynome of degree < 1
//so must have sum w^i = 1

// Quadrature on tetrahedron
real[int, int] qq3 = [
    [1./4., 0, 0, 0],
    [1./4., 1, 0, 0],
    [1./4., 0, 1, 0],
    [1./4., 0, 0, 1]
];

QF3 qf3(1, qq3); //def of quadrature formulae qf3 on get
//remark:
//1 is the order of the quadrature exact for polynome of degree < 1)

// Verification in 1d and 2d
mesh Th = square(10, 10);

real I1 = int1d(Th, qfe=qf1)(x^2);
real I1l = int1d(Th, qfe=qf1pElump)(x^2);

real I2 = int2d(Th, qft=qf2)(x^2);
real I2l = int2d(Th, qft=qf1pTlump)(x^2);

cout << I1 << " == " << I1l << endl;
cout << I2 << " == " << I2l << endl;
assert( abs(I1-I1l) < 1e-10 );
assert( abs(I2-I2l) < 1e-10 );

The output is

1
2
1.67 == 1.67
0.335 == 0.335

Variational Form, Sparse Matrix, PDE Data Vector#

In FreeFem++ it is possible to define variational forms, and use them to build matrices and vectors, and store them to speed-up the script (4 times faster here).

For example let us solve the Thermal Conduction problem. The variational formulation is in L^2(0,T;H^1(\Omega)); we shall seek u^n satisfying

\forall w \in V_{0}; \qquad \int_\Omega \frac{u^n-u^{n-1}}{\delta t} w + \kappa\n u^n\n w) +\int_\Gamma\alpha(u^n-u_{ue})w=0

where V_0 = \{w\in H^1(\Omega)/ w_{|\Gamma_{24}}=0\}.

So to code the method with the matrices A=(A_{ij}), M=(M_{ij}), and the vectors u^n, b^n, b',b", b_{cl} (notation if w is a vector then w_i is a component of the vector).

\begin{equation} u^n = A^{-1} b^n, \quad \quad b' = b_0 + M u^{n-1}, \quad b"= \frac{1}{\varepsilon} \; b_{cl}, \quad b^n_i = \left\{ \begin{array}{cl} b"_i & \mbox{if }\ i \in \Gamma_{24} \\ b'_i & \mbox{else if } \not\in \Gamma_{24} \end{array}\right. \end{equation}

Where with \frac{1}{\varepsilon} = \mathtt{tgv} = 10^{30} :

\begin{eqnarray} A_{ij} &=& \left\{\begin{array}{cl} \frac{1}{\varepsilon} & \mbox{if } i \in \Gamma_{24}, \mbox{and} j=i \\ \displaystyle \int_{\Omega} w_j w_i / dt + k (\nabla w_j. \nabla w_i ) + \int_{\Gamma_{13}} \alpha w_j w_i & \mbox{else if } i \not\in \Gamma_{24}, \mbox{or} j\ne i \end{array}\right.\\ M_{ij} &=& \left\{\begin{array}{cl} \frac{1}{\varepsilon} & \mbox{if } i \in \Gamma_{24}, \mbox{and} j=i\\ \displaystyle \int_{\Omega} w_j w_i / dt & \mbox{else if }i \not\in \Gamma_{24}, \mbox{or} j\ne i \end{array}\right. \\ b_{0,i} &=& \int_{\Gamma_{13}} \alpha u_{ue} w_i \\ b_{cl} &=& u^{0} \quad \mbox{the initial data} \end{eqnarray}
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
// Parameters
func fu0 = 10 + 90*x/6;
func k = 1.8*(y<0.5) + 0.2;
real ue = 25.;
real alpha = 0.25;
real T = 5;
real dt = 0.1 ;

// Mesh
mesh Th = square(30, 5, [6*x, y]);

// Fespace
fespace Vh(Th, P1);
Vh u0 = fu0, u = u0;

Create three variational formulation, and build the matrices A,M.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
// Problem
varf vthermic (u, v)
    = int2d(Th)(
          u*v/dt
        + k*(dx(u)*dx(v) + dy(u)*dy(v))
    )
    + int1d(Th, 1, 3)(
          alpha*u*v
    )
    + on(2, 4, u=1)
    ;

varf vthermic0 (u, v)
    = int1d(Th, 1, 3)(
          alpha*ue*v
    )
    ;

varf vMass (u, v)
    = int2d(Th)(
          u*v/dt
    )
    + on(2, 4, u=1)
    ;

real tgv = 1e30;
matrix A = vthermic(Vh, Vh, tgv=tgv, solver=CG);
matrix M = vMass(Vh, Vh);

Now, to build the right hand size we need 4 vectors.

1
2
3
4
real[int] b0 = vthermic0(0, Vh); //constant part of the RHS
real[int] bcn = vthermic(0, Vh); //tgv on Dirichlet boundary node ( !=0 )
//we have for the node i : i in Gamma_24 -> bcn[i] != 0
real[int] bcl = tgv*u0[]; //the Dirichlet boundary condition part

Note

The boundary condition is implemented by penalization and vector bcn contains the contribution of the boundary condition u=1, so to change the boundary condition, we have just to multiply the vector bc[] by the current value f of the new boundary condition term by term with the operator .*. Uzawa model gives a real example of using all this features.

And the new version of the algorithm is now:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
// Time loop
ofstream ff("thermic.dat");
for(real t = 0; t < T; t += dt){
    // Update
    real[int] b = b0; //for the RHS
    b += M*u[]; //add the the time dependent part
    //lock boundary part:
    b = bcn ? bcl : b; //do forall i: b[i] = bcn[i] ? bcl[i] : b[i]

    // Solve
    u[] = A^-1*b;

    // Save
    ff << t << " " << u(3, 0.5) << endl;

    // Plot
    plot(u);
}

// Display
for(int i = 0; i < 20; i++)
    cout << dy(u)(6.0*i/20.0, 0.9) << endl;

// Plot
plot(u, fill=true, wait=true);

Note

The functions appearing in the variational form are formal and local to the varf definition, the only important thing is the order in the parameter list, like in

1
2
varf vb1([u1, u2], q) = int2d(Th)((dy(u1) + dy(u2))*q) + int2d(Th)(1*q);
varf vb2([v1, v2], p) = int2d(Th)((dy(v1) + dy(v2))*p) + int2d(Th)(1*p);

To build matrix A from the bilinear part the variational form a of type varf simply write:

1
2
3
4
A = a(Vh, Wh , []...]);
// where
//Vh is "fespace" for the unknown fields with a correct number of component
//Wh is "fespace" for the test fields with a correct number of component

Possible named parameters in , [...] are

  • solver= LU, CG, Crout, Cholesky, GMRES, sparsesolver, UMFPACK ...

    The default solver is GMRES.

    The storage mode of the matrix of the underlying linear system depends on the type of solver chosen; for LU the matrix is sky-line non symmetric, for Crout the matrix is sky-line symmetric, for Cholesky the matrix is sky-line symmetric positive definite, for CG the matrix is sparse symmetric positive, and for GMRES, sparsesolver or UMFPACK the matrix is just sparse.

  • factorize = If true then do the matrix factorization for LU, Cholesky or Crout, the default value is false.

  • eps= A real expression. \varepsilon sets the stopping test for the iterative methods like CG. Note that if \varepsilon is negative then the stopping test is: if it is positive then the stopping test is

  • precon= Name of a function (for example P) to set the preconditioner.

    The prototype for the function P must be

    1
    func real[int] P(real[int] & xx) ;
    
  • tgv= Huge value (10^{30}) used to implement Dirichlet boundary conditions.

  • tolpivot= Set the tolerance of the pivot in UMFPACK (10^-1) and, LU, Crout, Cholesky factorization (10^{-20}).

  • tolpivotsym= Set the tolerance of the pivot sym in UMFPACK

  • strategy= Set the integer UMFPACK strategy (0 by default).

Note

The line of the matrix corresponding to the space Wh and the column of the matrix corresponding to the space Vh.

To build the dual vector b (of type real[int]) from the linear part of the variational form a do simply

1
2
real b(Vh.ndof);
b = a(0, Vh);

A first example to compute the area of each triangle K of mesh Th, just do:

1
2
3
4
fespace Nh(Th, P0); //the space function constant / triangle
Nh areaK;
varf varea (unused, chiK) = int2d(Th)(chiK);
etaK[] = varea(0, Ph);

Effectively, the basic functions of space Nh, are the characteristic function of the element of Th, and the numbering is the numeration of the element, so by construction:

Now, we can use this to compute error indicators like in example Adaptation using residual error indicator.

First to compute a continuous approximation to the function h "density mesh size" of the mesh Th.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
fespace Vh(Th, P1);
Vh h ;
real[int] count(Th.nv);
varf vmeshsizen (u, v) = intalledges(Th, qfnbpE=1)(v);
varf vedgecount (u, v) = intalledges(Th, qfnbpE=1)(v/lenEdge);

// Computation of the mesh size
count = vedgecount(0, Vh); //number of edge / vertex
h[] = vmeshsizen(0, Vh); //sum length edge / vertex
h[] = h[]./count; //mean length edge / vertex

To compute error indicator for Poisson equation : where h_K is size of the longest edge (hTriangle), h_e is the size of the current edge (lenEdge), n the normal.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
fespace Nh(Th, P0); // the space function constant / triangle
Nh etak;
varf vetaK (unused, chiK)
    = intalledges(Th)(
          chiK*lenEdge*square(jump(N.x*dx(u) + N.y*dy(u)))
    )
    + int2d(Th)(
          chiK*square(hTriangle*(f + dxx(u) + dyy(u)))
    )
    ;

 etak[] = vetaK(0, Ph);

We add automatic expression optimization by default, if this optimization creates problems, it can be removed with the keyword optimize as in the following example:

1
2
3
4
5
6
7
8
varf a (u1, u2)
    = int2d(Th, optimize=0)(
          dx(u1)*dx(u2)
        + dy(u1)*dy(u2)
    )
    + on(1, 2, 4, u1=0)
    + on(3, u1=1)
    ;

or you can also do optimization and remove the check by setting optimize=2.

Remark, it is all possible to build interpolation matrix, like in the following example:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
mesh TH = square(3, 4);
mesh th = square(2, 3);
mesh Th = square(4, 4);

fespace VH(TH, P1);
fespace Vh(th, P1);
fespace Wh(Th, P1);

matrix B = interpolate(VH, Vh); //build interpolation matrix Vh->VH
matrix BB = interpolate(Wh, Vh); //build interpolation matrix Vh->Wh

and after some operations on sparse matrices are available for example:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
int N = 10;
real [int, int] A(N, N); //a full matrix
real [int] a(N), b(N);
A = 0;
for (int i = 0; i < N; i++){
    A(i, i) = 1+i;
    if (i+1 < N) A(i, i+1) = -i;
    a[i] = i;
}
b = A*b;
matrix sparseA = A;
cout << sparseA << endl;
sparseA = 2*sparseA + sparseA';
sparseA = 4*sparseA + sparseA*5;
matrix sparseB = sparseA + sparseA + sparseA; ;
cout << "sparseB = " << sparseB(0,0) << endl;

Interpolation matrix#

It is also possible to store the matrix of a linear interpolation operator from a finite element space V_h to another W_h to interpolate(W_h,V_h,...) a function.

Note that the continuous finite functions are extended by continuity outside of the domain.

The named parameters of function interpolate are:

  • inside= set true to create zero-extension.
  • t= set true to get the transposed matrix
  • op= set an integer written below

    • 0 the default value and interpolate of the function
    • 1 interpolate the \p_x
    • 2 interpolate the \p_y
    • 3 interpolate the \p_z
  • U2Vc= set the which is the component of W_h come in V_h in interpolate process in a int array so the size of the array is number of component of W_h, if the put -1 then component is set to 0, like in the following example: (by default the component number is unchanged).

    1
    2
    3
    4
    5
    6
    7
    fespace V4h(Th4, [P1, P1, P1, P1]);
    fespace V3h(Th, [P1, P1, P1]);
    int[int] u2vc = [1, 3, -1]; //-1 -> put zero on the component
    matrix IV34 = interpolate(V3h, V4h, inside=0, U2Vc=u2vc); //V3h <- V4h
    V4h [a1, a2, a3, a4] = [1, 2, 3, 4];
    V3h [b1, b2, b3] = [10, 20, 30];
    b1[] = IV34*a1[];
    

    So here we have:

    1
    b1 == 2, b2 == 4, b3 == 0 ...
    

Matrix interpolation

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
// Mesh
mesh Th = square(4, 4);
mesh Th4 = square(2, 2, [x*0.5, y*0.5]);
plot(Th, Th4, wait=true);

// Fespace
fespace Vh(Th, P1);
Vh v, vv;
fespace Vh4(Th4, P1);
Vh4 v4=x*y;

fespace Wh(Th, P0);
fespace Wh4(Th4, P0);

// Interpolation
matrix IV = interpolate(Vh, Vh4); //here the function is exended by continuity
cout << "IV Vh<-Vh4 " << IV << endl;

v=v4;
vv[]= IV*v4[]; //here v == vv

real[int] diff= vv[] - v[];
cout << "|| v - vv || = " << diff.linfty << endl;
assert( diff.linfty<= 1e-6);

matrix IV0 = interpolate(Vh, Vh4, inside=1); //here the function is exended by zero
cout << "IV Vh<-Vh4 (inside=1) " << IV0 << endl;

matrix IVt0 = interpolate(Vh, Vh4, inside=1, t=1);
cout << "IV Vh<-Vh4^t (inside=1) " << IVt0 << endl;

matrix IV4t0 = interpolate(Vh4, Vh);
cout << "IV Vh4<-Vh^t " << IV4t0 << endl;

matrix IW4 = interpolate(Wh4, Wh);
cout << "IV Wh4<-Wh " << IW4 << endl;

matrix IW4V = interpolate(Wh4, Vh);
cout << "IV Wh4<-Vh " << IW4 << endl;

Build interpolation matrix A at a array of points (xx[j],yy[j]), i = 0, 2 here where w_i is the basic finite element function, c the component number, dop the type of diff operator like in op def.

1
2
3
4
5
6
7
real[int] xx = [.3, .4], yy = [.1, .4];
int c = 0, dop = 0;
matrix Ixx = interpolate(Vh, xx, yy, op=dop, composante=c);
cout << Ixx << endl;
Vh ww;
real[int] dd = [1, 2];
ww[] = Ixx*dd;

Schwarz

The following shows how to implement with an interpolation matrix a domain decomposition algorithm based on Schwarz method with Robin conditions.

Given a non-overlapping partition \bar\Omega=\bar\Omega_0\cup\bar\Omega_1 with \Omega_0\cap\Omega_1=\emptyset, \Sigma:=\bar\Omega_0\cap\bar\Omega_1 the algorithm is :

The same in variational form is: To discretized with the P^1 triangular Lagrangian finite element space V_h simply replace H^1_0(\Omega) by V_h(\Omega_0)\cup V_h(\Omega_1).

Then difficulty is to compute \int_{\Omega_j} \nabla u_j\cdot\nabla v when v is a basis function of V_h(\Omega_i), i\ne j.

It is done as follows (with \Gamma=\partial\Omega) :

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
// Parameters
int n = 30;
int Gamma = 1;
int Sigma = 2;

func f = 1.;
real alpha = 1.;

int Niter = 50;

// Mesh
mesh[int] Th(2);
int[int] reg(2);

border a0(t=0, 1){x=t; y=0; label=Gamma;}
border a1(t=1, 2){x=t; y=0; label=Gamma;}
border b1(t=0, 1){x=2; y=t; label=Gamma;}
border c1(t=2, 1){x=t; y=1; label=Gamma;}
border c0(t=1, 0){x=t; y=1; label=Gamma;}
border b0(t=1, 0){x=0; y=t; label=Gamma;}
border d(t=0, 1){x=1; y=t; label=Sigma;}
plot(a0(n) + a1(n) + b1(n) + c1(n) + c0(n) + b0(n) + d(n));
mesh TH = buildmesh(a0(n) + a1(n) + b1(n) + c1(n) + c0(n) + b0(n) + d(n));

reg(0) = TH(0.5, 0.5).region;
reg(1) = TH(1.5, 0.5).region;

for(int i = 0; i < 2; i++) Th[i] = trunc(TH, region==reg(i));

// Fespace
fespace Vh0(Th[0], P1);
Vh0 u0 = 0;

fespace Vh1(Th[1], P1);
Vh1 u1 = 0;

// Macro
macro grad(u) [dx(u), dy(u)] //

// Problem
int i;
varf a (u, v)
    = int2d(Th[i])(
          grad(u)'*grad(v)
    )
    + int1d(Th[i], Sigma)(
          alpha*u*v
    )
    + on(Gamma, u=0)
    ;

varf b (u, v)
    = int2d(Th[i])(
          f*v
    )
    + on(Gamma, u=0)
    ;

varf du1dn (u, v)
    =-int2d(Th[1])(
          grad(u1)'*grad(v)
        - f*v
    )
    + int1d(Th[1], Sigma)(
          alpha*u1*v
    )
    +on(Gamma, u=0)
    ;

varf du0dn (u, v)
    =-int2d(Th[0])(
          grad(u0)'*grad(v)
        - f*v
    )
    + int1d(Th[0], Sigma)(
          alpha*u0*v
    )
    +on(Gamma, u=0)
    ;

matrix I01 = interpolate(Vh1, Vh0);
matrix I10 = interpolate(Vh0, Vh1);

matrix[int] A(2);
i = 0; A[i] = a(Vh0, Vh0);
i = 1; A[i] = a(Vh1, Vh1);

// Solving loop
for(int iter = 0; iter < Niter; iter++){
    // Solve on Th[0]
    {
        i = 0;
        real[int] b0 = b(0, Vh0);
        real[int] Du1dn = du1dn(0, Vh1);
        real[int] Tdu1dn(Vh0.ndof); Tdu1dn = I01'*Du1dn;
        b0 += Tdu1dn;
        u0[] = A[0]^-1*b0;
    }
    // Solve on Th[1]
    {
        i = 1;
        real[int] b1 = b(0, Vh1);
        real[int] Du0dn = du0dn(0, Vh0);
        real[int] Tdu0dn(Vh1.ndof); Tdu0dn = I10'*Du0dn;
        b1 += Tdu0dn;
        u1[] = A[1]^-1*b1;
    }
    plot(u0, u1, cmm="iter="+iter);
}

Finite elements connectivity#

Here, we show how to get informations on a finite element space W_h({\cal T}_n,*), where "*" may be P1, P2, P1nc, etc.

  • Wh.nt gives the number of element of W_h
  • Wh.ndof gives the number of degrees of freedom or unknown
  • Wh.ndofK gives the number of degrees of freedom on one element
  • Wh(k,i) gives the number of ith degrees of freedom of element k.

See the following example:

Finite element connectivity

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
// Mesh
mesh Th = square(5, 5);

// Fespace
fespace Wh(Th, P2);

cout << "Number of degree of freedom = " << Wh.ndof << endl;
cout << "Number of degree of freedom / ELEMENT = " << Wh.ndofK << endl;

int k = 2, kdf = Wh.ndofK; //element 2
cout << "Degree of freedom of element " << k << ":" << endl;
for (int i = 0; i < kdf; i++)
    cout << Wh(k,i) << " ";
cout << endl;

The output is:

1
2
3
4
Number of degree of freedom = 121
Number of degree of freedom / ELEMENT = 6
Degree of freedom of element 2:
78 95 83 87 79 92

References#

[BERNADOU1980] BERNADOU, Michel, BOISSERIE, Jean-Marie, et HASSAN, Kamal. Sur l'implémentation des éléments finis de Hsieh-Clough-Tocher complet et réduit. 1980. Thèse de doctorat. INRIA.

[BERNARDI1985] BERNARDI, Christine et RAUGEL, Genevieve. Analysis of some finite elements for the Stokes problem. Mathematics of Computation, 1985, p. 71-79.

[THOMASSET2012] THOMASSET, François. Implementation of finite element methods for Navier-Stokes equations. Springer Science & Business Media, 2012.

[CROUZEIX1984] CROUZEIX, Michel et MIGNOT, Alain L. Analyse numérique des équations différentielles. Masson, 1984.

[TAYLOR2005] TAYLOR, Mark A., WINGATE, Beth A., et BOS, Len P. Several new quadrature formulas for polynomial integration in the triangle. arXiv preprint math/0501496, 2005.