Integrating along one dimension of a finite elment solution

[Nils Schween]

Dear deal.II community,

I am looking for a way to integrate u_h(x,y,z) = \sum \alpha_i
\phi_i(x,y,z) over a single dimension, say z, i.e. a way to perform
\tilde{u}_h = \int u_h(x,y,z) dz. u_h is the finite element solution of
the partial differential equation that I am solving. As usual, in my
application u_h is represented in the form of a solution vector that
contains the coefficients \alpha_i corresponding to the basis functions
\phi_i of the finite element space I chose for the spatial
discretization. The integration can be considered a post processing
step.

My application is using the MPI capabilities of the deal.II library and,
hence, the solution vector is distributed among many cores/nodes. It
seems to me that this complicates the computation of the average, i.e.
the above integral, because when integrating along one direction, I
need information located on many different cores.

I looked through the tutorials, the FAQ, the google user group and
searched different places in the library, if there are functions that
integrate out one direction/dimension. However, I was unlucky.

I would like to ask, if I overlooked something, if anyone did implement
something like this already or if anyone could give a rough strategy how
he/she would go about implementing this average.

I am looking forward to your answers and thank you in advance,
Nils

--
Nils Schween

Phone: +49 6221 516 557
Mail: nils.s...@mpi-hd.mpg.de
PGP-Key: 4DD3DCC0532EE96DB0C1F8B5368DBFA14CB81849

Max Planck Institute for Nuclear Physics
Astrophysical Plasma Theory (APT)
Saupfercheckweg 1, D-69117 Heidelberg
https://www.mpi-hd.mpg.de/mpi/en/research/scientific-divisions-and-groups/independent-research-groups/apt

[Luca Heltai]

Dear Nils,

I think you want to take a look at the functionalities of Utilities::MPI::RemotePointEvaluation

https://www.dealii.org/current/doxygen/deal.II/classUtilities_1_1MPI_1_1RemotePointEvaluation.html

What you want to do is to generate a one dimensional quadrature formula of your choice, collecting all the quadrature points where you need to evaluate your functions, and then call


RemotePointEvaluation rpe;
rpe.reinit(qpoints, tria, mapping);
auto v = rpe.evaluate_and_process<double>(…);

Take a look at the directory

https://github.com/dealii/dealii/blob/master/tests/remote_point_evaluation/

to see examples of how to use this.

L.

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> To view this discussion visit https://groups.google.com/d/msgid/dealii/87o6w7p3pm.fsf%40mpi-hd.mpg.de.

[Magdalena]

Dear Nils,

in addition to Luca's response, I would recommend going through step-87, which gives an overview on point evaluation functionalities for distributed meshes.

Also, please find below related code snippets from our user code. We implemented a prototype to compute averaged values across an interface thickness by integrating along a line (normal to the interface).

The main steps within the code snippet are:

1. Generate points along a line (in our case it was along the normal vector arround a narrow band in a implicitly given interface)
2. Setup `RemotePointEvaluation` with the given point cloud and use `VectorTools::point_values` to get the values at the request points
3. Perform the numerical integration along the line
4. Fill a DoF vector with the averaged values.

Feel free to reach out if you'd like to discuss anything in more detail — I’d be happy to help!

Code snippets:
 

```cpp
  /**
   * 1d numerical integration of @p vals given at @p points using a trapezoidal rule.
   */
  template <int dim, typename number>
  number
  integrate_over_line(const std::vector<number>             &vals,
                      const std::vector<dealii::Point<dim>> &points)
  {
    number result = 0;
    for (unsigned int i = 0; i < vals.size() - 1; ++i)
      result += (vals[i] + vals[i + 1]) / 2. * (points[i + 1].distance(points[i]));
    return result;
  }  

template <int dim, typename number>
  void
  EvaporationMassFluxOperatorThicknessIntegration<dim, number>::compute_evaporative_mass_flux(
    VectorType       &evaporative_mass_flux,
    const VectorType &temperature) const
  {
    Assert(dim > 1, ExcMessage("The requested operation can only be performed for dim>1."));

    if constexpr (dim > 1)
      {
        scratch_data.initialize_dof_vector(evaporative_mass_flux, evapor_mass_flux_dof_idx);
       
/*
         * 1) generate point cloud normal to interface
         */
        std::vector<Point<dim>>   global_points_normal_to_interface;
        std::vector<unsigned int> global_points_normal_to_interface_pointer;

        const auto thickness_integration_band =
          5 * reinit_data.compute_interface_thickness_parameter_epsilon(
                scratch_data.get_min_cell_size());

        Assert(thickness_integration_band > 0.0,
               ExcMessage("Thickness for interface integration not set."));

        LevelSet::Tools::generate_points_along_normal<dim>(
          global_points_normal_to_interface,
          global_points_normal_to_interface_pointer,
          scratch_data.get_dof_handler(ls_hanging_nodes_dof_idx),
          fe_dim,
          scratch_data.get_mapping(),
          level_set_as_heaviside,
          normal_vector,
          thickness_integration_band,
          thickness_integration_data.subdivisions_per_side,
          /* bidirectional */ true,
          /* contour_value */ 0.5,
          thickness_integration_data.subdivisions_MCA);

        const bool update_ghosts = !level_set_as_heaviside.has_ghost_elements();
        if (update_ghosts)
          level_set_as_heaviside.update_ghost_values();
        const bool heat_update_ghosts = !temperature.has_ghost_elements();
        if (heat_update_ghosts)
          temperature.update_ghost_values();

        /*
         * 2) setup remote point evaluation and evaluate result at points along line
         */
        Utilities::MPI::RemotePointEvaluation<dim, dim> remote_point_evaluation(
          1e-6 /*tolerance*/, true /*unique mapping*/);

        remote_point_evaluation.reinit(global_points_normal_to_interface,
                                       scratch_data.get_triangulation(),
                                       scratch_data.get_mapping());


        const auto temperature_evaluation_values =
          dealii::VectorTools::point_values<1>(remote_point_evaluation,
                                               scratch_data.get_dof_handler(
                                                 heat_hanging_nodes_dof_idx),
                                               temperature);

        const std::vector<Tensor<1, dim, number>> level_set_gradient_values =
          dealii::VectorTools::point_gradients<1>(remote_point_evaluation,
                                                  scratch_data.get_dof_handler(
                                                    ls_hanging_nodes_dof_idx),
                                                  level_set_as_heaviside);
        /*
         * 3) perform line integral
         */
        std::vector<number> mass_flux_val;
        mass_flux_val.resize(global_points_normal_to_interface_pointer.back());

        // loop over points located at the interface and determined by means of MCA
        for (unsigned int i = 0; i < global_points_normal_to_interface_pointer.size() - 1; ++i)
          {
            const auto start = global_points_normal_to_interface_pointer[i];
            const auto end   = global_points_normal_to_interface_pointer[i + 1];
            const auto size  = end - start;

            std::vector<number> func_vals(size);

            // loop over all points along normal at one MC point
            for (unsigned int l = 0; l < size; ++l)
              {
                func_vals[l] = evaporation_model.local_compute_evaporative_mass_flux(
                                 temperature_evaluation_values[start + l]) *
                               level_set_gradient_values[start + l].norm();
              }

            number mass_flux_average = UtilityFunctions::integrate_over_line<dim>(
              func_vals,
              std::vector(global_points_normal_to_interface.begin() + start,
                          global_points_normal_to_interface.begin() + end));

            /*
             * set for each point along the interface the value equal to the one determined
             * by the line integral
             */
            for (unsigned int l = 0; l < size; ++l)
              mass_flux_val[start + l] = mass_flux_average;
          }

        /*
         * 4) Broadcast values determined along the normal direction to nodal points by means
         * of a simple averaging procedure. Subsequently, a DoF-vector is filled.
         *
         * @note We also tried the extrapolation via a radial basis function, however it was
         * too sensitive with respect to the inherent distance parameter.
         */
        VectorType vector_multiplicity;
        vector_multiplicity.reinit(evaporative_mass_flux);

        const auto broadcast_function = [&](const auto &values, const auto &cell_data) {
          // loop over all cells where points along normal are interior
          for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
            {
              typename DoFHandler<dim>::active_cell_iterator cell = {
                &remote_point_evaluation.get_triangulation(),
                cell_data.cells[i].first,
                cell_data.cells[i].second,
                &scratch_data.get_dof_handler(ls_hanging_nodes_dof_idx)};

              // values at the points along normal in reference coordinates
              const ArrayView<const number> heat_vals(values.data() +
                                                        cell_data.reference_point_ptrs[i],
                                                      cell_data.reference_point_ptrs[i + 1] -
                                                        cell_data.reference_point_ptrs[i]);

              // extract dof indices of cell
              std::vector<types::global_dof_index> local_dof_indices(
                cell->get_fe().n_dofs_per_cell());
              cell->get_dof_indices(local_dof_indices);

              /*
               * Compute mean value from the values of points within the cell and set the values at
               * the cell nodes equal to the latter.
               */
              Vector<number> nodal_values(cell->get_fe().n_dofs_per_cell());

              const number average =
                std::accumulate(heat_vals.begin(), heat_vals.end(), 0.0) / heat_vals.size();

              for (auto &n : nodal_values)
                n = average;

              scratch_data.get_constraint(evapor_mass_flux_dof_idx)
                .distribute_local_to_global(nodal_values, local_dof_indices, evaporative_mass_flux);

              // count the number of written values (multiplicity)
              for (auto &val : nodal_values)
                val = 1.0;

              scratch_data.get_constraint(ls_hanging_nodes_dof_idx)
                .distribute_local_to_global(nodal_values, local_dof_indices, vector_multiplicity);
            }
        };

        // fill dof vector
        std::vector<number> buffer;

        remote_point_evaluation.template process_and_evaluate<number>(mass_flux_val,
                                                                      buffer,
                                                                      broadcast_function);

        // take care dofs shared by multiple geometric entities
        evaporative_mass_flux.compress(VectorOperation::add);
        vector_multiplicity.compress(VectorOperation::add);

        for (unsigned int i = 0; i < vector_multiplicity.locally_owned_size(); ++i)
          if (vector_multiplicity.local_element(i) > 1.0)
            evaporative_mass_flux.local_element(i) /= vector_multiplicity.local_element(i);

        scratch_data.get_constraint(evapor_mass_flux_dof_idx).distribute(evaporative_mass_flux);

        if (heat_update_ghosts)
          temperature.zero_out_ghost_values();

        if (update_ghosts)
          level_set_as_heaviside.zero_out_ghost_values();
      }
  }
```  

[Nils Schween]

Dear Luca, Dear Magdalena,

thank you very, very much. This is extremely useful.
@Luca thanks for pointing me to the tests. @Magdalena thanks for sharing
your code. I particularly like the idea to fill a DoF Vector. This
allows me to use the DataOut routines. It's likely that I will come back to
you with more questions while implementing our version; thanks for the
offer. When I am done, I will share our implementation as well.

Once more a big thank you,
Nils

--
Nils Schween
PhD Student