1006 * For the iterative solution procedure of the closest-
point projection,
1007 * the maximum number of iterations and the tolerance
for the maximum
1008 * absolute acceptable change in the distance in one iteration are set.
1011 * pcout <<
" - determine closest point iteratively" << std::endl;
1012 *
constexpr int max_iter = 30;
1013 *
constexpr double tol_distance = 1
e-6;
1017 * Now, we are ready to perform the algorithm by setting an
initial guess
1018 *
for the projection points simply corresponding to the collected support
1019 * points. We collect the global indices of the support points and the
1020 * total number of points that need to be processed and
do not
1021 * fulfill the required tolerance. Those will be gradually reduced
1022 * upon the iterative process.
1025 * std::vector<Point<dim>> closest_points = support_points;
1027 * std::vector<unsigned int> unmatched_points_idx(closest_points.size());
1028 * std::iota(unmatched_points_idx.begin(), unmatched_points_idx.end(), 0);
1030 *
int n_unmatched_points =
1036 * start the
loop for the fix-
point iteration. We update the list of points
1037 * that still need to be processed and subsequently pass
this information
1039 * illustration, we
export the coordinates of the points to be updated
for
1040 * each iteration to a CSV file. Next, we can evaluate the
signed distance
1041 * function and the
gradient at those points to update the current solution
1042 *
for the closest points. We perform the update only
if the
signed
1043 * distance of the closest
point is not already within the tolerance
1044 * and
register those points that still need to be processed.
1049 *
for (
int it = 0; it < max_iter && n_unmatched_points > 0; ++it)
1051 * pcout <<
" - iteration " << it <<
": " << n_unmatched_points;
1053 * std::vector<Point<dim>> unmatched_points(unmatched_points_idx.size());
1054 *
for (
unsigned int i = 0; i < unmatched_points_idx.size(); ++i)
1055 * unmatched_points[i] = closest_points[unmatched_points_idx[i]];
1057 *
const auto all_unmatched_points =
1058 * Utilities::MPI::reduce<std::vector<Point<dim>>>(
1059 * unmatched_points, MPI_COMM_WORLD, [](
const auto &a,
const auto &
b) {
1061 * result.insert(result.end(),
b.begin(),
b.end());
1067 * std::ofstream file(
"example_2_" + std::to_string(it) +
".csv");
1068 *
for (
const auto &p : all_unmatched_points)
1069 * file << p <<
std::endl;
1073 * rpe.reinit(unmatched_points, tria, mapping);
1076 * ExcMessage(
"Processed point is outside domain."));
1078 *
const auto eval_values =
1079 * VectorTools::point_values<1>(rpe, dof_handler, signed_distance);
1081 *
const auto eval_gradient =
1082 * VectorTools::point_gradients<1>(rpe, dof_handler, signed_distance);
1084 * std::vector<unsigned int> unmatched_points_idx_next;
1086 *
for (
unsigned int i = 0; i < unmatched_points_idx.size(); ++i)
1087 *
if (
std::abs(eval_values[i]) > tol_distance)
1089 * closest_points[unmatched_points_idx[i]] -=
1090 * eval_values[i] * eval_gradient[i];
1092 * unmatched_points_idx_next.emplace_back(unmatched_points_idx[i]);
1095 * unmatched_points_idx.swap(unmatched_points_idx_next);
1097 * n_unmatched_points =
1100 * pcout <<
" -> " << n_unmatched_points << std::endl;
1105 * We print a warning message
if we exceed the maximum number of allowed
1106 * iterations and
if there are still projection points with a distance
1107 *
value exceeding the tolerance.
1110 *
if (n_unmatched_points > 0)
1111 * pcout <<
"WARNING: The tolerance of " << n_unmatched_points
1112 * <<
" points is not yet attained." << std::endl;
1116 * As a result, we obtain a list of support points and corresponding
1117 * closest points at the zero-isosurface
level set. This information
1118 * can be used to update the
signed distance function, i.e., the
1119 * reinitialization the
values of the
level-set function to maintain
1120 * the
signed distance
property @cite henri2022geometrical.
1123 * pcout <<
" - determine distance in narrow band" << std::endl;
1125 * solution_distance.
reinit(solution);
1127 *
for (
unsigned int i = 0; i < closest_points.size(); ++i)
1128 * solution_distance[support_points_idx[i]] =
1129 * support_points[i].distance(closest_points[i]);
1133 * In addition, we use the information of the closest
point to
1135 * set isosurface, to the support points within the narrow band.
1136 * This might be helpful to improve accuracy,
e.g.,
for
1137 * diffuse
interface fluxes where certain quantities are only
1138 * accurately determined at the interface (
e.g. curvature
1139 *
for surface tension @cite coquerelle2016fourth).
1142 * pcout <<
" - perform extrapolation in narrow band" << std::endl;
1143 * rpe.reinit(closest_points, tria, mapping);
1144 * solution.update_ghost_values();
1145 *
const auto vals = VectorTools::point_values<1>(rpe, dof_handler, solution);
1148 * solution_extrapolated.
reinit(solution);
1150 *
for (
unsigned int i = 0; i < closest_points.size(); ++i)
1151 * solution_extrapolated[support_points_idx[i]] = vals[i];
1155 * Finally, we output the results to a VTU file.
1158 * pcout <<
" - output results" << std::endl;
1160 * data_out.
add_data_vector(dof_handler, signed_distance,
"signed_distance");
1161 * data_out.add_data_vector(dof_handler, solution,
"solution");
1162 * data_out.add_data_vector(dof_handler,
1163 * solution_distance,
1164 *
"solution_distance");
1165 * data_out.add_data_vector(dof_handler,
1166 * solution_extrapolated,
1167 *
"solution_extrapolated");
1168 * data_out.build_patches(mapping);
1169 * data_out.write_vtu_in_parallel(
"example_2.vtu", MPI_COMM_WORLD);
1171 * pcout << std::endl;
1177 * <a name=
"step_87-Miniexample3Sharpinterfacemethodontheexampleofsurfacetensionforfronttracking"></a>
1178 * <h3>Mini example 3: Sharp
interface method on the example of surface tension for front tracking</h3>
1182 * The
final mini example presents a basic implementation of
1183 * front tracking @cite peskin1977numerical, @cite unverdi1992front
1184 * of a surface mesh @f$\mathbb{T}_\Gamma@f$ immersed
1185 * in a Eulerian background fluid mesh @f$\mathbb{T}_\Omega@f$.
1189 * We assume that the immersed surface is transported according to a
1190 * prescribed velocity field from the background mesh. Subsequently,
1191 * we perform a sharp computation of the surface-tension force:
1193 * (\boldsymbol v_i (\boldsymbol{x}), \boldsymbol F_S
1194 * (\boldsymbol{x}))_{\Omega}
1196 * \left( \boldsymbol{v}_i (\boldsymbol{x}), \sigma (\boldsymbol{x}) \kappa
1197 * (\boldsymbol{x}) \boldsymbol{n} (\boldsymbol{x}) \right)_\Gamma \approx
1198 * \sum_{q\in\mathbb{T}_\Gamma} \boldsymbol{v}_i^T (\boldsymbol{x}_q)
1199 * \sigma (\boldsymbol{x}_q) \kappa (\boldsymbol{x}_q) \boldsymbol{n}
1200 * (\boldsymbol{x}_q) |J(\boldsymbol{x}_q)|
w(\boldsymbol{x}_q) \quad \forall
1201 * i\in\mathbb{T}_\Omega
1204 * We decompose
this operation into two steps. In the
first step, we evaluate
1205 * the force contributions @f$\sigma (\boldsymbol{x}_q) \kappa
1206 * (\boldsymbol{x}_q) \boldsymbol{n}
1207 * (\boldsymbol{x}_q)@f$ at the quadrature points defined on the immersed mesh
1208 * and multiply them with the mapped quadrature weight @f$|J(\boldsymbol{x}_q)|
1211 * \boldsymbol{
F}_S (\boldsymbol{x}_q) \gets \sigma (\boldsymbol{x}_q) \kappa
1212 * (\boldsymbol{x}_q) \boldsymbol{n} (\boldsymbol{x}_q) |J(\boldsymbol{x}_q)|
1213 * w_q \quad \forall q\in\mathbb{T}_\Gamma.
1215 * In the
second step, we compute the discretized weak form by multiplying
1216 * with test
functions on the background mesh:
1218 * (\boldsymbol v_i (\boldsymbol{x}), \boldsymbol F_S
1219 * (\boldsymbol{x}))_{\Omega} \gets \sum_q \boldsymbol{v}_i^T
1220 * (\boldsymbol{x}_q) \boldsymbol{
F}_S
1221 * (\boldsymbol{x}_q)
1222 * \quad \forall i\in\mathbb{T}_\Omega
1225 * Obviously, we need to communicate between the two steps. The
second step
1226 * can be handled completely by
Utilities::
MPI::RemotePointEvaluation, which
1227 * provides the function
1228 *
Utilities::
MPI::RemotePointEvaluation::process_and_evaluate() for this
1233 * We start with setting the parameters consisting of the polynomial degree of
1234 * the shape
functions, the dimension of the background mesh, the time-step
1235 * size to be considered for transporting the surface mesh and the number of
1244 * constexpr unsigned
int degree = 3;
1245 *
constexpr unsigned int dim = 2;
1246 *
const double dt = 0.01;
1247 *
const unsigned int n_time_steps = 200;
1251 * This program is intended to be executed in 2D or 3D.
1254 *
static_assert(dim == 2 || dim == 3,
"Only implemented for 2D or 3D.");
1260 * pcout <<
"Running: example 3" << std::endl;
1264 * Next, we create the standard objects necessary
for the finite element
1265 * representation of the background mesh
1268 * pcout <<
" - creating background mesh" << std::endl;
1269 * DistributedTriangulation<dim> tria_background(MPI_COMM_WORLD);
1271 * tria_background.refine_global(5);
1276 * dof_handler_background.distribute_dofs(fe_background);
1280 * and, similarly,
for the immersed surface mesh.
1281 * We use a sphere with radius @f$r=0.75@f$ which is
1282 * placed in the
center of the top half of the cubic background domain.
1285 * pcout <<
" - creating immersed mesh" << std::endl;
1287 *
Point<dim>(0.5, 0.75, 0.5));
1288 *
const double radius = 0.15;
1290 * DistributedTriangulation<dim - 1, dim> tria_immersed(MPI_COMM_WORLD);
1292 * tria_immersed.refine_global(4);
1296 * Two different mappings are considered
for the immersed
1297 * surface mesh: one
valid for the
initial configuration and one
1298 * that is updated in every time step according to the nodal
1299 * displacements. Two
types of finite elements are used to
1303 *
const MappingQ<dim - 1, dim> mapping_immersed_base(3);
1305 * mapping_immersed.initialize(mapping_immersed_base, tria_immersed);
1306 *
const QGauss<dim - 1> quadrature_immersed(degree + 1);
1308 *
const FE_Q<dim - 1, dim> fe_scalar_immersed(degree);
1309 *
const FESystem<dim - 1, dim> fe_immersed(fe_scalar_immersed, dim);
1310 *
DoFHandler<dim - 1, dim> dof_handler_immersed(tria_immersed);
1311 * dof_handler_immersed.distribute_dofs(fe_immersed);
1315 * We renumber the DoFs related to the vector-valued problem to
1316 * simplify access to the individual components.
1323 * We fill a DoF vector on the background mesh with an analytical
1324 * velocity field considering the Rayleigh-Kothe vortex flow and
1325 * initialize a DoF vector
for the weak form of the surface-tension force.
1329 * velocity.
reinit(dof_handler_background.locally_owned_dofs(),
1331 * dof_handler_background),
1336 * dof_handler_background.locally_owned_dofs(),
1342 * Next, we collect the real positions @f$\boldsymbol{x}_q@f$ of the quadrature
1343 * points of the surface mesh in a vector.
1347 * collect_support_points(mapping_immersed,
1348 * dof_handler_immersed,
1349 * immersed_support_points);
1354 * the time
loop. For any other step than the
initial one, we
first move the
1355 * support points of the surface mesh according to the fluid velocity of the
1356 * background mesh. Thereto, we
first update the time of the velocity
1357 * function. Then, we update the
internal data structures of the
1359 * points of the immersed mesh. We
throw an exception
if one of the points
1360 * cannot be found within the domain of the background mesh. Next, we
1361 * evaluate the velocity at the surface-mesh support points and compute the
1362 * resulting update of the coordinates. Finally, we update the mapping of
1363 * the immersed surface mesh to the current position.
1367 *
double time = 0.0;
1368 *
for (
unsigned int it = 0; it <= n_time_steps; ++it, time += dt)
1370 * pcout <<
"time: " << time << std::endl;
1373 * pcout <<
" - move support points (immersed mesh)" << std::endl;
1374 * vortex.set_time(time);
1376 * dof_handler_background,
1379 * rpe.reinit(convert<dim>(immersed_support_points),
1381 * mapping_background);
1385 *
"Immersed domain leaves background domain!"));
1387 * velocity.update_ghost_values();
1388 *
const auto immersed_velocity =
1389 * VectorTools::point_values<dim>(rpe,
1390 * dof_handler_background,
1392 * velocity.zero_out_ghost_values();
1394 *
for (
unsigned int i = 0, c = 0;
1395 * i < immersed_support_points.locally_owned_size() / dim;
1397 *
for (
unsigned int d = 0;
d < dim; ++
d, ++c)
1398 * immersed_support_points.local_element(c) +=
1399 * immersed_velocity[i][
d] * dt;
1401 * mapping_immersed.initialize(mapping_immersed_base,
1402 * dof_handler_immersed,
1403 * immersed_support_points,
1409 * Next, we
loop over all locally owned cells of the immersed mesh and
1410 * its quadrature points to compute the
value for the local surface
1411 * tension force contribution @f$\boldsymbol{
F}_S(\boldsymbol{x}_q)@f$. We
1412 * store the real coordinates of the quadrature points and the
1413 * corresponding force contributions in two individual vectors. For
1414 * computation of the latter, the normal vector
1415 * @f$\boldsymbol{n}(\boldsymbol{x}_q)@f$ can be directly extracted from the
1416 * surface mesh via
FEValues and,
for the curvature, we use the
1417 * following approximation:
1419 * \kappa(\boldsymbol{x}_q)
1421 * \nabla \cdot \boldsymbol{n}(\boldsymbol{x}_q)
1423 * \text{tr}\left({\nabla \boldsymbol{n}(\boldsymbol{x}_q)}\right)
1425 * \text{tr}\left({\nabla \sum_i \boldsymbol{N}_i (\boldsymbol{x}_q)
1426 * \boldsymbol n_i}\right)
1428 * \sum_i\text{tr}\left({\nabla \boldsymbol{N}_i (\boldsymbol{x}_q)
1429 * \boldsymbol n_i}\right)
1430 * \;\text{with}\; i\in[0,n_{\text{dofs\_per\_cell}}),
1432 * which we can apply since the immersed mesh is consistently
1433 * orientated. The surface tension coefficient is set to 1
for the
1434 * sake of demonstration.
1437 * pcout <<
" - compute to be tested values (immersed mesh)" << std::endl;
1440 * std::vector<Point<dim>> integration_points;
1441 * std::vector<value_type> integration_values;
1443 *
FEValues<dim - 1, dim> fe_values(mapping_immersed,
1445 * quadrature_immersed,
1450 *
FEValues<dim - 1, dim> fe_values_co(
1452 * fe_scalar_immersed,
1453 * fe_scalar_immersed.get_unit_support_points(),
1456 * std::vector<unsigned int> component_to_system_index(
1457 * fe_immersed.n_dofs_per_cell());
1459 *
for (
unsigned int i = 0, c = 0;
1460 * i < fe_scalar_immersed.n_dofs_per_cell();
1462 *
for (
unsigned int d = 0;
d < dim; ++
d, ++c)
1463 * component_to_system_index[c] =
1464 * fe_immersed.component_to_system_index(d, i);
1466 *
for (
const auto &cell : tria_immersed.active_cell_iterators() |
1469 * fe_values.
reinit(cell);
1470 * fe_values_co.reinit(cell);
1472 *
for (
const auto &q : fe_values.quadrature_point_indices())
1474 * integration_points.emplace_back(fe_values.quadrature_point(q));
1476 *
const auto sigma = 1.0;
1478 *
const auto normal = fe_values.normal_vector(q);
1479 *
double curvature = 0;
1480 *
for (
unsigned int i = 0, c = 0;
1481 * i < fe_scalar_immersed.n_dofs_per_cell();
1483 *
for (
unsigned int d = 0;
d < dim; ++
d, ++c)
1484 * curvature += fe_values.shape_grad_component(
1485 * component_to_system_index[c], q, d)[
d] *
1486 * fe_values_co.normal_vector(i)[
d];
1488 *
const auto FxJxW =
1489 * sigma * curvature * normal * fe_values.JxW(q);
1491 * integration_values.emplace_back(FxJxW);
1497 * Before we evaluate the weak form of the surface-tension force, the
1499 * set up from the quadrature points of the immersed mesh, determining
1500 * the surrounding cells on the background mesh.
1503 * pcout <<
" - test values (background mesh)" << std::endl;
1505 * rpe.
reinit(integration_points, tria_background, mapping_background);
1509 * In preparation
for utilizing
1512 * multiplication with the test function, we set up a callback function
1513 * that contains the operation on the
intersected cells of the
1514 * background mesh. Within
this function, we initialize a
1516 * arbitrary points within a cell. We
loop over the cells that surround
1517 * quadrature points of the immersed mesh -- provided by the callback
1518 * function. From the provided
CellData object, we retrieve the unit
1519 * points, i.e., the quadrature points of the immersed mesh that lie
1520 * within the current cell and a pointer to the stored
values on the
1521 * current cell (local surface-tension force)
for convenience. We
1523 * according to the unit points. Next, we
loop over the quadrature
1524 * points attributed to the cell and submit the local surface-tension
1527 * multiplied by the
values of the test function and a summation over
1528 * all given points is performed. Subsequently, the contributions are
1529 * assembled into the global vector containing the weak form of the
1530 * surface-tension force.
1533 *
const auto integration_function = [&](
const auto &
values,
1534 *
const auto &cell_data) {
1539 * std::vector<double> local_values;
1540 * std::vector<types::global_dof_index> local_dof_indices;
1542 *
for (
const auto cell : cell_data.cell_indices())
1544 * const auto cell_dofs =
1545 * cell_data.get_active_cell_iterator(cell)
1546 * ->as_dof_handler_iterator(dof_handler_background);
1548 *
const auto unit_points = cell_data.get_unit_points(cell);
1549 *
const auto FxJxW = cell_data.get_data_view(cell, values);
1551 * phi_force.reinit(cell_dofs, unit_points);
1553 *
for (
const auto q : phi_force.quadrature_point_indices())
1554 * phi_force.submit_value(FxJxW[q], q);
1556 * local_values.resize(cell_dofs->get_fe().n_dofs_per_cell());
1559 * local_dof_indices.resize(cell_dofs->get_fe().n_dofs_per_cell());
1560 * cell_dofs->get_dof_indices(local_dof_indices);
1562 * local_values, local_dof_indices, force_vector);
1568 * The callback function is passed together with the vector holding the
1569 * surface-tension force contribution at each quadrature
point of the
1572 * missing step is to
compress the distributed force vector.
1575 * rpe.process_and_evaluate<value_type>(integration_values,
1576 * integration_function);
1581 * After every tenth step or at the beginning/
end of the time
loop, we
1582 * output the force vector and the velocity of the background mesh to
1583 * a VTU file. In addition, we also
export the geometry of the
1584 * (deformed) immersed surface mesh to a separate VTU file.
1587 * if (it % 10 == 0 || it == n_time_steps)
1591 * vector_component_interpretation(
1593 * pcout <<
" - write data (background mesh)" << std::endl;
1597 * data_out_background.set_flags(flags_background);
1598 * data_out_background.add_data_vector(
1599 * dof_handler_background,
1601 * std::vector<std::string>(dim,
"force"),
1602 * vector_component_interpretation);
1603 * data_out_background.add_data_vector(
1604 * dof_handler_background,
1606 * std::vector<std::string>(dim,
"velocity"),
1607 * vector_component_interpretation);
1608 * data_out_background.build_patches(mapping_background, 3);
1609 * data_out_background.write_vtu_in_parallel(
"example_3_background_" +
1610 * std::to_string(it) +
1614 * pcout <<
" - write mesh (immersed mesh)" << std::endl;
1615 *
DataOut<dim - 1, dim> data_out_immersed;
1617 * data_out_immersed.build_patches(mapping_immersed, 3);
1618 * data_out_immersed.write_vtu_in_parallel(
"example_3_immersed_" +
1619 * std::to_string(it) +
1623 * pcout << std::endl;
1630 * <a name=
"step_87-Utilityfunctionsdefinition"></a>
1631 * <h3>Utility
functions (definition)</h3>
1634 *
template <
int spacedim>
1635 * std::vector<Point<spacedim>>
1638 *
const unsigned int n_subdivisions)
1640 *
Assert(n_subdivisions >= 1, ExcInternalError());
1642 * std::vector<Point<spacedim>> points;
1643 * points.reserve(n_subdivisions + 1);
1645 * points.emplace_back(p0);
1646 *
for (
unsigned int i = 1; i < n_subdivisions; ++i)
1647 * points.emplace_back(p0 + (p1 - p0) *
static_cast<double>(i) /
1648 *
static_cast<double>(n_subdivisions));
1649 * points.emplace_back(p1);
1654 *
template <
int spacedim,
typename T0,
typename T1>
1655 *
void print_along_line(
const std::string &file_name,
1657 *
const std::vector<T0> &values_0,
1658 *
const std::vector<T1> &values_1)
1661 * (values_1.size() == points.size() || values_1.empty()),
1662 * ExcMessage(
"The provided vectors must have the same length."));
1664 * std::ofstream file(file_name);
1666 *
for (
unsigned int i = 0; i < points.size(); ++i)
1668 * file << std::fixed << std::right << std::setw(5) << std::setprecision(3)
1669 * << points[0].distance(points[i]);
1671 *
for (
unsigned int d = 0;
d < spacedim; ++
d)
1672 * file << std::fixed << std::right << std::setw(10)
1673 * << std::setprecision(3) << points[i][
d];
1675 * file << std::fixed << std::right << std::setw(10)
1676 * << std::setprecision(3) << values_0[i];
1678 *
if (!values_1.empty())
1679 * file << std::fixed << std::right << std::setw(10)
1680 * << std::setprecision(3) << values_1[i];
1681 * file << std::endl;
1685 *
template <
int dim,
int spacedim>
1686 *
void collect_support_points(
1691 * support_points.reinit(dof_handler.locally_owned_dofs(),
1693 * dof_handler.get_mpi_communicator());
1695 *
const auto &fe = dof_handler.get_fe();
1700 * fe.base_element(0).get_unit_support_points(),
1703 * std::vector<types::global_dof_index> local_dof_indices(
1704 * fe.n_dofs_per_cell());
1706 * std::vector<unsigned int> component_to_system_index(
1707 * fe_values.n_quadrature_points * spacedim);
1709 *
for (
unsigned int q = 0, c = 0; q < fe_values.n_quadrature_points; ++q)
1710 *
for (
unsigned int d = 0;
d < spacedim; ++
d, ++c)
1711 * component_to_system_index[c] = fe.component_to_system_index(d, q);
1713 *
for (
const auto &cell : dof_handler.active_cell_iterators() |
1716 * fe_values.
reinit(cell);
1717 * cell->get_dof_indices(local_dof_indices);
1719 *
for (
unsigned int q = 0, c = 0; q < fe_values.n_quadrature_points; ++q)
1720 *
for (
unsigned int d = 0;
d < spacedim; ++
d, ++c)
1721 * support_points[local_dof_indices[component_to_system_index[c]]] =
1722 * fe_values.quadrature_point(q)[
d];
1726 *
template <
int dim,
int spacedim,
typename T>
1727 * std::tuple<std::vector<Point<spacedim>>, std::vector<types::global_dof_index>>
1728 * collect_support_points_with_narrow_band(
1733 *
const double narrow_band_threshold)
1736 * ExcMessage(
"The narrow band threshold"
1737 *
" must be larger than or equal to 0."));
1738 *
const auto &tria = dof_handler_signed_distance.get_triangulation();
1741 * .get_unit_support_points());
1744 * dof_handler_signed_distance.get_fe(),
1749 * dof_handler_support_points.get_fe(),
1753 * std::vector<T> temp_distance(quad.size());
1754 * std::vector<types::global_dof_index> local_dof_indices(
1755 * dof_handler_support_points.get_fe().n_dofs_per_cell());
1757 * std::vector<Point<dim>> support_points;
1758 * std::vector<types::global_dof_index> support_points_idx;
1760 *
const bool has_ghost_elements = signed_distance.has_ghost_elements();
1762 *
const auto &locally_owned_dofs_req =
1763 * dof_handler_support_points.locally_owned_dofs();
1764 * std::vector<bool> flags(locally_owned_dofs_req.n_elements(),
false);
1766 *
if (has_ghost_elements ==
false)
1767 * signed_distance.update_ghost_values();
1769 *
for (
const auto &cell :
1772 * const auto cell_distance =
1773 * cell->as_dof_handler_iterator(dof_handler_signed_distance);
1774 * distance_values.reinit(cell_distance);
1775 * distance_values.get_function_values(signed_distance, temp_distance);
1777 *
const auto cell_req =
1778 * cell->as_dof_handler_iterator(dof_handler_support_points);
1779 * req_values.reinit(cell_req);
1780 * cell_req->get_dof_indices(local_dof_indices);
1782 *
for (
const auto q : req_values.quadrature_point_indices())
1783 * if (
std::
abs(temp_distance[q]) < narrow_band_threshold)
1785 * const auto idx = local_dof_indices[q];
1787 *
if (locally_owned_dofs_req.is_element(idx) ==
false ||
1788 * flags[locally_owned_dofs_req.index_within_set(idx)])
1791 * flags[locally_owned_dofs_req.index_within_set(idx)] =
true;
1793 * support_points_idx.emplace_back(idx);
1794 * support_points.emplace_back(req_values.quadrature_point(q));
1798 *
if (has_ghost_elements ==
false)
1799 * signed_distance.zero_out_ghost_values();
1801 *
return {support_points, support_points_idx};
1804 *
template <
int spacedim>
1805 * std::vector<Point<spacedim>> convert(
1808 *
const unsigned int n_points =
1809 * support_points_unrolled.locally_owned_size() / spacedim;
1811 * std::vector<Point<spacedim>> points(n_points);
1813 *
for (
unsigned int i = 0, c = 0; i < n_points; ++i)
1814 *
for (
unsigned int d = 0;
d < spacedim; ++
d, ++c)
1815 * points[i][d] = support_points_unrolled.local_element(c);
1825 * <a name=
"step_87-Driver"></a>
1830 * Finally, the driver of the program executes the four mini examples.
1833 *
int main(
int argc,
char **argv)
1835 *
using namespace dealii;
1837 * std::cout.precision(5);
1840 * Step87::example_0();
1842 * Step87::example_1();
1843 * Step87::example_2();
1844 * Step87::example_3();
1847<a name=
"step_87-Results"></a><h1>Results</h1>
1850<a name=
"step_87-Miniexample0"></a><h3>Mini example 0</h3>
1853We present a part of the terminal output. It shows,
for each point, the
1854determined cell and reference position. Also, one can see that
1860 - Found point with real coordinates: 0 0.5
1861 - in cell with
vertices: (0 0.4) (0.2 0.4) (0 0.6) (0.2 0.6)
1862 - with coordinates on the unit cell: (0 0.5)
1868 - Found point with real coordinates: 0.05 0.5
1869 - in cell with
vertices: (0 0.4) (0.2 0.4) (0 0.6) (0.2 0.6)
1870 - with coordinates on the unit cell: (0.25 0.5)
1878 - Found point with real coordinates: 1 0.5
1879 - in cell with
vertices: (0.8 0.4) (1 0.4) (0.8 0.6) (1 0.6)
1880 - with coordinates on the unit cell: (1 0.5)
1889The CSV output is as follows and contains, in the
1890first column, the distances with respect to the
first point,
1891the
second and the third column represent the coordinates
1892of the points and the fourth column the evaluated solution
1893values at those points.
18960.000 0.000 0.500 0.250
18970.050 0.050 0.500 0.200
18980.100 0.100 0.500 0.150
18990.150 0.150 0.500 0.100
19000.200 0.200 0.500 0.050
19010.250 0.250 0.500 0.000
19020.300 0.300 0.500 -0.050
19030.350 0.350 0.500 -0.100
19040.400 0.400 0.500 -0.149
19050.450 0.450 0.500 -0.200
19060.500 0.500 0.500 -0.222
19070.550 0.550 0.500 -0.200
19080.600 0.600 0.500 -0.149
19090.650 0.650 0.500 -0.100
19100.700 0.700 0.500 -0.050
19110.750 0.750 0.500 0.000
19120.800 0.800 0.500 0.050
19130.850 0.850 0.500 0.100
19140.900 0.900 0.500 0.150
19150.950 0.950 0.500 0.200
19161.000 1.000 0.500 0.250
1919<a name="step_87-Miniexample1"></a><h3>Mini example 1</h3>
1922We show the terminal output.
1929The CSV output is as follows and identical to the results
1930of the serial case presented in mini example 0.
1931The fifth column shows the
1932user quantity evaluated additionally in this mini example.
19350.000 0.000 0.500 0.250 0.000
19360.050 0.050 0.500 0.200 0.050
19370.100 0.100 0.500 0.150 0.100
19380.150 0.150 0.500 0.100 0.150
19390.200 0.200 0.500 0.050 0.200
19400.250 0.250 0.500 0.000 0.250
19410.300 0.300 0.500 -0.050 0.300
19420.350 0.350 0.500 -0.100 0.350
19430.400 0.400 0.500 -0.149 0.400
19440.450 0.450 0.500 -0.200 0.450
19450.500 0.500 0.500 -0.222 0.500
19460.550 0.550 0.500 -0.200 0.550
19470.600 0.600 0.500 -0.149 0.600
19480.650 0.650 0.500 -0.100 0.650
19490.700 0.700 0.500 -0.050 0.700
19500.750 0.750 0.500 0.000 0.750
19510.800 0.800 0.500 0.050 0.800
19520.850 0.850 0.500 0.100 0.850
19530.900 0.900 0.500 0.150 0.900
19540.950 0.950 0.500 0.200 0.950
19551.000 1.000 0.500 0.250 1.000
1959<a name="step_87-Miniexample2"></a><h3>Mini example 2</h3>
1962We show the terminal output.
1966 - determine narrow band
1967 - determine closest point iteratively
1968 - iteration 0: 7076 -> 7076
1969 - iteration 1: 7076 -> 104
1970 - iteration 2: 104 -> 0
1971 - determine distance in narrow band
1972 - perform extrapolation in narrow band
1976The following three plots, representing the performed iterations of the
1977closest-point projection, show the current position of the closest
1978points exceeding the required tolerance of the discrete interface
1979of the circle and still need to
1981It can be seen that the
numbers of points to be processed decrease
1982from iteration to iteration.
1983<table align="
center" class="doxtable">
1997The output visualized in Paraview looks like the following: On the
1998left, the original distance function is shown as the light gray surface.
1999In addition, the contour values refer to the distance values determined
2000from calculation of the distance to the closest points at the interface
2001in the narrow band. It can be seen that the two functions coincide.
2002Similarly, on the right, the original solution and the extrapolated
2003solution from the interface is shown.
2005<table align="
center" class="doxtable">
2016<a name="step_87-Miniexample3"></a><h3>Mini example 3</h3>
2019We show a shortened version of the terminal output.
2023 - creating background mesh
2024 - creating immersed mesh
2026 - compute to be tested values (immersed mesh)
2027 - test values (background mesh)
2028 - write data (background mesh)
2029 - write mesh (immersed mesh)
2032 - move support points (immersed mesh)
2033 - compute to be tested values (immersed mesh)
2034 - test values (background mesh)
2037 - move support points (immersed mesh)
2038 - compute to be tested values (immersed mesh)
2039 - test values (background mesh)
2044 - move support points (immersed mesh)
2045 - compute to be tested values (immersed mesh)
2046 - test values (background mesh)
2047 - write data (background mesh)
2048 - write mesh (immersed mesh)
2051The output visualized in Paraview looks like the following: The deformation of
2052the immersed mesh by the reversible vortex flow can be seen. Due to
2053discretization errors, the shape is not exactly circular at the end, illustrated
2054in the right figure. The sharp nature of the surface-tension force vector, shown
2055as vector plots, can be seen by its restriction to cells that are intersected by
2058<table align="
center" class="doxtable">
2072<a name="step_87-Possibilitiesforextension"></a><h3>Possibilities for extension</h3>
2075This program highlights some of the main capabilities
2076of the distributed evaluation routines in deal.II. However, there are many
2077related topics worth mentioning:
2078- Performing a distributed search is an expensive step. That is why we suggest
2079to provide hints to
Utilities::
MPI::RemotePointEvaluation and to reuse
2081instances in the case that the communication pattern has not changed.
2082Furthermore, there are instances where no search is needed and the points are
2083already sorted into the right cells. This is the case if the points are
2084generated on the cell
level (see @ref step_85 "step-85"; CutFEM) or the points are
2085automatically sorted into the correct (neighboring) cell (see @ref step_68 "step-68"; PIC with
2086Particles::ParticleHandler). Having said that, the
2087Particles::ParticleHandler::insert_global_particles() function uses
2088the described infrastructure to perform the initial sorting of particles into
2090- We concentrated on parallelization aspects in this tutorial. However, we would
2091like to point out the need for fast evaluation on cell
level.
2094\hat{u}(\hat{\boldsymbol{x}}) = \sum_i \hat{N}_i(\hat{\boldsymbol{x}}) \hat{u}_i
2096to derive fast implementations,
e.g.,
for tensor-product elements
2098\hat{u}(\hat{x}_0, \hat{x}_1, \hat{x}_2) =
2099\sum_k \hat{N}^{\text{1D}}_k(\hat{x}_2)
2100\sum_j \hat{N}^{\text{1D}}_j(\hat{x}_1)
2101\sum_i \hat{N}^{\text{1D}}_i(\hat{x}_0)
2104Since only 1D shape
functions are queried and are re-used in re-occurring terms,
2105this formulation is faster than without exploitation of the structure.
2107The
class DataOutResample allows to output results on a different mesh than
2108the computational mesh. This is useful if one wants to output the results
2109on a coarser mesh or one does not want to output 3D results but instead 2D
2111and restrict residuals between two
independent meshes. By passing a sequence
2113non-nested multigrid can be computed.
2115grids via surfaces (example: fluid-structure interaction, independently created
2116grids). The evaluation points can come from any side (pointwise interpolation)
2117or are defined on
intersected meshes (Nitsche-type mortaring
2118@cite heinz2022high). Concerning the creation of such
intersected meshes and the
2119corresponding evaluation points, see
2120GridTools::internal::distributed\_compute_intersection_locations().
2122preCICE @cite bungartz2016precice @cite chourdakis2021precice can be used. The
2123code-gallery program "Laplace equation coupled to an external simulation
2124program" shows how to use this library with deal.II.
2127<a name="step_87-PlainProg"></a>
2128<h1> The plain program</h1>
2129@include
"step-87.cc"
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
void attach_triangulation(const Triangulation< dim, spacedim > &)
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
void test_and_sum(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
void reinit(const size_type size, const bool omit_zeroing_entries=false)
Abstract base class for mapping classes.
void process_and_evaluate(const std::vector< T > &input, std::vector< T > &buffer, const std::function< void(const ArrayView< const T > &, const CellData &)> &evaluation_function, const bool sort_data=true) const
void reinit(const std::vector< Point< spacedim > > &points, const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping)
#define Assert(cond, exc)
#define AssertThrow(cond, exc)
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())
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
DataComponentInterpretation
@ component_is_part_of_vector
void support_point_wise(DoFHandler< dim, spacedim > &dof_handler)
void hyper_sphere(Triangulation< spacedim - 1, spacedim > &tria, const Point< spacedim > ¢er=Point< spacedim >(), const double radius=1.)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
@ valid
Iterator points to a valid object.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::string compress(const std::string &input)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
bool write_higher_order_cells