372 const std::pair<unsigned int, unsigned int>
face_range) {
385 for (
unsigned int q = 0;
q <
phi_m.n_q_points; ++
q)
447 * <a name=
"step_89-CommProg"></a>
451 * <a name=
"step_89-Includefiles"></a>
503 * <a name=
"step_89-Initialconditionsforvibratingmembrane"></a>
511 * (components 1
to 1 + dim).
535 *
const double modes)
539 *
static_assert(dim == 2,
"Only implemented for dim==2");
545 *
const unsigned int comp)
const
564 * <a name=
"step_89-Gausspulse"></a>
594 *
static_assert(dim == 2,
"Only implemented for dim==2");
599 *
const unsigned int comp)
const
611 * <a name=
"step_89-Helperfunctions"></a>
642 *
template <
int dim,
typename Number,
typename VectorType>
648 *
matrix_free.get_dof_handler(),
670 *
template <
typename VectorType,
int dim>
671 *
void write_vtu(
const VectorType &solution,
674 *
const unsigned int degree,
679 *
flags.write_higher_order_cells =
true;
680 *
data_out.set_flags(flags);
682 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
685 *
std::vector<std::string> names(dim + 1,
"U");
690 *
data_out.add_data_vector(dof_handler, solution, names,
interpretation);
693 *
data_out.write_vtu_in_parallel(
name_prefix +
".vtu",
694 *
dof_handler.get_mpi_communicator());
701 * <a name=
"step_89-Materialparameterdescription"></a>
731 *
template <
typename Number>
743 *
ExcMessage(
"No materials given to CellwiseMaterialData."));
752 *
const auto n_cell_batches =
753 *
matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches();
756 *
density.resize(n_cell_batches);
758 *
for (
unsigned int cell = 0; cell < n_cell_batches; ++cell)
762 *
for (
unsigned int v = 0;
763 *
v < matrix_free.n_active_entries_per_cell_batch(cell);
767 *
matrix_free.get_cell_iterator(cell, v)->material_id();
873 *
void reinit_face(
const unsigned int face)
917 * <a name=
"step_89-Boundaryconditions"></a>
928 *
template <
int dim,
typename Number>
937 *
get_value(
const unsigned int q)
const
948 *
We don't have to apply boundary conditions for the velocity, i.e.
949 * @f$\mathbf{u}_h^+=\mathbf{u}_h^-@f$.
952 * template <int dim, typename Number>
953 * class BCEvaluationU
956 * BCEvaluationU(const FEFaceEvaluation<dim, -1, 0, dim, Number> &velocity_m)
957 * : velocity_m(velocity_m)
960 * typename FEFaceEvaluation<dim, -1, 0, dim, Number>::value_type
961 * get_value(const unsigned int q) const
963 * return velocity_m.get_value(q);
967 * const FEFaceEvaluation<dim, -1, 0, dim, Number> &velocity_m;
973 * <a name="step_89-Acousticoperator"></a>
974 * <h3>Acoustic operator</h3>
978 * The following class then defines the acoustic operator. The class is
979 * heavily based on matrix-free methods. For a better understanding in
980 * matrix-free methods please refer to @ref step_67 "step-67".
984 * At the top we define a flag that decides whether we want to use
985 * mortaring. If the remote evaluators are set up with a
986 * VectorizedArray we are using point-to-point interpolation;
987 * otherwise we make use of Nitsche-type mortaring. The decision is
988 * made using `std::is_floating_point_v`, which is a variable that
989 * is true if the template argument is a floating point type (such
990 * as `float` or `double`) and false otherwise (in particular if the
991 * template argument is a vectorized array type).
994 * template <int dim, typename Number, typename remote_value_type>
995 * class AcousticOperator
997 * static constexpr bool use_mortaring =
998 * std::is_floating_point_v<remote_value_type>;
1003 * In case of Nitsche-type mortaring, `NonMatching::MappingInfo` has to
1004 * be provided in the constructor.
1008 * const MatrixFree<dim, Number> &matrix_free,
1009 * std::shared_ptr<CellwiseMaterialData<Number>> material_data,
1010 * const std::set<types::boundary_id> &remote_face_ids,
1011 * std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
1013 * std::shared_ptr<FERemoteEvaluation<dim, dim, remote_value_type>>
1015 * std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>> c_r_eval,
1016 * std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>> rho_r_eval,
1017 * std::shared_ptr<NonMatching::MappingInfo<dim, dim, Number>> nm_info =
1019 * : matrix_free(matrix_free)
1020 * , material_data(material_data)
1021 * , remote_face_ids(remote_face_ids)
1022 * , pressure_r_eval(pressure_r_eval)
1023 * , velocity_r_eval(velocity_r_eval)
1024 * , c_r_eval(c_r_eval)
1025 * , rho_r_eval(rho_r_eval)
1026 * , nm_mapping_info(nm_info)
1028 * if (use_mortaring)
1031 * "In case of Nitsche-type mortaring NonMatching::MappingInfo \
1032 * has to be provided."));
1037 * The following function then evaluates the acoustic operator.
1038 * It first updates the precomputed values in corresponding the
1039 * FERemoteEvaluation objects. The material parameters do not change and
1040 * thus, we do not have to update precomputed values in @c c_r_eval and @c
1045 * It then either performs a matrix-free loop with Nitsche-type
1046 * mortaring at non-matching faces (if `use_mortaring` is true) or
1047 * with point-to-point interpolation at non-matching faces (in the
1048 * `else` branch). The difference is only in the third argument to
1049 * the loop function, denoting the function object that is
1050 * executed at boundary faces.
1053 * template <typename VectorType>
1054 * void evaluate(VectorType &dst, const VectorType &src) const
1056 * pressure_r_eval->gather_evaluate(src, EvaluationFlags::values);
1057 * velocity_r_eval->gather_evaluate(src, EvaluationFlags::values);
1059 * if constexpr (use_mortaring)
1062 * &AcousticOperator::local_apply_cell<VectorType>,
1063 * &AcousticOperator::local_apply_face<VectorType>,
1064 * &AcousticOperator::local_apply_boundary_face_mortaring<VectorType>,
1069 * MatrixFree<dim, Number>::DataAccessOnFaces::values,
1070 * MatrixFree<dim, Number>::DataAccessOnFaces::values);
1075 * &AcousticOperator::local_apply_cell<VectorType>,
1076 * &AcousticOperator::local_apply_face<VectorType>,
1077 * &AcousticOperator::local_apply_boundary_face_point_to_point<
1083 * MatrixFree<dim, Number>::DataAccessOnFaces::values,
1084 * MatrixFree<dim, Number>::DataAccessOnFaces::values);
1090 * In the `private` section of the class, we then define the
1091 * functions that evaluate volume, interior face, and boundary
1092 * face integrals. The concrete terms these functions evaluate are
1093 * stated in the documentation at the top of this tutorial
1094 * program. Each of these functions has its own object of type
1095 * `MaterialEvaluation` that provides access to the material at
1096 * each cell or face.
1100 * template <typename VectorType>
1101 * void local_apply_cell(
1102 * const MatrixFree<dim, Number> &matrix_free,
1104 * const VectorType &src,
1105 * const std::pair<unsigned int, unsigned int> &cell_range) const
1107 * FEEvaluation<dim, -1, 0, 1, Number> pressure(matrix_free, 0, 0, 0);
1108 * FEEvaluation<dim, -1, 0, dim, Number> velocity(matrix_free, 0, 0, 1);
1110 * MaterialEvaluation material(matrix_free, *material_data);
1112 * for (unsigned int cell = cell_range.first; cell < cell_range.second;
1115 * velocity.reinit(cell);
1116 * pressure.reinit(cell);
1118 * pressure.gather_evaluate(src, EvaluationFlags::gradients);
1119 * velocity.gather_evaluate(src, EvaluationFlags::gradients);
1123 * Get the materials on the corresponding cell. Since we
1124 * introduced @c MaterialEvaluation we can write the code
1125 * independent of whether the material is homogeneous or
1129 * material.reinit_cell(cell);
1130 * const auto c = material.get_speed_of_sound();
1131 * const auto rho = material.get_density();
1133 * for (const unsigned int q : pressure.quadrature_point_indices())
1135 * pressure.submit_value(rho * c * c * velocity.get_divergence(q),
1137 * velocity.submit_value(1.0 / rho * pressure.get_gradient(q), q);
1140 * pressure.integrate_scatter(EvaluationFlags::values, dst);
1141 * velocity.integrate_scatter(EvaluationFlags::values, dst);
1147 * This next function evaluates the fluxes at faces between cells with the
1148 * same material. If boundary faces are under consideration fluxes into
1149 * neighboring faces do not have to be considered which is enforced via
1150 * `weight_neighbor=false`. For non-matching faces the fluxes into
1151 * neighboring faces are not considered as well. This is because we iterate
1152 * over each side of the non-matching face separately (similar to a cell
1157 * In this and following functions, we also introduce the factors
1158 * @f$\tau@f$ and @f$\gamma@f$ that are computed from @f$\rho@f$ and @f$c@f$ along
1159 * interfaces and that appear in the bilinear forms. Their
1160 * definitions are provided in the introduction.
1163 * template <bool weight_neighbor,
1164 * typename InternalFaceIntegratorPressure,
1165 * typename InternalFaceIntegratorVelocity,
1166 * typename ExternalFaceIntegratorPressure,
1167 * typename ExternalFaceIntegratorVelocity>
1168 * inline DEAL_II_ALWAYS_INLINE void evaluate_face_kernel(
1169 * InternalFaceIntegratorPressure &pressure_m,
1170 * InternalFaceIntegratorVelocity &velocity_m,
1171 * ExternalFaceIntegratorPressure &pressure_p,
1172 * ExternalFaceIntegratorVelocity &velocity_p,
1173 * const typename InternalFaceIntegratorPressure::value_type c,
1174 * const typename InternalFaceIntegratorPressure::value_type rho) const
1176 * const auto tau = 0.5 * rho * c;
1177 * const auto gamma = 0.5 / (rho * c);
1179 * for (const unsigned int q : pressure_m.quadrature_point_indices())
1181 * const auto n = pressure_m.normal_vector(q);
1182 * const auto pm = pressure_m.get_value(q);
1183 * const auto um = velocity_m.get_value(q);
1185 * const auto pp = pressure_p.get_value(q);
1186 * const auto up = velocity_p.get_value(q);
1190 * Compute homogeneous local Lax-Friedrichs fluxes and submit the
1191 * corresponding values to the integrators.
1194 * const auto momentum_flux =
1195 * 0.5 * (pm + pp) + 0.5 * tau * (um - up) * n;
1196 * velocity_m.submit_value(1.0 / rho * (momentum_flux - pm) * n, q);
1197 * if constexpr (weight_neighbor)
1198 * velocity_p.submit_value(1.0 / rho * (momentum_flux - pp) * (-n), q);
1200 * const auto mass_flux = 0.5 * (um + up) + 0.5 * gamma * (pm - pp) * n;
1201 * pressure_m.submit_value(rho * c * c * (mass_flux - um) * n, q);
1202 * if constexpr (weight_neighbor)
1203 * pressure_p.submit_value(rho * c * c * (mass_flux - up) * (-n), q);
1209 * This function evaluates the fluxes at faces between cells with different
1210 * materials. This can only happen over non-matching interfaces. Therefore,
1211 * it is implicitly known that `weight_neighbor=false` and we can omit the
1215 * template <typename InternalFaceIntegratorPressure,
1216 * typename InternalFaceIntegratorVelocity,
1217 * typename ExternalFaceIntegratorPressure,
1218 * typename ExternalFaceIntegratorVelocity,
1219 * typename MaterialIntegrator>
1220 * void evaluate_face_kernel_inhomogeneous(
1221 * InternalFaceIntegratorPressure &pressure_m,
1222 * InternalFaceIntegratorVelocity &velocity_m,
1223 * const ExternalFaceIntegratorPressure &pressure_p,
1224 * const ExternalFaceIntegratorVelocity &velocity_p,
1225 * const typename InternalFaceIntegratorPressure::value_type c,
1226 * const typename InternalFaceIntegratorPressure::value_type rho,
1227 * const MaterialIntegrator &c_r,
1228 * const MaterialIntegrator &rho_r) const
1230 * const auto tau_m = 0.5 * rho * c;
1231 * const auto gamma_m = 0.5 / (rho * c);
1233 * for (const unsigned int q : pressure_m.quadrature_point_indices())
1235 * const auto c_p = c_r.get_value(q);
1236 * const auto rho_p = rho_r.get_value(q);
1237 * const auto tau_p = 0.5 * rho_p * c_p;
1238 * const auto gamma_p = 0.5 / (rho_p * c_p);
1239 * const auto tau_sum_inv = 1.0 / (tau_m + tau_p);
1240 * const auto gamma_sum_inv = 1.0 / (gamma_m + gamma_p);
1242 * const auto n = pressure_m.normal_vector(q);
1243 * const auto pm = pressure_m.get_value(q);
1244 * const auto um = velocity_m.get_value(q);
1246 * const auto pp = pressure_p.get_value(q);
1247 * const auto up = velocity_p.get_value(q);
1252 * Compute inhomogeneous fluxes and submit the corresponding values
1253 * to the integrators.
1256 * const auto momentum_flux =
1257 * pm - tau_m * tau_sum_inv * (pm - pp) +
1258 * tau_m * tau_p * tau_sum_inv * (um - up) * n;
1259 * velocity_m.submit_value(1.0 / rho * (momentum_flux - pm) * n, q);
1262 * const auto mass_flux =
1263 * um - gamma_m * gamma_sum_inv * (um - up) +
1264 * gamma_m * gamma_p * gamma_sum_inv * (pm - pp) * n;
1266 * pressure_m.submit_value(rho * c * c * (mass_flux - um) * n, q);
1272 * This function evaluates the inner face integrals.
1275 * template <typename VectorType>
1276 * void local_apply_face(
1277 * const MatrixFree<dim, Number> &matrix_free,
1279 * const VectorType &src,
1280 * const std::pair<unsigned int, unsigned int> &face_range) const
1282 * FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
1283 * matrix_free, true, 0, 0, 0);
1284 * FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_p(
1285 * matrix_free, false, 0, 0, 0);
1286 * FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
1287 * matrix_free, true, 0, 0, 1);
1288 * FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_p(
1289 * matrix_free, false, 0, 0, 1);
1291 * MaterialEvaluation material(matrix_free, *material_data);
1293 * for (unsigned int face = face_range.first; face < face_range.second;
1296 * velocity_m.reinit(face);
1297 * velocity_p.reinit(face);
1299 * pressure_m.reinit(face);
1300 * pressure_p.reinit(face);
1302 * pressure_m.gather_evaluate(src, EvaluationFlags::values);
1303 * pressure_p.gather_evaluate(src, EvaluationFlags::values);
1305 * velocity_m.gather_evaluate(src, EvaluationFlags::values);
1306 * velocity_p.gather_evaluate(src, EvaluationFlags::values);
1308 * material.reinit_face(face);
1309 * evaluate_face_kernel<true>(pressure_m,
1313 * material.get_speed_of_sound(),
1314 * material.get_density());
1316 * pressure_m.integrate_scatter(EvaluationFlags::values, dst);
1317 * pressure_p.integrate_scatter(EvaluationFlags::values, dst);
1318 * velocity_m.integrate_scatter(EvaluationFlags::values, dst);
1319 * velocity_p.integrate_scatter(EvaluationFlags::values, dst);
1327 * <a name="step_89-Matrixfreeboundaryfunctionforpointtopointinterpolation"></a>
1328 * <h4>Matrix-free boundary function for point-to-point interpolation</h4>
1332 * The following function then evaluates the boundary face integrals and the
1333 * non-matching face integrals using point-to-point interpolation.
1336 * template <typename VectorType>
1337 * void local_apply_boundary_face_point_to_point(
1338 * const MatrixFree<dim, Number> &matrix_free,
1340 * const VectorType &src,
1341 * const std::pair<unsigned int, unsigned int> &face_range) const
1343 * FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
1344 * matrix_free, true, 0, 0, 0);
1345 * FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
1346 * matrix_free, true, 0, 0, 1);
1348 * BCEvaluationP pressure_bc(pressure_m);
1349 * BCEvaluationU velocity_bc(velocity_m);
1351 * MaterialEvaluation material(matrix_free, *material_data);
1355 * Remote evaluators.
1358 * auto pressure_r = pressure_r_eval->get_data_accessor();
1359 * auto velocity_r = velocity_r_eval->get_data_accessor();
1360 * auto c_r = c_r_eval->get_data_accessor();
1361 * auto rho_r = rho_r_eval->get_data_accessor();
1363 * for (unsigned int face = face_range.first; face < face_range.second;
1366 * velocity_m.reinit(face);
1367 * pressure_m.reinit(face);
1369 * pressure_m.gather_evaluate(src, EvaluationFlags::values);
1370 * velocity_m.gather_evaluate(src, EvaluationFlags::values);
1372 * if (HelperFunctions::is_non_matching_face(
1373 * remote_face_ids, matrix_free.get_boundary_id(face)))
1377 * If @c face is non-matching we have to query values via the
1378 * FERemoteEvaluaton objects. This is done by passing the
1379 * corresponding FERemoteEvaluaton objects to the function that
1380 * evaluates the kernel. As mentioned above, each side of the
1381 * non-matching interface is traversed separately and we do not
1382 * have to consider the neighbor in the kernel. Note, that the
1383 * values in the FERemoteEvaluaton objects are already updated at
1388 * For point-to-point interpolation we simply use the
1389 * corresponding FERemoteEvaluaton objects in combination with the
1390 * standard FEFaceEvaluation objects.
1393 * velocity_r.reinit(face);
1394 * pressure_r.reinit(face);
1396 * material.reinit_face(face);
1400 * If we are considering a homogeneous material, do not use the
1401 * inhomogeneous fluxes. While it would be possible
1402 * to use the inhomogeneous fluxes they are more expensive to
1406 * if (material.is_homogeneous())
1408 * evaluate_face_kernel<false>(pressure_m,
1412 * material.get_speed_of_sound(),
1413 * material.get_density());
1418 * rho_r.reinit(face);
1419 * evaluate_face_kernel_inhomogeneous(
1424 * material.get_speed_of_sound(),
1425 * material.get_density(),
1434 * If @c face is a standard boundary face, evaluate the integral
1435 * as usual in the matrix free context. To be able to use the same
1436 * kernel as for inner faces we pass the boundary condition
1437 * objects to the function that evaluates the kernel. As detailed
1438 * above `weight_neighbor=false`.
1441 * material.reinit_face(face);
1442 * evaluate_face_kernel<false>(pressure_m,
1446 * material.get_speed_of_sound(),
1447 * material.get_density());
1450 * pressure_m.integrate_scatter(EvaluationFlags::values, dst);
1451 * velocity_m.integrate_scatter(EvaluationFlags::values, dst);
1458 * <a name="step_89-MatrixfreeboundaryfunctionforNitschetypemortaring"></a>
1459 * <h4>Matrix-free boundary function for Nitsche-type mortaring</h4>
1463 * This function evaluates the boundary face integrals and the
1464 * non-matching face integrals using Nitsche-type mortaring.
1467 * template <typename VectorType>
1468 * void local_apply_boundary_face_mortaring(
1469 * const MatrixFree<dim, Number> &matrix_free,
1471 * const VectorType &src,
1472 * const std::pair<unsigned int, unsigned int> &face_range) const
1474 * FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
1475 * matrix_free, true, 0, 0, 0);
1476 * FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
1477 * matrix_free, true, 0, 0, 1);
1481 * For Nitsche-type mortaring we are evaluating the integrals over
1482 * intersections. This is why, quadrature points are arbitrarily
1483 * distributed on every face. Thus, we can not make use of face batches
1484 * and FEFaceEvaluation but have to consider each face individually and
1485 * make use of @c FEFacePointEvaluation to evaluate the integrals in the
1486 * arbitrarily distributed quadrature points.
1487 * Since the setup of FEFacePointEvaluation is more expensive than that of
1488 * FEEvaluation we do the setup only once. For this we are using the
1489 * helper function @c get_thread_safe_fe_face_point_evaluation_object().
1492 * FEFacePointEvaluation<1, dim, dim, Number> &pressure_m_mortar =
1493 * get_thread_safe_fe_face_point_evaluation_object<1>(
1494 * thread_local_pressure_m_mortar, 0);
1495 * FEFacePointEvaluation<dim, dim, dim, Number> &velocity_m_mortar =
1496 * get_thread_safe_fe_face_point_evaluation_object<dim>(
1497 * thread_local_velocity_m_mortar, 1);
1499 * BCEvaluationP pressure_bc(pressure_m);
1500 * BCEvaluationU velocity_bc(velocity_m);
1502 * MaterialEvaluation material(matrix_free, *material_data);
1504 * auto pressure_r_mortar = pressure_r_eval->get_data_accessor();
1505 * auto velocity_r_mortar = velocity_r_eval->get_data_accessor();
1506 * auto c_r = c_r_eval->get_data_accessor();
1507 * auto rho_r = rho_r_eval->get_data_accessor();
1509 * for (unsigned int face = face_range.first; face < face_range.second;
1512 * if (HelperFunctions::is_non_matching_face(
1513 * remote_face_ids, matrix_free.get_boundary_id(face)))
1515 * material.reinit_face(face);
1519 * First fetch the DoF values with standard FEFaceEvaluation
1523 * pressure_m.reinit(face);
1524 * velocity_m.reinit(face);
1526 * pressure_m.read_dof_values(src);
1527 * velocity_m.read_dof_values(src);
1531 * Project the internally stored values into the face DoFs
1532 * of the current face.
1535 * pressure_m.project_to_face(EvaluationFlags::values);
1536 * velocity_m.project_to_face(EvaluationFlags::values);
1540 * For mortaring, we have to consider every face from the face
1541 * batches separately and have to use the FEFacePointEvaluation
1542 * objects to be able to evaluate the integrals with the
1543 * arbitrarily distributed quadrature points.
1546 * for (unsigned int v = 0;
1547 * v < matrix_free.n_active_entries_per_face_batch(face);
1550 * constexpr unsigned int n_lanes =
1551 * VectorizedArray<Number>::size();
1552 * velocity_m_mortar.reinit(face * n_lanes + v);
1553 * pressure_m_mortar.reinit(face * n_lanes + v);
1557 * Evaluate using FEFacePointEvaluation. As buffer,
1558 * simply use the internal buffers from the
1559 * FEFaceEvaluation objects.
1562 * velocity_m_mortar.evaluate_in_face(
1563 * &velocity_m.get_scratch_data().begin()[0][v],
1564 * EvaluationFlags::values);
1566 * pressure_m_mortar.evaluate_in_face(
1567 * &pressure_m.get_scratch_data().begin()[0][v],
1568 * EvaluationFlags::values);
1570 * velocity_r_mortar.reinit(face * n_lanes + v);
1571 * pressure_r_mortar.reinit(face * n_lanes + v);
1575 * As above, if we are considering a homogeneous
1576 * material, do not use the inhomogeneous
1577 * fluxes. Since we are operating on face @c v we
1578 * call @c material.get_density()[v].
1581 * if (material.is_homogeneous())
1583 * evaluate_face_kernel<false>(
1584 * pressure_m_mortar,
1585 * velocity_m_mortar,
1586 * pressure_r_mortar,
1587 * velocity_r_mortar,
1588 * material.get_speed_of_sound()[v],
1589 * material.get_density()[v]);
1593 * c_r.reinit(face * n_lanes + v);
1594 * rho_r.reinit(face * n_lanes + v);
1596 * evaluate_face_kernel_inhomogeneous(
1597 * pressure_m_mortar,
1598 * velocity_m_mortar,
1599 * pressure_r_mortar,
1600 * velocity_r_mortar,
1601 * material.get_speed_of_sound()[v],
1602 * material.get_density()[v],
1609 * Integrate using FEFacePointEvaluation. As buffer,
1610 * simply use the internal buffers from the
1611 * FEFaceEvaluation objects.
1614 * velocity_m_mortar.integrate_in_face(
1615 * &velocity_m.get_scratch_data().begin()[0][v],
1616 * EvaluationFlags::values);
1618 * pressure_m_mortar.integrate_in_face(
1619 * &pressure_m.get_scratch_data().begin()[0][v],
1620 * EvaluationFlags::values);
1625 * Collect the contributions from the face DoFs to
1626 * the internal cell DoFs to be able to use the
1627 * member function @c distribute_local_to_global().
1630 * pressure_m.collect_from_face(EvaluationFlags::values);
1631 * velocity_m.collect_from_face(EvaluationFlags::values);
1633 * pressure_m.distribute_local_to_global(dst);
1634 * velocity_m.distribute_local_to_global(dst);
1638 * velocity_m.reinit(face);
1639 * pressure_m.reinit(face);
1641 * pressure_m.gather_evaluate(src, EvaluationFlags::values);
1642 * velocity_m.gather_evaluate(src, EvaluationFlags::values);
1644 * material.reinit_face(face);
1645 * evaluate_face_kernel<false>(pressure_m,
1649 * material.get_speed_of_sound(),
1650 * material.get_density());
1652 * pressure_m.integrate_scatter(EvaluationFlags::values, dst);
1653 * velocity_m.integrate_scatter(EvaluationFlags::values, dst);
1658 * const MatrixFree<dim, Number> &matrix_free;
1660 * const std::shared_ptr<CellwiseMaterialData<Number>> material_data;
1662 * const std::set<types::boundary_id> remote_face_ids;
1666 * FERemoteEvaluation objects are stored as shared pointers. This way,
1667 * they can also be used for other operators without precomputing the values
1671 * const std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
1673 * const std::shared_ptr<FERemoteEvaluation<dim, dim, remote_value_type>>
1676 * const std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
1678 * const std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
1681 * const std::shared_ptr<NonMatching::MappingInfo<dim, dim, Number>>
1686 * We store FEFacePointEvaluation objects as members in a thread local
1687 * way, since its creation is more expensive compared to FEEvaluation
1691 * mutable Threads::ThreadLocalStorage<
1692 * std::unique_ptr<FEFacePointEvaluation<1, dim, dim, Number>>>
1693 * thread_local_pressure_m_mortar;
1695 * mutable Threads::ThreadLocalStorage<
1696 * std::unique_ptr<FEFacePointEvaluation<dim, dim, dim, Number>>>
1697 * thread_local_velocity_m_mortar;
1701 * Helper function to create and get FEFacePointEvaluation objects in a
1702 * thread safe way. On each thread, FEFacePointEvaluation is created if it
1703 * has not been created by now. After that, simply return the object
1704 * corresponding to the thread under consideration.
1707 * template <int n_components>
1708 * FEFacePointEvaluation<n_components, dim, dim, Number> &
1709 * get_thread_safe_fe_face_point_evaluation_object(
1710 * Threads::ThreadLocalStorage<
1711 * std::unique_ptr<FEFacePointEvaluation<n_components, dim, dim, Number>>>
1712 * &fe_face_point_eval_thread_local,
1713 * unsigned int fist_selected_comp) const
1715 * if (fe_face_point_eval_thread_local.get() == nullptr)
1717 * fe_face_point_eval_thread_local = std::make_unique<
1718 * FEFacePointEvaluation<n_components, dim, dim, Number>>(
1720 * matrix_free.get_dof_handler().get_fe(),
1722 * fist_selected_comp);
1724 * return *fe_face_point_eval_thread_local.get();
1731 * <a name="step_89-Inversemassoperator"></a>
1732 * <h3>Inverse mass operator</h3>
1736 * For the time stepping methods below, we need the inverse mass
1737 * operator. We apply it via a loop over all (batches of) cells as
1738 * always, where the contribution of each cell is computed in a
1742 * template <int dim, typename Number>
1743 * class InverseMassOperator
1746 * InverseMassOperator(const MatrixFree<dim, Number> &matrix_free)
1747 * : matrix_free(matrix_free)
1750 * template <typename VectorType>
1751 * void apply(VectorType &dst, const VectorType &src) const
1753 * dst.zero_out_ghost_values();
1754 * matrix_free.cell_loop(&InverseMassOperator::local_apply_cell<VectorType>,
1761 * template <typename VectorType>
1762 * void local_apply_cell(
1763 * const MatrixFree<dim, Number> &mf,
1765 * const VectorType &src,
1766 * const std::pair<unsigned int, unsigned int> &cell_range) const
1768 * FEEvaluation<dim, -1, 0, dim + 1, Number> phi(mf);
1769 * MatrixFreeOperators::CellwiseInverseMassMatrix<dim, -1, dim + 1, Number>
1772 * for (unsigned int cell = cell_range.first; cell < cell_range.second;
1776 * phi.read_dof_values(src);
1777 * minv.apply(phi.begin_dof_values(), phi.begin_dof_values());
1778 * phi.set_dof_values(dst);
1782 * const MatrixFree<dim, Number> &matrix_free;
1788 * <a name="step_89-RungeKuttatimestepping"></a>
1789 * <h3>Runge-Kutta time-stepping</h3>
1793 * This class implements a Runge-Kutta scheme of order 2.
1796 * template <int dim, typename Number, typename remote_value_type>
1799 * using VectorType = LinearAlgebra::distributed::Vector<Number>;
1803 * const std::shared_ptr<InverseMassOperator<dim, Number>>
1804 * inverse_mass_operator,
1805 * const std::shared_ptr<AcousticOperator<dim, Number, remote_value_type>>
1806 * acoustic_operator)
1807 * : inverse_mass_operator(inverse_mass_operator)
1808 * , acoustic_operator(acoustic_operator)
1813 * The `run()` function of this class contains the time loop. It
1814 * starts by initializing some member variables (such as
1815 * short-cuts to objects that describe the finite element, its
1816 * properties, and the mapping; an by initializing vectors). It
1817 * then computes the time step size via minimum element edge
1818 * length. Assuming non-distorted elements, we can compute the
1819 * edge length as the distance between two vertices. From this,
1820 * we can then obtain the time step size via the CFL condition.
1823 * void run(const MatrixFree<dim, Number> &matrix_free,
1825 * const double end_time,
1826 * const double speed_of_sound,
1827 * const Function<dim> &initial_condition,
1828 * const std::string &vtk_prefix)
1830 * const auto &dof_handler = matrix_free.get_dof_handler();
1831 * const auto &mapping = *matrix_free.get_mapping_info().mapping;
1832 * const auto degree = dof_handler.get_fe().degree;
1834 * VectorType solution;
1835 * matrix_free.initialize_dof_vector(solution);
1836 * VectorType solution_temp;
1837 * matrix_free.initialize_dof_vector(solution_temp);
1839 * HelperFunctions::set_initial_condition(matrix_free,
1840 * initial_condition,
1843 * double h_local_min = std::numeric_limits<double>::max();
1844 * for (const auto &cell : dof_handler.active_cell_iterators())
1846 * std::min(h_local_min,
1847 * (cell->vertex(1) - cell->vertex(0)).norm_square());
1848 * h_local_min = std::sqrt(h_local_min);
1849 * const double h_min =
1850 * Utilities::MPI::min(h_local_min, dof_handler.get_mpi_communicator());
1853 * cr * HelperFunctions::compute_dt_cfl(h_min, degree, speed_of_sound);
1857 * We can then perform the time integration loop:
1860 * double time = 0.0;
1861 * unsigned int timestep = 0;
1862 * while (time < end_time)
1864 * HelperFunctions::write_vtu(solution,
1865 * matrix_free.get_dof_handler(),
1868 * "step_89-" + vtk_prefix +
1869 * std::to_string(timestep));
1871 * std::swap(solution, solution_temp);
1874 * perform_time_step(dt, solution, solution_temp);
1880 * The main work of this class is done by a `private` member
1881 * function that performs one Runge-Kutta 2 time step. Recall that
1882 * this method requires two sub-steps ("stages") computing
1883 * intermediate values `k1` and `k2`, after which the intermediate
1884 * values are summed with weights to obtain the new solution at
1885 * the end of the time step. The RK2 method allows for the
1886 * elimination of the intermediate vector `k2` by utilizing the
1887 * `dst` vector as temporary storage.
1892 * perform_time_step(const double dt, VectorType &dst, const VectorType &src)
1894 * VectorType k1 = src;
1896 * /* First stage. */
1897 * evaluate_stage(k1, src);
1899 * /* Second stage. */
1900 * k1.sadd(0.5 * dt, 1.0, src);
1901 * evaluate_stage(dst, k1);
1903 * /* Summing things into the output vector. */
1904 * dst.sadd(dt, 1.0, src);
1909 * Evaluating a single Runge-Kutta stage is a straightforward step
1910 * that really only requires applying the operator once, and then
1911 * applying the inverse of the mass matrix.
1914 * void evaluate_stage(VectorType &dst, const VectorType &src)
1916 * acoustic_operator->evaluate(dst, src);
1918 * inverse_mass_operator->apply(dst, dst);
1921 * const std::shared_ptr<InverseMassOperator<dim, Number>>
1922 * inverse_mass_operator;
1923 * const std::shared_ptr<AcousticOperator<dim, Number, remote_value_type>>
1924 * acoustic_operator;
1931 * <a name="step_89-Constructionofnonmatchingtriangulations"></a>
1932 * <h3>Construction of non-matching triangulations</h3>
1936 * Let us now make our way to the higher-level functions of this program.
1940 * The first of these functions creates a two dimensional square
1941 * triangulation that spans from @f$(0,0)@f$ to @f$(1,1)@f$. It consists of
1942 * two sub-domains. The left sub-domain spans from @f$(0,0)@f$ to
1943 * @f$(0.525,1)@f$. The right sub-domain spans from @f$(0.525,0)@f$ to
1944 * @f$(1,1)@f$. The left sub-domain has elements that are three times
1945 * smaller compared to the ones for the right sub-domain.
1949 * At non-matching interfaces, we need to provide different boundary
1950 * IDs for the cells that make up the two parts (because, while they
1951 * may be physically adjacent, they are not logically neighbors
1952 * given that the faces of cells on both sides do not match, and the
1953 * Triangulation class will therefore treat the interface between
1954 * the two parts as a "boundary"). These boundary IDs have to differ
1955 * because later on RemotePointEvaluation has to search for remote
1956 * points for each face, that are defined in the same mesh (since we
1957 * merge the mesh) but not on the same side of the non-matching
1958 * interface. As a consequence, we declare at the top symbolic names
1959 * for these boundary indicators, and ensure that we return a set
1960 * with these values to the caller for later use.
1964 * The actual mesh is then constructed by first constructing the
1965 * left and right parts separately (setting material ids to zero and
1966 * one, respectively), and using the appropriate boundary ids for
1967 * all parts of the mesh. We then use
1968 * GridGenerator::merge_triangulations() to combine them into one
1969 * (non-matching) mesh. We have to pay attention that should the two
1970 * sub-triangulations have vertices at the same locations, that they
1971 * are not merged (connecting the two triangulations logically)
1972 * since we want the interface to be an actual boundary. We achieve
1973 * this by providing a tolerance of zero for the merge, see the
1974 * documentation of the function
1975 * GridGenerator::merge_triangulations().
1978 * template <int dim>
1979 * void build_non_matching_triangulation(
1980 * Triangulation<dim> &tria,
1981 * std::set<types::boundary_id> &non_matching_faces,
1982 * const unsigned int refinements)
1984 * const double length = 1.0;
1986 * const types::boundary_id non_matching_id_left = 98;
1987 * const types::boundary_id non_matching_id_right = 99;
1989 * non_matching_faces = {non_matching_id_left, non_matching_id_right};
1991 * /* Construct left part of mesh. */
1992 * Triangulation<dim> tria_left;
1993 * const unsigned int subdiv_left = 11;
1994 * GridGenerator::subdivided_hyper_rectangle(tria_left,
1995 * {subdiv_left, 2 * subdiv_left},
1997 * {0.525 * length, length});
1999 * for (const auto &cell : tria_left.active_cell_iterators())
2000 * cell->set_material_id(0);
2001 * for (const auto &face : tria_left.active_face_iterators())
2002 * if (face->at_boundary())
2004 * face->set_boundary_id(0);
2005 * if (face->center()[0] > 0.525 * length - 1e-6)
2006 * face->set_boundary_id(non_matching_id_left);
2009 * /* Construct right part of mesh. */
2010 * Triangulation<dim> tria_right;
2011 * const unsigned int subdiv_right = 4;
2012 * GridGenerator::subdivided_hyper_rectangle(tria_right,
2013 * {subdiv_right, 2 * subdiv_right},
2014 * {0.525 * length, 0.0},
2015 * {length, length});
2017 * for (const auto &cell : tria_right.active_cell_iterators())
2018 * cell->set_material_id(1);
2019 * for (const auto &face : tria_right.active_face_iterators())
2020 * if (face->at_boundary())
2022 * face->set_boundary_id(0);
2023 * if (face->center()[0] < 0.525 * length + 1e-6)
2024 * face->set_boundary_id(non_matching_id_right);
2027 * /* Merge triangulations. */
2028 * GridGenerator::merge_triangulations(tria_left,
2032 * /*copy_manifold_ids*/ false,
2033 * /*copy_boundary_ids*/ true);
2035 * /* Refine the result. */
2036 * tria.refine_global(refinements);
2042 * <a name="step_89-Setupandrunningofthetwoschemes"></a>
2043 * <h3>Set up and running of the two schemes</h3>
2048 * <a name="step_89-Setupandrunningofthepointtopointinterpolationscheme"></a>
2049 * <h4>Setup and running of the point-to-point interpolation scheme</h4>
2053 * We are now at the two functions that run the overall schemes (the
2054 * point-to-point and the mortaring schemes). The first of these
2055 * functions fills a FERemoteEvaluationCommunicator object that is
2056 * needed for point-to-point interpolation. Additionally, the
2057 * corresponding remote evaluators are set up using this remote
2058 * communicator. Eventually, the operators are handed to the time
2059 * integrator that runs the simulation.
2062 * template <int dim, typename Number>
2063 * void run_with_point_to_point_interpolation(
2064 * const MatrixFree<dim, Number> &matrix_free,
2065 * const std::set<types::boundary_id> &non_matching_faces,
2066 * const std::map<types::material_id, std::pair<double, double>> &materials,
2067 * const double end_time,
2068 * const Function<dim> &initial_condition,
2069 * const std::string &vtk_prefix)
2071 * const auto &dof_handler = matrix_free.get_dof_handler();
2072 * const auto &tria = dof_handler.get_triangulation();
2076 * Communication objects know about the communication pattern. That is,
2077 * they know about the cells and quadrature points that have to be
2078 * evaluated at remote faces. This information is given via
2079 * RemotePointEvaluation. Additionally, the communication objects
2080 * have to be able to match the quadrature points of the remote
2081 * points (that provide exterior information) to the quadrature points
2082 * defined at the interior cell. In case of point-to-point interpolation
2083 * a vector of pairs with face batch Ids and the number of faces in the
2084 * batch is needed. @c FERemoteCommunicationObjectEntityBatches
2085 * is a container to store this information.
2089 * The information is filled outside of the actual class since in some cases
2090 * the information is available from some heuristic and
2091 * it is possible to skip some expensive operations. This is for example
2092 * the case for sliding rotating interfaces with equally spaced elements on
2093 * both sides of the non-matching interface @cite duerrwaechter2021an.
2097 * For the standard case of point to point-to-point interpolation without
2098 * any heuristic we make use of the utility function
2099 * Utilities::compute_remote_communicator_faces_point_to_point_interpolation().
2100 * Please refer to the documentation of this function to see how to manually
2101 * set up the remote communicator from outside.
2105 * Among the inputs for the remote communicator we need a list of
2106 * boundary ids for the non-matching faces, along with a function
2107 * object for each boundary id that returns a vector of true/false
2108 * flags in which exactly the vertices of cells of the
2109 * triangulation are marked that have a face at the boundary id in
2114 * std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
2115 * non_matching_faces_marked_vertices;
2116 * for (const auto &nm_face : non_matching_faces)
2118 * auto marked_vertices = [&nm_face, &tria]() -> std::vector<bool> {
2119 * std::vector<bool> mask(tria.n_vertices(), true);
2121 * for (const auto &cell : tria.active_cell_iterators())
2122 * for (auto const &f : cell->face_indices())
2123 * if (cell->face(f)->at_boundary() &&
2124 * cell->face(f)->boundary_id() == nm_face)
2125 * for (const auto v : cell->vertex_indices())
2126 * mask[cell->vertex_index(v)] = false;
2131 * non_matching_faces_marked_vertices.emplace_back(
2132 * std::make_pair(nm_face, marked_vertices));
2135 * auto remote_communicator =
2136 * Utilities::compute_remote_communicator_faces_point_to_point_interpolation(
2137 * matrix_free, non_matching_faces_marked_vertices);
2141 * We are using point-to-point interpolation and can therefore
2142 * easily access the corresponding data at face batches. This
2143 * is why we use a @c VectorizedArray as @c remote_value_type:
2146 * using remote_value_type = VectorizedArray<Number>;
2150 * We then set up FERemoteEvaluation objects that access the
2151 * pressure and velocity at remote faces, along with an object to
2152 * describe cell-wise material data.
2155 * const auto pressure_r =
2156 * std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
2157 * remote_communicator, dof_handler, /*first_selected_component*/ 0);
2159 * const auto velocity_r =
2160 * std::make_shared<FERemoteEvaluation<dim, dim, remote_value_type>>(
2161 * remote_communicator, dof_handler, /*first_selected_component*/ 1);
2163 * const auto material_data =
2164 * std::make_shared<CellwiseMaterialData<Number>>(matrix_free, materials);
2168 * If we have an inhomogeneous problem, we have to set up the
2169 * material handler that accesses the materials at remote faces.
2173 * std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
2174 * remote_communicator,
2175 * matrix_free.get_dof_handler().get_triangulation(),
2176 * /*first_selected_component*/ 0);
2177 * const auto rho_r =
2178 * std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
2179 * remote_communicator,
2180 * matrix_free.get_dof_handler().get_triangulation(),
2181 * /*first_selected_component*/ 0);
2185 * If the domain is not homogeneous, i.e., if material parameters
2186 * change from cell to cell, we initialize and fill DoF vectors
2187 * that contain the material properties. Materials do not change
2188 * during the simulation, therefore there is no need to ever
2189 * compute the values after the first @c gather_evaluate() (at the
2190 * end of the following block) again.
2193 * if (!material_data->is_homogeneous())
2196 * matrix_free.get_dof_handler().get_triangulation().n_active_cells());
2197 * Vector<Number> rho(
2198 * matrix_free.get_dof_handler().get_triangulation().n_active_cells());
2200 * for (const auto &cell : matrix_free.get_dof_handler()
2201 * .get_triangulation()
2202 * .active_cell_iterators())
2204 * c[cell->active_cell_index()] =
2205 * materials.at(cell->material_id()).first;
2206 * rho[cell->active_cell_index()] =
2207 * materials.at(cell->material_id()).second;
2210 * c_r->gather_evaluate(c, EvaluationFlags::values);
2211 * rho_r->gather_evaluate(rho, EvaluationFlags::values);
2217 * Next, we set up the inverse mass operator and the acoustic
2218 * operator. Using `remote_value_type=VectorizedArray<Number>`
2219 * makes the operator use point-to-point interpolation. These two
2220 * objects are then used to create a `RungeKutta2` object to
2221 * perform the time integration.
2225 * We also compute the maximum speed of sound, needed for the
2226 * computation of the time-step size, and then run the time integrator.
2227 * For the examples considered here, we found a limiting Courant number of
2228 * @f$\mathrm{Cr}\approx 0.36@f$ to maintain stability. To ensure, the
2229 * error of the temporal discretization is small, we use a considerably
2230 * smaller Courant number of @f$0.2@f$.
2233 * const auto inverse_mass_operator =
2234 * std::make_shared<InverseMassOperator<dim, Number>>(matrix_free);
2236 * const auto acoustic_operator =
2237 * std::make_shared<AcousticOperator<dim, Number, remote_value_type>>(
2240 * non_matching_faces,
2246 * RungeKutta2<dim, Number, remote_value_type> time_integrator(
2247 * inverse_mass_operator, acoustic_operator);
2249 * double speed_of_sound_max = 0.0;
2250 * for (const auto &mat : materials)
2251 * speed_of_sound_max = std::max(speed_of_sound_max, mat.second.first);
2254 * time_integrator.run(matrix_free,
2257 * speed_of_sound_max,
2258 * initial_condition,
2265 * <a name="step_89-SetupandrunningoftheNitschetypemortaringscheme"></a>
2266 * <h4>Setup and running of the Nitsche-type mortaring scheme</h4>
2270 * The alternative to the previous function is to use the mortaring
2271 * scheme -- implemented in the following function. This function
2272 * can only be run when deal.II is configured using CGAL (but the
2273 * previous function can still be used without CGAL), and so errors
2274 * out if CGAL is not available.
2277 * template <int dim, typename Number>
2278 * void run_with_nitsche_type_mortaring(
2279 * const MatrixFree<dim, Number> &matrix_free,
2280 * const std::set<types::boundary_id> &non_matching_faces,
2281 * const std::map<types::material_id, std::pair<double, double>> &materials,
2282 * const double end_time,
2283 * const Function<dim> &initial_condition,
2284 * const std::string &vtk_prefix)
2286 * #ifndef DEAL_II_WITH_CGAL
2287 * (void)matrix_free;
2288 * (void)non_matching_faces;
2291 * (void)initial_condition;
2294 * ConditionalOStream pcout(
2295 * std::cout, (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0));
2297 * pcout << "In this function, mortars are computed using CGAL. "
2298 * "Configure deal.II with DEAL_II_WITH_CGAL to run this function.\n";
2303 * const auto &dof_handler = matrix_free.get_dof_handler();
2304 * const auto &tria = dof_handler.get_triangulation();
2305 * const auto &mapping = *matrix_free.get_mapping_info().mapping;
2306 * const auto n_quadrature_pnts = matrix_free.get_quadrature().size();
2309 * std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
2310 * non_matching_faces_marked_vertices;
2312 * for (const auto &nm_face : non_matching_faces)
2314 * auto marked_vertices = [&]() {
2315 * std::vector<bool> mask(tria.n_vertices(), true);
2317 * for (const auto &cell : tria.active_cell_iterators())
2318 * for (auto const &f : cell->face_indices())
2319 * if (cell->face(f)->at_boundary() &&
2320 * cell->face(f)->boundary_id() == nm_face)
2321 * for (const auto v : cell->vertex_indices())
2322 * mask[cell->vertex_index(v)] = false;
2327 * non_matching_faces_marked_vertices.emplace_back(
2328 * std::make_pair(nm_face, marked_vertices));
2333 * The only parts in this function that are functionally different
2334 * from the previous one follow here.
2338 * First, quadrature points are arbitrarily distributed on each non-matching
2339 * face. Therefore, we have to make use of FEFacePointEvaluation.
2340 * FEFacePointEvaluation needs NonMatching::MappingInfo to work at the
2341 * correct quadrature points that are in sync with used FERemoteEvaluation
2343 * `compute_remote_communicator_faces_nitsche_type_mortaring()` to reinit
2344 * NonMatching::MappingInfo ensures this. In the case of mortaring, we have
2345 * to use the weights provided by the quadrature rules that are used to set
2346 * up NonMatching::MappingInfo. Therefore we set the flag @c
2347 * use_global_weights.
2350 * typename NonMatching::MappingInfo<dim, dim, Number>::AdditionalData
2352 * additional_data.use_global_weights = true;
2356 * Set up NonMatching::MappingInfo with needed update flags and
2357 * @c additional_data.
2360 * auto nm_mapping_info =
2361 * std::make_shared<NonMatching::MappingInfo<dim, dim, Number>>(
2363 * update_values | update_JxW_values | update_normal_vectors |
2364 * update_quadrature_points,
2367 * auto remote_communicator =
2368 * Utilities::compute_remote_communicator_faces_nitsche_type_mortaring(
2370 * non_matching_faces_marked_vertices,
2371 * n_quadrature_pnts,
2373 * nm_mapping_info.get());
2377 * Second, since quadrature points are arbitrarily distributed we
2378 * have to consider each face in a batch separately and can not
2379 * make use of @c VecorizedArray.
2382 * using remote_value_type = Number;
2386 * The rest of the code is then identical to what we had in the
2387 * previous function (though it functions differently because of
2388 * the difference in `remote_value_type`).
2391 * const auto pressure_r =
2392 * std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
2393 * remote_communicator, dof_handler, /*first_selected_component*/ 0);
2395 * const auto velocity_r =
2396 * std::make_shared<FERemoteEvaluation<dim, dim, remote_value_type>>(
2397 * remote_communicator, dof_handler, /*first_selected_component*/ 1);
2399 * const auto material_data =
2400 * std::make_shared<CellwiseMaterialData<Number>>(matrix_free, materials);
2403 * std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
2404 * remote_communicator,
2405 * matrix_free.get_dof_handler().get_triangulation(),
2406 * /*first_selected_component*/ 0);
2407 * const auto rho_r =
2408 * std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
2409 * remote_communicator,
2410 * matrix_free.get_dof_handler().get_triangulation(),
2411 * /*first_selected_component*/ 0);
2413 * if (!material_data->is_homogeneous())
2416 * matrix_free.get_dof_handler().get_triangulation().n_active_cells());
2417 * Vector<Number> rho(
2418 * matrix_free.get_dof_handler().get_triangulation().n_active_cells());
2420 * for (const auto &cell : matrix_free.get_dof_handler()
2421 * .get_triangulation()
2422 * .active_cell_iterators())
2424 * c[cell->active_cell_index()] =
2425 * materials.at(cell->material_id()).first;
2426 * rho[cell->active_cell_index()] =
2427 * materials.at(cell->material_id()).second;
2430 * c_r->gather_evaluate(c, EvaluationFlags::values);
2431 * rho_r->gather_evaluate(rho, EvaluationFlags::values);
2434 * const auto inverse_mass_operator =
2435 * std::make_shared<InverseMassOperator<dim, Number>>(matrix_free);
2437 * const auto acoustic_operator =
2438 * std::make_shared<AcousticOperator<dim, Number, remote_value_type>>(
2441 * non_matching_faces,
2448 * RungeKutta2<dim, Number, remote_value_type> time_integrator(
2449 * inverse_mass_operator, acoustic_operator);
2451 * double speed_of_sound_max = 0.0;
2452 * for (const auto &mat : materials)
2453 * speed_of_sound_max = std::max(speed_of_sound_max, mat.second.first);
2455 * time_integrator.run(matrix_free,
2458 * speed_of_sound_max,
2459 * initial_condition,
2463 * } // namespace Step89
2469 * <a name="step_89-main"></a>
2474 * Finally, the `main()` function executes the different versions of handling
2475 * non-matching interfaces.
2479 * Similar to @ref step_87 "step-87", the minimum requirement of this tutorial is MPI.
2480 * The parallel::distributed::Triangulation class is used if deal.II is
2481 * configured with p4est. Otherwise parallel::shared::Triangulation is used.
2484 * int main(int argc, char *argv[])
2486 * using namespace dealii;
2487 * constexpr int dim = 2;
2488 * using Number = double;
2490 * Utilities::MPI::MPI_InitFinalize mpi(argc, argv);
2491 * std::cout.precision(5);
2492 * ConditionalOStream pcout(std::cout,
2493 * (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) ==
2496 * const unsigned int refinements = 1;
2497 * const unsigned int degree = 3;
2499 * #ifdef DEAL_II_WITH_P4EST
2500 * parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
2502 * parallel::shared::Triangulation<dim> tria(MPI_COMM_WORLD);
2505 * pcout << "Create non-matching grid..." << std::endl;
2507 * std::set<types::boundary_id> non_matching_faces;
2508 * Step89::build_non_matching_triangulation(tria,
2509 * non_matching_faces,
2512 * pcout << " - Refinement level: " << refinements << std::endl;
2513 * pcout << " - Number of cells: " << tria.n_cells() << std::endl;
2515 * pcout << "Create DoFHandler..." << std::endl;
2516 * DoFHandler<dim> dof_handler(tria);
2517 * dof_handler.distribute_dofs(FESystem<dim>(FE_DGQ<dim>(degree) ^ (dim + 1)));
2518 * pcout << " - Number of DoFs: " << dof_handler.n_dofs() << std::endl;
2520 * AffineConstraints<Number> constraints;
2521 * constraints.close();
2523 * pcout << "Set up MatrixFree..." << std::endl;
2524 * typename MatrixFree<dim, Number>::AdditionalData data;
2525 * data.mapping_update_flags = update_gradients | update_values;
2526 * data.mapping_update_flags_inner_faces = update_values;
2527 * data.mapping_update_flags_boundary_faces =
2528 * update_quadrature_points | update_values;
2530 * MatrixFree<dim, Number> matrix_free;
2531 * matrix_free.reinit(
2532 * MappingQ1<dim>(), dof_handler, constraints, QGauss<dim>(degree + 1), data);
2538 * <a name="step_89-RunvibratingmembranetestcaseHomogeneouspressure"></a>
2539 * <h4>Run vibrating membrane test case} Homogeneous pressure</h4>
2540 * Dirichlet boundary conditions are applied for
2541 * simplicity. Therefore, modes can not be chosen arbitrarily.
2544 * pcout << "Run vibrating membrane test case..." << std::endl;
2545 * const double modes = 10.0;
2546 * std::map<types::material_id, std::pair<double, double>> homogeneous_material;
2547 * homogeneous_material[numbers::invalid_material_id] = std::make_pair(1.0, 1.0);
2548 * const auto initial_solution_membrane =
2549 * Step89::InitialConditionVibratingMembrane<dim>(modes);
2551 * /* Run vibrating membrane test case using point-to-point interpolation: */
2552 * pcout << " - Point-to-point interpolation: " << std::endl;
2553 * Step89::run_with_point_to_point_interpolation(
2555 * non_matching_faces,
2556 * homogeneous_material,
2557 * 8.0 * initial_solution_membrane.get_period_duration(
2558 * homogeneous_material.begin()->second.first),
2559 * initial_solution_membrane,
2562 * /* Run vibrating membrane test case using Nitsche-type mortaring: */
2563 * pcout << " - Nitsche-type mortaring: " << std::endl;
2564 * Step89::run_with_nitsche_type_mortaring(
2566 * non_matching_faces,
2567 * homogeneous_material,
2568 * 8.0 * initial_solution_membrane.get_period_duration(
2569 * homogeneous_material.begin()->second.first),
2570 * initial_solution_membrane,
2576 * <a name="step_89-Runtestcasewithinhomogeneousmaterial"></a>
2577 * <h4>Run test case with inhomogeneous material</h4>
2578 * Run simple test case with inhomogeneous material and Nitsche-type
2582 * pcout << "Run test case with inhomogeneous material..." << std::endl;
2583 * std::map<types::material_id, std::pair<double, double>>
2584 * inhomogeneous_material;
2585 * inhomogeneous_material[0] = std::make_pair(1.0, 1.0);
2586 * inhomogeneous_material[1] = std::make_pair(3.0, 1.0);
2587 * Step89::run_with_nitsche_type_mortaring(matrix_free,
2588 * non_matching_faces,
2589 * inhomogeneous_material,
2591 * Step89::GaussPulse<dim>(0.3, 0.5),
2598<a name="step_89-Results"></a><h1>Results</h1>
2601<a name="step_89-VibratingmembranePointtopointinterpolationvsNitschetypemortaring"></a><h3>Vibrating membrane: Point-to-point interpolation vs. Nitsche-type mortaring</h3>
2604We compare the results of the simulations after the last time step, i.e. at @f$t=8T@f$.
2605The @f$y@f$-component of the velocity field using Nitsche-type mortaring is depicted on the left.
2606The same field using point-to-point interpolation is depicted on the right.
2608<table align="center" class="doxtable">
2611 @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_mortaring_velocity_Y.png "" width=60%
2614 @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_point_to_point_velocity_Y.png "" width=60%
2619Besides this, the results for the pressure and the velocity in @f$y@f$ direction
2620are plotted along the horizontal line that spans from (0,0.69) to (1,0.69).
2622<table align="center" class="doxtable">
2625 @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_mortaring_vs_point_to_point_pressure.svg "" width=150%
2628 @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_mortaring_vs_point_to_point_velocity_Y.svg "" width=150%
2633While the results of the pressure are similar, @f$u_y@f$ clearly differs. At certain
2634positions we can see aliasing errors for the point-to-point interpolation.
2635For different mesh configurations and/or longer run times, the aliasing effects
2636of the point-to-point simulation accumulate and the simulation becomes instable.
2637To keep the tutorial short we have chosen one mesh that can be used for all
2638examples. For a configuration that yields instable results for a wide range of
2639polynomial degrees, see @cite heinz2022high.
2641<a name="step_89-Wavepropagationthroughinhomogeneousfluid"></a><h3>Wave propagation through in-homogeneous fluid</h3>
2644The example that follows is just one example in which non-matching discretizations can be efficiently
2645used to reduce the number of DoFs. The example is nice, since results for a similar
2646test case are shown in multiple publications. As before, we slightly adapted the
2647test case to be able to use the same mesh for all simulations. The pressure field
2648at @f$t=0.3@f$ is depicted below.
2650@image html https://www.dealii.org/images/steps/developer/step_89_inhomogenous_test_case_pressure.png "" width=30%
2652As expected, we can easily see that the wave length in the right domain is roughly
2653three times times the wave length in the left domain. Hence, the wave can be
2654resolved with a coarser discretization.
2656Using the same element size in the whole domain, we can compute a reference result.
2657The displayed reference result is obtained by choosing the same subdivision level
2658for both sub-domains, i.e. @c subdiv_right = 11. In this particular example the
2659reference result uses @f$92,928@f$ DoFs, while the non-matching result uses @f$52,608@f$ DoFs.
2660The pressure result is plotted along the horizontal line that spans from (0,0.5) to (1,0.5).
2662@image html https://www.dealii.org/images/steps/developer/step_89_inhomogenous_test_case_conforming_vs_nonmatching.svg "" width=60%
2664The results computed using the non-matching discretization are clearly in good agreement with
2665the reference result.
2667<a name="step_89-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
2670All the implementations are done with overlapping triangulations in mind. In particular the
2671intersections in the mortaring case are constructed such that they are computed correctly
2672for overlapping triangulations. For this the intersection requests are of dimension @f$dim-1@f$.
2673The cells which are intersected with the intersection requests are of dimension @f$dim@f$. For the
2674simple case of non-conforming interfaces it would be sufficient to compute the intersections
2675between @f$dim-1@f$ and @f$dim-1@f$ entities. Furthermore, the lambda could be adapted, such that cells are
2676only marked if they are connected to a certain boundary ID (in this case, e.g., 99) instead of
2677marking every cell that is <i>not</i> connected to the opposite boundary ID (in this case, e.g., 98).
2678Marking fewer cells can reduce the setup costs significantly.
2680Note that the use of inhomogeneous materials in this procedure is questionable, since it is not clear which
2681material is present in the overlapping part of the mesh.
2684<a name="step_89-PlainProg"></a>
2685<h1> The plain program</h1>
2686@include "step-89.cc"
T read_cell_data(const AlignedVector< T > &array) const
std::conditional_t< rank_==1, Number, Tensor< rank_ - 1, dim, Number > > value_type
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
std::vector< index_type > data
@ component_is_part_of_vector
void write_vtu(const std::vector< Patch< dim, spacedim > > &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation > > &nonscalar_data_ranges, const VtkFlags &flags, std::ostream &out)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int one
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)