--- jupytext: main_language: python text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.16.7 --- # Divergence theorem This demo is implemented in {download}`demo_div_thm.py`. It illustrates: - The basic concepts of unfitted domain discretizations. - How to create an {py:class}`unfitted implicit domain ` described by an {py:class}`implicit function `. - How to compute integrals inside the unfitted domain and on its boundary using [FEniCSx](https://fenicsproject.org). ## Theorem and problem definition Given $\Omega \subset \mathbb{R}^n$ an open bounded domain with boundary $\partial \Omega = \Gamma$, and $\mathbf{F}:\bar\Omega\to\mathbb{R}^n$ a continuously differentiable vector field, where $\bar\Omega$ is the closure of $\Omega$, the divergence theorem states: $$ \begin{align} \int_{\Omega} \text{div} \textbf{F} \text{d} x = \int_{\Gamma} \textbf{F} \cdot \textbf{n} {\rm d} s \end{align} $$ where $\textbf{n}$ is the outward pointing unit normal at almost each point on the boundary $\Gamma$. Check the [wikipedia entry](https://en.wikipedia.org/wiki/Divergence_theorem) for further details. Let now consider that the domain $\Omega$ is embedded into a larger domain $\Omega^\ast\supset\Omega$ that is discretized with a mesh $\mathcal{T}(\Omega^\ast)$. The aim of this demo is to verify the divergence theorem by computing both integrals, and comparing them, considering $\Omega$ (and $\Gamma$) to be immersed in the mesh $\mathcal{T}(\Omega^\ast)$. For that purpose, we will use the quadrature rules generated by QUGaR in combination with FEniCSx capabilities. Hereinafter, and for the sake of simplicity, we will consider that the mesh $\mathcal{T}(\Omega^\ast)$ corresponds exactly to the domain $\Omega^\ast$. This is generally true, as $\Omega^\ast$ is usually chosen to be a $d$-dimensional hypercube that can be easily discretized with a Cartesian mesh. We also drop the depencency on $\Omega^\ast$ in the notation. The mesh $\mathcal{T}$ is naturally split into three disjoint families of cells as $\mathcal{T}=\mathcal{T}_{\text{cut}}\cup\mathcal{T}_{\text{full}}\cup\mathcal{T}_{\text{empty}}$, where $\mathcal{T}_{\text{cut}}$, $\mathcal{T}_{\text{full}}$, and $\mathcal{T}_{\text{empty}}$ are the submeshes containing the cut, full, and empty cells, respectively, and are defined as - $\mathcal{T}_{\text{cut}} =\left\lbrace\tau\in\mathcal{T}: \tau\cap\Omega\neq\emptyset\,\text{and}\,\tau\cap\Omega\neq\tau\right\rbrace$ - $\mathcal{T}_{\text{full}} =\left\lbrace\tau\in\mathcal{T}: \tau\cap\Omega=\tau\right\rbrace$ - $\mathcal{T}_{\text{empty}} =\left\lbrace\tau\in\mathcal{T}: \tau\cap\Omega=\emptyset\right\rbrace$ Therefore, it is easy to realize that $\mathcal{T}_{\text{full}}\subseteq\Omega$, $\mathcal{T}_{\text{cut}}\cap\Omega\neq\emptyset$, and $\mathcal{T}_{\text{empty}}\not\subset\Omega$. Thus, the volumetric integral over $\Omega$ will be approximated as the integral over $\mathcal{T}_{\text{full}}$ and $\mathcal{T}_{\text{cut}}$ as $$ \begin{align} \int_{\Omega} \text{div} \textbf{F} \text{d} x \approx \sum_{\tau\in\mathcal{T}_{\scriptsize\mbox{full}}}\int_{\tau} \text{div} \textbf{F} \text{d} x + \sum_{\tau\in\mathcal{T}_{\scriptsize\mbox{cut}}}\int_{\tau\cap\Omega} {\rm div} \textbf{F} \text{d} x \end{align} $$ The integrals over the full elements in $\mathcal{T}_{\text{full}}$ are computed using DOLFINx standard capabilities. However, the integrals over the active part of the cut elements ($\tau\cap\Omega,\,\forall\tau\in\mathcal{T}_{\text{cut}}$) will be computed using QUGaR generated quadratures in combination with new custom forms that allow DOLFINx to use quadrature rules defined at runtime. Regarding the right-hand-side of the divergence theorem, for computing the surface integral we must consider that the boundary $\Gamma$ is composed of two parts as $\Gamma=\Gamma_{\text{ext}}\cup\Gamma_{\text{int}}$, where $\Gamma_{\text{ext}}=\Gamma\cap\partial\Omega^\ast$ and $\Gamma_{\text{int}}=\Gamma\setminus\partial\Omega^\ast$ The boundary $\Gamma_{\text{ext}}$ is discretized through a subset of external facets (those that belong to a single cell) of the mesh $\mathcal{T}$ (that we denote as $\mathcal{F}$). Let's introduce $\mathcal{F}_{\text{ext}}\subset\mathcal{F}$, the subset of external facets of $\mathcal{T}$ that intersects $\Gamma_{\text{ext}}$, i.e., $\mathcal{F}_{\text{ext}} = \left\lbrace\zeta\in\mathcal{F}: \zeta\cap\Gamma\neq\emptyset\right\rbrace$. Then, $\mathcal{F}_{\text{ext}}$ can be split into two disjoint families of facets as $\mathcal{F}_{\text{ext}}=\mathcal{F}_{\text{cut}}\cup\mathcal{F}_{\text{full}}$ where where - $\mathcal{F}_{\text{cut}} =\left\lbrace\zeta\in\mathcal{F}: \zeta\cap\Gamma\neq\emptyset\,\text{and}\,\zeta\cap\Gamma\neq\zeta\right\rbrace$ - $\mathcal{F}_{\text{full}} =\left\lbrace\zeta\in\mathcal{F}: \zeta\cap\Gamma=\zeta\right\rbrace$ For computing the integral over $\Gamma_{\text{int}}$, the part of the $\Gamma$ that is immersed in the mesh $\mathcal{T}$ but does not intersect $\mathcal{F}$, we will special quadrature rules defined over the cut cells $\mathcal{T}_{\text{cut}}$. Thus, the surface integral of the divergence theorem will be approximated as $$ \begin{align} \int_{\Gamma} \textbf{F} \cdot \textbf{n} {\rm d} s \approx \sum_{\zeta\in\mathcal{F}_{\scriptsize\mbox{full}}}\int_{\zeta} \textbf{F} \cdot \textbf{n} {\rm d} s + \sum_{\zeta\in\mathcal{F}_{\scriptsize\mbox{cut}}}\int_{\zeta\cap\mathcal{F}} \textbf{F} \cdot \textbf{n} {\rm d} s + \sum_{\tau\in\mathcal{T}_{\scriptsize\mbox{cut}}}\int_{\tau\cap\Gamma_{\text{int}}} \textbf{F} \cdot \textbf{n} {\rm d} s \end{align} $$ In this demo, we will use QUGaR (and FEniCSx) to compute all the integrals defined above. For that purpose, we consider the following setup: - $\Omega^\ast = [0,1]^d$ (with $d=2,3$) (a square in 2D, a cube in 3D) - $\mathcal{T}$ is created with $16\times16$ quadrangles in 2D, and and $16\times16\times16$ hexahedra in 3D. - $\Omega=\left\lbrace x\in\Omega^\ast: \lVert x-x_0\rVert < R\right\rbrace$ (a disk in 2D, a sphere in 3D), with - $x_0 = (-0.25, 0.15)$ in 2D, and $x_0 = (0.25, 0.35, 0.45)$ in 3D - $R=0.6$ in 2D, and $R=0.7$ in 3D - $\mathbf{F} = \left(\sin(x),\cos(y)\right)$ in 2D, and $\mathbf{F} = \left(\sin(x),\cos(y),\sin(z)\right)$ in 3D. ## Implementation First we check that FEniCSx installation is available. ```python import qugar.utils if not qugar.utils.has_FEniCSx: raise ValueError("FEniCSx installation not found is required.") ``` Then the modules and functions to be used are imported: ```python from typing import cast from mpi4py import MPI import dolfinx.fem import numpy as np import ufl import qugar import qugar.cpp import qugar.impl from qugar.dolfinx import CustomForm, dx_bdry_unf, form_custom, mapped_normal from qugar.mesh import create_Cartesian_mesh from qugar.quad import create_quadrature_generator ``` We define the floating point type to be used in the computations. So far, QUGaR supports 32 and 64 bits real floating point types. In the near future support for 64 and 128 bits complex types will be also available. ```python dtype = np.float64 ``` We first define the geometry to be considered: a {py:class}`disk` in 2D and a {py:class}`sphere` in 3D. Uncomment the one you want to use. ```python # radius_2D = dtype(0.6) # center_2D = np.array([-0.25, 0.15], dtype=dtype) # impl_func = qugar.impl.create_disk(center=center_2D, radius=radius_2D) # name = "disk" radius_3D = dtype(0.7) center_3D = np.array([0.25, 0.35, 0.45], dtype=dtype) impl_func = qugar.impl.create_sphere(center=center_3D, radius=radius_3D) name = "sphere" print(f"{name.capitalize()} - divergence theorem verification") ``` Other possible geometries can be check in the [Implicit functions demo](demo_impl_funcs.md). +++ We create a {py:class}`Cartesian mesh` (corresponding to $\mathcal{T}$) in which we will embbed the domain $\Omega$, ```python dim = impl_func.dim n_cells = [16] * dim # xmin and xmax define the domain $\Omega^\ast$ xmin = np.zeros(dim, dtype) xmax = np.ones(dim, dtype) cart_mesh = create_Cartesian_mesh( MPI.COMM_WORLD, n_cells, xmin, xmax, ) dlf_mesh = cart_mesh.dolfinx_mesh ``` and then we create an unfitted domain that stores the partition $\mathcal{T}=\mathcal{T}_{\text{cut}}\cup\mathcal{T}_{\text{full}}\cup\mathcal{T}_{\text{empty}}$. ```python unf_domain = qugar.impl.create_unfitted_impl_domain(impl_func, cart_mesh) ``` We create the vector function $\mathbf{F}$ and its divergence. ```python x = ufl.SpatialCoordinate(dlf_mesh) if dim == 2: F = ufl.as_vector([ufl.sin(x[0]), ufl.cos(x[1])]) # type: ignore else: F = ufl.as_vector([ufl.sin(x[0]), ufl.cos(x[1]), ufl.sin(x[2])]) # type: ignore div_F = ufl.div(F) ``` In order to be able to integrate over the different families of cells and facets, we create tags for the cut and full cells and facets using QUGaR built-in functions. ```python cell_tags_vol = unf_domain.create_cell_tags(cut_tag=0, full_tag=0) cell_tags_int_srf = unf_domain.create_cell_tags(cut_tag=0) facet_tags = unf_domain.create_facet_tags(cut_tag=0, full_tag=0) ``` Note that the same tag can be used for two different families of cells or facets. +++ For computing volumetric integrals, we use standard UFL measures with the defined cell tags. ```python dx = ufl.dx( subdomain_id=0, domain=dlf_mesh, subdomain_data=cell_tags_vol, ) ufl_form_vol = div_F * dx ``` While in the case of the surface integrals, the procedure is twofold: for $\Gamma_{\text{int}}$ we use the {py:class}`dx_bdry_unf` measure introduced in QUGaR (note that we integrate over the cut cells only) ```python ds_int = dx_bdry_unf(subdomain_id=0, domain=dlf_mesh, subdomain_data=cell_tags_int_srf) ``` and the standard UFL external facet measure for (both parts of) $\Gamma_{\text{ext}}$. ```python ds = ufl.ds(subdomain_id=0, domain=dlf_mesh, subdomain_data=facet_tags) ``` In the same way, the unit outward normal vector at the boundary $\Gamma_{\text{int}}$ requires the {py:class}`mapped_normal` function introduced in QUGaR. ```python bound_normal = mapped_normal(dlf_mesh) ``` while the normal at the external facets is defined using the standard UFL function. ```python facet_normal = ufl.FacetNormal(dlf_mesh) ``` Using these measures and normals we define the integrals over the surface as ```python ufl_form_srf = ufl.inner(F, bound_normal) * ds_int + ufl.inner(F, facet_normal) * ds ``` Finally, we create the DOLFINx forms for the volumetric and surface integrals using the {py:class}`form_custom` functionality from QUGaR ```python form_vol = cast(CustomForm, form_custom(ufl_form_vol, dtype=dtype)) form_srf = cast(CustomForm, form_custom(ufl_form_srf, dtype=dtype)) ``` and assemble them using the custom quadrature rules generated by QUGaR and fed to DOLFINx at runtime through the coefficients. ```python quad_gen = create_quadrature_generator(unf_domain) vol_intgr = dolfinx.fem.assemble_scalar(form_vol, coeffs=form_vol.pack_coefficients(quad_gen)) srf_intgr = dolfinx.fem.assemble_scalar(form_srf, coeffs=form_srf.pack_coefficients(quad_gen)) ``` Finally, we compare both integrals ```python print(f" - Volumetric integral: {vol_intgr}") print(f" - Surface integral: {srf_intgr}") ```