FreeFem++, a generic tools to solve mechanical PDE with finite elements F. Hecht Laboratoire Jacques-Louis Lions Universit´ e Pierre et Marie Curie Paris, France with O. Pironneau, J. Morice, P. Ventura, O. Pantz http://www.freefem.org mailto:hecht@ann.jussieu.fr With the support of ANR (French gov.) ANR-07-CIS7-002-01 http://www.freefem.org/ff2a3/ http://www-anr-ci.cea.fr/ ECCM2010, Paris, Mai 2010 1
PLAN – Examples : – Introduction Freefem++ – Coupling ´ equation : Elasticy, Piezo, – some History Dielectic, and BEM, – Linear elastic Problem, Varia- – Vessel tional formulation – Cantilever – Eigen Value problem – Conclusion / Future – http://www.freefem.org/ ECCM2010, Paris, Mai 2010 2
Introduction FreeFem++ is a software to solve numerically partial differential equations R 2 ) and in I R 3 ) with finite elements methods. We used a user (PDE) in I language to set and control the problem. The FreeFem++ language allows for a quick specification of linear PDE’s, with the variational formulation of a linear steady state problem and the user can write they own script to solve no linear problem and time depend problem. You can solve coupled problem or problem with moving domain or eigenvalue problem, do mesh adaptation , compute error indicator, etc ... FreeFem++ is a freeware and this run on Mac, Unix and Window architecture, in parallel with MPI. ECCM2010, Paris, Mai 2010 3
History 1987 MacFem/PCFem les ancˆ etres (O. Pironneau en Pascal) payant. 1992 FreeFem r´ e´ ecriture de C++ (P1,P0 un maillage) O. Pironneau, D. Bernardi, F. Hecht , C. Prudhomme (adaptation Maillage, bamg). 1996 FreeFem+ r´ e´ ecriture de C++ (P1,P0 plusieurs maillages) O. Piron- neau, D. Bernardi, F. Hecht (alg` ebre de fonction). 1998 FreeFem++ r´ e´ ecriture avec autre noyau ´ el´ ement fini, et un autre lan- gage utilisateur ; F. Hecht, O. Pironneau, K.Ohtsuka. 1999 FreeFem 3d (S. Del Pino) , Une premi` ere version de freefem en 3d avec des m´ ethodes de domaine fictif. 2008 FreeFem++ v3 r´ e´ ecriture du noyau ´ el´ ement fini pour prendre en compte les cas multidimensionnels : 1d,2d,3d... ECCM2010, Paris, Mai 2010 4
For who, for what ! For what 1. R&D 2. Academic Research , 3. Teaching of FEM, PDE, Weak form and variational form 4. Algorithmes prototyping 5. Numerical experimentation 6. Scientific computing and Parallel computing For who : the researcher, engineer, professor, student... ECCM2010, Paris, Mai 2010 5
The main characteristics of FreeFem++ I/II (2D) – Wide range of finite elements : linear (2d,3d) and quadratic Lagrangian (2d,3d) elements, discontinuous P1 and Raviart-Thomas elements (2d,3d), 3d Edge element , vectorial element , mini-element( 2d, 3d), ... – Automatic interpolation of data from a mesh to an other one, so a finite element function is view as a function of ( x, y, z ) or as an array, with matrix construction if need. – Definition of the problem (complex or real value) with the variational form with access to the vectors and the matrix if needed – Discontinuous Galerkin formulation (only in 2d to day). ECCM2010, Paris, Mai 2010 6
The main characteristics of FreeFem++ II/II (2D) – Analytic description of boundaries, with specification by the user of the intersection of boundaries in 2d. – Automatic mesh generator, based on the Delaunay-Vorono ¨ ı algorithm. (2d,3d) – load and save Mesh, solution – Mesh adaptation based on metric, possibly anisotropic, with optional auto- matic computation of the metric from the Hessian of a solution. – LU, Cholesky, Crout, CG, GMRES, UMFPack, SuperLU, MUMPS, HIPS, HYPRE, SUPERLU DIST, PASTIX. ... sparse linear solver ; eigenvalue and eigenvector computation with ARPACK. – Online graphics with OpenGL/GLUT, C++ like syntax. – Link with other soft : modulef, emc2, medit, gnuplot, tetgen, superlu, mumps, ... – Dynamic linking to add plugin. – Wide range of of examples : Navier-Stokes 3d, elasticity 3d, fluid structure, eigenvalue problem, Schwarz’ domain decomposition algorithm, residual er- ror indicator, ... ECCM2010, Paris, Mai 2010 7
How to use on Unix build a ”yours.edp” file with your favorite editor : emacs, vi, nedit, etc. Enter FreeFem++ yours.edp or FreeFem++-nw yours.edp to execute your script. Remark, this application FreeFem++ must be in a directory of your PATH shell variable. on Window, MacOs X build a ”yours.edp” file with your favorite text editor (raw text, not word text) : emacs, winedit, wordpad, bbedit, ... and click on the icon of the application FreeFem++ and load you file via de open file dialog box or drag and drop the icon of your built file on the application FreeFem++ icon. ECCM2010, Paris, Mai 2010 8
Linear elasticty equation, weak form Let a domain Ω ⊂ R d with a partition of ∂ Ω in Γ d , Γ n . Find the displacement u field such that : − ∇ .σ ( u ) = f in Ω , u = 0 on Γ d , σ ( u ) . n = 0 on Γ n (1) Where ε ( u ) = 1 2 ( ∇ u + t ∇ u ) and σ ( u ) = A ε ( u ) with A the linear positif operator on symmetric d × d matrix corresponding to the material propriety. Denote V g = { v ∈ H 1 (Ω) d / v | Γ d = g } . The Basic displacement variational formulation with is : find u ∈ V 0 (Ω) , such that � � � Ω ε ( v ) : A ε ( u ) = Ω v . f + Γ (( A ε ( u )) .n ) .v, ∀ v ∈ V 0 (Ω) (2) ECCM2010, Paris, Mai 2010 9
Linear elasticty equation, in FreeFem++ The finite element method is just : replace V g with a finite element space, and the FreeFem++ code : load "medit" include "cube.idp" int[int] Nxyz=[20,5,5]; real [int,int] Bxyz=[[0.,5.],[0.,1.],[0.,1.]]; int [int,int] Lxyz=[[1,2],[2,2],[2,2]]; mesh3 Th=Cube(Nxyz,Bxyz,Lxyz); // Alu ... real rhoAlu = 2600, alu11= 1.11e11 , alu12 = 0.61e11 ; real alu44= (alu11-alu12)*0.5; func Aalu = [ [alu11, alu12,alu12, 0. ,0. ,0. ], [alu12, alu11,alu12, 0. ,0. ,0. ], [alu12, alu12,alu11, 0. ,0. ,0. ], [0. , 0. , 0. , alu44,0. ,0. ], [0. , 0. , 0. , 0. ,alu44,0. ], [0. , 0. , 0. , 0. ,0. ,alu44] ]; real gravity = -9.81; ECCM2010, Paris, Mai 2010 10
Linear elasticty equation, in FreeFem++ The problem code : fespace Vh(Th,[P1,P1,P1]); Vh [u1,u2,u3], [v1,v2,v3]; macro Strain(u1,u2,u3) [dx(u1),dy(u2),dz(u3),(dz(u2)+dy(u3)), (dz(u1)+dx(u3)),(dy(u1)+dx(u2))] // EOM macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM solve Lame([u1,u2,u3],[v1,v2,v3])= int3d(Th)( Strain(v1,v2,v3)’*(Aalu*Strain(u1,u2,u3)) ) - int3d(Th) ( rhoAlu*gravity*v3) + on(1,u1=0,u2=0,u3=0) ; real coef= 0.1/u1[].linfty; int[int] ref2=[1,0,2,0]; mesh3 Thm=movemesh3(Th,transfo=[x+u1*coef,y+u2*coef,z+u3*coef],label=ref2); plot(Th,Thm, wait=1,cmm="coef amplification = "+coef ); medit("Th-Thm",Th,Thm); ECCM2010, Paris, Mai 2010 11
Lame equation / figure Execute beam-A-3d.edp Execute fish.edp Execute EqPoisson.edp ECCM2010, Paris, Mai 2010 12
Linear elasticty equation Eigen Value The code to compute and show eigen vector of the Linear elasticty equation varf vLame([u1,u2,u3],[v1,v2,v3])= int3d(Th)( Strain(v1,v2,v3)’*(Aalu*Strain(u1,u2,u3)) )+ on(1,u1=0,u2=0,u3=0) ; varf vMass([u1,u2,u3],[v1,v2,v3])= int3d(Th)( [v1,v2,v3]’*[u1,u2,u3]*rhoalu ); matrix A= vLame(Vh,Vh,solver=sparsesolver); matrix B= vMass(Vh,Vh,solver=CG,eps=1e-20); int nev=20; // number of computed eigen value close to 0 real[int] ev(nev); // to store nev eigein value Vh[int] [eu1,eu2,eu3](nev); // to store nev eigen vector int k=EigenValue(A,B,sym=true,value=ev,vector=eu1,tol=1e-10); nev=min(k,nev); // some time the number of converged EV. can be greater than nev; for (int i=0;i<nev;i++) { real coef= 0.5/eu1[i][].linfty; mesh3 Thm=movemesh3(Th,transfo=[x+eu1[i]*coef,y+eu2[i]*coef,z+eu3[i]*coef]); medit("Thm-"+ev[i],Thm,wait=1);} Execute beam-EV-3d.edp ECCM2010, Paris, Mai 2010 13
Some Trick to build meshes The problem is to compute eigenvalue of a potential flow on the Chesapeake bay (Thank to Mme. Sonia Garcia, smg @ usna.edu) – Read the image in freefem, adaptmesh , trunc to build a first mesh of the bay and finally remove no connected component. We use : ξ > 0 . 9 || ξ || ∞ where ξ is solution of ∂ξ 10 − 5 ξ − ∆ ξ = 0 in Ω; ∂n = 1 on Γ . Remark, on each connect component ω of Ω, we have ∂ω 1 � ξ | ω ≃ 10 5 ω 1 . � Execute Chesapeake/Chesapeake-mesh.edp – Solve the eigen value, on this mesh. – Execute Chesapeake/Chesapeake-flow.edp ECCM2010, Paris, Mai 2010 14
A true problem with 3 physics : Piezo, Elastic, Dielectric In domain Ω =]0 , p [ × R ; split in 3 part Ω D , Ω E , Ω P (resp. Dielectric part, Elastic part, Pieze Electric part). We solve a follow problem : Compute u , φ in Ω such that : ∀ ( v , ψ ) � � ε ( v ) ∗ : A ε ( u ) ǫE ( φ ) ∗ E ( ψ ) + Ω D Ω E � − ǫE ( φ ) ∗ E ( ψ ) + ε ( v ) ∗ : A ε ( u ) − ε ( v ) ∗ : E E ( φ ) − E ( ψ ) ∗ E ∗ : ε ( u ) + Ω P � Ω ω 2 ρ v ∗ u − = 0 where E ( φ ) = −∇ φ , ∗ is the transpose and complex conjugation, ε ( u ) = 1 2 ( ∇ u + t ∇ u ) , A the linear positif operator on symmetric d × d matrix corresponding to the material elactic propriety E is the piezo/elastic coupling tensor, ω is the pulse, ρ the density, ǫ the permitivity material. Boundary condition are : α p -periodic (i.e : u ( x + p, . ) = αu ( x, . ) and φ given on the elastic boundary. ECCM2010, Paris, Mai 2010 15
Recommend
More recommend
