// @sect3{Initial conditions for vibrating membrane}
//
- // Function that provides the initial condition for the vibrating membrane
- // test case.
+ // In the following, let us first define a function that provides
+ // the initial condition for the vibrating membrane test case. It
+ // implementes both the initial pressure (component 0) and velocity
+ // (components 1 to 1 + dim).
+ //
+ // There is also a function that computes the duration of one
+ // oscillation.
template <int dim>
class InitialConditionVibratingMembrane : public Function<dim>
{
public:
InitialConditionVibratingMembrane(const double modes);
- // Function that the gives the initial pressure (comp 0) and velocity (comp
- // 1 to 1 + dim).
double value(const Point<dim> &p, const unsigned int comp) const final;
- // Function that calculates the duration of one oscillation.
double get_period_duration(const double speed_of_sound) const;
private:
// @sect3{Gauss pulse}
//
- // Function that provides the values of a pressure Gauss pulse.
+ // Next up is a function that provides the values of a pressure
+ // Gauss pulse. As with the previous function, it implements both
+ // the initial pressure (component 0) and velocity (components 1 to
+ // 1 + dim).
template <int dim>
class GaussPulse : public Function<dim>
{
public:
GaussPulse(const double shift_x, const double shift_y);
- // Function that the gives the initial pressure (comp 0) and velocity (comp
- // 1 to 1 + dim).
double value(const Point<dim> &p, const unsigned int comp) const final;
private:
static_assert(dim == 2, "Only implemented for dim==2");
}
- // Function that the gives the initial pressure (comp 0) and velocity (comp 1
- // to 1 + dim).
template <int dim>
double GaussPulse<dim>::value(const Point<dim> &p,
const unsigned int comp) const
{
// Helper function to check if a boundary ID is related to a non-matching
// face. A @c std::set that contains all non-matching boundary IDs is
- // handed over additionally to the face ID under question. This function
- // could certainly also be defined inline but this way the code is more easy
- // to read.
+ // handed over in addition to the face ID under question.
bool is_non_matching_face(
const std::set<types::boundary_id> &non_matching_face_ids,
const types::boundary_id face_id)
}
} // namespace HelperFunctions
- //@sect3{Material access}
+ // @sect3{Material parameter description}
+ //
+ // The following class stores the information if the fluid is
+ // homogeneous as well as the material properties at every cell.
+ // This class helps access the correct values without accessing a
+ // large vector of materials in the homogeneous case.
//
- // This class stores the information if the fluid is homogeneous
- // as well as the material properties at every cell.
- // This class helps to access the correct values without accessing
- // a large vector of materials in the homogeneous case.
+ // The class is provided a map from material ids to material
+ // properties (given as a pair of values for the speed of sound and
+ // the density). If the map has only one entry, the material is
+ // homogeneous -- using the same values everywhere -- and we can
+ // remember that fact in the `homogeneous` member variable and use it
+ // to optimize some code paths below. If the material is not
+ // homogeneous, we will fill a vector with the correct materials for
+ // each batch of cells; this information can then be access via
+ // FEEvaluationData::read_cell_data().
+ //
+ // As is usual when working with the MatrixFree framework, we will
+ // not only access material parameters at a single site, but for
+ // whole "batches" of cells. As a consequence, the functions below
+ // returned VectorizedArray objects for a batch at a time.
template <typename Number>
class CellwiseMaterialData
{
const MatrixFree<dim, Number, VectorizedArray<Number>> &matrix_free,
const std::map<types::material_id, std::pair<double, double>>
&material_id_map)
- // If the map is of size 1, the material is constant in every cell.
: homogeneous(material_id_map.size() == 1)
{
Assert(material_id_map.size() > 0,
- ExcMessage("No materials given to CellwiseMaterialData"));
+ ExcMessage("No materials given to CellwiseMaterialData."));
if (homogeneous)
{
- // In the homogeneous case we know the materials in the whole domain.
speed_of_sound_homogeneous = material_id_map.begin()->second.first;
density_homogeneous = material_id_map.begin()->second.second;
}
else
{
- // In the in-homogeneous case materials vary between cells. We are
- // filling a vector with the correct materials, that can be processed
- // via
- // @c read_cell_data().
const auto n_cell_batches =
matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches();
const AlignedVector<VectorizedArray<Number>> &get_speed_of_sound() const
{
- Assert(!homogeneous, ExcMessage("Use get_homogeneous_speed_of_sound()"));
+ Assert(!homogeneous, ExcMessage("Use get_homogeneous_speed_of_sound()."));
return speed_of_sound;
}
const AlignedVector<VectorizedArray<Number>> &get_density() const
{
- Assert(!homogeneous, ExcMessage("Use get_homogeneous_density()"));
+ Assert(!homogeneous, ExcMessage("Use get_homogeneous_density()."));
return density;
}
VectorizedArray<Number> get_homogeneous_speed_of_sound() const
{
- Assert(homogeneous, ExcMessage("Use get_speed_of_sound()"));
+ Assert(homogeneous, ExcMessage("Use get_speed_of_sound()."));
return speed_of_sound_homogeneous;
}
VectorizedArray<Number> get_homogeneous_density() const
{
- Assert(homogeneous, ExcMessage("Use get_density()"));
+ Assert(homogeneous, ExcMessage("Use get_density()."));
return density_homogeneous;
}
private:
const bool homogeneous;
- // Materials in the in-homogeneous case.
+ /* Materials in the inhomogeneous case. */
AlignedVector<VectorizedArray<Number>> speed_of_sound;
AlignedVector<VectorizedArray<Number>> density;
- // Materials in the homogeneous case.
+ /* Materials in the homogeneous case. */
VectorizedArray<Number> speed_of_sound_homogeneous;
VectorizedArray<Number> density_homogeneous;
};
- // To be able to access the material data in every cell in a thread safe way
- // @c MaterialEvaluation is used. Similar to @c FEEvaluation, every thread
- // creates its own instance and thus, there are no race conditions. For
- // in-homogeneous materials, a @c reinit_cell() or @c reinit_face() function
- // is used to set the correct material at the current cell batch. In the
- // homogeneous case the @c _reinit() functions don't have to reset the
- // materials.
+ // To be able to access the material data in every cell in a thread
+ // safe way, the following class @c MaterialEvaluation is
+ // used. Similar to @c FEEvaluation, functions below will create
+ // their own instances of this class; thus, there can be no race
+ // conditions even if these functions run multiple times in
+ // parallel. For inhomogeneous materials, a @c reinit_cell() or @c
+ // reinit_face() function is used to set the correct material at the
+ // current cell batch. In the homogeneous case the @c _reinit()
+ // functions don't have to reset the materials.
template <int dim, typename Number>
class MaterialEvaluation
{
{
if (material_data.is_homogeneous())
{
- // Set the material that is used in every cell.
speed_of_sound = material_data.get_homogeneous_speed_of_sound();
density = material_data.get_homogeneous_density();
}
return material_data.is_homogeneous();
}
- // Update the cell data, given a cell batch index.
+ // The following two functions update the data for the current
+ // cell or face, given a cell batch index. If the material is
+ // homogeneous, there is nothing to do. Otherwise, we reinit the
+ // FEEvaluation object and store the data for the current object.
void reinit_cell(const unsigned int cell)
{
- // In the homogeneous case we do not have to reset the cell data.
if (!material_data.is_homogeneous())
{
- // Reinit the FEEvaluation object and set the cell data.
phi.reinit(cell);
speed_of_sound =
phi.read_cell_data(material_data.get_speed_of_sound());
}
}
- // Update the cell data, given a face batch index.
void reinit_face(const unsigned int face)
{
- // In the homogeneous case we do not have to reset the cell data.
if (!material_data.is_homogeneous())
{
- // Reinit the FEFaceEvaluation object and set the cell data.
phi_face.reinit(face);
speed_of_sound =
phi_face.read_cell_data(material_data.get_speed_of_sound());
}
}
- // Return the speed of sound at the current cell batch.
+ // The following functions then return the speed of sound and
+ // density for the current cell batch.
VectorizedArray<Number> get_speed_of_sound() const
{
return speed_of_sound;
}
- // Return the density at the current cell batch.
VectorizedArray<Number> get_density() const
{
return density;
}
private:
- // Members needed for the in-homogeneous case.
+ /* Members needed for the inhomogeneous case. */
FEEvaluation<dim, -1, 0, 1, Number> phi;
FEFaceEvaluation<dim, -1, 0, 1, Number> phi_face;
- // Material defined at every cell.
+ /* Material defined at every cell. */
const CellwiseMaterialData<Number> &material_data;
- // Materials at current cell.
+ /* Materials at current cell. */
VectorizedArray<Number> speed_of_sound;
VectorizedArray<Number> density;
};
- //@sect3{Boundary conditions}
+ // @sect3{Boundary conditions}
//
// To be able to use the same kernel, for all face integrals we define
// a class that returns the needed values at boundaries. In this tutorial
const FEFaceEvaluation<dim, -1, 0, dim, Number> &velocity_m;
};
- //@sect3{Acoustic operator}
+ // @sect3{Acoustic operator}
+ //
+ // The following class then defines the acoustic operator. The class is
+ // heavily based on matrix-free methods. For a better understanding in
+ // matrix-free methods please refer to step-67.
//
- // Class that defines the acoustic operator. The class is heavily based on
- // matrix-free methods. For a better understanding in matrix-free methods
- // please refer to step-67.
+ // At the top we define a flag that decides whether we want to use
+ // mortaring. If the remote evaluators are set up with a
+ // VectorizedArray we are using point-to-point interpolation;
+ // otherwise we make use of Nitsche-type mortaring. The decision is
+ // made using `std::is_floating_point_v`, which is a variable that
+ // is true if the template argument is a floating point type (such
+ // as `float` or `double`) and false otherwise (in particular if the
+ // template argument is a vectorized array type).
template <int dim, typename Number, typename remote_value_type>
class AcousticOperator
{
- // If the remote evaluators are set up with a VectorizedArray we are
- // using point-to-point interpolation. Otherwise we make use of
- // Nitsche-type mortaring.
static constexpr bool use_mortaring =
std::is_floating_point_v<remote_value_type>;
has to be provided."));
}
- // Function to evaluate the acoustic operator.
+ // The following function then evaluates the acoustic operator.
+ // It first updates the precomputed values in corresponding the
+ // FERemoteEvaluation objects. The material parameters do not change and
+ // thus, we do not have to update precomputed values in @c c_r_eval and @c
+ // rho_r_eval.
+ //
+ // It then either performs a matrix-free loop with Nitsche-type
+ // mortaring at non-matching faces (if `use_mortaring` is true) or
+ // with point-to-point interpolation at non-matching faces (in the
+ // `else` branch). The difference is only in the third argument to
+ // the loop function, denoting the function object that is
+ // executed at boundary faces.
template <typename VectorType>
void evaluate(VectorType &dst, const VectorType &src) const
{
- // Update the precomputed values in corresponding the FERemoteEvaluation
- // objects. The material parameters do not change and thus, we do
- // not have to update precomputed values in @c c_r_eval and @c rho_r_eval.
pressure_r_eval->gather_evaluate(src, EvaluationFlags::values);
velocity_r_eval->gather_evaluate(src, EvaluationFlags::values);
if constexpr (use_mortaring)
{
- // Perform matrix free loop with Nitsche-type mortaring at
- // non-matching faces.
matrix_free.loop(
&AcousticOperator::local_apply_cell<VectorType>,
&AcousticOperator::local_apply_face<VectorType>,
}
else
{
- // Perform matrix free loop with point-to-point interpolation at
- // non-matching faces.
matrix_free.loop(
&AcousticOperator::local_apply_cell<VectorType>,
&AcousticOperator::local_apply_face<VectorType>,
}
}
+ // In the `private` section of the class, we then define the
+ // functions that evaluate volume, interior face, and boundary
+ // face integrals. The concrete terms these functions evaluate are
+ // stated in the documentation at the top of this tutorial
+ // program. Each of these functions has its own object of type
+ // `MaterialEvaluation` that provides access to the material at
+ // each cell or face.
private:
- // This function evaluates the volume integrals.
template <typename VectorType>
void local_apply_cell(
const MatrixFree<dim, Number> &matrix_free,
FEEvaluation<dim, -1, 0, 1, Number> pressure(matrix_free, 0, 0, 0);
FEEvaluation<dim, -1, 0, dim, Number> velocity(matrix_free, 0, 0, 1);
- // Class that gives access to the material at each cell
MaterialEvaluation material(matrix_free, *material_data);
for (unsigned int cell = cell_range.first; cell < cell_range.second;
pressure.gather_evaluate(src, EvaluationFlags::gradients);
velocity.gather_evaluate(src, EvaluationFlags::gradients);
- // Get the materials at the corresponding cell. Since we
+ // Get the materials on the corresponding cell. Since we
// introduced @c MaterialEvaluation we can write the code
- // independent if the material is homogeneous or in-homogeneous.
+ // independent of whether the material is homogeneous or
+ // inhomogeneous.
material.reinit_cell(cell);
const auto c = material.get_speed_of_sound();
const auto rho = material.get_density();
- for (unsigned int q : pressure.quadrature_point_indices())
+
+ for (const unsigned int q : pressure.quadrature_point_indices())
{
pressure.submit_value(rho * c * c * velocity.get_divergence(q),
q);
}
}
- // This function evaluates the fluxes at faces between cells with the same
- // material. If boundary faces are under consideration fluxes into
+ // This next function evaluates the fluxes at faces between cells with the
+ // same material. If boundary faces are under consideration fluxes into
// neighboring faces do not have to be considered which is enforced via
// `weight_neighbor=false`. For non-matching faces the fluxes into
// neighboring faces are not considered as well. This is because we iterate
// over each side of the non-matching face separately (similar to a cell
// centric loop).
+ //
+ // In this and following functions, we also introduce the factors
+ // $\tau$ and $\gamma$ that are computed from $\rho$ and $c$ along
+ // interfaces and that appear in the bilinear forms. Their
+ // definitions are provided in the introduction.
template <bool weight_neighbor,
typename InternalFaceIntegratorPressure,
typename InternalFaceIntegratorVelocity,
const typename InternalFaceIntegratorPressure::value_type c,
const typename InternalFaceIntegratorPressure::value_type rho) const
{
- // Compute penalty parameters from material parameters.
const auto tau = 0.5 * rho * c;
const auto gamma = 0.5 / (rho * c);
- for (unsigned int q : pressure_m.quadrature_point_indices())
+ for (const unsigned int q : pressure_m.quadrature_point_indices())
{
const auto n = pressure_m.normal_vector(q);
const auto pm = pressure_m.get_value(q);
const MaterialIntegrator &c_r,
const MaterialIntegrator &rho_r) const
{
- // Interior material information is constant over quadrature points
const auto tau_m = 0.5 * rho * c;
const auto gamma_m = 0.5 / (rho * c);
- for (unsigned int q : pressure_m.quadrature_point_indices())
+ for (const unsigned int q : pressure_m.quadrature_point_indices())
{
- // The material at the neighboring face might vary in every quadrature
- // point.
const auto c_p = c_r.get_value(q);
const auto rho_p = rho_r.get_value(q);
const auto tau_p = 0.5 * rho_p * c_p;
FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_p(
matrix_free, false, 0, 0, 1);
- // Class that gives access to the material at each cell
MaterialEvaluation material(matrix_free, *material_data);
for (unsigned int face = face_range.first; face < face_range.second;
- face++)
+ ++face)
{
velocity_m.reinit(face);
velocity_p.reinit(face);
}
- //@sect4{Matrix-free boundary function for point-to-point interpolation}
+ // @sect4{Matrix-free boundary function for point-to-point interpolation}
//
- // This function evaluates the boundary face integrals and the
+ // The following function then evaluates the boundary face integrals and the
// non-matching face integrals using point-to-point interpolation.
template <typename VectorType>
void local_apply_boundary_face_point_to_point(
const VectorType &src,
const std::pair<unsigned int, unsigned int> &face_range) const
{
- // Standard face evaluators.
FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
matrix_free, true, 0, 0, 0);
FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
matrix_free, true, 0, 0, 1);
- // Classes that return the correct BC values.
BCEvaluationP pressure_bc(pressure_m);
BCEvaluationU velocity_bc(velocity_m);
- // Class that gives access to the material at each cell
MaterialEvaluation material(matrix_free, *material_data);
// Remote evaluators.
auto rho_r = rho_r_eval->get_data_accessor();
for (unsigned int face = face_range.first; face < face_range.second;
- face++)
+ ++face)
{
velocity_m.reinit(face);
pressure_m.reinit(face);
material.reinit_face(face);
+ // If we are considering a homogeneous material, do not use the
+ // inhomogeneous fluxes. While it would be possible
+ // to use the inhomogeneous fluxes they are more expensive to
+ // compute.
if (material.is_homogeneous())
{
- // If homogeneous material is considered do not use the
- // inhomogeneous fluxes. While it would be possible
- // to use the inhomogeneous fluxes they are more expensive to
- // compute.
evaluate_face_kernel<false>(pressure_m,
velocity_m,
pressure_r,
}
else
{
- // If inhomogeneous material is considered use the
- // in-homogeneous fluxes.
c_r.reinit(face);
rho_r.reinit(face);
evaluate_face_kernel_inhomogeneous(
}
}
- //@sect4{Matrix-free boundary function for Nitsche-type mortaring}
+ // @sect4{Matrix-free boundary function for Nitsche-type mortaring}
//
// This function evaluates the boundary face integrals and the
// non-matching face integrals using Nitsche-type mortaring.
const VectorType &src,
const std::pair<unsigned int, unsigned int> &face_range) const
{
- // Standard face evaluators for BCs.
FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
matrix_free, true, 0, 0, 0);
FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
velocity_r_mortar.reinit(face * n_lanes + v);
pressure_r_mortar.reinit(face * n_lanes + v);
+ // As above, if we are considering a homogeneous
+ // material, do not use the inhomogeneous
+ // fluxes. Since we are operating on face @c v we
+ // call @c material.get_density()[v].
if (material.is_homogeneous())
{
- // If homogeneous material is considered do not use the
- // inhomogeneous fluxes. While it would be possible
- // to use the inhomogeneous fluxes they are more
- // expensive to compute. Since we are operating on face @c
- // v we call @c material.get_density()[v].
evaluate_face_kernel<false>(
pressure_m_mortar,
velocity_m_mortar,
}
else
{
- // Same as in @c local_apply_boundary_face_point_to_point().
velocity_m.reinit(face);
pressure_m.reinit(face);
const MatrixFree<dim, Number> &matrix_free;
- // CellwiseMaterialData is stored as shared pointer with the same
- // argumentation.
const std::shared_ptr<CellwiseMaterialData<Number>> material_data;
const std::set<types::boundary_id> remote_face_ids;
}
};
- //@sect3{Inverse mass operator}
+ // @sect3{Inverse mass operator}
//
- // Class to apply the inverse mass operator.
+ // For the time stepping methods below, we need the inverse mass
+ // operator. We apply it via a loop over all (batches of) cells as
+ // always, where the contribution of each cell is computed in a
+ // matrix-free way:
template <int dim, typename Number>
class InverseMassOperator
{
: matrix_free(matrix_free)
{}
- // Function to apply the inverse mass operator.
template <typename VectorType>
void apply(VectorType &dst, const VectorType &src) const
{
}
private:
- // Apply the inverse mass operator onto every cell batch.
template <typename VectorType>
void local_apply_cell(
const MatrixFree<dim, Number> &mf,
const MatrixFree<dim, Number> &matrix_free;
};
- //@sect3{Runge-Kutta time-stepping}
+ // @sect3{Runge-Kutta time-stepping}
//
// This class implements a Runge-Kutta scheme of order 2.
template <int dim, typename Number, typename remote_value_type>
, acoustic_operator(acoustic_operator)
{}
- // Set up and run time loop.
+ // The `run()` function of this class contains the time loop. It
+ // starts by initializing some member variables (such as
+ // short-cuts to objects that describe the finite element, its
+ // properties, and the mapping; an by initializing vectors). It
+ // then computes the time step size via minimum element edge
+ // length. Assuming non-distorted elements, we can compute the
+ // edge length as the distance between two vertices. From this,
+ // we can then obtain the time step size via the CFL condition.
void run(const MatrixFree<dim, Number> &matrix_free,
const double cr,
const double end_time,
const Function<dim> &initial_condition,
const std::string &vtk_prefix)
{
- // Get needed members of matrix free.
const auto &dof_handler = matrix_free.get_dof_handler();
const auto &mapping = *matrix_free.get_mapping_info().mapping;
const auto degree = dof_handler.get_fe().degree;
- // Initialize needed Vectors.
VectorType solution;
matrix_free.initialize_dof_vector(solution);
VectorType solution_temp;
matrix_free.initialize_dof_vector(solution_temp);
- // Set the initial condition.
HelperFunctions::set_initial_condition(matrix_free,
initial_condition,
solution);
- // Compute time step size: Compute minimum element edge length.
- // We assume non-distorted elements, therefore we only compute
- // the distance between two vertices
double h_local_min = std::numeric_limits<double>::max();
for (const auto &cell : dof_handler.active_cell_iterators())
h_local_min =
const double h_min =
Utilities::MPI::min(h_local_min, dof_handler.get_communicator());
- // Compute constant time step size via the CFL condition.
const double dt =
cr * HelperFunctions::compute_dt_cfl(h_min, degree, speed_of_sound);
- // Perform time integration loop.
+ // We can then perform the time integration loop:
double time = 0.0;
unsigned int timestep = 0;
while (time < end_time)
{
- // Write output.
HelperFunctions::write_vtu(solution,
matrix_free.get_dof_handler(),
mapping,
"step_89-" + vtk_prefix +
std::to_string(timestep));
- // Perform a single time step.
std::swap(solution, solution_temp);
time += dt;
- timestep++;
+ ++timestep;
perform_time_step(dt, solution, solution_temp);
}
}
+ // The main work of this class is done by a `private` member
+ // function that performs one Runge-Kutta 2 time step. Recall that
+ // this method requires two sub-steps ("stages") computing
+ // intermediate values `k1` and `k2`, after which the intermediate
+ // values are summed with weights to obtain the new solution at
+ // the end of the time step. The RK2 method allows for the
+ // elimination of the intermediate vector `k2` by utilizing the
+ // `dst` vector as temporary storage.
private:
- // Perform one Runge-Kutta 2 time step.
void
perform_time_step(const double dt, VectorType &dst, const VectorType &src)
{
VectorType k1 = src;
- // First stage.
+ /* First stage. */
evaluate_stage(k1, src);
- // Second stage.
+ /* Second stage. */
k1.sadd(0.5 * dt, 1.0, src);
evaluate_stage(dst, k1);
+
+ /* Summing things into the output vector. */
dst.sadd(dt, 1.0, src);
}
- // Evaluate a single Runge-Kutta stage.
+ // Evaluating a single Runge-Kutta stage is a straightforward step
+ // that really only requires applying the operator once, and then
+ // applying the inverse of the mass matrix.
void evaluate_stage(VectorType &dst, const VectorType &src)
{
- // Evaluate the stage
acoustic_operator->evaluate(dst, src);
dst *= -1.0;
inverse_mass_operator->apply(dst, dst);
}
- // Needed operators.
const std::shared_ptr<InverseMassOperator<dim, Number>>
inverse_mass_operator;
const std::shared_ptr<AcousticOperator<dim, Number, remote_value_type>>
// @sect3{Construction of non-matching triangulations}
//
- // This function creates a two dimensional squared triangulation
- // that spans from (0,0) to (1,1). It consists of two sub-domains.
- // The left sub-domain spans from (0,0) to (0.525,1). The right
- // sub-domain spans from (0.525,0) to (1,1). The left sub-domain has
- // three times smaller elements compared to the right sub-domain.
+ // Let us now make our way to the higher-level functions of this program.
+ //
+ // The first of these functions creates a two dimensional square
+ // triangulation that spans from $(0,0)$ to $(1,1)$. It consists of
+ // two sub-domains. The left sub-domain spans from $(0,0)$ to
+ // $(0.525,1)$. The right sub-domain spans from $(0.525,0)$ to
+ // $(1,1)$. The left sub-domain has elements that are three times
+ // smaller compared to the ones for the right sub-domain.
+ //
+ // At non-matching interfaces, we need to provide different boundary
+ // IDs for the cells that make up the two parts (because, while they
+ // may be physically adjacent, they are not logically neighbors
+ // given that the faces of cells on both sides do not match, and the
+ // Triangulation class will therefore treat the interface between
+ // the two parts as a "boundary"). These boundary IDs have to differ
+ // because later on RemotePointEvaluation has to search for remote
+ // points for each face, that are defined in the same mesh (since we
+ // merge the mesh) but not on the same side of the non-matching
+ // interface. As a consequence, we declare at the top symbolic names
+ // for these boundary indicators, and ensure that we return a set
+ // with these values to the caller for later use.
+ //
+ // The actual mesh is then constructed by first constructing the
+ // left and right parts separately (setting material ids to zero and
+ // one, respectively), and using the appropriate boundary ids for
+ // all parts of the mesh. We then use
+ // GridGenerator::merge_triangulations() to combine them into one
+ // (non-matching) mesh. We have to pay attention that should the two
+ // sub-triangulations have vertices at the same locations, that they
+ // are not merged (connecting the two triangulations logically)
+ // since we want the interface to be an actual boundary. We achieve
+ // this by providing a tolerance of zero for the merge, see the
+ // documentation of the function
+ // GridGenerator::merge_triangulations().
template <int dim>
void build_non_matching_triangulation(
Triangulation<dim> &tria,
{
const double length = 1.0;
- // At non-matching interfaces, we provide different boundary
- // IDs. These boundary IDs have to differ because later on
- // RemotePointEvaluation has to search for remote points for
- // each face, that are defined in the same mesh (since we merge
- // the mesh) but not on the same side of the non-matching interface.
const types::boundary_id non_matching_id_left = 98;
const types::boundary_id non_matching_id_right = 99;
- // Provide this information to the caller.
- non_matching_faces.insert(non_matching_id_left);
- non_matching_faces.insert(non_matching_id_right);
+ non_matching_faces = {non_matching_id_left, non_matching_id_right};
- // Construct left part of mesh.
+ /* Construct left part of mesh. */
Triangulation<dim> tria_left;
const unsigned int subdiv_left = 11;
GridGenerator::subdivided_hyper_rectangle(tria_left,
{0.0, 0.0},
{0.525 * length, length});
- // The left part of the mesh has the material ID 0.
for (const auto &cell : tria_left.active_cell_iterators())
cell->set_material_id(0);
-
- // The right face is non-matching. All other boundary IDs
- // are set to 0.
for (const auto &face : tria_left.active_face_iterators())
if (face->at_boundary())
{
face->set_boundary_id(non_matching_id_left);
}
- // Construct right part of mesh.
+ /* Construct right part of mesh. */
Triangulation<dim> tria_right;
const unsigned int subdiv_right = 4;
GridGenerator::subdivided_hyper_rectangle(tria_right,
{0.525 * length, 0.0},
{length, length});
- // The right part of the mesh has the material ID 1.
for (const auto &cell : tria_right.active_cell_iterators())
cell->set_material_id(1);
-
- // The left face is non-matching. All other boundary IDs
- // are set to 0.
for (const auto &face : tria_right.active_face_iterators())
if (face->at_boundary())
{
face->set_boundary_id(non_matching_id_right);
}
- // Merge triangulations with tolerance 0 to ensure no vertices
- // are merged, see the documentation of the function
- // @c merge_triangulations().
+ /* Merge triangulations. */
GridGenerator::merge_triangulations(tria_left,
tria_right,
tria,
/*tolerance*/ 0.,
/*copy_manifold_ids*/ false,
/*copy_boundary_ids*/ true);
+
+ /* Refine the result. */
tria.refine_global(refinements);
}
- // @sect3{Set up and run point-to-point interpolation}
- //
- // The main purpose of this function is to fill a
- // `FERemoteEvaluationCommunicator` object that is needed for point-to-point
- // interpolation. Additionally, the corresponding remote evaluators are set up
- // using this remote communicator. Eventually, the operators are handed to the
- // time integrator that runs the simulation.
+ // @sect3{Set up and running of the two schemes}
+
+ // @sect4{Setup and running of the point-to-point interpolation scheme}
//
+ // We are now at the two functions that run the overall schemes (the
+ // point-to-point and the mortaring schemes). The first of these
+ // functions fills a FERemoteEvaluationCommunicator object that is
+ // needed for point-to-point interpolation. Additionally, the
+ // corresponding remote evaluators are set up using this remote
+ // communicator. Eventually, the operators are handed to the time
+ // integrator that runs the simulation.
template <int dim, typename Number>
void run_with_point_to_point_interpolation(
const MatrixFree<dim, Number> &matrix_free,
const auto &dof_handler = matrix_free.get_dof_handler();
const auto &tria = dof_handler.get_triangulation();
- // Communication objects know about the communication pattern. I.e.,
+ // Communication objects know about the communication pattern. That is,
// they know about the cells and quadrature points that have to be
// evaluated at remote faces. This information is given via
// RemotePointEvaluation. Additionally, the communication objects
//
// For the standard case of point to point-to-point interpolation without
// any heuristic we make use of the utility function
- // @c compute_remote_communicator_faces_point_to_point_interpolation().
- // Please refer to this function to see how to manually set up the
- // remote communicator from outside.
-
+ // Utilities::compute_remote_communicator_faces_point_to_point_interpolation().
+ // Please refer to the documentation of this function to see how to manually
+ // set up the remote communicator from outside.
+ //
+ // Among the inputs for the remote communicator we need a list of
+ // boundary ids for the non-matching faces, along with a function
+ // object for each boundary id that returns a vector of true/false
+ // flags in which exactly the vertices of cells of the
+ // triangulation are marked that have a face at the boundary id in
+ // question.
std::vector<
std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
non_matching_faces_marked_vertices;
-
for (const auto &nm_face : non_matching_faces)
{
- // Sufficient lambda, that rules out all cells connected to the current
- // side of the non-matching interface to avoid self intersections.
- auto marked_vertices = [&]() {
- // only search points at cells that are not connected to
- // @c nm_face
+ auto marked_vertices = [&nm_face, &tria]() -> std::vector<bool> {
std::vector<bool> mask(tria.n_vertices(), true);
for (const auto &cell : tria.active_cell_iterators())
// We are using point-to-point interpolation and can therefore
// easily access the corresponding data at face batches. This
- // is why we use a @c VectorizedArray as @c remote_value_type
+ // is why we use a @c VectorizedArray as @c remote_value_type:
using remote_value_type = VectorizedArray<Number>;
- // Set up FERemoteEvaluation object that accesses the pressure
- // at remote faces.
+ // We then set up FERemoteEvaluation objects that access the
+ // pressure and velocity at remote faces, along with an object to
+ // describe cell-wise material data.
const auto pressure_r =
std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
remote_communicator, dof_handler, /*first_selected_component*/ 0);
- // Set up FERemoteEvaluation object that accesses the velocity
- // at remote faces.
const auto velocity_r =
std::make_shared<FERemoteEvaluation<dim, dim, remote_value_type>>(
remote_communicator, dof_handler, /*first_selected_component*/ 1);
- // Set up cell-wise material data.
const auto material_data =
std::make_shared<CellwiseMaterialData<Number>>(matrix_free, materials);
matrix_free.get_dof_handler().get_triangulation(),
/*first_selected_component*/ 0);
+ // If the domain is not homogeneous, i.e., if material parameters
+ // change from cell to cell, we initialize and fill DoF vectors
+ // that contain the material properties. Materials do not change
+ // during the simulation, therefore there is no need to ever
+ // compute the values after the first @c gather_evaluate() (at the
+ // end of the following block) again.
if (!material_data->is_homogeneous())
{
- // Initialize and fill DoF vectors that contain the materials.
Vector<Number> c(
matrix_free.get_dof_handler().get_triangulation().n_active_cells());
Vector<Number> rho(
materials.at(cell->material_id()).second;
}
- // Materials do not change during the simulation, therefore
- // there is no need to precompute the values after
- // the first @c gather_evaluate() again.
c_r->gather_evaluate(c, EvaluationFlags::values);
rho_r->gather_evaluate(rho, EvaluationFlags::values);
}
- // Set up inverse mass operator.
+ // Next, we set up the inverse mass operator and the acoustic
+ // operator. Using `remote_value_type=VectorizedArray<Number>`
+ // makes the operator use point-to-point interpolation. These two
+ // objects are then used to create a `RungeKutta2` object to
+ // perform the time integration.
+ //
+ // We also compute the maximum speed of sound, needed for the
+ // computation of the time-step size, and then run the time integrator.
+ // For the examples considered here, we found a limiting Courant number of
+ // $\mathrm{Cr}\approx 0.36$ to maintain stability. To ensure, the
+ // error of the temporal discretization is small, we use a considerably
+ // smaller Courant number of $0.2$.
const auto inverse_mass_operator =
std::make_shared<InverseMassOperator<dim, Number>>(matrix_free);
- // Set up the acoustic operator. Using
- // `remote_value_type=VectorizedArray<Number>` makes the operator use
- // point-to-point interpolation.
const auto acoustic_operator =
std::make_shared<AcousticOperator<dim, Number, remote_value_type>>(
matrix_free,
c_r,
rho_r);
- // Compute the the maximum speed of sound, needed for the computation of
- // the time-step size.
+ RungeKutta2<dim, Number, remote_value_type> time_integrator(
+ inverse_mass_operator, acoustic_operator);
+
double speed_of_sound_max = 0.0;
for (const auto &mat : materials)
speed_of_sound_max = std::max(speed_of_sound_max, mat.second.first);
- // Set up time integrator.
- RungeKutta2<dim, Number, remote_value_type> time_integrator(
- inverse_mass_operator, acoustic_operator);
- // For considered examples, we found a limiting Courant number of
- // $\mathrm{Cr}\approx 0.36$ to maintain stability. To ensure, the
- // error of the temporal discretization is small, we use a considerably
- // smaller Courant number of $0.2$.
time_integrator.run(matrix_free,
/*Cr*/ 0.2,
end_time,
vtk_prefix);
}
- // @sect3{Set up and run Nitsche-type mortaring}
- //
- // The main purpose of this function is to fill a
- // `FERemoteEvaluationCommunicator` object that is needed for Nitsche-type
- // mortaring. Additionally, the corresponding remote evaluators are set up
- // using this remote communicator. Eventually, the operators are handed to the
- // time integrator that runs the simulation.
+ // @sect4{Setup and running of the Nitsche-type mortaring scheme}
//
+ // The alternative to the previous function is to use the mortaring
+ // scheme -- implemented in the following function. This function
+ // can only be run when deal.II is configured using CGAL (but the
+ // previous function can still be used without CGAL), and so errors
+ // out if CGAL is not available.
template <int dim, typename Number>
void run_with_nitsche_type_mortaring(
const MatrixFree<dim, Number> &matrix_free,
const auto &mapping = *matrix_free.get_mapping_info().mapping;
const auto n_quadrature_pnts = matrix_free.get_quadrature().size();
- // In case of Nitsche-type mortaring a vector of pairs with cell iterator
- // and face number is needed as communication object.
- // @c FERemoteCommunicationObjectFaces is a container to store this
- // information.
- //
- // For the standard case of Nitsche-type mortaring without
- // any heuristic we make use of the utility function
- // @c compute_remote_communicator_faces_nitsche_type_mortaring().
- // Please refer to this function to see how to manually set up the
- // remote communicator from outside and how to reinit
- // NonMatching::MappingInfo.
-
std::vector<
std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
non_matching_faces_marked_vertices;
for (const auto &nm_face : non_matching_faces)
{
- // Sufficient lambda, that rules out all cells connected to the current
- // side of the non-matching interface to avoid self intersections.
auto marked_vertices = [&]() {
- // only search points at cells that are not connected to
- // @c nm_face
std::vector<bool> mask(tria.n_vertices(), true);
for (const auto &cell : tria.active_cell_iterators())
std::make_pair(nm_face, marked_vertices));
}
- // Quadrature points are arbitrarily distributed on each non-matching
+ // The only parts in this function that are functionally different
+ // from the previous one follow here.
+ //
+ // First, quadrature points are arbitrarily distributed on each non-matching
// face. Therefore, we have to make use of FEFacePointEvaluation.
// FEFacePointEvaluation needs NonMatching::MappingInfo to work at the
// correct quadrature points that are in sync with used FERemoteEvaluation
0,
nm_mapping_info.get());
- // Same as above but since quadrature points are aribtrarily distributed
- // we have to consider each face in a batch separately and can not make
- // use of @c VecorizedArray.
+ // Second, since quadrature points are arbitrarily distributed we
+ // have to consider each face in a batch separately and can not
+ // make use of @c VecorizedArray.
using remote_value_type = Number;
+ // The rest of the code is then identical to what we had in the
+ // previous function (though it functions differently because of
+ // the difference in `remote_value_type`).
const auto pressure_r =
std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
remote_communicator, dof_handler, /*first_selected_component*/ 0);
rho_r->gather_evaluate(rho, EvaluationFlags::values);
}
- // Set up inverse mass operator.
const auto inverse_mass_operator =
std::make_shared<InverseMassOperator<dim, Number>>(matrix_free);
- // Set up the acoustic operator. Using `remote_value_type=Number`
- // makes the operator use Nitsche-type mortaring.
const auto acoustic_operator =
std::make_shared<AcousticOperator<dim, Number, remote_value_type>>(
matrix_free,
rho_r,
nm_mapping_info);
- // Compute the the maximum speed of sound, needed for the computation of
- // the time-step size.
+ RungeKutta2<dim, Number, remote_value_type> time_integrator(
+ inverse_mass_operator, acoustic_operator);
+
double speed_of_sound_max = 0.0;
for (const auto &mat : materials)
speed_of_sound_max = std::max(speed_of_sound_max, mat.second.first);
-
- // Set up time integrator.
- RungeKutta2<dim, Number, remote_value_type> time_integrator(
- inverse_mass_operator, acoustic_operator);
-
- // Run time loop with Courant number $0.2$.
time_integrator.run(matrix_free,
/*Cr*/ 0.2,
end_time,
//
// Finally, the `main()` function executes the different versions of handling
// non-matching interfaces.
+//
+// Similar to step-87, the minimum requirement of this tutorial is MPI.
+// The parallel::distributed::Triangulation class is used if deal.II is
+// configured with p4est. Otherwise parallel::shared::Triangulation is used.
int main(int argc, char *argv[])
{
using namespace dealii;
const unsigned int refinements = 1;
const unsigned int degree = 3;
- // Construct non-matching triangulation and fill non-matching boundary IDs.
-
- // Similar to step-87, the minimum requirement of this tutorial is MPI.
- // The parallel::distributed::Triangulation class is used if deal.II is
- // configured with p4est. Otherwise parallel::shared::Triangulation is used.
#ifdef DEAL_II_WITH_P4EST
parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
#else
pcout << " - Refinement level: " << refinements << std::endl;
pcout << " - Number of cells: " << tria.n_cells() << std::endl;
- // Set up MatrixFree.
-
pcout << "Create DoFHandler..." << std::endl;
DoFHandler<dim> dof_handler(tria);
dof_handler.distribute_dofs(FESystem<dim>(FE_DGQ<dim>(degree) ^ (dim + 1)));
MappingQ1<dim>(), dof_handler, constraints, QGauss<dim>(degree + 1), data);
- //@sect4{Run vibrating membrane test case}
+ // @sect4{Run vibrating membrane test case} Homogeneous pressure}
+ // Dirichlet boundary conditions are applied for
+ // simplicity. Therefore, modes can not be chosen arbitrarily.
pcout << "Run vibrating membrane test case..." << std::endl;
- // Vibrating membrane test case:
- //
- // Homogeneous pressure DBCs are applied for simplicity. Therefore,
- // modes can not be chosen arbitrarily.
const double modes = 10.0;
std::map<types::material_id, std::pair<double, double>> homogeneous_material;
homogeneous_material[numbers::invalid_material_id] = std::make_pair(1.0, 1.0);
const auto initial_solution_membrane =
Step89::InitialConditionVibratingMembrane<dim>(modes);
+ /* Run vibrating membrane test case using point-to-point interpolation: */
pcout << " - Point-to-point interpolation: " << std::endl;
- // Run vibrating membrane test case using point-to-point interpolation:
-
Step89::run_with_point_to_point_interpolation(
matrix_free,
non_matching_faces,
initial_solution_membrane,
"vm-p2p");
+ /* Run vibrating membrane test case using Nitsche-type mortaring: */
pcout << " - Nitsche-type mortaring: " << std::endl;
- // Run vibrating membrane test case using Nitsche-type mortaring:
Step89::run_with_nitsche_type_mortaring(
matrix_free,
non_matching_faces,
initial_solution_membrane,
"vm-nitsche");
- //@sect4{Run test case with in-homogeneous material}
- pcout << "Run test case with in-homogeneous material..." << std::endl;
- // In-homogeneous material test case:
- //
- // Run simple test case with in-homogeneous material and Nitsche-type
+ // @sect4{Run test case with inhomogeneous material}
+ // Run simple test case with inhomogeneous material and Nitsche-type
// mortaring:
+ pcout << "Run test case with inhomogeneous material..." << std::endl;
std::map<types::material_id, std::pair<double, double>>
inhomogeneous_material;
inhomogeneous_material[0] = std::make_pair(1.0, 1.0);
projects / dealii.git / commitdiff