Reference documentation for deal.II version 9.6.0
\(\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
1003);
1004 *  
1005 * @endcode
1006 *
1007 * For the iterative solution procedure of the closest-point projection,
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 *
1018 * Now, we are ready to perform the algorithm by setting an initial guess
1019 * for the projection points simply corresponding to the collected support
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
1022 * fulfill the required tolerance. Those will be gradually reduced
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 =
1032 *   Utilities::MPI::sum(unmatched_points_idx.size(), MPI_COMM_WORLD);
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
1038 * that still need to be processed and subsequently pass this information
1039 * to the Utilities::MPI::RemotePointEvaluation object. For the sake of
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)
1056 *   unmatched_points[i] = closest_points[unmatched_points_idx[i]];
1057 *  
1058 *   const auto all_unmatched_points =
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 *  
1066 *   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
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 *   {
1090 *   closest_points[unmatched_points_idx[i]] -=
1091 *   eval_values[i] * eval_gradient[i];
1092 *  
1093 *   unmatched_points_idx_next.emplace_back(unmatched_points_idx[i]);
1094 *   }
1095 *  
1096 *   unmatched_points_idx.swap(unmatched_points_idx_next);
1097 *  
1098 *   n_unmatched_points =
1099 *   Utilities::MPI::sum(unmatched_points_idx.size(), MPI_COMM_WORLD);
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
1118 * closest points at the zero-isosurface level set. This information
1119 * can be used to update the signed distance function, i.e., the
1120 * reinitialization the values of the level-set function to maintain
1121 * the signed distance property @cite henri2022geometrical.
1122 *
1123 * @code
1124 *   pcout << " - determine distance in narrow band" << std::endl;
1125 *   LinearAlgebra::distributed::Vector<double> solution_distance;
1126 *   solution_distance.reinit(solution);
1127 *  
1128 *   for (unsigned int i = 0; i < closest_points.size(); ++i)
1129 *   solution_distance[support_points_idx[i]] =
1130 *   support_points[i].distance(closest_points[i]);
1131 *  
1132 * @endcode
1133 *
1134 * In addition, we use the information of the closest point to
1135 * extrapolate values from the interface, i.e., the zero-level
1136 * set isosurface, to the support points within the narrow band.
1137 * This might be helpful to improve accuracy, e.g., for
1138 * diffuse interface fluxes where certain quantities are only
1139 * accurately determined at the interface (e.g. curvature
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 *  
1148 *   LinearAlgebra::distributed::Vector<double> solution_extrapolated;
1149 *   solution_extrapolated.reinit(solution);
1150 *  
1151 *   for (unsigned int i = 0; i < closest_points.size(); ++i)
1152 *   solution_extrapolated[support_points_idx[i]] = vals[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,
1164 *   solution_distance,
1165 *   "solution_distance");
1166 *   data_out.add_data_vector(dof_handler,
1167 *   solution_extrapolated,
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 *
1183 * The final mini example presents a basic implementation of
1184 * front tracking @cite peskin1977numerical, @cite unverdi1992front
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 *
1190 * We assume that the immersed surface is transported according to a
1191 * prescribed velocity field from the background mesh. Subsequently,
1192 * we perform a sharp computation of the surface-tension force:
1193 * @f[
1194 * (\boldsymbol v_i (\boldsymbol{x}), \boldsymbol F_S
1195 * (\boldsymbol{x}))_{\Omega}
1196 * =
1197 * \left( \boldsymbol{v}_i (\boldsymbol{x}), \sigma (\boldsymbol{x}) \kappa
1198 * (\boldsymbol{x}) \boldsymbol{n} (\boldsymbol{x}) \right)_\Gamma \approx
1199 * \sum_{q\in\mathbb{T}_\Gamma} \boldsymbol{v}_i^T (\boldsymbol{x}_q)
1200 * \sigma (\boldsymbol{x}_q) \kappa (\boldsymbol{x}_q) \boldsymbol{n}
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
1206 * the force contributions @f$\sigma (\boldsymbol{x}_q) \kappa
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[
1212 * \boldsymbol{F}_S (\boldsymbol{x}_q) \gets \sigma (\boldsymbol{x}_q) \kappa
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]
1216 * In the second step, we compute the discretized weak form by multiplying
1217 * with test functions on the background mesh:
1218 * @f[
1219 * (\boldsymbol v_i (\boldsymbol{x}), \boldsymbol F_S
1220 * (\boldsymbol{x}))_{\Omega} \gets \sum_q \boldsymbol{v}_i^T
1221 * (\boldsymbol{x}_q) \boldsymbol{F}_S
1222 * (\boldsymbol{x}_q)
1223 * \quad \forall i\in\mathbb{T}_\Omega
1224 * .
1225 * @f]
1226 * Obviously, we need to communicate between the two steps. The second step
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,
1258 *   Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) ==
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;
1270 *   DistributedTriangulation<dim> tria_background(MPI_COMM_WORLD);
1271 *   GridGenerator::hyper_cube(tria_background);
1272 *   tria_background.refine_global(5);
1273 *  
1274 *   MappingQ1<dim> mapping_background;
1275 *   const FESystem<dim> fe_background(FE_Q<dim>(degree), dim);
1276 *   DoFHandler<dim> dof_handler_background(tria_background);
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 *  
1291 *   DistributedTriangulation<dim - 1, dim> tria_immersed(MPI_COMM_WORLD);
1292 *   GridGenerator::hyper_sphere(tria_immersed, center, radius);
1293 *   tria_immersed.refine_global(4);
1294 *  
1295 * @endcode
1296 *
1297 * Two different mappings are considered for the immersed
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
1301 * represent scalar and vector-valued DoF values.
1302 *
1303 * @code
1304 *   MappingQ<dim - 1, dim> mapping_immersed_base(3);
1305 *   MappingQCache<dim - 1, dim> mapping_immersed(3);
1306 *   mapping_immersed.initialize(mapping_immersed_base, tria_immersed);
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);
1311 *   DoFHandler<dim - 1, dim> dof_handler_immersed(tria_immersed);
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
1320 *   DoFRenumbering::support_point_wise(dof_handler_immersed);
1321 *  
1322 * @endcode
1323 *
1324 * We fill a DoF vector on the background mesh with an analytical
1325 * velocity field considering the Rayleigh-Kothe vortex flow and
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(),
1332 *   dof_handler_background),
1333 *   MPI_COMM_WORLD);
1335 *  
1337 *   dof_handler_background.locally_owned_dofs(),
1338 *   DoFTools::extract_locally_active_dofs(dof_handler_background),
1339 *   MPI_COMM_WORLD);
1340 *  
1341 * @endcode
1342 *
1343 * Next, we collect the real positions @f$\boldsymbol{x}_q@f$ of the quadrature
1344 * points of the surface mesh in a vector.
1345 *
1346 * @code
1347 *   LinearAlgebra::distributed::Vector<double> immersed_support_points;
1348 *   collect_support_points(mapping_immersed,
1349 *   dof_handler_immersed,
1350 *   immersed_support_points);
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
1358 * function. Then, we update the internal data structures of the
1359 * Utilities::MPI::RemotePointEvaluation object with the collected support
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
1363 * resulting update of the coordinates. Finally, we update the mapping of
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);
1376 *   VectorTools::interpolate(mapping_background,
1377 *   dof_handler_background,
1378 *   vortex,
1379 *   velocity);
1380 *   rpe.reinit(convert<dim>(immersed_support_points),
1381 *   tria_background,
1382 *   mapping_background);
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 =
1391 *   dof_handler_background,
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 *  
1402 *   mapping_immersed.initialize(mapping_immersed_base,
1403 *   dof_handler_immersed,
1404 *   immersed_support_points,
1405 *   false);
1406 *   }
1407 *  
1408 * @endcode
1409 *
1410 * Next, we loop over all locally owned cells of the immersed mesh and
1411 * its quadrature points to compute the value for the local surface
1412 * tension force contribution @f$\boldsymbol{F}_S(\boldsymbol{x}_q)@f$. We
1413 * store the real coordinates of the quadrature points and the
1414 * corresponding force contributions in two individual vectors. For
1415 * computation of the latter, the normal vector
1416 * @f$\boldsymbol{n}(\boldsymbol{x}_q)@f$ can be directly extracted from the
1417 * surface mesh via FEValues and, for the curvature, we use the
1418 * following approximation:
1419 * @f[
1420 * \kappa(\boldsymbol{x}_q)
1421 * =
1422 * \nabla \cdot \boldsymbol{n}(\boldsymbol{x}_q)
1423 * =
1424 * \text{tr}\left({\nabla \boldsymbol{n}(\boldsymbol{x}_q)}\right)
1425 * \approx
1426 * \text{tr}\left({\nabla \sum_i \boldsymbol{N}_i (\boldsymbol{x}_q)
1427 * \boldsymbol n_i}\right)
1428 * =
1429 * \sum_i\text{tr}\left({\nabla \boldsymbol{N}_i (\boldsymbol{x}_q)
1430 * \boldsymbol n_i}\right)
1431 * \;\text{with}\; i\in[0,n_{\text{dofs\_per\_cell}}),
1432 * @f]
1433 * which we can apply since the immersed mesh is consistently
1434 * orientated. The surface tension coefficient is set to 1 for the
1435 * sake of demonstration.
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,
1446 *   quadrature_immersed,
1450 *  
1451 *   FEValues<dim - 1, dim> fe_values_co(
1452 *   mapping_immersed,
1453 *   fe_scalar_immersed,
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
1499 * communication pattern of Utilities::MPI::RemotePointEvaluation is
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 *  
1506 *   rpe.reinit(integration_points, tria_background, mapping_background);
1507 *  
1508 * @endcode
1509 *
1510 * In preparation for utilizing
1512 * performs the
1513 * multiplication with the test function, we set up a callback function
1514 * that contains the operation on the intersected cells of the
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
1523 * reinitialize the data structure of FEPointEvaluation on every cell
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
1527 * FEPointEvaluation::test_and_sum(), the submitted values are
1528 * multiplied by the values of the test function and a summation over
1529 * all given points is performed. Subsequently, the contributions are
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) {
1536 *   FEPointEvaluation<dim, dim> phi_force(mapping_background,
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 *   {
1545 *   const auto cell_dofs =
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());
1558 *   phi_force.test_and_sum(local_values, EvaluationFlags::values);
1559 *  
1560 *   local_dof_indices.resize(cell_dofs->get_fe().n_dofs_per_cell());
1561 *   cell_dofs->get_dof_indices(local_dof_indices);
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,
1577 *   integration_function);
1578 *   force_vector.compress(VectorOperation::add);
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<
1592 *   vector_component_interpretation(
1594 *   pcout << " - write data (background mesh)" << std::endl;
1595 *   DataOut<dim> data_out_background;
1596 *   DataOutBase::VtkFlags flags_background;
1597 *   flags_background.write_higher_order_cells = true;
1598 *   data_out_background.set_flags(flags_background);
1599 *   data_out_background.add_data_vector(
1600 *   dof_handler_background,
1601 *   force_vector,
1602 *   std::vector<std::string>(dim, "force"),
1603 *   vector_component_interpretation);
1604 *   data_out_background.add_data_vector(
1605 *   dof_handler_background,
1606 *   velocity,
1607 *   std::vector<std::string>(dim, "velocity"),
1608 *   vector_component_interpretation);
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>
1632 * <h3>Utility functions (definition)</h3>
1633 *
1634 * @code
1635 *   template <int spacedim>
1636 *   std::vector<Point<spacedim>>
1637 *   create_points_along_line(const Point<spacedim> &p0,
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>
1687 *   void collect_support_points(
1688 *   const Mapping<dim, spacedim> &mapping,
1689 *   const DoFHandler<dim, spacedim> &dof_handler,
1690 *   LinearAlgebra::distributed::Vector<double> &support_points)
1691 *   {
1692 *   support_points.reinit(dof_handler.locally_owned_dofs(),
1693 *   DoFTools::extract_locally_active_dofs(dof_handler),
1694 *   dof_handler.get_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>>
1729 *   collect_support_points_with_narrow_band(
1730 *   const Mapping<dim, spacedim> &mapping,
1731 *   const DoFHandler<dim, spacedim> &dof_handler_signed_distance,
1732 *   const LinearAlgebra::distributed::Vector<T> &signed_distance,
1733 *   const DoFHandler<dim, spacedim> &dof_handler_support_points,
1734 *   const double narrow_band_threshold)
1735 *   {
1736 *   AssertThrow(narrow_band_threshold >= 0,
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 *  
1744 *   FEValues<dim> distance_values(mapping,
1745 *   dof_handler_signed_distance.get_fe(),
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 *   {
1773 *   const auto cell_distance =
1774 *   cell->as_dof_handler_iterator(dof_handler_signed_distance);
1775 *   distance_values.reinit(cell_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())
1784 *   if (std::abs(temp_distance[q]) < narrow_band_threshold)
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 *  
1802 *   return {support_points, support_points_idx};
1803 *   }
1804 *  
1805 *   template <int spacedim>
1806 *   std::vector<Point<spacedim>> convert(
1807 *   const LinearAlgebra::distributed::Vector<double> &support_points_unrolled)
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 *
1831 * Finally, the driver of the program executes the four mini examples.
1832 *
1833 * @code
1834 *   int main(int argc, char **argv)
1835 *   {
1836 *   using namespace dealii;
1837 *   Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
1838 *   std::cout.precision(5);
1839 *  
1840 *   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
1841 *   Step87::example_0(); // only run on root process
1842 *  
1843 *   Step87::example_1();
1844 *   Step87::example_2();
1845 *   Step87::example_3();
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
1856the values evaluated with FEValues, FEPointEvaluation, and
1857VectorTools::point_values() are identical, as expected.
1858
1859@verbatim
1860Running: example 0
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
1888@endverbatim
1889
1890The CSV output is as follows and contains, in the
1891first column, the distances with respect to the first point,
1892the second and the third column represent the coordinates
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
1918@endverbatim
1919
1920<a name="step_87-Miniexample1"></a><h3>Mini example 1</h3>
1921
1922
1923We show the terminal output.
1924
1925@verbatim
1926Running: example 1
1927 - writing csv file
1928@endverbatim
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
1933user quantity evaluated additionally in this mini example.
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
1957@endverbatim
1958
1959
1960<a name="step_87-Miniexample2"></a><h3>Mini example 2</h3>
1961
1962
1963We show the terminal output.
1964@verbatim
1965Running: example 2
1966 - create system
1967 - determine narrow band
1968 - determine closest point iteratively
1969 - iteration 0: 7076 -> 7076
1970 - iteration 1: 7076 -> 104
1971 - iteration 2: 104 -> 0
1972 - determine distance in narrow band
1973 - perform extrapolation in narrow band
1974 - output results
1975@endverbatim
1976
1977The following three plots, representing the performed iterations of the
1978closest-point projection, show the current position of the closest
1979points exceeding the required tolerance of the discrete interface
1980of the circle and still need to
1981be corrected.
1982It can be seen that the numbers of points to be processed decrease
1983from iteration to iteration.
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
1998The output visualized in Paraview looks like the following: On the
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
2023Running: example 3
2024 - creating background mesh
2025 - creating immersed 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)
2050@endverbatim
2051
2052The output visualized in Paraview looks like the following: The deformation of
2053the immersed mesh by the reversible vortex flow can be seen. Due to
2054discretization errors, the shape is not exactly circular at the end, illustrated
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
2057the immersed mesh.
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
2076This program highlights some of the main capabilities
2077of the distributed evaluation routines in deal.II. However, there are many
2078related topics worth mentioning:
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
2082instances in the case that the communication pattern has not changed.
2083Furthermore, there are instances where no search is needed and the points are
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
2089the described infrastructure to perform the initial sorting of particles into
2090cells.
2091- We concentrated on parallelization aspects in this tutorial. However, we would
2092like to point out the need for fast evaluation on cell level.
2093The task for this in deal.II is FEPointEvaluation. It exploits the structure of
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]
2105Since only 1D shape functions are queried and are re-used in re-occurring terms,
2106this formulation is faster than without exploitation of the structure.
2107- Utilities::MPI::RemotePointEvaluation is used in multiple places in deal.II.
2108The class DataOutResample allows to output results on a different mesh than
2109the computational mesh. This is useful if one wants to output the results
2110on a coarser mesh or one does not want to output 3D results but instead 2D
2111slices. In addition, MGTwoLevelTransferNonNested allows to prolongate solutions
2112and restrict residuals between two independent meshes. By passing a sequence
2113of such two-level transfer operators to MGTransferMF and, finally, to Multigrid,
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
2120corresponding evaluation points, see
2121GridTools::internal::distributed\_compute_intersection_locations().
2122- Alternatively to the coupling via Utilities::MPI::RemotePointEvaluation,
2123preCICE @cite bungartz2016precice @cite chourdakis2021precice can be used. The
2124code-gallery program "Laplace equation coupled to an external simulation
2125program" shows how to use this library with deal.II.
2126 *
2127 *
2128<a name="step_87-PlainProg"></a>
2129<h1> The plain program</h1>
2130@include "step-87.cc"
2131*/
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)
Definition fe_q.h:554
void reinit(const size_type size, const bool omit_zeroing_entries=false)
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
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< 3 > center
Point< 3 > vertices[4]
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
unsigned int level
Definition grid_out.cc:4626
__global__ void set(Number *val, const Number s, const size_type N)
#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.
const Event initial
Definition event.cc:64
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.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
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:107
T reduce(const T &local_value, const MPI_Comm comm, const std::function< T(const T &, const T &)> &combiner, const unsigned int root_process=0)
std::string compress(const std::string &input)
Definition utilities.cc:389
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={})
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)
std::vector< typename FEPointEvaluation< n_components, dim, spacedim, typename VectorType::value_type >::gradient_type > point_gradients(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