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();
}
}
```
[email protected] schrieb am Montag, 5. Mai 2025 um 21:27:39 UTC+2:
> 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.
>
> > On 5 May 2025, at 15:49, Nils Schween <[email protected]> wrote:
> >
> > 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 <+49%206221%20516557>
> > Mail: [email protected]
> > 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
> >
> > --
> > 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 [email protected].
> > To view this discussion visit
> https://groups.google.com/d/msgid/dealii/87o6w7p3pm.fsf%40mpi-hd.mpg.de.
>
>
--
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 [email protected].
To view this discussion visit
https://groups.google.com/d/msgid/dealii/6162529a-4275-47d9-8d7e-583dbde8938an%40googlegroups.com.
