1 Plan Friday Lecture 1 Introduction in 2D History Lecture 2: - - PowerPoint PPT Presentation

▶
1 plan
SMART_READER_LITE
LIVE PREVIEW

1 Plan Friday Lecture 1 Introduction in 2D History Lecture 2: - - PowerPoint PPT Presentation

Apr 09, 2024 237 likes •634 views

1 Plan Friday Lecture 1 Introduction in 2D History Lecture 2: Mesh, 3D cases Generalities Installation procedures Tutorial 1: simple cases Syntax Saturday Example 1: Laplaces equation Lecture 4: Fluid mechanics

slide-1
SLIDE 1

1

slide-2
SLIDE 2

Plan

  • History
  • Generalities
  • Installation procedures
  • Syntax
  • Example 1: Laplace’s equation
  • Example 2: Black-Scholes equation
  • Example 3: Schwarz’ algorithm
  • Example 4: fluid-Structure

interaction

  • Example 5: a nonlinear problem
  • Example 6: an optimization problem.

2

Friday Lecture 1 Introduction in 2D à Lecture 2: Mesh, 3D cases Tutorial 1: simple cases Saturday Lecture 4: Fluid mechanics Lecture 5: Advanced topics Tutorial 2: Your own problem

slide-3
SLIDE 3

Introduction to freefem++

Olivier Pironneau LJLL-UPMC

https://dl.dropbox.com/u/6801560/freefemConfUofH.zip

slide-4
SLIDE 4

4

History

freefem has always been an interpreter of a language for partial differential equations

  • 1986: MacFEM – PCFEM written in Pascal
  • 1990: Syntax analyzer (+ D. Bernardi) freefem
  • 1995: freefem+ (OP + Hecht) in C++
  • 1999: freefem++ (spaghetti => Hecht rewrites all)
  • 2000: freefem3D (DelPino, Havé , Pironneau) à
  • 2003: integrated environment freefem++-cs

(Leyaric)

  • 2005: a new documentation (+ Ohtsuka)
  • 2009: freefem++ does 3D à
  • A web site www.freefem.org

– 3D Visualisation by medit (P. Frey)

  • 2011: Parallel version with MPI

Since 2000 freefem’s only author is Frédéric Hecht

slide-5
SLIDE 5
  • Follow the math => variational formulation
  • Algorithms are the user’s responsibility

built from blocks such as elliptic solver and convection operator using Finite Element Methods

  • Automatic mesh generation with adaptivity
  • Improvements follows the research front: if it is FEM it can be done with freefem++

Leading ideas

solve A(u,w) = int2d(th)(u*w/dt + nu*(dx(u)*dx(w)+dy(u)*dy(w))

  • int2d(th)(convect(v,[a1,a2],-dt)/dt + f*w)

+ on(bdy, u=0);

slide-6
SLIDE 6

Example 1: A Dirichlet Problem

Border bdy(t=0, 2*pi) { x = cos(t); y = sin(t); } mesh th = buildmesh(bdy(70)); fespace Vh(th,P2); func f = exp(x) * sin(y) ; Vh u,w; solve A(u,w)=int2d(th)(dx(u)*dx(w)+ dy(u)*dy(w))-int2d(th)(f*w) + on(bdy,u=0); plot(u, ps= “membrane.eps“);

Variational formulation Approximation

try

About speed: freefem++ interprets the formulae to prepare the data for the direct

  • solver. So against a C++ program speed is lost in the formula interpretation only.

About autonomy: freefem++ does not generate an autonomous module!

Leading ideas

slide-7
SLIDE 7

7

Begin with : freefem++-cs

slide-8
SLIDE 8

8

Build your own Integrated Development Environment

  • Download freefem++ and use it like a compiler/interpreter
  • Use your favorite editor + terminal window
  • Fraise (Mac); Crimpson or Notepad++ (Windows)
slide-9
SLIDE 9

linux: use freefem++-cs or recompile the source code

e.g. Port on ubuntu with openGL graphics by L. Dumont

slide-10
SLIDE 10

10

Other freefem++ Tools

  • Documentation in PDF à
  • Web site: www.freefem.org/ff++
  • Mailing list
  • Wiki
  • Examples

Third Edition, Version 3.20 http://www.freefem.org/ff++

F . Hecht, O. Pir

Laboratoire Jacques-Louis Lions, Universit´ e Pierre et Marie Curie, Paris

Note that although it is possible to solve hyperbolic systems there is no special module for it.

slide-11
SLIDE 11

Web site

11

slide-12
SLIDE 12

How to Install freefem++

  • Beginners: Download freefem++-cs
  • Intermediate: Mac OSX > 10.6, Windows > XP:

Download + installer freefem++ from archive Install your own text editor (Fraise, Notepad++, crimpson) Use freefem as a compiler/interpreter

  • Advance: Linux: download+unzip+recompile

http://www.ann.jussieu.fr/~lehyaric/ffcs

slide-13
SLIDE 13

2D Mesh Generation

  • 1. mesh th = square(5,5);
  • //unit square:
  • 2. mesh Th = square(5,10,[x-1.5, 10*y]);//(-0.5,0.5)x(0,10)
  • 1. border a(t=0,2*pi){ x = cos(t); y = sin(t);label=2;}
  • 2. border b(t=0,2*pi){ x =0.5+0.3*cos(-t); y =0.2*sin(-t);}
  • 3. mesh th1 = buildmesh( a(20) + b(10));
  • 4. mesh th2 = buildmesh( a(20) + b(-10));
  • 5. mesh th3 = movemesh(th2,[x+1,y+2]);
  • 6. mesh th4 = readmesh(“mymesh.msh“);
  • 7. func f = sin(x+1);
  • 8. mesh th4 = adaptmesh(th1,f, err=1e-4);

Rule 1: The domain is on the left of its oriented boundary Rule 2: Borders are defined piecewise analytically but must define continuous and closed curves. Rule 3: borders must not overlap nor cross each other. Rule 4: Each border is assigned a number but can be referred by names also. Unless overwritten the number is the order of appearance of the key word «border».

3 4 2 1 try

slide-14
SLIDE 14

Finite Element Spaces in 2D

fespace Vh(th,P1dc); Vh v,vh; varf A(v,vh) = int2d(th)(v*vh/dt/2); varf B(vh,w) =intalledges(th)(vh*mean(w)*(N.x*u1+N.y*u2))

  • int2d(th)( w*(u1*dx(vh)+u2*dy(vh)));

[N.x,N.y]= normal vector Mean(v)=(v+ + v-)/2

P0, P1, P2, P3, P1nc, P1dc,P2dc, P1b,RT0 P03d, P13d, P23d, RT03d, Edge03d, P1b3d

slide-15
SLIDE 15

Boundary Conditions

  • Dirichlet cond by using on(thebdylabel,u=z)
  • Neumann cond are in the variational formulation: int1d(th,2)(nu*g*w)
  • Periodic conditions are within the space definition
  • mesh Th=square(10,15);
  • fespace Vh(Th,P1,periodic= [2,y],[4,y]);
  • Conditions on RT0 elements can be tricky to formulate:
slide-16
SLIDE 16

The language: freefem script

  • Reserved words

– for variables: x,y,z, N.x, N.y, pi … – For objects: int, real, boundary, mesh, fespace, problem, func … – For actions: plot, solve, movemesh, adaptmesh, movemesh … – For programming: do, while, for, if, cout, cin, ofstream … – For mathematical functions and operators: sin(x), dx(), int2d() Tips:

  • Freefem syntax follows C++. A large part of C++ language is implemented.
  • Several external librairies can be used : solvers (UMFPACS, MUMPS..),
  • ptimization: IPOPT, CMAES…)
  • Can launch external programs like gnuplot

16

slide-17
SLIDE 17

Operators

  • fespace Vh(th,P2);
  • Vh u;
  • dx(u), dy(u), dxx(u), dyy(u), dxy(u)
  • convect(u,[a_1,a_2],dt), mean(u), jump(u)
  • You can make your own
  • macro div(u,v) ( dx(u)+dy(v) ) //
  • sin(u), exp(u), ...
  • int2d(u), u[].max, ...

Rule: freefem is an interpreter so expressions are evaluated pointwise as needed. Example:

  • real I = intalledge(th)(sin(dx(u))^2);

is computed as the sum of the values of the integrand at the quadrature points of the edges in a loop over all triangles.

slide-18
SLIDE 18

Quadrature formulae and Solvers

mesh th = square(15,15); fespace Vh(th,P1); Vh u, w, g = x+y; solve a(u,w,solver=LU) = int2d(th)(dx(u)*dx(w)+dy(u)*dy(w)) + int1d(th,qfe=qf1pE) (1.e10* u*w)

  • int1d (th, qfe=qf1pE)(1.e10*g*w) ;
  • Because the quadrature points are the vertices, it is the same as
  • solve a(u,w) = int2d(th)(grad(u)*grad(w))
  • int2d(f,w) + on(th,1,2,3,4,u=0);
  • solver

= LU,CG,Crout,Cholesky,GMRES,sparsesolver, UMFPACK, MUMPS

try

slide-19
SLIDE 19

Syntax: an incomplete extension of C++

mesh Th = square(5,5); fespace Vh(Th,P1); Vh u=0; Vh<complex> uc = x+2i*y ; //complex FE function or array

  • int i = 0 ;

real a=2/5 ;

  • // quiz? value of a?

bool b=(a<2) ; real[int] aa(10) ;

  • // a real array of 10 value

cout<<u(.5,0.6)<<endl ; //value of FE function u at (.5,.6) if(u<1.0) a=2; else a=1; // wrong

  • Vh au = (u<1.0) ? 2.0 : 1.0;
  • fstream ff("myfile.txt");

for(i=0;i<Th.nt;i++) // also while, break, continue for(int j=0;j<3;j++) cout<<Th[i][j].x<<"\t"<<Th[i][j].y<<"\t"<<u[Vh(i,j)]<<endl;

  • for (int i=0 ;i<u[].n ;++i) { u[][i]=1 ;

plot(u,wait=true,dim=3,fill=1,cmm=" v"+i) ; u[][i]=0;}

  • try
slide-20
SLIDE 20

Plots

  • border a(t=0,pi){ x=cos(t); y=sin(t);}
  • border b(t=pi,2*pi){ x=cos(t); y=sin(t);}
  • plot(a(20)+b(40),wait=true);
  • mesh th=buildmesh(a(20)+b(40));
  • plot(th, wait=1,ps=“th.eps“);
  • fespace Vh(th,P2); Vh u=sin(x*y), v=x*exp(-y);
  • plot(u,v,wait=1, value=1,fill=1,dim=3);
  • plot([u,v],wait=1);

More: cut, link with gnuplot and medit

slide-21
SLIDE 21

Summary

A freefem script is xxx.edp and contains

  • a mesh definition: mesh =
  • FEM space déclaration: fespace Vh
  • pde definition : solve a(u,w)= or problem a(u,w)= containing a bilinear part and

a linear part (and boundary conditions) or matrix &right hand side definitions by varf

  • some kind of output, plot, cout…

Example: Darcy flow: -div(A grad p)=0, p= x+ y on top and left sides, p=0 on bottom side and n.(A grad p) = 0 on the right side.

  • 1. mesh th = square(10,10);
  • 2. fespace Vh(th,P1);
  • 3. Vh p,q;
  • 4. real A11=0.5, A22=1.5;
  • 5. solve a(p,q)=int2d(th)(A11*dx(p)*dx(q)+A22*dy(p )*dy(q))

+on(1,p=0)+on(3,4,p=x+y);

  • 1. plot(p, ps= "darcy.eps",wait=true);
  • 2. fespace Wh(th,P0);
  • 3. Wh u=-A11*dx(p ), v= -A22*dy( p );
  • 4. plot([u,v],wait=true,value=true);

21

try

slide-22
SLIDE 22

Time dependent 1d problem: Black-Scholes equation

real Lx=4, Ly=0.1; mesh Th = square(150,2,[Lx*x,Ly*y]); fespace Vh(Th,P1, periodic=[[1,x],[3,x]]); Vh v,u, uold=max(1.-x,0.); plot(Th,wait=1); real dt=0.01, sigmax=0.5, r=0.01, T=1; problem basket(u,v) = int2d(Th) ( u*v*(r+1/dt) + x*x*sigmax^2*dx(u)*dx(v)/2 –(r-sigmax^2)*dx(u)*v)

  • int2d(Th)(uold*v/dt) + on(2,u=0);

for(real t=0;t<T;t+=dt){ basket; uold = u; } plot(u); ofstream fich("u.txt"); for(real t=0;t<=1;t+=0.05){ real c=t*Lx, BS=0.5*c*(1+erf(log(c)/(sigmax*sqrt(2*T))+sigmax*sqrt(T/8.)))

  • 0.5*(1+erf(log(c)/(sigmax*sqrt(2*T))-sigmax*sqrt(T/8.)))+1-c;

real error10 = 10*(abs(u(c,0)-BS); fich<<c<<" "<< u(c,0)<<" "<<BS<<" "<<error10<<endl; } // display by gnuplot: plot "u.txt"using 1:2 w l, ”u.txt"using 1:4 w l

try

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.5 1 1.5 2 2.5 3 3.5 4 "u.txt"using 1:2 "u.txt"using 1:3 "u.txt" using 1:4

slide-23
SLIDE 23

freefem++ and Fluid Dynamics

Olivier Pironneau LJLL-UPMC

https://dl.dropbox.com/u/6801560/OP_ffDay2.zip

slide-24
SLIDE 24

Interpolation: e.g. Schwarz DDM

border a(t=-pi/2,pi/2) { x=cos(t); y=sin(t);} border b(t=pi/2,3*pi/2){ x=cos(t); y=sin(t);} mesh th1=buildmesh(a(10)+b(10)); mesh th2=square(5,5,[2*x,2*y-1]); plot(th1,th2, wait=1); fespace Vh1(th1,P2); Vh1 w1,u1=0; fespace Vh2(th2,P1); Vh2 w2,u2=0; macro Grad(u) [dx(u),dy(u)] // func f=1;

  • int i;
  • // to avoid refactorizing LU

problem L1(u1,w1,solver=LU,init=i)=int2d(th1)(Grad(u1)'*Grad(w1))

  • int2d(th1)(f*w1) + on(b,u1=0)+on(a,u1=u2);

problem L2(u2,w2,solver=LU, init=i) = int2d(th2)(Grad(u2)'*Grad(w2))

  • int2d(th2)(f*w2)
  • + on(4,u2=u1) +on(1,2,3,u2=0);

for(i=0;i<5;i++){ L2; L1; plot(u1,u2,wait=1);}

  • Rule: pointwise evaluation when needed

try

slide-25
SLIDE 25

Stokes Flow: the driven cavity

//file “stokesCavity.edp” mesh Th=square(40,40); fespace Vh(Th,P2) fespace Qh(Th,P1); Vh u,v,uu,vv; Qh p,pp; problem Stokes([u,v,p],[uu,vv,pp]) = int2d(Th)( nu*( dx(u)*dx(uu) + dy(u)*dy(uu) + dx(v)*dx(vv) + dy(v)*dy(vv))

  • p*pp*1.e-6
  • p*(dx(uu) +dy(vv))- pp*(dx(u)+dy(v))

) + on(1,2,4,u=0,v=0) + on(3,u=1,v=0);

  • Stokes;

plot([u,v],p,value=true, wait=true);

  • try

Th = adaptmesh(Th, u,v,p,err=1e-4); stokes; plot( [u,v],p, value=true, wait=true);

3 4 2 1

slide-26
SLIDE 26

Time Dependent Stokes Flow

//file timeStokesCavity.edp mesh Th=square(40,40); fespace Vh(Th,P2), Qh(Th,P1); Vh u,v,uu,vv, uold,vold; Qh p,pp; real nu=0.01, T=1., dt = 0.1; int m, M= T/dt;

  • problem Stokes([u,v,p],[uu,vv,pp], init=m, solver=LU) =

int2d(Th)( (u*uu+v*vv)/dt + nu*(dx(u)*dx(uu) + dy(u)*dy(uu) + dx(v)*dx(vv) + dy(v)*dy(vv))

  • p*pp*1.e-6 - p*(dx(uu) +dy(vv))- pp*(dx(u)+dy(v))
  • int2d(Th)((uold*uu+vold*vv)/dt)

+ on(1,2,4,u=0,v=0) + on(3,u=1,v=0); for(m=0;m<M;m++){ Stokes; uold=u; vold=v; plot(p,value=true, fill=true, wait=false); }

  • try

3 4 2 1

slide-27
SLIDE 27

Navier - Stokes Flow (I)

//file NavierStokesCavity.edp mesh Th=square(40,40); fespace Vh(Th,P2), Qh(Th,P1); Vh u,v,uu,vv, uold,vold; Qh p,pp; real nu=0.01, T=1., dt = 0.1; int m, M= T/dt;

  • problem NavierStokes([u,v,p],[uu,vv,pp], init=m, solver=LU) =

int2d(Th)( (u*uu+v*vv)/dt + nu*(dx(u)*dx(uu) + dy(u)*dy(uu) + dx(v)*dx(vv) + dy(v)*dy(vv))

  • p*pp*1.e-6 - p*(dx(uu) +dy(vv))- pp*(dx(u)+dy(v))
  • + (uold*dx(u)+vold*dy(u))*uu
  • + (uold*dx(v)+vold*dy(v))*vv

) - int2d(Th)((uold*uu+vold*vv)/dt) + on(1,2,4,u=0,v=0) + on(3,u=1,v=0); for(m=0;m<M;m++){ NavierStokes; uold=u; vold=v; plot(p,value=true, fill=true, wait=false); }

  • try

3 4 2 1

slide-28
SLIDE 28

Navier - Stokes Flow (II)

//file NSCavity.edp mesh Th=square(40,40); fespace Vh(Th,P2), Qh(Th,P1); Vh u,v,uu,vv, uold,vold; Qh p,pp; real nu=0.01, T=1., dt = 0.1; int m, M= T/dt;

  • problem NS([u,v,p],[uu,vv,pp], init=m, solver=Crout) =

int2d(Th)( (u*uu+v*vv)/dt + nu*(dx(u)*dx(uu) + dy(u)*dy(uu) + dx(v)*dx(vv) + dy(v)*dy(vv))

  • p*pp*1.e-6 - p*(dx(uu) +dy(vv))- pp*(dx(u)+dy(v)))
  • int2d(Th) (( uold(x-uold(x,y)*dt,y-uold(x,y)*dt)*uu
  • + vold(x-uold(x,y)*dt,y-uold(x,y)*dt)*vv )/dt)

+ on(1,2,4,u=0,v=0) + on(3,u=1,v=0); for(m=0;m<M;m++){ NS; uold=u; vold=v; plot(p,value=true, fill=true, wait=false); }

  • try

3 4 2 1

slide-29
SLIDE 29

Navier - Stokes Flow (III)

//file NSCylinder.edp mesh Th=square(40,40); fespace Vh(Th,P2), Qh(Th,P1); Vh u,v,uu,vv, uold,vold; Qh p,pp; real nu=0.01, T=1., dt = 0.1; int m, M= T/dt;

  • problem NS([u,v,p],[uu,vv,pp], init=m, solver=Crout) =

int2d(Th)( (u*uu+v*vv)/dt + nu*(dx(u)*dx(uu) + dy(u)*dy(uu) + dx(v)*dx(vv) + dy(v)*dy(vv))

  • p*pp*1.e-6 - p*(dx(uu) +dy(vv))- pp*(dx(u)+dy(v)))
  • int2d(Th) ((convect([uold,vold],-dt,uold)*uu

+ convect([uold,vold],-dt,vold)*vv )/dt ) + on(1,2,3,u=0,v=0) + on(4,u=1,v=0); for(m=0;m<M;m++){ NS; uold=u; vold=v; plot(p,value=true, fill=true, wait=false); }

try 4 1 1 2 3

slide-30
SLIDE 30

Vh r,rr, rold; real kappadt= 0.1*dt; problem NS([u,v,p],[uu,vv,pp], init=m, solver=Crout) = int2d(Th)( (u*uu+v*vv)/dt + nu*(dx(u)*dx(uu) + dy(u)*dy(uu) + dx(v)*dx(vv) + dy(v)*dy(vv))

  • p*pp*1.e-6 - p*(dx(uu) +dy(vv))- pp*(dx(u)+dy(v)))
  • int2d(Th) ((convect([uold,vold],-dt,uold)*uu

+ convect([uold,vold],-dt,vold)*vv )/dt –g*r) + on(1,2,3,4,u=0,v=0);

  • for(m=0;m<M;m++){

NS; thermal; uold=u; vold=v; rold=r; plot(r,value=true, fill=true, wait=false); }

  • Boussinesq

try 3 4 2 1

problem thermal(r,rr,solver=Crout, init=m) = int2d(Th)(r*rr + kappadt*(dx(r)*dx(rr)+dy(r)*dy(rr)))

  • int2d(Th)(convect([uold,vold],-dt,rold)*rr)

+ on(2,r=-1)+on(4,r=3);

slide-31
SLIDE 31

//file floatincavity.edp real x0=0.25,y0=0.25, R=0.15; border c(t=0,2*pi){x=x0+R*cos(t);y=y0+R*sin(t);} mesh Th=buildmesh(a1(n)+a2(n)+a3(n)+a4(n)+c(-n)); fespace Vh(Th,P1b); Vh u,uu, v,vv, uold=0, vold=0,U,V; fespace Ph(Th,P1); Ph p=0,pp; real u0=0,v0=0; for(int m=0;m<M;m++){ real Fx=-nu*int1d(Th,c)((2*dx(u)-p/nu)*N.x+(dy(u)+dx(v))*N.y); real Fy=-nu*int1d(Th,c)((dy(u)+dx(v))*N.x+(2*dy(v)-p/nu)*N.y); real om=nu*int1d(Th,c)(((2*dx(u)-p/nu)*N.x+(dy(u)

  • +dx(v))*N.y)*(y-y0)
  • ((dy(u)+dx(v))*N.x+(2*dy(v)-p/nu)*N.y)*(x-x0));

u0+=dt*Fx/mu; v0+=dt*Fy/mu; x0+=dt*u0; y0+=dt*v0; U=u0+om*(y-y0); V=v0-om*(x-x0); solve NS([u,v,p],[uu,vv,pp]) = int2d(Th)(... +on(c,u=U,v=V); Th=buildmesh(a1(n)+a2(n)+a3(n)+a4(n)+c(-n)); uold=u;vold=v; }

Navier - Stokes Flow (III) try 1 1 2 c 1

slide-32
SLIDE 32

32

border circle(t=0,2*pi){x=cos(t); y=sin(t);label=1;} mesh Thx = buildmesh (circle(20)); int[int] rup=[0,3], rdown=[0,2], rmid=[1,1 ]; mesh3 th=buildlayers(Thx,30*n,zbound=[0,5],labelmid=rmid,labelup=rup, labeldown = rdown,transfo=[x,(4+y)*cos(z*pi/10)-4,(4+y)*sin(z*pi/10)]); fespace Vh(th,P2); fespace Qh(th,P1); fespace Zh(Thx,P2); Zh w2; Vh u,v,w,uh,vh,wh,uold=0,vold=0,wold=0; Qh p,ph; real p1=1, dt=0.1, nu=0.005; int n=-1; problem a([u,v,w,p],[uh,vh,wh,ph],init=n)= int3d(th)((u*uh+v*vh+w*wh)/dt + nu*( dx(u)*dx(uh)+dy(u)*dy(uh)+dz(u)*dz(uh) + dx(v)*dx(vh)+dy(v)*dy(vh)+dz(v)*dz(vh) + dx(w)*dx(wh)+dy(w)*dy(wh)+dz(w)*dz(wh) )

  • (dx(uh)+dy(vh)+dz(wh))*p - (dx(u)+dy(v)+dz(w))*ph + p*ph*1e-5 )
  • int3d(th)(( convect([uold,vold,wold],-dt,uold)*uh
  • +convect([uold,vold,wold],-dt,vold)*vh
  • +convect([uold,vold,wold],-dt,wold)*wh )/dt)
  • int2d(th,2)(p1*wh) - int2d(th,3)(p1*vh)
  • + on(1,u=0,v=0,w=0)+ on(2,v=0,u=0)+ on(3,u=0,w=0);

for(real t=0;t<2.6;t+=dt){ p1=cos(t); n+=1; a; uold=u; vold=v; wold=w;}

  • Blood Flow

try

slide-33
SLIDE 33

Wave equation try

mesh Th = readmesh("chesapeake.msh"); real dt=3, dt2=dt*dt, nu=0.01; fespace Vh(Th,P1); Vh u,v,uold,uoldold;

  • u= 1+2*exp(-((y-400)^2+(x-300)^2)/300); uold= u; uoldold = u;
  • problem wave(u,v)=

int2d(Th)( u*v/dt2 + nu*( dx(u)*dx(v)+dy(u)*dy(v) ))

  • int2d(Th)((2*uold-uoldold)*v/dt2 );// + on(2,u=0);
  • real[int] visual(20);

for(int k=0;k<visual.n;k++) visual[k]=2.0*k/visual.n;

  • for(int k=0;k<500;k++) {

wave; uoldold=uold; uold=u; plot(u,wait=0,value=1,fill=1, viso=visual, nbiso=visual.n); }

slide-34
SLIDE 34

Shallow-Water equations try

mesh Th = readmesh("chesapeake.msh"); real dt=1, nu=0.01; fespace Wh(Th,P1); fespace Vh(Th,P1); Wh u,v,u1,v1,uh,vh; Vh r,rh,r1, fro;

  • macro dn(u) (N.x*dx(u)+N.y*dy(u) ) // normal derivative

r1= 1+2*exp(-((y-400)^2+(x-300)^2)/300); v1= 0; u1 = 0;

  • problem euler(u,v,r,uh,vh,rh)=

int2d(Th)( (u*uh+v*vh+r*rh)/dt + dx(r)*uh+dy(r)*vh + r*rh*(dx(u1)+dy(v1)) +nu*( dx(u)*dx(uh)+dy(u)*dy(uh)+dx(v)*dx(vh+dy(u)*dy(uh)))

  • int2d(Th)((rh*convect([u1,v1],-dt,r1)
  • + uh*convect([u1,v1],-dt,u1
  • + vh*convect([u1,v1],-dt,v1))/dt)

+on(2,r=1,u=0,v=0); for(int k=0;k<100;k++) { euler; u1=u; v1=v; r1=r;}

slide-35
SLIDE 35

Perspectives

  • FreeFem++ is easy to use for simple problems and hard
  • n complex problem (the no-free lunch theorem)
  • Now 3D but speed is an issue: parallel version
  • Sensitivity, optimisation, eigenvalues, matrix form,
  • ptimal control, mesh adaptivity, etc?
slide-36
SLIDE 36

Multi-Physics

  • mesh th=square(20,10,[x,y/4+1]);

fespace Vh(th,P2); Vh u,v,uu,vv; mesh Th=square(20,20); fespace Uh(Th,P1b); Uh uf,vf,uuf,vvf; fespace Ph(Th,P1); Ph p,pp; solve stokes([uf,vf,p],[uuf,vvf,pp]) = int2d(Th)(dx(uf)*dx(uuf)+dy(uf)*dy(uuf) + dx(vf)*dx(vvf)+ dy(vf)*dy(vvf) + dx(p)*uuf + dy(p)*vvf + pp*(dx(uf)+dy(vf))) + on(1,2,4,uf=0,vf=0) + on(3,uf=1,vf=0);

  • macro epsilon(u1,u2) [dx(u1),dy(u2),(dy(u1)+dx(u2))/sqrt(2.)] //

macro div(u,v) ( dx(u)+dy(v) ) // real E=21e5, nu=0.28, mu=E/(2*(1+nu)), lambda=E*nu/((1+nu)*(1-2*nu)), f=-1; solve lame([u,v],[uu,vv])= int2d(th)( lambda*div(u,v)*div(uu,vv) +2*mu*( epsilon(u,v)'*epsilon(uu,vv) ) ) - int2d(th)(f*vv) +int1d(th,1)(50*p*vv) + on(2,3,4,u=0,v=0); th = movemesh(th,[x,y+400*v]); Th = movemesh(Th, [x, y+400*y*v(x,1.0)]); u=u; v=v; uf=uf;vf=vf; plot(v,[uf,vf],wait=0);

try

slide-37
SLIDE 37

Non-linear problem Try the fixed point scheme:

mesh th = square(10,10); fespace Vh(th,P1); func f= exp(x) * sin(y) ; Vh u,w, uold=x*(x-1)*y*(y-1);

  • problem A(u,w,solver=LU)

=int2d(th)((1+uold^3) *(dx(u)*dx(w)+ dy(u)*dy(w)))

  • int2d(th)(f*w) + on(1,2,3,4,u=0);
  • for(int m=0;m<5;m++){

A; w=u-uold; plot(w,wait=true,value=true); uold=u; }

try

slide-38
SLIDE 38

Optimization e.g.

A better method is to solve, with q=p+1 : mesh th = square(10,10); fespace Vh(th,P1); fespace Ph(th,P0); func f=1; func real F(real v){return (1+v^2)^4;}//v ->|grad(u)| func real dF(real v){return 8*(1+v^2)^3;} func real J(real[int] & u) { Vh w; w[]=u; // copy array u in the FEM function w return int2d(th)(F( dx(w)^2 + dy(w)^2 ) - 8*f*w) ; } func real[int] dJ(real[int] & u) { Vh w;w[]=u; Ph rho=dF( dx(w)^2 + dy(w)^2); varf au(uh,vh)=int2d(th)(rho*(dx(w)*dx(vh)+dy(w)*dy(vh))

  • 8*f*vh) + on(1,2,3,4,uh=0);

u= au(0,Vh); //above with vh replaced by the ith hat function return u; } real[int] u(th.nv); u=0; BFGS(J,dJ,u,eps=1.e-6,nbiter=10,nbiterline=10); Vh w; w[]=u; plot(w);

  • try