deal.II version GIT relicensing-1721-g8100761196 2024-08-31 12:30:00+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);
1003 *  
1004 * @endcode
1005 *
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.
1009 *
1010 * @code
1011 *   pcout << " - determine closest point iteratively" << std::endl;
1012 *   constexpr int max_iter = 30;
1013 *   constexpr double tol_distance = 1e-6;
1014 *  
1015 * @endcode
1016 *
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.
1023 *
1024 * @code
1025 *   std::vector<Point<dim>> closest_points = support_points; // initial guess
1026 *  
1027 *   std::vector<unsigned int> unmatched_points_idx(closest_points.size());
1028 *   std::iota(unmatched_points_idx.begin(), unmatched_points_idx.end(), 0);
1029 *  
1030 *   int n_unmatched_points =
1031 *   Utilities::MPI::sum(unmatched_points_idx.size(), MPI_COMM_WORLD);
1032 *  
1033 * @endcode
1034 *
1035 * Now, we create a Utilities::MPI::RemotePointEvaluation cache object and
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
1038 * to the Utilities::MPI::RemotePointEvaluation object. For the sake of
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.
1045 *
1046 * @code
1048 *  
1049 *   for (int it = 0; it < max_iter && n_unmatched_points > 0; ++it)
1050 *   {
1051 *   pcout << " - iteration " << it << ": " << n_unmatched_points;
1052 *  
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]];
1056 *  
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) {
1060 *   auto result = a;
1061 *   result.insert(result.end(), b.begin(), b.end());
1062 *   return result;
1063 *   });
1064 *  
1065 *   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
1066 *   {
1067 *   std::ofstream file("example_2_" + std::to_string(it) + ".csv");
1068 *   for (const auto &p : all_unmatched_points)
1069 *   file << p << std::endl;
1070 *   file.close();
1071 *   }
1072 *  
1073 *   rpe.reinit(unmatched_points, tria, mapping);
1074 *  
1075 *   AssertThrow(rpe.all_points_found(),
1076 *   ExcMessage("Processed point is outside domain."));
1077 *  
1078 *   const auto eval_values =
1079 *   VectorTools::point_values<1>(rpe, dof_handler, signed_distance);
1080 *  
1081 *   const auto eval_gradient =
1082 *   VectorTools::point_gradients<1>(rpe, dof_handler, signed_distance);
1083 *  
1084 *   std::vector<unsigned int> unmatched_points_idx_next;
1085 *  
1086 *   for (unsigned int i = 0; i < unmatched_points_idx.size(); ++i)
1087 *   if (std::abs(eval_values[i]) > tol_distance)
1088 *   {
1089 *   closest_points[unmatched_points_idx[i]] -=
1090 *   eval_values[i] * eval_gradient[i];
1091 *  
1092 *   unmatched_points_idx_next.emplace_back(unmatched_points_idx[i]);
1093 *   }
1094 *  
1095 *   unmatched_points_idx.swap(unmatched_points_idx_next);
1096 *  
1097 *   n_unmatched_points =
1098 *   Utilities::MPI::sum(unmatched_points_idx.size(), MPI_COMM_WORLD);
1099 *  
1100 *   pcout << " -> " << n_unmatched_points << std::endl;
1101 *   }
1102 *  
1103 * @endcode
1104 *
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.
1108 *
1109 * @code
1110 *   if (n_unmatched_points > 0)
1111 *   pcout << "WARNING: The tolerance of " << n_unmatched_points
1112 *   << " points is not yet attained." << std::endl;
1113 *  
1114 * @endcode
1115 *
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.
1121 *
1122 * @code
1123 *   pcout << " - determine distance in narrow band" << std::endl;
1124 *   LinearAlgebra::distributed::Vector<double> solution_distance;
1125 *   solution_distance.reinit(solution);
1126 *  
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]);
1130 *  
1131 * @endcode
1132 *
1133 * In addition, we use the information of the closest point to
1134 * extrapolate values from the interface, i.e., the zero-level
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).
1140 *
1141 * @code
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);
1146 *  
1147 *   LinearAlgebra::distributed::Vector<double> solution_extrapolated;
1148 *   solution_extrapolated.reinit(solution);
1149 *  
1150 *   for (unsigned int i = 0; i < closest_points.size(); ++i)
1151 *   solution_extrapolated[support_points_idx[i]] = vals[i];
1152 *  
1153 * @endcode
1154 *
1155 * Finally, we output the results to a VTU file.
1156 *
1157 * @code
1158 *   pcout << " - output results" << std::endl;
1159 *   DataOut<dim> data_out;
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);
1170 *  
1171 *   pcout << std::endl;
1172 *   }
1173 *  
1174 * @endcode
1175 *
1176 *
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>
1179 *
1180
1181 *
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$.
1186 *
1187
1188 *
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:
1192 * @f[
1193 * (\boldsymbol v_i (\boldsymbol{x}), \boldsymbol F_S
1194 * (\boldsymbol{x}))_{\Omega}
1195 * =
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
1202 * .
1203 * @f]
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)|
1209 * w_q@f$:
1210 * @f[
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.
1214 * @f]
1215 * In the second step, we compute the discretized weak form by multiplying
1216 * with test functions on the background mesh:
1217 * @f[
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
1223 * .
1224 * @f]
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
1229 * purpose.
1230 *
1231
1232 *
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
1236 * time steps.
1237 *
1238
1239 *
1240 *
1241 * @code
1242 *   void example_3()
1243 *   {
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;
1248 *  
1249 * @endcode
1250 *
1251 * This program is intended to be executed in 2D or 3D.
1252 *
1253 * @code
1254 *   static_assert(dim == 2 || dim == 3, "Only implemented for 2D or 3D.");
1255 *  
1256 *   ConditionalOStream pcout(std::cout,
1257 *   Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) ==
1258 *   0);
1259 *  
1260 *   pcout << "Running: example 3" << std::endl;
1261 *  
1262 * @endcode
1263 *
1264 * Next, we create the standard objects necessary for the finite element
1265 * representation of the background mesh
1266 *
1267 * @code
1268 *   pcout << " - creating background mesh" << std::endl;
1269 *   DistributedTriangulation<dim> tria_background(MPI_COMM_WORLD);
1270 *   GridGenerator::hyper_cube(tria_background);
1271 *   tria_background.refine_global(5);
1272 *  
1273 *   const MappingQ1<dim> mapping_background;
1274 *   const FESystem<dim> fe_background(FE_Q<dim>(degree), dim);
1275 *   DoFHandler<dim> dof_handler_background(tria_background);
1276 *   dof_handler_background.distribute_dofs(fe_background);
1277 *  
1278 * @endcode
1279 *
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.
1283 *
1284 * @code
1285 *   pcout << " - creating immersed mesh" << std::endl;
1286 *   const Point<dim> center((dim == 2) ? Point<dim>(0.5, 0.75) :
1287 *   Point<dim>(0.5, 0.75, 0.5));
1288 *   const double radius = 0.15;
1289 *  
1290 *   DistributedTriangulation<dim - 1, dim> tria_immersed(MPI_COMM_WORLD);
1291 *   GridGenerator::hyper_sphere(tria_immersed, center, radius);
1292 *   tria_immersed.refine_global(4);
1293 *  
1294 * @endcode
1295 *
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
1300 * represent scalar and vector-valued DoF values.
1301 *
1302 * @code
1303 *   const MappingQ<dim - 1, dim> mapping_immersed_base(3);
1304 *   MappingQCache<dim - 1, dim> mapping_immersed(3);
1305 *   mapping_immersed.initialize(mapping_immersed_base, tria_immersed);
1306 *   const QGauss<dim - 1> quadrature_immersed(degree + 1);
1307 *  
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);
1312 *  
1313 * @endcode
1314 *
1315 * We renumber the DoFs related to the vector-valued problem to
1316 * simplify access to the individual components.
1317 *
1318 * @code
1319 *   DoFRenumbering::support_point_wise(dof_handler_immersed);
1320 *  
1321 * @endcode
1322 *
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.
1326 *
1327 * @code
1329 *   velocity.reinit(dof_handler_background.locally_owned_dofs(),
1331 *   dof_handler_background),
1332 *   MPI_COMM_WORLD);
1334 *  
1336 *   dof_handler_background.locally_owned_dofs(),
1337 *   DoFTools::extract_locally_active_dofs(dof_handler_background),
1338 *   MPI_COMM_WORLD);
1339 *  
1340 * @endcode
1341 *
1342 * Next, we collect the real positions @f$\boldsymbol{x}_q@f$ of the quadrature
1343 * points of the surface mesh in a vector.
1344 *
1345 * @code
1346 *   LinearAlgebra::distributed::Vector<double> immersed_support_points;
1347 *   collect_support_points(mapping_immersed,
1348 *   dof_handler_immersed,
1349 *   immersed_support_points);
1350 *  
1351 * @endcode
1352 *
1353 * We initialize a Utilities::MPI::RemotePointEvaluation object and start
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
1358 * Utilities::MPI::RemotePointEvaluation object with the collected support
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.
1364 *
1365 * @code
1367 *   double time = 0.0;
1368 *   for (unsigned int it = 0; it <= n_time_steps; ++it, time += dt)
1369 *   {
1370 *   pcout << "time: " << time << std::endl;
1371 *   if (it > 0)
1372 *   {
1373 *   pcout << " - move support points (immersed mesh)" << std::endl;
1374 *   vortex.set_time(time);
1375 *   VectorTools::interpolate(mapping_background,
1376 *   dof_handler_background,
1377 *   vortex,
1378 *   velocity);
1379 *   rpe.reinit(convert<dim>(immersed_support_points),
1380 *   tria_background,
1381 *   mapping_background);
1382 *  
1383 *   AssertThrow(rpe.all_points_found(),
1384 *   ExcMessage(
1385 *   "Immersed domain leaves background domain!"));
1386 *  
1387 *   velocity.update_ghost_values();
1388 *   const auto immersed_velocity =
1389 *   VectorTools::point_values<dim>(rpe,
1390 *   dof_handler_background,
1391 *   velocity);
1392 *   velocity.zero_out_ghost_values();
1393 *  
1394 *   for (unsigned int i = 0, c = 0;
1395 *   i < immersed_support_points.locally_owned_size() / dim;
1396 *   ++i)
1397 *   for (unsigned int d = 0; d < dim; ++d, ++c)
1398 *   immersed_support_points.local_element(c) +=
1399 *   immersed_velocity[i][d] * dt;
1400 *  
1401 *   mapping_immersed.initialize(mapping_immersed_base,
1402 *   dof_handler_immersed,
1403 *   immersed_support_points,
1404 *   false);
1405 *   }
1406 *  
1407 * @endcode
1408 *
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:
1418 * @f[
1419 * \kappa(\boldsymbol{x}_q)
1420 * =
1421 * \nabla \cdot \boldsymbol{n}(\boldsymbol{x}_q)
1422 * =
1423 * \text{tr}\left({\nabla \boldsymbol{n}(\boldsymbol{x}_q)}\right)
1424 * \approx
1425 * \text{tr}\left({\nabla \sum_i \boldsymbol{N}_i (\boldsymbol{x}_q)
1426 * \boldsymbol n_i}\right)
1427 * =
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}}),
1431 * @f]
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.
1435 *
1436 * @code
1437 *   pcout << " - compute to be tested values (immersed mesh)" << std::endl;
1438 *   using value_type = Tensor<1, dim, double>;
1439 *  
1440 *   std::vector<Point<dim>> integration_points;
1441 *   std::vector<value_type> integration_values;
1442 *  
1443 *   FEValues<dim - 1, dim> fe_values(mapping_immersed,
1444 *   fe_immersed,
1445 *   quadrature_immersed,
1449 *  
1450 *   FEValues<dim - 1, dim> fe_values_co(
1451 *   mapping_immersed,
1452 *   fe_scalar_immersed,
1453 *   fe_scalar_immersed.get_unit_support_points(),
1455 *  
1456 *   std::vector<unsigned int> component_to_system_index(
1457 *   fe_immersed.n_dofs_per_cell());
1458 *  
1459 *   for (unsigned int i = 0, c = 0;
1460 *   i < fe_scalar_immersed.n_dofs_per_cell();
1461 *   ++i)
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);
1465 *  
1466 *   for (const auto &cell : tria_immersed.active_cell_iterators() |
1467 *   IteratorFilters::LocallyOwnedCell())
1468 *   {
1469 *   fe_values.reinit(cell);
1470 *   fe_values_co.reinit(cell);
1471 *  
1472 *   for (const auto &q : fe_values.quadrature_point_indices())
1473 *   {
1474 *   integration_points.emplace_back(fe_values.quadrature_point(q));
1475 *  
1476 *   const auto sigma = 1.0; // surface tension coefficient
1477 *  
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();
1482 *   ++i)
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];
1487 *  
1488 *   const auto FxJxW =
1489 *   sigma * curvature * normal * fe_values.JxW(q);
1490 *  
1491 *   integration_values.emplace_back(FxJxW);
1492 *   }
1493 *   }
1494 *  
1495 * @endcode
1496 *
1497 * Before we evaluate the weak form of the surface-tension force, the
1498 * communication pattern of Utilities::MPI::RemotePointEvaluation is
1499 * set up from the quadrature points of the immersed mesh, determining
1500 * the surrounding cells on the background mesh.
1501 *
1502 * @code
1503 *   pcout << " - test values (background mesh)" << std::endl;
1504 *  
1505 *   rpe.reinit(integration_points, tria_background, mapping_background);
1506 *  
1507 * @endcode
1508 *
1509 * In preparation for utilizing
1511 * performs the
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
1515 * FEPointEvaluation object that allows us to integrate values at
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
1522 * reinitialize the data structure of FEPointEvaluation on every cell
1523 * according to the unit points. Next, we loop over the quadrature
1524 * points attributed to the cell and submit the local surface-tension
1525 * force to the FEPointEvaluation object. Via
1526 * FEPointEvaluation::test_and_sum(), the submitted values are
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.
1531 *
1532 * @code
1533 *   const auto integration_function = [&](const auto &values,
1534 *   const auto &cell_data) {
1535 *   FEPointEvaluation<dim, dim> phi_force(mapping_background,
1536 *   fe_background,
1537 *   update_values);
1538 *  
1539 *   std::vector<double> local_values;
1540 *   std::vector<types::global_dof_index> local_dof_indices;
1541 *  
1542 *   for (const auto cell : cell_data.cell_indices())
1543 *   {
1544 *   const auto cell_dofs =
1545 *   cell_data.get_active_cell_iterator(cell)
1546 *   ->as_dof_handler_iterator(dof_handler_background);
1547 *  
1548 *   const auto unit_points = cell_data.get_unit_points(cell);
1549 *   const auto FxJxW = cell_data.get_data_view(cell, values);
1550 *  
1551 *   phi_force.reinit(cell_dofs, unit_points);
1552 *  
1553 *   for (const auto q : phi_force.quadrature_point_indices())
1554 *   phi_force.submit_value(FxJxW[q], q);
1555 *  
1556 *   local_values.resize(cell_dofs->get_fe().n_dofs_per_cell());
1557 *   phi_force.test_and_sum(local_values, EvaluationFlags::values);
1558 *  
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);
1563 *   }
1564 *   };
1565 *  
1566 * @endcode
1567 *
1568 * The callback function is passed together with the vector holding the
1569 * surface-tension force contribution at each quadrature point of the
1570 * immersed mesh to
1572 * missing step is to compress the distributed force vector.
1573 *
1574 * @code
1575 *   rpe.process_and_evaluate<value_type>(integration_values,
1576 *   integration_function);
1577 *   force_vector.compress(VectorOperation::add);
1578 *  
1579 * @endcode
1580 *
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.
1585 *
1586 * @code
1587 *   if (it % 10 == 0 || it == n_time_steps)
1588 *   {
1589 *   std::vector<
1591 *   vector_component_interpretation(
1593 *   pcout << " - write data (background mesh)" << std::endl;
1594 *   DataOut<dim> data_out_background;
1595 *   DataOutBase::VtkFlags flags_background;
1596 *   flags_background.write_higher_order_cells = true;
1597 *   data_out_background.set_flags(flags_background);
1598 *   data_out_background.add_data_vector(
1599 *   dof_handler_background,
1600 *   force_vector,
1601 *   std::vector<std::string>(dim, "force"),
1602 *   vector_component_interpretation);
1603 *   data_out_background.add_data_vector(
1604 *   dof_handler_background,
1605 *   velocity,
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) +
1611 *   ".vtu",
1612 *   MPI_COMM_WORLD);
1613 *  
1614 *   pcout << " - write mesh (immersed mesh)" << std::endl;
1615 *   DataOut<dim - 1, dim> data_out_immersed;
1616 *   data_out_immersed.attach_triangulation(tria_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) +
1620 *   ".vtu",
1621 *   MPI_COMM_WORLD);
1622 *   }
1623 *   pcout << std::endl;
1624 *   }
1625 *   }
1626 *  
1627 * @endcode
1628 *
1629 *
1630 * <a name="step_87-Utilityfunctionsdefinition"></a>
1631 * <h3>Utility functions (definition)</h3>
1632 *
1633 * @code
1634 *   template <int spacedim>
1635 *   std::vector<Point<spacedim>>
1636 *   create_points_along_line(const Point<spacedim> &p0,
1637 *   const Point<spacedim> &p1,
1638 *   const unsigned int n_subdivisions)
1639 *   {
1640 *   Assert(n_subdivisions >= 1, ExcInternalError());
1641 *  
1642 *   std::vector<Point<spacedim>> points;
1643 *   points.reserve(n_subdivisions + 1);
1644 *  
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);
1650 *  
1651 *   return points;
1652 *   }
1653 *  
1654 *   template <int spacedim, typename T0, typename T1>
1655 *   void print_along_line(const std::string &file_name,
1656 *   const std::vector<Point<spacedim>> &points,
1657 *   const std::vector<T0> &values_0,
1658 *   const std::vector<T1> &values_1)
1659 *   {
1660 *   AssertThrow(points.size() == values_0.size() &&
1661 *   (values_1.size() == points.size() || values_1.empty()),
1662 *   ExcMessage("The provided vectors must have the same length."));
1663 *  
1664 *   std::ofstream file(file_name);
1665 *  
1666 *   for (unsigned int i = 0; i < points.size(); ++i)
1667 *   {
1668 *   file << std::fixed << std::right << std::setw(5) << std::setprecision(3)
1669 *   << points[0].distance(points[i]);
1670 *  
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];
1674 *  
1675 *   file << std::fixed << std::right << std::setw(10)
1676 *   << std::setprecision(3) << values_0[i];
1677 *  
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;
1682 *   }
1683 *   }
1684 *  
1685 *   template <int dim, int spacedim>
1686 *   void collect_support_points(
1687 *   const Mapping<dim, spacedim> &mapping,
1688 *   const DoFHandler<dim, spacedim> &dof_handler,
1689 *   LinearAlgebra::distributed::Vector<double> &support_points)
1690 *   {
1691 *   support_points.reinit(dof_handler.locally_owned_dofs(),
1692 *   DoFTools::extract_locally_active_dofs(dof_handler),
1693 *   dof_handler.get_mpi_communicator());
1694 *  
1695 *   const auto &fe = dof_handler.get_fe();
1696 *  
1697 *   FEValues<dim, spacedim> fe_values(
1698 *   mapping,
1699 *   fe,
1700 *   fe.base_element(0).get_unit_support_points(),
1702 *  
1703 *   std::vector<types::global_dof_index> local_dof_indices(
1704 *   fe.n_dofs_per_cell());
1705 *  
1706 *   std::vector<unsigned int> component_to_system_index(
1707 *   fe_values.n_quadrature_points * spacedim);
1708 *  
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);
1712 *  
1713 *   for (const auto &cell : dof_handler.active_cell_iterators() |
1714 *   IteratorFilters::LocallyOwnedCell())
1715 *   {
1716 *   fe_values.reinit(cell);
1717 *   cell->get_dof_indices(local_dof_indices);
1718 *  
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];
1723 *   }
1724 *   }
1725 *  
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(
1729 *   const Mapping<dim, spacedim> &mapping,
1730 *   const DoFHandler<dim, spacedim> &dof_handler_signed_distance,
1731 *   const LinearAlgebra::distributed::Vector<T> &signed_distance,
1732 *   const DoFHandler<dim, spacedim> &dof_handler_support_points,
1733 *   const double narrow_band_threshold)
1734 *   {
1735 *   AssertThrow(narrow_band_threshold >= 0,
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();
1739 *   const Quadrature<dim> quad(dof_handler_support_points.get_fe()
1740 *   .base_element(0)
1741 *   .get_unit_support_points());
1742 *  
1743 *   FEValues<dim> distance_values(mapping,
1744 *   dof_handler_signed_distance.get_fe(),
1745 *   quad,
1746 *   update_values);
1747 *  
1748 *   FEValues<dim> req_values(mapping,
1749 *   dof_handler_support_points.get_fe(),
1750 *   quad,
1752 *  
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());
1756 *  
1757 *   std::vector<Point<dim>> support_points;
1758 *   std::vector<types::global_dof_index> support_points_idx;
1759 *  
1760 *   const bool has_ghost_elements = signed_distance.has_ghost_elements();
1761 *  
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);
1765 *  
1766 *   if (has_ghost_elements == false)
1767 *   signed_distance.update_ghost_values();
1768 *  
1769 *   for (const auto &cell :
1770 *   tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
1771 *   {
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);
1776 *  
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);
1781 *  
1782 *   for (const auto q : req_values.quadrature_point_indices())
1783 *   if (std::abs(temp_distance[q]) < narrow_band_threshold)
1784 *   {
1785 *   const auto idx = local_dof_indices[q];
1786 *  
1787 *   if (locally_owned_dofs_req.is_element(idx) == false ||
1788 *   flags[locally_owned_dofs_req.index_within_set(idx)])
1789 *   continue;
1790 *  
1791 *   flags[locally_owned_dofs_req.index_within_set(idx)] = true;
1792 *  
1793 *   support_points_idx.emplace_back(idx);
1794 *   support_points.emplace_back(req_values.quadrature_point(q));
1795 *   }
1796 *   }
1797 *  
1798 *   if (has_ghost_elements == false)
1799 *   signed_distance.zero_out_ghost_values();
1800 *  
1801 *   return {support_points, support_points_idx};
1802 *   }
1803 *  
1804 *   template <int spacedim>
1805 *   std::vector<Point<spacedim>> convert(
1806 *   const LinearAlgebra::distributed::Vector<double> &support_points_unrolled)
1807 *   {
1808 *   const unsigned int n_points =
1809 *   support_points_unrolled.locally_owned_size() / spacedim;
1810 *  
1811 *   std::vector<Point<spacedim>> points(n_points);
1812 *  
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);
1816 *  
1817 *   return points;
1818 *   }
1819 *  
1820 *   } // namespace Step87
1821 *  
1822 * @endcode
1823 *
1824 *
1825 * <a name="step_87-Driver"></a>
1826 * <h3>Driver</h3>
1827 *
1828
1829 *
1830 * Finally, the driver of the program executes the four mini examples.
1831 *
1832 * @code
1833 *   int main(int argc, char **argv)
1834 *   {
1835 *   using namespace dealii;
1836 *   Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
1837 *   std::cout.precision(5);
1838 *  
1839 *   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
1840 *   Step87::example_0(); // only run on root process
1841 *  
1842 *   Step87::example_1();
1843 *   Step87::example_2();
1844 *   Step87::example_3();
1845 *   }
1846 * @endcode
1847<a name="step_87-Results"></a><h1>Results</h1>
1848
1849
1850<a name="step_87-Miniexample0"></a><h3>Mini example 0</h3>
1851
1852
1853We present a part of the terminal output. It shows, for each point, the
1854determined cell and reference position. Also, one can see that
1855the values evaluated with FEValues, FEPointEvaluation, and
1856VectorTools::point_values() are identical, as expected.
1857
1858@verbatim
1859Running: example 0
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)
1863 - Values at point:
1864 - 0.25002 (w. FEValues)
1865 - 0.25002 (w. FEPointEvaluation)
1866 - 0.25002 (w. VectorTools::point_value())
1867
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)
1871 - Values at point:
1872 - 0.20003 (w. FEValues)
1873 - 0.20003 (w. FEPointEvaluation)
1874 - 0.20003 (w. VectorTools::point_value())
1875
1876...
1877
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)
1881 - Values at point:
1882 - 0.25002 (w. FEValues)
1883 - 0.25002 (w. FEPointEvaluation)
1884 - 0.25002 (w. VectorTools::point_value())
1885
1886 - writing csv file
1887@endverbatim
1888
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.
1894
1895@verbatim
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
1917@endverbatim
1918
1919<a name="step_87-Miniexample1"></a><h3>Mini example 1</h3>
1920
1921
1922We show the terminal output.
1923
1924@verbatim
1925Running: example 1
1926 - writing csv file
1927@endverbatim
1928
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.
1933
1934@verbatim
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
1956@endverbatim
1957
1958
1959<a name="step_87-Miniexample2"></a><h3>Mini example 2</h3>
1960
1961
1962We show the terminal output.
1963@verbatim
1964Running: example 2
1965 - create system
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
1973 - output results
1974@endverbatim
1975
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
1980be corrected.
1981It can be seen that the numbers of points to be processed decrease
1982from iteration to iteration.
1983<table align="center" class="doxtable">
1984 <tr>
1985 <td>
1986 @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_p_0.png
1987 </td>
1988 <td>
1989 @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_p_1.png
1990 </td>
1991 <td>
1992 @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_p_2.png
1993 </td>
1994 </tr>
1995</table>
1996
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.
2004
2005<table align="center" class="doxtable">
2006 <tr>
2007 <td>
2008 @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_res_0.png
2009 </td>
2010 <td>
2011 @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_res_1.png
2012 </td>
2013 </tr>
2014</table>
2015
2016<a name="step_87-Miniexample3"></a><h3>Mini example 3</h3>
2017
2018
2019We show a shortened version of the terminal output.
2020
2021@verbatim
2022Running: example 3
2023 - creating background mesh
2024 - creating immersed mesh
2025time: 0
2026 - compute to be tested values (immersed mesh)
2027 - test values (background mesh)
2028 - write data (background mesh)
2029 - write mesh (immersed mesh)
2030
2031time: 0.01
2032 - move support points (immersed mesh)
2033 - compute to be tested values (immersed mesh)
2034 - test values (background mesh)
2035
2036time: 0.02
2037 - move support points (immersed mesh)
2038 - compute to be tested values (immersed mesh)
2039 - test values (background mesh)
2040
2041...
2042
2043time: 2
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)
2049@endverbatim
2050
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
2056the immersed mesh.
2057
2058<table align="center" class="doxtable">
2059 <tr>
2060 <td>
2061 @image html https://www.dealii.org/images/steps/developer/step_87_ex_3_force.0000.png
2062 </td>
2063 <td>
2064 @image html https://www.dealii.org/images/steps/developer/step_87_ex_3_force.0010.png
2065 </td>
2066 <td>
2067 @image html https://www.dealii.org/images/steps/developer/step_87_ex_3_force.0020.png
2068 </td>
2069 </tr>
2070</table>
2071
2072<a name="step_87-Possibilitiesforextension"></a><h3>Possibilities for extension</h3>
2073
2074
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
2080Utilities::MPI::RemotePointEvaluation
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
2089cells.
2090- We concentrated on parallelization aspects in this tutorial. However, we would
2091like to point out the need for fast evaluation on cell level.
2092The task for this in deal.II is FEPointEvaluation. It exploits the structure of
2093@f[
2094\hat{u}(\hat{\boldsymbol{x}}) = \sum_i \hat{N}_i(\hat{\boldsymbol{x}}) \hat{u}_i
2095@f]
2096to derive fast implementations, e.g., for tensor-product elements
2097@f[
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)
2102\hat{u}_{ijk}.
2103@f]
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.
2106- Utilities::MPI::RemotePointEvaluation is used in multiple places in deal.II.
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
2110slices. In addition, MGTwoLevelTransferNonNested allows to prolongate solutions
2111and restrict residuals between two independent meshes. By passing a sequence
2112of such two-level transfer operators to MGTransferMF and, finally, to Multigrid,
2113non-nested multigrid can be computed.
2114- Utilities::MPI::RemotePointEvaluation can be used to couple non-matching
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().
2121- Alternatively to the coupling via Utilities::MPI::RemotePointEvaluation,
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.
2125 *
2126 *
2127<a name="step_87-PlainProg"></a>
2128<h1> The plain program</h1>
2129@include "step-87.cc"
2130*/
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
void reinit(const std::vector< Point< spacedim > > &points, const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping)
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
#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
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={}, 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