step-87.h

Go to the documentation of this file.

);
 *   
 * @endcode
 * 
 * For the iterative solution procedure of the closest-point projection,
 * the maximum number of iterations and the tolerance for the maximum
 * absolute acceptable change in the distance in one iteration are set.
 * 
 * @code
 *       pcout << "  - determine closest point iteratively" << std::endl;
 *       constexpr int    max_iter     = 30;
 *       constexpr double tol_distance = 1e-6;
 *   
 * @endcode
 * 
 * Now, we are ready to perform the algorithm by setting an initial guess
 * for the projection points simply corresponding to the collected support
 * points. We collect the global indices of the support points and the
 * total number of points that need to be processed and do not
 * fulfill the required tolerance. Those will be gradually reduced
 * upon the iterative process.
 * 
 * @code
 *       std::vector<Point<dim>> closest_points = support_points; // initial guess
 *   
 *       std::vector<unsigned int> unmatched_points_idx(closest_points.size());
 *       std::iota(unmatched_points_idx.begin(), unmatched_points_idx.end(), 0);
 *   
 *       int n_unmatched_points =
 *         Utilities::MPI::sum(unmatched_points_idx.size(), MPI_COMM_WORLD);
 *   
 * @endcode
 * 
 * Now, we create a Utilities::MPI::RemotePointEvaluation cache object and
 * start the loop for the fix-point iteration. We update the list of points
 * that still need to be processed and subsequently pass this information
 * to the Utilities::MPI::RemotePointEvaluation object. For the sake of
 * illustration, we export the coordinates of the points to be updated for
 * each iteration to a CSV file. Next, we can evaluate the signed distance
 * function and the gradient at those points to update the current solution
 * for the closest points. We perform the update only if the signed
 * distance of the closest point is not already within the tolerance
 * and register those points that still need to be processed.
 * 
 * @code
 *       Utilities::MPI::RemotePointEvaluation<dim, dim> rpe;
 *   
 *       for (int it = 0; it < max_iter && n_unmatched_points > 0; ++it)
 *         {
 *           pcout << "    - iteration " << it << ": " << n_unmatched_points;
 *   
 *           std::vector<Point<dim>> unmatched_points(unmatched_points_idx.size());
 *           for (unsigned int i = 0; i < unmatched_points_idx.size(); ++i)
 *             unmatched_points[i] = closest_points[unmatched_points_idx[i]];
 *   
 *           const auto all_unmatched_points =
 *             Utilities::MPI::reduce<std::vector<Point<dim>>>(
 *               unmatched_points, MPI_COMM_WORLD, [](const auto &a, const auto &b) {
 *                 auto result = a;
 *                 result.insert(result.end(), b.begin(), b.end());
 *                 return result;
 *               });
 *   
 *           if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
 *             {
 *               std::ofstream file("example_2_" + std::to_string(it) + ".csv");
 *               for (const auto &p : all_unmatched_points)
 *                 file << p << std::endl;
 *               file.close();
 *             }
 *   
 *           rpe.reinit(unmatched_points, tria, mapping);
 *   
 *           AssertThrow(rpe.all_points_found(),
 *                       ExcMessage("Processed point is outside domain."));
 *   
 *           const auto eval_values =
 *             VectorTools::point_values<1>(rpe, dof_handler, signed_distance);
 *   
 *           const auto eval_gradient =
 *             VectorTools::point_gradients<1>(rpe, dof_handler, signed_distance);
 *   
 *           std::vector<unsigned int> unmatched_points_idx_next;
 *   
 *           for (unsigned int i = 0; i < unmatched_points_idx.size(); ++i)
 *             if (std::abs(eval_values[i]) > tol_distance)
 *               {
 *                 closest_points[unmatched_points_idx[i]] -=
 *                   eval_values[i] * eval_gradient[i];
 *   
 *                 unmatched_points_idx_next.emplace_back(unmatched_points_idx[i]);
 *               }
 *   
 *           unmatched_points_idx.swap(unmatched_points_idx_next);
 *   
 *           n_unmatched_points =
 *             Utilities::MPI::sum(unmatched_points_idx.size(), MPI_COMM_WORLD);
 *   
 *           pcout << " -> " << n_unmatched_points << std::endl;
 *         }
 *   
 * @endcode
 * 
 * We print a warning message if we exceed the maximum number of allowed
 * iterations and if there are still projection points with a distance
 * value exceeding the tolerance.
 * 
 * @code
 *       if (n_unmatched_points > 0)
 *         pcout << "WARNING: The tolerance of " << n_unmatched_points
 *               << " points is not yet attained." << std::endl;
 *   
 * @endcode
 * 
 * As a result, we obtain a list of support points and corresponding
 * closest points at the zero-isosurface level set. This information
 * can be used to update the signed distance function, i.e., the
 * reinitialization the values of the level-set function to maintain
 * the signed distance property @cite henri2022geometrical.
 * 
 * @code
 *       pcout << "  - determine distance in narrow band" << std::endl;
 *       LinearAlgebra::distributed::Vector<double> solution_distance;
 *       solution_distance.reinit(solution);
 *   
 *       for (unsigned int i = 0; i < closest_points.size(); ++i)
 *         solution_distance[support_points_idx[i]] =
 *           support_points[i].distance(closest_points[i]);
 *   
 * @endcode
 * 
 * In addition, we use the information of the closest point to
 * extrapolate values from the interface, i.e., the zero-level
 * set isosurface, to the support points within the narrow band.
 * This might be helpful to improve accuracy, e.g., for
 * diffuse interface fluxes where certain quantities are only
 * accurately determined at the interface (e.g. curvature
 * for surface tension @cite coquerelle2016fourth).
 * 
 * @code
 *       pcout << "  - perform extrapolation in narrow band" << std::endl;
 *       rpe.reinit(closest_points, tria, mapping);
 *       solution.update_ghost_values();
 *       const auto vals = VectorTools::point_values<1>(rpe, dof_handler, solution);
 *   
 *       LinearAlgebra::distributed::Vector<double> solution_extrapolated;
 *       solution_extrapolated.reinit(solution);
 *   
 *       for (unsigned int i = 0; i < closest_points.size(); ++i)
 *         solution_extrapolated[support_points_idx[i]] = vals[i];
 *   
 * @endcode
 * 
 * Finally, we output the results to a VTU file.
 * 
 * @code
 *       pcout << "  - output results" << std::endl;
 *       DataOut<dim> data_out;
 *       data_out.add_data_vector(dof_handler, signed_distance, "signed_distance");
 *       data_out.add_data_vector(dof_handler, solution, "solution");
 *       data_out.add_data_vector(dof_handler,
 *                                solution_distance,
 *                                "solution_distance");
 *       data_out.add_data_vector(dof_handler,
 *                                solution_extrapolated,
 *                                "solution_extrapolated");
 *       data_out.build_patches(mapping);
 *       data_out.write_vtu_in_parallel("example_2.vtu", MPI_COMM_WORLD);
 *   
 *       pcout << std::endl;
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="step_87-Miniexample3Sharpinterfacemethodontheexampleofsurfacetensionforfronttracking"></a> 
 * <h3>Mini example 3: Sharp interface method on the example of surface tension for front tracking</h3>
 *   
 
 * 
 * The final mini example presents a basic implementation of
 * front tracking @cite peskin1977numerical, @cite unverdi1992front
 * of a surface mesh @f$\mathbb{T}_\Gamma@f$ immersed
 * in a Eulerian background fluid mesh @f$\mathbb{T}_\Omega@f$.
 *   
 
 * 
 * We assume that the immersed surface is transported according to a
 * prescribed velocity field from the background mesh. Subsequently,
 * we perform a sharp computation of the surface-tension force:
 * @f[
 * (\boldsymbol v_i (\boldsymbol{x}), \boldsymbol F_S
 * (\boldsymbol{x}))_{\Omega}
 * =
 * \left( \boldsymbol{v}_i (\boldsymbol{x}), \sigma (\boldsymbol{x}) \kappa
 * (\boldsymbol{x}) \boldsymbol{n} (\boldsymbol{x}) \right)_\Gamma \approx
 * \sum_{q\in\mathbb{T}_\Gamma} \boldsymbol{v}_i^T (\boldsymbol{x}_q)
 * \sigma (\boldsymbol{x}_q) \kappa (\boldsymbol{x}_q) \boldsymbol{n}
 * (\boldsymbol{x}_q) |J(\boldsymbol{x}_q)| w(\boldsymbol{x}_q) \quad \forall
 * i\in\mathbb{T}_\Omega
 * .
 * @f]
 * We decompose this operation into two steps. In the first step, we evaluate
 * the force contributions @f$\sigma (\boldsymbol{x}_q) \kappa
 * (\boldsymbol{x}_q) \boldsymbol{n}
 * (\boldsymbol{x}_q)@f$ at the quadrature points defined on the immersed mesh
 * and multiply them with the mapped quadrature weight @f$|J(\boldsymbol{x}_q)|
 * w_q@f$:
 * @f[
 * \boldsymbol{F}_S (\boldsymbol{x}_q) \gets \sigma (\boldsymbol{x}_q) \kappa
 * (\boldsymbol{x}_q) \boldsymbol{n} (\boldsymbol{x}_q) |J(\boldsymbol{x}_q)|
 * w_q \quad \forall q\in\mathbb{T}_\Gamma.
 * @f]
 * In the second step, we compute the discretized weak form by multiplying
 * with test functions on the background mesh:
 * @f[
 * (\boldsymbol v_i (\boldsymbol{x}), \boldsymbol F_S
 * (\boldsymbol{x}))_{\Omega} \gets \sum_q \boldsymbol{v}_i^T
 * (\boldsymbol{x}_q) \boldsymbol{F}_S
 * (\boldsymbol{x}_q)
 * \quad \forall i\in\mathbb{T}_\Omega
 * .
 * @f]
 * Obviously, we need to communicate between the two steps. The second step
 * can be handled completely by Utilities::MPI::RemotePointEvaluation, which
 * provides the function
 * Utilities::MPI::RemotePointEvaluation::process_and_evaluate() for this
 * purpose.
 *   
 
 * 
 * We start with setting the parameters consisting of the polynomial degree of
 * the shape functions, the dimension of the background mesh, the time-step
 * size to be considered for transporting the surface mesh and the number of
 * time steps.
 * 
 
 * 
 * 
 * @code
 *     void example_3()
 *     {
 *       constexpr unsigned int degree       = 3;
 *       constexpr unsigned int dim          = 2;
 *       const double           dt           = 0.01;
 *       const unsigned int     n_time_steps = 200;
 *   
 * @endcode
 * 
 * This program is intended to be executed in 2D or 3D.
 * 
 * @code
 *       static_assert(dim == 2 || dim == 3, "Only implemented for 2D or 3D.");
 *   
 *       ConditionalOStream pcout(std::cout,
 *                                Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) ==
 *                                  0);
 *   
 *       pcout << "Running: example 3" << std::endl;
 *   
 * @endcode
 * 
 * Next, we create the standard objects necessary for the finite element
 * representation of the background mesh
 * 
 * @code
 *       pcout << "  - creating background mesh" << std::endl;
 *       DistributedTriangulation<dim> tria_background(MPI_COMM_WORLD);
 *       GridGenerator::hyper_cube(tria_background);
 *       tria_background.refine_global(5);
 *   
 *       const MappingQ1<dim> mapping_background;
 *       const FESystem<dim>  fe_background(FE_Q<dim>(degree), dim);
 *       DoFHandler<dim>      dof_handler_background(tria_background);
 *       dof_handler_background.distribute_dofs(fe_background);
 *   
 * @endcode
 * 
 * and, similarly, for the immersed surface mesh.
 * We use a sphere with radius @f$r=0.75@f$ which is
 * placed in the center of the top half of the cubic background domain.
 * 
 * @code
 *       pcout << "  - creating immersed mesh" << std::endl;
 *       const Point<dim> center((dim == 2) ? Point<dim>(0.5, 0.75) :
 *                                            Point<dim>(0.5, 0.75, 0.5));
 *       const double     radius = 0.15;
 *   
 *       DistributedTriangulation<dim - 1, dim> tria_immersed(MPI_COMM_WORLD);
 *       GridGenerator::hyper_sphere(tria_immersed, center, radius);
 *       tria_immersed.refine_global(4);
 *   
 * @endcode
 * 
 * Two different mappings are considered for the immersed
 * surface mesh: one valid for the initial configuration and one
 * that is updated in every time step according to the nodal
 * displacements. Two types of finite elements are used to
 * represent scalar and vector-valued DoF values.
 * 
 * @code
 *       const MappingQ<dim - 1, dim> mapping_immersed_base(3);
 *       MappingQCache<dim - 1, dim>  mapping_immersed(3);
 *       mapping_immersed.initialize(mapping_immersed_base, tria_immersed);
 *       const QGauss<dim - 1> quadrature_immersed(degree + 1);
 *   
 *       const FE_Q<dim - 1, dim>     fe_scalar_immersed(degree);
 *       const FESystem<dim - 1, dim> fe_immersed(fe_scalar_immersed, dim);
 *       DoFHandler<dim - 1, dim>     dof_handler_immersed(tria_immersed);
 *       dof_handler_immersed.distribute_dofs(fe_immersed);
 *   
 * @endcode
 * 
 * We renumber the DoFs related to the vector-valued problem to
 * simplify access to the individual components.
 * 
 * @code
 *       DoFRenumbering::support_point_wise(dof_handler_immersed);
 *   
 * @endcode
 * 
 * We fill a DoF vector on the background mesh with an analytical
 * velocity field considering the Rayleigh-Kothe vortex flow and
 * initialize a DoF vector for the weak form of the surface-tension force.
 * 
 * @code
 *       LinearAlgebra::distributed::Vector<double> velocity;
 *       velocity.reinit(dof_handler_background.locally_owned_dofs(),
 *                       DoFTools::extract_locally_active_dofs(
 *                         dof_handler_background),
 *                       MPI_COMM_WORLD);
 *       Functions::RayleighKotheVortex<dim> vortex(2);
 *   
 *       LinearAlgebra::distributed::Vector<double> force_vector(
 *         dof_handler_background.locally_owned_dofs(),
 *         DoFTools::extract_locally_active_dofs(dof_handler_background),
 *         MPI_COMM_WORLD);
 *   
 * @endcode
 * 
 * Next, we collect the real positions @f$\boldsymbol{x}_q@f$ of the quadrature
 * points of the surface mesh in a vector.
 * 
 * @code
 *       LinearAlgebra::distributed::Vector<double> immersed_support_points;
 *       collect_support_points(mapping_immersed,
 *                              dof_handler_immersed,
 *                              immersed_support_points);
 *   
 * @endcode
 * 
 * We initialize a Utilities::MPI::RemotePointEvaluation object and start
 * the time loop. For any other step than the initial one, we first move the
 * support points of the surface mesh according to the fluid velocity of the
 * background mesh. Thereto, we first update the time of the velocity
 * function. Then, we update the internal data structures of the
 * Utilities::MPI::RemotePointEvaluation object with the collected support
 * points of the immersed mesh. We throw an exception if one of the points
 * cannot be found within the domain of the background mesh. Next, we
 * evaluate the velocity at the surface-mesh support points and compute the
 * resulting update of the coordinates. Finally, we update the mapping of
 * the immersed surface mesh to the current position.
 * 
 * @code
 *       Utilities::MPI::RemotePointEvaluation<dim> rpe;
 *       double                                     time = 0.0;
 *       for (unsigned int it = 0; it <= n_time_steps; ++it, time += dt)
 *         {
 *           pcout << "time: " << time << std::endl;
 *           if (it > 0)
 *             {
 *               pcout << "  - move support points (immersed mesh)" << std::endl;
 *               vortex.set_time(time);
 *               VectorTools::interpolate(mapping_background,
 *                                        dof_handler_background,
 *                                        vortex,
 *                                        velocity);
 *               rpe.reinit(convert<dim>(immersed_support_points),
 *                          tria_background,
 *                          mapping_background);
 *   
 *               AssertThrow(rpe.all_points_found(),
 *                           ExcMessage(
 *                             "Immersed domain leaves background domain!"));
 *   
 *               velocity.update_ghost_values();
 *               const auto immersed_velocity =
 *                 VectorTools::point_values<dim>(rpe,
 *                                                dof_handler_background,
 *                                                velocity);
 *               velocity.zero_out_ghost_values();
 *   
 *               for (unsigned int i = 0, c = 0;
 *                    i < immersed_support_points.locally_owned_size() / dim;
 *                    ++i)
 *                 for (unsigned int d = 0; d < dim; ++d, ++c)
 *                   immersed_support_points.local_element(c) +=
 *                     immersed_velocity[i][d] * dt;
 *   
 *               mapping_immersed.initialize(mapping_immersed_base,
 *                                           dof_handler_immersed,
 *                                           immersed_support_points,
 *                                           false);
 *             }
 *   
 * @endcode
 * 
 * Next, we loop over all locally owned cells of the immersed mesh and
 * its quadrature points to compute the value for the local surface
 * tension force contribution @f$\boldsymbol{F}_S(\boldsymbol{x}_q)@f$. We
 * store the real coordinates of the quadrature points and the
 * corresponding force contributions in two individual vectors. For
 * computation of the latter, the normal vector
 * @f$\boldsymbol{n}(\boldsymbol{x}_q)@f$ can be directly extracted from the
 * surface mesh via FEValues and, for the curvature, we use the
 * following approximation:
 * @f[
 * \kappa(\boldsymbol{x}_q)
 * =
 * \nabla \cdot \boldsymbol{n}(\boldsymbol{x}_q)
 * =
 * \text{tr}\left({\nabla \boldsymbol{n}(\boldsymbol{x}_q)}\right)
 * \approx
 * \text{tr}\left({\nabla \sum_i \boldsymbol{N}_i (\boldsymbol{x}_q)
 * \boldsymbol n_i}\right)
 * =
 * \sum_i\text{tr}\left({\nabla \boldsymbol{N}_i (\boldsymbol{x}_q)
 * \boldsymbol n_i}\right)
 * \;\text{with}\; i\in[0,n_{\text{dofs\_per\_cell}}),
 * @f]
 * which we can apply since the immersed mesh is consistently
 * orientated. The surface tension coefficient is set to 1 for the
 * sake of demonstration.
 * 
 * @code
 *           pcout << "  - compute to be tested values (immersed mesh)" << std::endl;
 *           using value_type = Tensor<1, dim, double>;
 *   
 *           std::vector<Point<dim>> integration_points;
 *           std::vector<value_type> integration_values;
 *   
 *           FEValues<dim - 1, dim> fe_values(mapping_immersed,
 *                                            fe_immersed,
 *                                            quadrature_immersed,
 *                                            update_JxW_values | update_gradients |
 *                                              update_normal_vectors |
 *                                              update_quadrature_points);
 *   
 *           FEValues<dim - 1, dim> fe_values_co(
 *             mapping_immersed,
 *             fe_scalar_immersed,
 *             fe_scalar_immersed.get_unit_support_points(),
 *             update_JxW_values | update_normal_vectors);
 *   
 *           std::vector<unsigned int> component_to_system_index(
 *             fe_immersed.n_dofs_per_cell());
 *   
 *           for (unsigned int i = 0, c = 0;
 *                i < fe_scalar_immersed.n_dofs_per_cell();
 *                ++i)
 *             for (unsigned int d = 0; d < dim; ++d, ++c)
 *               component_to_system_index[c] =
 *                 fe_immersed.component_to_system_index(d, i);
 *   
 *           for (const auto &cell : tria_immersed.active_cell_iterators() |
 *                                     IteratorFilters::LocallyOwnedCell())
 *             {
 *               fe_values.reinit(cell);
 *               fe_values_co.reinit(cell);
 *   
 *               for (const auto &q : fe_values.quadrature_point_indices())
 *                 {
 *                   integration_points.emplace_back(fe_values.quadrature_point(q));
 *   
 *                   const auto sigma = 1.0; // surface tension coefficient
 *   
 *                   const auto normal    = fe_values.normal_vector(q);
 *                   double     curvature = 0;
 *                   for (unsigned int i = 0, c = 0;
 *                        i < fe_scalar_immersed.n_dofs_per_cell();
 *                        ++i)
 *                     for (unsigned int d = 0; d < dim; ++d, ++c)
 *                       curvature += fe_values.shape_grad_component(
 *                                      component_to_system_index[c], q, d)[d] *
 *                                    fe_values_co.normal_vector(i)[d];
 *   
 *                   const auto FxJxW =
 *                     sigma * curvature * normal * fe_values.JxW(q);
 *   
 *                   integration_values.emplace_back(FxJxW);
 *                 }
 *             }
 *   
 * @endcode
 * 
 * Before we evaluate the weak form of the surface-tension force, the
 * communication pattern of Utilities::MPI::RemotePointEvaluation is
 * set up from the quadrature points of the immersed mesh, determining
 * the surrounding cells on the background mesh.
 * 
 * @code
 *           pcout << "  - test values (background mesh)" << std::endl;
 *   
 *           rpe.reinit(integration_points, tria_background, mapping_background);
 *   
 * @endcode
 * 
 * In preparation for utilizing
 * Utilities::MPI::RemotePointEvaluation::process_and_evaluate that
 * performs the
 * multiplication with the test function, we set up a callback function
 * that contains the operation on the intersected cells of the
 * background mesh. Within this function, we initialize a
 * FEPointEvaluation object that allows us to integrate values at
 * arbitrary points within a cell. We loop over the cells that surround
 * quadrature points of the immersed mesh -- provided by the callback
 * function. From the provided CellData object, we retrieve the unit
 * points, i.e., the quadrature points of the immersed mesh that lie
 * within the current cell and a pointer to the stored values on the
 * current cell (local surface-tension force) for convenience. We
 * reinitialize the data structure of FEPointEvaluation on every cell
 * according to the unit points. Next, we loop over the quadrature
 * points attributed to the cell and submit the local surface-tension
 * force to the FEPointEvaluation object. Via
 * FEPointEvaluation::test_and_sum(), the submitted values are
 * multiplied by the values of the test function and a summation over
 * all given points is performed. Subsequently, the contributions are
 * assembled into the global vector containing the weak form of the
 * surface-tension force.
 * 
 * @code
 *           const auto integration_function = [&](const auto &values,
 *                                                 const auto &cell_data) {
 *             FEPointEvaluation<dim, dim> phi_force(mapping_background,
 *                                                   fe_background,
 *                                                   update_values);
 *   
 *             std::vector<double>                  local_values;
 *             std::vector<types::global_dof_index> local_dof_indices;
 *   
 *             for (const auto cell : cell_data.cell_indices())
 *               {
 *                 const auto cell_dofs =
 *                   cell_data.get_active_cell_iterator(cell)
 *                     ->as_dof_handler_iterator(dof_handler_background);
 *   
 *                 const auto unit_points = cell_data.get_unit_points(cell);
 *                 const auto FxJxW       = cell_data.get_data_view(cell, values);
 *   
 *                 phi_force.reinit(cell_dofs, unit_points);
 *   
 *                 for (const auto q : phi_force.quadrature_point_indices())
 *                   phi_force.submit_value(FxJxW[q], q);
 *   
 *                 local_values.resize(cell_dofs->get_fe().n_dofs_per_cell());
 *                 phi_force.test_and_sum(local_values, EvaluationFlags::values);
 *   
 *                 local_dof_indices.resize(cell_dofs->get_fe().n_dofs_per_cell());
 *                 cell_dofs->get_dof_indices(local_dof_indices);
 *                 AffineConstraints<double>().distribute_local_to_global(
 *                   local_values, local_dof_indices, force_vector);
 *               }
 *           };
 *   
 * @endcode
 * 
 * The callback function is passed together with the vector holding the
 * surface-tension force contribution at each quadrature point of the
 * immersed mesh to
 * Utilities::MPI::RemotePointEvaluation::process_and_evaluate. The only
 * missing step is to compress the distributed force vector.
 * 
 * @code
 *           rpe.process_and_evaluate<value_type>(integration_values,
 *                                                integration_function);
 *           force_vector.compress(VectorOperation::add);
 *   
 * @endcode
 * 
 * After every tenth step or at the beginning/end of the time loop, we
 * output the force vector and the velocity of the background mesh to
 * a VTU file. In addition, we also export the geometry of the
 * (deformed) immersed surface mesh to a separate VTU file.
 * 
 * @code
 *           if (it % 10 == 0 || it == n_time_steps)
 *             {
 *               std::vector<
 *                 DataComponentInterpretation::DataComponentInterpretation>
 *                 vector_component_interpretation(
 *                   dim, DataComponentInterpretation::component_is_part_of_vector);
 *               pcout << "  - write data (background mesh)" << std::endl;
 *               DataOut<dim>          data_out_background;
 *               DataOutBase::VtkFlags flags_background;
 *               flags_background.write_higher_order_cells = true;
 *               data_out_background.set_flags(flags_background);
 *               data_out_background.add_data_vector(
 *                 dof_handler_background,
 *                 force_vector,
 *                 std::vector<std::string>(dim, "force"),
 *                 vector_component_interpretation);
 *               data_out_background.add_data_vector(
 *                 dof_handler_background,
 *                 velocity,
 *                 std::vector<std::string>(dim, "velocity"),
 *                 vector_component_interpretation);
 *               data_out_background.build_patches(mapping_background, 3);
 *               data_out_background.write_vtu_in_parallel("example_3_background_" +
 *                                                           std::to_string(it) +
 *                                                           ".vtu",
 *                                                         MPI_COMM_WORLD);
 *   
 *               pcout << "  - write mesh (immersed mesh)" << std::endl;
 *               DataOut<dim - 1, dim> data_out_immersed;
 *               data_out_immersed.attach_triangulation(tria_immersed);
 *               data_out_immersed.build_patches(mapping_immersed, 3);
 *               data_out_immersed.write_vtu_in_parallel("example_3_immersed_" +
 *                                                         std::to_string(it) +
 *                                                         ".vtu",
 *                                                       MPI_COMM_WORLD);
 *             }
 *           pcout << std::endl;
 *         }
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="step_87-Utilityfunctionsdefinition"></a> 
 * <h3>Utility functions (definition)</h3>
 * 
 * @code
 *     template <int spacedim>
 *     std::vector<Point<spacedim>>
 *     create_points_along_line(const Point<spacedim> &p0,
 *                              const Point<spacedim> &p1,
 *                              const unsigned int     n_subdivisions)
 *     {
 *       Assert(n_subdivisions >= 1, ExcInternalError());
 *   
 *       std::vector<Point<spacedim>> points;
 *       points.reserve(n_subdivisions + 1);
 *   
 *       points.emplace_back(p0);
 *       for (unsigned int i = 1; i < n_subdivisions; ++i)
 *         points.emplace_back(p0 + (p1 - p0) * static_cast<double>(i) /
 *                                    static_cast<double>(n_subdivisions));
 *       points.emplace_back(p1);
 *   
 *       return points;
 *     }
 *   
 *     template <int spacedim, typename T0, typename T1>
 *     void print_along_line(const std::string                  &file_name,
 *                           const std::vector<Point<spacedim>> &points,
 *                           const std::vector<T0>              &values_0,
 *                           const std::vector<T1>              &values_1)
 *     {
 *       AssertThrow(points.size() == values_0.size() &&
 *                     (values_1.size() == points.size() || values_1.empty()),
 *                   ExcMessage("The provided vectors must have the same length."));
 *   
 *       std::ofstream file(file_name);
 *   
 *       for (unsigned int i = 0; i < points.size(); ++i)
 *         {
 *           file << std::fixed << std::right << std::setw(5) << std::setprecision(3)
 *                << points[0].distance(points[i]);
 *   
 *           for (unsigned int d = 0; d < spacedim; ++d)
 *             file << std::fixed << std::right << std::setw(10)
 *                  << std::setprecision(3) << points[i][d];
 *   
 *           file << std::fixed << std::right << std::setw(10)
 *                << std::setprecision(3) << values_0[i];
 *   
 *           if (!values_1.empty())
 *             file << std::fixed << std::right << std::setw(10)
 *                  << std::setprecision(3) << values_1[i];
 *           file << std::endl;
 *         }
 *     }
 *   
 *     template <int dim, int spacedim>
 *     void collect_support_points(
 *       const Mapping<dim, spacedim>               &mapping,
 *       const DoFHandler<dim, spacedim>            &dof_handler,
 *       LinearAlgebra::distributed::Vector<double> &support_points)
 *     {
 *       support_points.reinit(dof_handler.locally_owned_dofs(),
 *                             DoFTools::extract_locally_active_dofs(dof_handler),
 *                             dof_handler.get_mpi_communicator());
 *   
 *       const auto &fe = dof_handler.get_fe();
 *   
 *       FEValues<dim, spacedim> fe_values(
 *         mapping,
 *         fe,
 *         fe.base_element(0).get_unit_support_points(),
 *         update_quadrature_points);
 *   
 *       std::vector<types::global_dof_index> local_dof_indices(
 *         fe.n_dofs_per_cell());
 *   
 *       std::vector<unsigned int> component_to_system_index(
 *         fe_values.n_quadrature_points * spacedim);
 *   
 *       for (unsigned int q = 0, c = 0; q < fe_values.n_quadrature_points; ++q)
 *         for (unsigned int d = 0; d < spacedim; ++d, ++c)
 *           component_to_system_index[c] = fe.component_to_system_index(d, q);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators() |
 *                                 IteratorFilters::LocallyOwnedCell())
 *         {
 *           fe_values.reinit(cell);
 *           cell->get_dof_indices(local_dof_indices);
 *   
 *           for (unsigned int q = 0, c = 0; q < fe_values.n_quadrature_points; ++q)
 *             for (unsigned int d = 0; d < spacedim; ++d, ++c)
 *               support_points[local_dof_indices[component_to_system_index[c]]] =
 *                 fe_values.quadrature_point(q)[d];
 *         }
 *     }
 *   
 *     template <int dim, int spacedim, typename T>
 *     std::tuple<std::vector<Point<spacedim>>, std::vector<types::global_dof_index>>
 *     collect_support_points_with_narrow_band(
 *       const Mapping<dim, spacedim>                &mapping,
 *       const DoFHandler<dim, spacedim>             &dof_handler_signed_distance,
 *       const LinearAlgebra::distributed::Vector<T> &signed_distance,
 *       const DoFHandler<dim, spacedim>             &dof_handler_support_points,
 *       const double                                 narrow_band_threshold)
 *     {
 *       AssertThrow(narrow_band_threshold >= 0,
 *                   ExcMessage("The narrow band threshold"
 *                              " must be larger than or equal to 0."));
 *       const auto &tria = dof_handler_signed_distance.get_triangulation();
 *       const Quadrature<dim> quad(dof_handler_support_points.get_fe()
 *                                    .base_element(0)
 *                                    .get_unit_support_points());
 *   
 *       FEValues<dim> distance_values(mapping,
 *                                     dof_handler_signed_distance.get_fe(),
 *                                     quad,
 *                                     update_values);
 *   
 *       FEValues<dim> req_values(mapping,
 *                                dof_handler_support_points.get_fe(),
 *                                quad,
 *                                update_quadrature_points);
 *   
 *       std::vector<T>                       temp_distance(quad.size());
 *       std::vector<types::global_dof_index> local_dof_indices(
 *         dof_handler_support_points.get_fe().n_dofs_per_cell());
 *   
 *       std::vector<Point<dim>>              support_points;
 *       std::vector<types::global_dof_index> support_points_idx;
 *   
 *       const bool has_ghost_elements = signed_distance.has_ghost_elements();
 *   
 *       const auto &locally_owned_dofs_req =
 *         dof_handler_support_points.locally_owned_dofs();
 *       std::vector<bool> flags(locally_owned_dofs_req.n_elements(), false);
 *   
 *       if (has_ghost_elements == false)
 *         signed_distance.update_ghost_values();
 *   
 *       for (const auto &cell :
 *            tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
 *         {
 *           const auto cell_distance =
 *             cell->as_dof_handler_iterator(dof_handler_signed_distance);
 *           distance_values.reinit(cell_distance);
 *           distance_values.get_function_values(signed_distance, temp_distance);
 *   
 *           const auto cell_req =
 *             cell->as_dof_handler_iterator(dof_handler_support_points);
 *           req_values.reinit(cell_req);
 *           cell_req->get_dof_indices(local_dof_indices);
 *   
 *           for (const auto q : req_values.quadrature_point_indices())
 *             if (std::abs(temp_distance[q]) < narrow_band_threshold)
 *               {
 *                 const auto idx = local_dof_indices[q];
 *   
 *                 if (locally_owned_dofs_req.is_element(idx) == false ||
 *                     flags[locally_owned_dofs_req.index_within_set(idx)])
 *                   continue;
 *   
 *                 flags[locally_owned_dofs_req.index_within_set(idx)] = true;
 *   
 *                 support_points_idx.emplace_back(idx);
 *                 support_points.emplace_back(req_values.quadrature_point(q));
 *               }
 *         }
 *   
 *       if (has_ghost_elements == false)
 *         signed_distance.zero_out_ghost_values();
 *   
 *       return {support_points, support_points_idx};
 *     }
 *   
 *     template <int spacedim>
 *     std::vector<Point<spacedim>> convert(
 *       const LinearAlgebra::distributed::Vector<double> &support_points_unrolled)
 *     {
 *       const unsigned int n_points =
 *         support_points_unrolled.locally_owned_size() / spacedim;
 *   
 *       std::vector<Point<spacedim>> points(n_points);
 *   
 *       for (unsigned int i = 0, c = 0; i < n_points; ++i)
 *         for (unsigned int d = 0; d < spacedim; ++d, ++c)
 *           points[i][d] = support_points_unrolled.local_element(c);
 *   
 *       return points;
 *     }
 *   
 *   } // namespace Step87
 *   
 * @endcode
 * 
 * 
 * <a name="step_87-Driver"></a> 
 * <h3>Driver</h3>
 * 
 
 * 
 * Finally, the driver of the program executes the four mini examples.
 * 
 * @code
 *   int main(int argc, char **argv)
 *   {
 *     using namespace dealii;
 *     Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
 *     std::cout.precision(5);
 *   
 *     if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
 *       Step87::example_0(); // only run on root process
 *   
 *     Step87::example_1();
 *     Step87::example_2();
 *     Step87::example_3();
 *   }
 * @endcode
<a name="step_87-Results"></a><h1>Results</h1>
 
 
<a name="step_87-Miniexample0"></a><h3>Mini example 0</h3>
 
 
We present a part of the terminal output. It shows, for each point, the
determined cell and reference position. Also, one can see that
the values evaluated with FEValues, FEPointEvaluation, and
VectorTools::point_values() are identical, as expected.
 
@verbatim
Running: example 0
 - Found point with real coordinates: 0 0.5
   - in cell with vertices: (0 0.4) (0.2 0.4) (0 0.6) (0.2 0.6)
   - with coordinates on the unit cell: (0 0.5)
 - Values at point:
  - 0.25002 (w. FEValues)
  - 0.25002 (w. FEPointEvaluation)
  - 0.25002 (w. VectorTools::point_value())
 
 - Found point with real coordinates: 0.05 0.5
   - in cell with vertices: (0 0.4) (0.2 0.4) (0 0.6) (0.2 0.6)
   - with coordinates on the unit cell: (0.25 0.5)
 - Values at point:
  - 0.20003 (w. FEValues)
  - 0.20003 (w. FEPointEvaluation)
  - 0.20003 (w. VectorTools::point_value())
 
...
 
 - Found point with real coordinates: 1 0.5
   - in cell with vertices: (0.8 0.4) (1 0.4) (0.8 0.6) (1 0.6)
   - with coordinates on the unit cell: (1 0.5)
 - Values at point:
  - 0.25002 (w. FEValues)
  - 0.25002 (w. FEPointEvaluation)
  - 0.25002 (w. VectorTools::point_value())
 
 - writing csv file
@endverbatim
 
The CSV output is as follows and contains, in the
first column, the distances with respect to the first point,
the second and the third column represent the coordinates
of the points and the fourth column the evaluated solution
values at those points.
 
@verbatim
0.000     0.000     0.500     0.250
0.050     0.050     0.500     0.200
0.100     0.100     0.500     0.150
0.150     0.150     0.500     0.100
0.200     0.200     0.500     0.050
0.250     0.250     0.500     0.000
0.300     0.300     0.500    -0.050
0.350     0.350     0.500    -0.100
0.400     0.400     0.500    -0.149
0.450     0.450     0.500    -0.200
0.500     0.500     0.500    -0.222
0.550     0.550     0.500    -0.200
0.600     0.600     0.500    -0.149
0.650     0.650     0.500    -0.100
0.700     0.700     0.500    -0.050
0.750     0.750     0.500     0.000
0.800     0.800     0.500     0.050
0.850     0.850     0.500     0.100
0.900     0.900     0.500     0.150
0.950     0.950     0.500     0.200
1.000     1.000     0.500     0.250
@endverbatim
 
<a name="step_87-Miniexample1"></a><h3>Mini example 1</h3>
 
 
We show the terminal output.
 
@verbatim
Running: example 1
 - writing csv file
@endverbatim
 
The CSV output is as follows and identical to the results
of the serial case presented in mini example 0.
The fifth column shows the
user quantity evaluated additionally in this mini example.
 
@verbatim
0.000     0.000     0.500     0.250     0.000
0.050     0.050     0.500     0.200     0.050
0.100     0.100     0.500     0.150     0.100
0.150     0.150     0.500     0.100     0.150
0.200     0.200     0.500     0.050     0.200
0.250     0.250     0.500     0.000     0.250
0.300     0.300     0.500    -0.050     0.300
0.350     0.350     0.500    -0.100     0.350
0.400     0.400     0.500    -0.149     0.400
0.450     0.450     0.500    -0.200     0.450
0.500     0.500     0.500    -0.222     0.500
0.550     0.550     0.500    -0.200     0.550
0.600     0.600     0.500    -0.149     0.600
0.650     0.650     0.500    -0.100     0.650
0.700     0.700     0.500    -0.050     0.700
0.750     0.750     0.500     0.000     0.750
0.800     0.800     0.500     0.050     0.800
0.850     0.850     0.500     0.100     0.850
0.900     0.900     0.500     0.150     0.900
0.950     0.950     0.500     0.200     0.950
1.000     1.000     0.500     0.250     1.000
@endverbatim
 
 
<a name="step_87-Miniexample2"></a><h3>Mini example 2</h3>
 
 
We show the terminal output.
@verbatim
Running: example 2
  - create system
  - determine narrow band
  - determine closest point iteratively
    - iteration 0: 7076 -> 7076
    - iteration 1: 7076 -> 104
    - iteration 2: 104 -> 0
  - determine distance in narrow band
  - perform extrapolation in narrow band
  - output results
@endverbatim
 
The following three plots, representing the performed iterations of the
closest-point projection, show the current position of the closest
points exceeding the required tolerance of the discrete interface
of the circle and still need to
be corrected.
It can be seen that the numbers of points to be processed decrease
from iteration to iteration.
<table align="center" class="doxtable">
  <tr>
    <td>
        @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_p_0.png
    </td>
    <td>
        @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_p_1.png
    </td>
    <td>
        @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_p_2.png
    </td>
  </tr>
</table>
 
The output visualized in Paraview looks like the following: On the
left, the original distance function is shown as the light gray surface.
In addition, the contour values refer to the distance values determined
from calculation of the distance to the closest points at the interface
in the narrow band. It can be seen that the two functions coincide.
Similarly, on the right, the original solution and the extrapolated
solution from the interface is shown.
 
<table align="center" class="doxtable">
  <tr>
    <td>
        @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_res_0.png
    </td>
    <td>
        @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_res_1.png
    </td>
  </tr>
</table>
 
<a name="step_87-Miniexample3"></a><h3>Mini example 3</h3>
 
 
We show a shortened version of the terminal output.
 
@verbatim
Running: example 3
  - creating background mesh
  - creating immersed mesh
time: 0
  - compute to be tested values (immersed mesh)
  - test values (background mesh)
  - write data (background mesh)
  - write mesh (immersed mesh)
 
time: 0.01
  - move support points (immersed mesh)
  - compute to be tested values (immersed mesh)
  - test values (background mesh)
 
time: 0.02
  - move support points (immersed mesh)
  - compute to be tested values (immersed mesh)
  - test values (background mesh)
 
...
 
time: 2
  - move support points (immersed mesh)
  - compute to be tested values (immersed mesh)
  - test values (background mesh)
  - write data (background mesh)
  - write mesh (immersed mesh)
@endverbatim
 
The output visualized in Paraview looks like the following: The deformation of
the immersed mesh by the reversible vortex flow can be seen. Due to
discretization errors, the shape is not exactly circular at the end, illustrated
in the right figure. The sharp nature of the surface-tension force vector, shown
as vector plots, can be seen by its restriction to cells that are intersected by
the immersed mesh.
 
<table align="center" class="doxtable">
  <tr>
    <td>
        @image html https://www.dealii.org/images/steps/developer/step_87_ex_3_force.0000.png
    </td>
    <td>
        @image html https://www.dealii.org/images/steps/developer/step_87_ex_3_force.0010.png
    </td>
    <td>
        @image html https://www.dealii.org/images/steps/developer/step_87_ex_3_force.0020.png
    </td>
  </tr>
</table>
 
<a name="step_87-Possibilitiesforextension"></a><h3>Possibilities for extension</h3>
 
 
This program highlights some of the main capabilities
of the distributed evaluation routines in deal.II. However, there are many
related topics worth mentioning:
- Performing a distributed search is an expensive step. That is why we suggest
to provide hints to Utilities::MPI::RemotePointEvaluation and to reuse
Utilities::MPI::RemotePointEvaluation
instances in the case that the communication pattern has not changed.
Furthermore, there  are instances where no search is needed and the points are
already sorted into the right cells. This is the case if the points are
generated on the cell level (see @ref step_85 "step-85"; CutFEM) or the points are
automatically sorted into the correct (neighboring) cell (see @ref step_68 "step-68"; PIC with
Particles::ParticleHandler). Having said that, the
Particles::ParticleHandler::insert_global_particles() function uses
the described infrastructure to perform the initial sorting of particles into
cells.
- We concentrated on parallelization aspects in this tutorial. However, we would
like to point out the need for fast evaluation on cell level.
The task for this in deal.II is FEPointEvaluation. It exploits the structure of
@f[
\hat{u}(\hat{\boldsymbol{x}}) = \sum_i \hat{N}_i(\hat{\boldsymbol{x}}) \hat{u}_i
@f]
to derive fast implementations, e.g., for tensor-product elements
@f[
\hat{u}(\hat{x}_0, \hat{x}_1, \hat{x}_2) =
\sum_k \hat{N}^{\text{1D}}_k(\hat{x}_2)
\sum_j \hat{N}^{\text{1D}}_j(\hat{x}_1)
\sum_i \hat{N}^{\text{1D}}_i(\hat{x}_0)
\hat{u}_{ijk}.
@f]
Since only 1D shape functions are queried and are re-used in re-occurring terms,
this formulation is faster than without exploitation of the structure.
- Utilities::MPI::RemotePointEvaluation is used in multiple places in deal.II.
The class DataOutResample allows to output results on a different mesh than
the computational mesh. This is useful if one wants to output the results
on a coarser mesh or one does not want to output 3D results but instead 2D
slices. In addition, MGTwoLevelTransferNonNested allows to prolongate solutions
and restrict residuals between two independent meshes. By passing a sequence
of such two-level transfer operators to MGTransferMF and, finally, to Multigrid,
non-nested multigrid can be computed.
- Utilities::MPI::RemotePointEvaluation can be used to couple non-matching
grids via surfaces (example: fluid-structure interaction, independently created
grids). The evaluation points can come from any side (pointwise interpolation)
or are defined on intersected meshes (Nitsche-type mortaring
@cite heinz2022high). Concerning the creation of such intersected meshes and the
corresponding evaluation points, see
GridTools::internal::distributed\_compute_intersection_locations().
- Alternatively to the coupling via Utilities::MPI::RemotePointEvaluation,
preCICE @cite bungartz2016precice @cite chourdakis2021precice can be used. The
code-gallery program "Laplace equation coupled to an external simulation
program" shows how to use this library with deal.II.
 *
 *
<a name="step_87-PlainProg"></a>
<h1> The plain program</h1>
@include "step-87.cc"
*/