layout: true class: nord-dark, typo, typo-selection --- count: false class: center, middle # Curved Grid
.letter-spacing-20[A Dune module for surface parametrization]
by Florian Stenger and
Simon Praetorius
--- layout: true class: nord-light, typo, typo-selection --- # Initial Example ```c++ using namespace Dune; // 1. Construct a reference grid auto refGrid = GmshReader< FoamGrid<2,3> >::read("sphere.msh"); // 2. Define the geometry mapping auto sphere = [](const auto& x) { return x / x.two_norm(); }; auto sphereGF = analyticDiscreteFunction(sphere, *refGrid, order); // 3. Wrap the reference grid to build a curved grid CurvedGrid grid{*refGrid, sphereGF}; ``` --- # Initial Example ```c++ using namespace Dune; // 1. Construct a reference grid * auto refGrid = GmshReader< FoamGrid<2,3> >::read("sphere.msh"); // 2. Define the geometry mapping auto sphere = [](const auto& x) { return x / x.two_norm(); }; auto sphereGF = analyticDiscreteFunction(sphere, *refGrid, order); // 3. Wrap the reference grid to build a curved grid CurvedGrid grid{*refGrid, sphereGF}; ``` `1.` Read a grid from GMsh file and construct a FoamGrid with grid dimension 2 and world dimension 3. --- # Initial Example ```c++ using namespace Dune; // 1. Construct a reference grid auto refGrid = GmshReader< FoamGrid<2,3> >::read("sphere.msh"); // 2. Define the geometry mapping * auto sphere = [](const auto& x) { return x / x.two_norm(); }; * auto sphereGF = analyticDiscreteFunction(sphere, *refGrid, order); // 3. Wrap the reference grid to build a curved grid CurvedGrid grid{*refGrid, sphereGF}; ``` `2.` The surface parametrization is described by a (closest-point) projection that needs to be transformed into a *grid-function*. --- # Initial Example ```c++ using namespace Dune; // 1. Construct a reference grid auto refGrid = GmshReader< FoamGrid<2,3> >::read("sphere.msh"); // 2. Define the geometry mapping auto sphere = [](const auto& x) { return x / x.two_norm(); }; auto sphereGF = analyticDiscreteFunction(sphere, *refGrid, order); // 3. Wrap the reference grid to build a curved grid * CurvedGrid grid{*refGrid, sphereGF}; ``` `3.` The CurvedGrid wrapper takes a grid and a grid-function defined on that grid. The grid-function is used to describe the geometry mappings --- # A Reference Surface - Consider an element \(e\) in a reference surface \(\Gamma_h\) - (...) - (...) .center[
] --- # A Reference Surface - Consider an element \(e\) in a reference surface \(\Gamma_h\) - Mapping from reference element \(\hat{e}\) to element \(e\) denoted by \(\mu_e\) - (...) .center[
] --- # Surface Parametrization - Consider an element \(e\) in a reference surface \(\Gamma_h\) - Mapping from reference element \(\hat{e}\) to element \(e\) denoted by \(\mu_e\) - Curved surface \(\Gamma\) parametrized by \(\mathbf{X}\), or reference mapping \(\mathbf{X}_e\) .center[
] --- # `CurvedGeometry` - A Wrapper Class - **Input**: Surface \(\Gamma\) representation
Geometry parametrization: \(\mathbf{X}_e:\hat{e}\to\Gamma\subset\mathbb{R}^n\)
- **Output**: Geometry class with member functions:
`global(local)`: \(\mathbf{X}_e:\hat{e}\to\Gamma\)
--- # `CurvedGeometry` - A Wrapper Class - **Input**: Surface \(\Gamma\) representation
(closest-point) projection: \(\mathbf{X}:e\to\Gamma\subset\mathbb{R}^n\)
reference-mapping: \(\mu_e:\hat{e}\to e\)
- **Output**: Geometry class with member functions:
`global(local)`: \(\mathbf{X}_e:\hat{e}\to\Gamma;\;\mathbf{X}_e=\mathbf{X}\circ\mu_e\)
--- # `CurvedGeometry` - A Wrapper Class - **Input**: Surface \(\Gamma\) representation
(closest-point) projection: \(\mathbf{X}:e\to\Gamma\subset\mathbb{R}^n\)
Jacobian of projection: \(D\mathbf{X}:e\to\mathbb{R}^{n\times n}\)
reference-mapping: \(\mu_e:\hat{e}\to e,\; D^T\mu_e:\hat{e}\to\mathbb{R}^{m\times n}\)
- **Output**: Geometry class with member functions:
`global(local)`: \(\mathbf{X}_e:\hat{e}\to\Gamma;\;\mathbf{X}_e=\mathbf{X}\circ\mu_e\)
`jacobianTransposed(local)`: \(D^T X_e:\hat{e}\to \mathbb{R}^{m\times n}\)
\(D^T X_e=D^T\mu_e\cdot(DX\circ\mu_e)^T\)
`jacobianInverseTransposed(local)`: \(D^{-T} X_e:\hat{e}\to \mathbb{R}^{n\times m}\)
--- # `CurvedGeometry` - A Wrapper Class - Mapping of sub-entities \(s\preceq e\) or intersections \(s=\bar{e}_1\cap \bar{e}_2\) with \(dim(s)=l\leq m\) - (...) - (...)
--- # `CurvedGeometry` - A Wrapper Class - Mapping of sub-entities \(s\preceq e\) or intersections \(s=\bar{e}_1\cap \bar{e}_2\) with \(dim(s)=l\leq m\) - Mapping of sub-entity coordinates to element coordinates: \(\eta_{s,e}:\hat{s}\to\hat{e}\) - (...)
--- # `CurvedGeometry` - A Wrapper Class - Mapping of sub-entities \(s\preceq e\) or intersections \(s=\bar{e}_1\cap \bar{e}_2\) with \(dim(s)=l\leq m\) - Mapping of sub-entity coordinates to element coordinates: \(\eta_{s,e}:\hat{s}\to\hat{e}\) - Diagram commutes, i.e, \(\mu_s=\mu_e\circ\eta_{s,e}\), if grid is *twist-free*
--- # `CurvedGeometry` - A Wrapper Class - **Input**: Surface \(\Gamma\) representation
(closest-point) projection: \(\mathbf{X}:e\to\Gamma\subset\mathbb{R}^n\)
Jacobian of projection: \(D\mathbf{X}:e\to\Gamma\subset\mathbb{R}^{n\times n}\)
reference-mapping: \(\mu_e:\hat{e}\to e,\; D^T\mu_e:\hat{e}\to\mathbb{R}^{m\times n}\)
sub-entity geometry: \(\eta_{s,e}:\hat{s}\to\hat{e},\; D^T\eta_{s,e}:\hat{s}\to\mathbb{R}^{l\times m}\)
- **Output**: Geometry class with member functions:
`global(local)`: \(\mathbf{X}_s:\hat{s}\to\Gamma;\;\mathbf{X}_s=\mathbf{X}\circ\mu_e\circ\eta_{s,e}\)
`jacobianTransposed(local)`: \(D^T X_s:\hat{s}\to \mathbb{R}^{l\times n}\)
\(D^T X_s=D^T\eta_{s,e}\cdot(D^T\mu_e\circ\eta_{s,e})\cdot(DX\circ\mu_e\circ\eta_{s,e})^T\)
`jacobianInverseTransposed(local)`: \(D^{-T} X_s:\hat{s}\to \mathbb{R}^{n\times l}\)
--- # Two Implementations 1. `ParametrizedGeometry` 2. `LocalFunctionGeometry` --- # 1. `ParametrizedGeometry` Interpolation \(\mathbb{I}_h^{(k)}\) into Lagrange basis functions \(\{\phi_j^{(k)}\}\) of order \(k\) \[ \mathbf{X}^{(k)} := \mathbb{I}_h^{(k)}\mathbf{X}\text{ with } \\ \mathbf{X}^{(k)}(\mu_e(\hat{x})) = \mathbf{X}_e^{(k)}(\hat{x}) := \sum_{j=1}^{n_k} X_j^{(k)}\cdot \phi^{(k)}_j(\hat{x}),\quad\text{for }\hat{x}\in\hat{e}\,, \] with \(X_j^{(k)} := \mathbf{X}_e(\hat{x}^{(k)}_j)\in\Gamma\) for \(\{\hat{x}^{(k)}_j\}_{j=1\ldots n_k}\) local Lagrange nodes associated to \(\phi_j\). .center[
] --- # 1. `ParametrizedGeometry` - The local parametrization \(\mathbf{X}_e^{(k)}\) is a k-th order approximation \(\Gamma_h^{(k)}\) of \(\Gamma\) - Local basis functions define the Jacobian of the parametrization by \[ D\mathbf{X}_e^{(k)} = \sum_j^{n_k} X_j^{(k)}\otimes D\phi_j^{(k)}\,. \] .font-m[ ```c++ //
template
class ParametrizedGeometry; // constructors: template
ParametrizedGeometry(ReferenceElement, const LocalFiniteElement&, const Mapping& X_e); ParametrizedGeometry(ReferenceElement, const LocalFiniteElement&, std::vector
Xj); ``` ] --- # `ParametrizedGeometry` - Example ```c++ // define the geometry projection auto X = SphereProjection
{1.0}; // Construct a local finite-element cache, of order k LagrangeLFECache
lfeCache(k); // traverse the reference grid for (const auto& e : elements(refGrid->leafGridView())) { // projection from local coordinates auto X_e = [&X, mu_e=e.geometry()](const auto& local) { return X(mu_e.global(local)); }; // construct the CurvedGeometry from local parametrization const auto& lfe = lfeCache.get(e.type()); ParametrizedGeometry curvedGeo(e.type(), lfe, X_e); } ``` --- # 2. `LocalFunctionGeometry` - Function \(f:e\to\mathcal{R}\) defined on grid elements and localizable to elements is called *grid-function*. - Local-function \(f_e\) defined by grid-function as \(f_e:\hat{e}\to\mathcal{R};\;f_e=f\circ\mu_e\), - or grid-function \(f\) defined by local-function \(f_e\) as \(f|_e=f_e\circ\mu_e^{-1}\) -- - If \(f\) is differentiable, so is \(f_e\) with localized Jacobian \[ (Df)_e:\hat{e}\to L(e,\mathcal{R})\cong\mathbb{R}^{dim(\mathcal{R})\times n};\; (Df)_e=Df\circ\mu_e\,, \] -- - For sub-entity mappings, we need \(\eta_{s,e}\) and its Jacobian, called `LocalGeometry`. --- # 2. `LocalFunctionGeometry` Implementation class parametrized with the type of \(f_e\) and \(\eta_{s,e}\): ```c++ //
template
class LocalFunctionGeometry; template
using ElementLocalFunctionGeometry = LocalFunctionGeometry
, LocalFunctionGeometryTraits
>; ``` --- # `LocalFunctionGeometry` - Example ```c++ // Define the geometry grid-function auto X = SphereProjection
{1.0}; auto f = analyticGridFunction
>(X); // Define a local-function from the grid-function auto f_e = localFunction(f); auto eta = DefaultLocalGeometry
{}; // traverse the reference grid for (const auto& e : elements(refGrid->leafGridView())) { // bind local-function to element f_e.bind(e); // construct the Geometry from f_e and eta LocalFunctionGeometry localFctGeo(e.type(), f_e, eta); } ``` --- # Examples of local-functions - `AnalyticGridFunction`, expects a differentiable function \(X\) as input. Wrapps \(X\) and \(DX=derivative(X)\) and builds grid-functions and local-function interface. - `AnalyticDiscreteFunction`, expects a function \(X\) and a polynomial order \(k\). Performs local interpolation into Lagrange space on `localFunction.bind(element)`. Similar to `LagrangeCurvedGeometry`. - `DiscreteGridViewFunction`, a global basis and global coefficient vector that performs evaluation of local basis functions and linear combination with localized coefficients when bound to an element. Similar to `DiscreteGlobalBasisFunction` from dune-functions (but differentiable) and to `DiscreteFunction` from AMDiS (or `valueOf(DOFVector)` from old AMDiS). --- # Examples of Geometries - Some example geometries with projection \(\mathbf{X}(x)\), normal vector \(n(x)\), mean curvature \(H(x)\), surface area, and some provide the derivative \(D\mathbf{X}(x)\) and gaussian curvature \(K(x)\). - `SphereProjection`, `EllipsoidProjection`, and `TorusProjection` - `ExplicitSurfaceProjection`: use the closest-point projection to a triangulated grid \(\Gamma_h^\text{fine}\approx\Gamma\). Implemented by fast tree-based search for nearest neighbor coordinates, using a KDTree implementation, and a linear interpolation in the fine grid. .pure-g[ .pure-u-1-3.right[
] .pure-u-1-3.center[ \(\Gamma_h\overset{\mathbf{X}}{\longrightarrow}\Gamma_h^\text{fine}\) ] .pure-u-1-3.left[
] ] --- # Examples of Geometries - `[Simple]ImplicitSurfaceProjection`: If surface \(\Gamma\) is described by zero-level set of function \(\psi:\mathbb{R}^n\to\mathbb{R}\), i.e., \( \Gamma = \{x\in\mathbb{R}^n\,:\,\psi(x)=0\} \) an approximation of the closest-point projection of a point \(x_0\in\Gamma_h\) can be obtained iteratively, by (1) .footnote1[(1) I. Nitschke. Diskretes Äußeres Kalkül (DEC) auf Oberflächen ohne Rand, 2014] \[ x_{i+1} = x_i - \nabla\psi(x_i) \frac{\psi(x_i)}{\|\nabla\psi(x_i)\|^2}, \] or with correction of the distance (2) .footnote3[(2) A. Demlow and G. Dziuk. An adaptive finite element method for the laplace–beltrami operator on implicitly defined surfaces., 2007] \[ \tilde{x}_{i+1} = x_i - \nabla\psi(x_i) \frac{\psi(x_i)}{\|\nabla\psi(x_i)\|^2},\\ \text{dist} = \operatorname{sign}(\psi(x_0))\|\tilde{x}_{i+1}-x_0\|, \\ x_{i+1} = x_0 - \nabla\psi(\tilde{x}_{i+1}) \frac{\text{dist}}{\|\nabla\psi(\tilde{x}_{i+1})\|} \] --- # Examples of Geometries Surface of genus 2: \[ \psi(x,y,z) = 2y(y^{2}-3x^{2})(1-z^{2})+(x^{2}+y^{2})^{2}-(9z^{2}-1)(1-z^{2}) \] .center[
] --- # `CurvedGrid` - A Wrapper Class - A Grid provides a traversal over a (subset of) elements. - Each element can return a geometry - The `CurvedGrid` traverses a reference grid but returns on each element a curved geometry: ```c++ //
template
class CurvedGrid; ``` It is parametrized with a grid-function. --- # `CurvedGrid` - A Wrapper Class - Construct either a `ParametrizedGeometry` with Lagrange basis function (if order `k > 0`)... ```c++ // Construct a reference grid auto refGrid = GmshReader< FoamGrid<2,3> >::read("sphere.msh"); // Define the geometry mapping auto sphere = [](const auto& x) { return x / x.two_norm(); }; template
using Order_t = std::integral_constant
; // Wrap the reference grid to build a curved grid CurvedGrid grid{*refGrid, sphere, Order_t
{}}; ``` --- # `CurvedGrid` - A Wrapper Class - or construct a `LocalFunctionGeometry` if the grid-function is differentiable ```c++ auto gv = refGrid->leafGridView(); // Define a discrete grid-function on the reference grid // with dim(range) = 3 and polynomial order k auto f = discreteGridViewFunction<3>(gv, k); // Interpolate the parametrization into the grid-function Functions::interpolate(f.basis(), f.coefficients(), sphere); // Wrap the reference grid and the grid-function CurvedGrid grid{*refGrid, f}; ``` --- # Validation of Geometrie - Geometric error bounds: for \(h\) small enough an \(\Gamma\) smooth: .footnote[A. Demlow. Higher-order finite element methods and pointwise error estimates for elliptic problems on surfaces., 2009] \[ \|\mathbf{X} - \mathbf{X}^{(k)}\|_{L^\infty(\Gamma_h^{(k)})} \leq C h^{k+1}, \\ \|n\circ \mathbf{X} - n_h^k\|_{L^\infty(\Gamma_h^{(k)})} \leq C h^k, \\ \|H\circ \mathbf{X} - H_h^k\|_{L^\infty(\Gamma_h^{(k)})} \leq C h^{k-1} \] .center[
] --- # Integration into AMDiS - AMDiS defines a problem class to store a PDE description, the finite-element spaces and the discrete solution (upon others) - Simply exchange the grid with `CurvedGrid` ```c++ using namespace Dune; using namespace AMDiS; using namespace Dune::Functions::BasisFactory; // Read parameters for grid creation from init-file auto refGrid = GmshReader< AlbertaGrid<2,3> >::read("sphere.msh"); // Construct the curved grid auto X = SphereProjection
{1.0}; CurvedGrid grid(*refGrid, X); // Construct the ProblemStat ProblemStat prob("prob", grid, lagrange<2>()); ``` --- # **Example 1**: Helmholtz equation Tangential vector (tensor) fields on curved surfaces described by surface Helmholtz equation: Find \(\mathbf{u}\in H^1(\Gamma,T^{(r,s)}(\Gamma))\), s.t. \[ \big(\nabla \mathbf{u},\,\nabla \mathbf{v}\big)_{\Gamma} + \big(\mathbf{u},\,\mathbf{v}\big)_{\Gamma} = \big(\mathbf{f},\, \mathbf{v}\big)_{\Gamma} \] for all \(\mathbf{v}\in H^1(\Gamma,T^{(r,s)}(\Gamma))\) with \(\nabla\mathbf{u}\) the covariant derivative. -- **Discretization**: - Tensor extended to ambient space \(\Rightarrow T^{(r,s)}(\mathbb{R}^n)\cong[\mathbb{R}^n]^{(r+s)}\) (non-tangential) with tensor-product finite-element space \(V_h^{(r,s)}:=[V_h]^{n^{(r+s)}}\). - Discrete function space \(V_h\) for scalar fields on \(\Gamma_h^{(k_g)}\) \[ V_h := \{v\in C^0(\Gamma_h^{(k_g)})\,:\,v\circ \mathbf{X}^{(k_g)}_e \in\mathbb{P}_{k_u}(\hat{e}),\;\forall e\in\mathcal{G}_h\} \] .footnote1[Michael Nestler, Ingo Nitschke and Axel Voigt, *A finite element approach for vector- and tensor-valued surface PDEs*, (2019)] .footnote2[Hanne Hardering and Simon Praetorius, *Higher-Order Surface Finite-Elements for the Tensor-Laplacian*, (in preparation)] --- # **Example 1**: Helmholtz equation - Enforce tangentiality by penalty term \(\big(Q \mathbf{u},\,Q \mathbf{v}\big)_{\Gamma}\) - Tangential projection \(P:=P^1\otimes\!\ldots\!\otimes P^1\otimes P_1\otimes\!\ldots\!\otimes P_1\), normal projection \(Q = Id-P\) - Modified Helmholtz equation: Find \(\mathbf{u}\in V_h^{(r,s)}\), s.t. .font-md[ \[ \big(\nabla P \mathbf{u},\,\nabla P \mathbf{v}\big)_{\Gamma_h^{(k_g)}} + \big(P \mathbf{u},\,P \mathbf{v}\big)_{\Gamma_h^{(k_g)}} + \beta h^{-2\alpha}\big(\tilde{Q} \mathbf{u},\,\tilde{Q} \mathbf{v}\big)_{\Gamma_h^{(k_g)}} = \big(\mathbf{f}^e,\, P \mathbf{v}\big)_{\Gamma_h^{(k_g)}} \] ] for all \(\mathbf{v}\in V_h^{(r,s)}\) where projections \(P^1,P_1\) are defined w.r.t. surface normals \(n_h\perp\Gamma_h^{(k_g)}\) and \(\tilde{Q}\) is defined w.r.t. "better" surface normal \(\tilde{n}_h\) of order \(k_p\). - The surface covariant derivative on \(\Gamma_h^{(k_g)}\) is given by \[ \nabla P \mathbf{u} = P D\mathbf{u} - P DQ \mathbf{u} \] --- # **Example 1**: Helmholtz equation Convergence orders for \(k:=k_u=k_g\) with \(k_p=k+1\)
--- # **Example 1**: Helmholtz equation - Convergence orders: \(k_u=k_g=1\) normal order \(k_p=2\) .font-sm.mb-xs[ | level | h | err(L2) | eoc(L2) | err(Ah) | eoc(Ah) | | ----- | -------- | -------- | ------- | --------- | -------- | | 0 | 7.07e-01 | 3.40e-01 | | 5.20e-01 | | | 2 | 2.05e-01 | 3.43e-02 | 2.07 | 8.18e-02 | 1.45 | | 4 | 5.20e-02 | 2.20e-03 | 2.00 | 1.72e-02 | 1.07 | ] - Convergence orders: \(k_u=k_g=4\) normal order \(k_p=5\) .font-sm.mb-xs[ | level | h | err(L2) | eoc(L2) | err(Ah) | eoc(Ah) | | ----- | -------- | -------- | -------- | -------- | --------- | | 0 | 7.46e-01 | 4.24e-02 | | 2.19e-01 | | | 2 | 2.06e-01 | 7.60e-05 | 4.04 | 1.32e-04 | 3.50 | | 4 | 5.20e-02 | 4.40e-08 | 5.48 | 3.62e-06 | 4.31 | ] --- # **Example 2**: Mean Curvature Flow Parametrization of a family of surfaces \(\Gamma(t)\subset\mathbb{R}^3\) over \(\Gamma_h\) \[ \Gamma(t) = \left\{\mathbf{X}(p,t)\,:\,p\in\Gamma_h\right\}\,. \] Evolution of coordinates by *mean-curvature-flow* \[ \partial_t \mathbf{X} = -H(\mathbf{X}) n(\mathbf{X}) = \Delta_\Gamma \mathbf{X}\,. \] Weak formulation: For all \(t\in [0,T]\), find \(\mathbf{X}(\cdot,t)\in [H^1(\Gamma_h)]^3\) s.t. \[ \int_{\Gamma(t)} \partial_t x(t)\cdot y\,\textrm{d}S = -\int_{\Gamma(t)}\sum_{i}\nabla_{\Gamma(t)}x_i(t)\cdot\nabla_{\Gamma(t)}y_i\,\textrm{d}S,\;\forall y\in [H^1(\Gamma(t))]^3 \] with \(x=\mathbf{X}(p,t)\) for \(p\in\Gamma_h\). --- # **Example 2**: Mean Curvature Flow Discretization by surface finite-elements \[ V_h := \{v\in C^0(\Gamma_h)\,:\,v\circ \mu_e \in\mathbb{P}_k(\hat{e}),\;\forall e\in\mathcal{G}_h\} \] Define (curved) surface \(\Gamma_t:=\mathbf{X}_t(\Gamma_h)\) by grid-function \(\mathbf{X}_t\) Weak formulation: For all \(t\in [0,T]\), find \(\mathbf{X}_t\in [V_h]^3\) s.t. \[ \int_{\Gamma_h} \partial_t \mathbf{X}_t\cdot \mathbf{Y}\,\textrm{d}S(\Gamma_t) = -\int_{\Gamma_h}\sum_{i}\nabla_{\Gamma_t}X_t^i\cdot\nabla_{\Gamma_t}Y^i\,\textrm{d}S(\Gamma_t),\;\forall \mathbf{Y}\in [V_h]^3 \] .center[
] --- # Summary and Outlook
Summary
Curved geometries implemented in
ParametrizedGeometry
and
LocalFunctionGeometry
Mapping from intersection, sub-entities, and elements to curved surface
Curved grid wrapper implemented in
CurvedGrid
Can be used as replacement for any Dune grid.
--- # Summary and Outlook
Summary
Curved geometries implemented in
ParametrizedGeometry
and
LocalFunctionGeometry
Mapping from intersection, sub-entities, and elements to curved surface
Curved grid wrapper implemented in
CurvedGrid
Can be used as replacement for any Dune grid.
Outlook
Volume grid with only curved boundaries.
Use B-Spline basis for parametrization.
Automatic differentiation in implicit geometries.
--- layout: true class: nord-dark, typo, typo-selection --- count: false class: center, middle # Find the code on
[gitlab.mn.tu-dresden.de/...](https://gitlab.mn.tu-dresden.de) CurvedGeometry: [...iwr/dune-curvedgeometry](https://gitlab.mn.tu-dresden.de/iwr/dune-curvedgeometry) CurvedGrid: [...iwr/dune-curvedgrid](https://gitlab.mn.tu-dresden.de/iwr/dune-curvedgrid) AMDiS Integration: [...amdis/amdis-extensions](https://gitlab.mn.tu-dresden.de/amdis/amdis-extensions)