In this example, we show the use of the matrix-free framework for the cut finite element method (CutFEM) and the cut discontinuous Galerkin method (CutDG) in deal.II. This tutorial is based on step-85. Please read the introduction there for details about the level-set discretization, quadrature generation (with NonMatching::QuadratureGenerator) and mesh classification (with NonMatching::MeshClassifier) before starting with this tutorial. We choose the same problem like step-85, Poisson's equation:
\begin{align*}
-\Delta u &= f \qquad && \text{in }\, \Omega,
\\
u &= u_D \qquad && \text{on }\, \Gamma_D = \partial \Omega,
\end{align*}
where we choose \(f(x) = 4\) and \(u_D(x) = 0\).
With CutFEM/CutDG this leads to the following weak formulation: Find \(u_h \in V_\Omega^h\) such that
\begin{equation*}
a_h(u_h, v_h) = L_h(v_h), \quad \forall v_h \in V_\Omega^h,
\end{equation*}
\begin{align*}
a_{h,CG}(u_h, v_h) &= (\nabla u_h, \nabla v_h)_\Omega
- (u_h, \partial_n v_h)_{\Gamma_D}
- \left (\partial_n u_h - \gamma_D h^{-1} u_h, v_h \right )_{\Gamma_D}
+ g_h(u_h, v_h),
\\
L_h(v_h) &= (f,v)_\Omega
+ \left (u_D, \gamma_D h^{-1} v_h -\partial_n v_h \right )_{\Gamma_D},
\end{align*}
where \(g_h\) is the ghost penalty stabilization that circumvents the small cut cell problem.
For the discontinuous Galerkin variant (CutDG) we get an additional SIPG flux term on interior faces:
\begin{equation*}
a_{h,DG}(u_h, v_h) = a_{h,CG}(u_h, v_h)
- (\jump{u_h}, \average{\partial_n v_h})_{\Gamma_{int}} -
(\average{\partial_n u_h} - \gamma_{SIPG}h^{-1}\jump{u_h}, \jump{v_h})_{\Gamma_{int}}
\end{equation*}
\begin{equation*}
g_h(u_h,v_h) = \sum_{F \in \mathcal{F}_h} \sum_{j=0}^p \gamma_{GP,j} \left(h_F^{2j-1}[\partial_n^j u_h], [\partial_n^j v_h] \right)_F
\end{equation*}
over faces \(F\in\mathcal{F}_h\) between intersected cells and intersected and inside cells. In this tutorial it is implemented for degree \(p\leq 2\).
\begin{equation*}
y=A(u)=\sum_{cell=1}^{n_{cells}}A_{cell}(u)+\sum_{face=1}^{n_{faces}}A_{face}(u)
\end{equation*}
The first include files have all been treated in previous examples.
And most important the header for flexible matrix-free evaluation.
We define helper functions to determine if a cell (or a batch of cells) is inside or intersected depending on the active_fe_index.
We define 3 types frequently used throughout the tutorial keep the definition of the Number type at a central place.
We also define types for the matrix-free evaluation classes for cells and faces and for structured and unstructured quadrature to keep the definition of the template arguments for those at a central place.
Finally, we introduce a shortcut for the width of the detected SIMD vectorization.
For the SIPG flux, the Nitsche term, and the ghost-penalty terms, we need a scaling by a characteristic element length. To quickly access these values, we precompute the characteristic element length and store it vectorized in a vector.
This function is the interface function for linear solvers that applies the operator evaluation, which is an abstraction of a matrix-vector product by executing loops over cells and faces.
The right-hand-side vector is assembled with a loop over the cells.
This functions partitions a vector with matrix-free functionality.
With this function, we can assemble the diagonal of the operator instead of applying the matrix-vector product. It uses the same function as the operator evaluation for the local operation on cells and faces but sets the template argument assemble to true.
We define two helper functions needed for matrix-free assembly operations.
This function implements the local cell operation.
We define the cell_operation for inside cells as the standard Poisson term.
This function implements the local face operation.
We start with the face operations for the DG case. We define the face_operation for an inside face as the SIPG term.
For the CG case we only need to define the ghost penalty term for mixed and intersected faces.
For cases where we need to call evaluate twice on the same cell DoFs, we introduce an option to save the DoF values to a buffer such that we can reuse them. This way, we don't have to read them from the global vector again.
We do not need a local operation of the fitted boundary for this tutorial as the whole Dirichlet boundary is immersed in the volume of the domain. However, this tutorial is easily extendable to the case with fitted and unfitted boundaries.
This is the actual implementation of the quadrature point operation of the Poisson term in the weak form. It is templated over the EvaluatorType to be usable by FEEvaluation as well as FEPointEvaluation.
The implementation of the Nitsche term.
The implementation of the face-based ghost penalty term (up to degree 2 / normal hessians). For continuous elements, we have value_m = value_p thanks to continuity, which allows us to skip this term.
The implementation of the right-hand-side term evaluating the rhs function (unfortunately, we cannot evaluate a Function object with vectorized types directly, so we have to reshuffle the quadrature point data).
Helper function which computes the interior penalty parameter for the SIPG flux.
We generate a background mesh with perfectly Cartesian cells. Our domain is a unit disc centered at the origin, so we need to make the background mesh a bit larger than \([-1, 1]^{\text{dim}}\) to completely cover \(\Omega\).
The discrete level set function is defined on the whole background mesh. Thus, to set up the DoFHandler for the level set function, we distribute DoFs over all elements in \(\mathcal{T}_h\). We then set up the discrete level set function by interpolating onto this finite element space.
We set up the level set vector with all locally relevant DoFs. This is currently required by the NonMatching::MeshClassifier.
In case of DG, we also have to compute mapping data for cut faces, so we again fill the containers in matrix-free ordering.
Here we create the Jacobi preconditioner, which assembles the diagonal on construction. Then, we call the preconditioned conjugate gradient solver.
To test that the implementation works as expected, we want to compute the error in the solution in the \(L^2\)-norm. The analytical solution to the Poisson problem stated in the introduction reads
\begin{align*}
u(x) = - \frac{2}{\text{dim}}(\| x \|^2 - 1) , \qquad x \in
\overline{\Omega}.
\end{align*}
Of course, the analytical solution, and thus also the error, is only defined in \(\overline{\Omega}\). Thus, to compute the \(L^2\)-error we must proceed in the same way as for the operator.
Finally, we do a convergence study to check that the \(L^2\)-error decreases with the expected rate. We refine the background mesh a few times. In each refinement cycle, we solve the problem, compute the error, and add the \(L^2\)-error and the mesh size to a ConvergenceTable.
The results of the convergence studies are shown in the tables below.
We see that the \(L^2\) error decreases as we refine and that the estimated order of convergence, EOC, is close to 2.
If we increase the degree of the finite element for both the level-set and the solution we observe EOC close to 3.
We observe optimal order convergence also for DG.
This tutorial could be easily extended to time dependent problems like the heat equation by adding a mass operator and a time integration scheme in the solver class.
Another option would be to implement a different ghost penalty like the volume based ghost penalty to allow stabilization for arbitrary polynomial degree of the shape functions.
Also, the preconditioner could be extended to a more robust scheme like multigrid.
If you are interested in the mentioned possibilities for extensions we recommend our paper "High-performance matrix-free
unfitted finite element operator evaluation" [32] where volume ghost-penalty and multigrid preconditioning already have been addressed.
#include <fstream>
#include <vector>
namespace Step95
{
inline bool is_inside(const unsigned int active_fe_index)
{
return active_fe_index ==
}
inline bool is_intersected(const unsigned int active_fe_index)
{
return active_fe_index == static_cast<unsigned int>(
}
template <int dim, typename Number>
const unsigned int component = 0)
{
for (unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v)
{
for (
unsigned int d = 0;
d < dim; ++
d)
p[d] = p_vectorized[d][v];
result[v] = function.
value(p, component);
}
return result;
}
template <int dim>
class PoissonOperator
{
using CellEvaluator =
using FaceEvaluator =
using GenericCellEvaluator =
using GenericFaceEvaluator =
static constexpr unsigned int n_lanes = VectorizedArrayType::size();
public:
void
&mapping_info_cell_in,
&mapping_info_surface_in,
&mapping_info_faces_in,
const bool is_dg_in);
void vmult(VectorType &dst, const VectorType &src) const;
void initialize_dof_vector(VectorType &vec) const;
private:
template <typename EvaluatorType>
void create_zero_basis(EvaluatorType &evaluator) const;
template <typename EvaluatorType>
void create_standard_basis(const unsigned int j,
EvaluatorType &evaluator) const;
template <bool assemble = false>
void local_apply_cell(
VectorType &dst,
const VectorType &src,
const std::pair<unsigned int, unsigned int> &cell_range) const;
template <bool assemble = false>
void local_apply_face(
VectorType &dst,
const VectorType &src,
const std::pair<unsigned int, unsigned int> &face_range) const;
void local_apply_boundary_face(
VectorType &,
const VectorType &,
const std::pair<unsigned int, unsigned int> &) const;
template <typename EvaluatorType>
void do_poisson_cell_term(EvaluatorType &evaluator,
const unsigned int q) const;
template <typename EvaluatorType, typename Number2>
void do_flux_term(EvaluatorType &evaluator_m,
EvaluatorType &evaluator_p,
const Number2 &tau,
const unsigned int q) const;
template <typename EvaluatorType, typename Number2>
void do_boundary_flux_term_homogeneous(EvaluatorType &evaluator_m,
const Number2 &tau,
const unsigned int q) const;
template <bool do_normal_hessians, bool do_values, typename EvaluatorType>
void do_gp_face_term(
EvaluatorType &evaluator_m,
EvaluatorType &evaluator_p,
const typename EvaluatorType::NumberType &masked_factor_value,
const typename EvaluatorType::NumberType &masked_factor_gradient,
const typename EvaluatorType::NumberType &masked_factor_hessian,
const unsigned int q) const;
template <typename EvaluatorType>
void do_rhs_cell_term(EvaluatorType &evaluator,
const unsigned int q) const;
void do_local_apply_sipg_term(FaceEvaluator &evaluator_m,
FaceEvaluator &evaluator_p,
const VectorizedArrayType *dof_ptr_m,
const VectorizedArrayType *dof_ptr_p,
const VectorizedArrayType &diameter) const;
template <bool is_dg_>
void do_local_apply_gp_face_term(FaceEvaluator &evaluator_m,
FaceEvaluator &evaluator_p,
const VectorizedArrayType *dof_ptr_m,
const VectorizedArrayType *dof_ptr_p,
const VectorizedArrayType &diameter,
const bool sum_into_values) const;
bool
is_inside_face(std::pair<unsigned int, unsigned int> face_category) const;
bool
is_mixed_face(std::pair<unsigned int, unsigned int> face_category) const;
bool is_intersected_face(
std::pair<unsigned int, unsigned int> face_category) const;
VectorizedArrayType
compute_diameter_of_inner_face_batch(unsigned int face_batch_index) const;
VectorizedArrayType compute_interior_penalty_parameter(
const VectorizedArrayType &diameter) const;
matrix_free;
mapping_info_cell;
mapping_info_surface;
mapping_info_faces;
std::unique_ptr<GenericCellEvaluator> evaluator_cell;
std::unique_ptr<GenericCellEvaluator> evaluator_surface;
std::unique_ptr<GenericFaceEvaluator> evaluator_face_m;
std::unique_ptr<GenericFaceEvaluator> evaluator_face_p;
const unsigned int dof_index = 0;
const unsigned int quad_index = 0;
bool is_dg = false;
};
template <int dim>
void PoissonOperator<dim>::reinit(
&mapping_info_cell_in,
&mapping_info_surface_in,
&mapping_info_faces_in,
const bool is_dg_in)
{
is_dg = is_dg_in;
matrix_free = &matrix_free_in;
mapping_info_cell = &mapping_info_cell_in;
mapping_info_surface = &mapping_info_surface_in;
mapping_info_faces = is_dg ? (&mapping_info_faces_in) : nullptr;
const auto &fe = matrix_free->get_dof_handler().get_fe();
evaluator_cell =
std::make_unique<GenericCellEvaluator>(*mapping_info_cell, fe, 0, true);
evaluator_surface = std::make_unique<GenericCellEvaluator>(
*mapping_info_surface, fe, 0, true);
if (is_dg)
{
evaluator_face_m =
std::make_unique<GenericFaceEvaluator>(*mapping_info_faces, fe, true);
evaluator_face_p =
std::make_unique<GenericFaceEvaluator>(*mapping_info_faces,
fe,
false);
}
matrix_free->initialize_cell_data_vector(cell_diameter);
for (unsigned int cell_batch_index = 0;
cell_batch_index <
matrix_free->n_cell_batches() + matrix_free->n_ghost_cell_batches();
++cell_batch_index)
{
auto &
diameter = cell_diameter[cell_batch_index];
for (unsigned int v = 0;
v < matrix_free->n_active_entries_per_cell_batch(cell_batch_index);
++v)
{
const auto cell_accessor_inside =
matrix_free->get_cell_iterator(cell_batch_index, v);
diameter[v] = cell_accessor_inside->minimum_vertex_distance();
}
}
}
template <int dim>
void PoissonOperator<dim>::vmult(PoissonOperator::VectorType &dst,
const PoissonOperator::VectorType &src) const
{
matrix_free->loop(&PoissonOperator::local_apply_cell,
&PoissonOperator::local_apply_face,
&PoissonOperator::local_apply_boundary_face,
this,
dst,
src,
true);
}
template <int dim>
void PoissonOperator<dim>::rhs(PoissonOperator::VectorType &rhs,
{
const unsigned int dummy = 0;
matrix_free->template cell_loop<VectorType, unsigned int>(
VectorType &dst,
const VectorType &,
const std::pair<unsigned int, unsigned int> &cell_range) {
CellEvaluator evaluator(*matrix_free, dof_index, quad_index);
const auto cell_range_category =
matrix_free->get_cell_range_category(cell_range);
if (is_inside(cell_range_category))
{
for (unsigned int cell_batch_index = cell_range.first;
cell_batch_index < cell_range.second;
++cell_batch_index)
{
evaluator.reinit(cell_batch_index);
for (unsigned int q : evaluator.quadrature_point_indices())
do_rhs_cell_term(evaluator, rhs_function, q);
evaluator.distribute_local_to_global(dst);
}
}
else if (is_intersected(cell_range_category))
{
const auto dofs_per_cell = evaluator.dofs_per_cell;
for (unsigned int cell_batch_index = cell_range.first;
cell_batch_index < cell_range.second;
++cell_batch_index)
{
evaluator.reinit(cell_batch_index);
for (unsigned int v = 0; v < n_lanes; ++v)
{
cell_batch_index * n_lanes + v;
for (const auto q :
evaluator_cell->quadrature_point_indices())
do_rhs_cell_term(*evaluator_cell, rhs_function, q);
evaluator_cell->integrate(
&evaluator.begin_dof_values()[0][v], dofs_per_cell),
}
evaluator.distribute_local_to_global(dst);
}
}
},
rhs,
dummy);
}
template <int dim>
void PoissonOperator<dim>::initialize_dof_vector(
PoissonOperator::VectorType &vec) const
{
matrix_free->initialize_dof_vector(vec);
}
template <int dim>
void PoissonOperator<dim>::compute_diagonal(
PoissonOperator::VectorType &diagonal) const
{
matrix_free->loop(&PoissonOperator::local_apply_cell<true>,
&PoissonOperator::local_apply_face<true>,
&PoissonOperator::local_apply_boundary_face,
this,
diagonal,
diagonal);
}
template <int dim>
template <typename EvaluatorType>
void PoissonOperator<dim>::create_zero_basis(EvaluatorType &evaluator) const
{
for (unsigned int i = 0; i < evaluator.dofs_per_cell; ++i)
evaluator.begin_dof_values()[i] = VectorizedArrayType(0.);
}
template <int dim>
template <typename EvaluatorType>
void
PoissonOperator<dim>::create_standard_basis(const unsigned int j,
EvaluatorType &evaluator) const
{
create_zero_basis(evaluator);
evaluator.begin_dof_values()[j] = VectorizedArrayType(1.);
}
template <int dim>
template <bool assemble>
void PoissonOperator<dim>::local_apply_cell(
PoissonOperator::VectorType &dst,
const PoissonOperator::VectorType &src,
const std::pair<unsigned int, unsigned int> &cell_range) const
{
CellEvaluator evaluator(*matrix_free, dof_index, quad_index);
const auto dofs_per_cell = evaluator.dofs_per_cell;
const unsigned int cell_range_category =
matrix_free->get_cell_range_category(cell_range);
const auto inside_cell_operation = [&](CellEvaluator &evaluator) {
for (unsigned int q : evaluator.quadrature_point_indices())
do_poisson_cell_term(evaluator, q);
};
const auto intersected_cell_operation = [&](CellEvaluator &evaluator) {
const unsigned int cell_batch_index =
evaluator.get_cell_or_face_batch_id();
const VectorizedArrayType
diameter = cell_diameter[cell_batch_index];
const auto tau = compute_interior_penalty_parameter(diameter);
for (unsigned int v = 0; v < n_lanes; ++v)
{
const unsigned int cell_index = cell_batch_index * n_lanes + v;
&evaluator.begin_dof_values()[0][v],
dofs_per_cell),
&evaluator.begin_dof_values()[0][v],
dofs_per_cell),
for (const auto q : evaluator_cell->quadrature_point_indices())
do_poisson_cell_term(*evaluator_cell, q);
for (const auto q : evaluator_surface->quadrature_point_indices())
do_boundary_flux_term_homogeneous(*evaluator_surface, tau[v], q);
&evaluator.begin_dof_values()[0][v],
dofs_per_cell),
&evaluator.begin_dof_values()[0][v],
dofs_per_cell),
true);
}
};
std::function<void(CellEvaluator &)> cell_operation;
if (is_inside(cell_range_category))
cell_operation = inside_cell_operation;
else if (is_intersected(cell_range_category))
cell_operation = intersected_cell_operation;
else
return;
for (unsigned int cell_batch_index = cell_range.first;
cell_batch_index < cell_range.second;
++cell_batch_index)
{
evaluator.reinit(cell_batch_index);
if (assemble)
{
for (unsigned int j = 0; j < evaluator.dofs_per_cell; ++j)
{
create_standard_basis(j, evaluator);
cell_operation(evaluator);
dof_buffer[j] = evaluator.begin_dof_values()[j];
}
for (unsigned int j = 0; j < evaluator.dofs_per_cell; ++j)
evaluator.begin_dof_values()[j] = dof_buffer[j];
}
else
{
evaluator.read_dof_values(src);
cell_operation(evaluator);
}
evaluator.distribute_local_to_global(dst);
}
}
template <int dim>
template <bool assemble>
void PoissonOperator<dim>::local_apply_face(
PoissonOperator::VectorType &dst,
const PoissonOperator::VectorType &src,
const std::pair<unsigned int, unsigned int> &face_range) const
{
FaceEvaluator evaluator_m(*matrix_free, true, dof_index, quad_index);
FaceEvaluator evaluator_p(*matrix_free, false, dof_index, quad_index);
const std::pair<unsigned int, unsigned int> face_range_category =
matrix_free->get_face_range_category(face_range);
const bool need_buffer = is_dg && !is_inside_face(face_range_category);
need_buffer ? evaluator_m.dofs_per_cell : 0);
need_buffer ? evaluator_p.dofs_per_cell : 0);
assemble ? evaluator_m.dofs_per_cell : 0);
assemble ? evaluator_p.dofs_per_cell : 0);
const auto inside_face_operation_dg = [&](FaceEvaluator &evaluator_m,
FaceEvaluator &evaluator_p) {
const unsigned int face_batch_index =
evaluator_m.get_cell_or_face_batch_id();
compute_diameter_of_inner_face_batch(face_batch_index);
do_local_apply_sipg_term(evaluator_m,
evaluator_p,
evaluator_m.begin_dof_values(),
evaluator_p.begin_dof_values(),
diameter);
};
const auto mixed_face_operation_dg = [&](FaceEvaluator &evaluator_m,
FaceEvaluator &evaluator_p) {
const unsigned int face_batch_index =
evaluator_m.get_cell_or_face_batch_id();
compute_diameter_of_inner_face_batch(face_batch_index);
do_local_apply_sipg_term(evaluator_m,
evaluator_p,
dof_buffer_m.data(),
dof_buffer_p.data(),
diameter);
do_local_apply_gp_face_term<true>(evaluator_m,
evaluator_p,
dof_buffer_m.data(),
dof_buffer_p.data(),
diameter,
true);
};
const auto intersected_face_operation_dg = [&](FaceEvaluator &evaluator_m,
FaceEvaluator &evaluator_p) {
const unsigned int face_batch_index =
evaluator_m.get_cell_or_face_batch_id();
compute_diameter_of_inner_face_batch(face_batch_index);
const auto tau = compute_interior_penalty_parameter(diameter);
evaluator_m.project_to_face(dof_buffer_m.data(),
evaluator_p.project_to_face(dof_buffer_p.data(),
for (unsigned int v = 0; v < n_lanes; ++v)
{
const unsigned int face_index = face_batch_index * n_lanes + v;
evaluator_face_m->reinit(face_index);
evaluator_face_p->reinit(face_index);
evaluator_face_m->evaluate_in_face(
&evaluator_m.get_scratch_data().begin()[0][v],
evaluator_face_p->evaluate_in_face(
&evaluator_p.get_scratch_data().begin()[0][v],
for (const auto q : evaluator_face_m->quadrature_point_indices())
do_flux_term(*evaluator_face_m, *evaluator_face_p, tau[v], q);
evaluator_face_m->integrate_in_face(
&evaluator_m.get_scratch_data().begin()[0][v],
evaluator_face_p->integrate_in_face(
&evaluator_p.get_scratch_data().begin()[0][v],
}
do_local_apply_gp_face_term<true>(evaluator_m,
evaluator_p,
dof_buffer_m.data(),
dof_buffer_p.data(),
diameter,
true);
};
const auto intersected_mixed_face_operation_cg =
[&](FaceEvaluator &evaluator_m, FaceEvaluator &evaluator_p) {
const unsigned int face_batch_index =
evaluator_m.get_cell_or_face_batch_id();
compute_diameter_of_inner_face_batch(face_batch_index);
do_local_apply_gp_face_term<false>(evaluator_m,
evaluator_p,
evaluator_m.begin_dof_values(),
evaluator_p.begin_dof_values(),
diameter,
false);
};
std::function<void(FaceEvaluator &, FaceEvaluator &)> face_operation;
bool buffer_dof_values = false;
if (is_dg)
{
if (is_inside_face(face_range_category))
face_operation = inside_face_operation_dg;
else if (is_mixed_face(face_range_category))
{
face_operation = mixed_face_operation_dg;
buffer_dof_values = true;
}
else if (is_intersected_face(face_range_category))
{
face_operation = intersected_face_operation_dg;
buffer_dof_values = true;
}
else
return;
}
else
{
if (is_mixed_face(face_range_category) ||
is_intersected_face(face_range_category))
face_operation = intersected_mixed_face_operation_cg;
else
return;
}
for (unsigned int face_batch_index = face_range.first;
face_batch_index < face_range.second;
++face_batch_index)
{
evaluator_m.reinit(face_batch_index);
evaluator_p.reinit(face_batch_index);
if (assemble)
{
for (unsigned int j = 0; j < evaluator_m.dofs_per_cell; ++j)
for (unsigned int p = 0; p < 2; ++p)
{
if (p == 0)
{
create_standard_basis(j, evaluator_m);
create_zero_basis(evaluator_p);
}
else
{
create_zero_basis(evaluator_m);
create_standard_basis(j, evaluator_p);
}
if (buffer_dof_values)
for (unsigned int i = 0; i < evaluator_m.dofs_per_cell; ++i)
{
dof_buffer_m[i] = evaluator_m.begin_dof_values()[i];
dof_buffer_p[i] = evaluator_p.begin_dof_values()[i];
}
face_operation(evaluator_m, evaluator_p);
if (p == 0)
local_diagonal_m[j] = evaluator_m.begin_dof_values()[j];
else
local_diagonal_p[j] = evaluator_p.begin_dof_values()[j];
}
for (unsigned int j = 0; j < evaluator_m.dofs_per_cell; ++j)
{
evaluator_m.begin_dof_values()[j] = local_diagonal_m[j];
evaluator_p.begin_dof_values()[j] = local_diagonal_p[j];
}
}
else
{
evaluator_m.read_dof_values(src);
evaluator_p.read_dof_values(src);
if (buffer_dof_values)
for (unsigned int i = 0; i < evaluator_m.dofs_per_cell; ++i)
{
dof_buffer_m[i] = evaluator_m.begin_dof_values()[i];
dof_buffer_p[i] = evaluator_p.begin_dof_values()[i];
}
face_operation(evaluator_m, evaluator_p);
}
evaluator_m.distribute_local_to_global(dst);
evaluator_p.distribute_local_to_global(dst);
}
}
template <int dim>
void PoissonOperator<dim>::local_apply_boundary_face(
PoissonOperator::VectorType &,
const PoissonOperator::VectorType &,
const std::pair<unsigned int, unsigned int> &) const
{}
template <int dim>
template <typename EvaluatorType>
void PoissonOperator<dim>::do_poisson_cell_term(EvaluatorType &evaluator,
const unsigned int q) const
{
evaluator.submit_gradient(evaluator.get_gradient(q), q);
}
template <int dim>
template <typename EvaluatorType, typename Number2>
void PoissonOperator<dim>::do_flux_term(EvaluatorType &evaluator_m,
EvaluatorType &evaluator_p,
const Number2 &tau,
const unsigned int q) const
{
const auto normal_gradient_m = evaluator_m.get_normal_derivative(q);
const auto normal_gradient_p = evaluator_p.get_normal_derivative(q);
const auto value_m = evaluator_m.get_value(q);
const auto value_p = evaluator_p.get_value(q);
const auto jump_value = value_m - value_p;
const auto central_flux_gradient =
0.5 * (normal_gradient_m + normal_gradient_p);
const auto value_terms = central_flux_gradient - tau * jump_value;
evaluator_m.submit_value(-value_terms, q);
evaluator_p.submit_value(value_terms, q);
const auto gradient_terms = -0.5 * jump_value;
evaluator_m.submit_normal_derivative(gradient_terms, q);
evaluator_p.submit_normal_derivative(gradient_terms, q);
}
template <int dim>
template <typename EvaluatorType, typename Number2>
void PoissonOperator<dim>::do_boundary_flux_term_homogeneous(
EvaluatorType &evaluator_m,
const Number2 &tau,
const unsigned int q) const
{
const auto value = evaluator_m.get_value(q);
const auto normal_gradient = evaluator_m.get_normal_derivative(q);
const auto value_term = 2. * tau *
value - normal_gradient;
const auto normal_gradient_term = -
value;
evaluator_m.submit_value(value_term, q);
evaluator_m.submit_normal_derivative(normal_gradient_term, q);
}
template <int dim>
template <bool do_normal_hessians, bool do_values, typename EvaluatorType>
void PoissonOperator<dim>::do_gp_face_term(
EvaluatorType &evaluator_m,
EvaluatorType &evaluator_p,
const typename EvaluatorType::NumberType &masked_factor_value,
const typename EvaluatorType::NumberType &masked_factor_gradient,
const typename EvaluatorType::NumberType &masked_factor_hessian,
const unsigned int q) const
{
if (do_values)
{
const auto value_m = evaluator_m.get_value(q);
const auto value_p = evaluator_p.get_value(q);
const auto jump_value = value_m - value_p;
const auto value_term = masked_factor_value * jump_value;
evaluator_m.submit_value(value_term, q);
evaluator_p.submit_value(-value_term, q);
}
{
const auto normal_gradient_m = evaluator_m.get_normal_derivative(q);
const auto normal_gradient_p = evaluator_p.get_normal_derivative(q);
const auto jump_normal_gradient = normal_gradient_m - normal_gradient_p;
const auto gradient_term = masked_factor_gradient * jump_normal_gradient;
evaluator_m.submit_normal_derivative(gradient_term, q);
evaluator_p.submit_normal_derivative(-gradient_term, q);
}
if (do_normal_hessians)
{
const auto normal_hessian_m = evaluator_m.get_normal_hessian(q);
const auto normal_hessian_p = evaluator_p.get_normal_hessian(q);
const auto jump_normal_hessian = normal_hessian_m - normal_hessian_p;
const auto hessian_term = masked_factor_hessian * jump_normal_hessian;
evaluator_m.submit_normal_hessian(hessian_term, q);
evaluator_p.submit_normal_hessian(-hessian_term, q);
}
}
template <int dim>
template <typename EvaluatorType>
void PoissonOperator<dim>::do_rhs_cell_term(EvaluatorType &evaluator,
const unsigned int q) const
{
evaluator.submit_value(evaluate_function(rhs_function,
evaluator.quadrature_point(q)),
q);
}
template <int dim>
void PoissonOperator<dim>::do_local_apply_sipg_term(
PoissonOperator::FaceEvaluator &evaluator_m,
PoissonOperator::FaceEvaluator &evaluator_p,
const PoissonOperator::VectorizedArrayType *dof_ptr_m,
const PoissonOperator::VectorizedArrayType *dof_ptr_p,
const PoissonOperator::VectorizedArrayType &diameter) const
{
const auto tau = compute_interior_penalty_parameter(diameter);
evaluator_m.evaluate(dof_ptr_m,
evaluator_p.evaluate(dof_ptr_p,
for (const auto q : evaluator_m.quadrature_point_indices())
do_flux_term(evaluator_m, evaluator_p, tau, q);
}
template <int dim>
template <bool is_dg_>
void PoissonOperator<dim>::do_local_apply_gp_face_term(
PoissonOperator::FaceEvaluator &evaluator_m,
PoissonOperator::FaceEvaluator &evaluator_p,
const PoissonOperator::VectorizedArrayType *dof_ptr_m,
const PoissonOperator::VectorizedArrayType *dof_ptr_p,
const PoissonOperator::VectorizedArrayType &diameter,
const bool sum_into_values) const
{
if (is_dg_)
const unsigned int degree =
matrix_free->get_dof_handler(dof_index).get_fe().degree;
const bool do_hessians = degree > 1;
if (do_hessians)
degree <= 2,
ExcMessage(
"Face-based ghost-penalty stabilization only implemented up to degree 2!"));
evaluator_m.evaluate(dof_ptr_m, evaluation_flags);
evaluator_p.evaluate(dof_ptr_p, evaluation_flags);
const VectorizedArrayType factor_values = 0.5 /
diameter;
const VectorizedArrayType factor_gradient = 0.5 *
diameter;
const VectorizedArrayType factor_hessians =
if (do_hessians)
for (const auto q : evaluator_m.quadrature_point_indices())
do_gp_face_term<true, is_dg_>(evaluator_m,
evaluator_p,
factor_values,
factor_gradient,
factor_hessians,
q);
else
for (const auto q : evaluator_m.quadrature_point_indices())
do_gp_face_term<false, is_dg_>(evaluator_m,
evaluator_p,
factor_values,
factor_gradient,
factor_hessians,
q);
evaluator_m.integrate(evaluation_flags, sum_into_values);
evaluator_p.integrate(evaluation_flags, sum_into_values);
}
template <int dim>
bool PoissonOperator<dim>::is_inside_face(
std::pair<unsigned int, unsigned int> face_category) const
{
return is_inside(face_category.first) && is_inside(face_category.second);
}
template <int dim>
bool PoissonOperator<dim>::is_mixed_face(
std::pair<unsigned int, unsigned int> face_category) const
{
return (is_inside(face_category.first) &&
is_intersected(face_category.second)) ||
(is_intersected(face_category.first) &&
is_inside(face_category.second));
}
template <int dim>
bool PoissonOperator<dim>::is_intersected_face(
std::pair<unsigned int, unsigned int> face_category) const
{
return is_intersected(face_category.first) &&
is_intersected(face_category.second);
}
template <int dim>
typename PoissonOperator<dim>::VectorizedArrayType
PoissonOperator<dim>::compute_diameter_of_inner_face_batch(
unsigned int face_batch_index) const
{
const auto &face_info = matrix_free->get_face_info(face_batch_index);
for (unsigned int v = 0;
v < matrix_free->n_active_entries_per_face_batch(face_batch_index);
++v)
{
const auto cell_batch_index_interior =
face_info.cells_interior[v] / n_lanes;
const auto cell_lane_index_interior =
face_info.cells_interior[v] % n_lanes;
const auto cell_batch_index_exterior =
face_info.cells_exterior[v] / n_lanes;
const auto cell_lane_index_exterior =
face_info.cells_exterior[v] % n_lanes;
cell_diameter[cell_batch_index_interior][cell_lane_index_interior],
cell_diameter[cell_batch_index_exterior][cell_lane_index_exterior]);
}
}
template <int dim>
typename PoissonOperator<dim>::VectorizedArrayType
PoissonOperator<dim>::compute_interior_penalty_parameter(
const PoissonOperator::VectorizedArrayType &diameter) const
{
const auto k = matrix_free->get_dof_handler(dof_index).get_fe().degree;
return 2. * Utilities::fixed_power<2>(k + 1) /
diameter;
}
template <int dim>
class JacobiPreconditioner
{
public:
JacobiPreconditioner(const PoissonOperator<dim> &poisson_operator);
void vmult(VectorType &dst, const VectorType &src) const;
private:
};
template <int dim>
JacobiPreconditioner<dim>::JacobiPreconditioner(
const PoissonOperator<dim> &poisson_operator)
{
poisson_operator.initialize_dof_vector(inverse_diagonal);
poisson_operator.compute_diagonal(inverse_diagonal);
for (unsigned int i = 0; i < inverse_diagonal.locally_owned_size(); ++i)
{
if (
std::abs(inverse_diagonal.local_element(i)) > 1.0e-10)
inverse_diagonal.local_element(i) =
1.0 / inverse_diagonal.local_element(i);
else
inverse_diagonal.local_element(i) = 1.0;
}
}
template <int dim>
void JacobiPreconditioner<dim>::vmult(
JacobiPreconditioner::VectorType &dst,
const JacobiPreconditioner::VectorType &src) const
{
{
dst.scale(inverse_diagonal);
}
else
{
for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
dst.local_element(i) =
inverse_diagonal.local_element(i) * src.local_element(i);
}
}
template <int dim>
class PoissonSolver
{
using CellEvaluator =
using GenericCellEvaluator =
static constexpr unsigned int n_lanes = VectorizedArrayType::size();
public:
PoissonSolver();
void run(
bool is_dg,
unsigned int fe_degree);
private:
void make_grid();
void setup_discrete_level_set();
void distribute_dofs();
void setup_matrix_free();
void setup_mapping_data();
void solve();
void output_results() const;
double compute_L2_error() const;
unsigned int fe_degree;
PoissonOperator<dim> poisson_operator;
std::unique_ptr<FE_Q<dim>> fe_level_set;
std::unique_ptr<NonMatching::MappingInfo<dim, dim, VectorizedArrayType>>
mapping_info_cell;
std::unique_ptr<NonMatching::MappingInfo<dim, dim, VectorizedArrayType>>
mapping_info_surface;
std::unique_ptr<NonMatching::MappingInfo<dim, dim, VectorizedArrayType>>
mapping_info_faces;
const unsigned int dof_index = 0;
const unsigned int quad_index = 0;
bool is_dg;
};
template <int dim>
PoissonSolver<dim>::PoissonSolver()
: fe_degree(1)
, rhs_function(4.0)
, triangulation(MPI_COMM_WORLD)
triangulation.get_mpi_communicator()) == 0)
, level_set_dof_handler(triangulation)
, dof_handler(triangulation)
, mesh_classifier(level_set_dof_handler, level_set)
, mapping(1)
, is_dg(false)
{}
template <int dim>
void PoissonSolver<dim>::make_grid()
{
pcout << "Creating background mesh" << std::endl;
triangulation.clear();
triangulation.refine_global(2);
}
template <int dim>
void PoissonSolver<dim>::setup_discrete_level_set()
{
pcout << "Setting up discrete level set function" << std::endl;
fe_level_set = std::make_unique<FE_Q<dim>>(fe_degree);
level_set_dof_handler.distribute_dofs(*fe_level_set);
level_set.reinit(level_set_dof_handler.locally_owned_dofs(),
level_set_dof_handler),
triangulation.get_mpi_communicator());
signed_distance_sphere,
level_set);
level_set.update_ghost_values();
}
template <int dim>
void PoissonSolver<dim>::distribute_dofs()
{
pcout << "Distributing degrees of freedom" << std::endl;
std::unique_ptr<FE_Poly<dim>> fe;
if (is_dg)
fe = std::make_unique<FE_DGQ<dim>>(fe_degree);
else
fe = std::make_unique<FE_Q<dim>>(fe_degree);
fe_collection.push_back(*fe);
for (const auto &cell : dof_handler.active_cell_iterators() |
{
mesh_classifier.location_to_level_set(cell);
cell->set_active_fe_index(
cell->set_active_fe_index(static_cast<unsigned int>(
cell->set_active_fe_index(static_cast<unsigned int>(
else
cell->set_active_fe_index(static_cast<unsigned int>(
}
dof_handler.distribute_dofs(fe_collection);
}
template <int dim>
void PoissonSolver<dim>::setup_matrix_free()
{
pcout << "Setting up matrix-free" << std::endl;
affine_constraints.
close();
additional_data;
if (dof_handler.get_fe().degree > 1)
matrix_free.reinit(
mapping, dof_handler, affine_constraints, quadrature, additional_data);
}
template <int dim>
void PoissonSolver<dim>::setup_mapping_data()
{
pcout << "Setting up non matching mapping info" << std::endl;
const auto is_intersected_cell =
return mesh_classifier.location_to_level_set(cell) ==
};
std::vector<Quadrature<dim>> quad_vec_cells;
quad_vec_cells.reserve(
(matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches()) *
n_lanes);
std::vector<NonMatching::ImmersedSurfaceQuadrature<dim>> quad_vec_surface;
quad_vec_surface.reserve(
(matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches()) *
n_lanes * n_lanes);
q_collection1D, level_set_dof_handler, level_set);
std::vector<typename DoFHandler<dim>::cell_iterator> vector_accessors;
vector_accessors.reserve(
(matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches()) *
n_lanes);
for (unsigned int cell_batch = 0;
cell_batch <
matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches();
++cell_batch)
for (unsigned int v = 0; v < n_lanes; ++v)
{
if (v < matrix_free.n_active_entries_per_cell_batch(cell_batch))
vector_accessors.push_back(
matrix_free.get_cell_iterator(cell_batch, v));
else
vector_accessors.push_back(
matrix_free.get_cell_iterator(cell_batch, 0));
const auto &cell = vector_accessors.back();
if (is_intersected_cell(cell))
{
quadrature_generator.generate(cell);
quad_vec_cells.push_back(
quadrature_generator.get_inside_quadrature());
quad_vec_surface.push_back(
quadrature_generator.get_surface_quadrature());
}
else
{
quad_vec_cells.emplace_back();
quad_vec_surface.emplace_back();
}
}
mapping_info_cell = std::make_unique<
mapping_info_cell->
reinit_cells(vector_accessors, quad_vec_cells);
mapping_info_surface = std::make_unique<
mapping,
mapping_info_surface->
reinit_surface(vector_accessors, quad_vec_surface);
if (is_dg)
{
face_quadrature_generator(q_collection1D,
level_set_dof_handler,
level_set);
quad_vec_faces.reserve((matrix_free.n_inner_face_batches() +
matrix_free.n_boundary_face_batches()) *
n_lanes);
std::vector<
std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>>
vector_face_accessors_m;
vector_face_accessors_m.reserve(
(matrix_free.n_inner_face_batches() +
matrix_free.n_boundary_face_batches()) *
n_lanes);
unsigned int face_batch = 0;
for (; face_batch < matrix_free.n_inner_face_batches(); ++face_batch)
{
for (unsigned int v = 0; v < n_lanes; ++v)
{
if (v < matrix_free.n_active_entries_per_face_batch(face_batch))
vector_face_accessors_m.push_back(
matrix_free.get_face_iterator(face_batch, v, true));
else
vector_face_accessors_m.push_back(
matrix_free.get_face_iterator(face_batch, 0, true));
const auto &cell_m = vector_face_accessors_m.back().first;
const unsigned int f = vector_face_accessors_m.back().second;
if (is_intersected_cell(cell_m))
{
face_quadrature_generator.generate(cell_m, f);
quad_vec_faces.push_back(
face_quadrature_generator.get_inside_quadrature());
}
else
quad_vec_faces.emplace_back();
}
}
for (; face_batch < (matrix_free.n_inner_face_batches() +
matrix_free.n_boundary_face_batches());
++face_batch)
{
for (unsigned int v = 0; v < n_lanes; ++v)
{
if (v < matrix_free.n_active_entries_per_face_batch(face_batch))
vector_face_accessors_m.push_back(
matrix_free.get_face_iterator(face_batch, v, true));
else
vector_face_accessors_m.push_back(
matrix_free.get_face_iterator(face_batch, 0, true));
const auto &cell_m = vector_face_accessors_m.back().first;
const unsigned int f = vector_face_accessors_m.back().second;
if (is_intersected_cell(cell_m))
{
face_quadrature_generator.generate(cell_m, f);
quad_vec_faces.push_back(
face_quadrature_generator.get_inside_quadrature());
}
else
quad_vec_faces.emplace_back();
}
}
mapping_info_faces = std::make_unique<
mapping,
quad_vec_faces);
}
}
template <int dim>
void PoissonSolver<dim>::solve()
{
pcout << "Solving system" << std::endl;
JacobiPreconditioner<dim> jacobi_preconditioner(poisson_operator);
const unsigned int max_iterations = solution.size();
solver.solve(poisson_operator, solution, rhs, jacobi_preconditioner);
}
template <int dim>
void PoissonSolver<dim>::output_results() const
{
pcout << "Writing vtu file" << std::endl;
return cell->is_active() && cell->is_locally_owned() &&
mesh_classifier.location_to_level_set(cell) !=
});
"step-95",
0,
triangulation.get_mpi_communicator());
}
template <int dim>
class AnalyticalSolution :
public Function<dim>
{
public:
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>
double PoissonSolver<dim>::compute_L2_error() const
{
pcout << "Computing L2 error" << std::endl;
const QGauss<1> quadrature_1D(fe_degree + 1);
AnalyticalSolution<dim> analytical_solution;
double error_L2_squared = 0;
auto l2_kernel = [](auto &evaluator,
const unsigned int q) {
const VectorizedArrayType
value =
evaluate_function(analytical_solution_function,
evaluator.quadrature_point(q));
const auto difference = evaluator.get_value(q) -
value;
evaluator.submit_value(difference * difference, q);
};
unsigned int dummy = 0;
matrix_free.template cell_loop<unsigned int, VectorType>(
unsigned int &,
const VectorType &src,
const std::pair<unsigned int, unsigned int> &cell_range) {
CellEvaluator evaluator(matrix_free, dof_index, quad_index);
GenericCellEvaluator evaluator_cell(*mapping_info_cell, fe_dgq);
const auto cell_range_category =
matrix_free.get_cell_range_category(cell_range);
if (is_inside(cell_range_category))
{
for (unsigned int cell_batch_index = cell_range.first;
cell_batch_index < cell_range.second;
++cell_batch_index)
{
evaluator.reinit(cell_batch_index);
evaluator.read_dof_values(src);
for (unsigned int q : evaluator.quadrature_point_indices())
l2_kernel(evaluator, analytical_solution, q);
for (unsigned int v = 0;
v < matrix_free.n_active_entries_per_cell_batch(
cell_batch_index);
++v)
error_L2_squared += evaluator.integrate_value()[v];
}
}
else if (is_intersected(cell_range_category))
{
const auto dofs_per_cell = evaluator.dofs_per_cell;
for (unsigned int cell_batch_index = cell_range.first;
cell_batch_index < cell_range.second;
++cell_batch_index)
{
evaluator.reinit(cell_batch_index);
evaluator.read_dof_values(src);
for (unsigned int v = 0;
v < matrix_free.n_active_entries_per_cell_batch(
cell_batch_index);
++v)
{
cell_batch_index * n_lanes + v;
evaluator_cell.evaluate(
&evaluator.begin_dof_values()[0][v], dofs_per_cell),
for (const auto q :
evaluator_cell.quadrature_point_indices())
l2_kernel(evaluator_cell, analytical_solution, q);
error_L2_squared += evaluator_cell.integrate_value();
}
}
}
},
dummy,
solution);
triangulation.get_mpi_communicator()));
}
template <int dim>
void PoissonSolver<dim>::run(const bool is_dg_in,
const unsigned int fe_degree_in)
{
is_dg = is_dg_in;
fe_degree = fe_degree_in;
if (is_dg)
pcout << "Run DG convergence study with degree " << fe_degree
<< std::endl;
else
pcout << "Run CG convergence study with degree " << fe_degree
<< std::endl;
const unsigned int n_refinements = 3;
make_grid();
for (unsigned int cycle = 0; cycle <= n_refinements; cycle++)
{
pcout << "Refinement cycle " << cycle << std::endl;
triangulation.refine_global(1);
setup_discrete_level_set();
pcout << "Classifying cells" << std::endl;
mesh_classifier.reclassify();
distribute_dofs();
setup_matrix_free();
setup_mapping_data();
poisson_operator.reinit(matrix_free,
*mapping_info_cell,
*mapping_info_surface,
is_dg ? *mapping_info_faces :
*mapping_info_cell,
is_dg);
matrix_free.initialize_dof_vector(solution);
matrix_free.initialize_dof_vector(rhs);
poisson_operator.rhs(rhs, rhs_function);
solve();
if (cycle == n_refinements)
output_results();
const double error_L2 = compute_L2_error();
const double cell_side_length =
triangulation.begin_active(triangulation.n_levels() - 1)
->minimum_vertex_distance();
convergence_table.
add_value(
"Mesh size", cell_side_length);
convergence_table.
add_value(
"L2-Error", error_L2);
pcout << std::endl;
triangulation.get_mpi_communicator()) == 0)
pcout << std::endl;
}
pcout <<
"wall time: " << timer.
wall_time() <<
"\n" << std::endl;
}
}
int main(
int argc,
char **argv)
{
constexpr int dim = 2;
const unsigned int fe_degree = 2;
Step95::PoissonSolver<dim> poisson_solver;
poisson_solver.run(false , fe_degree);
poisson_solver.run(true , fe_degree);
}
void evaluate_convergence_rates(const std::string &data_column_key, const std::string &reference_column_key, const RateMode rate_mode, const unsigned int dim=2)
void set_flags(const FlagType &flags)
std::string write_vtu_with_pvtu_record(const std::string &directory, const std::string &filename_without_extension, const unsigned int counter, const MPI_Comm mpi_communicator, const unsigned int n_digits_for_counter=numbers::invalid_unsigned_int, const unsigned int n_groups=0) const
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={})
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 write_text(std::ostream &out, const TextOutputFormat format=table_with_headers) const
void add_value(const std::string &key, const T value)
void set_scientific(const std::string &key, const bool scientific)
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
* * * RotationFunction< dim, Number >::RotationFunction Number(dim)
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)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
UpdateFlags mapping_update_flags_inner_faces
UpdateFlags mapping_update_flags_boundary_faces
1.9.8