Unfitted Finite Element Methods

Decoupling the mesh from the geometry

Eric Neiva - Warwick Applied Maths Seminar - 25-11-2022

 

Acknowledgement of EU funding
This material is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 949267).
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Mesh-based approximation of PDEs

Continuous Discrete
L(u)=f in Ω\mathcal{L}(u) = f \text{ in } \Omega Au=fA \boldsymbol{u} = \boldsymbol{f}
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Mesh-based approximation of PDEs

Continuous Discrete
L(u)=f in Ω\mathcal{L}(u) = f \text{ in } \Omega Au=fA {\color{blue} \boldsymbol{u}} = \boldsymbol{f}
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Unfitted aka Embedded aka Immersed boundary aka etc

Body-fitted Unfitted
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Large-scale body-fitted meshing is very inefficient

Meshing bottleneck

In Hughes et al. CMAME 194 ('05) 4135–4195:

The typical situation in engineering practice is that designs are encapsulated in CAD systems and meshes are generated from CAD data. [...] It is estimated that about 80% of overall analysis time is devoted to mesh generation in the automotive, aerospace, and ship building industries. In the automotive industry, a mesh for an entire vehicle takes about four months to create.

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Large-scale body-fitted meshing is very inefficient

Partitioning bottleneck
Body-fitted > expensive Cartesian grid > trivial
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Large-scale body-fitted meshing is very inefficient

Adaptive mesh refinement bottleneck

Tree-based are currently the most scalable AMR approach:

By Offermans et al.

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Unfitted FEM can overcome all these bottlenecks

STLCutters.jl in Badia, Martorell & Verdugo, '22

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Unfitted FEM can overcome all these bottlenecks

STLCutters.jl in Badia, Martorell & Verdugo, '22
julia> include("examples/LinearElasticity.jl")
julia> filename = "test/data/550964.stl"
julia> LinearElasticity.main(filename,n=50,force=(tand(5),0,-1),output="Eiffel")

julia> include("examples/Stokes.jl")
julia> filename = "test/data/47076.stl"
julia> Stokes.main(filename,n=10,output="ChichenItza")

Try other STLs from the Thingi10K dataset!

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Unfitted FEM can overcome all these bottlenecks

Uniform grids Adaptive grids
Verdugo, Martín & Badia, '19 Badia, N., Martín & Verdugo, '21
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Outline

  • Motivation
  • The Finite Element Method
  • Unfitted FE approximations
  • Numerical examples
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Poisson problem L(u)=−Δu\mathcal{L}(u) = - \Delta u

Differential equation or strong form

{−Δu=f in Ω,u=0 on ΓD,∇u⋅n=g on ΓN. \left\lbrace \begin{aligned} - \Delta u &= f \quad \ &\text{in} \ \Omega, \\ u &= 0 \ &\text{on} \ \Gamma_{\rm D}, \\ \nabla u \cdot \boldsymbol{n} &= g \ &\text{on} \ \Gamma_{\rm N}. \end{aligned} \right.

u∈C2(Ω‾), f∈C(Ω),\footnotesize u \in \mathcal{C}^2(\overline{\Omega}), \ f \in \mathcal{C}(\Omega),
g∈C(∂Ω∖ΓD) and n∈[C(∂Ω)]d\footnotesize g \in \mathcal{C}(\partial \Omega \setminus \Gamma_{\rm D}) \ \text{and} \ \boldsymbol{n} \in [\mathcal{C}(\partial \Omega)]^d

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

FEM approximates the variational or weak formulation

Let

H1(Ω)≐{v∈L2(Ω):∇v∈L2(Ω)}, andV≐{v∈H1(Ω):v=0 on ΓD}=HΓD1(Ω). \begin{aligned} H^1(\Omega) &\doteq \{v \in L_2(\Omega) : \nabla v \in L_2(\Omega)\}, \ \text{and} \\ V &\doteq \{v \in H^1(\Omega) : v = 0 \ \text{on} \ \Gamma_D\} = H^1_{\Gamma_{\rm D}}(\Omega). \end{aligned}

Multiply by v∈Vv \in V, average over Ω\Omega and integrate by parts:

∫Ωfv dΩ=∫Ω(−Δu)v dΩ=∫Ω∇u⋅∇v dΩ −∫ΓN(∇u⋅n)v dΩ==∫Ω∇u⋅∇v dΩ −∫ΓNgv dΓ,∀ v∈V \begin{aligned} \int_{\Omega} f v \ \mathrm{d}\Omega = \int_\Omega (- \Delta u)v \ \mathrm{d}\Omega &= \int_\Omega \nabla u \cdot \nabla v \ \mathrm{d}\Omega \ - \int_{\Gamma_{\rm N}} ( \nabla u \cdot \boldsymbol{n} ) v \ \mathrm{d}\Omega = \\ &= \int_\Omega \nabla u \cdot \nabla v \ \mathrm{d}\Omega \ - \int_{\Gamma_{\rm N}} g v \ \mathrm{d}\Gamma, \quad \forall \ v \in V \end{aligned}

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

FEM approximates the variational or weak formulation

Let

H1(Ω)≐{v∈L2(Ω):∇v∈L2(Ω)}, andV≐{v∈H1(Ω):v=0 on ΓD}=HΓD1(Ω). \begin{aligned} H^1(\Omega) &\doteq \{v \in L_2(\Omega) : \nabla v \in L_2(\Omega)\}, \ \text{and} \\ V &\doteq \{v \in H^1(\Omega) : {\color{yellow}v = 0} \ \text{on} \ \Gamma_D\} = H^1_{\Gamma_{\rm D}}(\Omega). \end{aligned}

Multiply by v∈Vv \in V, average over Ω\Omega and integrate by parts:

∫Ωfv dΩ=∫Ω(−Δu)v dΩ=∫Ω∇u⋅∇v dΩ −∫ΓN(∇u⋅n)v dΩ==∫Ω∇u⋅∇v dΩ −∫ΓNgv dΓ,∀ v∈V \begin{aligned} \int_{\Omega} f v \ \mathrm{d}\Omega = \int_\Omega (- \Delta u)v \ \mathrm{d}\Omega &= \int_\Omega \nabla u \cdot \nabla v \ \mathrm{d}\Omega \ - \int_{\Gamma_{\rm N}} ( \nabla u \cdot \boldsymbol{n} ) v \ \mathrm{d}\Omega = \\ &= \int_\Omega \nabla u \cdot \nabla v \ \mathrm{d}\Omega \ - {\color{yellow}\int_{\Gamma_{\rm N}} g v \ \mathrm{d}\Gamma}, \quad \forall \ v \in V \end{aligned}

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

FEM approximates the variational or weak formulation

Statement of the weak form

Find u∈V such that a(u,v)=l(v),∀v∈V, \text{Find} \ u \in V \ \text{such that} \ a(u,v) = l(v), \enskip \forall v \in V,

where

a(u,v)≐∫Ω∇u⋅∇v dΩbilinear symmetric (continuous)l(v)≐∫Ωfv dΩ+∫ΓNgv dΓlinear (bounded) \begin{aligned} a(u,v) &\doteq \int_\Omega \nabla u \cdot \nabla v \ \mathrm{d}\Omega &\text{bilinear symmetric (continuous)}\\ l(v) &\doteq \int_{\Omega} f v \ \mathrm{d}\Omega + \int_{\Gamma_{\rm N}} g v \ \mathrm{d}\Gamma &\text{linear (bounded)} \end{aligned}

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Variational problem is related to the minimisation problem

Total potential energy

J(v)≐12a(v,v)−ℓ(v) J(v) \doteq \frac{1}{2} a(v,v) - \ell(v)

Statement of the minimisation problem

u≐argminv∈V J(v)u \doteq \underset{v \in V}{\mathrm{argmin}} \ J(v)

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Variational problem is related to the minimisation problem

Perturbing around the global min. u∈Vu \in V with v∈Vv \in V and α∈R\alpha \in \mathbb{R},

0=lim⁡α→0J(u+αv)−J(u)α=lim⁡α→0αa(u,v)+α22a(v,v)−αℓ(v)α=a(u,v)−ℓ(v),∀v∈V,\begin{align*} 0 = \lim_{\alpha \to 0} \frac{J(u + \alpha v) - J(u)}{\alpha} &= \lim_{\alpha \to 0} \frac{\alpha a(u,v) + \frac{\alpha^2}{2} a(v,v) - \alpha \ell(v)}{\alpha} \\ &= a(u,v) - \ell(v), \quad \forall v \in V, \end{align*}

we recover the variational statement.

Remark

Equivalence between strong and weak in the sense of distributions

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Galerkin approximation

Goal: Map ∞\infty-dim. variational problem to a discrete one.

Given

u∈V:a(u,v)=l(v),∀v∈V,u \in V : a(u,v) = l(v), \quad \forall v \in V,

let us consider a finite dimensional vector subspace VN⊂VV_N \subset V, with dim(VN)=N\mathrm{dim}(V_N) = N. The Galerkin approximation of the problem reads

uN∈VN:a(uN,vN)=l(vN),∀vN∈VN.u_N \in V_N : a(u_N,v_N) = l(v_N), \quad \forall v_N \in V_N.

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Galerkin approximation as a linear system

Since VNV_N is a real vector space of dim <∞< \infty, we can pick a basis B={φ1,…,φN}\mathcal{B} = \{ \varphi_1, \ldots, \varphi_N \} of VNV_N. Thus, the Galerkin problem is equivalent to:

Find uN=∑i=1Nuiφiu_N = \sum_{i=1}^N u_i \varphi_i with u=(u1,…,un)t∈RN\boldsymbol{u} = (u_1,\ldots,u_n)^t \in \mathbb{R}^N such that

Au=f, whereA∈RN×N, Aij≐a(bj,bi),i,j∈{1,…,N},f∈RN, fi≐ℓ(bi),i∈{1,…,N}. \begin{aligned} \mathbf{A}& \boldsymbol{u} = \mathbf{f}, \ \text{where} \\ \mathbf{A}& \in \mathbb{R}^{N \times N}, \ A_{ij} \doteq a(b^j,b^i), \quad &i,j \in \{ 1, \ldots, N \}, \\ \mathbf{f}& \in \mathbb{R}^{N}, \ f_i \doteq \ell(b^i), \quad &i \in \{ 1, \ldots, N \}. \end{aligned}

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

The Finite Element Method

  • Mesh
    • Th\mathcal{T}_h is partition of Ωh≈Ω\Omega_h \approx \Omega into cells KK
    • h=max⁡K∈Thdiam(K)h = \max_{K \in \mathcal{T}_h} {\rm diam}(K) (mesh size)
  • Basis of VNV_N
    • Globally continuous functions
    • Polynomials inside KK

Vh={vh∈C(Ωh):vh∣K∈P(K),∀K∈Th} V_h = \{ v_h \in \mathcal{C}(\Omega_h) : v_h |_K \in \mathcal{P}(K), \forall K \in \mathcal{T}_h \}

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Lagrangian Finite Elements

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Numerical integration of the weak form

E.g. Aij=a(φj,φi)=∫Ωh∇φj⋅∇φi dΩ,=∑K∈Th∫K∇φj∣K⋅∇φi∣K dΩK,=∑K∈Th(∑g∈Qwg∇φj∣K(xg)⋅∇φi∣K(xg)) \footnotesize{} \begin{aligned} \text{E.g.} \ A_{ij} = a(\varphi_j,\varphi_i) &= \int_{\Omega_h} \nabla \varphi_j \cdot \nabla \varphi_i \ \mathrm{d}\Omega, \\ &= \sum_{K \in \mathcal{T}_h} \int_K \nabla \varphi_j |_K \cdot \nabla \varphi_i |_K \ \mathrm{d}{\Omega_K}, \\ &= \sum_{K \in \mathcal{T}_h} \left( \sum_{g \in Q} w_g \nabla \varphi_j |_K (x_g) \cdot \nabla \varphi_i |_K (x_g) \right) \end{aligned}

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Unfitted finite element methods

Body-fitted Unfitted
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Unfitted finite element methods

  • Implicit geometry representation is flexible

    • A level-set function ψ\psi such that ∂Ω≡{ψ=0}\partial \Omega \equiv \{ \psi = 0 \}
    • An STL file, scan data... the key is to tell inside from outside Ω\Omega

  • Challenges

    • Cut-cell integration
    • Imposing Dirichlet (essential) boundary conditions
    • Stability and conditioning w.r.t. small cut elements

  • See de Prenter, Verhoosel, van Brummelen, Larson & Badia, subm.

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Cut-cell integration

  • Now Th\mathcal{T}_h triangulates Ω~\tilde{\Omega} and Ω⊂Ω~\Omega \subset \tilde{\Omega}, so

E.g. Aij=a(φj,φi)=∑K∈Th∫Ω∩K∇φj∣K⋅∇φi∣K dΩK \text{E.g.} \ A_{ij} = a(\varphi_j,\varphi_i) = \sum_{K \in \mathcal{T}_h} \int_{\color{yellow} \Omega \cap K} \nabla \varphi_j |_K \cdot \nabla \varphi_i |_K \ \mathrm{d}{\Omega_K}

  • Must modify local integration on cut cells
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Cut-cell integration uses implicit geometry representation

  • See review in, e.g., Saye, '22
  • Techniques w/different accuracy/robustness/cost balances, e.g.:
Subtriangulation Height functions
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Imposing Dirichlet boundary conditions

Model problem Body-fitted Unfitted
Close-up Remove Dirichlet DoFs ???
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Imposing Dirichlet boundary conditions

aK(u,v)≐∫Ω∩K∇u⋅∇v dΩ+∫ΓD∩KβKuv−(∇u⋅n)v−(∇v⋅n)u dΩlK(v)≐∫Ω∩Kfv dΩ+∫ΓN∩Kgv dΓ+∫ΓD∩KβKuDv−(∇v⋅n)uD dΩ \begin{aligned} a_K(u,v) &\doteq \int_{\Omega \cap K} \nabla u \cdot \nabla v \ \mathrm{d}\Omega \color{yellow}+ \int_{\Gamma_{\rm D} \cap K} \beta_K uv - (\nabla u \cdot \boldsymbol{n} ) v - (\nabla v \cdot \boldsymbol{n} ) u \ \mathrm{d}\Omega \\ l_K(v) &\doteq \int_{\Omega \cap K} f v \ \mathrm{d}\Omega + \int_{\Gamma_{\rm N} \cap K} g v \ \mathrm{d}\Gamma \color{yellow}+ \int_{\Gamma_{\rm D} \cap K} \beta_K u_{\rm D}v - (\nabla v \cdot \boldsymbol{n} ) u_{\rm D} \ \mathrm{d}\Omega \end{aligned}

  • βK\color{yellow}\beta_K must be large enough to ensure coercivity (uniqueness)
    (In our model problem uD=0u_D = 0)
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Stability to small cut cells

  • If nothing done in unfitted, βK\beta_K can be arbitrarily large:
Body-fitted Naive unfitted
βK∼hK−1\beta_K \sim h_{K}^{-1} βK∼hΩ∩K−1\beta_K \sim h_{\Omega \cap K}^{-1}
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Ill-conditioning with small cut cells

  • If nothing done in unfitted, k2(A)k_2(A) can be arbitrarily large:
Body-fitted Naive unfitted
k2(A)∼h−2k_2(A) \sim h^{-2} k2(A)∼η−(2p+1−2/d)k_2(A) \sim \eta^{-(2p+1-2/d)}
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022
Stable and well-conditioned unfitted FEMs

Ghost penalty

Add stabilising term. Linear case:

s(u,v)=∑F∈Fg∫Fγsh ⟦nF⋅∇u⟧⟦nF⋅∇v⟧ s(u,v) = \sum_{F \in {\color{orange}\mathcal{F_g}}} \int_{F} \gamma^s h \ \llbracket \boldsymbol{n}_F \cdot \nabla u \rrbracket \llbracket \boldsymbol{n}_F \cdot \nabla v \rrbracket

where Fg{\color{orange}\mathcal{F_g}} is the ghost skeleton.

Order > 1: add high-order derivatives.

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022
Stable and well-conditioned unfitted FEMs

Aggregated FEM (AgFEM)

∙{\color{blue}\bullet} interior DOFs (I\mathcal{I}) and ×{\color{red}\boldsymbol{\times}} exterior DOFs (E\mathcal{E})

Extrapolate ×{\color{red}\boldsymbol{\times}} DOFs from ∙{\color{blue}\bullet} DOFs:

u×=∑∙∈K~(×)φ∙(x×)u∙,∀×∈E u_{{\color{red}\boldsymbol{\times}}} = \sum_{{\color{blue}\bullet}\in \tilde{K}({\color{red}\boldsymbol{\times}})} \varphi_{{\color{blue}\bullet}}(x_{\color{red}\boldsymbol{\times}})u_{{\color{blue}\bullet}}, \quad \forall {\color{red}\boldsymbol{\times}} \in \mathcal{E}

where K~(×)\tilde{K}({\color{red}\boldsymbol{\times}}) is an interior cell of Th\mathcal{T}_h.

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Stable and well-conditioned unfitted FEMs

β∼h−1\beta \sim h^{-1} and k2(A)∼h−2k_2(A) \sim h^{-2} holds for both methods ✅

Refer to

for

  • Mathematical analysis and numerical comparison
  • Links between AgFEM and ghost penalty
 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

AgFEM examples: (Navier-)Stokes flows

Badia, Verdugo and Martín, SISC 40(6) ('18), B1541-B1576

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

AgFEM examples: (Multi-)interface flows

N. and Badia, CMAME 380 ('21), 113769

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

AgFEM examples: High-order approximations

Badia, N. and Verdugo, Comput. Math. with Appl. 127 ('22), 105-126

 
Eric Neiva | CDF-CNRS | Unfitted FEM | 25-11-2022

Try unfitted FEM yourself at GridapEmbedded.jl

  • A very expressive API to solve complex PDEs with few lines of code
  • Write weak form with almost 1:1 to the mathematical notation
 

Thank you!!!


Special thanks to F. Verdugo and P.A. Martorell for providing some of their pictures for these slides.

This material is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 949267).

Except where otherwise noted, this work and its contents (texts and illustrations) are licensed under the Attribution 4.0 International (CC BY 4.0).

Maybe recall that mesh distributed among processors

- Bypassing functional analysis and existence and uniqueness theorems - The keys are that FEM approximates the variational problem and its relation to the minimisation problem

Go fast. Natural vs essential BCs

Go fast. Natural vs essential BCs

Comment on regularity

- Recasting the problem into a discrete vector subspace - If (V) is well-posed, then (G) is well-posed (previous theorem for finite dim, also)

The two key ingredients are:

- We exchanged geometric complexity for integration complexity - See [de Prenter, Verhoosel, van Brummelen, Larson & Badia, subm.](https://arxiv.org/abs/2208.08538)

| [Burman et al., 2015](https://onlinelibrary.wiley.com/doi/full/10.1002/nme.4823) | [Saye, 2015](https://epubs.siam.org/doi/10.1137/140966290) and [2022](https://www.sciencedirect.com/science/article/pii/S002199912100615X) |

Recasting as the graph of a height function