Finite element programming by FreeFem++
Atsushi Suzuki1
1 Cybermedia Center, Osaka University
atsushi.suzuki@cas.cmc.osaka-u.ac.jp

RIIT Tutorial, Kyushu Univ. 25, November 2016

1 / 70Numerical simulation with finite element method
I mathematical modeling
I discretization of time for evolution problem
I discretization scheme for the spacemesh generation / adaptive mesh refinement
stiffness matrix from finite elements and variational
formulation
I linear solver ⇐ CG, GMRES, direct solver: UMFPACK, MUMPS
I

I

FreeFem++ provides vast amounts of tools
I nonlinear solver
I optimization solver

2 / 70Outline I
Weak formulation for partial differential equations with 2ndorder derivatives
Poisson equation with mixed boundary conditions
Stokes equations with mixed boundary conditions
magneto-static problem with Dirichlet boundary conditions
P1/P2 finite elements and numerical quadrature
Finite element stiffness matrix from weak form
positivity of the matrix from coercivity of the bilinear forma penalty method to treat Dirichlet boundary data
direct method
CG / GMRES methods for symmetric/unsymmetric sparse
matrices
Nonlinear finite element solution by Newton method
stationary Navier-Stokes equations and a weak formulationdifferential calculus of nonlinear operator and Newton
iteration
3 / 70Outline II
Syntax and data structure of FreeFem++ language
control flow
data structure; array and FE object
sparse matrix
function macro
procedure to compute SpMV for built-in iterative solver
Schwarz algorithm as preconditioner for global Krylov iteration
overlapping subdomains and RAS/ASM
2-level algorithm with a coarse space
Time-dependent Navier-Stokes equations around a cylinder
boundary conditions of incompressible flow around acylinderThermal convection problem by Rayleigh-Bénard eqs.
governing equations by Boussinesq approximation
stationary solution by Newton method from numericalstationary solution
4 / 70Outline
Weak formulation for partial differential equations with 2ndorder derivatives
Poisson equation with mixed boundary conditions
Stokes equations with mixed boundary conditions
magneto-static problem with Dirichlet boundary conditions
P1/P2 finite elements and numerical quadrature
Finite element stiffness matrix from weak form
Nonlinear finite element solution by Newton method
Syntax and data structure of FreeFem++ language
Schwarz algorithm as preconditioner for global Krylov iteration
Time-dependent Navier-Stokes equations around a cylinder
5 / 70Poisson equation with mixed B.C. and a weak formulation: 1/2Ω ⊂ R2 , ∂Ω = ΓD ∪ ΓN
−4u = f in Ω,
u = g on ΓD ,∂u
= h on ΓN .
∂nweak formulation
V : function space, V (g) = {u ∈ V ; u = g on ΓD }.
V = C 1 (Ω) ∩ C 0 (Ω̄) ?
Find u ∈ V (g) s.t.
ZZ
−4u vdx =f vdx ∀v ∈ V (0)ΩΩ

Lemma (Gauss-Green’s
formula) nu, v ∈ V , n = 1 : outer normal to ∂ΩZ n2ZZ(∂i u)v dx = − u∂i v dx +ΩΩu ni v ds .

∂Ω

6 / 70Poisson equation with mixed B.C. and a weak formulation: 2/2ZZZ(∂1 u∂1 v+∂2 u∂2 v) dx − (∂1 u n1 +∂2 u n2 )v ds
Ω
Ω
∂Ω
Z
Z
=
∇u · ∇v dx −∇u · n v dsΩ
ΓD ∪ΓNZ
Z
v = 0 on ΓD ⇒ =∇u · ∇v dx −hv ds(−∂12 − ∂22 )u v dx =ΩΓNFind u ∈ V (g) s.t.
ZZZ∇u · ∇vdx =f vdx +ΩΩ

h vds ∀v ∈ V (0)ΓN

I a(·, ·) : V × V → R
I F (·) : V → R: bilinear form
: functional

Find u ∈ V (g) s.t.
a(u, v) = F (v) ∀v ∈ V (0)7 / 70FreeFem++ script to solve Poisson equation
finite element basis,
P span[ϕ1 , . . . , ϕN ] = Vh ⊂ V
uh ∈ Vh ⇒ uh = 1≤i≤N ui ϕiDirichlet data : u(Pj ) = g(Pj ) Pj ∈ ΓD
FindZuh ∈ Vh (g) s.t. Z
Z∇uh · ∇vh dx =Ωf vh dx +Ωh vh ds∀vh ∈ Vh (0).

ΓNexample1.edp
varf+matrix
mesh Th=square(20,20);
fespace Vh(Th,P1);
Vh uh,vh;func f=5.0/4.0*pi*pi*sin(pi*x)*sin(pi*y/2.0);
func g=sin(pi*x)*sin(pi*y/2.0);
func h=(-pi)/2.0*sin(pi*x);
solve poisson(uh,vh)=
int2d(Th)(dx(uh)*dx(vh)+dy(uh)*dy(vh))- int2d(Th)(f*vh) - int1d(Th,1)(h*vh)// boundary 1 : (x,0)+ on(2,3,4,uh=g);plot(uh);

8 / 70Stokes equations and a weak formulation : 1/3Ω = (0, 1) × (0, 1)−2∇ · D(u) + ∇p = f in Ω
∇ · u = 0 in Ω
u = g on ΓD ,2D(u)n − n p = h on ΓN .
strain rate tensor : [D(u)]i j = 21 (∂i uj + ∂j ui ).
I V (g) = {v ∈ H 1 (Ω)2 ; v = g on ΓD }, V = V (0)I Q = L2 (Ω)
weak formulation : Find (u, p) ∈ V (g) × Q s.t.
Z
ZZ
Z
2D(u) : D(v) dx − ∇ · v p dx = f · v dx + h · v dx ∀v ∈ V,
Ω
Ω
Ω
ΓN
Z
− ∇ · u q dx = 0 ∀q ∈ Q.
Ω

9 / 70Stokes equations and a weak formulation : 2/3

Lemma (Gauss-Green’s formula)u, v ∈ C 1 (Ω) ∩ C 0 (Ω̄), n : outer normal to ∂ΩZZZ(∂i u)v dx = − u∂i v dx +u ni v ds .
Ω

Ω

∂Ω

Z
−2

(∇ · D(u)) · v dx =
Ω

Z
−2
Ω

1
j ∂j 2 (∂i uj + ∂j ui )vi dx =

P P
i

Z

P

ZΩ

−
∂Ω

Z
Ω

P

i,j (∂i uj + ∂j ui )nj vi ds

Z
2D(u) : D(v) dx −=i,j (∂i uj + ∂j ui )∂j vi dx

2D(u) n · v ds∂Ω

from the symmetry of D(u)PPi,j (∂i uj + ∂j ui )∂j vi= i,j (∂i uj +∂j ui )(∂j vi +∂i vj )/2 = 2D(u):D(v).
10 / 70

Stokes equations and a weak formulation : 3/3
Z
Ω

Z

Z
P
p∂
v
dx
+
i i i
i p ni viΩ∂Ω
ZZ= − p∇ · v +pn · v

Pi (∂i p) vi dx = −

Ω

P

∂Ω

On the boundary ∂Ω = ΓD ∪ ΓN ,Z
(2D(u)n − n p) · v ds
ΓD ∪ΓNZ
Z(2D(u)n − n p) · v ds +(2D(u)n − n p) · v ds
=
ΓD
ΓN
Z
=
h · v ds. ⇐ V = {v ∈ H 1 (Ω)2 ; v = 0 on ∂Ω}.
ΓN

RemarkP
P−2[∇ · D(u)]i = − j ∂j (∂i uj + ∂j ui ) = − j ∂j2 ui = −[4u]i .
11 / 70

FreeFem++ script to solve Stokes equations by P2/P1
Find (u, p) ∈ Vh (g) × Qh s.t.
a(u, v) + b(v, p) + b(u, q) = (f, v) ∀(v, q) ∈ Vh × Qh .
example2.edp

fespace Vh(Th,P2),Qh(Th,P1);
func f1=5.0/8.0*pi*pi*sin(pi*x)*sin(pi*y/2.0)+2.0*x;
func f2=5.0/4.0*pi*pi*cos(pi*x)*cos(pi*y/2.0)+2.0*y;func g1=sin(pi*x)*sin(pi*y/2.0)/2.0;
func g2=cos(pi*x)*cos(pi*y/2.0);
func h1=3.0/4.0*pi*sin(pi*x)*cos(pi*y/2.0);
func h2=pi*cos(pi*x)*sin(pi*y/2.0)+x*x+y*y;
Vh u1,u2,v1,v2; Qh p,q;
macro d12(u1,u2) (dy(u1) + dx(u2))/2.0 //real epsln=1.0e-6;
solve stokes(u1,u2,p, v1,v2,q) =int2d(Th)(2.0*(dx(u1)*dx(v1)+dy(u2)*dy(v2)+
2.0*d12(u1,u2)*d12(v1,v2))-p*dx(v1)-p*dy(v2)-dx(u1)*q-dy(u2)*q)- int2d(Th)(f1*v1+f2*v1) - int1d(Th,1)(h1*v1+h2*v1)+ on(2,3,4,u1=g1,u2=g2);
plot([u1,u2],p,wait=1,value=true,coef=0.1);
12 / 70

stationary Navier-Stokes equations and a weak formulationΩ = (0, 1) × (0, 1)−2ν∇ · D(u) + u · ∇u + ∇p = f in Ω
∇ · u = 0 in Ω
u = g on ∂Ω
I V (g) = {v ∈ H 1 (Ω)2 ; v = g on ∂Ω}, V = V (0)R
I Q = L2 (Ω) = {p ∈ L2 (Ω) ;0Ω p dx = 0}

outflowbi/tri-linear
Z forms and weak formulation :a(u, v) = 2νD(u) : D(v) dx u, v ∈ H 1 (Ω)2Ω Z
a1 (u, v, w) = (u · ∇v) · w u, v, w ∈ H 1 (Ω)2Z Ωb(v, p) = − ∇ · v p dx v ∈ H 1 (Ω)2 , p ∈ L2 (Ω)Ω

Find (u, p) ∈ V (g) × Q s.t.
a(u, v) + a1 (u, u, v) + b(v, p) = (f, v) ∀v ∈ V,b(u, q) = 0∀q ∈ Q.
13 / 70magneto-static problem with Dirichlet b.c. : 1/2
constraints on the external force: ∇ · f = 0 in Ω ⊂ R3 .
∇ × (∇ × u) = f

in Ω,

∇·u=0

in Ω,

u×n=0

on ∂Ω.

H0 (curl ; Ω) = {u ∈ L2 (Ω)3 ; ∇ × u ∈ L2 (Ω)3 ; u × n = 0}find (u, p) ∈ H0 (curl ; Ω) × H01 (Ω)(∇ × u, ∇ × v) + (v, ∇p) = (f, v)(u, ∇q) = 0∀v ∈ H0 (curl ; Ω)∀v ∈ H01 (Ω)

has a unique solution.
I ∇ · f = 0 ⇒ p = 0.
I (∇×·, ∇×·) : coercive on W .
W = H0 (curl ; Ω) ∩ {u ∈ H(div; Ω) ; div u = 0}.
I H0 (curl ; Ω) = gradH 1 (Ω) ⊕ W .
0
14 / 70magneto-static problem with Dirichlet b.c. : 2/2
vector valued Sobolev space:H(curl ; Ω) = {u ∈ L2 (Ω)3 ; ∇ × u ∈ L2 (Ω)3 }.
with essential boundary conditions on ∂Ω:H0 (curl ; Ω) = {u ∈ L2 (Ω)3 ; ∇ × u ∈ L2 (Ω)3 ; u × n = 0}.
     
Z
Z
∂1
u1
v1
∂2  × u2  · v2  dx
∇ × u · vdx =
Ω
Ω
∂3
u3
v3


Z
(∂2 u3 − ∂3 u2 )v1 + (∂3 u1 − ∂1 u3 )v2 +  dx=Ω(∂1 u2 − ∂2 u1 )v3




Z
Z
(n2 u3 − n3 u2 )v1 +(u3 ∂2 v1 − u2 ∂3 v1 )+ (n3 u1 − n1 u3 )v2 +  ds
= −  (u1 ∂3 v2 − u3 ∂1 v2 )+  dx +∂ΩΩ
(u2 ∂1 v3 − u1 ∂2 v3 )(n1 u2 − n2 u1 )v3ZZ=u · ∇ × vdx +(n × u) · vdsΩ∂Ω

15 / 70stiffness matrix of magneto-static problem by FreeFem++load "msh3"
load "Dissection"
defaulttoDissection;
mesh3 Th=cube(20,20,20);
fespace VQh(Th, [Edge03d, P1]); // Nedelec element
VQh [u1, u2, u3, p], [v1, v2, v3, q];varf aa([u1, u2, u3, p], [v1, v2, v3, q]) =int3d(Th)((dy(u3)-dz(u2)) * (dy(v3)-dz(v2)) +(dz(u1)-dx(u3)) * (dz(v1)-dx(v3)) +(dx(u2)-dy(u1)) * (dx(v2)-dy(v1)) +dx(p)*v1+dy(p)*v2+dz(p)*v3 +dx(q)*u1+dy(q)*u2+dz(q)*u3);
matrix A = aa(VQh, VQh, solver=sparsesolver,
tolpivot=1.0e-2,strategy=102);

Nédélec element of degree 0 and P1
N0 (K) = (P0 (K))3 ⊕ [x × (P0 (K))3 ],P1 (K)

fespace VQh(Th,[Edge03d,P1]);
16 / 70

P1 finite element and sparse matrix
Th : triangulation of a domain Ω, triangular element K ∈ Th
piecewise linear element : ϕi |K (x1 , x2 ) = a0 + a1 x1 + a2 x2
ϕi |K (Pj ) = δi jZ
X Z
[A]i j = a(ϕj , ϕi ) =∇ϕj · ∇ϕi dx =∇ϕj · ∇ϕi dx.
Ω

K

K∈Th

1
15

18

15

18

12

16

19

19

11

11
16

12
6

17
6

13

7

14

8

1

14

8

1

9

9
10

2

17

13

7

3

10

2
3

4

5

4

5

A : sparse matrix, CRS (Compressed Row Storage) format tostore
17 / 70

P2 finite element
Th : triangulation of a domain Ω, triangular element K ∈ Th
piecewise quadratic element : 6 DOF on element K.
ϕi |K (x1 , x2 ) = a0 + a1 x1 + a2 x2 + a3 x21 + a4 x1 x2 + a5 x22
P1
ϕi |K (Pj ) = δi j
1
5

6

λ1

3
2

λ2

λ3

4

P3

P2

by using area coordinates {λ1 , λ2 , λ3 }, λ1 + λ2 + λ3 = 1.
   2  

ϕ1
1
−1 −1
λ1
λ1 (2λ1 − 1)
ϕ 2   1
 2  

−1
−1
  
  λ22  λ2 (2λ2 − 1)
ϕ 3    λ3  λ3 (2λ3 − 1)
1 −1 −1
 =

 

ϕ 4   λ2 λ3  =  4λ2 λ3 
4
  

 

ϕ 5   λ3 λ1   4λ3 λ1 
44
4λ1 λ2
ϕ6
λ1 λ2
fespace Vh(Th,P2);
18 / 70

numerical integration
Numerical quadrature:{Pi }i≤i≤m : integration points in K, {ωi }i≤i≤m : weights
m
X Z
X X
2
|(u − uh )(Pi )|2 ωi|u − uh |0,Ω =|u − uh |2 dx ∼K∈Th

KK∈Th i=1

formula : degree 5, 7 points, qf5pT,
P.C. Hammer, O.J. Marlowe, A.H. Stroud [1956]
area coordinates {λi }3i=1weight1 1
9
P1( 13 , √,)×1
3 3
40√|K|√√
15 6− 15 9+2 15
155− 15( 6−21, 21, 21√ )
|K| ×3 λ 3 P λ 2
1200
√
√
√
6+ 15 6+ 15 9−2 15
155+ 15( 21 , 21 , 21 )×31200 |K|
P2

λ1

P3

Remarkit is not good idea to use interpolation of continuous function tofinite element space, for verification of convergence order.
|Πh u − uh |1,Ω may be smaller (in extreme cases, super convergence)
19 / 70

treatment of Neumann data around mixed boundary
Neumann data is evaluated by line integral with FEM basis ϕi .
Z
h ϕi ds
ΓN

Dirichlet

NeumannQ

For given discrete
Neumann data, h is interpolated in FEM
Pspace, h = j hj ϕj |ΓN ,XjZ
ϕj ϕi ds.

hj
ΓN

On the node Q ∈ Γ̄D ∩ Γ̄N , both Dirichlet and Neumann arenecessary.
20 / 70

advantages of finite element formulation : 1/2
I weak formulation is obtained by integration by part withclear description on the boundary
I Dirichlet boundary condition, called as essential boundary

condition, is embedded in a functional space
I Neumann boundary condition, called as natural boundary

condition, is treated with surface/line integral by
Gauss-Green’s formula
I computation of local stiffness matrix is performed bynumerical quadrature formula

21 / 70

Outline
Weak formulation for partial differential equations with 2ndorder derivatives
Finite element stiffness matrix from weak form
positivity of the matrix from coercivity of the bilinear forma penalty method to treat Dirichlet boundary data
direct method
CG / GMRES methods for symmetric/unsymmetric sparse
matrices
Nonlinear finite element solution by Newton method
Syntax and data structure of FreeFem++ language
Schwarz algorithm as preconditioner for global Krylov iteration
Time-dependent Navier-Stokes equations around a cylinder 22 / 70

discretization and matrix formulation : 1/2
finite element basis,
P span[ϕ1 , . . . , ϕN ] = Vh ⊂ V
uh ∈ Vh ⇒ uh = 1≤i≤N ui ϕifinite element nodes {Pj }N
j=1 , ϕi (Pj ) = δi j Lagrange element
ΛD ⊂ Λ = {1, . . . , N } : index of node on the Dirichlet boundaryXVh (g) = {uh ∈ Vh ; uh =ui ϕi , uk = gk (k ∈ ΛD )}Find uh ∈ Vh (g) s.t.
a(uh , vh ) = F (vh ) ∀vh ∈ Vh (0).
Find {uj }, uk = gk (k ∈ ΛD ) s.t.
XXXa(u j ϕj ,vi ϕi ) = F (vi ϕi ) ∀{vi }, vk = 0(k ∈ ΛD )jii

Find {uj }j∈Λ s.t.
Xa(ϕj , ϕi )uj = F (ϕi )juk = gk∀i ∈ Λ \ ΛD∀k ∈ ΛD
23 / 70

discretization and matrix formulation : 2/2
FindX{uj }j∈Λ\ΛD s.t.
Xa(ϕj , ϕi )uj = F (ϕi ) −a(ϕk , ϕi )gk
j∈Λ\ΛDA = {a(ϕj , ϕi )}i,j∈Λ\ΛD∀i ∈ Λ \ ΛDk∈ΛDA ∈ Rn×n , f ∈ Rn , n = #(Λ \ ΛD )Lemma
A : positive definite (coercive) ⇔ (Au, u) > 0 ∀u 6= 0
⇒ Au = f has a unique solution.
A : bijective
I injective: Au = 0, 0 = (Au, u) > 0⇒ u = 0.

I surjective:
Rn = ImA ⊕ (ImA)⊥ , u ∈ (ImA)⊥ ⇒ (Av, u) = 0 ∀v ∈ Rn

by putting v = u, 0 = (Au, u) ⇒ u = 0(ImA)⊥ = {0} ⇒ ImA = Rn .
unymmetric ⇒ solution by LDU -factorization, GMRES method
symmetric ⇒ solution by LDLT -factorization, CG method
24 / 70

Sobolev space
P1 element element space does not belong to C 1 (Ω).
V = H 1 (Ω),
Z(u, v) =u v + ∇u · ∇v, ||u||21 = (u, u) < +∞ΩZ
2||u||0 =u u,ZΩ|u|21 =
∇u · ∇u.
Ω

H01 = {u ∈ H 1 (Ω) ; u = 0 on ∂Ω}.

Lemma (Poincaré’s inequality)∃C(Ω) u ∈ H01 ⇒ ||u||0 ≤ C(Ω)|u|1 .
RΩ ∇u · ∇v is coercive on V .
1a(u, u) = |u|21 ≥ 2 ||u||20 .

a(u, v) =C

25 / 70

FreeFem++ script to solve Poisson eq. using matrix
Find uh ∈ Vh (g) s.t. a(uh , vh ) = F (vh ) ∀vh ∈ Vh (0).
example3.edp

solve

Vh u,v;varf aa(u,v)=int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))+on(2,3,4,u=g);
varf external(u,v)=int2d(Th)(f*v)+int1d(Th,1)(h*v)+on(2,3,4,u=g);
real tgv=1.0e+30;
matrix A = aa(Vh,Vh,tgv=tgv,solver=CG);
real[int] ff = external(0,Vh,tgv=tgv);
u[] = A^-1 * ff;// u : fem unknown, u[] : vector
plot(u);useful liner solver; solver=CG
iterative solver for SPD matrix
GMRES
iterative solver for nonsingular matrix
UMFPACK
direct solver for nonsingular matrix
sparsesolver other solvers called by dynamic link
26 / 70

penalty method to solve inhomogeneous Dirichlet problem
modification of diagonal entries of A where index k ∈ ΛD
penalization parameter τ = 1/ε; tgvukττ gk , k ∈ Λ D
=

ui

fi

[A]i j = a(ϕj , ϕi )
τ uk +

X

ak j uj = τ gk ⇔ uk − gk = ε(−

j6=k

X

X

ak j uj ),j6=kai j uj = fi∀i ∈ {1, . . . , N } \ ΛD .

jkeeping symmetry of the matrix without changing index
numbering.
27 / 70

LDLT /LDU factorization

 


1
0
a11 β1T
a11 β1T=α1 A220I2α1 a−1
S22
11




1
0
a11 0 1 a−1
β1T
11=
0 I2 0α1 a−1
S22
I2
11T
Schur complement S22 = A22 − α1 a−1
11 β1

rank-1 update

LDU -factorization with symmetric pivoting : A = ΠT LDU Π
do k = 1, · · · , N
find k < l ≤ n max|A(l, l)|,exchange rows and columns : A(k, ∗) ↔ A(l, ∗), A(∗, k) ↔ A(∗, l).
rank-1 updatekl

28 / 70

conjugate gradient method
A~u = f~.
preconditioner Q ∼ A−1
Krylov subsp. : Kn (Q~r 0 , QA)=span[Q~r 0 , QAQ~r ,0. . . , (QA)n−1 Q~r 0 ]Find ~u n ∈ Kn (Q~r 0 , QA) + ~u 0 s.t.
(A~u n − f~, ~v ) = 0 ∀~v ∈ Kn (Q~r 0 , QA).
Preconditioned CG method~u 0 : initial step for CG.
~r 0 = f~ − A~u 0
p~ 0 = Q~r 0 .
loop n = 0, 1, . . .
αn = (Q~r n , ~r n )/(A~p n , p~ n ),
n+1
n
n
~u
= ~u + αn p~ ,n+1~r
= ~r n − αn A~
p n,n+1
if ||~r|| <  exit loop.
βn = (Q~r n+1 , ~r n+1 )/(Q~r n , ~r n ),p~ n+1 = Q~r n+1 + βn p~ n .
LinearCG(opA,u,f,precon=opQ,nbiter=100,eps=1.0e-10)
29 / 70

GMRES method : 1/2
Krylov subspace : Kn (~r 0 , A) = span[~r 0 , A~r 0 , . . . , An−1~r 0 ]Find ~u n ∈ Kn (~r 0 , A) + ~u 0 s.t.
||A~u n − f~|| ≤ ||A~v n − f~|| ∀~v ∈ Kn (~r 0 , A) + ~u 0 .
Vm : Arnoldi basis generated by Gram-Schmidt orthogonization
for Krylov vectors.
~u = Vm ~y , ~y ∈ Rm

J(~y ) := ||AVm ~y − ~r0 ||T= ||Vm+1(AVm ~y − ~r0 )||TT= ||(Vm+1
AVm )~y − (Vm+1~r0 )||= ||H̄m ~y − β~e1 || .
Find ~y ∈ Rm

(β = ||~r0 ||)J(~y ) ≤ J(~z) ∀~z ∈ Rm .

minimization problem with Hessenberg matrix H̄m ∈ R(m+1)×m
is solved by Givens rotation.
LinearGMRES(opA,u,f,precon=opQ,nbiter=100,
eps=1.0e-10)
30 / 70

GMRES method : 2/2
Arnoldi method (Gram-Schmidt method on Krylov subspace)||~v1 || = 1; do j = 1, 2, . . . , m
do i = 1, 2, . . . , jPw~ j := A~vj − 1≤i≤j hi j ~vi ,~vj+1 := w
~ j /hj+1 j ,hi j := (A ~vj , ~vi )hj+1 j := ||w~ j ||

Givens rotation matrices Ωi ∈ R(m+1)×(m+1)
Ii−1 csii, c1 := √ h1 1
Ωi := , s1 := √ 2h2 1 2 .
h21 1 +h22 2
h1 1 +h2 2−si ci I
m−i
Qm := Ωm Ωm−1 · · · Ω1 ∈ R(m+1)×(m+1) ,R̄m := Qm H̄m : upper triangular,Tḡm := Qm (βe1 ) = [ γ1 γ2 · · · γm+1 ] , 

Rm
gm
m×m
R̄m :=(Rm ∈ R), ḡm :=(gm ∈ Rm ).
0 ··· 0γm+1min||βe1 − H̄m y|| = min||ḡm − R̄m y|| = |γm+1 | = |s1 s2 · · · sm |β .
−1 g attains the minimum.
ym = Rmm−1 (1 ≤ m ≤ M ) for all non-singular matrix A .
Remark : ∃Rm
31 / 70

advantages of finite element formulation 2/2
I solvability of linear system is inherited from solvability ofcontinuous weak formulation
better to learn for efficient computation
I treatment of Dirichlet boundary conditions in FreeFem++with explicit usage of matrix and linear solver
I Direct solver like UMFPACK is efficient for 2D problems,

but for 3D iterative solvers such as CG/GMERES withgood preconditioner (e.g. additive Schwarz method) arenecessary.

32 / 70

Outline
Weak formulation for partial differential equations with 2ndorder derivatives
Finite element stiffness matrix from weak form
Nonlinear finite element solution by Newton method
stationary Navier-Stokes equations and a weak formulationdifferential calculus of nonlinear operator and Newton
iteration
Syntax and data structure of FreeFem++ language
Schwarz algorithm as preconditioner for global Krylov iteration
Time-dependent Navier-Stokes equations around a cylinder
Thermal convection problem by Rayleigh-Bénard eqs.

33 / 70

stationary Navier-Stokes equations and a weak formulationΩ = (0, 1) × (0, 1)−2ν∇ · D(u) + u · ∇u + ∇p = f in Ω
∇ · u = 0 in Ω
u = g on ∂Ω
I V (g) = {v ∈ H 1 (Ω)2 ; v = g on ∂Ω}, V = V (0)R
I Q = L2 (Ω) = {p ∈ L2 (Ω) ;0Ω p dx = 0}

outflowbi/tri-linear
Z forms and weak formulation :a(u, v) = 2νD(u) : D(v) dx u, v ∈ H 1 (Ω)2Ω Z
a1 (u, v, w) = (u · ∇v) · w u, v, w ∈ H 1 (Ω)2Z Ωb(v, p) = − ∇ · v p dx v ∈ H 1 (Ω)2 , p ∈ L2 (Ω)Ω

Find (u, p) ∈ V (g) × Q s.t.
a(u, v) + a1 (u, u, v) + b(v, p) = (f, v) ∀v ∈ V,b(u, q) = 0∀q ∈ Q.
34 / 70

nonlinear system of the stationary solution
A(u, p ; v, q) = a(u, v) + a1 (u, u, v) + b(v, p) + b(u, q)nonlinear problem:
Find (u, p) ∈ V (g) × Q s.t. A(u, p ; v, q) = (f, v) ∀(v, q) ∈ V × Q.
a1 (·, ·, ·) : trilinear form,
a1 (u + δu,u + δu, v) = a1 (u, u + δu, v) + a1 (δu, u + δu, v)= a1 (u, u, v) + a1 (u, δu, v) + a1 (δu, u, v) + a1 (δu, δu, v)
A(u + δu, p + δp ; v, q) − A(u, p ; v, q)=a(u + δu, v) − a(u, v)+ b(v, p + δp) − b(v, p) + b(u + δu, q) − b(u, q)+ a1 (u + δu, u + δu, v) − a1 (u, u, v)=a(δu, v) + b(v, δp) + b(δu, q) + a1 (δu, u, v) + a1 (u, δu, v) + O(||δu||2 )Find (δu, δp) ∈ V × Q s.t.
a(δu, v) + b(v, δp) + b(δu, q)+a1 (δu, u, v) + a1 (u, δu, v) =− A(u, p ; v, q) ∀(v, q) ∈ V × Q
35 / 70

Newton iteration(u0 , p0 ) ∈ V (g) × Q
loop n = 0, 1 . . .
Find (δu, δp) ∈ V × Q s.t.
a(δu, v) + b(v, δp) + b(δu, q) + a1 (δu, un , v) + a1 (un , δu, v) =A(un , pn ; v, q) ∀(v, q) ∈ V × Q
if ||(δu, δp)||V ×Q ≤ ε then breakun+1 = un − δu,pn+1 = pn − δp.
loop end.
example4.edp(u(0) , p(0) ) ∈ V (g) × Q : solution of the Stokes eqs., ν = 1.
while (ν > νmin )Newton iteration (u(k+1) , p(k+1) ) ∈ V (g) × Q from (u(k) , p(k) ).
ν = ν/2, k++.
while end.
initial guess from the stationary state of lower Reynolds number
36 / 70mesh adaptation
fespace XXMh(Th,[P2,P2,P1]);
XXMh [u1,u2,p];real lerr=0.01;
Th=adaptmesh(Th, [u1,u2], p, err=lerr, nbvx=100000);[u1,u2,p]=[u1,u2,p]; // interpolation on the new mesh⇒err : P1 interpolation error levelnbvx : maximum number of vertexes to be generated.
37 / 70

FreeFem++ for non-linear problem
I Jacobian of Newton method needs to be obtained fromdifferential calculus
I good initial guess is necessary to achieve fast

convergence of Newton method
I combination with mesh adaptation is an effective techniqueto improve initial conditions for Newton method

38 / 70

Outline
Weak formulation for partial differential equations with 2ndorder derivatives
Finite element stiffness matrix from weak form
Nonlinear finite element solution by Newton method
Syntax and data structure of FreeFem++ language
control flow
data structure; array and FE object
sparse matrix
function macro
procedure to compute SpMV for built-in iterative solver
Schwarz algorithm as preconditioner for global Krylov iteration
Time-dependent Navier-Stokes equations around a cylinder 39 / 70

syntax of FreeFem++ script
for loopfor (int i=0; i<10; i++) {...
if (err < 1.0e-6) break;}while loop
int i = 0;
while (i < 10) {...
if (err < 1.0e-6) break;i++;}procedure (function)
func real[int] ff(real[int] &pp) { // C++ reference
real [int] qq;
qq = ... pp ;// calculate qq from pp
return qq;// return data real[int]}40 / 70array and finite element object

mesh Th=square(20,20);
fespace Xh(Th, P1); // P1 finite element on Th
Xh uh;// P1 finite element object
int n = Xh.ndof;// DOF of FE space Xhreal[int] u(n);//u = uh[];// copy data from FE object to arrayfor (int i = 0; i < n; i++)// index-wise copy, but slowu(i)=uh[](i);real[int] u;// not yet allocatedu.resize(10);// same as C++ STL vector
real[int] vv = v; // allocated as same size of v.na(2)=0.0 ;// set value of 3rd index
a += b;// a(i) = a(i) + b(i)a = b .* c ;// a(i) = b(i) * c(i); element-wisea = b < c ? b : c // a(i) = min(b(i), c(i)); C-syntax
a.sum;// sum a(i);a.n;// size of array

There are other operations such as `1 , `2 , `∞ -norms, max, min.
cf. Finite element analysis by mathematical programming language
FreeFem++, Ohtsuka-Takaishi [2014].
41 / 70

vector-valued finite element object and 2D array
vector valued FE object
mesh Th=square(20,20);
fespace Xh(Th, [P1, P1]); // P1-P1 finite element on Th
Xh [uh1, uh2];// P1 finite element object

real[int] u(uh1[].n+uh2[].n);// allocation of array
real[int] v(Xh,nodf);//
u = [uh1[], uh2[]]; // copy of vector-valued object
v = uh1[];// u1h[] shows a pointer to an array// in [uh1, uh2]

2-D arraymesh Th=square(20,20);
fespace Xh(Th, P1); // P1 finite element on Th
Xh[int] uh(10);// 10 sized array of FE object
real[int,int] uu(Xh.ndof, 10); // 2D arrayfor (int m = 0; m < 10; m++) {uu(:,m)=uh[m][](:);// copy each array}
42 / 70

variational form and sparse matrix
mesh Th=square(20,20);
fespace Xh(Th, P1);
Xh u, v;func fnc = sin(pi * x + pi * y);varf a(u, v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))+ on(1,u=0.0);
varf f(u, v) = int2d(Th)(fnc * v);
matrix AA = a(Xh, Xh, solver=UMFPACK);
real[int] x(Xh.ndof), b(X.ndof);
b = f(Xh, 0);// RHS from external force
x = AA^-1 * b;// solution of linear system
int n = 10;// example to generate
real[int, int] mm(n, n); // sparse matrix from 2D array
Xh [int] uh(n);for (int i = 0; i < n; i++)uh[i] = sin(pi * real(i+1) * x);
for (int j = 0; j < n; j++)for (int i = 0; i < n; i++)mm(i, j)=uh[i][]’*uh[j][]; //’ inner productmatrix aa = mm;// copy into sparse matrix
set (aa, solver=UMFPACK);// set sparse solver
43 / 70

function macro and FEM object
func fnc = sin(pi*x)*sin(pi*y); // function with x,ymesh Th=square(20,20);
fespace Vh(Th, P2);Vh fh;
fh = fnc;// P2 interpolation using mesh data
real sum0 = int2d(Th)(fh * fh); //numerical quadrature// applied to P2 FE object
real sum1 = int2d(Th)(fnc * fnc);//numerical quadrature// applied to (x,y)-func.

quadrature qf5pT

44 / 70

procedure and CG built-in function : 1/2

int n = 20;
example5.edp
varf+matrix
mesh Th=square(n,n);
fespace Vh(Th,P1); Vh u,v;
func f = 5.0/4.0*pi*pi*sin(pi*x)*sin(pi*y/2.0);
func h = (-pi)/2.0*sin(pi*x);
func g = sin(pi*x)*sin(pi*y/2.0);
varf aa(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))+ on(2,3,4,u=1.0);
varf external(u,v) = int2d(Th)(f*v) + int1d(Th,1)(h*v);
real tgv=1.0e+30;
matrix A;
real[int] bc = aa(0, Vh, tgv=tgv);
func real[int] opA(real[int] &pp){// SpMV operation only for inner
pp = bc ? 0.0 : pp;
real[int] qq = A * pp;// node without Dirichlet bdry.
pp = bc ? 0.0 : qq;
return pp;}
45 / 70

procedure and CG built-in function : 2/2
func real[int] opQ(real[int] &pp) // diagonal scaling{for (int i = 0; i < pp.n; i++) {pp(i) = pp(i) / A(i, i);}pp = bc ? 0.0 : pp;
return pp;}
A = aa(Vh, Vh, tgv=tgv, solver=sparsesolver);real[int] ff = external(0, Vh);// g is valid only on the bdry 1
v = g;
u[] = bc ? v[] : 0.0; // u[] has Dirichlet data w/o tgv
v[] = A * u[];// v[] = A_{12}*gv[] = bc ? 0.0 : v[];// ff_{1} -= A_{12}*u_{2}ff -= v[];
ff = bc ? 0.0 : ff;
LinearCG(opA, u[], ff, precon=opQ, nbiter=200,
eps=1.0e-10,verbosity=50);
46 / 70

Outline
Weak formulation for partial differential equations with 2ndorder derivatives
Finite element stiffness matrix from weak form
Nonlinear finite element solution by Newton method
Syntax and data structure of FreeFem++ language
Schwarz algorithm as preconditioner for global Krylov iteration
overlapping subdomains and RAS/ASM
2-level algorithm with a coarse space
Time-dependent Navier-Stokes equations around a cylinder
Thermal convection problem by Rayleigh-Bénard eqs.
47 / 70

Schwarz algorithm : 1/2Ω2Ω1
Jacobi-Schwarz algorithm
u01 , u02 : given
loop n = 0, 1, 2, · · ·−4un+1= f in Ω1
1−4un+1= f in Ω22

un+1= 0 on ∂Ω1 ∩ ∂Ω
1

un+1
= 0 on ∂Ω2 ∩ ∂Ω
2

un+1
= un2 on ∂Ω1 ∩ Ω̄2
1

un+1= un1 on ∂Ω2 ∩ Ω̄1
2

I parallel computation, extendable to more than two

subdomains
I overlapping subdomain

48 / 70

Restricted Additive Schwarz algorithm : 1/2
partition
P of unityu = 2i=1 Ei (χi ui )I Ei : extension from Ωi to ΩI χi : partition of unity function in Ω̄iloop n = 0, 1, 2, · · ·−4win+1 = f in Ωiwin+1 = 0 on ∂Ωi ∩ ∂Ω
win+1 = un on ∂Ωi ∩ Ω̄j
P2

n+1
)
i=1 Ei (χi wi
by substituting 4un that satisfies un = 0 on ∂Ω

un+1 =i ∩ ∂Ω
n+1
nn−4wi + 4u = f + 4u in Ωiwin+1 − un = 0 on ∂Ωi ∩ ∂Ω,
win+1 − un = 0 on ∂Ωi ∩ Ω̄jPn+1) − 2i=1 Ei (χi uni )i=1 Ei (χi wi
P= 2i=1 Ei (χi (win+1 − un ))un+1 − un =P2

49 / 70

Restricted Additive Schwarz algorithm : 2/2
Restricted additive Schwarz (RAS) algorithm
u01 , u02 : given
loop n = 0, 1, 2, · · ·rn = f + 4un−4vin+1 = rn in Ωivin+1 = 0 on ∂Ωi
u

n+1

n

=u +

2
X

Ei (χi vin+1 )i=1RAS ⇒ Jacobi-Schwarz method
proof by induction−4(un + vin ) = −4(un ) + rn = f in Ωi
un + vin = un on ∂Ωi ∩ ∂Ω .
un = E1 (χ1 un1 ) + E2 (χ2 un2 )
un = E1 (0 · un1 ) + E2 (1 · un2 ) = un2 on ∂Ωi ∩ Ω̄2 .
50 / 70

Additive Schwarz algorithm
loop n = 0, 1, 2, · · ·−4win+1 = f in Ωiwin+1 = 0 on ∂Ωi ∩ ∂Ω
win+1 = un on ∂Ωi ∩ Ω̄jun+1 =

P2

n+1)
i=1 Ei (wiAdditive Schwarz Method (ASM)
u01 , u02 : given
loop n = 0, 1, 2, · · ·rn = f + 4un−4vin+1 = rn in Ωivin+1 = 0 on ∂Ωiun+1 = un +2X

Ei (vin+1 )i=151 / 70

Schwarz methods as preconditioner
ASM preconditioner
−1
=
MASM

M
X

RpT (Rp ARpT )−1 Rp

p=1−1
ASM does not converge as fixed point iteration, but MASM
issymmetric and works well as a preconditioner for CG method.

RAS preconditioner
−1
MRAS
=

M
X

RpT Dp (Rp ARpT )−1 Rp

p=1−1
is not symmetric and then works
RAS does converge but MRAS
as a preconditioner for GMRES method.
convergence of ASM, RAS: slow for many subdomains⇒ coarse space

52 / 702-level Schwarz methods with a coarse spacecoarse space by Nicolaides
Dp : discrete representation
the partition of unity
PM of
TD R = I ,RN
p=1 p p p{~zp } ⊂ RN : basis of coarse space, Z = [~z1 , · · · , ~zM ].
~zp = RpT Dp Rp~1,
R0 = Z T .
2-level ASM preconditioner−1
MASM,2= R0T (R0 AR0T )−1 R0 +M
X

RpT (Rp ARpT )−1 Rp

p=1

2-level RAS preconditioner−1
MRAS,2= R0T (R0 AR0T )−1 R0 +MX

RpT Dp (Rp ARpT )−1 Rp

p=1

53 / 70

FreeFem++ script for a 3D domain computation
Subdomain problems : solved by a direct solver
can be distributed in paralle machine.

example6.edp

I domain decomposition by METIS

metisdual()
I overlapping subdomain from non-overlapping one using

numerical diffusion operator
func bool AddLayers()
I partition of unity

func bool SubdomainsPartitionUnity()
cf. V. Dolean, P Jolivet, F. Nataf, An Introduction to Domain
Decomposition Methods – Algorithms, Theory, and Parallel
Implementation, SIAM, 2015. ISBN 978-1-611974-05-8

54 / 70

Outline
Weak formulation for partial differential equations with 2ndorder derivatives
Finite element stiffness matrix from weak form
Nonlinear finite element solution by Newton method
Syntax and data structure of FreeFem++ language
Schwarz algorithm as preconditioner for global Krylov iteration
Time-dependent Navier-Stokes equations around a cylinder
boundary conditions of incompressible flow around acylinderThermal convection problem by Rayleigh-Bénard eqs.
55 / 70

incompressible flow around a cylinder : boundary conditionsΩ = (−1, 9) × (−1, 1)

Γ4

Γ3

ω

Γ2

Γ1
∂u
+ u · ∇u − 2ν∇ · D(u) + ∇p = 0 in Ω∂t
∇ · u = 0 in Ω
u = g on ∂Ωboundary conditions:
Poiseuille flow on Γ4 : u = (1 − y 2 , 0).u·n=0(2νD(u)n − np) · t = 0
no-slip boundary condition on ω : u = 0
outflow boundary condition on Γ2 : 2νD(u)n − np = 0
slip boundary condition on Γ1 ∪ Γ3 :56 / 70slip boundary conditions and function space
Γ3
Γ4

ω

Γ2


Γ1
slip boundary condition on Γ1 ∪ Γ3 :

Z

u·n=0
(2ν D(u)n − n p) · t = 0

I V (g) = {v ∈ H 1 (Ω)2 ; v = g on Γ4 ∪ ω, v · n = 0 on Γ1 ∪ Γ3 },I Q = L2 (Ω).
non-slipZ(2νD(u)n − n p) · vds =Γ1 ∪Γ3(2νD(u)n − n p) · (vn n + vt t)ds
Γ1 ∪Γ3

Z(2νD(u)n − n p) · (v · n)n ds=
Γ1 ∪Γ3

Z(2νD(u)n − n p) · tvt ds = 0+
Γ1 ∪Γ357 / 70characteristic line and material derivativeu(x1 , x2 , t) : Ω × (0, T ] → R2 , given velocity field.
φ(x1 , x2 , t) : Ω × (0, T ] → R.
X(t) : (0, T ] → R2 , characteristic curve :dX(t) = u(X(t), t), X(0) = X0dtd∂dφ(X(t), t) =∇φ(X(t), t) · X(t) + φ(X(t), t)
dt
dt
∂t
∂
=∇φ(X(t), t) · u(X(t), t) + φ(X(t), t)∂t
Dφ∂= φ + u · ∇φ.
Dt∂tapproximation by differenceDφ(X(t), t)φ(X(t), t) − φ(X(t − ∆t), t − ∆t)∼Dt∆tmaterial derivative :58 / 70characteristic Galerkin method to discretized material derivative
approximation by Euler method :tn < tn+1 , tn+1 = ∆t + tn .
X(tn+1 ) = x

L2
L3
K
n

X1,h(K)

L4

X(tn ) = X n (x) + O(∆t2 )

L1

uh

L
X n (x) = x − u(x, tn )L
φ(x, tn+1 ) − φ(X n , tn )Dφ(X(tn+1 ), tn+1 )=+ O(∆t)Dt∆tφn+1 − φn ◦ X n∼.
∆t
un : obtained in the previous time step.
Find (un+1 , pn+1 ) ∈ V (g) × Q s.t.

 n+1
u
− un ◦ X n, v + a(un+1 , v) + b(v, pn+1 ) = 0 ∀v ∈ V,∆t5

6b(un+1 , q) = 0∀q ∈ Q.
59 / 70

FreeFem++ script using characteristic Galerkin method
FreeFem++ provides convect to compute (un ◦ X n , ·).
example7.edp
real nu=1.0/Re;
real alpha=1.0/dt;
int i;
problem NS([u1,u2,p],[v1,v2,q],solver=UMFPACK,init=i) =int2d(Th)(alpha*(u1*v1 + u2*v2)+2.0*nu*(dx(u1)*dx(v1)+2.0*d12(u1,u2)*d12(v1,v2)+dy(u2)*dy(v2))- p * div(v1, v2) - q * div(u1, u2))
- int2d(Th)(alpha*( convect([up1,up2],-dt,up1)*v1+convect([up1,up2],-dt,up2)*v2) )+ on(1,3,u2=0)+on(4,u1=1.0-y*y,u2=0)+on(5,u1=0,u2=0);

for (i = 0; i <= timestepmax; i++) {up1 = u1; up2 = u2; pp = p;NS;// factorization is called when i=0
plot([up1,up2],pp,wait=0,value=true,coef=0.1);}
60 / 70

FreeFem++ script for mesh generation around a cylinder
Delaunay triangulation from nodes given on the boundary
boundary segments are oriented and should be connected.
int n1 = 30;
int n2 = 60;
border ba(t=0,1.0){x=t*10.0-1.0;y=-1.0;label=1;};
border bb(t=0,1.0){x=9.0;y=2.0*t-1.0;label=2;};
border bc(t=0,1.0){x=9.0-10.0*t;y=1.0;label=3;};
border bd(t=0,1.0){x=-1.0;y=1.0-2.0*t;label=4;};
border cc(t=0,2*pi){x=cos(t)*0.25+0.75;
y=sin(t)*0.25;label=5;};
mesh Th=buildmesh(ba(n2)+bb(n1)+bc(n2)+bd(n1)+cc(-n1));plot(Th);

61 / 70

stream line for visualization of flow around a cylinder : 1/2
∂2 ψ
stream function ψ : Ω → R, u =.
−∂1 ψboundary conditions for the stream line:
Z y
Z y
y3
u1 (x1 , t)dt =∂2 ψ(x1 , t)dt = ψ(x1 , y) − ψ(x1 , 0)
inlet: y −
=
3
0
0
y3
y3
ψ(x1 , y) = ψ(x1 , 0) + y −
=y− .
3
3
Z
Z
x

x

u2 (t, ± 1)dt =slip: 0 =x1−∂1 ψ(t, ±1)dt = ψ(x1 , ±1) − ψ(x, ±1)
x1

ψ(x, ±1) = ψ(x1 , ±1) = ψ(x1 , 0) ± 2/3 = ±2/3.
Z θZ θ−∂1 ψ r sin θ + ∂2 ψ r cos θcylinder: 0 =u · n dθ =
−π

−π

Z θ
=

∂
ψ(r, θ)dθ.
∂θ−πcenter from inlet: u2 = 0 ⇒ same as slip wall ,ψ|ω = ψ(x1 , 0) = 0.

62 / 70

stream line for visualization of flow around a cylinder : 2/2
Γ3
Γ4

ω

Γ2

Γ1
slip boundary condition on Γ1 ∪ Γ3 , outflow on Γ2 .

63 / 70

Outline
Weak formulation for partial differential equations with 2ndorder derivatives
Finite element stiffness matrix from weak form
Nonlinear finite element solution by Newton method
Syntax and data structure of FreeFem++ language
Schwarz algorithm as preconditioner for global Krylov iteration
Time-dependent Navier-Stokes equations around a cylinder
Thermal convection problem by Rayleigh-Bénard eqs.
governing equations by Boussinesq approximation
stationary solution by Newton method from numericalstationary solution

64 / 70thermal convection in a box : 1/2
Γ3 : θ = θ 0 , u 2 = 0
Γ4 : ∂n θ = 0, u1 = 0

Γ2 : ∂n θ = 0, u1 = 0

Γ1 : θ = θ0 + ∆θ, u2 = 0
Rayleigh-Bénard equations∂u+ u · ∇u − 2∇ · µ0 D(u) + ∇p = −ρ g~e2 in Ω ,ρ∂t∇ · u = 0 in Ω ,∂θ+ u · ∇θ − ∇ · (κθ) = 0 in Ω .
∂t~e2 : unit normal of y-direction
d : height of the box, g : gravity acceleration,
κ : thermal diffusivity, µ0 : viscosity
65 / 70thermal convection in a box : 2/2
Boussinesq approximation : ρ = ρ0 {1 − α(θ − θ0 )}, θ0 = 0.
ρ0 : representative density, α : thermal expansion coefficient.
non-dimensional Rayleigh-Bénard equations
1∂u+ u · ∇u − 2∇ · D(u) + ∇p = Raθ~e2 in Ω ,P r ∂t∇ · u = 0 in Ω ,∂θ+ u · ∇θ − 4θ = 0 in Ω∂tu · n = 0 on ∂Ω ,
θ = 1 on Γ1 ,
θ = 0 on Γ3 ,∂n θ = 0 on Γ2 ∪ Γ4 .
µ0
: Prandtl number,κρ0
ρ0 gα∆θd3
I Ra =: Rayleigh number.
κµ0
I Pr =

66 / 70

a weak form to solve time-dependent Rayleigh-Bénard eqs.
I velocity : V = {v ∈ H 1 (Ω)2 ; v · n = 0 on ∂Ω},
R
I pressure : Q = L2 (Ω) = {p ∈ L2 (Ω) ;0Ω p dx = 0},
I temperature : ΨD = {θ ∈ H 1 (Ω) ; θ = 1 on Γ1 , θ = 0 on Γ3 }.

bilinear forms:Ra0 (u, v) = RΩ 2D(u) : D(v),c0 (θ, ψ) = Ω ∇θ · ∇ψ.

Rb(v, p) = − Ω ∇ · v p,

example8.edp
using Characteristic Galerkin method:nn(u , θ ) ∈ V × ΨD : from previous time step
Find (un+1 , pn+1 , θn+1 ) ∈ V × Q × ΨD s.t.
 n+1

1
u
− un ◦ X n, v + a0 (un+1 , v) + b(v, pn+1 ) = Ra(θn~e2 , v)
Pr∆t∀v ∈ V,b(un+1 , q) = 0∀q ∈ Q, n+1θ− θn ◦ X n, ψ + c0 (θn+1 , ψ) = 0∆t∀ψ ∈ Ψ0 .
67 / 70

a weak form to solve stationary Rayleigh-Bénard eqs.
trilinear forms and bilinear form for the Navier-Stokes eqs.
R
I a1 (u, v, w) = 1
P r Ω (u · ∇v) · w
RI c1 (u, θ, ψ) =Ω (u · ∇θ) · ψ
I A(u, p ; v, q) = a0 (u, v) + a1 (u, u, v) + b(v, p) + b(u, q)
example9.edp
Newton iteration (u0 , p0 , θ0 ) ∈ V × Q × ΨDloop n = 0, 1 . . .
Find (δu, δp, δθ) ∈ V × Q × Ψ0 s.t.
a0 (δu, v) + b(v, δp) + b(δu, q) + a1 (δu,un , v) + a1 (un , δu, v)−Ra(δθ~e2 , v) = A(un , pn ; v, q) − Ra(θn~e2 , v) ∀(v, q) ∈ V × Q
c0 (δθ, ψ) + c1 (un , δθ, ψ) + c1 (δu, θn , ψ) = c0 (θn , ψ) + c1 (un , θn , ψ)∀ψ ∈ Ψ0
if ||(δu, δp, δθ)||V ×Q×Ψ ≤ ε then breakun+1 = un − δu, pn+1 = pn − δp, θn+1 = θn − δθ.
loop end.
initial data ⇐ stationary solution by time-dependent problem.
68 / 70

Application of finite element method to fluid problems
Time-dependent Navier-Stokes equations
I material derivative is approximated by Characteristic

Galerkin method
I functional space of pressure depends on boundary

conditions of flow, e.g., inflow, non-slip, slip, and outflow.
Thermal convection problem by Rayleigh-Bénard equations
I time-dependent problems for fluid and temperature by

convection are solved by Characteristic Galerkin method
I stationary solution is obtained by Newton iteration using an

initial value obtained from time-evolutionary solution

69 / 70References
FreeFem++:
I F. Hecht, FreeFem++ manual, 3rd ed., 2015.
I K. Ohtsuka, T. Takaishi, Finite element analysis bymathematical programming language FreeFem++ (inJapanese), Industrial and Applied Mathematics Series
Vol.4, Kyoritsu, 2014.
I I. Danaila, F. Hecht, O. Pironneau, Simulation numériqueen C++, Dunod, 2003.
Finite element theory:
I D. Braess, Finite elements – Theory, fast solvers andapplication in solid mechanics, 3rd ed., Cambridge Univ.
Press, 2007.
I A. Ern, J.-L. Guermond, Theory and practice of finite
elements, Springer Verlag, New-York, 2004.
I M. Tabata, Numerical solution of partial differential
equations II (in Japanese), Iwanami Shoten, 1994.
70 / 70

example1.edp 1/1// example 1 : poisson-mixedBC.edp[slide page 8]// finite element solution of Poisson equation with mixed boundary condition// for RIIT Tutorial at Kyushu University, 25 Nov.2016, Atsushi Suzuki
int n = 20;
mesh Th=square(n,n);
fespace Vh(Th,P1);
Vh uh,vh;real err, hh;
func f = 5.0/4.0 * pi * pi * sin(pi * x) * sin(pi * y / 2.0);
func h = (-pi)/2.0 * sin(pi * x);
func g = sin(pi * x) * sin(pi * y / 2.0);// for error estimation
func sol = sin(pi * x) * sin(pi * y / 2.0);
func solx = pi * cos(pi * x) * sin(pi * y / 2.0);
func soly = (pi / 2.0) * sin(pi * x) * cos(pi * y / 2.0);
solve poisson(uh,vh) =int2d(Th)( dx(uh)*dx(vh)+dy(uh)*dy(vh) )- int2d(Th)( f*vh )- int1d(Th,1) (h * vh)+ on(2,3,4,uh=g);
hh = 1.0 / real(n) * sqrt(2.0);// int2d uses qf5pT : 5th order integration quadrature
err = int2d(Th)( (dx(uh) - solx) * (dx(uh) - solx) +(dy(uh) - soly) * (dy(uh) - soly) +(uh - sol) * (uh - sol));err = sqrt(err);
cout << "DOF=" << uh[].n << "\t h=" << hh << " err-H1=" << err << endl;
plot(uh,wait=1);

example2.edp 1/2// example 2 : stokes-mixedBC.edp[slide page 12]// error estimation for finite element solution of Stokes equations// with P2/P1 or P1b/P1 element// for RITT Tutorial at Kyushu University, 25 Nov.2016, Atsushi Suzuki
int n1 = 20;
int n2 = n1 * 2;mesh Th1=square(n1,n1);mesh Th2=square(n2,n2);// finite element spaces
fespace Vh1(Th1,P2),Qh1(Th1,P1);
fespace Vh2(Th2,P2),Qh2(Th2,P1);// fespace Vh1(Th1,P1),Qh1(Th1,P1);// fespace Vh2(Th2,P1),Qh2(Th2,P1);// fespace Vh1(Th1,P2),Qh1(Th1,P0);// fespace Vh2(Th2,P2),Qh2(Th2,P0);// external force
func f1 = 5.0/8.0 * pi * pi * sin(pi * x) * sin(pi * y / 2.0) + 2.0 * x;
func f2 = 5.0/4.0 * pi * pi * cos(pi * x) * cos(pi * y / 2.0) + 2.0 * y;func h1 = 3.0/4.0 * pi * sin(pi * x) * cos(pi * y / 2.0);
func h2 = pi * cos(pi * x) * sin(pi * y / 2.0) + x * x + y * y;func sol1 = sin(pi * x) * sin(pi * y / 2.0) / 2.0;
func sol2 = cos(pi * x) * cos(pi * y / 2.0);
func solp = x * x + y * y;func sol1x = pi * cos(pi * x) * sin(pi * y / 2.0) / 2.0;
func sol1y = pi / 2.0 * sin(pi * x) * cos(pi * y / 2.0) / 2.0;
func sol2x = (-pi) * sin(pi * x) * cos(pi * y / 2.0);
func sol2y = (-pi / 2.0) * cos(pi * x) * sin(pi * y / 2.0);// finite element solutions and test functions
Vh1 u11,u12,v11,v12;
Qh1 p1,q1;
Vh2 u21,u22,v21,v22;
Qh2 p2,q2;
real epsln = 1.0e-6;// definitions of macros for strain rate tensor
macro d11(u1)dx(u1) //macro d22(u2)dy(u2) //macro d12(u1,u2) (dy(u1) + dx(u2))/2.0 //// stokes problem
solve stokes1(u11,u12,p1, v11,v12,q1) =int2d(Th1)(
2.0*(d11(u11)*d11(v11)+2.0*d12(u11,u12)*d12(v11,v12)+d22(u12)*d22(v12))- p1*dx(v11) - p1*dy(v12)- dx(u11)*q1 - dy(u12)*q1)- int2d(Th1)(f1*v11+f2*v12) - int1d(Th1,1)(h1*v11+h2*v12)+ on(2,3,4,u11=sol1,u12=sol2);
real meanp,err1,err2,hh1,hh2;plot([u11,u12],p1,wait=1,value=true,coef=0.1);
solve stokes2(u21,u22,p2, v21,v22,q2) =int2d(Th2)(
2.0*(d11(u21)*d11(v21)+2.0*d12(u21,u22)*d12(v21,v22)+d22(u22)*d22(v22))- p2*dx(v21) - p2*dy(v22)- dx(u21)*q2 - dy(u22)*q2)- int2d(Th2)(f1*v21+f2*v22) - int1d(Th2,1)(h1*v21+h2*v22)+ on(2,3,4,u21=sol1,u22=sol2);
plot([u21,u22],p2,wait=1,value=true,coef=0.1);

example2.edp 2/2
hh1 = 1.0 / n1 * sqrt(2.0);
hh2 = 1.0 / n2 * sqrt(2.0);
err1 = int2d(Th1)( (dx(u11) - sol1x) * (dx(u11) - sol1x)+ (dy(u11) - sol1y) * (dy(u11) - sol1y)+ (u11 - sol1) * (u11 - sol1) +(dx(u12) - sol2x) * (dx(u12) - sol2x)+ (dy(u12) - sol2y) * (dy(u12) - sol2y)+ (u12 - sol2) * (u12 - sol2) +(p1 - solp) * (p1 - solp));
err1 = sqrt(err1);
err2 = int2d(Th2)( (dx(u21) - sol1x) * (dx(u21) - sol1x)+ (dy(u21) - sol1y) * (dy(u21) - sol1y)+ (u21 - sol1) * (u21 - sol1) +(dx(u22) - sol2x) * (dx(u22) - sol2x)+ (dy(u22) - sol2y) * (dy(u22) - sol2y)+ (u22 - sol2) * (u22 - sol2) +(p2 - solp) * (p2 - solp));
err2 = sqrt(err2);

cout << "coarse mesh: h=" << hh1 << " err-H1/L2=" << err1 << endl;
cout << "fine mesh:h=" << hh2 << " err-H1/L2=" << err2 << endl;
cout << "O(h)=" << log(err1/err2)/log(hh1/hh2) << endl;

example3.edp 1/1// example 3 : poisson-matrix.edp[slide page 26]// finite element solution of Poisson equation with mixed boundary condition// with P1 element// for RIIT Tutorial at Kyushu University, 25 Nov.2016, Atsushi Suzuki
int n = 20;
mesh Th=square(n,n);
fespace Vh(Th,P1);
Vh u,v;real err, hh;
func f = 5.0/4.0 * pi * pi * sin(pi * x) * sin(pi * y / 2.0);
func h = (-pi)/2.0 * sin(pi * x);
func g = sin(pi * x) * sin(pi * y / 2.0);// for error estimation
func sol = sin(pi * x) * sin(pi * y / 2.0);
func solx = pi * cos(pi * x) * sin(pi * y / 2.0);
func soly = (pi / 2.0) * sin(pi * x) * cos(pi * y / 2.0);
varf aa(u,v) = int2d(Th)( dx(u)*dx(v)+dy(u)*dy(v) )+ on(2,3,4,u=g); // u=1 is enough to say which is Dirichlet node
varf external(u,v) = int2d(Th)( f*v ) + int1d(Th,1) (h * v)
+ on(2,3,4,u=g); // inhomogenoues Dirichlet data are given
real tgv=1.0e+30;
matrix A = aa(Vh,Vh,tgv=tgv,solver=CG);
real[int] ff = external(0,Vh,tgv=tgv); // containing Dirichlet data with tgv
u[] = A^-1 * ff;
hh = 1.0 / real(n) * sqrt(2.0);// int2d uses qf5pT : 5th order integration quadrature
err = int2d(Th)( (dx(u) - solx) * (dx(u) - solx) +(dy(u) - soly) * (dy(u) - soly) +(u - sol) * (u - sol));err = sqrt(err);
cout << "DOF=" << u[].n << "\t h=" << hh << " err-H1=" << err << endl;

example4.edp 1/2// example 4 : cavityNewton.edp[slide page 36]// stational Navier-Stokes equations in a cavity by Newton iteration// P2/P1 element// for RIIT Tutorial at Kyushu University, 25 Nov.2016, Atsushi Suzuki// based on examples++-tutorial/cavityNewtow.edp
mesh Th=square(40,40);
fespace Xh(Th,P2);
fespace Mh(Th,P1);
fespace XXMh(Th,[P2,P2,P1]);
XXMh [u1,u2,p], [v1,v2,q];
macro d11(u1)dx(u1) //macro d22(u2)dy(u2) //macro d12(u1,u2) (dy(u1) + dx(u2))/2.0 //macro div(u1,u2) (dx(u1) + dy(u2))//
macro grad(u1,u2) [dx(u1), dy(u2)]//
macro ugrad(u1,u2,v) (u1*dx(v) + u2*dy(v)) //macro Ugrad(u1,u2,v1,v2) [ugrad(u1,u2,v1), ugrad(u1,u2,v2)]//
real epsln = 1.0e-6;
solve Stokes ([u1,u2,p],[v1,v2,q],solver=UMFPACK) =int2d(Th)(2.0*(d11(u1)*d11(v1) + 2.0*d12(u1,u2)*d12(v1,v2) + d22(u2)*d22(v2))- p * div(v1,v2) - q * div(u1,u2)- p * q * epsln)+ on(3,u1=4*x*(1-x),u2=0)+ on(1,2,4,u1=0,u2=0);
plot(coef=0.2,cmm=" [u1,u2] and p",p,[u1,u2],wait=1);

Mh psi,phi;problem streamlines(psi,phi,solver=UMFPACK) =int2d(Th)( dx(psi)*dx(phi) + dy(psi)*dy(phi))+ int2d(Th)( -phi*(dy(u1)-dx(u2)))+ on(1,2,3,4,psi=0);
streamlines;
plot(psi,wait=1);
real nu=1.0;
XXMh [up1,up2,pp];varf vDNS ([u1,u2,p],[v1,v2,q]) =int2d(Th)(nu * 2.0*(d11(u1)*d11(v1)+2.0*d12(u1,u2)*d12(v1,v2)+d22(u2)*d22(v2))- p * div(v1, v2) - q * div(u1, u2)- p * q * epsln+ Ugrad(u1,u2,up1,up2)’*[v1,v2] //’+ Ugrad(up1,up2,u1,u2)’*[v1,v2]) //’+ on(1,2,3,4,u1=0,u2=0);// [up1, up2, pp] are obtained from the previous step
varf vNS ([u1,u2,p],[v1,v2,q]) =int2d(Th)(nu * 2.0*(d11(up1)*d11(v1)+2.0*d12(up1,up2)*d12(v1,v2)+
d22(up2)*d22(v2))- pp * div(v1, v2) - q * div(up1, up2)- pp * q * epsln+ Ugrad(up1,up2,up1,up2)’*[v1,v2]) //’+ on(1,2,3,4,u1=0,u2=0);
Xh uu1=u1, uu2=u2; // initial condition is given by Stokes eqs : Re=0.
up1[] = 0.0; // initialize for [up1, up2, pp]
real reyini = 100.0;
real reymax = 12800.0;
real re = reyini;
int kreymax = log(reymax / reyini)/log(2.0) * 2.0;
for(int k = 0; k < kreymax; k++) {re *= sqrt(2.0);
real lerr=0.02; // parameter to control mesh refinmente

example4.edp 2/2
if(re>8000) lerr=0.01;
if(re>10000) lerr=0.005;
for(int step= 0 ;step < 2; step++) {Th=adaptmesh(Th, [u1,u2], p, err=lerr, nbvx=100000);[u1, u2, p]=[u1, u2, p];[up1, up2, pp]=[up1, up2, pp];plot(Th, wait=1);
for (int i = 0 ; i < 20; i++) {nu = 1.0 / re;up1[] = u1[];
real[int] b = vNS(0, XXMh);matrix Ans = vDNS(XXMh, XXMh, solver=UMFPACK);
real[int] w = Ans^-1*b;u1[] -= w;
cout << "iter = "<< i << " " << w.l2 << " Reynolds number = " << re<< endl;if(w.l2 < 1.0e-6) break;} // loop : i}// loop : step
streamlines;
plot(psi,nbiso=30,wait=1);// extract velocity component from [u1, u2, p]
uu1 = u1;
uu2 = u2;plot(coef=0.2,cmm="rey="+re+" [u1,u2] and p ",psi,[uu1,uu2],wait=1,nbiso=20);} // loop : re

example5.edp 1/1// example 5 : poisson-LinearCG.edp[slide page 45]// finite element solution of Poisson equation with mixed boundary condition// in matrix form and solved by LinearCG with diagonal preconditioner// for RIIT Tutorial at Kyushu University, 25 Nov.2016, Atsushi Suzuki
int n = 20;
mesh Th=square(n,n);
fespace Vh(Th,P1);
Vh u,v;real err, hh;
func f = 5.0/4.0 * pi * pi * sin(pi * x) * sin(pi * y / 2.0);
func h = (-pi)/2.0 * sin(pi * x);
func g = sin(pi * x) * sin(pi * y / 2.0);// for error estimation
func sol = sin(pi * x) * sin(pi * y / 2.0);
func solx = pi * cos(pi * x) * sin(pi * y / 2.0);
func soly = (pi / 2.0) * sin(pi * x) * cos(pi * y / 2.0);
varf aa(u,v) = int2d(Th)( dx(u)*dx(v)+dy(u)*dy(v) )+ on(2,3,4,u=1.0);
varf external(u,v) = int2d(Th)( f*v ) + int1d(Th,1) (h * v);
real tgv=1.0e+30;
matrix A;
real[int] bc = aa(0, Vh, tgv=tgv);
func real[int] opA(real[int] &pp){
pp = bc ? 0.0 : pp;// SpMV operation only for node in interior ofreal[int] qq = A*pp;// the domain without Dirichlet nodes
pp = bc ? 0.0 : qq;
return pp;}
func real[int] opQ(real[int] &pp){
for (int i = 0; i < pp.n; i++) {pp(i) = pp(i) / A(i, i);}pp = bc ? 0.0 : pp;
return pp;}
A = aa(Vh, Vh, tgv=tgv, solver=sparsesolver);real[int] ff = external(0, Vh);v = g;// g is valid only on the boundary 1u[] = bc ? v[] : 0.0; // u[] has Dirichlet data without tgv
v[] = A * u[];// v[] = A_{12}*gv[] = bc ? 0.0 : v[];
ff -= v[];// ff_{1} -= A_{12}*u_{2}ff = bc ? 0.0 : ff;
LinearCG(opA, u[], ff, precon=opQ, nbiter=200, eps=1.0e-10,verbosity=50);
hh = 1.0 / real(n) * sqrt(2.0);// int2d uses qf5pT : 5th order integration quadrature
err = int2d(Th)( (dx(u) - solx) * (dx(u) - solx) +(dy(u) - soly) * (dy(u) - soly) +(u - sol) * (u - sol));err = sqrt(err);
cout << "DOF=" << u[].n << "\t h=" << hh << " err-H1=" << err << endl;

Schwarz3d.edp 1/5// example6 : Schwarz3d[slide page 54]// Schwarz preconiditoner for Krylov method// for RIIT Tutorial at Kyushu University, 25 Nov.2016, Atsushi Suzuki// based on FreeFem++ manual, section 10 and// scripts in An Introduction to Domain Decomposition Methods -//Algorithms, Theory, and Parallel Implementation,//
V. Dolean, P Jolivet, F. Nataf, SIAM 2015, ISBN 978-1-611974-05-8

load "medit";load "metis";include "cube.idp"
bool flagRAS=false;
int sizeoverlaps=1;// size of overlap

int[int] NN=[50,50,50]; //bool withmetis=true;
int npart= 8;// using numerical diffusion of mass matrix : FreeFem++ manual Examples 11.9
func bool AddLayers(mesh3 & Th,real[int] &ssd,int n){
fespace Vh(Th,P1);
fespace Ph(Th,P0);Ph s;s[]= ssd;
Vh u;varf vM(u,v)=int3d(Th,qforder=1)(u*v/volume);
matrix M=vM(Ph,Vh);for(int i=0;i<n;++i) {u[]= M*s[];u = u > 0.1;s[]= M’*u[]; //’;
s = s > 0.1;}
ssd=s[];
return true;}
func bool SubdomainsPartitionUnity(mesh3 & Th, int nnpart, real[int] & partdof,
int sizeoverlaps,
mesh3[int] & Tha,matrix[int] & Rih, matrix[int] & Dih){fespace Vh(Th,P1);
fespace Ph(Th,P0);
mesh3 Thi=Th;
fespace Vhi(Thi,P1); // FreeFem++ trick, formal definition
Vhi[int] pun(nnpart), dum(nnpart); // local fem functions Vh sun=0, unssd=s0;
Vh sun = 0, unssd = 0, demo = 0;
Ph part;part[]=partdof;
for(int i=0;i<nnpart;++i) {Ph suppi= abs(part - i) < 0.1; // boolean 1 in the subdomain 0 elsewhere
AddLayers(Th,suppi[],sizeoverlaps);// partitions by adding layers
Thi=Tha[i]=trunc(Th,suppi>0,label=10,split=1); // mesh interfaces label 10
Rih[i]=interpolate(Vhi,Vh,inside=true); // Restriction operator : Vh to Vhi
pun[i][] = 1.0;
sun[] += Rih[i]’*pun[i][]; // ’}for(int i=0;i<nnpart;++i) {Thi=Tha[i];//medit("Thi"+i,Thi);
pun[i]= pun[i]/sun;Dih[i]=pun[i][]; //diagonal matrix built from a vector if(verbosity > 1)
dum[i] = (pun[i] == 1.0 ? 0.0 : pun[i]);demo[] += Rih[i]’*dum[i][];}Schwarz3d.edp 2/5
plot(demo,cmm="overlapped skelton",wait=1);
return true;}
real [int,int] BB=[[0,1],[0,1],[0,1]];
int [int,int] L=[[1,1],[1,1],[1,1]];
mesh3 Thg=Cube(NN,BB,L);//medit("Th",Thg);
fespace Ph(Thg,P0);
fespace Vh(Thg,P1);//fespace Xh(Thg,[P1,P1,P1]);
Ph part;Vh sun=0,unssd=0;
Ph xx=x,yy=y;if (withmetis){
int[int] nupart(Thg.nt);
metisdual(nupart,Thg,npart);
for(int n=0;n<nupart.n; n++)part[][n]=nupart[n];}//plot(part,fill=1,cmm="subdomains",wait=false);
mesh3[int] aTh(npart);
matrix[int] Rih(npart);
matrix[int] Dih(npart);
matrix[int] aA(npart);
real[int] partdof(npart);Vh[int] Z(npart);
matrix E;// coarse space : only used as set of arrays

SubdomainsPartitionUnity(Thg, npart, part[], sizeoverlaps, aTh, Rih, Dih);//plot(part,fill=1,cmm="subdomains",wait=1);
macro Grad(u) [dx(u),dy(u),dz(u)] // EOM
func f = 1;// external force
func g = 0 ;// homogeneous Dirichlet data
func kappa = 30.; // viscosity
func eta = 0;
Vh rhsglobal,uglob; // rhs and solution of the global problem
varf vaglobal(u,v) = int3d(Thg)((1.0 + (kappa - 1.0) * x * y * z) *Grad(u)’*Grad(v)) //’+on(1,u=1.0);
varf vexternal(u,v)= int3d(Thg)(f*v);
matrix Aglobal;
real tgv=1.e+30;
real[int] bc = vaglobal(0, Vh, tgv=tgv);
Aglobal = vaglobal(Vh,Vh,tgv=tgv,solver = CG); // global matrix
rhsglobal[] = vexternal(0,Vh); // global rhs
rhsglobal[] = bc ? 0.0 : rhsglobal[];//uglob[] = Aglobal^-1*rhsglobal[];//plot(uglob,value=1,fill=1,wait=1,cmm="Solution by a direct method",dim=3);
for(int n = 0; n < npart; n++) {matrix aT = Aglobal*Rih[n]’; // ’aA[n] = Rih[n]*aT;set(aA[n], solver = sparsesolver);}func real[int] opA(real[int] &v){
v = bc ? 0.0 : v;
real[int] s = Aglobal * v;
s = bc ? 0.0 : s;return s;}
func real[int] opScale(real[int] &v)

Schwarz3d.edp 3/5{
v = bc ? 0.0 : v;
real[int] diag(v.n);diag = Aglobal.diag;real[int] s = v ./ diag;s = bc ? 0.0 : s;return s;// division on each element of array}func bool CoarseSpace(matrix &EE){
for (int n = 0; n < npart; n++) {Z[n] = 1.0;
real[int] zit = Rih[n] * Z[n][];real[int] zitemp = Dih[n] * zit;
Z[n][] = Rih[n]’*zitemp; //’}real[int,int] Ef(npart,npart); //for(int m = 0; m < npart; m++) {real[int] zz = opA(Z[m][]);
for(int n = 0; n < npart; n++) {Ef(m, n) = Z[n][]’*zz; // ’}}
EE = Ef;set(EE,solver=UMFPACK);
return true;}
func real[int] opQ(real[int] &v){
real [int] s(v.n);s = 0.0;
real[int] vv(npart);v = bc ? 0.0 : v;
for(int n = 0; n < npart; n++) {vv[n]= Z[n][]’*v; // ’}real[int] zz = E^-1 * vv;
for(int n = 0; n < npart; n++) {s +=zz[n] * Z[n][];}return s;}
func real[int] opASM(real[int] &v){
real [int] s(v.n);s = 0.0;v = bc ? 0.0 : v;
for (int n = 0; n < npart; n++) {real[int] bi = Rih[n]*v;
real[int] ui = aA[n]^-1*bi; // local solve
if (flagRAS) {bi = Dih[n]*ui;}
else {
bi = ui;// ASM is appropriate for preconditioner}
s += Rih[n]’*bi; // ’}s = bc ? 0.0 : s;return s;}
func real[int] opASM2(real[int] &v){
real[int] s(v.n);s = opQ(v);s += opASM(v);Schwarz3d.edp 4/5
return s;}
func real[int] opP(real[int] &v){
real[int] s(v.n);real[int] s2(v.n);s = opA(v);s2 = opQ(s);s = v - s2;
return s;}// A-orthogonal projection onto Z^T : I - AQ

func real[int] opPt(real[int] &v){
real[int] s(v.n);real[int] s2(v.n);s = opQ(v);s2 = opA(s);s = v - s2;
return s;}
func real[int] opASMQ(real[int] &v){
real[int] s(v.n);s = opQ(v);real[int] ss(v.n);
ss = opPt(v);real[int] sss(v.n);
sss = opASM(ss);ss = opP(sss);
s += ss;
return s;}
Vh un;
CoarseSpace(E);// coarse space Vh[int] Z is also generated here
if (flagRAS) {un = 0.0;
LinearGMRES(opA, un[], rhsglobal[], nbiter=400, precon=opScale,
eps=1.0e-10,verbosity=50);
un = 0.0;
LinearGMRES(opA, un[], rhsglobal[], nbiter=200, precon=opASM,
eps=1.0e-10,verbosity=50);
un = 0;LinearGMRES(opA, un[], rhsglobal[], nbiter=200, precon=opASM2,
eps=1.0e-10,verbosity=50);
un= 0;
LinearGMRES(opA, un[], rhsglobal[], nbiter=200, precon=opASMQ,
eps=1.0e-10,verbosity=50);
plot(un,value=1,fill=1,wait=1,cmm="solution by GMRES-RAS",dim=3);
un = un - uglob;
plot(un,value=1,fill=1,wait=1,cmm="error of GMRES-RAS",dim=3);}
else {
un = 0.0;
LinearCG(opA, un[], rhsglobal[], nbiter=400, precon=opScale,
eps=1.0e-10,verbosity=50);
un = 0.0;
LinearCG(opA, un[], rhsglobal[], nbiter=200, precon=opASM,
eps=1.0e-10,verbosity=50);
un = 0.0;
LinearCG(opA, un[], rhsglobal[], nbiter=200, precon=opASM2,
eps=1.0e-10,verbosity=50);
un = 0.0;

Schwarz3d.edp 5/5
LinearCG(opA, un[], rhsglobal[], nbiter=200, precon=opASMQ,
eps=1.0e-10,verbosity=50);
plot(un,value=1,fill=1,wait=1,cmm="solution by CG-ASM",dim=3);
un = un - uglob;
plot(un,value=1,fill=1,wait=1,cmm="error of CG-ASM",dim=3);}

example7.edp 1/2// example 7 : NS-cylinder.edp[slide page 60]// Navier-Stokes flow around a cylinder// P2/P1 element with Characteristic Galerkin// for RIIT Tutorial at Kyushu University, 25 Nov.2016, Atsushi Suzuki
int n1 = 30;
int n2 = 60;
real nu = 1.0/400.0;
real dt = 0.05;
real alpha = 1.0/dt;
int timestepmax = 400;
border ba(t=0,1.0){x=t*10.0-1.0;y=-1.0;label=1;};
border bb(t=0,1.0){x=9.0;y=2.0*t-1.0;label=2;};
border bc(t=0,1.0){x=9.0-10.0*t;y=1.0;label=3;};
border bd(t=0,1.0){x=-1.0;y=1.0-2.0*t;label=4;};
border cc(t=0,2*pi){x=cos(t)*0.25+0.75;y=sin(t)*0.25;label=5;};
mesh Th=buildmesh(ba(n2)+bb(n1)+bc(n2)+bd(n1)+cc(-n1));plot(Th);
fespace Xh(Th,[P2,P2,P1]);
fespace Vh(Th,P2);
fespace Qh(Th,P1);

Xh [u1,u2,p], [v1,v2,q];
Vh up1, up2;
Qh pp;
macro d11(u1)dx(u1) //macro d22(u2)dy(u2) //macro d12(u1,u2) (dy(u1) + dx(u2))/2.0 //macro div(u1,u2) (dx(u1) + dy(u2)) //Qh psi,phi;func stinflow=y-y*y*y/3.0;
problem streamlines(psi,phi,solver=UMFPACK) =int2d(Th)( dx(psi)*dx(phi) + dy(psi)*dy(phi))+ int2d(Th)( phi*(dy(u1)-dx(u2)))+ on(1,psi=(-2.0/3.0))+ on(4,psi=stinflow)+ on(3,psi=(2.0/3.0))+ on(5,psi=0.0);
streamlines;
plot(psi,wait=1);
problem Stokes([u1,u2,p],[v1,v2,q],solver=UMFPACK) =int2d(Th)(+ 2.0*nu * (d11(u1)*d11(v1)+2.0*d12(u1,u2)*d12(v1,v2)+d22(u2)*d22(v2))- p * div(v1, v2) - q * div(u1, u2)) //+ on(1,3,u2=0)+ on(4,u1=1.0-y*y,u2=0)+ on(5,u1=0,u2=0);
int i;
problem NS([u1,u2,p],[v1,v2,q],solver=UMFPACK,init=i) =int2d(Th)( alpha * (u1*v1 + u2*v2)+ 2.0*nu * (d11(u1)*d11(v1)+2.0*d12(u1,u2)*d12(v1,v2)+d22(u2)*d22(v2))- p * div(v1, v2) - q * div(u1, u2))
- int2d(Th)( alpha * (convect([up1,up2],-dt,up1)*v1+convect([up1,up2],-dt,up2)*v2) )+ on(1,3,u2=0)+ on(4,u1=1.0-y*y,u2=0)+ on(5,u1=0,u2=0);
plot(Th,wait=1);
Stokes;for (i = 0; i < timestepmax; i++) {up1 = u1;

example7.edp 2/2
up2 = u2;
pp = p;NS;
streamlines;
plot(psi, nbiso=30,wait=0);
if (i % 20 == 0) {plot([up1,up2],pp,wait=0,value=true,coef=0.1);}}

example8.edp 1/2// example 8 : RayleighBenard.edp[slide page 69]// Rayleigh-Benard thermal convection in a box// P2/P1/P2 element with Characteristic Galerkin// time evolution data will be stored in "rb.data" for Rayleigh-Benard-stat.edp// for RIIT Tutorial at Kyushu University, 25 Nov.2016, Atsushi Suzuki
int n1 = 80;
int n2 = 20;
real Pr = 0.71;
real Ra = 1500.0;
real dt = 0.01;
real alpha = 1.0/(dt * Pr);
int timestepmax = 600;
mesh Th=square(n1,n2,[x*4.0,y]);fespace Xh(Th,[P2,P2,P1]);
fespace Vh(Th,P2);
fespace Qh(Th,P1);
Xh [u1,u2,p], [v1, v2, q];
Vh up1, up2, th, thp, psi;
Qh pp, ss, rr;macro d11(u1)dx(u1) //macro d22(u2)dy(u2) //macro d12(u1,u2) (dy(u1) + dx(u2))/2.0 //macro div(u1,u2) (dx(u1) + dy(u2)) //real epsln = 1.0e-6;// penalization parameter to avoid pressure ambiguity
int i;
problem NS([u1,u2,p],[v1,v2,q],solver=UMFPACK,init=i) =int2d(Th)(alpha * (u1*v1 + u2*v2)+ 2.0* (d11(u1)*d11(v1)+2.0*d12(u1,u2)*d12(v1,v2)+d22(u2)*d22(v2))- p * div(v1, v2) - q * div(u1, u2)- p * q * epsln)- int2d(Th)( alpha * (convect([up1,up2],-dt,up1)*v1+convect([up1,up2],-dt,up2)*v2) )- int2d(Th) (Ra * thp * v2)+ on(1,3,u2=0)+ on(2,4,u1=0);
problem Heat(th,psi,solver=UMFPACK,init=i) =int2d(Th)(alpha * (th * psi)
+ dx(th) * dx(psi) + dy(th) * dy(psi))- int2d(Th)(alpha * convect([up1, up2], -dt, thp) * psi)+ on(1,th=1)+ on(3,th=0);
problem streamlines(ss,rr,solver=UMFPACK) =int2d(Th)( dx(ss)*dx(rr) + dy(ss)*dy(rr))+ int2d(Th)( rr*(dy(u1)-dx(u2)))+ on(1,2,3,4,ss=0.0);
plot(Th,wait=1);u1[] = 0.0;// impulsive startth = (1.0-y); // conductive solutionfor (i = 0; i < timestepmax; i++) {up1 = u1;
up2 = u2;
pp = p;thp = th;NS;
Heat;plot(th,value=true);streamlines;
plot(ss,nbiso=30);
if (i % 20 == 0) {plot ([u1,u2],value=true,wait=0,coef=0.1);}

example8.edp 2/2}// write time-evolution data for RayleighBenrad-stat.edp
up1=u1;up2=u2;pp=p;
thp=th;{ofstream file("rb.data", binary);
file.precision(16);for (int i = 0; i < up1[].n; i++) {file << up1[](i) << " ";}file << endl;
for (int i = 0; i < up2[].n; i++) {file << up2[](i) << " ";}file << endl;
for (int i = 0; i < pp[].n; i++) {file << pp[](i) << " " ;}file << endl;
for (int i = 0; i < thp[].n; i++) {file << thp[](i) << " ";}file << endl;}

example9.edp 1/2// example 9 : RayleighBenard-stat.edp[slide page 70]// Navier-Stokes flow around a cylinder// P2/P1 element with Characteristic Galerkin// initial data for Newton iteration needs to be prepared in "rb.data"// by Rayleigh-Benard-stat.edp// for RIIT Tutorial at Kyushu University, 25 Nov.2016, Atsushi Suzuki
int n1 = 80;
int n2 = 20;
real Pr = 0.71;
real Ra = 1500.0;
int timestepmax = 600;
mesh Th=square(n1,n2,[x*4.0,y]);fespace Wh(Th,[P2,P2,P1,P2]);
fespace Vh(Th,P2);
fespace Qh(Th,P1);
Wh [u1,u2,p,th], [v1, v2, q, psi];
Vh up1, up2, thp;
Qh pp, ss, rr;macro d11(u1)dx(u1) //macro d22(u2)dy(u2) //macro d12(u1,u2) (dy(u1) + dx(u2))/2.0 //macro div(u1,u2) (dx(u1) + dy(u2)) //macro ugrad(u1,u2,v) (u1*dx(v) + u2*dy(v)) //macro Ugrad(u1,u2,v1,v2) [ugrad(u1,u2,v1), ugrad(u1,u2,v2)]//
real epsln = 1.0e-6;// penalization parameter to avoid pressure ambiguity

varf vDRB ([u1,u2,p,th],[v1,v2,q,psi]) =int2d(Th)( 2.0*(d11(u1)*d11(v1)+2.0*d12(u1,u2)*d12(v1,v2)+d22(u2)*d22(v2))- p * div(v1, v2) - q * div(u1, u2)- p * q * epsln+ (Ugrad(u1,u2,up1,up2)’*[v1,v2] // ’+ Ugrad(up1,up2,u1,u2)’*[v1,v2]) / (2.0 * Pr) //’)- int2d(Th)( Ra * th * v2 )+ int2d(Th)( dx(th) * dx(psi) + dy(th) * dy(psi) )+ int2d(Th)( ugrad(up1,up2,th)*psi + ugrad(u1,u2,thp)*psi)+ on(1,3,u2=0)+ on(2,4,u1=0)+ on(1,th=0)+ on(3,th=0);// [up1, up2, pp] are obtained from the previous step
varf vRB ([u1,u2,p,th],[v1,v2,q,psi]) =int2d(Th)( 2.0*(d11(up1)*d11(v1)+2.0*d12(up1,up2)*d12(v1,v2)+
d22(up2)*d22(v2))- pp * div(v1, v2) - q * div(up1, up2)- pp * q * epsln+ Ugrad(up1,up2,up1,up2)’*[v1,v2] / (2.0 * Pr) ) //’- int2d(Th)( Ra * thp * v2 )+ int2d(Th)( dx(thp) * dx(psi) + dy(thp) * dy(psi) )+ int2d(Th)( ugrad(up1,up2,thp)*psi )+ on(1,3,u2=0)+ on(2,4,u1=0)+ on(1,th=0) // force homogeneous Dirichlet B.C. to be compatible with+ on(3,th=0); // Jacobian matrix of Newton iteration
problem streamlines(ss,rr,solver=UMFPACK) =int2d(Th)( dx(ss)*dx(rr) + dy(ss)*dy(rr))+ int2d(Th)( rr*(dy(up1)-dx(up2)))+ on(1,2,3,4,ss=0.0);
plot(Th,wait=1);// read stationary data computed by RayleighBenrad.edp{
ifstream file("rb.data", binary);example9.edp 2/2
for (int i = 0; i < up1[].n; i++) {file >> up1[](i);}for (int i = 0; i < up2[].n; i++) {file >> up2[](i);}for (int i = 0; i < pp[].n; i++) {file >> pp[](i);}for (int i = 0; i < thp[].n; i++) {file >> thp[](i);}}[u1, u2, p, th] = [up1, up2, pp, thp];plot(th,cmm="temperature is read from the file",value=true, wait=1);
cout << "start newton" << endl;
for (int n = 0; n < 20; n++) {up1 = u1;
up2 = u2;
pp = p;thp = th;plot(thp,cmm="newton n="+n,wait=0,value=true);real[int] b = vRB(0, Wh);
cout << "||b||_l2 " << b.l2 << endl;
matrix A = vDRB(Wh, Wh, solver=UMFPACK);
real[int] w = A^-1 * b;u1[] -= w;
cout << "iter = "<< n << " " << w.l2 << endl;
if (w.l2 < 1.e-6) break;}plot(thp,value=true);up1 = u1;
up2 = u2;
streamlines;
plot(ss,nbiso=30);