Divergence theorem

This demo is implemented in demo_div_thm.py. It illustrates:

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 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.

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:

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.

dtype = np.float64

We first define the geometry to be considered: a disk in 2D and a sphere in 3D. Uncomment the one you want to use.

# 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.

We create a Cartesian mesh (corresponding to \(\mathcal{T}\)) in which we will embbed the domain \(\Omega\),

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}}\).

unf_domain = qugar.impl.create_unfitted_impl_domain(impl_func, cart_mesh)

We create the vector function \(\mathbf{F}\) and its divergence.

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.

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.

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 dx_bdry_unf measure introduced in QUGaR (note that we integrate over the cut cells only)

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}}\).

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 mapped_normal function introduced in QUGaR.

bound_normal = mapped_normal(dlf_mesh)

while the normal at the external facets is defined using the standard UFL function.

facet_normal = ufl.FacetNormal(dlf_mesh)

Using these measures and normals we define the integrals over the surface as

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 form_custom functionality from QUGaR

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.

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

print(f"  - Volumetric integral: {vol_intgr}")
print(f"  - Surface integral: {srf_intgr}")