{ "cells": [ { "cell_type": "markdown", "source": [ "# Tutorial 19: Unfitted Poisson" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "In this **tutorial**, we will learn\n", "\n", " - How to use [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl)\n", " - How to formulate an unfitted finite element method\n", " - How to construct an unfitted boundary geometry setup\n", " - How to use the Aggregated Finite Element Method\n", " - How to impose Dirichlet boundary conditions _weakly_\n", "\n", "## Problem statement\n", "\n", "We want to solve the **Poisson equation** on the crescent moon 2D shape shown in the next figure\n", "\n", "![fig0a](../assets/unfitted_poisson/crescent_moon.png)\n", "\n", "We prescribe the linear solution to the problem $u(x) = x - y$ and Dirichlet boundary conditions on the whole boundary. Thus, the problem to solve is: find the scalar field $u$ such that\n", "\n", "$$\n", "\\left\\lbrace\n", "\\begin{aligned}\n", "-\\Delta u &= 0 \\ &\\text{in} \\ \\Omega,\\\\\n", "u &= x - y \\ &\\text{on}\\ \\partial\\Omega\n", "\\end{aligned}\n", "\\right.\n", "$$\n", "\n", "## Numerical scheme\n", "\n", "To solve this PDE, we will use an unfitted technique. The idea behind _unfitted_, aka _embedded_ or _immersed boundary_, methods is quite simple: Instead of relying on a mesh that fits to the domain boundary, we embed the domain in a background grid, upon which we will define our finite elements.\n", "\n", "| Body-fitted | Unfitted |\n", "| ---- | ---- |\n", "| ![fig0b](../assets/unfitted_poisson/fig_2_body-fitted_a_bis.png) | ![fig0c](../assets/unfitted_poisson/fig_3_unfitted.png)\n", "\n", "Unfitted methods allow one to circumvent the meshing bottleneck when dealing with very complex geometries. But, as usual, this comes at a price. In unfitted methods, integration of the weak form and imposition of Dirichlet boundary conditions becomes more involved. In addition, we have to be careful with the small cut cell problem, which refers to stability and ill-conditioning issues that appear in presence of very small intersections between the background cells and the domain.\n", "\n", "Fortunately, there are plenty of solutions available to take care of all these challenges. For more details, we refer to this recent review:\n", "\n", "> F. de Prenter, C. Verhoosel, H. van Brummelen, M. Larson, S. Badia, _Stability and conditioning of immersed finite element methods: analysis and remedies_ arXiv preprint [arXiv:2208.08538](https://arxiv.org/pdf/2208.08538.pdf)\n", "\n", "In this tutorial, we will use the functionality available at [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl) to formulate and solve the unfitted problem above. We will illustrate a typical unfitted FE simulation pipeline in [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl), including how we setup the embedded geometry, how we deal with the integration of the weak form and the imposition of Dirichlet boundary conditions.\n", "\n", "We will discretise the problem with the aggregated finite element method, which is robust to the small cut cell problem. For more details, we refer to:\n", "> S. Badia, A.F. Martín, F. Verdugo, _[The aggregated unfitted finite element method for elliptic problems](https://www.sciencedirect.com/science/article/pii/S0045782518301476)_, Computer Methods in Applied Mechanics and Engineering 336 (2018), 533-553.\n", "\n", "We recommend to read both references above to grasp the details of this tutorial.\n", "\n", "The first step is to load the [Gridap](https://github.com/gridap/Gridap.jl) and [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl) packages." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "using Gridap\n", "using GridapEmbedded" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "### Implicit boundary representation\n", "\n", "The next step is to generate the embedded geometry. The most straightforward way to represent the geometry is with the level-set method. It amounts to represent a closed curve (2D) or surface (3D) $\\Gamma$ using an auxiliary function $\\psi$, called the level-set function. $\\Gamma$ is represented as the zero-level-set of $\\psi$ by\n", "\n", "$$\n", "\\Gamma = \\{ (x,y) \\ | \\ \\psi(x,y) = 0 \\}\n", "$$\n", "\n", "and the level-set method manipulates $\\Gamma$ _implicitly_, through the function $\\psi$. This function $\\psi$ is assumed to take negative values inside the region delimited by $\\Gamma$ and positive values outside.\n", "\n", "[GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl) offers out-of-the-box level set surface descriptions for the `disk`, `cylinder`, `sphere`, `square`, `tube`, among others. You can also implement your own level-set function.\n", "\n", "Interestingly, you can also construct complex domains from simple level-set descriptions. This is done by translating Boolean set operations and geometric transformations into simple manipulations of the level-set functions. For instance, given two level-set functions $\\psi_1$ and $\\psi_2$, representing the domains $\\Omega_1$ and $\\Omega_2$, the level-set function $\\psi = \\min(\\psi_1,\\psi_2)$ gives the domain $\\Omega_1 \\cup \\Omega_2$.\n", "\n", "We will use the former properties to generate the crescent moon shape with [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl) as follows. First, we note that the crescent moon shape is the outcome of subtracting to a disk of radius $R=1/2$ another disk of radius R, whose center is at a relative position $(-R/2,R/2)$.\n", "\n", "![fig1](../assets/unfitted_poisson/crescent_moon_setup.png)\n", "\n", "Defining this geometry in [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl) is done as follows:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "R = 0.5\n", "L = 0.5*R\n", "p1 = Point(0.0,0.0)\n", "p2 = p1 + VectorValue(-L,L)\n", "\n", "geo1 = disk(R,x0=p1)\n", "geo2 = disk(R,x0=p2)\n", "geo3 = setdiff(geo1,geo2)" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "Where we have used the `setdiff` function to subtract the disk `geo2` to the disk `geo1`. This is a very simple example of constructive solid geometry. [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl) allows to create more complex CSG geometries as shown [here](https://github.com/gridap/GridapEmbedded.jl#constructive-solid-geometry-csg). For more general shapes, [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl) admits an STL description of the implicit geometry, see package [STLCutters](https://github.com/gridap/STLCutters.jl) for more details.\n", "\n", "### Unfitted triangulations\n", "\n", "As explained, the rationale of unfitted methods is to embed the domain of interest in a background Cartesian grid (or, more generally, any easy-to-generate mesh). In our case, we generate a 30x30 Cartesian grid model containing the geometry as" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "t = 1.01\n", "pmin = p1-t*R\n", "pmax = p1+t*R\n", "\n", "n = 30\n", "partition = (n,n)\n", "bgmodel = CartesianDiscreteModel(pmin,pmax,partition)\n", "dp = pmax - pmin" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "Unfitted FE formulations in [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl) hang on two different type of triangulations: \"Active\" and \"Physical\". In order to generate them, we need first to cut the embedded geometry against the model with the function `cut`." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "cutgeo = cut(bgmodel,geo3)" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "`cut` generates an `EmbeddedDiscretization` instance (here, `cutgeo`). An `EmbeddedDiscretization` classifies all cells and facets from the background model into inside, outside or cutting the embedded geometry. Cells and facets of the background model outside the embedded geometry play no role in the unfitted formulation, they are _inactive_, thus it is convenient to remove them and restrict the triangulations to the _active_ portion of the model (i.e., cut and interior cells and facets). In [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl) we can do this using the classification from `cutgeo` and the `ACTIVE` keyword." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "Ω_act = Triangulation(cutgeo,ACTIVE)" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "To illustrate this concept, we can plot both the background and active triangulations to compare them." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "Ω_bg = Triangulation(bgmodel)\n", "writevtk(Ω_bg,\"bg_trian\")\n", "writevtk(Ω_act,\"act_trian\")" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "In the picture below of the background grid, white cells are _inactive_, whereas gray cells are _active_.\n", "\n", "![fig2](../assets/unfitted_poisson/fig_active_triangulation.png)\n", "\n", "As we will see later, we define our unfitted FE spaces on _active_ triangulations.\n", "\n", "An `EmbeddedDiscretization` instance (here, `cutgeo`) also generates subtriangulations on each cut cells to represent the portion of the cell which is inside the domain of analysis. We use these subtriangulations to generate the so called _physical_ triangulations. Physical triangulations are nothing other than a body-fitted mesh of our domain $\\Omega$, but _we only use them to integrate the weak form_ of the problem in $\\Omega$, we won't define FE spaces and assign DoFs on top of them. In [GridapEmbedded](https://github.com/gridap/GridapEmbedded.jl) we build physical triangulations using the `PHYSICAL` keyword." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "Ω = Triangulation(cutgeo,PHYSICAL)\n", "writevtk(Ω,\"phys_trian\")" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "Once again, we can combine plots of the physical and active triangulations to illustrate these concepts. In the first plot, we show the physical triangulation within the background one.\n", "\n", "![fig3](../assets/unfitted_poisson/fig_physical_trian_1.png)\n", "\n", "Note that cut cells are subtriangulated to approximate the embedded geometry, as exposed in the following close-up plot.\n", "\n", "![fig4](../assets/unfitted_poisson/fig_physical_trian_2.png)\n", "\n", "In the third plot, we show the region represented by physical triangulation (shaded in gray) embedded in the active grid (black-contoured cells).\n", "\n", "![fig5](../assets/unfitted_poisson/fig_physical_trian_3.png\" width=\"420\")\n", "\n", "**In a nutshell,** to define the unfitted FE formulation of the problem we need the \\\"active\\\" and \\\"physical\\\" triangulations of the domain. The former triangulation is used to define the FE spaces, whereas the latter is used to integrate the weak form. We use a level-set function to derive both of them.\n", "\n", "| Active | Physical |\n", "| ---- | ---- |\n", "| _For_ FE spaces | _For_ measures |\n", "| ![fig6](../assets/unfitted_poisson/fig_active_triangulation.png) | ![fig7](../assets/unfitted_poisson/fig_physical_trian_1.png)\n", "\n", "### Unfitted FE spaces\n", "\n", "Standard FE methods, i.e. the usual ones for body-fitted meshes, are exposed to the small cut cell problem, which leads to stability and conditioning issues due to the presence of arbitrarily small cut cells. In unfitted FEM, it is impractical to have control over how the mesh intersects the geometry, so we need to treat these numerical issues.\n", "\n", "As said, in this tutorial, we leverage the aggregated unfitted finite element method (AgFEM), see the figure below. The main idea is to remove DoFs on cut cells, whose basis functions potentially have very small local support. In order to do that, we constrain exterior DoFs ${\\color{red}\\times}$ in terms of interior DOFs ${\\color{blue}\\bullet}$ (see figure below) as follows:\n", "1. We map every exterior DoF ${\\color{red}\\times}$ to an interior cell $\\tilde{K}({\\color{red}\\times})$ of the active triangulation using a cell aggregation scheme (thus the name of the method).\n", "2. We extrapolate the value at the exterior DoF ${\\color{red}\\times}$ with the local FE basis at the interior cell $\\tilde{K}({\\color{red}\\times})$. This leads to:\n", "\n", "$$\n", " u_{{\\color{red}\\times}} = \\sum_{{\\color{blue}\\bullet}\\in \\tilde{K}({\\color{red}\\times})} \\varphi_{{\\color{blue}\\bullet}}(x_{{\\color{red}\\times}})u_{{\\color{blue}\\bullet}}, \\quad \\forall {\\color{red}\\times}\n", "$$\n", "\n", "Hence, we end up extrapolating exterior DoFs ${\\color{red}\\times}$, in terms of interior ones ${\\color{blue}\\bullet}$.\n", "\n", "![fig8](../assets/unfitted_poisson/fig_agfemspace.png)\n", "\n", "We define the AgFEM space of our problem as follows. First, we generate a standard linear FE space **on the `ACTIVE` triangulation**. We do not prescribe Dirichlet boundaries on the standard FE space for reasons that will be clear later in the tutorial." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "order = 1\n", "reffe = ReferenceFE(lagrangian,Float64,order)\n", "Vstd = TestFESpace(Ω_act,reffe,conformity=:H1)" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "After this, we construct the aggregates which will be used to build the map of exterior DoFs to interior ones." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "strategy = AggregateAllCutCells()\n", "aggregates = aggregate(strategy,cutgeo)" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "Note the following remark:\n", "\n", "> Here, we aggregate all cut cells. Alternatively, we can leave out large enough cut cells using `AggregateCutCellsByThreshold(t)`, where `t` is a real number. Given a cut cell, if the ratio of the cut region volume and the background cell volume is larger than `t`, then this cell is _not_ aggregated.\n", "\n", "We can color the aggregates and easily inspect them on the background mesh." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "colors = color_aggregates(aggregates,bgmodel)\n", "Ω_bg = Triangulation(bgmodel)\n", "writevtk(Ω_bg,\"aggs_on_bg_trian\",celldata=[\"aggregate\"=>aggregates,\"color\"=>colors])" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "Finally, we use the aggregates to constrain the exterior DoFs of `Vstd` in terms of the interior ones. This leads to the AgFEM space." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "V = AgFEMSpace(Vstd,aggregates)\n", "U = TrialFESpace(V)" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "Note that `V` is a test FE space and the trial FE space `U` does not receive Dirichlet functions, because we have not prescribed Dirichlet boundaries in the standard test FE space `Vstd`.\n", "\n", "### Unfitted measures\n", "\n", "We define next the integration measures **on the `PHYSICAL` triangulation**." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "degree = 2*order\n", "dΩ = Measure(Ω,degree)" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "Although `degree = 2*order` is enough for this particular linear case. In general, we need to enforce `degree = 2*dim*order` on cut cells. This is a detail. We just leave a brief explanation of the reason below.\n", "\n", "> Since interior cells are rectangular, Fubini theorem applies and we can separate the integrals of the weak form by dimensions. Thus, we can take as quadrature degree `2*order` as usual. However, in the case of cut cells, Fubini theorem does not hold. So, to take into account that we cannot separate the integrals, the quadrature degree on cut cells must be premultiplied by the problem dimension, i.e. `degree = 2*dim*order`.\n", "\n", "### Imposing Dirichlet boundary conditions\n", "\n", "In unfitted FEMs, we cannot impose Dirichlet boundary contions on the Poisson problem as in body-fitted ones. The reason is apparent in the figure below.\n", "\n", "We assume a Poisson problem defined on a 2D shape with Dirichlet $\\Gamma_D$ and Neumann $\\Gamma_N$ BCs and a close-up on the Dirichlet boundary.\n", "\n", "| Model problem | Body-fitted | Unfitted |\n", "| ---- | ---- | ---- |\n", "| ![fig9](../assets/unfitted_poisson/fig_cp_bcs.png) | ![fig10](../assets/unfitted_poisson/fig_bf_bcs.png) | ![fig11](../assets/unfitted_poisson/fig_ud_bcs.png) |\n", "\n", "Clearly, on body-fitted meshes, we can impose Dirichlet BCs _strongly_, i.e. remove the basis functions associated to Dirichlet DoFs from the FE space. On unfitted meshes, we must enforce these conditions weakly. In order to do that we resort to Nitsche's method:\n", "\n", "$$\n", "\\begin{aligned}\n", "a_K(u,v) &\\doteq \\int_{\\Omega \\cap K} \\nabla u \\cdot \\nabla v \\ \\mathrm{d}\\Omega + \\int_{\\Gamma_{\\rm D} \\cap K} \\beta_K uv - (\\nabla u \\cdot \\boldsymbol{n} ) v - (\\nabla v \\cdot \\boldsymbol{n} ) u \\ \\mathrm{d}\\Omega \\\\\n", "l_K(v) &\\doteq \\int_{\\Omega \\cap K} f v \\ \\mathrm{d}\\Omega + \\int_{\\Gamma_{\\rm N} \\cap K} g v \\ \\mathrm{d}\\Gamma + \\int_{\\Gamma_{\\rm D} \\cap K} \\beta_K u_{\\rm D}v - (\\nabla v \\cdot \\boldsymbol{n} ) u_{\\rm D} \\ \\mathrm{d}\\Omega\n", "\\end{aligned}\n", "$$\n", "\n", "where $\\beta_K$ must be large enough to ensure coercivity of the problem.\n", "\n", "This means we need to integrate boundary tems along the embedded boundary. For this purpose we generate the embedded boundary triangulation and measure." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "Γ = EmbeddedBoundary(cutgeo)\n", "n_Γ = get_normal_vector(Γ)\n", "dΓ = Measure(Γ,degree)" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "We can now write the weak form of the problem." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "u(x) = x[1] - x[2] # Solution of the problem\n", "const γd = 10.0 # Nitsche coefficient\n", "const h = dp[1]/n # Mesh size according to the parameters of the background grid\n", "\n", "a(u,v) =\n", " ∫( ∇(v)⋅∇(u) )dΩ +\n", " ∫( (γd/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u )dΓ\n", "\n", "l(v) = ∫( (γd/h)*v*u - (n_Γ⋅∇(v))*u )dΓ" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "Once we have done all this, the rest of the pipeline follows the steps of a standard FE simulation." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "op = AffineFEOperator(a,l,U,V)\n", "uh = solve(op)\n", "\n", "e = u - uh\n", "\n", "l2(u) = sqrt(sum( ∫( u*u )*dΩ ))\n", "h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ ))\n", "\n", "el2 = l2(e)\n", "eh1 = h1(e)\n", "ul2 = l2(uh)\n", "uh1 = h1(uh)\n", "\n", "using Test\n", "@test el2/ul2 < 1.e-8\n", "@test eh1/uh1 < 1.e-7\n", "\n", "writevtk(Ω,\"results.vtu\",cellfields=[\"uh\"=>uh])" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "![fig10](../assets/unfitted_poisson/fig_results.png)" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "**Tutorial done!**" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "---\n", "\n", "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*" ], "metadata": {} } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.7" }, "kernelspec": { "name": "julia-1.6", "display_name": "Julia 1.6.7", "language": "julia" } }, "nbformat": 4 }