This tutorial depends on step-85.
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/convergence_table.h>
#include <deal.II/base/function.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/quadrature_lib.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_update_flags.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/grid_generator.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/dynamic_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/lac/trilinos_sparse_matrix.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_immersed_values.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
{
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.
{
lagrange = 0,
nothing = 1
};
unsigned short int fe_index
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:
TamarindShape()
{}
const unsigned int component = 0) const override
{
(void)component;
1.0;
}
const unsigned int component = 0) const override
{
(void)component;
grad[0] = 0.5 *
point[0] +
grad[1] = 2.0 *
point[1];
grad[2] = (2.0) * 4.0 * point[2] *
return grad;
}
const unsigned int component = 0) const override
{
(void)component;
0.5 +
2) +
}
};
virtual SymmetricTensor< 2, dim, RangeNumberType > hessian(const Point< dim > &p, const unsigned int component=0) const
virtual Tensor< 1, dim, RangeNumberType > gradient(const Point< dim > &p, const unsigned int component=0) const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
constexpr T pow(const T base, const int iexp)
static constexpr double PI
::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:
const TamarindShape<dim> tamarind;
public:
AnalyticalSolution()
{}
const unsigned int component = 0) const override;
const unsigned int component = 0) const override;
};
template <int dim>
double AnalyticalSolution<dim>::value(
const Point<dim> &point,
const unsigned int component) const
{
(void)component;
}
template <int dim>
AnalyticalSolution<dim>::gradient(
const Point<dim> &point,
const unsigned int component) const
{
(void)component;
projector_first_column[0] += 1.0;
projector_second_column[1] += 1.0;
point[1] * projector_first_column +
point[0] * projector_second_column;
return surface_gradient;
}
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>
{
const TamarindShape<dim> tamarind;
public:
RightHandSide()
{}
const unsigned int component = 0) const override;
};
template <int dim>
double RightHandSide<dim>::value(
const Point<dim> &point,
const unsigned int component) const
{
(void)component;
Assert(dim == 3, ExcNotImplemented());
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];
mean_curv /= grad.norm();
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
{
const VectorType &level_set,
: fe_values(
mapping,
fe_collection[static_cast<
types::fe_index>(ActiveFEIndex::lagrange)],
quadrature,
cell_update_flags)
, region_update_flags(nonmatching_update_flags)
, quadrature_1D(quadrature_edge)
, fe_collection(fe_collection)
, mesh_classifier(mesh_classifier)
, level_set_dof_handler(level_set_dof_handler)
, level_set(level_set)
, level_set_fe_values(mapping,
level_set_dof_handler.get_fe(),
quadrature,
cell_update_flags)
, non_matching_fe_values(fe_collection,
quadrature_edge,
nonmatching_update_flags,
mesh_classifier,
level_set_dof_handler,
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_dof_handler(scratch_data.level_set_dof_handler)
, 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,
quadrature_1D,
region_update_flags,
mesh_classifier,
level_set_dof_handler,
level_set)
{}
Abstract base class for mapping classes.
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
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.
const VectorType &level_set;
};
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
{
std::vector<types::global_dof_index> local_dof_indices;
{
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>
struct CopyDataError
{
double cell_L2_error_sqr;
double cell_H1_error_sqr;
double cell_stab_sqr;
{
cell_L2_error_sqr = 0.0;
cell_H1_error_sqr = 0.0;
cell_stab_sqr = 0.0;
}
};
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
unsigned int n_dofs_per_cell() const
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\).
template <int dim>
class NormalGradientVolumeStabilization
{
public:
NormalGradientVolumeStabilization()
: stabilization_parameter(1.0)
, stabilization_exponent(-1.0)
{}
bool needs_cell_worker() const
{
return true;
}
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.
void assemble_cell_worker(
VectorType &level_set,
ScratchData<dim> &scratch_data,
CopyData<dim> ©_data) const
{
scratch_data.level_set_fe_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 =
stabilization_parameter *
std::pow(cell->minimum_vertex_distance(), stabilization_exponent);
for (const unsigned int q : fe_values.quadrature_point_indices())
{
const
Tensor<1, dim> &normal =
grad_level_set[q] / grad_level_set[q].norm();
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];
}
}
void evaluate_cell_worker(
VectorType &solution,
VectorType &level_set,
ScratchData<dim> &scratch_data,
CopyDataError<dim> ©_data) const
{
double cell_stab_sqr = 0.0;
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 =
stabilization_parameter *
std::pow(cell->minimum_vertex_distance(), stabilization_exponent);
for (const unsigned int q : fe_values.quadrature_point_indices())
{
level_set_grad[q] / level_set_grad[q].
norm();
const double stabilization_at_point = normal * sol_grad[q];
cell_stab_sqr +=
}
copy_data.cell_stab_sqr = cell_stab_sqr;
}
private:
const double stabilization_parameter;
const double stabilization_exponent;
};
const std::vector< double > & get_JxW_values() const
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
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>
class LaplaceBeltramiSolver
{
public:
LaplaceBeltramiSolver();
void run();
private:
void make_grid();
void localize_surface();
void setup_discrete_level_set();
void distribute_dofs();
void initialize_matrices();
void assemble_system();
void solve();
void mark_intersected();
void refine_grid();
void compute_errors();
void output_level_set(unsigned int);
void output_solution();
The surface of interest corresponds to the zero contour of the following exact level set function:
const TamarindShape<dim> tamarind;
The manufactured solution to the Laplace–Beltrami problem and the corresponding right-hand side:
const AnalyticalSolution<dim> analytical_solution;
const RightHandSide<dim> right_hand_side;
There is a single triangulation which is shared by the discretizations of the solution and of the level set.
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
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;
const unsigned int level_set_fe_degree;
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.
VectorType completely_distributed_solution;
VectorType locally_relevant_solution;
VectorType locally_relevant_exact;
VectorType level_set;
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.
Any TraceFEM needs a stabilization, and we choose the normal-gradient, volume stabilization.
const NormalGradientVolumeStabilization<dim> stabilization_scheme;
Discrete right-hand side and the final matrix corresponding to dof_handler.
VectorType global_rhs;
MatrixType global_matrix;
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>
LaplaceBeltramiSolver<dim>::LaplaceBeltramiSolver()
: mpi_communicator(MPI_COMM_WORLD)
, tamarind()
, analytical_solution()
, right_hand_side()
, computing_timer(mpi_communicator,
pcout,
, fe_degree(1)
, level_set_fe_degree(1)
, level_set_fe(level_set_fe_degree)
, mapping()
, mesh_classifier(level_set_dof_handler, level_set)
, cell_quadrature(fe_degree + 1)
, quadrature_1D(fe_degree + 1)
, stabilization_scheme()
{
surface_update_flags.surface =
}
@ update_normal_vectors
Normal vectors.
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>
void LaplaceBeltramiSolver<dim>::make_grid()
{
pcout << "Creating background mesh..."
<< "\n"
<< std::flush;
const double cube_side = 2.008901281;
}
template <int dim>
void LaplaceBeltramiSolver<dim>::setup_discrete_level_set()
{
pcout
<< "Setting up discrete level set function and reclassifying cells... "
<< "\n"
<< std::flush;
level_set_dof_handler.distribute_dofs(level_set_fe);
level_set_constraints.clear();
const IndexSet level_set_locally_relevant_dofs =
level_set_constraints.reinit(level_set_locally_relevant_dofs);
level_set_constraints);
level_set_constraints.close();
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
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_constraints.distribute(tmp_sol);
level_set.reinit(level_set_locally_relevant_dofs,
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() |
{
if (mesh_classifier.location_to_level_set(cell) ==
cell->set_active_fe_index(
static_cast<
types::fe_index>(ActiveFEIndex::lagrange));
else
cell->set_active_fe_index(
}
}
The method fills in the indicator telling which cells are intersected. It is used in the adaptive refinement near the surface.
template <int dim>
void LaplaceBeltramiSolver<dim>::mark_intersected()
{
pcout << "Determining cells with active FE index..."
<< "\n"
<< std::flush;
for (const auto &cell : dof_handler.active_cell_iterators() |
{
if (mesh_classifier.location_to_level_set(cell) ==
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>
void LaplaceBeltramiSolver<dim>::refine_grid()
{
pcout << "Refining near surface..."
<< "\n"
<< std::flush;
}
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>
void LaplaceBeltramiSolver<dim>::distribute_dofs()
{
pcout << "Distributing degrees of freedom... "
<< "\n"
<< std::flush;
dof_handler.distribute_dofs(fe_collection);
locally_owned_dofs = dof_handler.locally_owned_dofs();
locally_relevant_dofs =
completely_distributed_solution.reinit(dof_handler.locally_owned_dofs(),
mpi_communicator);
locally_relevant_solution.reinit(locally_owned_dofs,
locally_relevant_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 =
level_set_dof_handler.n_dofs();
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>
void LaplaceBeltramiSolver<dim>::initialize_matrices()
{
pcout << "Initializing the matrix... "
<< "\n"
<< std::flush;
dof_handler.n_dofs(),
locally_relevant_dofs);
constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
constraints.close();
locally_owned_dofs,
mpi_communicator,
locally_relevant_dofs);
global_matrix.reinit(locally_owned_dofs,
locally_owned_dofs,
dsp,
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)
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>
void LaplaceBeltramiSolver<dim>::assemble_system()
{
pcout << "Assembling... "
<< "\n"
<< std::flush;
const auto cell_worker =
ScratchData<dim> &scratch_data,
CopyData<dim> ©_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>>
&surface_fe_values =
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) *
right_hand_side.value(
surface_fe_values->quadrature_point(q)) *
JxW_surface[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)) *
JxW_surface[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) *
JxW_surface[q];
}
}
}
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,
global_matrix,
global_rhs);
};
ScratchData<dim> scratch_data(mapping,
fe_collection,
mesh_classifier,
level_set_dof_handler,
level_set,
surface_update_flags,
cell_quadrature,
quadrature_1D);
CopyData<dim> copy_data;
dof_handler.end(),
cell_worker,
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)
In the following, we solve the resulting linear system of equations. We either use a direct solver or AMG.
template <int dim>
void LaplaceBeltramiSolver<dim>::solve()
{
bool apply_direct_solver = false;
const double relative_error = 1e-9 * global_rhs.l2_norm();
unsigned int n_iterations = 0;
if (apply_direct_solver)
{
pcout << "Solving directly... " << '\n' << std::flush;
trilinos.solve(global_matrix,
completely_distributed_solution,
global_rhs);
}
else
{
pcout << "Solving with AMG... "
<< "\n"
<< std::flush;
const unsigned int max_iterations = 500;
const std::vector<std::vector<bool>> 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;
preconditioner_stiffness.initialize(global_matrix);
cg.solve(global_matrix,
completely_distributed_solution,
global_rhs,
preconditioner_stiffness);
n_iterations = solver_control.last_step();
}
constraints.distribute(completely_distributed_solution);
locally_relevant_solution = completely_distributed_solution;
convergence_table.add_value("Iterations", n_iterations);
}
std::vector< index_type > data
std::vector< std::vector< bool > > constant_modes
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>
void LaplaceBeltramiSolver<dim>::compute_errors()
{
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 ©_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>>
&surface_fe_values =
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);
surface_fe_values->get_function_values(locally_relevant_solution,
sol);
std::vector<Tensor<1, dim>> sol_grad(n_q_points);
surface_fe_values->get_function_gradients(locally_relevant_solution,
sol_grad);
for (const unsigned int q :
surface_fe_values->quadrature_point_indices())
{
const
Point<dim> &point = surface_fe_values->quadrature_point(q);
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));
cell_L2_error_sqr +=
cell_H1_error_sqr +=
grad_error_at_point * grad_error_at_point * JxW_surface[q];
}
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(
locally_relevant_solution,
level_set,
cell,
scratch_data,
copy_data);
}
}
};
const auto copier = [&](const auto ©_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_dof_handler,
level_set,
surface_update_flags,
cell_quadrature,
quadrature_1D);
CopyDataError<dim> copy_data;
dof_handler.end(),
cell_worker,
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);
}
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;
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>
void LaplaceBeltramiSolver<dim>::output_solution()
{
pcout << "Writing vtu file... " << std::flush;
locally_relevant_solution,
"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(
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);
}
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
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>
void LaplaceBeltramiSolver<dim>::localize_surface()
{
unsigned int preliminary_levels = 3;
for (unsigned int localization_cycle = 0;
localization_cycle < preliminary_levels;
++localization_cycle)
{
pcout << std::endl
<< "Preliminary refinement #" << localization_cycle << std::endl;
setup_discrete_level_set();
mark_intersected();
output_level_set(localization_cycle);
refine_grid();
}
computing_timer.reset();
}
template <int dim>
void LaplaceBeltramiSolver<dim>::run()
{
make_grid();
localize_surface();
const unsigned int convergence_levels = 3;
for (unsigned int cycle = 0; cycle < convergence_levels; ++cycle)
{
pcout << std::endl << "Convergence refinement #" << cycle << std::endl;
setup_discrete_level_set();
distribute_dofs();
initialize_matrices();
assemble_system();
solve();
compute_errors();
convergence_table.write_text(pcout.get_stream());
computing_timer.print_summary();
computing_timer.reset();
if (cycle < convergence_levels - 1)
{
mark_intersected();
refine_grid();
}
else
output_solution();
computing_timer.print_summary();
computing_timer.reset();
}
}
}
int main(int argc, char *argv[])
{
try
{
using namespace Step90;
LaplaceBeltramiSolver<3> LB_solver;
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)
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
namespace Step90
{
{
lagrange = 0,
nothing = 1
};
template <int dim>
class TamarindShape :
public Function<dim>
{
public:
TamarindShape()
{}
const unsigned int component = 0) const override
{
(void)component;
Assert(dim == 3, ExcNotImplemented());
1.0;
}
const unsigned int component = 0) const override
{
(void)component;
Assert(dim == 3, ExcNotImplemented());
grad[0] = 0.5 *
point[0] +
grad[1] = 2.0 *
point[1];
grad[2] = (2.0) * 4.0 * point[2] *
return grad;
}
const unsigned int component = 0) const override
{
(void)component;
Assert(dim == 3, ExcNotImplemented());
0.5 +
2) +
(-8.0) * point[2] *
}
};
template <int dim>
class AnalyticalSolution :
public Function<dim>
{
private:
const TamarindShape<dim> tamarind;
public:
AnalyticalSolution()
{}
const unsigned int component = 0) const override;
const unsigned int component = 0) const override;
};
template <int dim>
double AnalyticalSolution<dim>::value(
const Point<dim> &point,
const unsigned int component) const
{
(void)component;
}
template <int dim>
AnalyticalSolution<dim>::gradient(
const Point<dim> &point,
const unsigned int component) const
{
(void)component;
projector_first_column[0] += 1.0;
projector_second_column[1] += 1.0;
point[1] * projector_first_column +
point[0] * projector_second_column;
return surface_gradient;
}
template <int dim>
class RightHandSide :
public Function<dim>
{
const TamarindShape<dim> tamarind;
public:
RightHandSide()
{}
const unsigned int component = 0) const override;
};
template <int dim>
double RightHandSide<dim>::value(
const Point<dim> &point,
const unsigned int component) const
{
(void)component;
Assert(dim == 3, ExcNotImplemented());
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];
mean_curv /= grad.
norm();
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
{
const VectorType &level_set,
: fe_values(
mapping,
fe_collection[static_cast<
types::
fe_index>(ActiveFEIndex::lagrange)],
quadrature,
cell_update_flags)
, region_update_flags(nonmatching_update_flags)
, quadrature_1D(quadrature_edge)
, fe_collection(fe_collection)
, mesh_classifier(mesh_classifier)
, level_set_dof_handler(level_set_dof_handler)
, level_set(level_set)
, level_set_fe_values(mapping,
level_set_dof_handler.get_fe(),
quadrature,
cell_update_flags)
, non_matching_fe_values(fe_collection,
quadrature_edge,
nonmatching_update_flags,
mesh_classifier,
level_set_dof_handler,
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_dof_handler(scratch_data.level_set_dof_handler)
, 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,
quadrature_1D,
region_update_flags,
mesh_classifier,
level_set_dof_handler,
level_set)
{}
};
template <int dim>
struct CopyData
{
std::vector<types::global_dof_index> local_dof_indices;
{
cell_rhs.
reinit(dofs_per_cell);
local_dof_indices.resize(dofs_per_cell);
cell->get_dof_indices(local_dof_indices);
}
};
template <int dim>
struct CopyDataError
{
double cell_L2_error_sqr;
double cell_H1_error_sqr;
double cell_stab_sqr;
{
cell_L2_error_sqr = 0.0;
cell_H1_error_sqr = 0.0;
cell_stab_sqr = 0.0;
}
};
template <int dim>
class NormalGradientVolumeStabilization
{
public:
NormalGradientVolumeStabilization()
: stabilization_parameter(1.0)
, stabilization_exponent(-1.0)
{}
bool needs_cell_worker() const
{
return true;
}
void assemble_cell_worker(
VectorType &level_set,
ScratchData<dim> &scratch_data,
CopyData<dim> ©_data) const
{
scratch_data.level_set_fe_values;
std::vector<Tensor<1, dim>> grad_level_set(
const double factor =
stabilization_parameter *
std::pow(cell->minimum_vertex_distance(), stabilization_exponent);
for (const unsigned int q : fe_values.quadrature_point_indices())
{
grad_level_set[q] / grad_level_set[q].
norm();
for (const unsigned int i : fe_values.dof_indices())
for (const unsigned
int j : fe_values.dof_indices())
factor * (normal * fe_values.shape_grad(i, q)) *
(normal * fe_values.shape_grad(j, q)) * JxW_cell[q];
}
}
void evaluate_cell_worker(
VectorType &solution,
VectorType &level_set,
ScratchData<dim> &scratch_data,
CopyDataError<dim> ©_data) const
{
double cell_stab_sqr = 0.0;
scratch_data.level_set_fe_values;
std::vector<Tensor<1, dim>> level_set_grad(n_q_points);
std::vector<Tensor<1, dim>> sol_grad(n_q_points);
const double factor =
stabilization_parameter *
std::pow(cell->minimum_vertex_distance(), stabilization_exponent);
for (const unsigned int q : fe_values.quadrature_point_indices())
{
level_set_grad[q] / level_set_grad[q].
norm();
const double stabilization_at_point = normal * sol_grad[q];
cell_stab_sqr +=
}
copy_data.cell_stab_sqr = cell_stab_sqr;
}
private:
const double stabilization_parameter;
const double stabilization_exponent;
};
template <int dim>
class LaplaceBeltramiSolver
{
public:
LaplaceBeltramiSolver();
private:
void make_grid();
void localize_surface();
void setup_discrete_level_set();
void distribute_dofs();
void initialize_matrices();
void assemble_system();
void solve();
void mark_intersected();
void refine_grid();
void compute_errors();
void output_level_set(unsigned int);
void output_solution();
const TamarindShape<dim> tamarind;
const AnalyticalSolution<dim> analytical_solution;
const RightHandSide<dim> right_hand_side;
const unsigned int fe_degree;
const unsigned int level_set_fe_degree;
const NormalGradientVolumeStabilization<dim> stabilization_scheme;
};
template <int dim>
LaplaceBeltramiSolver<dim>::LaplaceBeltramiSolver()
: mpi_communicator(MPI_COMM_WORLD)
, tamarind()
, analytical_solution()
, right_hand_side()
, computing_timer(mpi_communicator,
pcout,
, fe_degree(1)
, level_set_fe_degree(1)
, level_set_fe(level_set_fe_degree)
, mapping()
, mesh_classifier(level_set_dof_handler, level_set)
, cell_quadrature(fe_degree + 1)
, quadrature_1D(fe_degree + 1)
, stabilization_scheme()
{
surface_update_flags.surface =
}
template <int dim>
void LaplaceBeltramiSolver<dim>::make_grid()
{
pcout << "Creating background mesh..."
<< "\n"
<< std::flush;
const double cube_side = 2.008901281;
}
template <int dim>
void LaplaceBeltramiSolver<dim>::setup_discrete_level_set()
{
pcout
<< "Setting up discrete level set function and reclassifying cells... "
<< "\n"
<< std::flush;
level_set_dof_handler.distribute_dofs(level_set_fe);
level_set_constraints.clear();
const IndexSet level_set_locally_relevant_dofs =
level_set_constraints.reinit(level_set_locally_relevant_dofs);
level_set_constraints);
level_set_constraints.close();
VectorType tmp_sol(level_set_dof_handler.locally_owned_dofs(),
mpi_communicator);
level_set_constraints.distribute(tmp_sol);
level_set.reinit(level_set_locally_relevant_dofs,
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() |
{
if (mesh_classifier.location_to_level_set(cell) ==
cell->set_active_fe_index(
else
cell->set_active_fe_index(
}
}
template <int dim>
void LaplaceBeltramiSolver<dim>::mark_intersected()
{
pcout << "Determining cells with active FE index..."
<< "\n"
<< std::flush;
for (const auto &cell : dof_handler.active_cell_iterators() |
{
if (mesh_classifier.location_to_level_set(cell) ==
active_fe_indicator[cell->active_cell_index()] = 1.0;
}
}
template <int dim>
void LaplaceBeltramiSolver<dim>::refine_grid()
{
pcout << "Refining near surface..."
<< "\n"
<< std::flush;
}
template <int dim>
void LaplaceBeltramiSolver<dim>::distribute_dofs()
{
pcout << "Distributing degrees of freedom... "
<< "\n"
<< std::flush;
dof_handler.distribute_dofs(fe_collection);
locally_owned_dofs = dof_handler.locally_owned_dofs();
locally_relevant_dofs =
completely_distributed_solution.reinit(dof_handler.locally_owned_dofs(),
mpi_communicator);
locally_relevant_solution.reinit(locally_owned_dofs,
locally_relevant_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 =
level_set_dof_handler.n_dofs();
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>
void LaplaceBeltramiSolver<dim>::initialize_matrices()
{
pcout << "Initializing the matrix... "
<< "\n"
<< std::flush;
dof_handler.n_dofs(),
locally_relevant_dofs);
constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
constraints.close();
locally_owned_dofs,
mpi_communicator,
locally_relevant_dofs);
global_matrix.reinit(locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
}
template <int dim>
void LaplaceBeltramiSolver<dim>::assemble_system()
{
pcout << "Assembling... "
<< "\n"
<< std::flush;
const auto cell_worker =
ScratchData<dim> &scratch_data,
CopyData<dim> ©_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>>
&surface_fe_values =
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())
{
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) *
right_hand_side.value(
surface_fe_values->quadrature_point(q)) *
JxW_surface[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)) *
JxW_surface[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) *
JxW_surface[q];
}
}
}
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,
global_matrix,
global_rhs);
};
ScratchData<dim> scratch_data(mapping,
fe_collection,
mesh_classifier,
level_set_dof_handler,
level_set,
surface_update_flags,
cell_quadrature,
quadrature_1D);
CopyData<dim> copy_data;
dof_handler.end(),
cell_worker,
copier,
scratch_data,
copy_data,
}
template <int dim>
void LaplaceBeltramiSolver<dim>::solve()
{
bool apply_direct_solver = false;
const double relative_error = 1
e-9 * global_rhs.l2_norm();
unsigned int n_iterations = 0;
if (apply_direct_solver)
{
pcout << "Solving directly... " << '\n' << std::flush;
trilinos.solve(global_matrix,
completely_distributed_solution,
global_rhs);
}
else
{
pcout << "Solving with AMG... "
<< "\n"
<< std::flush;
const unsigned int max_iterations = 500;
const std::vector<std::vector<bool>> constant_modes =
preconditioner_stiffness.
initialize(global_matrix);
cg.solve(global_matrix,
completely_distributed_solution,
global_rhs,
preconditioner_stiffness);
n_iterations = solver_control.last_step();
}
constraints.distribute(completely_distributed_solution);
locally_relevant_solution = completely_distributed_solution;
convergence_table.add_value("Iterations", n_iterations);
}
template <int dim>
void LaplaceBeltramiSolver<dim>::compute_errors()
{
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 ©_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>>
&surface_fe_values =
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);
surface_fe_values->get_function_values(locally_relevant_solution,
sol);
std::vector<Tensor<1, dim>> sol_grad(n_q_points);
surface_fe_values->get_function_gradients(locally_relevant_solution,
sol_grad);
for (const unsigned int q :
surface_fe_values->quadrature_point_indices())
{
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));
cell_L2_error_sqr +=
cell_H1_error_sqr +=
grad_error_at_point * grad_error_at_point * JxW_surface[q];
}
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(
locally_relevant_solution,
level_set,
cell,
scratch_data,
copy_data);
}
}
};
const auto copier = [&](const auto ©_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_dof_handler,
level_set,
surface_update_flags,
cell_quadrature,
quadrature_1D);
CopyDataError<dim> copy_data;
dof_handler.end(),
cell_worker,
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;
mpi_communicator);
}
template <int dim>
void LaplaceBeltramiSolver<dim>::output_solution()
{
pcout << "Writing vtu file... " << std::flush;
locally_relevant_solution,
"solution");
return cell->is_active() && cell->is_locally_owned() &&
mesh_classifier.location_to_level_set(cell) ==
});
}
template <int dim>
void LaplaceBeltramiSolver<dim>::localize_surface()
{
unsigned int preliminary_levels = 3;
for (unsigned int localization_cycle = 0;
localization_cycle < preliminary_levels;
++localization_cycle)
{
pcout << std::endl
<< "Preliminary refinement #" << localization_cycle << std::endl;
setup_discrete_level_set();
mark_intersected();
output_level_set(localization_cycle);
refine_grid();
}
computing_timer.reset();
}
template <int dim>
void LaplaceBeltramiSolver<dim>::run()
{
make_grid();
localize_surface();
const unsigned int convergence_levels = 3;
for (unsigned int cycle = 0; cycle < convergence_levels; ++cycle)
{
pcout << std::endl << "Convergence refinement #" << cycle << std::endl;
setup_discrete_level_set();
distribute_dofs();
initialize_matrices();
assemble_system();
solve();
compute_errors();
convergence_table.write_text(pcout.get_stream());
computing_timer.print_summary();
computing_timer.reset();
if (cycle < convergence_levels - 1)
{
mark_intersected();
refine_grid();
}
else
output_solution();
computing_timer.print_summary();
computing_timer.reset();
}
}
}
int main(int argc, char *argv[])
{
try
{
using namespace Step90;
LaplaceBeltramiSolver<3> LB_solver;
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 write_vtu_in_parallel(const std::string &filename, const MPI_Comm comm) const
virtual void build_patches(const unsigned int n_subdivisions=0)
void set_cell_selection(const std::function< cell_iterator(const Triangulation< dim, spacedim > &)> &first_cell, const std::function< cell_iterator(const Triangulation< dim, spacedim > &, const cell_iterator &)> &next_cell)
void reinit(const Triangulation< dim, spacedim > &tria)
const std::vector< Point< spacedim > > & get_quadrature_points() const
void get_function_gradients(const ReadVector< Number > &fe_function, std::vector< Tensor< 1, spacedim, Number > > &gradients) const
const Quadrature< dim > & get_quadrature() const
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
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.)
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)
double aggregation_threshold
unsigned int smoother_sweeps
bool higher_order_elements
1.9.8