Gypsilab : a M ATLAB toolbox for FEM-BEM coupling Franois Alouges, - PowerPoint PPT Presentation

Gypsilab : a M ATLAB toolbox for FEM-BEM coupling François Alouges, (joint work with Matthieu Aussal) Workshop on Numerical methods for wave propagation and applications Sept. 1st 2017

Facts New numerical techniques in FEM (e.g. DDM) or BEM (compression techniques). Difficult to test them on different applications (equations might be quite complex). Numericians might not be specialists in fluid/elasticity/electromagnetic/acoustic/etc. domains. Difficult to transfer one technique that worked in a domain to another one. Difficult to transfer the new techniques to other people (e.g. companies).

FEM and BEM specificities FEM leads to sparse matrices N unknowns O ( N ) storage (2D-3D) FEM softwares are numerous. Some are free (e.g. FreeFem++, GetDP , Xlife++, FENICS), others are commercial products (COMSOL, Fluent, etc.). BEM leads to O ( N 2 ) storage (usually on the boundary 2D) BEM softwares are more rare (BEM++, Xlife++) BEM methods often need compression techniques (FMM, H -Matrices, SCSD) and singular integration techniques. This is often very technical. FEM-BEM coupling is even more difficult (coupling sparse and H -Matrices, direct vs iterative methods, etc.).

Prototyping languages When one creates a new numerical technique it is practical to have a prototyping environment (Matlab, Python, Julia). FreeFem++ (Nice for FEM, no BEM yet) FENICS/FireDrake (FEM in Python, BEM++ is being interfaced) Julia ? Matlab : everybody has his own library for FEM, some people implement BEM Need for a versatile easy to use toolbox “à la FreeFem++” for FEM-BEM

Gypsilab On-going work Written in full Matlab Compatible with the H -Matrices toolbox developed by M. Aussal Goals : Simple to express variational formulations that may involve convolution kernels. Efficient enough to prototype moderately large problems (say O ( 10 6 ) unknowns in FEM and O ( 10 5 ) unknowns in BEM). larger problems in progress...

Formalism Assembling matrices (LHS) and vectors (RHS) coming from variational formulations Leave the solve to the user (or eigenvalue, or non linear solving, etc.) Example 1: Laplace problem � − ∆ u + u = f in Ω ∂ u ∂ n = 0 on ∂ Ω

Laplace example % Finite elements u = fem(’P1’); v = fem(’P1’); % Matrix and RHS K = integral(Omega,grad(v),grad(u)) ... + integral(Omega,v,u); F = integral(Omega, v, f); % Resolution uh = K \ F; % erreur en norme L2 et H1 errL2 = u.diff(Omega, uh, Uex, ’L2’); errH1 = u.diff(Omega, uh, Uex, ’H1’); Syntax: integral(domain, test, fn, unknown)

Full program

Result error L2, CB Neumann error H1, CB Neumann −1 0 10 10 −2 10 −1 10 −3 10 −2 10 −4 10 EF P1 EF P1 EF P2 EF P2 slope 2 slope 1 slope 3 slope 2 −5 −3 10 10 −0.8 −0.3 −0.8 −0.3 10 10 10 10 h h

Formalism Example 2: Laplace problem with Fourier BC � − ∆ u + u = 0 in Ω ∂ u ∂ n + u = 1 on ∂ Ω � � � ( ∇ u · ∇ v + uv ) dx + uv ds = 1 · v ds Ω ∂ Ω ∂ Ω

Full program

Result 0.7 0.68 0.66 0.64 0.62 0.6 0.58 0.56 0.54 0.52 1 0.5 1 0.5 0 0 −0.5 −0.5 −1 −1

Performances in FEM � � ( ∇ u · ∇ v + uv ) dx + uv ds Ω ∂ Ω P 2 elements, 7 Gauss points per triangle. T assembling ( s ) T solve ( s ) N dof 41298 2.0 0.35 81082 4.4 0.76 203727 12.3 2.0 404116 26.0 4.6 806447 52.9 11.5 4016167 304 94.3

Eigenvalue problem Example 3 � − ∆ u = λ u in Ω u = 0 on ∂ Ω % Finite elements u = fem(’P1’, mesh, meshb); v = fem(’P1’, mesh, meshb); % Mass matrix M = integral(mesh,u,v); % Rigidity matrix K = integral(mesh,grad(u),grad(v)); % Find eigenvalues [V,EV] = eigs(K,M,2*Neig,’SM’);

Result

Performances � � ∇ u · ∇ v dx = λ uv dx Ω Ω P 1 elements, 4 Gauss points per tetrahedron 20 eigenvalues/eigenvectors. N v N dof T assembling ( s ) T solve ( s ) 10000 8448 2.27 0.94 50000 44688 13.6 12.9 100000 94608 29.3 48.5

Inside the box - FEM Express the matrix as products of sparse matrices. � = A i , j φ i φ j dx Ω � ∼ ω k φ i ( x k ) φ j ( x k ) k Setting B kj = φ j ( x k ) and W kk = ω k (all sparse matrices), we get A = B t WB .

Laplace operator � A i , j = ∇ φ i · ∇ φ j dx Ω ∂φ i ∂ x ( x k ) ∂φ j � ∼ ∂ x ( x k ) ω k k ∂φ i ∂ y ( x k ) ∂φ j � + ∂ y ( x k ) ω k k ∂φ i ∂ z ( x k ) ∂φ j � + ∂ z ( x k ) ω k k Setting C x kj = ∂φ j ∂ x ( x k ) (resp. C y and C z ) and W kk = ω k , we get A = C t x WC x + C t y WC y + C t z WC z

Inside the box - FEM A = integral(mesh,opr(fem_test),func,opr(fem_unk)) Builds on the fly the integration points and weights ( x k , ω k ) Builds on the fly the sparse matrices C 1 = opr ( fem _ test )( x k ) , C 2 = opr ( fem _ unk )( x k ) Builds on the fly the diagonal sparse matrix W kk = ω k ∗ func ( x k ) Returns A = C t 1 WC 2

Generalization to BEM Galerkin formalism � � A i , j = φ i ( x ) G ( x , y ) φ j ( y ) dx dy Ω 1 Ω 2 � � ∼ φ i ( x k ) ω k G ( x k , y l ) ω l φ j ( y l ) k l Setting B 1 ki = φ i ( x k ) , B 2 lj = φ j ( y l ) , W 1 kk = ω k , W 2 ll = ω l , and G kl = G ( x k , y l ) we get A = B t 1 W 1 GW 2 B 2 Same matrices W and B as before (on possibly different meshes). Only multiplication by the full matrix G .

BEM Formalism Galerkin % G(x,y) = exp(ik|x-y|)/|x-y| Gxy = @(X,Y) femGreenKernel(X,Y,’[exp(ikr)/r]’,k); % Finite elements u = fem(’P1’); v = fem(’P1’); % int_Sx int_Sy psi(x)’ G(x,y) psi(y) dx dy LHS = 1/(4*pi) * integral(bnd,bnd,v,Gxy,u); Collocation (or radiation) LHS = 1/(4*pi) * integral(X,bnd,Gxy,v);

Single layer potential

Single layer potential Total field solution 1.5 1.6 3 1.4 2 1 1 1.2 0 0.5 1 −1 0.8 −2 0 0.6 −3 3 0.4 2 −0.5 1 0 0.2 −1 −2 0.5 −1 0 −0.5 0 1 2 3 4 5 6 7 −3

Performances in BEM � � exp ( ik | x − y | ) u ( x ) v ( y ) dx dy | x − y | Σ Σ P 1 elements, 3 Gauss points per triangle. N dof T assembling T solve 1000 18.0 1.48 1500 33.9 5.12 2000 59.4 11.4 2500 170 22 3000 (swap pb) 676 34

Interfacing with OpenHMX OpenHMX an open source H -matrix toolbox written by Matthieu Aussal. Includes the whole H -matrix algebra, including LU factorization and resolution. Object oriented (new class of matrices on which one can use classical operations such as +, *, \, spy , lu , etc. In principle the matrices that we consider are of the form A = C t 1 GC 2 where C 1 and C 2 (dof to quad points) are sparse and G is dense (kernel computed at integration points). G could be compressed using H -matrices...

This is not what we want.. G is very big. Multiplying H -matrices with sparse matrices lead to a mixed structure... Build directly the H -matrix from the finite element structure. Subdivide using the dof structure. Apply the matrices C 1 , C 2 when building the leaf structure and the ACA compression. LHS = 1/(4*pi) * integral(bnd,bnd,v,Gxy,u,tol);

G Ψ lab - Plane wave scattering by submarine 100.000 ddl at 200 Hz Neumann boundary condition Assembling H-Matrix : 1720 sec and 4.5 G0 (7 cores) LU factorisation : 900 sec Resolution : 3.5 sec Radiation : 764 sec (7 cores)

Classical CFIE for PEC materials � TJ · J ′ ds = ik � � � � J ′ ( x ) · G ( x − y ) J ( y ) dx dy − i / k div ( J ′ )( x ) · G ( x − y ) div ( J )( y ) dx dy � nxKJ · J ′ ds = − � � det ( n × J ′ ( x ) , ∇ y G ( x − y ) , J ( y )) dx dy

Validation and performances Sphere scattering, CFIE (10 operators), 30000 dof Assembling (LHS, RHS) 1000 s, solving 510 s 20 10 15 8 10 6 5 4 0 2 −5 0 −10 −2 0 2 4 6 8 0 2 4 6 8

FEM-BEM coupling In progress... Unified formalism for FEM and BEM Example of electromagnetic of PEC+dielectric materials Domain Ω = Ω PEC ∪ Ω d Current J on ∂ Ω , Electric field E tr in Ω d ( Γ d = ∂ Ω d , k d = k √ ε r µ r ) Equations: TJ · J ′ ds + 1 ( K + 1 � � � 2 n × )( E tr × n ) · J ′ ds = E inc · J ′ ds √ ε r ∂ Ω Γ d ∂ Ω � d E tr · E ′ dx − i √ µ r k d � J · E ′ ds = 0 curl ( E tr ) · curl ( E ′ ) − k 2 Ω d Γ d

Discretization and subtleties J and J ′ are RWG surfacic functions defined on the boundary ∂ Ω E tr and E ′ are taken to be Nedelec (edge) volumic functions defined on Ω d (tetrahedral mesh) Surfacic and volumic integrals + coupling done in a trace sense on Γ d The global matrix of the system has a structure � A JJ � B JE M = C EJ D EE that we wish to store as a H -matrix for LU decomposition.

Conclusion G Ψ lab is a matlab toolbox for FEM-BEM programming Easy to use interface Includes H -matrix toolbox G Ψ lab is open source (GPL like license), downloadable at www.cmap.polytechnique.fr/˜aussal/gypsilab