deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-87.h
Go to the documentation of this file.
1);
1004 *  
1005 * @endcode
1006 *
1008 * the maximum number of iterations and the tolerance for the maximum
1009 * absolute acceptable change in the distance in one iteration are set.
1010 *
1011 * @code
1012 *   pcout << " - determine closest point iteratively" << std::endl;
1013 *   constexpr int max_iter = 30;
1014 *   constexpr double tol_distance = 1e-6;
1015 *  
1016 * @endcode
1017 *
1020 * points. We collect the global indices of the support points and the
1021 * total number of points that need to be processed and do not
1023 * upon the iterative process.
1024 *
1025 * @code
1026 *   std::vector<Point<dim>> closest_points = support_points; // initial guess
1027 *  
1028 *   std::vector<unsigned int> unmatched_points_idx(closest_points.size());
1029 *   std::iota(unmatched_points_idx.begin(), unmatched_points_idx.end(), 0);
1030 *  
1031 *   int n_unmatched_points =
1033 *  
1034 * @endcode
1035 *
1036 * Now, we create a Utilities::MPI::RemotePointEvaluation cache object and
1037 * start the loop for the fix-point iteration. We update the list of points
1040 * illustration, we export the coordinates of the points to be updated for
1041 * each iteration to a CSV file. Next, we can evaluate the signed distance
1042 * function and the gradient at those points to update the current solution
1043 * for the closest points. We perform the update only if the signed
1044 * distance of the closest point is not already within the tolerance
1045 * and register those points that still need to be processed.
1046 *
1047 * @code
1049 *  
1050 *   for (int it = 0; it < max_iter && n_unmatched_points > 0; ++it)
1051 *   {
1052 *   pcout << " - iteration " << it << ": " << n_unmatched_points;
1053 *  
1054 *   std::vector<Point<dim>> unmatched_points(unmatched_points_idx.size());
1055 *   for (unsigned int i = 0; i < unmatched_points_idx.size(); ++i)
1057 *  
1058 *   const auto all_unmatched_points =
1059 *   Utilities::MPI::reduce<std::vector<Point<dim>>>(
1060 *   unmatched_points, MPI_COMM_WORLD, [](const auto &a, const auto &b) {
1061 *   auto result = a;
1062 *   result.insert(result.end(), b.begin(), b.end());
1063 *   return result;
1064 *   });
1065 *  
1067 *   {
1068 *   std::ofstream file("example_2_" + std::to_string(it) + ".csv");
1069 *   for (const auto &p : all_unmatched_points)
1070 *   file << p << std::endl;
1071 *   file.close();
1072 *   }
1073 *  
1074 *   rpe.reinit(unmatched_points, tria, mapping);
1075 *  
1076 *   AssertThrow(rpe.all_points_found(),
1077 *   ExcMessage("Processed point is outside domain."));
1078 *  
1079 *   const auto eval_values =
1080 *   VectorTools::point_values<1>(rpe, dof_handler, signed_distance);
1081 *  
1082 *   const auto eval_gradient =
1083 *   VectorTools::point_gradients<1>(rpe, dof_handler, signed_distance);
1084 *  
1085 *   std::vector<unsigned int> unmatched_points_idx_next;
1086 *  
1087 *   for (unsigned int i = 0; i < unmatched_points_idx.size(); ++i)
1088 *   if (std::abs(eval_values[i]) > tol_distance)
1089 *   {
1091 *   eval_values[i] * eval_gradient[i];
1092 *  
1094 *   }
1095 *  
1097 *  
1100 *  
1101 *   pcout << " -> " << n_unmatched_points << std::endl;
1102 *   }
1103 *  
1104 * @endcode
1105 *
1106 * We print a warning message if we exceed the maximum number of allowed
1107 * iterations and if there are still projection points with a distance
1108 * value exceeding the tolerance.
1109 *
1110 * @code
1111 *   if (n_unmatched_points > 0)
1112 *   pcout << "WARNING: The tolerance of " << n_unmatched_points
1113 *   << " points is not yet attained." << std::endl;
1114 *  
1115 * @endcode
1116 *
1117 * As a result, we obtain a list of support points and corresponding
1119 * can be used to update the signed distance function, i.e., the
1121 * the signed distance property @cite henri2022geometrical.
1122 *
1123 * @code
1124 *   pcout << " - determine distance in narrow band" << std::endl;
1126 *   solution_distance.reinit(solution);
1127 *  
1128 *   for (unsigned int i = 0; i < closest_points.size(); ++i)
1130 *   support_points[i].distance(closest_points[i]);
1131 *  
1132 * @endcode
1133 *
1136 * set isosurface, to the support points within the narrow band.
1137 * This might be helpful to improve accuracy, e.g., for
1140 * for surface tension @cite coquerelle2016fourth).
1141 *
1142 * @code
1143 *   pcout << " - perform extrapolation in narrow band" << std::endl;
1144 *   rpe.reinit(closest_points, tria, mapping);
1145 *   solution.update_ghost_values();
1146 *   const auto vals = VectorTools::point_values<1>(rpe, dof_handler, solution);
1147 *  
1149 *   solution_extrapolated.reinit(solution);
1150 *  
1151 *   for (unsigned int i = 0; i < closest_points.size(); ++i)
1153 *  
1154 * @endcode
1155 *
1156 * Finally, we output the results to a VTU file.
1157 *
1158 * @code
1159 *   pcout << " - output results" << std::endl;
1160 *   DataOut<dim> data_out;
1161 *   data_out.add_data_vector(dof_handler, signed_distance, "signed_distance");
1162 *   data_out.add_data_vector(dof_handler, solution, "solution");
1163 *   data_out.add_data_vector(dof_handler,
1165 *   "solution_distance");
1166 *   data_out.add_data_vector(dof_handler,
1168 *   "solution_extrapolated");
1169 *   data_out.build_patches(mapping);
1170 *   data_out.write_vtu_in_parallel("example_2.vtu", MPI_COMM_WORLD);
1171 *  
1172 *   pcout << std::endl;
1173 *   }
1174 *  
1175 * @endcode
1176 *
1177 *
1178 * <a name="step_87-Miniexample3Sharpinterfacemethodontheexampleofsurfacetensionforfronttracking"></a>
1179 * <h3>Mini example 3: Sharp interface method on the example of surface tension for front tracking</h3>
1180 *
1181
1182 *
1185 * of a surface mesh @f$\mathbb{T}_\Gamma@f$ immersed
1186 * in a Eulerian background fluid mesh @f$\mathbb{T}_\Omega@f$.
1187 *
1188
1189 *
1193 * @f[
1195 * (\boldsymbol{x}))_{\Omega}
1196 * =
1198 * (\boldsymbol{x}) \boldsymbol{n} (\boldsymbol{x}) \right)_\Gamma \approx
1201 * (\boldsymbol{x}_q) |J(\boldsymbol{x}_q)| w(\boldsymbol{x}_q) \quad \forall
1202 * i\in\mathbb{T}_\Omega
1203 * .
1204 * @f]
1205 * We decompose this operation into two steps. In the first step, we evaluate
1207 * (\boldsymbol{x}_q) \boldsymbol{n}
1208 * (\boldsymbol{x}_q)@f$ at the quadrature points defined on the immersed mesh
1209 * and multiply them with the mapped quadrature weight @f$|J(\boldsymbol{x}_q)|
1210 * w_q@f$:
1211 * @f[
1213 * (\boldsymbol{x}_q) \boldsymbol{n} (\boldsymbol{x}_q) |J(\boldsymbol{x}_q)|
1214 * w_q \quad \forall q\in\mathbb{T}_\Gamma.
1215 * @f]
1217 * with test functions on the background mesh:
1218 * @f[
1220 * (\boldsymbol{x}))_{\Omega} \gets \sum_q \boldsymbol{v}_i^T
1222 * (\boldsymbol{x}_q)
1223 * \quad \forall i\in\mathbb{T}_\Omega
1224 * .
1225 * @f]
1227 * can be handled completely by Utilities::MPI::RemotePointEvaluation, which
1228 * provides the function
1229 * Utilities::MPI::RemotePointEvaluation::process_and_evaluate() for this
1230 * purpose.
1231 *
1232
1233 *
1234 * We start with setting the parameters consisting of the polynomial degree of
1235 * the shape functions, the dimension of the background mesh, the time-step
1236 * size to be considered for transporting the surface mesh and the number of
1237 * time steps.
1238 *
1239
1240 *
1241 *
1242 * @code
1243 *   void example_3()
1244 *   {
1245 *   constexpr unsigned int degree = 3;
1246 *   constexpr unsigned int dim = 2;
1247 *   const double dt = 0.01;
1248 *   const unsigned int n_time_steps = 200;
1249 *  
1250 * @endcode
1251 *
1252 * This program is intended to be executed in 2D or 3D.
1253 *
1254 * @code
1255 *   static_assert(dim == 2 || dim == 3, "Only implemented for 2D or 3D.");
1256 *  
1257 *   ConditionalOStream pcout(std::cout,
1259 *   0);
1260 *  
1261 *   pcout << "Running: example 3" << std::endl;
1262 *  
1263 * @endcode
1264 *
1265 * Next, we create the standard objects necessary for the finite element
1266 * representation of the background mesh
1267 *
1268 * @code
1269 *   pcout << " - creating background mesh" << std::endl;
1272 *   tria_background.refine_global(5);
1273 *  
1275 *   const FESystem<dim> fe_background(FE_Q<dim>(degree), dim);
1277 *   dof_handler_background.distribute_dofs(fe_background);
1278 *  
1279 * @endcode
1280 *
1281 * and, similarly, for the immersed surface mesh.
1282 * We use a sphere with radius @f$r=0.75@f$ which is
1283 * placed in the center of the top half of the cubic background domain.
1284 *
1285 * @code
1286 *   pcout << " - creating immersed mesh" << std::endl;
1287 *   const Point<dim> center((dim == 2) ? Point<dim>(0.5, 0.75) :
1288 *   Point<dim>(0.5, 0.75, 0.5));
1289 *   const double radius = 0.15;
1290 *  
1292 *   GridGenerator::hyper_sphere(tria_immersed, center, radius);
1293 *   tria_immersed.refine_global(4);
1294 *  
1295 * @endcode
1296 *
1298 * surface mesh: one valid for the initial configuration and one
1299 * that is updated in every time step according to the nodal
1300 * displacements. Two types of finite elements are used to
1302 *
1303 * @code
1304 *   const MappingQ<dim - 1, dim> mapping_immersed_base(3);
1305 *   MappingQCache<dim - 1, dim> mapping_immersed(3);
1307 *   const QGauss<dim - 1> quadrature_immersed(degree + 1);
1308 *  
1309 *   const FE_Q<dim - 1, dim> fe_scalar_immersed(degree);
1310 *   const FESystem<dim - 1, dim> fe_immersed(fe_scalar_immersed, dim);
1312 *   dof_handler_immersed.distribute_dofs(fe_immersed);
1313 *  
1314 * @endcode
1315 *
1316 * We renumber the DoFs related to the vector-valued problem to
1317 * simplify access to the individual components.
1318 *
1319 * @code
1321 *  
1322 * @endcode
1323 *
1324 * We fill a DoF vector on the background mesh with an analytical
1326 * initialize a DoF vector for the weak form of the surface-tension force.
1327 *
1328 * @code
1330 *   velocity.reinit(dof_handler_background.locally_owned_dofs(),
1333 *   MPI_COMM_WORLD);
1335 *  
1337 *   dof_handler_background.locally_owned_dofs(),
1339 *   MPI_COMM_WORLD);
1340 *  
1341 * @endcode
1342 *
1344 * points of the surface mesh in a vector.
1345 *
1346 * @code
1351 *  
1352 * @endcode
1353 *
1354 * We initialize a Utilities::MPI::RemotePointEvaluation object and start
1355 * the time loop. For any other step than the initial one, we first move the
1356 * support points of the surface mesh according to the fluid velocity of the
1357 * background mesh. Thereto, we first update the time of the velocity
1360 * points of the immersed mesh. We throw an exception if one of the points
1361 * cannot be found within the domain of the background mesh. Next, we
1362 * evaluate the velocity at the surface-mesh support points and compute the
1364 * the immersed surface mesh to the current position.
1365 *
1366 * @code
1368 *   double time = 0.0;
1369 *   for (unsigned int it = 0; it <= n_time_steps; ++it, time += dt)
1370 *   {
1371 *   pcout << "time: " << time << std::endl;
1372 *   if (it > 0)
1373 *   {
1374 *   pcout << " - move support points (immersed mesh)" << std::endl;
1375 *   vortex.set_time(time);
1378 *   vortex,
1379 *   velocity);
1383 *  
1384 *   AssertThrow(rpe.all_points_found(),
1385 *   ExcMessage(
1386 *   "Immersed domain leaves background domain!"));
1387 *  
1388 *   velocity.update_ghost_values();
1389 *   const auto immersed_velocity =
1390 *   VectorTools::point_values<dim>(rpe,
1392 *   velocity);
1393 *   velocity.zero_out_ghost_values();
1394 *  
1395 *   for (unsigned int i = 0, c = 0;
1396 *   i < immersed_support_points.locally_owned_size() / dim;
1397 *   ++i)
1398 *   for (unsigned int d = 0; d < dim; ++d, ++c)
1399 *   immersed_support_points.local_element(c) +=
1400 *   immersed_velocity[i][d] * dt;
1401 *  
1405 *   false);
1406 *   }
1407 *  
1408 * @endcode
1409 *
1411 * its quadrature points to compute the value for the local surface
1413 * store the real coordinates of the quadrature points and the
1415 * computation of the latter, the normal vector
1417 * surface mesh via FEValues and, for the curvature, we use the
1419 * @f[
1420 * \kappa(\boldsymbol{x}_q)
1421 * =
1423 * =
1424 * \text{tr}\left({\nabla \boldsymbol{n}(\boldsymbol{x}_q)}\right)
1425 * \approx
1427 * \boldsymbol n_i}\right)
1428 * =
1430 * \boldsymbol n_i}\right)
1431 * \;\text{with}\; i\in[0,n_{\text{dofs\_per\_cell}}),
1432 * @f]
1434 * orientated. The surface tension coefficient is set to 1 for the
1436 *
1437 * @code
1438 *   pcout << " - compute to be tested values (immersed mesh)" << std::endl;
1439 *   using value_type = Tensor<1, dim, double>;
1440 *  
1441 *   std::vector<Point<dim>> integration_points;
1442 *   std::vector<value_type> integration_values;
1443 *  
1444 *   FEValues<dim - 1, dim> fe_values(mapping_immersed,
1445 *   fe_immersed,
1450 *  
1451 *   FEValues<dim - 1, dim> fe_values_co(
1454 *   fe_scalar_immersed.get_unit_support_points(),
1456 *  
1457 *   std::vector<unsigned int> component_to_system_index(
1458 *   fe_immersed.n_dofs_per_cell());
1459 *  
1460 *   for (unsigned int i = 0, c = 0;
1461 *   i < fe_scalar_immersed.n_dofs_per_cell();
1462 *   ++i)
1463 *   for (unsigned int d = 0; d < dim; ++d, ++c)
1464 *   component_to_system_index[c] =
1465 *   fe_immersed.component_to_system_index(d, i);
1466 *  
1467 *   for (const auto &cell : tria_immersed.active_cell_iterators() |
1468 *   IteratorFilters::LocallyOwnedCell())
1469 *   {
1470 *   fe_values.reinit(cell);
1471 *   fe_values_co.reinit(cell);
1472 *  
1473 *   for (const auto &q : fe_values.quadrature_point_indices())
1474 *   {
1475 *   integration_points.emplace_back(fe_values.quadrature_point(q));
1476 *  
1477 *   const auto sigma = 1.0; // surface tension coefficient
1478 *  
1479 *   const auto normal = fe_values.normal_vector(q);
1480 *   double curvature = 0;
1481 *   for (unsigned int i = 0, c = 0;
1482 *   i < fe_scalar_immersed.n_dofs_per_cell();
1483 *   ++i)
1484 *   for (unsigned int d = 0; d < dim; ++d, ++c)
1485 *   curvature += fe_values.shape_grad_component(
1486 *   component_to_system_index[c], q, d)[d] *
1487 *   fe_values_co.normal_vector(i)[d];
1488 *  
1489 *   const auto FxJxW =
1490 *   sigma * curvature * normal * fe_values.JxW(q);
1491 *  
1492 *   integration_values.emplace_back(FxJxW);
1493 *   }
1494 *   }
1495 *  
1496 * @endcode
1497 *
1498 * Before we evaluate the weak form of the surface-tension force, the
1500 * set up from the quadrature points of the immersed mesh, determining
1501 * the surrounding cells on the background mesh.
1502 *
1503 * @code
1504 *   pcout << " - test values (background mesh)" << std::endl;
1505 *  
1507 *  
1508 * @endcode
1509 *
1512 * performs the
1513 * multiplication with the test function, we set up a callback function
1515 * background mesh. Within this function, we initialize a
1516 * FEPointEvaluation object that allows us to integrate values at
1517 * arbitrary points within a cell. We loop over the cells that surround
1518 * quadrature points of the immersed mesh -- provided by the callback
1519 * function. From the provided CellData object, we retrieve the unit
1520 * points, i.e., the quadrature points of the immersed mesh that lie
1521 * within the current cell and a pointer to the stored values on the
1522 * current cell (local surface-tension force) for convenience. We
1524 * according to the unit points. Next, we loop over the quadrature
1525 * points attributed to the cell and submit the local surface-tension
1526 * force to the FEPointEvaluation object. Via
1528 * multiplied by the values of the test function and a summation over
1530 * assembled into the global vector containing the weak form of the
1531 * surface-tension force.
1532 *
1533 * @code
1534 *   const auto integration_function = [&](const auto &values,
1535 *   const auto &cell_data) {
1537 *   fe_background,
1538 *   update_values);
1539 *  
1540 *   std::vector<double> local_values;
1541 *   std::vector<types::global_dof_index> local_dof_indices;
1542 *  
1543 *   for (const auto cell : cell_data.cell_indices())
1544 *   {
1546 *   cell_data.get_active_cell_iterator(cell)
1547 *   ->as_dof_handler_iterator(dof_handler_background);
1548 *  
1549 *   const auto unit_points = cell_data.get_unit_points(cell);
1550 *   const auto FxJxW = cell_data.get_data_view(cell, values);
1551 *  
1552 *   phi_force.reinit(cell_dofs, unit_points);
1553 *  
1554 *   for (const auto q : phi_force.quadrature_point_indices())
1555 *   phi_force.submit_value(FxJxW[q], q);
1556 *  
1557 *   local_values.resize(cell_dofs->get_fe().n_dofs_per_cell());
1559 *  
1560 *   local_dof_indices.resize(cell_dofs->get_fe().n_dofs_per_cell());
1561 *   cell_dofs->get_dof_indices(local_dof_indices);
1562 *   AffineConstraints<double>().distribute_local_to_global(
1563 *   local_values, local_dof_indices, force_vector);
1564 *   }
1565 *   };
1566 *  
1567 * @endcode
1568 *
1569 * The callback function is passed together with the vector holding the
1570 * surface-tension force contribution at each quadrature point of the
1571 * immersed mesh to
1573 * missing step is to compress the distributed force vector.
1574 *
1575 * @code
1576 *   rpe.process_and_evaluate<value_type>(integration_values,
1579 *  
1580 * @endcode
1581 *
1582 * After every tenth step or at the beginning/end of the time loop, we
1583 * output the force vector and the velocity of the background mesh to
1584 * a VTU file. In addition, we also export the geometry of the
1585 * (deformed) immersed surface mesh to a separate VTU file.
1586 *
1587 * @code
1588 *   if (it % 10 == 0 || it == n_time_steps)
1589 *   {
1590 *   std::vector<
1594 *   pcout << " - write data (background mesh)" << std::endl;
1597 *   flags_background.write_higher_order_cells = true;
1599 *   data_out_background.add_data_vector(
1601 *   force_vector,
1602 *   std::vector<std::string>(dim, "force"),
1604 *   data_out_background.add_data_vector(
1606 *   velocity,
1607 *   std::vector<std::string>(dim, "velocity"),
1609 *   data_out_background.build_patches(mapping_background, 3);
1610 *   data_out_background.write_vtu_in_parallel("example_3_background_" +
1611 *   std::to_string(it) +
1612 *   ".vtu",
1613 *   MPI_COMM_WORLD);
1614 *  
1615 *   pcout << " - write mesh (immersed mesh)" << std::endl;
1616 *   DataOut<dim - 1, dim> data_out_immersed;
1617 *   data_out_immersed.attach_triangulation(tria_immersed);
1618 *   data_out_immersed.build_patches(mapping_immersed, 3);
1619 *   data_out_immersed.write_vtu_in_parallel("example_3_immersed_" +
1620 *   std::to_string(it) +
1621 *   ".vtu",
1622 *   MPI_COMM_WORLD);
1623 *   }
1624 *   pcout << std::endl;
1625 *   }
1626 *   }
1627 *  
1628 * @endcode
1629 *
1630 *
1631 * <a name="step_87-Utilityfunctionsdefinition"></a>
1633 *
1634 * @code
1635 *   template <int spacedim>
1636 *   std::vector<Point<spacedim>>
1638 *   const Point<spacedim> &p1,
1639 *   const unsigned int n_subdivisions)
1640 *   {
1641 *   Assert(n_subdivisions >= 1, ExcInternalError());
1642 *  
1643 *   std::vector<Point<spacedim>> points;
1644 *   points.reserve(n_subdivisions + 1);
1645 *  
1646 *   points.emplace_back(p0);
1647 *   for (unsigned int i = 1; i < n_subdivisions; ++i)
1648 *   points.emplace_back(p0 + (p1 - p0) * static_cast<double>(i) /
1649 *   static_cast<double>(n_subdivisions));
1650 *   points.emplace_back(p1);
1651 *  
1652 *   return points;
1653 *   }
1654 *  
1655 *   template <int spacedim, typename T0, typename T1>
1656 *   void print_along_line(const std::string &file_name,
1657 *   const std::vector<Point<spacedim>> &points,
1658 *   const std::vector<T0> &values_0,
1659 *   const std::vector<T1> &values_1)
1660 *   {
1661 *   AssertThrow(points.size() == values_0.size() &&
1662 *   (values_1.size() == points.size() || values_1.empty()),
1663 *   ExcMessage("The provided vectors must have the same length."));
1664 *  
1665 *   std::ofstream file(file_name);
1666 *  
1667 *   for (unsigned int i = 0; i < points.size(); ++i)
1668 *   {
1669 *   file << std::fixed << std::right << std::setw(5) << std::setprecision(3)
1670 *   << points[0].distance(points[i]);
1671 *  
1672 *   for (unsigned int d = 0; d < spacedim; ++d)
1673 *   file << std::fixed << std::right << std::setw(10)
1674 *   << std::setprecision(3) << points[i][d];
1675 *  
1676 *   file << std::fixed << std::right << std::setw(10)
1677 *   << std::setprecision(3) << values_0[i];
1678 *  
1679 *   if (!values_1.empty())
1680 *   file << std::fixed << std::right << std::setw(10)
1681 *   << std::setprecision(3) << values_1[i];
1682 *   file << std::endl;
1683 *   }
1684 *   }
1685 *  
1686 *   template <int dim, int spacedim>
1688 *   const Mapping<dim, spacedim> &mapping,
1689 *   const DoFHandler<dim, spacedim> &dof_handler,
1691 *   {
1692 *   support_points.reinit(dof_handler.locally_owned_dofs(),
1694 *   dof_handler.get_mpi_communicator());
1695 *  
1696 *   const auto &fe = dof_handler.get_fe();
1697 *  
1698 *   FEValues<dim, spacedim> fe_values(
1699 *   mapping,
1700 *   fe,
1701 *   fe.base_element(0).get_unit_support_points(),
1703 *  
1704 *   std::vector<types::global_dof_index> local_dof_indices(
1705 *   fe.n_dofs_per_cell());
1706 *  
1707 *   std::vector<unsigned int> component_to_system_index(
1708 *   fe_values.n_quadrature_points * spacedim);
1709 *  
1710 *   for (unsigned int q = 0, c = 0; q < fe_values.n_quadrature_points; ++q)
1711 *   for (unsigned int d = 0; d < spacedim; ++d, ++c)
1712 *   component_to_system_index[c] = fe.component_to_system_index(d, q);
1713 *  
1714 *   for (const auto &cell : dof_handler.active_cell_iterators() |
1715 *   IteratorFilters::LocallyOwnedCell())
1716 *   {
1717 *   fe_values.reinit(cell);
1718 *   cell->get_dof_indices(local_dof_indices);
1719 *  
1720 *   for (unsigned int q = 0, c = 0; q < fe_values.n_quadrature_points; ++q)
1721 *   for (unsigned int d = 0; d < spacedim; ++d, ++c)
1722 *   support_points[local_dof_indices[component_to_system_index[c]]] =
1723 *   fe_values.quadrature_point(q)[d];
1724 *   }
1725 *   }
1726 *  
1727 *   template <int dim, int spacedim, typename T>
1728 *   std::tuple<std::vector<Point<spacedim>>, std::vector<types::global_dof_index>>
1730 *   const Mapping<dim, spacedim> &mapping,
1732 *   const LinearAlgebra::distributed::Vector<T> &signed_distance,
1734 *   const double narrow_band_threshold)
1735 *   {
1737 *   ExcMessage("The narrow band threshold"
1738 *   " must be larger than or equal to 0."));
1739 *   const auto &tria = dof_handler_signed_distance.get_triangulation();
1740 *   const Quadrature<dim> quad(dof_handler_support_points.get_fe()
1741 *   .base_element(0)
1742 *   .get_unit_support_points());
1743 *  
1746 *   quad,
1747 *   update_values);
1748 *  
1749 *   FEValues<dim> req_values(mapping,
1750 *   dof_handler_support_points.get_fe(),
1751 *   quad,
1753 *  
1754 *   std::vector<T> temp_distance(quad.size());
1755 *   std::vector<types::global_dof_index> local_dof_indices(
1756 *   dof_handler_support_points.get_fe().n_dofs_per_cell());
1757 *  
1758 *   std::vector<Point<dim>> support_points;
1759 *   std::vector<types::global_dof_index> support_points_idx;
1760 *  
1761 *   const bool has_ghost_elements = signed_distance.has_ghost_elements();
1762 *  
1763 *   const auto &locally_owned_dofs_req =
1764 *   dof_handler_support_points.locally_owned_dofs();
1765 *   std::vector<bool> flags(locally_owned_dofs_req.n_elements(), false);
1766 *  
1767 *   if (has_ghost_elements == false)
1768 *   signed_distance.update_ghost_values();
1769 *  
1770 *   for (const auto &cell :
1771 *   tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
1772 *   {
1774 *   cell->as_dof_handler_iterator(dof_handler_signed_distance);
1776 *   distance_values.get_function_values(signed_distance, temp_distance);
1777 *  
1778 *   const auto cell_req =
1779 *   cell->as_dof_handler_iterator(dof_handler_support_points);
1780 *   req_values.reinit(cell_req);
1781 *   cell_req->get_dof_indices(local_dof_indices);
1782 *  
1783 *   for (const auto q : req_values.quadrature_point_indices())
1785 *   {
1786 *   const auto idx = local_dof_indices[q];
1787 *  
1788 *   if (locally_owned_dofs_req.is_element(idx) == false ||
1789 *   flags[locally_owned_dofs_req.index_within_set(idx)])
1790 *   continue;
1791 *  
1792 *   flags[locally_owned_dofs_req.index_within_set(idx)] = true;
1793 *  
1794 *   support_points_idx.emplace_back(idx);
1795 *   support_points.emplace_back(req_values.quadrature_point(q));
1796 *   }
1797 *   }
1798 *  
1799 *   if (has_ghost_elements == false)
1800 *   signed_distance.zero_out_ghost_values();
1801 *  
1803 *   }
1804 *  
1805 *   template <int spacedim>
1806 *   std::vector<Point<spacedim>> convert(
1808 *   {
1809 *   const unsigned int n_points =
1810 *   support_points_unrolled.locally_owned_size() / spacedim;
1811 *  
1812 *   std::vector<Point<spacedim>> points(n_points);
1813 *  
1814 *   for (unsigned int i = 0, c = 0; i < n_points; ++i)
1815 *   for (unsigned int d = 0; d < spacedim; ++d, ++c)
1816 *   points[i][d] = support_points_unrolled.local_element(c);
1817 *  
1818 *   return points;
1819 *   }
1820 *  
1821 *   } // namespace Step87
1822 *  
1823 * @endcode
1824 *
1825 *
1826 * <a name="step_87-Driver"></a>
1827 * <h3>Driver</h3>
1828 *
1829
1830 *
1832 *
1833 * @code
1834 *   int main(int argc, char **argv)
1835 *   {
1836 *   using namespace dealii;
1838 *   std::cout.precision(5);
1839 *  
1841 *   Step87::example_0(); // only run on root process
1842 *  
1846 *   }
1847 * @endcode
1848<a name="step_87-Results"></a><h1>Results</h1>
1849
1850
1851<a name="step_87-Miniexample0"></a><h3>Mini example 0</h3>
1852
1853
1854We present a part of the terminal output. It shows, for each point, the
1855determined cell and reference position. Also, one can see that
1858
1859@verbatim
1861 - Found point with real coordinates: 0 0.5
1862 - in cell with vertices: (0 0.4) (0.2 0.4) (0 0.6) (0.2 0.6)
1863 - with coordinates on the unit cell: (0 0.5)
1864 - Values at point:
1865 - 0.25002 (w. FEValues)
1866 - 0.25002 (w. FEPointEvaluation)
1867 - 0.25002 (w. VectorTools::point_value())
1868
1869 - Found point with real coordinates: 0.05 0.5
1870 - in cell with vertices: (0 0.4) (0.2 0.4) (0 0.6) (0.2 0.6)
1871 - with coordinates on the unit cell: (0.25 0.5)
1872 - Values at point:
1873 - 0.20003 (w. FEValues)
1874 - 0.20003 (w. FEPointEvaluation)
1875 - 0.20003 (w. VectorTools::point_value())
1876
1877...
1878
1879 - Found point with real coordinates: 1 0.5
1880 - in cell with vertices: (0.8 0.4) (1 0.4) (0.8 0.6) (1 0.6)
1881 - with coordinates on the unit cell: (1 0.5)
1882 - Values at point:
1883 - 0.25002 (w. FEValues)
1884 - 0.25002 (w. FEPointEvaluation)
1885 - 0.25002 (w. VectorTools::point_value())
1886
1887 - writing csv file
1889
1890The CSV output is as follows and contains, in the
1891first column, the distances with respect to the first point,
1893of the points and the fourth column the evaluated solution
1894values at those points.
1895
1896@verbatim
18970.000 0.000 0.500 0.250
18980.050 0.050 0.500 0.200
18990.100 0.100 0.500 0.150
19000.150 0.150 0.500 0.100
19010.200 0.200 0.500 0.050
19020.250 0.250 0.500 0.000
19030.300 0.300 0.500 -0.050
19040.350 0.350 0.500 -0.100
19050.400 0.400 0.500 -0.149
19060.450 0.450 0.500 -0.200
19070.500 0.500 0.500 -0.222
19080.550 0.550 0.500 -0.200
19090.600 0.600 0.500 -0.149
19100.650 0.650 0.500 -0.100
19110.700 0.700 0.500 -0.050
19120.750 0.750 0.500 0.000
19130.800 0.800 0.500 0.050
19140.850 0.850 0.500 0.100
19150.900 0.900 0.500 0.150
19160.950 0.950 0.500 0.200
19171.000 1.000 0.500 0.250
1919
1920<a name="step_87-Miniexample1"></a><h3>Mini example 1</h3>
1921
1922
1923We show the terminal output.
1924
1925@verbatim
1927 - writing csv file
1929
1930The CSV output is as follows and identical to the results
1931of the serial case presented in mini example 0.
1932The fifth column shows the
1934
1935@verbatim
19360.000 0.000 0.500 0.250 0.000
19370.050 0.050 0.500 0.200 0.050
19380.100 0.100 0.500 0.150 0.100
19390.150 0.150 0.500 0.100 0.150
19400.200 0.200 0.500 0.050 0.200
19410.250 0.250 0.500 0.000 0.250
19420.300 0.300 0.500 -0.050 0.300
19430.350 0.350 0.500 -0.100 0.350
19440.400 0.400 0.500 -0.149 0.400
19450.450 0.450 0.500 -0.200 0.450
19460.500 0.500 0.500 -0.222 0.500
19470.550 0.550 0.500 -0.200 0.550
19480.600 0.600 0.500 -0.149 0.600
19490.650 0.650 0.500 -0.100 0.650
19500.700 0.700 0.500 -0.050 0.700
19510.750 0.750 0.500 0.000 0.750
19520.800 0.800 0.500 0.050 0.800
19530.850 0.850 0.500 0.100 0.850
19540.900 0.900 0.500 0.150 0.900
19550.950 0.950 0.500 0.200 0.950
19561.000 1.000 0.500 0.250 1.000
1958
1959
1960<a name="step_87-Miniexample2"></a><h3>Mini example 2</h3>
1961
1962
1963We show the terminal output.
1964@verbatim
1966 - create system
1969 - iteration 0: 7076 -> 7076
1970 - iteration 1: 7076 -> 104
1971 - iteration 2: 104 -> 0
1972 - determine distance in narrow band
1974 - output results
1976
1978closest-point projection, show the current position of the closest
1979points exceeding the required tolerance of the discrete interface
1984<table align="center" class="doxtable">
1985 <tr>
1986 <td>
1987 @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_p_0.png
1988 </td>
1989 <td>
1990 @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_p_1.png
1991 </td>
1992 <td>
1993 @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_p_2.png
1994 </td>
1995 </tr>
1996</table>
1997
1999left, the original distance function is shown as the light gray surface.
2000In addition, the contour values refer to the distance values determined
2001from calculation of the distance to the closest points at the interface
2002in the narrow band. It can be seen that the two functions coincide.
2003Similarly, on the right, the original solution and the extrapolated
2004solution from the interface is shown.
2005
2006<table align="center" class="doxtable">
2007 <tr>
2008 <td>
2009 @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_res_0.png
2010 </td>
2011 <td>
2012 @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_res_1.png
2013 </td>
2014 </tr>
2015</table>
2016
2017<a name="step_87-Miniexample3"></a><h3>Mini example 3</h3>
2018
2019
2020We show a shortened version of the terminal output.
2021
2022@verbatim
2024 - creating background mesh
2026time: 0
2027 - compute to be tested values (immersed mesh)
2028 - test values (background mesh)
2029 - write data (background mesh)
2030 - write mesh (immersed mesh)
2031
2032time: 0.01
2033 - move support points (immersed mesh)
2034 - compute to be tested values (immersed mesh)
2035 - test values (background mesh)
2036
2037time: 0.02
2038 - move support points (immersed mesh)
2039 - compute to be tested values (immersed mesh)
2040 - test values (background mesh)
2041
2042...
2043
2044time: 2
2045 - move support points (immersed mesh)
2046 - compute to be tested values (immersed mesh)
2047 - test values (background mesh)
2048 - write data (background mesh)
2049 - write mesh (immersed mesh)
2051
2055in the right figure. The sharp nature of the surface-tension force vector, shown
2056as vector plots, can be seen by its restriction to cells that are intersected by
2058
2059<table align="center" class="doxtable">
2060 <tr>
2061 <td>
2062 @image html https://www.dealii.org/images/steps/developer/step_87_ex_3_force.0000.png
2063 </td>
2064 <td>
2065 @image html https://www.dealii.org/images/steps/developer/step_87_ex_3_force.0010.png
2066 </td>
2067 <td>
2068 @image html https://www.dealii.org/images/steps/developer/step_87_ex_3_force.0020.png
2069 </td>
2070 </tr>
2071</table>
2072
2073<a name="step_87-Possibilitiesforextension"></a><h3>Possibilities for extension</h3>
2074
2075
2079- Performing a distributed search is an expensive step. That is why we suggest
2080to provide hints to Utilities::MPI::RemotePointEvaluation and to reuse
2081Utilities::MPI::RemotePointEvaluation
2084already sorted into the right cells. This is the case if the points are
2085generated on the cell level (see @ref step_85 "step-85"; CutFEM) or the points are
2086automatically sorted into the correct (neighboring) cell (see @ref step_68 "step-68"; PIC with
2087Particles::ParticleHandler). Having said that, the
2088Particles::ParticleHandler::insert_global_particles() function uses
2090cells.
2092like to point out the need for fast evaluation on cell level.
2094@f[
2095\hat{u}(\hat{\boldsymbol{x}}) = \sum_i \hat{N}_i(\hat{\boldsymbol{x}}) \hat{u}_i
2096@f]
2097to derive fast implementations, e.g., for tensor-product elements
2098@f[
2099\hat{u}(\hat{x}_0, \hat{x}_1, \hat{x}_2) =
2100\sum_k \hat{N}^{\text{1D}}_k(\hat{x}_2)
2101\sum_j \hat{N}^{\text{1D}}_j(\hat{x}_1)
2102\sum_i \hat{N}^{\text{1D}}_i(\hat{x}_0)
2103\hat{u}_{ijk}.
2104@f]
2108The class DataOutResample allows to output results on a different mesh than
2110on a coarser mesh or one does not want to output 3D results but instead 2D
2112and restrict residuals between two independent meshes. By passing a sequence
2114non-nested multigrid can be computed.
2115- Utilities::MPI::RemotePointEvaluation can be used to couple non-matching
2116grids via surfaces (example: fluid-structure interaction, independently created
2117grids). The evaluation points can come from any side (pointwise interpolation)
2118or are defined on intersected meshes (Nitsche-type mortaring
2119@cite heinz2022high). Concerning the creation of such intersected meshes and the
2121GridTools::internal::distributed\_compute_intersection_locations().
2126 *
2127 *
2128<a name="step_87-PlainProg"></a>
2129<h1> The plain program</h1>
2130@include "step-87.cc"
2131*/
void test_and_sum(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
Definition fe_q.h:554
Definition point.h:113
friend class Tensor
Definition tensor.h:865
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
Point< 2 > second
Definition grid_out.cc:4630
Point< 2 > first
Definition grid_out.cc:4629
unsigned int level
Definition grid_out.cc:4632
#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())
Definition loop.h:564
@ 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.
std::vector< index_type > data
Definition mpi.cc:740
std::size_t size
Definition mpi.cc:739
const Event initial
Definition event.cc:68
void support_point_wise(DoFHandler< dim, spacedim > &dof_handler)
IndexSet extract_locally_active_dofs(const DoFHandler< dim, spacedim > &dof_handler)
void extrapolate(const DoFHandler< dim, spacedim > &dof1, const InVector &z1, const DoFHandler< dim, spacedim > &dof2, OutVector &z2)
void hyper_sphere(Triangulation< spacedim - 1, spacedim > &tria, const Point< spacedim > &center=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.
constexpr char N
constexpr char T
constexpr types::blas_int zero
constexpr types::blas_int one
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:193
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)
Definition mpi.cc:112
std::string compress(const std::string &input)
Definition utilities.cc:393
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={}, const unsigned int level=numbers::invalid_unsigned_int)
std::vector< typename FEPointEvaluation< n_components, dim, spacedim, typename VectorType::value_type >::value_type > point_values(const Mapping< dim > &mapping, const MeshType< dim, spacedim > &mesh, const VectorType &vector, const std::vector< Point< spacedim > > &evaluation_points, Utilities::MPI::RemotePointEvaluation< dim, spacedim > &cache, const EvaluationFlags::EvaluationFlags flags=EvaluationFlags::avg, const unsigned int first_selected_component=0)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
STL namespace.
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition types.h:32