deal.II version GIT relicensing-2652-g41e7496fd4 2025-02-17 19:10:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
The step-90 tutorial program

This tutorial depends on step-85.

Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program

This program was contributed by Vladimir Yushutin and Timo Heister, Clemson University, 2023.

This material is based upon work partly supported by the National Science Foundation Award DMS-2028346, OAC-2015848, EAR-1925575, and by the Computational Infrastructure in Geodynamics initiative (CIG), through the NSF under Award EAR-0949446, EAR-1550901, EAR-2149126 via the University of California – Davis.

Introduction

In this tutorial, we implement the trace finite element method (TraceFEM) in deal.II. TraceFEM solves PDEs posed on a possibly evolving \((dim-1)\)-dimensional surface \(\Gamma\) employing a fixed uniform background mesh of a \(dim\)-dimensional domain in which the surface is embedded. Such surface PDEs arise in problems involving material films with complex properties and in other situations in which a non-trivial condition is imposed on either a stationary or a moving interface. Here we consider a steady, complex, non-trivial surface and the prototypical Laplace-Beltrami equation which is a counterpart of the Poisson problem on flat domains.

Being an unfitted method, TraceFEM allows to circumvent the need of remeshing of an evolving surface if it is implicitly given by the zero contour of a level-set function. At the same time, it easily provides with an extension of the discrete solution to a neighborhood of the surface which turns out to be very handy in case of non-stationary interfaces and films. Certainly, this flexibility comes with a price: one needs to design the nodes and weights for a quadrature customized for each implicit intersection of the zero level-set and the background mesh. Moreover, these intersections may be of arbitrary shape and size manifesting in the so-called "small cut" problem and requiring addition of a stabilization form that restores well-conditioning of the problem.

Two aspects are of our focus. First, the surface approximation is separated from the discretization of the surface PDE, e.g., a \(Q_2\) discrete level-set and a \(Q_1\) solution are possible on the same bulk triangulation. Second, we make sure that the performance of TraceFEM in the parallel implementation corresponds to that of a classical fitted FEM for a two-dimensional problem. We demonstrate how to achieve both goals by using a combination of MeshWorker and NonMatching capabilities.

A natural alternative to TraceFEM in solving surface PDEs is the parametric surface finite element method. The latter method relies on an explicit parametrization of the surface which may be not feasible especially for evolving interfaces with an unknown in advance shape - in this sense, TraceFEM is a technique inspired by the level-set description of interfaces. However, the parametric surface finite element method, when applicable, enjoys many well-known properties of fitted methods on flat domains provided the geometric errors - which a present for both methods - are taken under control.

A non-trivial surface

A fitted FEM on a flat two-dimensional domain, if discretized by piecewise linears with \(N\) degrees of freedom, typically results in \(O(h)=O(N^{-1/2})\) convergence rate of the energy error; requires \(O(N)\) storage for the degrees of freedom; and, more importantly, takes \(O(N)\) of construction time to create them, i.e. to mesh the domain. TraceFEM, although solving a two-dimensional problem, relies on the inherently three-dimensional mesh on which the level-set function must be defined and, if implemented naively, suffers from the increased storage and the increased construction time in terms of the active degrees of freedom \(N_a\) that actually enters the scheme with, hopefully, \(O(N_a^{-1/2})\) error. To combat these possible bottlenecks, we create iteratively a mesh which is localized near the zero contour line of the level set function, i.e near the surface, to restore the aforementioned two-dimensional performance typical for fitted FEM, see the first three typical iterations of this methodology below.

Iterative localization of the zero contour of a typical level set

The cells colored by red cary the active degrees of freedom (total number \(N_a\)) as the level set is not sign-definite at support points. Notice also that the mesh is graded: any cell has at most 4 neighbors adjacent to each of 6 faces.

Once a desired geometry approximation \(\Gamma_h\) is achieved using the iterative approach above, we can start forming the linear system using the constructed normals and quadratures. For the purposes of the tutorial we choose a non-trivial surface \(\Gamma\) given by

\begin{equation*} \frac{x^2}{4}+ y^2 + \frac{4 z^2} {(1 + 0.5 \sin(\pi x))^{2}} = 1 \end{equation*}

The OY and OX views of this tamarind-shaped, exact surface \(\Gamma\) are shown below along with the mesh after three iterations (the approximation \(\Gamma_h\) is not shown).

OY(left) and OZ(right) cross-sections of the background mesh along with the exact surface

Model problem

We would like to solve the simplest possible problem defined on a surface, namely the Laplace–Beltrami equation,

\begin{equation*} -\Delta_\Gamma u + c u = f \qquad \text{in }\, \Gamma, \end{equation*}

where we take \(c=1\) for concreteness. We added the term \(cu\) to the left-hand side so the problem becomes well-posed in the absence of any boundary; an alternative could be to take \(c=0\) but impose the zero mean condition.

Manufactured exact solution

We choose the test solution and the right-hand side forcing as the restriction to \(\Gamma\) of

\begin{equation*} u(x,y,z)=xy\,,\quad f(x,y,z)=xy + 2.0\,\mathbf{n}_x \mathbf{n}_y + \kappa (y \mathbf{n}_x + x\mathbf{n}_y), \end{equation*}

where the latter is manufactured using the exact normal \(\mathbf{n}\), the exact Hessian \(\nabla^2\mathbf{n}\) and the mean curvature, \(\kappa=\mathrm{div} n\) of the surface. Note that we do not need to impose any boundary conditions as the surface \(\Gamma\) is closed.

The Trace Finite Element Method

TraceFEM is an unfitted method: the surface \(\Gamma\) is immersed into a regular, uniform background mesh that stays fixed even if the surface would be evolving. To solve Laplace–Beltrami equation, we first construct a surface approximation \(\Gamma_h\) by intersecting implicitly the cells of the background mesh with the iso surface of an approximation of the level-set field. We note that we never actually create any two-dimensional meshes for the surface but only compute approximate quadrature points and surface normals. Next we distribute degrees of freedom over a thin subdomain \(\Omega_h\) that completely covers \(\Gamma_h\) and that consists of the intersected cells \(\mathcal{T}_\Gamma^h\),

\begin{equation*} \mathcal{T}_\Gamma^h = \{ T \in \mathcal{T}^{h} : T \cap \Gamma_h \neq \emptyset \}. \end{equation*}

The finite element space where we want to find our numerical solution, \(u_h\), is now

\begin{equation*} V_h = \{ v \in C(\Omega_h) : v \in Q_p(T), \, T \in \mathcal{T}_\Gamma^h \}, \end{equation*}

where \(\Omega_h\) is the union of all intersected cells from \(\bigcup_{T \in \mathcal{T}_\Gamma^h} \overline{T}\).

To create \(V_h\), we first add an FE_Q and an FE_Nothing element to an hp::FECollection. We then iterate over each cell \(T\) and, depending on whether \(T\) belongs to \(\mathcal{T}_\Gamma^h\) or not, we set the active_fe_index to either 0 or 1. To determine whether a cell is intersected or not, we use the class NonMatching::MeshClassifier.

A natural candidate for a weak formulation involves the following (bi)linear forms

\begin{align*} a_h(u_h, v_h) = (\nabla_{\Gamma_h} u_h, \nabla_{\Gamma_h} v_h)_{\Gamma_h}+(u_h, v_h)_{\Gamma_h}\,,\qquad L_h(v_h) = (f^e,v_h)_{\Gamma_h}. \end{align*}

where \(f^e\) is an extension (non-necessarily the the so-called normal extension) of \(f\) from \(\Gamma\) to \(\Omega_h\). Note that the right-hand side \(f\) of the Laplace-Beltrami problem is defined on the exact surface \(\Gamma\) only and we need to specify how to evaluate its action on the perturbed approximate geometry \(\Gamma_h\) which is immersed in \(\Omega_h\). For the purposes of this test, the forcing \(f\) is manufactured using \(u=xy\) and the level-set function and, therefore, is a function of Cartesian coordinates \(x\), \(y\), \(z\). The latter is identified with \(f^e\) on \(\Gamma_h\) and it is not the normal extension of the function \(f\).

However, the so-called "small-cut problem" may arise and one should introduce the stabilized version of TraceFEM: Find \(u_h \in V_h\) such that

\begin{equation*} a_h(u_h,v_h) + s_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h. \end{equation*}

Here the normal-gradient stabilization \(s_h\) involves the three-dimensional integration over whole (but intersected) cells and is given by

\begin{equation*} s_h(u_h,v_h) = h^{-1}(\mathbf{n}_h\cdot\nabla u_h, \mathbf{n}_h\cdot\nabla v_h)_{\Omega_h}, \end{equation*}

Note that the \(h^{-1}\) scaling may be relaxed for sufficiently smooth solutions such as the manufactured one, but we choose the strong scaling to demonstrate the extreme case [171].

Discrete Level Set Function

In TraceFEM we construct the approximation \(\Gamma_h\) using the interpolant \(\psi_h\) of the exact level-set function on the bulk triangulation:

\begin{align*} \Gamma_h &= \{x \in \mathbb{R}^{\text{3}} : \psi_h(x) = 0 \}. \end{align*}

The exact normal vector \(\mathbf{n}\) is approximated by \(\mathbf{n}_h=\nabla\psi_h/\|\nabla\psi_h\|\) which, together with approximate quadrature for the integration over \(\Gamma_h\), leads to the so-called "geometrical error". Luckily, one can show [171] that the method converges optimally for the model problem if the same element space \(V_h\) is employed for the discrete functions and for the interpolation of the level set function as if the exact domain would have been used. Furthermore, deal.II allows to choose independently the discrete space for the solution and a higher-order discrete space for the level set function for a more accurate geometric approximation.

The commented program

  #include <deal.II/base/function.h>
  #include <deal.II/base/numbers.h>
  #include <deal.II/base/point.h>
  #include <deal.II/base/quadrature.h>
  #include <deal.II/base/tensor.h>
  #include <deal.II/base/timer.h>
  #include <deal.II/distributed/grid_refinement.h>
  #include <deal.II/distributed/tria.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/fe/fe_nothing.h>
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/grid/grid_tools.h>
  #include <deal.II/grid/tria.h>
  #include <deal.II/hp/fe_collection.h>
  #include <deal.II/lac/affine_constraints.h>
  #include <deal.II/lac/precondition.h>
  #include <deal.II/lac/solver_control.h>
  #include <deal.II/lac/sparsity_pattern.h>
  #include <deal.II/lac/trilinos_vector.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/meshworker/mesh_loop.h>
  #include <deal.II/meshworker/scratch_data.h>
  #include <deal.II/non_matching/fe_values.h>
  #include <deal.II/non_matching/mesh_classifier.h>
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/numerics/error_estimator.h>
  #include <deal.II/numerics/vector_tools.h>
  namespace Step90
  {
  using namespace dealii;
  using VectorType = TrilinosWrappers::MPI::Vector;
Definition hp.h:117

The parallelization in this tutorial relies on the Trilinos library. We will grant to some cells empty finite element spaces FE_Nothing as done in step-85, but this time active DoFs will be only assigned to cell which are intersected by the surface approximation.

  {
  nothing = 1
  };

Exact surface

The following class defines the surface using the implicit level set representation. The exact surface normal uses the Cartesian gradient of the level set function. The exact Hessian is needed for the construction of the test case only.

  template <int dim>
  class TamarindShape : public Function<dim>
  {
  public:
  : Function<dim>()
  {}
  double value(const Point<dim> &point,
  const unsigned int component = 0) const override
  {
  AssertIndexRange(component, this->n_components);
  (void)component;
  return 0.25 * Utilities::pow(point[0], 2) + Utilities::pow(point[1], 2) +
  4.0 * Utilities::pow(point[2], 2) *
  std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -2) -
  1.0;
  }
  const unsigned int component = 0) const override
  {
  AssertIndexRange(component, this->n_components);
  (void)component;
  grad[0] = 0.5 * point[0] +
  (-2.0) * 4.0 * Utilities::pow(point[2], 2) *
  std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -3) *
  (0.5 * numbers::PI * std::cos(numbers::PI * point[0]));
  grad[1] = 2.0 * point[1];
  grad[2] = (2.0) * 4.0 * point[2] *
  std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -2);
  return grad;
  }
  hessian(const Point<dim> &point,
  const unsigned int component = 0) const override
  {
  AssertIndexRange(component, this->n_components);
  (void)component;
  hessian[0][0] =
  0.5 +
  8.0 * Utilities::pow(point[2], 2) *
  (3.0 * std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -4) *
  std::cos(numbers::PI * point[0]),
  2) +
  std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -3) * 0.5 *
  hessian[0][1] = 0.0;
  hessian[0][2] =
  (-8.0) * point[2] *
  std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -3) *
  hessian[1][1] = 2.0;
  hessian[1][2] = 0.0;
  hessian[2][2] =
  8.0 * std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -2);
  return hessian;
  }
  };
friend class Tensor
Definition tensor.h:865
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:193
constexpr T pow(const T base, const int iexp)
Definition utilities.h:967
constexpr double PI
Definition numbers.h:241
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)

Exact solution

The following class defines the chosen exact solution and its surface gradient. The exact solution we try to reproduce is \(u=xy\) and it may be evaluated away from \(\Gamma\) as any other function of Cartesian points. Also note that the gradient() method returns the surface gradient \(\nabla_\Gamma u\) of the exact solution.

  template <int dim>
  class AnalyticalSolution : public Function<dim>
  {
  private:
  public:
  : Function<dim>()
  {}
  double value(const Point<dim> &point,
  const unsigned int component = 0) const override;
  Tensor<1, dim> gradient(const Point<dim> &point,
  const unsigned int component = 0) const override;
  };
  template <int dim>
  const unsigned int component) const
  {
  AssertIndexRange(component, this->n_components);
  (void)component;
  return point[0] * point[1];
  }
  template <int dim>
  const unsigned int component) const
  {
  AssertIndexRange(component, this->n_components);
  (void)component;
  const Tensor<1, dim> grad = tamarind.gradient(point, component);
  const Tensor<1, dim> normal = grad / grad.norm();
  }
numbers::NumberTraits< Number >::real_type norm() const

Exact forcing

We choose the right hand side equal to the evaluation of the surface Laplacian for a manufactured solution \(u\). This corresponds to the exact forcing \(f=-\Delta_\Gamma u+u\):

  template <int dim>
  class RightHandSide : public Function<dim>
  {
  public:
  : Function<dim>()
  {}
  virtual double value(const Point<dim> &p,
  const unsigned int component = 0) const override;
  };
  template <int dim>
  const unsigned int component) const
  {
  AssertIndexRange(component, this->n_components);
  (void)component;
  Assert(dim == 3, ExcNotImplemented());
  const Tensor<1, dim> grad = tamarind.gradient(point, component);
  const Tensor<1, dim> normal = grad / grad.norm();
  const SymmetricTensor<2, dim> hessian = tamarind.hessian(point, component);
  double mean_curv = 0.0;
  for (int j = 0; j < 3; j++)
  for (int k = 0; k < 3; k++)
  mean_curv += ((j == k ? 1 : 0) - normal[j] * normal[k]) * hessian[j][k];
  return point[0] * point[1] + 2.0 * normal[0] * normal[1] +
  mean_curv * (point[1] * normal[0] + point[0] * normal[1]);
  }

Scratch and Copy objects for TraceFEM

Since the assembly procedure will be performed via MeshWorker, we need a Scratch object that handles the Non-Matching FEValues effectively. The input arguments of its constructor are discussed in the solver class below.

  template <int dim>
  struct ScratchData
  {
  ScratchData(const Mapping<dim> &mapping,
  const hp::FECollection<dim> &fe_collection,
  const NonMatching::MeshClassifier<dim> &mesh_classifier,
  const VectorType &level_set,
  const Quadrature<dim> &quadrature,
  const UpdateFlags cell_update_flags = update_values |
  : fe_values(
  mapping,
  fe_collection[static_cast<types::fe_index>(ActiveFEIndex::lagrange)],
  quadrature,
  cell_update_flags)
  , region_update_flags(nonmatching_update_flags)
  , fe_collection(fe_collection)
  , mesh_classifier(mesh_classifier)
  , level_set(level_set)
  quadrature,
  cell_update_flags)
  , non_matching_fe_values(fe_collection,
  mesh_classifier,
  level_set)
  {}
  ScratchData(const ScratchData<dim> &scratch_data)
  : fe_values(scratch_data.fe_values.get_mapping(),
  scratch_data.fe_values.get_fe(),
  scratch_data.fe_values.get_quadrature(),
  scratch_data.fe_values.get_update_flags())
  , region_update_flags(scratch_data.region_update_flags)
  , fe_collection(scratch_data.fe_collection)
  , mesh_classifier(scratch_data.mesh_classifier)
  , level_set(scratch_data.level_set)
  , level_set_fe_values(scratch_data.level_set_fe_values.get_mapping(),
  scratch_data.level_set_fe_values.get_fe(),
  scratch_data.level_set_fe_values.get_quadrature(),
  scratch_data.level_set_fe_values.get_update_flags())
  , non_matching_fe_values(fe_collection,
  region_update_flags,
  mesh_classifier,
  level_set)
  {}
UpdateFlags
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
Definition types.h:32

The following FEValues object is used for the standard quadrature on cells involving the FE space of the solution. In TraceFEM, we need this quadrature due to the stabilization term. In addition, a cell quadrature for the FE space of the level set is defined.

The MeshWorker framework also requires a "copy" data structure that is filled by the worker function working on a cell or face, and whose contents are then later copied into global matrices and vectors. This CopyData object is customized for TraceFEM. In particular, the implementation of the normal-gradient volume stabilization relies on it.

  template <int dim>
  struct CopyData
  {
  FullMatrix<double> cell_matrix;
  std::vector<types::global_dof_index> local_dof_indices;
  void reinit(const typename DoFHandler<dim>::active_cell_iterator &cell)
  {
  const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
  cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
  cell_rhs.reinit(dofs_per_cell);
  local_dof_indices.resize(dofs_per_cell);
  cell->get_dof_indices(local_dof_indices);
  }
  };
  template <int dim>
  {
  unsigned int cell_index;
  {
  cell_index = cell->active_cell_index();
  }
  };
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
unsigned int cell_index
typename ActiveSelector::active_cell_iterator active_cell_iterator
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)

Normal-gradient stabilization form of TraceFEM

The following class corresponds to the stabilization form, its contribution to the global matrix and to the error. More specifically, the method needs_cell_worker() indicates whether the bilinear form of the stabilization, unlike the main bilinear form of Laplace-Beltrami operator, needs the bulk cell quadratures. The cell worker which is useful in an accumulation by MeshWorkers is provided by the assemble_cell_worker() method. The remaining method evaluate_cell_worker() computes the stabilization error for the solution \(u_h\), i.e \(s_h(u_h,u_h)\). Also note that the method needs_cell_worker() indicates that the assembly and the evaluation of the form does require a bulk cell quadrature. This methodology may be utilized in the MeshWorker. The stabilization scaling is specified by \(\mathrm{stabilization\_parameter}\cdot h^\mathrm{stabilization\_exponent}\). For elliptic problems with smooth solutions we can choose any \(-1\leq \mathrm{stabilization\_exponent} \leq 1\) and a sufficiently large \(\mathrm{stabilization\_parameter}\) that depends of \(\Gamma\).

We define the stabilization form here assuming that ScratchData and CopyData arguments are initialized properly. The local contribution of the stabilization from this cell to the global matrix is given in assemble_cell_worker() and, later in evaluate_cell_worker(), the local bilinear form of the stabilization is evaluated on the solution. Note the gradients of the discrete level set are computed in the bulk cell quadrature points, which, upon normalization, give the discrete normal vector in a bulk cell.

  VectorType &level_set,
  ScratchData<dim> &scratch_data,
  CopyData<dim> &copy_data) const
  {
  const FEValues<dim> &fe_values = scratch_data.fe_values;
  scratch_data.level_set_fe_values;
  const std::vector<double> &JxW_cell = fe_values.get_JxW_values();
  std::vector<Tensor<1, dim>> grad_level_set(
  level_set_fe_values.get_quadrature().size());
  level_set_fe_values.get_function_gradients(level_set, grad_level_set);
  const double factor =
  std::pow(cell->minimum_vertex_distance(), stabilization_exponent);
  for (const unsigned int q : fe_values.quadrature_point_indices())
  {
  const Tensor<1, dim> &normal =
  for (const unsigned int i : fe_values.dof_indices())
  for (const unsigned int j : fe_values.dof_indices())
  copy_data.cell_matrix(i, j) +=
  factor * (normal * fe_values.shape_grad(i, q)) *
  (normal * fe_values.shape_grad(j, q)) * JxW_cell[q];
  }
  }
  VectorType &solution,
  VectorType &level_set,
  ScratchData<dim> &scratch_data,
  CopyDataError<dim> &copy_data) const
  {
  double cell_stab_sqr = 0.0;
  const FEValues<dim> &fe_values = scratch_data.fe_values;
  const std::vector<double> &JxW_cell = fe_values.get_JxW_values();
  const unsigned int n_q_points = fe_values.get_quadrature_points().size();
  scratch_data.level_set_fe_values;
  std::vector<Tensor<1, dim>> level_set_grad(n_q_points);
  level_set_fe_values.get_function_gradients(level_set, level_set_grad);
  std::vector<Tensor<1, dim>> sol_grad(n_q_points);
  fe_values.get_function_gradients(solution, sol_grad);
  const double factor =
  std::pow(cell->minimum_vertex_distance(), stabilization_exponent);
  for (const unsigned int q : fe_values.quadrature_point_indices())
  {
  const Tensor<1, dim> normal =
  const double stabilization_at_point = normal * sol_grad[q];
  }
  copy_data.cell_stab_sqr = cell_stab_sqr;
  }
  private:
  };
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471

Laplace–Beltrami solver

The main class whose method run() performs the computation. One may adjust main parameters of TraceFEM in the constructor. The other methods are discussed below.

  template <int dim>
  {
  public:
  void run();
  private:
  void make_grid();
  void distribute_dofs();
  void initialize_matrices();
  void solve();
  void refine_grid();
  void output_level_set(unsigned int);
  MPI_Comm mpi_communicator;

The surface of interest corresponds to the zero contour of the following exact level set function:

The manufactured solution to the Laplace–Beltrami problem and the corresponding right-hand side:

There is a single triangulation which is shared by the discretizations of the solution and of the level set.

We need two separate FE spaces. The first manages the TraceFEM space which is active on intersected elements. The second manages the discrete level set function that describes the geometry of the surface. Also, the degrees of the FE spaces and the corresponding DoFHandler objects are given in the following:

  const unsigned int fe_degree;
  hp::FECollection<dim> fe_collection;
  DoFHandler<dim> dof_handler;
  const unsigned int level_set_fe_degree;
  const MappingQ1<dim> mapping;

Since we will adaptively refine the bulk triangulation, two constraints are needed: one for the solution space and another for the level set space.

Discrete vectors initialized with dof_handler and level_set_dof_handler.

The following NonMatching::MeshClassifier object is used to separate intersected elements and non-intersected ones. We will then use different finite elements from an hp::FECollection for these two categories:

The first bulk quadrature is required for the for TraceFEM stabilization, while the integration over implicit surface is based on the last, one-dimensional rule.

  const QGauss<dim> cell_quadrature;

Any TraceFEM needs a stabilization, and we choose the normal-gradient, volume stabilization.

Discrete right-hand side and the final matrix corresponding to dof_handler.

  VectorType global_rhs;
  MatrixType global_matrix;
  SparsityPattern sparsity_pattern;
  IndexSet locally_owned_dofs;

Depending on the type of the quadrature, surface, face or volume, we need to define different update flags.

The following variables are used to display the results of the convergence test:

  };
  template <int dim>
  : mpi_communicator(MPI_COMM_WORLD)
  , triangulation(mpi_communicator)
  (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
  , computing_timer(mpi_communicator,
  TimerOutput::never,
  TimerOutput::wall_times)
  , fe_degree(1)
  , fe_collection(FE_Q<dim>(fe_degree), FE_Nothing<dim>())
  , dof_handler(triangulation)
  , mapping()
  , mesh_classifier(level_set_dof_handler, level_set)
  , cell_quadrature(fe_degree + 1)
  , quadrature_1D(fe_degree + 1)
  {
  }
Definition fe_q.h:554
@ update_normal_vectors
Normal vectors.
STL namespace.

Geometric approximation

Let us start with a function that creates the background mesh, using a domain size chosen to avoid situations in which level set function vanishes at mesh vertices. The initial refinement helps the level set to approximate the surface meaningfully.

In following next method we construct the discrete level set and determine which cells are intersected. Note that all cells, intersected and non-intersected, have a corresponding active_fe_indicator. Similarly, the exact level set function is approximated on the whole triangulation and postprocessed afterward resulting on a surface approximation with no gaps.

  template <int dim>
  {
  pcout << "Creating background mesh..."
  << "\n"
  << std::flush;
  const double cube_side = 2.008901281;
  triangulation.refine_global(3);
  }
  template <int dim>
  {
  << "Setting up discrete level set function and reclassifying cells... "
  << "\n"
  << std::flush;
  TimerOutput::Scope t(computing_timer, "setup_level_set");
  active_fe_indicator.reinit(triangulation.n_active_cells());
constexpr void clear()
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)

Here is where the geometric information enters the code. Next, using the discrete level set, we mark the cell which are intersected by its zero contour. Finally, once the triangulation's cells are classified, we determine which cells are active.

  VectorType tmp_sol(level_set_dof_handler.locally_owned_dofs(),
  mpi_communicator);
  level_set_dof_handler.locally_owned_dofs(),
  mpi_communicator);
  level_set = tmp_sol;
  mesh_classifier.reclassify();
  for (const auto &cell : dof_handler.active_cell_iterators() |
  IteratorFilters::LocallyOwnedCell())
  {
  if (mesh_classifier.location_to_level_set(cell) ==
  NonMatching::LocationToLevelSet::intersected)
  cell->set_active_fe_index(
  else
  cell->set_active_fe_index(
  static_cast<types::fe_index>(ActiveFEIndex::nothing));
  }
  }
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={}, const unsigned int level=numbers::invalid_unsigned_int)

The method fills in the indicator telling which cells are intersected. It is used in the adaptive refinement near the surface.

  template <int dim>
  {
  pcout << "Determining cells with active FE index..."
  << "\n"
  << std::flush;
  for (const auto &cell : dof_handler.active_cell_iterators() |
  IteratorFilters::LocallyOwnedCell())
  {
  if (mesh_classifier.location_to_level_set(cell) ==
  NonMatching::LocationToLevelSet::intersected)
  active_fe_indicator[cell->active_cell_index()] = 1.0;
  }
  }

We refine only intersected cells with active_fe_indicator=1. We are calling GridRefinement::refine_and_coarsen_fixed_fraction() instead of the GridRefinement::refine_and_coarsen_fixed_number() function called in most other tutorials because the number of non-intersected cells also grows interfering with the number of active, intersected cells.

  template <int dim>
  {
  pcout << "Refining near surface..."
  << "\n"
  << std::flush;
  triangulation.execute_coarsening_and_refinement();
  }
void refine_and_coarsen_fixed_fraction(::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_error, const double bottom_fraction_of_error, const VectorTools::NormType norm_type=VectorTools::L1_norm)

As the surface is properly approximated by several adaptive steps, we may now distribute the degrees of freedom on cells which are intersected by the discrete approximation. Next, we initialize matrices for active DoFs and apply the constraints for the solution.

  template <int dim>
  {
  pcout << "Distributing degrees of freedom... "
  << "\n"
  << std::flush;
  dof_handler.distribute_dofs(fe_collection);
  locally_owned_dofs = dof_handler.locally_owned_dofs();
  completely_distributed_solution.reinit(dof_handler.locally_owned_dofs(),
  mpi_communicator);
  locally_relevant_solution.reinit(locally_owned_dofs,
  mpi_communicator);
  global_rhs.reinit(locally_owned_dofs, mpi_communicator);
  const unsigned int dof_handler_size = dof_handler.n_dofs();
  const unsigned int level_set_dof_handler_size =
  convergence_table.add_value("LevelSet dofs", level_set_dof_handler_size);
  convergence_table.evaluate_convergence_rates(
  convergence_table.add_value("Active dofs", dof_handler_size);
  convergence_table.evaluate_convergence_rates(
  }
  template <int dim>
  {
  pcout << "Initializing the matrix... "
  << "\n"
  << std::flush;
  DynamicSparsityPattern dsp(dof_handler.n_dofs(),
  dof_handler.n_dofs(),
  constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
  DoFTools::make_hanging_node_constraints(dof_handler, constraints);
  constraints.close();
  DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
  locally_owned_dofs,
  mpi_communicator,
  global_matrix.reinit(locally_owned_dofs,
  locally_owned_dofs,
  mpi_communicator);
  }
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)

Assembly and surface accumulation

We use a MeshWorker to assemble the linear problem efficiently. This cell worker does not do anything for non-intersected cells.

  template <int dim>
  {
  pcout << "Assembling... "
  << "\n"
  << std::flush;
  const auto cell_worker =
  [&](const typename DoFHandler<dim>::active_cell_iterator &cell,
  ScratchData<dim> &scratch_data,
  CopyData<dim> &copy_data) {
  if (mesh_classifier.location_to_level_set(cell) ==
  cell->is_locally_owned())
  {

Once we know that the cell is intersected, we construct the unfitted quadratures for the solutions FE space on the cell.

  scratch_data.non_matching_fe_values.reinit(cell);
  copy_data.reinit(cell);
  copy_data.cell_matrix = 0;
  copy_data.cell_rhs = 0;
  const std::optional<NonMatching::FEImmersedSurfaceValues<dim>>
  scratch_data.non_matching_fe_values.get_surface_fe_values();
  const std::vector<double> &JxW_surface =
  surface_fe_values->get_JxW_values();

The accumulation of the surface integrals, including the forcing, is performed here.

  for (unsigned int q : surface_fe_values->quadrature_point_indices())
  {
  const Tensor<1, dim> &normal =
  surface_fe_values->normal_vector(q);
  for (const unsigned int i : surface_fe_values->dof_indices())
  {
  copy_data.cell_rhs(i) +=
  surface_fe_values->shape_value(i, q) *
  surface_fe_values->quadrature_point(q)) *
  for (const unsigned int j :
  surface_fe_values->dof_indices())
  {
  copy_data.cell_matrix(i, j) +=
  (surface_fe_values->shape_value(i, q) *
  surface_fe_values->shape_value(j, q)) *
  copy_data.cell_matrix(i, j) +=
  (surface_fe_values->shape_grad(i, q) -
  (normal * surface_fe_values->shape_grad(i, q)) *
  normal) *
  (surface_fe_values->shape_grad(j, q) -
  (normal * surface_fe_values->shape_grad(j, q)) *
  normal) *
  }
  }
  }

The normal-gradient volume stabilization form needs a bulk cell integration while other types of stabilization may need face quadratures, for example. So we check it first. The cell was provided by the solution's DoFHandler, so we recast it as a level set's DoFHandler cell. However, it is the same geometric entity of the common triangulation.

  if (stabilization_scheme.needs_cell_worker())
  {
  cell->as_dof_handler_iterator(level_set_dof_handler);
  scratch_data.fe_values.reinit(cell);
  scratch_data.level_set_fe_values.reinit(level_set_cell);
  stabilization_scheme.assemble_cell_worker(level_set,
  cell,
  scratch_data,
  copy_data);
  }
  }
  };

Next, the copier worker distributes the local contributions from the CopyData taking into account the constraints. Finally, the MeshWorker goes over all cells provided by the solutions' DoFHandler. Note that this includes non-intersected cells as well, but the cell worker does nothing on them.

  const auto copier = [&](const CopyData<dim> &c) {
  constraints.distribute_local_to_global(c.cell_matrix,
  c.cell_rhs,
  c.local_dof_indices,
  };
  ScratchData<dim> scratch_data(mapping,
  fe_collection,
  mesh_classifier,
  level_set,
  cell_quadrature,
  CopyData<dim> copy_data;
  MeshWorker::mesh_loop(dof_handler.begin_active(),
  dof_handler.end(),
  copier,
  scratch_data,
  copy_data,
  }
void mesh_loop(const CellIteratorType &begin, const CellIteratorType &end, const CellWorkerFunctionType &cell_worker, const CopierType &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const BoundaryWorkerFunctionType &boundary_worker=BoundaryWorkerFunctionType(), const FaceWorkerFunctionType &face_worker=FaceWorkerFunctionType(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition mesh_loop.h:281

In the following, we solve the resulting linear system of equations. We either use a direct solver or AMG.

  template <int dim>
  {
  bool apply_direct_solver = false;
  const double relative_error = 1e-9 * global_rhs.l2_norm();
  unsigned int n_iterations = 0;
  {
  pcout << "Solving directly... " << '\n' << std::flush;
  SolverControl solver_control(100, relative_error);
  }
  else
  {
  Timer timer;
  pcout << "Solving with AMG... "
  << "\n"
  << std::flush;
  const unsigned int max_iterations = 500;
  SolverControl solver_control(max_iterations, relative_error);
  const std::vector<std::vector<bool>> constant_modes =
  Amg_data.constant_modes = constant_modes;
  Amg_data.elliptic = true;
  Amg_data.higher_order_elements = false;
  Amg_data.smoother_sweeps = 2;
  Amg_data.aggregation_threshold = 0.02;
  Amg_data.output_details = true;
  SolverCG<VectorType> cg(solver_control);
  cg.solve(global_matrix,
  n_iterations = solver_control.last_step();
  }
  constraints.distribute(completely_distributed_solution);
  convergence_table.add_value("Iterations", n_iterations);
  }
Definition timer.h:117
std::vector< index_type > data
Definition mpi.cc:740
std::vector< std::vector< bool > > extract_constant_modes(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask={})

Similarly to what we do in the assembly() function, a MeshWorker is used to accumulate errors including the stabilization term. At the end, we collect the results, and print them out.

  template <int dim>
  {
  pcout << "Evaluating errors on the surface..."
  << "\n"
  << std::flush;
  double error_L2_sqr = 0.0;
  double error_H1_sqr = 0.0;
  double error_stab_sqr = 0.0;
  const auto cell_worker = [&](const auto &cell,
  auto &scratch_data,
  auto &copy_data) {
  if (mesh_classifier.location_to_level_set(cell) ==
  cell->is_locally_owned())
  {
  double cell_L2_error_sqr = 0.0;
  double cell_H1_error_sqr = 0.0;
  copy_data.reinit(cell);
  scratch_data.non_matching_fe_values.reinit(cell);
  const std::optional<NonMatching::FEImmersedSurfaceValues<dim>>
  scratch_data.non_matching_fe_values.get_surface_fe_values();
  const std::vector<double> &JxW_surface =
  surface_fe_values->get_JxW_values();
  const unsigned int n_q_points =
  surface_fe_values->n_quadrature_points;
  std::vector<double> sol(n_q_points);
  std::vector<Tensor<1, dim>> sol_grad(n_q_points);
  for (const unsigned int q :
  surface_fe_values->quadrature_point_indices())
  {
  const Point<dim> &point = surface_fe_values->quadrature_point(q);
  const Tensor<1, dim> &normal =
  surface_fe_values->normal_vector(q);
  const double error_at_point =
  sol.at(q) - analytical_solution.value(point);
  (sol_grad.at(q) - (normal * sol_grad.at(q)) * normal -
  analytical_solution.gradient(point));
  }
  copy_data.cell_L2_error_sqr = cell_L2_error_sqr;
  copy_data.cell_H1_error_sqr = cell_H1_error_sqr;
  if (stabilization_scheme.needs_cell_worker())
  {
  cell->as_dof_handler_iterator(level_set_dof_handler);
  scratch_data.fe_values.reinit(cell);
  scratch_data.level_set_fe_values.reinit(level_set_cell);
  stabilization_scheme.evaluate_cell_worker(
  level_set,
  cell,
  scratch_data,
  copy_data);
  }
  }
  };
  const auto copier = [&](const auto &copy_data) {
  if (copy_data.cell_index < active_fe_indicator.size())
  {
  error_L2_sqr += copy_data.cell_L2_error_sqr;
  error_H1_sqr += copy_data.cell_H1_error_sqr;
  error_stab_sqr += copy_data.cell_stab_sqr;
  }
  };
  ScratchData<dim> scratch_data(mapping,
  fe_collection,
  mesh_classifier,
  level_set,
  cell_quadrature,
  MeshWorker::mesh_loop(dof_handler.begin_active(),
  dof_handler.end(),
  copier,
  scratch_data,
  copy_data,
  const double error_L2 =
  const double error_semiH1 =
  const double error_stab =
  convergence_table.add_value("L2 Error", error_L2);
  convergence_table.evaluate_convergence_rates(
  convergence_table.set_scientific("L2 Error", true);
  convergence_table.add_value("H1 error", error_semiH1);
  convergence_table.evaluate_convergence_rates(
  convergence_table.set_scientific("H1 error", true);
  convergence_table.add_value("Stab norm", error_stab);
  convergence_table.evaluate_convergence_rates(
  convergence_table.set_scientific("Stab norm", true);
  }
Definition point.h:113
T sum(const T &t, const MPI_Comm mpi_communicator)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)

The following two methods perform VTK output of the preliminary mesh refinements for the geometry approximation and of the TraceFEM solution. The important difference between the two is that the non-intersected cells are excluded from the output saving considerable amount of time and storage.

  template <int dim>
  void LaplaceBeltramiSolver<dim>::output_level_set(const unsigned int cycle)
  {
  pcout << "Writing vtu file for surface... " << '\n' << std::flush;
  TimerOutput::Scope t(computing_timer, "output_level_set");
  DataOut<dim> data_out;
  data_out.add_data_vector(level_set_dof_handler, level_set, "level_set");
  data_out.add_data_vector(active_fe_indicator, "ref_indicator");
  data_out.build_patches();
  data_out.write_vtu_in_parallel("surface_" + std::to_string(cycle) + ".vtu",
  mpi_communicator);
  }
  template <int dim>
  {
  pcout << "Writing vtu file... " << std::flush;
  TimerOutput::Scope t(computing_timer, "output_solution");
  Vector<double> exact(dof_handler.locally_owned_dofs().size());
  DataOut<dim> data_out;
  data_out.add_data_vector(dof_handler,
  "solution");
  data_out.add_data_vector(dof_handler, exact, "exact");
  data_out.add_data_vector(level_set_dof_handler, level_set, "level_set");
  data_out.set_cell_selection(
  [this](const typename Triangulation<dim>::cell_iterator &cell) {
  return cell->is_active() && cell->is_locally_owned() &&
  mesh_classifier.location_to_level_set(cell) ==
  });
  data_out.build_patches();
  data_out.write_vtu_in_parallel("solution.vtu", mpi_communicator);
  }

The method localize_surface() generates iteratively a surface approximation as described above. Once the surface approximation is constructed, the main logic of the solver is executed as presented in the method run().

  template <int dim>
  {
  unsigned int preliminary_levels = 3;
  for (unsigned int localization_cycle = 0;
  {
  pcout << std::endl
  << "Preliminary refinement #" << localization_cycle << std::endl;
  refine_grid();
  }
  }
  template <int dim>
  {
  const unsigned int convergence_levels = 3;
  for (unsigned int cycle = 0; cycle < convergence_levels; ++cycle)
  {
  pcout << std::endl << "Convergence refinement #" << cycle << std::endl;
  distribute_dofs();
  initialize_matrices();
  solve();
  if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
  convergence_table.write_text(pcout.get_stream());
  computing_timer.print_summary();
  if (cycle < convergence_levels - 1)
  {
  refine_grid();
  }
  else
  computing_timer.print_summary();
  }
  }
  } // namespace Step90
  int main(int argc, char *argv[])
  {
  try
  {
  using namespace dealii;
  using namespace Step90;
  LB_solver.run();
  }
  catch (std::exception &exc)
  {
  std::cerr << std::endl
  << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Exception on processing: " << std::endl
  << exc.what() << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
  catch (...)
  {
  std::cerr << std::endl
  << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Unknown exception!" << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
  return 0;
  }
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:112

Results

The numerical solution \(u_h\) for a very fine mesh \(\Gamma_h\) is shown below by plotting in Paraview the zero contour of the approximate level set \(\psi_h\) and restricting the discrete solution \(u_h\) to the resulting surface approximation \(\Gamma_h\).

Next, we demonstrate the corresponding set of intersected cells with active degrees of freedom. Note that not all cells are of the same refinement level which is attributed to the insufficiently fine initial uniform grid.

Convergence test

The results of the convergence study are shown in the following table. The experimental orders of convergence (EOC) are reported for the surface errors and the stabilization.

Cycle DOFS Rate Iterations \(L^2\)-Error EOC \(H^1\)-Error EOC \(s_h^{1/2}(u_h)\) EOC
0 12370 - 15 7.6322e-02 - 3.6212e-01 - 2.2423e-01 -
1 49406 2.00 18 1.1950e-02 2.68 1.4752e-01 1.30 1.1238e-01 1.00
2 196848 1.99 19 1.7306e-03 2.79 7.4723e-02 0.98 6.1131e-02 0.88
3 785351 2.00 22 3.6276e-04 2.25 3.9329e-02 0.93 3.0185e-02 1.02
4 3136501 2.00 25 7.5910e-05 2.26 1.9694e-02 1.00 1.4875e-02 1.02
5 12536006 2.00 26 1.7279e-05 2.14 9.8443e-03 1.00 7.4067e-03 1.01
6 50122218 2.00 30 4.3891e-06 1.98 4.9219e-03 1.00 3.7042e-03 1.00

In this test we refine the mesh near the surface and, as a result, the number of degrees of freedom scales in the two-dimensional fashion. The optimal rates of error convergence in \(L^2(\Gamma)\) and \(H^1(\Gamma)\) norms are clearly observable. We also note the first order convergence of the stabilization \(s_h^{1/2}(u_h)=\sqrt{s_h(u_h, u_h)}\) evaluated at the solution \(u_h\).

Parallel scalability

The weak and strong scalability test results are shown in the following figure. Clearly, the refine() method is responsible for the certain lack of parallel scalability.

The plain program

/* ------------------------------------------------------------------------
*
* SPDX-License-Identifier: LGPL-2.1-or-later
* Copyright (C) 2024 by the deal.II authors
*
* This file is part of the deal.II library.
*
* Part of the source code is dual licensed under Apache-2.0 WITH
* LLVM-exception OR LGPL-2.1-or-later. Detailed license information
* governing the source code and code contributions can be found in
* LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
*
* ------------------------------------------------------------------------
*
* This program was contributed by Vladimir Yushutin and Timo Heister, Clemson
* University, 2023.
*/
namespace Step90
{
using namespace dealii;
using VectorType = TrilinosWrappers::MPI::Vector;
using MatrixType = TrilinosWrappers::SparseMatrix;
{
lagrange = 0,
nothing = 1
};
template <int dim>
class TamarindShape : public Function<dim>
{
public:
: Function<dim>()
{}
double value(const Point<dim> &point,
const unsigned int component = 0) const override
{
AssertIndexRange(component, this->n_components);
(void)component;
Assert(dim == 3, ExcNotImplemented());
return 0.25 * Utilities::pow(point[0], 2) + Utilities::pow(point[1], 2) +
4.0 * Utilities::pow(point[2], 2) *
std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -2) -
1.0;
}
const unsigned int component = 0) const override
{
AssertIndexRange(component, this->n_components);
(void)component;
Assert(dim == 3, ExcNotImplemented());
grad[0] = 0.5 * point[0] +
(-2.0) * 4.0 * Utilities::pow(point[2], 2) *
std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -3) *
(0.5 * numbers::PI * std::cos(numbers::PI * point[0]));
grad[1] = 2.0 * point[1];
grad[2] = (2.0) * 4.0 * point[2] *
std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -2);
return grad;
}
hessian(const Point<dim> &point,
const unsigned int component = 0) const override
{
AssertIndexRange(component, this->n_components);
(void)component;
Assert(dim == 3, ExcNotImplemented());
hessian[0][0] =
0.5 +
8.0 * Utilities::pow(point[2], 2) *
(3.0 * std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -4) *
std::cos(numbers::PI * point[0]),
2) +
std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -3) * 0.5 *
hessian[0][1] = 0.0;
hessian[0][2] =
(-8.0) * point[2] *
std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -3) *
hessian[1][1] = 2.0;
hessian[1][2] = 0.0;
hessian[2][2] =
8.0 * std::pow(1.0 + 0.5 * std::sin(numbers::PI * point[0]), -2);
return hessian;
}
};
template <int dim>
class AnalyticalSolution : public Function<dim>
{
private:
public:
: Function<dim>()
{}
double value(const Point<dim> &point,
const unsigned int component = 0) const override;
const unsigned int component = 0) const override;
};
template <int dim>
const unsigned int component) const
{
AssertIndexRange(component, this->n_components);
(void)component;
return point[0] * point[1];
}
template <int dim>
const unsigned int component) const
{
AssertIndexRange(component, this->n_components);
(void)component;
const Tensor<1, dim> grad = tamarind.gradient(point, component);
const Tensor<1, dim> normal = grad / grad.norm();
Tensor<1, dim> projector_first_column = -normal[0] * normal;
Tensor<1, dim> projector_second_column = -normal[1] * normal;
}
template <int dim>
class RightHandSide : public Function<dim>
{
public:
: Function<dim>()
{}
virtual double value(const Point<dim> &p,
const unsigned int component = 0) const override;
};
template <int dim>
const unsigned int component) const
{
AssertIndexRange(component, this->n_components);
(void)component;
Assert(dim == 3, ExcNotImplemented());
const Tensor<1, dim> grad = tamarind.gradient(point, component);
const Tensor<1, dim> normal = grad / grad.norm();
const SymmetricTensor<2, dim> hessian = tamarind.hessian(point, component);
double mean_curv = 0.0;
for (int j = 0; j < 3; j++)
for (int k = 0; k < 3; k++)
mean_curv += ((j == k ? 1 : 0) - normal[j] * normal[k]) * hessian[j][k];
return point[0] * point[1] + 2.0 * normal[0] * normal[1] +
mean_curv * (point[1] * normal[0] + point[0] * normal[1]);
}
template <int dim>
struct ScratchData
{
ScratchData(const Mapping<dim> &mapping,
const hp::FECollection<dim> &fe_collection,
const NonMatching::MeshClassifier<dim> &mesh_classifier,
const VectorType &level_set,
const Quadrature<dim> &quadrature,
const UpdateFlags cell_update_flags = update_values |
: fe_values(
mapping,
quadrature,
cell_update_flags)
, region_update_flags(nonmatching_update_flags)
, fe_collection(fe_collection)
, mesh_classifier(mesh_classifier)
, level_set(level_set)
quadrature,
cell_update_flags)
, non_matching_fe_values(fe_collection,
mesh_classifier,
level_set)
{}
ScratchData(const ScratchData<dim> &scratch_data)
: fe_values(scratch_data.fe_values.get_mapping(),
scratch_data.fe_values.get_fe(),
scratch_data.fe_values.get_quadrature(),
scratch_data.fe_values.get_update_flags())
, region_update_flags(scratch_data.region_update_flags)
, quadrature_1D(scratch_data.quadrature_1D)
, fe_collection(scratch_data.fe_collection)
, mesh_classifier(scratch_data.mesh_classifier)
, level_set(scratch_data.level_set)
, level_set_fe_values(scratch_data.level_set_fe_values.get_mapping(),
scratch_data.level_set_fe_values.get_fe(),
scratch_data.level_set_fe_values.get_quadrature(),
scratch_data.level_set_fe_values.get_update_flags())
, non_matching_fe_values(fe_collection,
region_update_flags,
mesh_classifier,
level_set)
{}
FEValues<dim> fe_values;
const NonMatching::RegionUpdateFlags region_update_flags;
const hp::FECollection<dim> &fe_collection;
const NonMatching::MeshClassifier<dim> &mesh_classifier;
const VectorType &level_set;
};
template <int dim>
struct CopyData
{
std::vector<types::global_dof_index> local_dof_indices;
void reinit(const typename DoFHandler<dim>::active_cell_iterator &cell)
{
const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
cell_rhs.reinit(dofs_per_cell);
local_dof_indices.resize(dofs_per_cell);
cell->get_dof_indices(local_dof_indices);
}
};
template <int dim>
{
unsigned int cell_index;
double cell_stab_sqr;
void reinit(const typename DoFHandler<dim>::active_cell_iterator &cell)
{
cell_index = cell->active_cell_index();
}
};
template <int dim>
{
public:
{}
bool needs_cell_worker() const
{
return true;
}
VectorType &level_set,
ScratchData<dim> &scratch_data,
CopyData<dim> &copy_data) const
{
const FEValues<dim> &fe_values = scratch_data.fe_values;
scratch_data.level_set_fe_values;
const std::vector<double> &JxW_cell = fe_values.get_JxW_values();
std::vector<Tensor<1, dim>> grad_level_set(
level_set_fe_values.get_quadrature().size());
level_set_fe_values.get_function_gradients(level_set, grad_level_set);
const double factor =
std::pow(cell->minimum_vertex_distance(), stabilization_exponent);
for (const unsigned int q : fe_values.quadrature_point_indices())
{
const Tensor<1, dim> &normal =
for (const unsigned int i : fe_values.dof_indices())
for (const unsigned int j : fe_values.dof_indices())
copy_data.cell_matrix(i, j) +=
factor * (normal * fe_values.shape_grad(i, q)) *
(normal * fe_values.shape_grad(j, q)) * JxW_cell[q];
}
}
VectorType &solution,
VectorType &level_set,
ScratchData<dim> &scratch_data,
CopyDataError<dim> &copy_data) const
{
double cell_stab_sqr = 0.0;
const FEValues<dim> &fe_values = scratch_data.fe_values;
const std::vector<double> &JxW_cell = fe_values.get_JxW_values();
const unsigned int n_q_points = fe_values.get_quadrature_points().size();
scratch_data.level_set_fe_values;
std::vector<Tensor<1, dim>> level_set_grad(n_q_points);
level_set_fe_values.get_function_gradients(level_set, level_set_grad);
std::vector<Tensor<1, dim>> sol_grad(n_q_points);
fe_values.get_function_gradients(solution, sol_grad);
const double factor =
std::pow(cell->minimum_vertex_distance(), stabilization_exponent);
for (const unsigned int q : fe_values.quadrature_point_indices())
{
const Tensor<1, dim> normal =
const double stabilization_at_point = normal * sol_grad[q];
}
copy_data.cell_stab_sqr = cell_stab_sqr;
}
private:
const double stabilization_exponent;
};
template <int dim>
{
public:
void run();
private:
void make_grid();
void distribute_dofs();
void initialize_matrices();
void solve();
void refine_grid();
void output_level_set(unsigned int);
MPI_Comm mpi_communicator;
const unsigned int fe_degree;
hp::FECollection<dim> fe_collection;
DoFHandler<dim> dof_handler;
const unsigned int level_set_fe_degree;
const MappingQ1<dim> mapping;
VectorType level_set;
const QGauss<dim> cell_quadrature;
SparsityPattern sparsity_pattern;
IndexSet locally_owned_dofs;
};
template <int dim>
: mpi_communicator(MPI_COMM_WORLD)
, triangulation(mpi_communicator)
(Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
, computing_timer(mpi_communicator,
TimerOutput::never,
TimerOutput::wall_times)
, fe_degree(1)
, fe_collection(FE_Q<dim>(fe_degree), FE_Nothing<dim>())
, dof_handler(triangulation)
, mapping()
, mesh_classifier(level_set_dof_handler, level_set)
, cell_quadrature(fe_degree + 1)
, quadrature_1D(fe_degree + 1)
{
}
template <int dim>
{
pcout << "Creating background mesh..."
<< "\n"
<< std::flush;
const double cube_side = 2.008901281;
triangulation.refine_global(3);
}
template <int dim>
{
<< "Setting up discrete level set function and reclassifying cells... "
<< "\n"
<< std::flush;
TimerOutput::Scope t(computing_timer, "setup_level_set");
active_fe_indicator.reinit(triangulation.n_active_cells());
mpi_communicator);
level_set_dof_handler.locally_owned_dofs(),
mpi_communicator);
level_set = tmp_sol;
mesh_classifier.reclassify();
for (const auto &cell : dof_handler.active_cell_iterators() |
IteratorFilters::LocallyOwnedCell())
{
if (mesh_classifier.location_to_level_set(cell) ==
cell->set_active_fe_index(
else
cell->set_active_fe_index(
static_cast<types::fe_index>(ActiveFEIndex::nothing));
}
}
template <int dim>
{
pcout << "Determining cells with active FE index..."
<< "\n"
<< std::flush;
for (const auto &cell : dof_handler.active_cell_iterators() |
IteratorFilters::LocallyOwnedCell())
{
if (mesh_classifier.location_to_level_set(cell) ==
active_fe_indicator[cell->active_cell_index()] = 1.0;
}
}
template <int dim>
{
pcout << "Refining near surface..."
<< "\n"
<< std::flush;
triangulation.execute_coarsening_and_refinement();
}
template <int dim>
{
pcout << "Distributing degrees of freedom... "
<< "\n"
<< std::flush;
dof_handler.distribute_dofs(fe_collection);
locally_owned_dofs = dof_handler.locally_owned_dofs();
completely_distributed_solution.reinit(dof_handler.locally_owned_dofs(),
mpi_communicator);
locally_relevant_solution.reinit(locally_owned_dofs,
mpi_communicator);
global_rhs.reinit(locally_owned_dofs, mpi_communicator);
const unsigned int dof_handler_size = dof_handler.n_dofs();
const unsigned int level_set_dof_handler_size =
convergence_table.add_value("LevelSet dofs", level_set_dof_handler_size);
convergence_table.evaluate_convergence_rates(
convergence_table.add_value("Active dofs", dof_handler_size);
convergence_table.evaluate_convergence_rates(
}
template <int dim>
{
pcout << "Initializing the matrix... "
<< "\n"
<< std::flush;
DynamicSparsityPattern dsp(dof_handler.n_dofs(),
dof_handler.n_dofs(),
constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
constraints.close();
DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
locally_owned_dofs,
mpi_communicator,
global_matrix.reinit(locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
}
template <int dim>
{
pcout << "Assembling... "
<< "\n"
<< std::flush;
const auto cell_worker =
[&](const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData<dim> &scratch_data,
CopyData<dim> &copy_data) {
if (mesh_classifier.location_to_level_set(cell) ==
cell->is_locally_owned())
{
scratch_data.non_matching_fe_values.reinit(cell);
copy_data.reinit(cell);
copy_data.cell_matrix = 0;
copy_data.cell_rhs = 0;
const std::optional<NonMatching::FEImmersedSurfaceValues<dim>>
scratch_data.non_matching_fe_values.get_surface_fe_values();
const std::vector<double> &JxW_surface =
surface_fe_values->get_JxW_values();
for (unsigned int q : surface_fe_values->quadrature_point_indices())
{
const Tensor<1, dim> &normal =
surface_fe_values->normal_vector(q);
for (const unsigned int i : surface_fe_values->dof_indices())
{
copy_data.cell_rhs(i) +=
surface_fe_values->shape_value(i, q) *
surface_fe_values->quadrature_point(q)) *
for (const unsigned int j :
surface_fe_values->dof_indices())
{
copy_data.cell_matrix(i, j) +=
(surface_fe_values->shape_value(i, q) *
surface_fe_values->shape_value(j, q)) *
copy_data.cell_matrix(i, j) +=
(surface_fe_values->shape_grad(i, q) -
(normal * surface_fe_values->shape_grad(i, q)) *
normal) *
(surface_fe_values->shape_grad(j, q) -
(normal * surface_fe_values->shape_grad(j, q)) *
normal) *
}
}
}
if (stabilization_scheme.needs_cell_worker())
{
cell->as_dof_handler_iterator(level_set_dof_handler);
scratch_data.fe_values.reinit(cell);
scratch_data.level_set_fe_values.reinit(level_set_cell);
stabilization_scheme.assemble_cell_worker(level_set,
cell,
scratch_data,
copy_data);
}
}
};
const auto copier = [&](const CopyData<dim> &c) {
constraints.distribute_local_to_global(c.cell_matrix,
c.cell_rhs,
c.local_dof_indices,
};
ScratchData<dim> scratch_data(mapping,
fe_collection,
mesh_classifier,
level_set,
cell_quadrature,
CopyData<dim> copy_data;
MeshWorker::mesh_loop(dof_handler.begin_active(),
dof_handler.end(),
copier,
scratch_data,
copy_data,
}
template <int dim>
{
bool apply_direct_solver = false;
const double relative_error = 1e-9 * global_rhs.l2_norm();
unsigned int n_iterations = 0;
{
pcout << "Solving directly... " << '\n' << std::flush;
SolverControl solver_control(100, relative_error);
}
else
{
Timer timer;
pcout << "Solving with AMG... "
<< "\n"
<< std::flush;
const unsigned int max_iterations = 500;
SolverControl solver_control(max_iterations, relative_error);
const std::vector<std::vector<bool>> constant_modes =
Amg_data.constant_modes = constant_modes;
Amg_data.elliptic = true;
Amg_data.higher_order_elements = false;
Amg_data.smoother_sweeps = 2;
Amg_data.aggregation_threshold = 0.02;
Amg_data.output_details = true;
SolverCG<VectorType> cg(solver_control);
cg.solve(global_matrix,
n_iterations = solver_control.last_step();
}
constraints.distribute(completely_distributed_solution);
convergence_table.add_value("Iterations", n_iterations);
}
template <int dim>
{
pcout << "Evaluating errors on the surface..."
<< "\n"
<< std::flush;
double error_L2_sqr = 0.0;
double error_H1_sqr = 0.0;
double error_stab_sqr = 0.0;
const auto cell_worker = [&](const auto &cell,
auto &scratch_data,
auto &copy_data) {
if (mesh_classifier.location_to_level_set(cell) ==
cell->is_locally_owned())
{
double cell_L2_error_sqr = 0.0;
double cell_H1_error_sqr = 0.0;
copy_data.reinit(cell);
scratch_data.non_matching_fe_values.reinit(cell);
const std::optional<NonMatching::FEImmersedSurfaceValues<dim>>
scratch_data.non_matching_fe_values.get_surface_fe_values();
const std::vector<double> &JxW_surface =
surface_fe_values->get_JxW_values();
const unsigned int n_q_points =
surface_fe_values->n_quadrature_points;
std::vector<double> sol(n_q_points);
sol);
std::vector<Tensor<1, dim>> sol_grad(n_q_points);
for (const unsigned int q :
surface_fe_values->quadrature_point_indices())
{
const Point<dim> &point = surface_fe_values->quadrature_point(q);
const Tensor<1, dim> &normal =
surface_fe_values->normal_vector(q);
const double error_at_point =
sol.at(q) - analytical_solution.value(point);
(sol_grad.at(q) - (normal * sol_grad.at(q)) * normal -
analytical_solution.gradient(point));
}
copy_data.cell_L2_error_sqr = cell_L2_error_sqr;
copy_data.cell_H1_error_sqr = cell_H1_error_sqr;
if (stabilization_scheme.needs_cell_worker())
{
cell->as_dof_handler_iterator(level_set_dof_handler);
scratch_data.fe_values.reinit(cell);
scratch_data.level_set_fe_values.reinit(level_set_cell);
stabilization_scheme.evaluate_cell_worker(
level_set,
cell,
scratch_data,
copy_data);
}
}
};
const auto copier = [&](const auto &copy_data) {
if (copy_data.cell_index < active_fe_indicator.size())
{
error_L2_sqr += copy_data.cell_L2_error_sqr;
error_H1_sqr += copy_data.cell_H1_error_sqr;
error_stab_sqr += copy_data.cell_stab_sqr;
}
};
ScratchData<dim> scratch_data(mapping,
fe_collection,
mesh_classifier,
level_set,
cell_quadrature,
CopyDataError<dim> copy_data;
MeshWorker::mesh_loop(dof_handler.begin_active(),
dof_handler.end(),
copier,
scratch_data,
copy_data,
const double error_L2 =
const double error_semiH1 =
const double error_stab =
convergence_table.add_value("L2 Error", error_L2);
convergence_table.evaluate_convergence_rates(
convergence_table.set_scientific("L2 Error", true);
convergence_table.add_value("H1 error", error_semiH1);
convergence_table.evaluate_convergence_rates(
convergence_table.set_scientific("H1 error", true);
convergence_table.add_value("Stab norm", error_stab);
convergence_table.evaluate_convergence_rates(
convergence_table.set_scientific("Stab norm", true);
}
template <int dim>
void LaplaceBeltramiSolver<dim>::output_level_set(const unsigned int cycle)
{
pcout << "Writing vtu file for surface... " << '\n' << std::flush;
TimerOutput::Scope t(computing_timer, "output_level_set");
DataOut<dim> data_out;
data_out.add_data_vector(level_set_dof_handler, level_set, "level_set");
data_out.add_data_vector(active_fe_indicator, "ref_indicator");
data_out.build_patches();
data_out.write_vtu_in_parallel("surface_" + std::to_string(cycle) + ".vtu",
mpi_communicator);
}
template <int dim>
{
pcout << "Writing vtu file... " << std::flush;
TimerOutput::Scope t(computing_timer, "output_solution");
Vector<double> exact(dof_handler.locally_owned_dofs().size());
DataOut<dim> data_out;
data_out.add_data_vector(dof_handler,
"solution");
data_out.add_data_vector(dof_handler, exact, "exact");
data_out.add_data_vector(level_set_dof_handler, level_set, "level_set");
data_out.set_cell_selection(
[this](const typename Triangulation<dim>::cell_iterator &cell) {
return cell->is_active() && cell->is_locally_owned() &&
mesh_classifier.location_to_level_set(cell) ==
});
data_out.build_patches();
data_out.write_vtu_in_parallel("solution.vtu", mpi_communicator);
}
template <int dim>
{
unsigned int preliminary_levels = 3;
for (unsigned int localization_cycle = 0;
{
pcout << std::endl
<< "Preliminary refinement #" << localization_cycle << std::endl;
refine_grid();
}
computing_timer.reset();
}
template <int dim>
{
const unsigned int convergence_levels = 3;
for (unsigned int cycle = 0; cycle < convergence_levels; ++cycle)
{
pcout << std::endl << "Convergence refinement #" << cycle << std::endl;
distribute_dofs();
initialize_matrices();
solve();
if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
convergence_table.write_text(pcout.get_stream());
computing_timer.print_summary();
computing_timer.reset();
if (cycle < convergence_levels - 1)
{
refine_grid();
}
else
computing_timer.print_summary();
computing_timer.reset();
}
}
} // namespace Step90
int main(int argc, char *argv[])
{
try
{
using namespace dealii;
using namespace Step90;
LB_solver.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
void reinit(const Triangulation< dim, spacedim > &tria)
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
unsigned short int fe_index
Definition types.h:68