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

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: Laplace’s equation Lecture 4: Fluid mechanics • Example 2: Black-Scholes equation Lecture 5: Advanced topics • Example 3: Schwarz’ algorithm Tutorial 2: Your own problem • Example 4: fluid-Structure interaction • Example 5: a nonlinear problem • Example 6: an optimization problem. 2

Introduction to freefem++ Olivier Pironneau LJLL-UPMC https://dl.dropbox.com/u/6801560/freefemConfUofH.zip

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 4

Leading ideas • 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++ 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); �

Leading ideas Example 1: A Dirichlet Problem Variational formulation Approximation 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 “ ); � 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! try

Begin with : freefem++-cs 7

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

linux: use freefem++-cs or recompile the source code e.g. Port on ubuntu with openGL graphics by L. Dumont

Other freefem++ Tools • Documentation in PDF à • Web site: www.freefem.org/ff++ • Mailing list Third Edition, Version 3.20 • Wiki http://www.freefem.org/ff++ • Examples 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. 10

Web site 11

How to Install freefem++ • Beginners : Download freefem++-cs http://www.ann.jussieu.fr/~lehyaric/ffcs • 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 �

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)); � 3 4. mesh th2 = buildmesh ( a(20) + b(-10)); � 5. mesh th3 = movemesh (th2,[x+1,y+2]); � 4 2 6. mesh th4 = readmesh ( “ mymesh.msh “ ); � 1 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 try 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».

Finite Element Spaces in 2D P0, P1, P2, P3, P1nc, P1dc,P2dc, P1b,RT0 � P03d, P13d, P23d, RT03d, Edge03d, P1b3d � 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 �

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:

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..), optimization: IPOPT, CMAES…) • Can launch external programs like gnuplot 16

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

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

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; � ofstream 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

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

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 ); 21 try 4. plot([u,v],wait=true,value=true);