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 Magdalena <[email protected]> writes: > 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 > > 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. -- Nils Schween PhD Student Phone: +49 6221 516 557 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/874ixy43e7.fsf%40mpi-hd.mpg.de.
smime.p7s
Description: S/MIME cryptographic signature
