FreeFEM Documentation on GitHub

stars - forks

# 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:

$C=\{(x,y);\;x=\cos t,\,y=\sin t,\,0\leq t\leq 2\pi \}$

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$

Similarly, for the edges AB and DC we have:

$-\mu n_y(x,y)\delta x,\quad\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{split}\begin{array}{rcccl} \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{array}\end{split}$

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:

$-\Delta u = f\quad \mbox{in }\Omega$

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:

$u=0\quad \mbox{on }\p \Omega$

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; 

In the 37th line, the $$L^2$$-error estimation between the exact solution $$u_e$$,

$\|u_h - u_e\|_{0,\Omega}=\left(\int_{\Omega}|u_h-u_e|^2\, \d x\d y\right)^{1/2}$

and in the following line, the $$H^1$$-error seminorm estimation:

$|u_h - u_e|_{1,\Omega}=\left(\int_{\Omega}|\nabla u_h-\nabla u_e|^2\, \d x\d y\right)^{1/2}$

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:

(42) $\begin{split}\begin{array}{rcl} \mathrm{div}\mathbf{E} &=& \rho/\epsilon\\ \mathrm{curl}\mathbf{E} &=& 0 \end{array}\end{split}$

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

From the equation (42) 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:

$-\Delta \varphi =0\quad \textrm{ in }\Omega$

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:

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

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. 138 Equipotential lines where $$C_1$$ is located in right hand side

## 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:

(43) $\Delta \varphi=0\quad\textrm{in }\Omega, \quad \varphi|_S=c,\quad \varphi|_{\Gamma_{\infty}}=u_{\infty 1x}-u_{\infty2x}$

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 (43) are linear, the solution $$\varphi_c$$ is a linear function of $$c$$

$\varphi_c = \varphi_0 + c\varphi_1$

where $$\varphi_0$$ is a solution of (43) 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:

$y = 0.17735\sqrt{x} - 0.075597x - 0.212836x^2 + 0.17363x^3 - 0.06254x^4$

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

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

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

$\begin{split}\begin{array}{ccccccc} -\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{array}\end{split}$

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

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

which can be programmed as:

$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)} .$
  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); 

## 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:

(44) $\begin{split}\begin{array}{rcl} |\nabla u - \nabla u_h|_{0,\Omega} &\le& C_1h \\ \|u - u_h\|_{0,\Omega} &\le& C_2h^2 \end{array}\end{split}$

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 (44). We will pick up numericall error if we use the numerical derivative, so we will use the following for (44).

$\begin{split}\begin{array}{rcl} \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{array}\end{split}$

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 << "("<

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:

$-\Delta u = sin(x+\pi/4.)*cos(y+\pi/4.)$

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. 141 The isovalue of solution $$u$$ with periodic boundary condition

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

Tip

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. 142 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

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

Tip

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 // 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) <

Tip

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 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); 

## 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{split}\begin{array}{rcll} -\Delta u &=& f & \textrm{ in }\Omega\\ u &=& g &\textrm{ on }\Gamma_D\\ \p u/\p n &=& 0 &\textrm{ on }\Gamma_N \end{array}\end{split}$

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:

$u=K_iu_S + u_R,\, u_R\in H^2(\textrm{near }\gamma_i),\, i=1,2$

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_e$$S 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{split}\begin{array}{rcl} \|u^0_h-u_e\|_{1,\Omega}/H1e&=&0.120421\\ \|u^a_h-u_e\|_{1,\Omega}/H1e&=&0.0150581 \end{array}\end{split}$

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{split}\begin{array}{rcll} -\Delta p &=& 1 & \textrm{ in }\Omega\\ p &=& g_d & \textrm{ on }\Gamma_D\\ \p p/\p n &=& g_n & \textrm{ on }\Gamma_N \end{array}\end{split}$

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{split}\begin{array}{rcll} \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{array}\end{split}$

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

The variational formulation is:

$\begin{split}\begin{array}{rcl} \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{array}\end{split}$

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. 145 and Fig. 146. Fig. 145 Density of the error indicator with isotropic $$P_{2}$$ metric Fig. 146 Density of the error indicator with isotropic $$P_{2}$$ metric

## 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); }