Skip to content

Static problems

Soap Film#

Our starting point here will be the mathematical model to find the shape of soap film which is glued to the ring on the xy-plane

\begin{equation*} C=\{(x,y);\;x=\cos t,\,y=\sin t,\,0\leq t\leq 2\pi \}. \end{equation*}

We assume the shape of the film is described by the graph (x,y,u(x,y)) of the vertical displacement u(x,y)\, (x^2+y^2<1) under a vertical pressure p in terms of force per unit area and an initial tension \mu in terms of force per unit length.

Consider the "small plane" ABCD, A:(x,y,u(x,y)), B:(x,y,u(x+\delta x,y)), C:(x,y,u(x+\delta x,y+\delta y)) and D:(x,y,u(x,y+\delta y)).

Denote by \vec{n}(x,y)=(n_x(x,y),n_y(x,y),n_z(x,y)) the normal vector of the surface z=u(x,y). We see that the vertical force due to the tension \mu acting along the edge AD is -\mu n_x(x,y)\delta y and the the vertical force acting along the edge AD is

\mu n_x(x+\delta x,y)\delta y\simeq \mu\left(n_x(x,y) +\frac{\p n_x}{\p x}\delta x\right)(x,y)\delta y.
htbp

Similarly, for the edges AB and DC we have

-\mu n_y(x,y)\delta x,\qquad \mu\left(n_y(x,y)+\p n_y/\p y\right)(x,y)\delta x.

The force in the vertical direction on the surface ABCD due to the tension \mu is given by

\mu\left(\p n_x/\p x\right)\delta x\delta y+T\left(\p n_y/\p y\right)\delta y\delta x.

Assuming small displacements, we have

\begin{eqnarray*} \nu_x&=&(\p u/\p x)/\sqrt{1+(\p u/\p x)^2+(\p u/\p y)^2}\simeq \p u/\p x,\\ \nu_y&=&(\p u/\p y)/\sqrt{1+(\p u/\p x)^2+(\p u/\p y)^2}\simeq \p u/\p y. \end{eqnarray*}

Letting \delta x\to dx,\, \delta y\to dy, we have the equilibrium of the vertical displacement of soap film on ABCD by p

\mu dx dy\p^2 u/\p x^2 + \mu dx dy\p^2 u/\p y^2 + p dx dy = 0.

Using the Laplace operator \Delta = \p^2 /\p x^2 + \p^2 /\p y^2, we can find the virtual displacement write the following

\begin{equation} -\Delta u = f\quad \mbox{in }\Omega \end{equation}

where f=p/\mu, \Omega =\{(x,y);\;x^{2}+y^{2}<1\}.

Poisson's equation appears also in electrostatics taking the form of $f=\rho /\epsilon $ where $\rho $ is the charge density, $\epsilon $ the dielectric constant and u is named as electrostatic potential.

The soap film is glued to the ring % \p \Omega =C, then we have the boundary condition

\begin{equation} u=0\quad \mbox{on }\p \Omega \end{equation}

If the force is gravity, for simplify, we assume that f=-1.

 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
// Parameters
int nn = 50;
func f = -1;
func ue = (x^2+y^2-1)/4; //ue: exact solution

// Mesh
border a(t=0, 2*pi){x=cos(t); y=sin(t); label=1;}
mesh disk = buildmesh(a(nn));
plot(disk);

// Fespace
fespace femp1(disk, P1);
femp1 u, v;

// Problem
problem laplace (u, v)
    = int2d(disk)( //bilinear form
          dx(u)*dx(v)
        + dy(u)*dy(v)
    )
    - int2d(disk)( //linear form
          f*v
    )
    + on(1, u=0) //boundary condition
    ;

// Solve
laplace;

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

// Error
femp1 err = u - ue;
plot(err, value=true, wait=true);

cout << "error L2 = " << sqrt( int2d(disk)(err^2) )<< endl;
cout << "error H10 = " << sqrt( int2d(disk)((dx(u)-x/2)^2) + int2d(disk)((dy(u)-y/2)^2) )<< endl;

/// Re-run with a mesh adaptation ///

// Mesh adaptation
disk = adaptmesh(disk, u, err=0.01);
plot(disk, wait=true);

// Solve
laplace;
plot (u, value=true, wait=true);

// Error
err = u - ue; //become FE-function on adapted mesh
plot(err, value=true, wait=true);

cout << "error L2 = " << sqrt( int2d(disk)(err^2) )<< endl;
cout << "error H10 = " << sqrt( int2d(disk)((dx(u)-x/2)^2) + int2d(disk)((dy(u)-y/2)^2) )<< endl;
Fig. 1: isovalue of u Fig 2: A side view of u
SoapFilmSol SoapFilm3D

In the 37th line, the L^2-error estimation between the exact solution u_e, and in the following line, the H^1-error seminorm estimation are done on the initial mesh. The results are \|u_h - u_e\|_{0,\Omega}=0.000384045,\, |u_h - u_e|_{1,\Omega}=0.0375506.

After the adaptation, we have \|u_h - u_e\|_{0,\Omega}=0.000109043,\, |u_h - u_e|_{1,\Omega}=0.0188411. So the numerical solution is improved by adaptation of mesh.

Electrostatics#

We assume that there is no current and a time independent charge distribution. Then the electric field \mathbf{E} satisfies

\begin{eqnarray} \mathrm{div}\mathbf{E} &=& \rho/\epsilon\\ \mathrm{curl}\mathbf{E} &=& 0 \label{eqn::Maxwell} \end{eqnarray}

where \rho is the charge density and \epsilon is called the permittivity of free space.

From the equation \eqref{eqn::Maxwell} We can introduce the electrostatic potential such that \mathbf{E}=-\nabla \phi. Then we have Poisson's equation -\Delta \phi=f, f=-\rho/\epsilon.

We now obtain the equipotential line which is the level curve of \phi, when there are no charges except conductors \{C_i\}_{1,\cdots,K}. Let us assume K conductors C_1,\cdots,C_K within an enclosure C_0.

Each one is held at an electrostatic potential \varphi_i. We assume that the enclosure C0 is held at potential 0. In order to know \varphi(x) at any point x of the domain \Omega, we must solve

\begin{equation} -\Delta \varphi =0\quad \textrm{ in }\Omega , \end{equation}

where \Omega is the interior of C_0 minus the conductors C_i, and \Gamma is the boundary of \Omega, that is \sum_{i=0}^N C_i.

Here g is any function of x equal to \varphi_i on C_i and to 0 on C_0. The boundary equation is a reduced form for:

\begin{equation} \varphi =\varphi_{i}\;\text{on }C_{i},\;i=1...N,\varphi =0\;\text{on }C_{0}. \end{equation}

First we give the geometrical informations; C_0=\{(x,y);\; x^2+y^2=5^2\}, C_1=\{(x,y):\;\frac{1}{0.3^2}(x-2)^2+\frac{1}{3^2}y^2=1\}, C_2=\{(x,y):\; \frac{1}{0.3^2}(x+2)^2+\frac{1}{3^2}y^2=1\}.

Let \Omega be the disk enclosed by C_0 with the elliptical holes enclosed by C_1 and C_2. Note that C_0 is described counterclockwise, whereas the elliptical holes are described clockwise, because the boundary must be oriented so that the computational domain is to its left.

 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
// Mesh
border C0(t=0, 2*pi){x=5*cos(t); y=5*sin(t);}
border C1(t=0, 2*pi){x=2+0.3*cos(t); y=3*sin(t);}
border C2(t=0, 2*pi){x=-2+0.3*cos(t); y=3*sin(t);}

mesh Th = buildmesh(C0(60) + C1(-50) + C2(-50));
plot(Th);

// Fespace
fespace Vh(Th, P1);
Vh uh, vh;

// Problem
problem Electro (uh, vh)
    = int2d(Th)( //bilinear
          dx(uh)*dx(vh)
        + dy(uh)*dy(vh)
    )
    + on(C0, uh=0) //boundary condition on C_0
    + on(C1, uh=1) //+1 volt on C_1
    + on(C2, uh=-1) //-1 volt on C_2
    ;

// Solve
Electro;
plot(uh);
Fig. 3: Disk with two elliptical holes Fig. 9.4: Equipotential lines, where C_1 is located in right hand side
ElectrostaticsMesh Electrostatics

Aerodynamics#

Let us consider a wing profile S in a uniform flow. Infinity will be represented by a large circle \Gamma_{\infty}. As previously, we must solve

\begin{equation} \Delta \varphi=0\quad\textrm{in }\Omega, \quad \varphi|_S=c,\quad \varphi|_{\Gamma_{\infty}}=u_{\infty 1x}-u_{\infty2x} \label{eqn:NACA-5-5} \end{equation}

where \Omega is the area occupied by the fluid, u_{\infty} is the air speed at infinity, c is a constant to be determined so that \p_n\varphi is continuous at the trailing edge P of S (so-called Kutta-Joukowski condition). Lift is proportional to c.

To find c we use a superposition method. As all equations in \eqref{eqn:NACA-5-5} are linear, the solution \varphi_c is a linear function of c

\begin{equation} \varphi_c = \varphi_0 + c\varphi_1, \end{equation}

where \varphi_0 is a solution of \eqref{eqn:NACA-5-5} with c = 0 and \varphi_1 is a solution with c = 1 and zero speed at infinity.

With these two fields computed, we shall determine c by requiring the continuity of \p \varphi /\p n at the trailing edge. An equation for the upper surface of a NACA0012 (this is a classical wing profile in aerodynamics; the rear of the wing is called the trailing edge) is:

\begin{equation} y = 0.17735\sqrt{x} - 0.075597x - 0.212836x^2 + 0.17363x^3 - 0.06254x^4. \end{equation}

Taking an incidence angle \alpha such that \tan \alpha = 0.1, we must solve

\begin{equation} -\Delta\varphi = 0\qquad \textrm{in }\Omega, \qquad \varphi|_{\Gamma_1} = y - 0.1x,\quad \varphi |_{\Gamma_2} = c, \end{equation}

where \Gamma_2 is the wing profile and \Gamma_1 is an approximation of infinity. One finds c by solving:

\begin{eqnarray} -\Delta\varphi_0 = 0 ~~\textrm{in }\Omega,\qquad \varphi_0|_{\Gamma_1} = y - 0.1x, \quad \varphi_0|_{\Gamma_2} = 0,\\ -\Delta\varphi_1 = 0 ~~\textrm{in }\Omega, \qquad \varphi_1|_{\Gamma_1} = 0, \quad \varphi_1|_{\Gamma_2} = 1. \end{eqnarray}

The solution \varphi = \varphi_0+c\varphi_1 allows us to find c by writing that \p_n\varphi has no jump at the trailing edge P = (1, 0).

We have \p n\varphi -(\varphi (P^+)-\varphi (P))/\delta where P^+ is the point just above P in the direction normal to the profile at a distance \delta. Thus the jump of \p_n\varphi is (\varphi_0|_{P^+} +c(\varphi_1|_{P^+} -1))+(\varphi_0|_{P^-} +c(\varphi_1|_{P^-} -1)) divided by \delta because the normal changes sign between the lower and upper surfaces. Thus

\begin{equation} c = -\frac{\varphi_0|_{P^+} + \varphi_0|_{P^-}} {(\varphi_1|_{P^+} + \varphi_1|_{P^-} - 2)} , \end{equation}

which can be programmed as:

\begin{equation} c = -\frac{\varphi_0(0.99, 0.01) + \varphi_0(0.99,-0.01)} {(\varphi_1(0.99, 0.01) + \varphi_1(0.99,-0.01) - 2)} . \end{equation}
 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
// Mesh
border a(t=0, 2*pi){x=5*cos(t); y=5*sin(t);}
border upper(t=0, 1) {
    x=t;
    y=0.17735*sqrt(t)-0.075597*t - 0.212836*(t^2) + 0.17363*(t^3) - 0.06254*(t^4);
}
border lower(t=1, 0) {
x=t;
y=-(0.17735*sqrt(t) - 0.075597*t - 0.212836*(t^2) + 0.17363*(t^3) - 0.06254*(t^4));
}
border c(t=0, 2*pi){x=0.8*cos(t)+0.5; y=0.8*sin(t);}

mesh Zoom = buildmesh(c(30) + upper(35) + lower(35));
mesh Th = buildmesh(a(30) + upper(35) + lower(35));

// Fespace
fespace Vh(Th, P2);
Vh psi0, psi1, vh;

fespace ZVh(Zoom, P2);

// Problem
solve Joukowski0(psi0, vh)
    = int2d(Th)(
          dx(psi0)*dx(vh)
        + dy(psi0)*dy(vh)
    )
    + on(a, psi0=y-0.1*x)
    + on(upper, lower, psi0=0)
    ;

plot(psi0);

solve Joukowski1(psi1,vh)
    = int2d(Th)(
          dx(psi1)*dx(vh)
        + dy(psi1)*dy(vh)
    )
    + on(a, psi1=0)
    + on(upper, lower, psi1=1);

plot(psi1);

//continuity of pressure at trailing edge
real beta = psi0(0.99,0.01) + psi0(0.99,-0.01);
beta = -beta / (psi1(0.99,0.01) + psi1(0.99,-0.01)-2);

Vh psi = beta*psi1 + psi0;
plot(psi);

ZVh Zpsi = psi;
plot(Zpsi, bw=true);

ZVh cp = -dx(psi)^2 - dy(psi)^2;
plot(cp);

ZVh Zcp = cp;
plot(Zcp, nbiso=40);
Fig. 5: isovalue of cp = -(\p_x\psi)^2 - (\p_y\psi)^2 Fig. 6: Zooming of cp
Aerodynamics1 Aerodynamics2

Error estimation#

There are famous estimation between the numerical result u_h and the exact solution u of the Poisson's problem:

If triangulations \{\mathcal{T}_h\}_{h\downarrow 0} is regular (see Regular Triangulation), then we have the estimates

\begin{eqnarray} |\nabla u - \nabla u_h|_{0,\Omega} &\le& C_1h \label{eqn::ErrorEstimatation1}\\ \|u - u_h\|_{0,\Omega} &\le& C_2h^2 \label{eqn::ErrorEstimatation2} \end{eqnarray}

with constants C_1,\, C_2 independent of h, if u is in H^2(\Omega). It is known that u\in H^2(\Omega) if \Omega is convex.

In this section we check \eqref{eqn::ErrorEstimatation1} and \eqref{eqn::ErrorEstimatation2}. We will pick up numericall error if we use the numerical derivative, so we will use the following for \eqref{eqn::ErrorEstimatation1}.

\begin{eqnarray*} \int_{\Omega}|\nabla u - \nabla u_h|^2\, \d x\d y &=&\int_{\Omega}\nabla u\cdot \nabla(u - 2u_h)\, \d x\d y+ \int_{\Omega}\nabla u_h\cdot \nabla u_h\, \d x\d y\\ &=&\int_{\Omega}f(u-2u_h)\, \d x\d y+\int_{\Omega}fu_h\, \d x\d y \end{eqnarray*}

The constants C_1,\, C_2 are depend on \mathcal{T}_h and f, so we will find them by FreeFem++.

In general, we cannot get the solution u as a elementary functions even if spetical functions are added. Instead of the exact solution, here we use the approximate solution u_0 in V_h(\mathcal{T}_h,P_2),\, h\sim 0.

 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
// Parameters
func f = x*y;

//Mesh
mesh Th0 = square(100, 100);

// Fespace
fespace V0h(Th0, P2);
V0h u0, v0;

// Problem
solve Poisson0 (u0, v0)
    = int2d(Th0)(
          dx(u0)*dx(v0)
        + dy(u0)*dy(v0)
    )
    - int2d(Th0)(
          f*v0
    )
    + on(1, 2, 3, 4, u0=0)
    ;
plot(u0);

// Error loop
real[int] errL2(10), errH1(10);
for (int i = 1; i <= 10; i++){
    // Mesh
    mesh Th = square(5+i*3,5+i*3);

    // Fespace
    fespace Vh(Th, P1);
    Vh u, v;
    fespace Ph(Th, P0);
    Ph h = hTriangle; //get the size of all triangles

    // Problem
    solve Poisson (u, v)
        = int2d(Th)(
              dx(u)*dx(v)
            + dy(u)*dy(v)
        )
        - int2d(Th)(
              f*v
        )
        + on(1, 2, 3, 4, u=0)
        ;

    // Error
    V0h uu = u; //interpolate solution on first mesh
    errL2[i-1] = sqrt( int2d(Th0)((uu - u0)^2) )/h[].max^2;
    errH1[i-1] = sqrt( int2d(Th0)(f*(u0 - 2*uu + uu)) )/h[].max;
}

// Display
cout << "C1 = " << errL2.max << "("<<errL2.min<<")" << endl;
cout << "C2 = " << errH1.max << "("<<errH1.min<<")" << endl;

We can guess that C_1=0.0179253(0.0173266) and C_2=0.0729566(0.0707543), where the numbers inside the parentheses are minimum in calculation.

Periodic Boundary Conditions#

We now solve the Poisson equation on the square ]0,2\pi[^2 under bi-periodic boundary condition u(0,y)=u(2\pi,y) for all y and u(x,0)=u(x,2\pi) for all x.

These boundary conditions are achieved from the definition of the periodic finite element space.

 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
// Parameters
func f = sin(x+pi/4.)*cos(y+pi/4.); //right hand side

// Mesh
mesh Th = square(10, 10, [2*x*pi, 2*y*pi]);

// Fespace
//defined the fespace with periodic condition
//label: 2 and 4 are left and right side with y the curve abscissa
//       1 and 2 are bottom and upper side with x the curve abscissa
fespace Vh(Th, P2, periodic=[[2, y], [4, y], [1, x], [3, x]]);
Vh uh, vh;

// Problem
problem laplace (uh, vh)
    = int2d(Th)(
          dx(uh)*dx(vh)
        + dy(uh)*dy(vh)
    )
    + int2d(Th)(
        - f*vh
    )
    ;

// Solve
laplace;

// Plot
plot(uh, value=true);
Fig. 7: The isovalue of solution u with periodic boundary condition
PeriodicBoundaryConditions

The periodic condition does not necessarily require parallel boundaries. The following example give such example.

Periodic boundary conditions - non-parallel boundaries

 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
// Parameters
int n = 10;
real r = 0.25;
real r2 = 1.732;
func f = (y+x+1)*(y+x-1)*(y-x+1)*(y-x-1);

// Mesh
border a(t=0, 1){x=-t+1; y=t; label=1;};
border b(t=0, 1){x=-t; y=1-t; label=2;};
border c(t=0, 1){x=t-1; y=-t; label=3;};
border d(t=0, 1){x=t; y=-1+t; label=4;};
border e(t=0, 2*pi){x=r*cos(t); y=-r*sin(t); label=0;};
mesh Th = buildmesh(a(n) + b(n) + c(n) + d(n) + e(n));
plot(Th, wait=true);

// Fespace
//warning for periodic condition:
//side a and c
//on side a (label 1) $ x \in [0,1] $ or $ x-y\in [-1,1] $
//on side c (label 3) $ x \in [-1,0]$ or $ x-y\in[-1,1] $
//so the common abscissa can be respectively $x$ and $x+1$
//or you can can try curviline abscissa $x-y$ and $x-y$
//1 first way
//fespace Vh(Th, P2, periodic=[[2, 1+x], [4, x], [1, x], [3, 1+x]]);
//2 second way
fespace Vh(Th, P2, periodic=[[2, x+y], [4, x+y], [1, x-y], [3, x-y]]);
Vh uh, vh;

// Problem
real intf = int2d(Th)(f);
real mTh = int2d(Th)(1);
real k =  intf / mTh;
problem laplace (uh, vh)
    = int2d(Th)(
          dx(uh)*dx(vh)
        + dy(uh)*dy(vh)
        )
    + int2d(Th)(
          (k-f)*vh
    )
    ;

// Solve
laplace;

// Plot
plot(uh, wait=true);
Fig. 8: The isovalue of solution u for \Delta u = ((y+x)^{2}+1)((y-x)^{2}+1) - k, in \Omega and \p_{n} u =0 on hole, and with two periodic boundary condition on external border
PeriodicBoundaryConditions2

An other example with no equal border, just to see if the code works.

Periodic boundary conditions - non-equal border

 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
// Macro
//irregular boundary condition to build border AB
macro LINEBORDER(A, B, lab)
    border A#B(t=0,1){
        real t1=1.-t;
        x=A#x*t1+B#x*t;
        y=A#y*t1+B#y*t;
        label=lab;
} //EOM
// compute ||AB|| A=(ax,ay) et B =(bx,by)
macro dist(ax, ay, bx, by) sqrt(square((ax)-(bx)) + square((ay)-(by))) //EOM
macro Grad(u) [dx(u), dy(u)] //EOM

// Parameters
int n = 10;
real Ax = 0.9, Ay = 1;
real Bx = 2, By = 1;
real Cx = 2.5, Cy = 2.5;
real Dx = 1, Dy = 2;
real gx = (Ax+Bx+Cx+Dx)/4.;
real gy = (Ay+By+Cy+Dy)/4.;

// Mesh
LINEBORDER(A,B,1)
LINEBORDER(B,C,2)
LINEBORDER(C,D,3)
LINEBORDER(D,A,4)
mesh Th=buildmesh(AB(n)+BC(n)+CD(n)+DA(n),fixedborder=1);

// Fespace
real l1 = dist(Ax,Ay,Bx,By);
real l2 = dist(Bx,By,Cx,Cy);
real l3 = dist(Cx,Cy,Dx,Dy);
real l4 = dist(Dx,Dy,Ax,Ay);
func s1 = dist(Ax,Ay,x,y)/l1; //absisse on AB = ||AX||/||AB||
func s2 = dist(Bx,By,x,y)/l2; //absisse on BC = ||BX||/||BC||
func s3 = dist(Cx,Cy,x,y)/l3; //absisse on CD = ||CX||/||CD||
func s4 = dist(Dx,Dy,x,y)/l4; //absisse on DA = ||DX||/||DA||
verbosity = 6; //to see the abscisse value of the periodic condition
fespace Vh(Th, P1, periodic=[[1, s1], [3, s3], [2, s2], [4, s4]]);
verbosity = 1; //reset verbosity
Vh u, v;

real cc = 0;
cc = int2d(Th)((x-gx)*(y-gy)-cc)/Th.area;
cout << "compatibility = " << int2d(Th)((x-gx)*(y-gy)-cc) <<endl;

// Problem
solve Poisson (u, v)
    = int2d(Th)(
          Grad(u)'*Grad(v)
        + 1e-10*u*v
    )
    -int2d(Th)(
          10*v*((x-gx)*(y-gy)-cc)
    )
    ;

// Plot
plot(u, value=true);

Periodic boundry conditions - Poisson cube-balloon

 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
load "msh3"
load "tetgen"
load "medit"

// Parameters
real hs = 0.1; //mesh size on sphere
int[int] N = [20, 20, 20];
real [int,int] B = [[-1, 1], [-1, 1], [-1, 1]];
int [int,int] L = [[1, 2], [3, 4], [5, 6]];

real x0 = 0.3, y0 = 0.4, z0 = 06;
func f = sin(x*2*pi+x0)*sin(y*2*pi+y0)*sin(z*2*pi+z0);

// Mesh
bool buildTh = 0;
mesh3 Th;
try { //a way to build one time the mesh or read it if the file exist
    Th = readmesh3("Th-hex-sph.mesh");
}
catch (...){
    buildTh = 1;
}

if (buildTh){
    include "MeshSurface.idp"

    // Surface Mesh
    mesh3 ThH = SurfaceHex(N, B, L, 1);
    mesh3 ThS = Sphere(0.5, hs, 7, 1);

    mesh3 ThHS = ThH + ThS;

    real voltet = (hs^3)/6.;
    real[int] domain = [0, 0, 0, 1, voltet, 0, 0, 0.7, 2, voltet];
    Th = tetg(ThHS, switch="pqaAAYYQ", nbofregions=2, regionlist=domain);

    savemesh(Th, "Th-hex-sph.mesh");
}

// Fespace
fespace Ph(Th, P0);
Ph reg = region;
cout << " centre = " << reg(0,0,0) << endl;
cout << " exterieur = " << reg(0,0,0.7) << endl;

verbosity = 50;
fespace Vh(Th, P1, periodic=[[3, x, z], [4, x, z], [1, y, z], [2, y, z], [5, x, y], [6, x, y]]);
verbosity = 1;
Vh uh,vh;

// Macro
macro Grad(u) [dx(u),dy(u),dz(u)] // EOM

// Problem
problem Poisson (uh, vh)
    = int3d(Th, 1)(
          Grad(uh)'*Grad(vh)*100
    )
    + int3d(Th, 2)(
          Grad(uh)'*Grad(vh)*2
    )
    + int3d(Th)(
          vh*f
    )
;

// Solve
Poisson;

// Plot
plot(uh, wait=true, nbiso=6);
medit("uh", Th, uh);
Fig. 9: View of the surface isovalue of periodic solution uh Fig. 10: View a the cut of the solution uh with ffmedit
PeriodicBoundaryConditionsPoisson1 PeriodicBoundaryConditionsPoisson2

Poisson Problems with mixed boundary condition#

Here we consider the Poisson equation with mixed boundary conditions: For given functions f and g, find u such that

\begin{eqnarray} -\Delta u &=& f & \textrm{ in }\Omega\\ u &=& g &\textrm{ on }\Gamma_D\\ \p u/\p n &=& 0 &\textrm{ on }\Gamma_N \end{eqnarray}

where \Gamma_D is a part of the boundary \Gamma and \Gamma_N=\Gamma\setminus \overline{\Gamma_D}.

The solution u has the singularity at the points \{\gamma_1,\gamma_2\}=\overline{\Gamma_D}\cap\overline{\Gamma_N}.

When \Omega=\{(x,y);\; -1<x<1,\, 0<y<1\}, \Gamma_N=\{(x,y);\; -1\le x<0,\, y=0\}, \Gamma_D=\p \Omega\setminus \Gamma_N, the singularity will appear at \gamma_1=(0,0),\, \gamma_2(-1,0), and u has the expression with a constants K_i.

Here u_S = r_j^{1/2}\sin(\theta_j/2) by the local polar coordinate (r_j,\theta_j at \gamma_j such that (r_1,\theta_1)=(r,\theta).

Instead of polar coordinate system (r,\theta), we use that r=sqrt(x^2+y^2) and \theta =atan2(y,x) in FreeFem++.

Assume that f=-2\times 30(x^2+y^2) and g=u_e=10(x^2+y^2)^{1/4}\sin\left([\tan^{-1}(y/x)]/2\right)+30(x^2y^2), where u_eS is the exact solution.

 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
// Parameters
func f = -2*30*(x^2+y^2); //given function
//the singular term of the solution is K*us (K: constant)
func us = sin(atan2(y,x)/2)*sqrt( sqrt(x^2+y^2) );
real K = 10.;
func ue = K*us + 30*(x^2*y^2);

// Mesh
border N(t=0, 1){x=-1+t; y=0; label=1;};
border D1(t=0, 1){x=t; y=0; label=2;};
border D2(t=0, 1){x=1; y=t; label=2;};
border D3(t=0, 2){x=1-t; y=1; label=2;};
border D4(t=0, 1){x=-1; y=1-t; label=2;};

mesh T0h = buildmesh(N(10) + D1(10) + D2(10) + D3(20) + D4(10));
plot(T0h, wait=true);

// Fespace
fespace V0h(T0h, P1);
V0h u0, v0;

//Problem
solve Poisson0 (u0, v0)
    = int2d(T0h)(
          dx(u0)*dx(v0)
        + dy(u0)*dy(v0)
    )
    - int2d(T0h)(
          f*v0
    )
    + on(2, u0=ue)
    ;

// Mesh adaptation by the singular term
mesh Th = adaptmesh(T0h, us);
for (int i = 0; i < 5; i++)
    mesh Th = adaptmesh(Th, us);

// Fespace
fespace Vh(Th, P1);
Vh u, v;

// Problem
solve Poisson (u, v)
    = int2d(Th)(
          dx(u)*dx(v)
        + dy(u)*dy(v)
    )
    - int2d(Th)(
          f*v
    )
    + on(2, u=ue)
    ;

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

// Error in H1 norm
Vh uue = ue;
real H1e = sqrt( int2d(Th)(dx(uue)^2 + dy(uue)^2 + uue^2) );
Vh err0 = u0 - ue;
Vh err = u - ue;
Vh H1err0 = int2d(Th)(dx(err0)^2 + dy(err0)^2 + err0^2);
Vh H1err = int2d(Th)(dx(err)^2 + dy(err)^2 + err^2);
cout << "Relative error in first mesh = "<< int2d(Th)(H1err0)/H1e << endl;
cout << "Relative error in adaptive mesh = "<< int2d(Th)(H1err)/H1e << endl;

From line 35 to 37, mesh adaptations are done using the base of singular term.

In line 61, H1e=\|u_e\|_{1,\Omega} is calculated.

In lines 64 and 65, the relative errors are calculated, that is,

\begin{eqnarray*} \|u^0_h-u_e\|_{1,\Omega}/H1e&=&0.120421\\ \|u^a_h-u_e\|_{1,\Omega}/H1e&=&0.0150581 \end{eqnarray*}

where u^0_h is the numerical solution in T0h and u^a_h is u in this program.

Poisson with mixed finite element#

Here we consider the Poisson equation with mixed boundary value problems:

For given functions f , g_d, g_n, find p such that

\begin{eqnarray} -\Delta p &=& 1 & \textrm{ in }\Omega\\ p &=& g_d & \textrm{ on }\Gamma_D\\ \p p/\p n &=& g_n & \textrm{ on }\Gamma_N \end{eqnarray}

where \Gamma_D is a part of the boundary \Gamma and \Gamma_N=\Gamma\setminus \overline{\Gamma_D}.

The mixed formulation is: find p and \mathbf{u} such that

\begin{eqnarray} \nabla p + \mathbf{u} &=& \mathbf{0} & \textrm{ in }\Omega\\ \nabla. \mathbf{u} &=& f & \textrm{ in }\Omega\\ p &=& g_d & \textrm{ on }\Gamma_D\\ \p u. n &=& \mathbf{g}_n.n & \textrm{ on }\Gamma_N \end{eqnarray}

where \mathbf{g}_n is a vector such that \mathbf{g}_n.n = g_n.

The variational formulation is,

\begin{eqnarray} \forall \mathbf{v} \in \mathbb{V}_0: & \int_\Omega p \nabla.v + \mathbf{v} \mathbf{v} &= \int_{\Gamma_d} g_d \mathbf{v}.n\\ \forall {q} \in \mathbb{P}: & \int_\Omega q \nabla.u &= \int_\Omega q f\nonumber\\ & \p u. n &= \mathbf{g}_n.n \quad \textrm{on }\Gamma_N \end{eqnarray}

where the functional space are:

\mathbb{P}= L^2(\Omega), \qquad\mathbb{V}= H(div)=\{\mathbf{v}\in L^2(\Omega)^2,\nabla.\mathbf{v}\in L^2(\Omega)\}

and

\mathbb{V}_0 = \{\mathbf{v}\in \mathbb{V};\quad\mathbf{v}. n = 0 \quad\mathrm{on }\;\;\Gamma_N\}

To write the FreeFem++ example, we have just to choose the finites elements spaces.

Here \mathbb{V} space is discretize with Raviart-Thomas finite element RT0 and \mathbb{P} is discretize by constant finite element P0.

Example 9.10 LaplaceRT.edp

 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
// Parameters
func gd = 1.;
func g1n = 1.;
func g2n = 1.;

// Mesh
mesh Th = square(10, 10);

// Fespace
fespace Vh(Th, RT0);
Vh [u1, u2];
Vh [v1, v2];

fespace Ph(Th, P0);
Ph p, q;

// Problem
problem laplaceMixte ([u1, u2, p], [v1, v2, q], solver=GMRES, eps=1.0e-10, tgv=1e30, dimKrylov=150)
    = int2d(Th)(
          p*q*1e-15 //this term is here to be sure
                    // that all sub matrix are inversible (LU requirement)
        + u1*v1
        + u2*v2
        + p*(dx(v1)+dy(v2))
        + (dx(u1)+dy(u2))*q
    )
    + int2d(Th) (
          q
    )
    - int1d(Th, 1, 2, 3)(
          gd*(v1*N.x +v2*N.y)
    )
    + on(4, u1=g1n, u2=g2n)
    ;

// Solve
laplaceMixte;

// Plot
plot([u1, u2], coef=0.1, wait=true, value=true);
plot(p, fill=1, wait=true, value=true);

Metric Adaptation and residual error indicator#

We do metric mesh adaption and compute the classical residual error indicator \eta_{T} on the element T for the Poisson problem.

First, we solve the same problem as in a previous example.

 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
// Parameters
real[int] viso(21);
for (int i = 0; i < viso.n; i++)
    viso[i] = 10.^(+(i-16.)/2.);
real error = 0.01;
func f = (x-y);

// Mesh
border ba(t=0, 1.0){x=t; y=0; label=1;}
border bb(t=0, 0.5){x=1; y=t; label=2;}
border bc(t=0, 0.5){x=1-t; y=0.5; label=3;}
border bd(t=0.5, 1){x=0.5; y=t; label=4;}
border be(t=0.5, 1){x=1-t; y=1; label=5;}
border bf(t=0.0, 1){x=0; y=1-t; label=6;}
mesh Th = buildmesh(ba(6) + bb(4) + bc(4) + bd(4) + be(4) + bf(6));

// Fespace
fespace Vh(Th, P2);
Vh u, v;

fespace Nh(Th, P0);
Nh rho;

// Problem
problem Probem1 (u, v, solver=CG, eps=1.0e-6)
    = int2d(Th, qforder=5)(
          u*v*1.0e-10
        + dx(u)*dx(v)
        + dy(u)*dy(v)
    )
    + int2d(Th, qforder=5)(
        - f*v
    )
    ;

Now, the local error indicator \eta_{T} is:

\eta_{T} =\left( h_{T}^{2} || f + \Delta u_{{h}} ||_{L^{2}(T)}^{2} +\sum_{e\in \mathcal{E}_{K}} h_{e} \,||\, [ \frac{\p u_{h}}{\p n_{k}}] \,||^{2}_{L^{2}(e)} \right)^{\frac{1}{2}}

where h_{T} is the longest edge of T, {\cal E}_T is the set of T edge not on \Gamma=\p \Omega, n_{T} is the outside unit normal to K, h_{e} is the length of edge e, [ g ] is the jump of the function g across edge (left value minus right value).

Of course, we can use a variational form to compute \eta_{T}^{2}, with test function constant function in each triangle.

 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
// Error
varf indicator2 (uu, 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)))
    )
    ;

// Mesh adaptation loop
for (int i = 0; i < 4; i++){
    // Solve
    Probem1;
    cout << u[].min << " " << u[].max << endl;
    plot(u, wait=true);

    // Error
    rho[] = indicator2(0, Nh);
    rho = sqrt(rho);
    cout << "rho = min " << rho[].min << " max=" << rho[].max << endl;
    plot(rho, fill=true, wait=true, cmm="indicator density", value=true, viso=viso, nbiso=viso.n);

    // Mesh adaptation
    plot(Th, wait=true, cmm="Mesh (before adaptation)");
    Th = adaptmesh(Th, [dx(u), dy(u)], err=error, anisomax=1);
    plot(Th, wait=true, cmm="Mesh (after adaptation)");
    u = u;
    rho = rho;
    error = error/2;
}

If the method is correct, we expect to look the graphics by an almost constant function \eta on your computer as in Fig. 11.

Fig. 11: Density of the error indicator with isotropic P_{2} metric
MetricAdaptation
MetricAdaptation2

Adaptation using residual error indicator#

In the previous example we compute the error indicator, now we use it, to adapt the mesh. The new mesh size is given by the following formulae:

h_{n+1}(x) = \frac{h_{n}(x)}{f_{n}(\eta_K(x))}

where \eta_n(x) is the level of error at point x given by the local error indicator, h_n is the previous "mesh size" field, and f_n is a user function define by f_n = min(3,max(1/3,\eta_n / \eta_n^* )) where $\eta_n^* =mean(\eta_n) c $, and c is an user coefficient generally close to one.

First a macro MeshSizecomputation is defined to get a P_1 mesh size as the average of edge length.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
// macro the get the current mesh size parameter
// in:
// Th the mesh
// Vh P1 fespace on Th
// out :
// h: the Vh finite element finite set to the current mesh size
macro MeshSizecomputation (Th, Vh, h)
{
    real[int] count(Th.nv);
    /*mesh size (lenEdge = integral(e) 1 ds)*/
    varf vmeshsizen (u, v) = intalledges(Th, qfnbpE=1)(v);
    /*number of edges per vertex*/
    varf vedgecount (u, v) = intalledges(Th, qfnbpE=1)(v/lenEdge);
    /*mesh size*/
    count = vedgecount(0, Vh);
    h[] = 0.;
    h[] = vmeshsizen(0, Vh);
    cout << "count min = " << count.min << " max = " << count.max << endl;
    h[] = h[]./count;
    cout << "-- bound meshsize = " << h[].min << " " << h[].max << endl;
} //

A second macro to re-mesh according to the new mesh size.

 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
// macro to remesh according the de residual indicator
// in:
// Th the mesh
// Ph P0 fespace on Th
// Vh P1 fespace on Th
// vindicator the varf to evaluate the indicator
// coef on etameam
macro ReMeshIndicator (Th, Ph, Vh, vindicator, coef)
{
    Vh h=0;
    /*evaluate the mesh size*/
    MeshSizecomputation(Th, Vh, h);
    Ph etak;
    etak[] = vindicator(0, Ph);
    etak[] = sqrt(etak[]);
    real etastar= coef*(etak[].sum/etak[].n);
    cout << "etastar = " << etastar << " sum = " << etak[].sum << " " << endl;

    /*etaK is discontinous*/
    /*we use P1 L2 projection with mass lumping*/
    Vh fn, sigma;
    varf veta(unused, v) = int2d(Th)(etak*v);
    varf vun(unused, v) = int2d(Th)(1*v);
    fn[] = veta(0, Vh);
    sigma[] = vun(0, Vh);
    fn[] = fn[]./ sigma[];
    fn = max(min(fn/etastar,3.),0.3333);

    /*new mesh size*/
    h = h / fn;
    /*build the mesh*/
    Th = adaptmesh(Th, IsMetric=1, h, splitpbedge=1, nbvx=10000);
} //
 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
// Parameters
real hinit = 0.2; //initial mesh size
func f=(x-y);

// Mesh
border ba(t=0, 1.0){x=t; y=0; label=1;}
border bb(t=0, 0.5){x=1; y=t; label=2;}
border bc(t=0, 0.5){x=1-t; y=0.5; label=3;}
border bd(t=0.5, 1){x=0.5; y=t; label=4;}
border be(t=0.5, 1){x=1-t; y=1; label=5;}
border bf(t=0.0, 1){x=0; y=1-t; label=6;}
mesh Th = buildmesh(ba(6) + bb(4) + bc(4) + bd(4) + be(4) + bf(6));

// Fespace
fespace Vh(Th, P1); //for the mesh size and solution
Vh h = hinit; //the FE function for the mesh size
Vh u, v;

fespace Ph(Th, P0); //for the error indicator

//Build a mesh with the given mesh size hinit
Th = adaptmesh(Th, h, IsMetric=1, splitpbedge=1, nbvx=10000);
plot(Th, wait=1);

// Problem
problem Poisson (u, v)
    = int2d(Th, qforder=5)(
          u*v*1.0e-10
        + dx(u)*dx(v)
        + dy(u)*dy(v)
    )
    - int2d(Th, qforder=5)(
          f*v
    )
    ;

 varf indicator2 (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)))
    )
    ;

// Mesh adaptation loop
for (int i = 0; i < 10; i++){
    u = u;

    // Solve
    Poisson;
    plot(Th, u, wait=true);

    real cc = 0.8;
    if (i > 5) cc=1;
    ReMeshIndicator(Th, Ph, Vh, indicator2, cc);
    plot(Th, wait=true);
}
Fig. 9.12: The error indicator with isotropic P_{1}, the mesh and isovalue of the solution
AdaptationResidualError
AdaptationResidualError2