Reference documentation for deal.II version GIT relicensing-487-ge9eb5ab491 2024-04-25 07:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-87.h
Go to the documentation of this file.
1);
1004 *  
1005 * @endcode
1006 *
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 =
1059 *   Utilities::MPI::reduce<std::vector<Point<dim>>>(
1060 *   unmatched_points, MPI_COMM_WORLD, [](const auto &a, const auto &b) {
1061 *   auto result = a;
1062 *   result.insert(result.end(), b.begin(), b.end());
1063 *   return result;
1064 *   });
1065 *  
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 *   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 *   QGauss<dim - 1> quadrature_immersed(degree + 1);
1308 *  
1309 *   FE_Q<dim - 1, dim> fe_scalar_immersed(degree);
1310 *   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 =
1390 *   VectorTools::point_values<dim>(rpe,
1391 *   dof_handler_background,
1392 *   velocity);
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_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:550
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
Triangulation< dim, spacedim > & get_triangulation()
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
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(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:442
@ 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)
spacedim & mapping
spacedim const Point< spacedim > & p
Definition grid_tools.h:990
const std::vector< bool > & used
const Triangulation< dim, spacedim > & tria
spacedim & mesh
Definition grid_tools.h:989
if(marked_vertices.size() !=0) for(auto it
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
@ 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)
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={})
spacedim const DoFHandler< dim, spacedim > const FullMatrix< double > & transfer
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
const Iterator const std_cxx20::type_identity_t< Iterator > & end
Definition parallel.h:610
const InputIterator OutputIterator const Function & function
Definition parallel.h:168
STL namespace.
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition types.h:32