From 538090ffffa77c93f5919952433c3be4509e8bd6 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 29 Jun 2025 22:00:49 -0600 Subject: [PATCH 01/10] Implement an interface that external surface deformation tools can use. --- .../external_tool_interface.h | 191 ++++++++++++++++++ .../external_tool_interface.cc | 159 +++++++++++++++ 2 files changed, 350 insertions(+) create mode 100644 include/aspect/mesh_deformation/external_tool_interface.h create mode 100644 source/mesh_deformation/external_tool_interface.cc diff --git a/include/aspect/mesh_deformation/external_tool_interface.h b/include/aspect/mesh_deformation/external_tool_interface.h new file mode 100644 index 00000000000..8133ec5c085 --- /dev/null +++ b/include/aspect/mesh_deformation/external_tool_interface.h @@ -0,0 +1,191 @@ +/* + Copyright (C) 2011 - 2024 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + + +#ifndef _aspect_mesh_deformation_external_tool_interface_h +#define _aspect_mesh_deformation_external_tool_interface_h + +#include + + +namespace aspect +{ + namespace MeshDeformation + { + /** + * This class provides some support for writing classes derived + * from MeshDeformation::Interface when implementing surface + * deformation models that are based on external tools such as + * Fastscape, Landlab, OpenLEM, etc. The primary complication of + * interacting with these external tools is that they generally + * use their own discretization of surface processes, using a mesh + * that, in most cases, will be different from the one used by + * ASPECT. (Of course, ASPECT's use is a *volume* mesh, whereas + * surface processes use *surface* meshes. ASPECT's own surface + * diffusion, for example, works on the surface faces of ASPECT's + * mesh, but in general, external tools use a mesh that is + * *entirely unrelated* to the surface cells of ASPECT's volume + * mesh.) In these cases, one needs to implement transfer + * operations to take ASPECT's solution and interpolate it to the + * external tool's mesh, run some surface evolution steps, and + * then interpolate back to the surface nodes of ASPECT's + * mesh. These interpolation operations are difficult to implement + * in their own right, but are made particularly cumbersome in + * parallel because then, the points at which the external tool + * wants to know the solution on any given MPI process will, in + * general, not be located on that part of ASPECT's mesh that is + * owned by that MPI process. As a consequence, implementing the + * interpolation requires finding the MPI process that owns the + * cell on which that point is located, and then communication + * back and forth. + * + * @sect3{Helper functions} + * + * This class implements the interpolation functionality for + * derived classes to use. It works through a two-stage approach: + * First, derived classes declare at which points they require + * ASPECT's solution to be evaluated, via the + * set_evaluation_points() function. This function then finds + * which process owns the cell around this point, and sets up + * communication structures that will make later evaluation + * efficient. Second, this class provides the + * evaluate_aspect_variables_at_points() function that uses these + * communication structures to evaluate ASPECT's current solution + * at the points previously set. The class also provides the + * interpolate_vertical_velocities_to_surface_points() function + * that takes a set of (vertical) velocity values at these points + * and uses them to interpolate the information back onto ASPECT's + * mesh. + * + * All three of these functions are `protected` member functions + * of this class, ready to be called by derived classes. + * + * + * @sect3{A high-level function} + * + * All classes derived from Interface need to implement the + * Interface::compute_velocity_constraints_on_boundary() + * function. Because the basic outline of this function looks + * essentially the same for all external tools, this class also + * provides a high-level implementation of this function is also + * provided as part of this class. In essence, it looks as + * follows (pseudo-code), relying on the + * `compute_updated_velocities_at_points()` function that gets + * ASPECT's solution at the evaluation points and returns velocity + * vectors at all of these evaluation points: + * @code + * // Interpolate ASPECT's solution at evaluation points: + * aspect_surface_velocities = evaluate_aspect_variables_at_points(); + * + * // Call derive class's method to compute updated velocities: + * external_surface_velocities = compute_updated_velocities_at_points(aspect_surface_velocities); + * + * // Turn the result into the constraints that + * // compute_velocity_constraints_on_boundary is supposed + * // to return. + * @endcode + */ + template + class ExternalToolInterface : public Interface, public SimulatorAccess + { + public: + virtual + void + compute_velocity_constraints_on_boundary(const DoFHandler &mesh_deformation_dof_handler, + AffineConstraints &mesh_velocity_constraints, + const std::set &boundary_ids) const override; + + virtual + std::vector> + compute_updated_velocities_at_points (const std::vector> ¤t_solution_at_points) const = 0; + + + protected: + // Helper functions to be used in derived classes: + + /** + * Declare at which points the external tool driven by a class + * derived from the current one needs to evaluate the ASPECT + * solution. Each process will call this function with the points + * it needs. Different processes will, in general, provide + * distinct sets of points. + * + * The points provided here are given in the undeformed + * configuration that corresponds to the initial mesh. For + * example, if the mesh describes a box geometry, then the + * points should like in the $x$-$y$-plane that forms the top + * surface, rather than on the currently deformed top surface + * that describes the current elevation map. + * + * @note This function sets up communication structures that + * encode, among other things, which process owns which + * points. This information becomes outdated when the + * ASPECT mesh changes for mesh refinement. As a + * consequence, the current function sets up a process that + * invalidates the communication structures upon mesh + * refinement. Derived classes therefore have to set up + * their own code to call set_evaluation_points() again when + * the ASPECT mesh changes, or perhaps simply check whether + * the information needs to be updated in an overload of the + * update() function called at the beginning of each time + * step. + */ + void + set_evaluation_points (const std::vector> &evaluation_points); + + /** + * Return the value of the ASPECT solution at the set of points + * previously set via set_evaluation_points(). + * + * For each point, the return object contains a vector with as + * many components as there are ASPECT solution components. + */ + std::vector> + evaluate_aspect_variables_at_points () const; + + /** + * Given velocities (typically vertical, but not always) at + * all of the points this process has previously set via + * set_evaluation_points(), compute a (global) finite element + * field vector that in the velocity components of surface + * nodes corresponds to an interpolation of the velocities + * provided. + * + * The output of this function can then be used as input for a + * function that implements the + * compute_velocity_constraints_on_boundary() function of the + * base class. + */ + LinearAlgebra::Vector + interpolate_velocities_to_surface_points (const std::vector> &vertical_velocities) const; + + + private: + + std::vector> evaluation_points; + + // Timo: Replace by whatever type you need here + class X {}; + std::unique_ptr remote_point_evaluator; + }; + } +} + +#endif diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc new file mode 100644 index 00000000000..1f3931183f8 --- /dev/null +++ b/source/mesh_deformation/external_tool_interface.cc @@ -0,0 +1,159 @@ +/* + Copyright (C) 2014 - 2024 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . + */ + + +#include + +namespace aspect +{ + namespace MeshDeformation + { + template + void + ExternalToolInterface:: + compute_velocity_constraints_on_boundary(const DoFHandler &mesh_deformation_dof_handler, + AffineConstraints &mesh_velocity_constraints, + const std::set &boundary_ids) const + { + // First compute a (global) vector that has the correct velocities + // set at all boundary nodes: + const std::vector> aspect_surface_velocities = evaluate_aspect_variables_at_points(); + + const std::vector> external_surface_velocities + = compute_updated_velocities_at_points(aspect_surface_velocities); + + const LinearAlgebra::Vector v_interpolated + = interpolate_velocities_to_surface_points(external_surface_velocities); + + // Turn v_interpolated into constraints. For this, loop over all + // boundary DoFs and if a boundary DoF is locally owned, create a + // constraint. We later make that consistent across processors to + // ensure we also know about the locally relevant DoFs' + // constraints: + std::vector face_dof_indices (mesh_deformation_dof_handler.get_fe().dofs_per_face); + for (const auto &cell : mesh_deformation_dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + for (const auto &face : cell->face_iterators()) + if (face->at_boundary() && + (boundary_ids.find(face->boundary_id()) != boundary_ids.end())) + { + face->get_dof_indices(face_dof_indices); + + for (types::global_dof_index i : face_dof_indices) + if (v_interpolated.locally_owned_elements().is_element(i)) + mesh_velocity_constraints.add_constraint (i, {}, v_interpolated(i)); + } + mesh_velocity_constraints.make_consistent_in_parallel(mesh_deformation_dof_handler.locally_owned_dofs(), + DoFTools::extract_locally_relevant_dofs(mesh_deformation_dof_handler), + this->get_mpi_communicator()); + + } + + + template + void + ExternalToolInterface:: + set_evaluation_points (const std::vector> &evaluation_points) + { + // First, save a copy of the points at which we need the solution, + // among other reasons so that we can track that input arguments + // for later function calls describe the same number of points. + this->evaluation_points = evaluation_points; + + // Then also invalidate the previous evaluator: + remote_point_evaluator.reset(); + + // Timo: + // TODO: Implement set-up phase via MPIRemotePointEvaluation + + + // Finally, also ensure that upon mesh refinement, all of the + // information set herein is invalidated: + this->get_signals().pre_refinement_store_user_data + .connect([this](typename parallel::distributed::Triangulation &) + { + this->evaluation_points.clear(); + this->remote_point_evaluator.reset(); + } + ); + } + + + + template + std::vector> + ExternalToolInterface:: + evaluate_aspect_variables_at_points () const + { + Assert (remote_point_evaluator != nullptr, + ExcMessage("You can only call this function if you have previously " + "set the evaluation points by calling set_evaluation_points(), " + "and if the evaluator has not been invalidated by a mesh " + "refinement step.")); + + std::vector> solution_at_points (evaluation_points.size(), + std::vector(this->introspection().n_components)); + // Timo: Implement via MPIRemoteEvaluation + + return solution_at_points; + } + + + + template + LinearAlgebra::Vector + ExternalToolInterface:: + interpolate_velocities_to_surface_points (const std::vector> &velocities) const + { + Assert (remote_point_evaluator != nullptr, + ExcMessage("You can only call this function if you have previously " + "set the evaluation points by calling set_evaluation_points(), " + "and if the evaluator has not been invalidated by a mesh " + "refinement step.")); + AssertDimension(velocities.size(), evaluation_points.size()); + + // Create the output vector. TODO: We need to get access to the + // locally owned DoFs index set. Look up how this is done in the + // other implementations of the Interface base class, if they do + // it this way at all. + + // LinearAlgebra::Vector vector_with_surface_velocities(this->get_mesh_deformation_handler().mesh_deformation_dof_handler.locally_owned_dofs(), + // mpi_communicator) + LinearAlgebra::Vector vector_with_surface_velocities; + + + // Timo: Implement via MPIRemoteEvaluation + + return vector_with_surface_velocities; + } + } + + + + namespace MeshDeformation + { +#define INSTANTIATE(dim) \ + template class ExternalToolInterface; + + ASPECT_INSTANTIATE(INSTANTIATE) + +#undef INSTANTIATE + } +} From a713a3debf3b245febacba96689b6465d5129ba3 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 30 Jun 2025 08:17:23 -0600 Subject: [PATCH 02/10] Use the same style to compute constraints as in all other mesh deformation plugins. --- .../external_tool_interface.cc | 35 +++++++++++-------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc index 1f3931183f8..27292688781 100644 --- a/source/mesh_deformation/external_tool_interface.cc +++ b/source/mesh_deformation/external_tool_interface.cc @@ -47,23 +47,28 @@ namespace aspect // constraint. We later make that consistent across processors to // ensure we also know about the locally relevant DoFs' // constraints: - std::vector face_dof_indices (mesh_deformation_dof_handler.get_fe().dofs_per_face); - for (const auto &cell : mesh_deformation_dof_handler.active_cell_iterators()) - if (cell->is_locally_owned()) - for (const auto &face : cell->face_iterators()) - if (face->at_boundary() && - (boundary_ids.find(face->boundary_id()) != boundary_ids.end())) + // now insert the relevant part of the solution into the mesh constraints + const IndexSet constrained_dofs = + DoFTools::extract_boundary_dofs(mesh_deformation_dof_handler, + ComponentMask(dim, true), + boundary_ids); + + for (unsigned int i = 0; i < constrained_dofs.n_elements(); ++i) + { + types::global_dof_index index = constrained_dofs.nth_index_in_set(i); + if (mesh_velocity_constraints.can_store_line(index)) + if (mesh_velocity_constraints.is_constrained(index)==false) { - face->get_dof_indices(face_dof_indices); - - for (types::global_dof_index i : face_dof_indices) - if (v_interpolated.locally_owned_elements().is_element(i)) - mesh_velocity_constraints.add_constraint (i, {}, v_interpolated(i)); +#if DEAL_II_VERSION_GTE(9,6,0) + mesh_velocity_constraints.add_constraint(index, + {}, + v_interpolated(index)); +#else + mesh_velocity_constraints.add_line(index); + mesh_velocity_constraints.set_inhomogeneity(index, v_interpolated(index)); +#endif } - mesh_velocity_constraints.make_consistent_in_parallel(mesh_deformation_dof_handler.locally_owned_dofs(), - DoFTools::extract_locally_relevant_dofs(mesh_deformation_dof_handler), - this->get_mpi_communicator()); - + } } From 869729046708ab87036735527b9c6538cf90a9cf Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 30 Jun 2025 08:23:15 -0600 Subject: [PATCH 03/10] Make index handling more efficient. --- source/mesh_deformation/external_tool_interface.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc index 27292688781..74225a52fc3 100644 --- a/source/mesh_deformation/external_tool_interface.cc +++ b/source/mesh_deformation/external_tool_interface.cc @@ -53,9 +53,8 @@ namespace aspect ComponentMask(dim, true), boundary_ids); - for (unsigned int i = 0; i < constrained_dofs.n_elements(); ++i) + for (const types::global_dof_index index : constrained_dofs) { - types::global_dof_index index = constrained_dofs.nth_index_in_set(i); if (mesh_velocity_constraints.can_store_line(index)) if (mesh_velocity_constraints.is_constrained(index)==false) { From 724f70ba6bf7f30f78ec5d95a05410a20cc50373 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 30 Jun 2025 11:07:08 -0600 Subject: [PATCH 04/10] Add necessary header include. --- source/mesh_deformation/external_tool_interface.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc index 74225a52fc3..0e72f3e2ab8 100644 --- a/source/mesh_deformation/external_tool_interface.cc +++ b/source/mesh_deformation/external_tool_interface.cc @@ -20,6 +20,8 @@ #include +#include + namespace aspect { @@ -144,6 +146,8 @@ namespace aspect // Timo: Implement via MPIRemoteEvaluation + (void)velocities; + return vector_with_surface_velocities; } From 5b983f57d6a6c598002737effd76134f7cddeaa1 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Tue, 1 Jul 2025 17:41:29 -0400 Subject: [PATCH 05/10] add RPE --- .../external_tool_interface.h | 12 +- .../external_tool_interface.cc | 391 +++++++++++++++++- tests/mesh_deformation_external_01.cc | 162 ++++++++ tests/mesh_deformation_external_01.prm | 96 +++++ .../screen-output | 159 +++++++ tests/mesh_deformation_external_01/statistics | 19 + tests/mesh_deformation_external_01_3d.cc | 21 + tests/mesh_deformation_external_01_3d.prm | 97 +++++ .../screen-output | 75 ++++ .../statistics | 13 + 10 files changed, 1020 insertions(+), 25 deletions(-) create mode 100644 tests/mesh_deformation_external_01.cc create mode 100644 tests/mesh_deformation_external_01.prm create mode 100644 tests/mesh_deformation_external_01/screen-output create mode 100644 tests/mesh_deformation_external_01/statistics create mode 100644 tests/mesh_deformation_external_01_3d.cc create mode 100644 tests/mesh_deformation_external_01_3d.prm create mode 100644 tests/mesh_deformation_external_01_3d/screen-output create mode 100644 tests/mesh_deformation_external_01_3d/statistics diff --git a/include/aspect/mesh_deformation/external_tool_interface.h b/include/aspect/mesh_deformation/external_tool_interface.h index 8133ec5c085..83dceec22c9 100644 --- a/include/aspect/mesh_deformation/external_tool_interface.h +++ b/include/aspect/mesh_deformation/external_tool_interface.h @@ -112,6 +112,11 @@ namespace aspect AffineConstraints &mesh_velocity_constraints, const std::set &boundary_ids) const override; + /** + * Given all solution variables at each surface point, compute velocities at these points. + * + * This needs to be implemented by the derived class and implement the "surface evolution". + */ virtual std::vector> compute_updated_velocities_at_points (const std::vector> ¤t_solution_at_points) const = 0; @@ -176,14 +181,11 @@ namespace aspect LinearAlgebra::Vector interpolate_velocities_to_surface_points (const std::vector> &vertical_velocities) const; - - private: + protected: std::vector> evaluation_points; - // Timo: Replace by whatever type you need here - class X {}; - std::unique_ptr remote_point_evaluator; + std::unique_ptr> remote_point_evaluator; }; } } diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc index 0e72f3e2ab8..c37b67e0414 100644 --- a/source/mesh_deformation/external_tool_interface.cc +++ b/source/mesh_deformation/external_tool_interface.cc @@ -22,11 +22,287 @@ #include #include +#include + +/* +TODO: +- PointDataOut: move to deal.II +- initial topography from external tool: for now compute_initial_deformation_on_boundary(), later refactor + to provide FE vector +*/ namespace aspect { namespace MeshDeformation { + + + + template + class PointDataOut : public dealii::DataOutInterface<0, spacedim> + { + public: + /** + * Default constructor. + */ + PointDataOut() = default; + + /** + * Default destructor. + */ + ~PointDataOut() = default; + + + /** + * Build the patches for a given set of points and optionally data in each point. + * + * @param [in] locations The point locations. + * @param [in] data A vector of data values for each point. + * @param [in] data_component_names An optional vector of strings that + * describe the properties of each datum. + * @param [in] data_component_interpretations An optional vector that + * controls if the properties are interpreted as scalars, vectors, + * or tensors. Has to be of the same length as @p data_component_names. + */ + void + build_patches(const std::vector> &locations, + const std::vector> &data = {}, + const std::vector &data_component_names = {}, + const std::vector< + DataComponentInterpretation::DataComponentInterpretation> + &data_component_interpretations_ = {}) + { + Assert( + data_component_names.size() == data_component_interpretations_.size(), + ExcMessage( + "When calling Particles::DataOut::build_patches with data component " + "names and interpretations you need to provide as many data component " + "names as interpretations. Provide the same name for components that " + "belong to a single vector or tensor.")); + + dataset_names.clear(); + dataset_names.emplace_back("id"); + dataset_names.insert(dataset_names.end(), + data_component_names.begin(), + data_component_names.end()); + + data_component_interpretations.clear(); + data_component_interpretations.emplace_back( + DataComponentInterpretation::component_is_scalar); + data_component_interpretations.insert( + data_component_interpretations.end(), + data_component_interpretations_.begin(), + data_component_interpretations_.end()); + + const unsigned int n_property_components = data_component_names.size(); + const unsigned int n_data_components = dataset_names.size(); + + patches.resize(locations.size()); + + for (unsigned int i = 0; i < locations.size(); ++i) + { + patches[i].vertices[0] = locations[i]; + patches[i].patch_index = i; + + // We have one more data components than dataset_names (the particle id) + patches[i].data.reinit(n_data_components, 1); + + patches[i].data(0, 0) = i; // id + + if (n_data_components > 1) + { + for (unsigned int property_index = 0; + property_index < n_property_components; + ++property_index) + patches[i].data(property_index + 1, 0) = + data[i][property_index]; + } + } + } + + protected: + /** + * Returns the patches built by the data_out class which was previously + * built using a particle handler + */ + virtual const std::vector> & + get_patches() const override + { + return patches; + } + + /** + * Virtual function through which the names of data sets are obtained from + * this class + */ + virtual std::vector + get_dataset_names() const override + { + return dataset_names; + } + + + /** + * Overload of the respective DataOutInterface::get_nonscalar_data_ranges() + * function. See there for a more extensive documentation. + * This function is a reimplementation of the function + * DataOut_DoFData::get_nonscalar_data_ranges(). + */ + virtual std::vector< + std::tuple> + get_nonscalar_data_ranges() const override + { + std::vector< + std::tuple> + ranges; + + // Make sure the data structures were set up correctly. Since they + // can only be filled by build_patches() above, they should have + // been checked already. + Assert(dataset_names.size() == data_component_interpretations.size(), + ExcInternalError()); + + // collect the ranges of particle data + const unsigned int n_output_components = + data_component_interpretations.size(); + unsigned int output_component = 0; + for (unsigned int i = 0; i < n_output_components; /* i is updated below */) + // see what kind of data we have here. note that for the purpose of the + // current function all we care about is vector data + switch (data_component_interpretations[i]) + { + case DataComponentInterpretation::component_is_scalar: + { + // Just move component forward by one + ++i; + ++output_component; + + break; + } + case DataComponentInterpretation::component_is_part_of_vector: + { + // ensure that there is a continuous number of next space_dim + // components that all deal with vectors + Assert( + i + spacedim <= n_output_components, + Exceptions::DataOutImplementation::ExcInvalidVectorDeclaration( + i, dataset_names[i])); + for (unsigned int dd = 1; dd < spacedim; ++dd) + Assert( + data_component_interpretations[i + dd] == + DataComponentInterpretation::component_is_part_of_vector, + Exceptions::DataOutImplementation:: + ExcInvalidVectorDeclaration(i, dataset_names[i])); + + // all seems right, so figure out whether there is a common + // name to these components. if not, leave the name empty and + // let the output format writer decide what to do here + std::string name = dataset_names[i]; + for (unsigned int dd = 1; dd < spacedim; ++dd) + if (name != dataset_names[i + dd]) + { + name = ""; + break; + } + + // Finally add a corresponding range. + // + // This sort of logic is also explained in some detail in + // DataOut::build_one_patch(). + ranges.emplace_back(std::forward_as_tuple( + output_component, + output_component + spacedim - 1, + name, + DataComponentInterpretation::component_is_part_of_vector)); + + // increase the 'component' counter by the appropriate amount, + // same for 'i', since we have already dealt with all these + // components + output_component += spacedim; + i += spacedim; + + break; + } + + case DataComponentInterpretation::component_is_part_of_tensor: + { + const unsigned int size = spacedim * spacedim; + // ensure that there is a continuous number of next + // spacedim*spacedim components that all deal with tensors + Assert( + i + size <= n_output_components, + Exceptions::DataOutImplementation::ExcInvalidTensorDeclaration( + i, dataset_names[i])); + for (unsigned int dd = 1; dd < size; ++dd) + Assert( + data_component_interpretations[i + dd] == + DataComponentInterpretation::component_is_part_of_tensor, + Exceptions::DataOutImplementation:: + ExcInvalidTensorDeclaration(i, dataset_names[i])); + + // all seems right, so figure out whether there is a common + // name to these components. if not, leave the name empty and + // let the output format writer decide what to do here + std::string name = dataset_names[i]; + for (unsigned int dd = 1; dd < size; ++dd) + if (name != dataset_names[i + dd]) + { + name = ""; + break; + } + + // Finally add a corresponding range. + ranges.emplace_back(std::forward_as_tuple( + output_component, + output_component + size - 1, + name, + DataComponentInterpretation::component_is_part_of_tensor)); + + // increase the 'component' counter by the appropriate amount, + // same for 'i', since we have already dealt with all these + // components + output_component += size; + i += size; + break; + } + + default: + Assert(false, ExcNotImplemented()); + } + return ranges; + } + + private: + /** + * This is a vector of patches that is created each time build_patches() is + * called. These patches are used in the output routines of the base + * classes. + */ + std::vector> patches; + + /** + * A vector of field names for all data components stored in patches. + */ + std::vector dataset_names; + + /** + * A vector that for each of the data components of the + * current data set indicates whether they are scalar fields, parts of a + * vector-field, or any of the other supported kinds of data. + */ + std::vector + data_component_interpretations; + }; + + + + template void ExternalToolInterface:: @@ -44,6 +320,8 @@ namespace aspect const LinearAlgebra::Vector v_interpolated = interpolate_velocities_to_surface_points(external_surface_velocities); + // TODO: need ghost values of v_interpolated? + // Turn v_interpolated into constraints. For this, loop over all // boundary DoFs and if a boundary DoF is locally owned, create a // constraint. We later make that consistent across processors to @@ -70,6 +348,8 @@ namespace aspect #endif } } + + // TODO: make consistent? } @@ -83,12 +363,11 @@ namespace aspect // for later function calls describe the same number of points. this->evaluation_points = evaluation_points; - // Then also invalidate the previous evaluator: - remote_point_evaluator.reset(); - - // Timo: - // TODO: Implement set-up phase via MPIRemotePointEvaluation - + // Set up RemotePointEvaluation with a Mapping of the undeformed mesh: + remote_point_evaluator = std::make_unique>(); + // TODO: does this need to be a higher order mapping for spherical problems? + static MappingQ mapping(1); + remote_point_evaluator->reinit(this->evaluation_points, this->get_triangulation(), mapping); // Finally, also ensure that upon mesh refinement, all of the // information set herein is invalidated: @@ -97,8 +376,7 @@ namespace aspect { this->evaluation_points.clear(); this->remote_point_evaluator.reset(); - } - ); + }); } @@ -114,9 +392,23 @@ namespace aspect "and if the evaluator has not been invalidated by a mesh " "refinement step.")); - std::vector> solution_at_points (evaluation_points.size(), - std::vector(this->introspection().n_components)); - // Timo: Implement via MPIRemoteEvaluation + const unsigned int n_components = this->introspection().n_components; + std::vector> solution_at_points (evaluation_points.size(), std::vector(n_components, 0.0)); + + // VectorTools::point_values can evaluate N components at a time, but this is a template argument and not a + // runtime argument. For now, we just evaluate them one component at the time. Of course it would be more + // efficient to branch and evaluate up to K at a time (for a reasonable number of K, say 10). Maybe something + // to put directly into deal.II... + for (unsigned int c=0; c values = VectorTools::point_values<1>(*this->remote_point_evaluator, + this->get_dof_handler(), + this->get_solution(), + dealii::VectorTools::EvaluationFlags::avg, + c); + for (unsigned int i=0; i out; + // const auto &mapping = this->get_mapping(); + std::vector> real_evaluation_points(evaluation_points.size()); + std::vector> data(evaluation_points.size(), std::vector(dim, 0.0)); + for (unsigned int i=0; iget_mesh_deformation_handler().mesh_deformation_dof_handler.locally_owned_dofs(), - // mpi_communicator) - LinearAlgebra::Vector vector_with_surface_velocities; + const std::vector data_component_names(dim, "velocity"); + const std::vector data_component_interpretations(dim, DataComponentInterpretation::component_is_part_of_vector); + out.build_patches(real_evaluation_points, data, data_component_names, data_component_interpretations); + + out.write_vtu_with_pvtu_record(this->get_output_directory(), "surf_points", output_no, this->get_mpi_communicator(), 4, 0); + + ++output_no; + } + + + + // Create the output vector. + const DoFHandler &mesh_dof_handler = this->get_mesh_deformation_handler().get_mesh_deformation_dof_handler(); + LinearAlgebra::Vector vector_with_surface_velocities(mesh_dof_handler.locally_owned_dofs(), + this->get_mpi_communicator()); + LinearAlgebra::Vector one_over_distance_vec(mesh_dof_handler.locally_owned_dofs(), + this->get_mpi_communicator()); + + + const unsigned int dofs_per_cell = mesh_dof_handler.get_fe().dofs_per_cell; + std::vector cell_dof_indices (dofs_per_cell); + + const auto eval_func = [&](const ArrayView< const Tensor<1,dim>> &values, + const typename Utilities::MPI::RemotePointEvaluation::CellData &cell_data) + { + for (const auto cell_index : cell_data.cell_indices()) + { + const auto cell_dofs = + cell_data.get_active_cell_iterator(cell_index)->as_dof_handler_iterator( + mesh_dof_handler); + + const ArrayView> unit_points = cell_data.get_unit_points(cell_index); + const auto local_values = cell_data.get_data_view(cell_index, values); + + cell_dofs->get_dof_indices(cell_dof_indices); + const std::vector< Point< dim >> &support_points = mesh_dof_handler.get_fe().get_unit_support_points(); + for (unsigned int i=0; i one_over_distance_vec(cell_dof_indices[j])) + { + one_over_distance_vec(cell_dof_indices[j]) = one_over_distance; + const unsigned int c = mesh_dof_handler.get_fe().system_to_component_index(j).first; + vector_with_surface_velocities(cell_dof_indices[j]) = local_values[i][c]; + } + } + } + } + }; - // Timo: Implement via MPIRemoteEvaluation - (void)velocities; + this->remote_point_evaluator->template process_and_evaluate,1>(velocities, eval_func, /*sort_data*/ true); + one_over_distance_vec.compress(VectorOperation::insert); + vector_with_surface_velocities.compress(VectorOperation::insert); return vector_with_surface_velocities; } diff --git a/tests/mesh_deformation_external_01.cc b/tests/mesh_deformation_external_01.cc new file mode 100644 index 00000000000..ed803861046 --- /dev/null +++ b/tests/mesh_deformation_external_01.cc @@ -0,0 +1,162 @@ +/* + Copyright (C) 2025 - 2025 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + +#include +#include + +#include + +using namespace aspect; + +#include +#include +#include +#include + + +namespace aspect +{ + namespace MeshDeformation + { + template + class TestExternalDeformation : public ExternalToolInterface + { + public: + TestExternalDeformation() = default; + + void initialize() override + { + this->get_pcout() << "initialize()" << std::endl; + } + + void + update () override + { + this->get_pcout() << "update()" << std::endl; + + if (!this->remote_point_evaluator) + { + std::vector> points; + this->get_pcout() << "\tsetting points" << std::endl; + + if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator())==0) + { + if constexpr (dim == 2) + { + points.emplace_back(0.3, 1.0); + points.emplace_back(0.5, 1.0); + } + else + { + points.emplace_back(0.3, 0.7, 1.0); + points.emplace_back(0.5, 0.5, 1.0); + points.emplace_back(0.7, 0.2, 1.0); + points.emplace_back(0.7, 0.91, 1.0); + } + } + else + { + if constexpr (dim == 2) + { + points.emplace_back(0.72, 1.0); + } + else + { + points.emplace_back(0.9, 0.9, 1.0); + } + } + this->set_evaluation_points (points); + } + } + + virtual + std::vector> + compute_updated_velocities_at_points (const std::vector> ¤t_solution_at_points) const override + { + this->get_pcout() << "compute_updated_velocities_at_points()" << std::endl; + + { + // Copy all data to rank 0 to print to the screen. + this->get_pcout() << "Solution at evaluation points:" << std::endl; + + std::vector>> locations_by_rank = Utilities::MPI::gather(this->get_mpi_communicator(), this->evaluation_points); + std::vector>> data_by_rank = Utilities::MPI::gather(this->get_mpi_communicator(), current_solution_at_points); + + const unsigned int rank = Utilities::MPI::this_mpi_process(this->get_mpi_communicator()); + const unsigned int size = Utilities::MPI::n_mpi_processes(this->get_mpi_communicator()); + + if (rank == 0) + { + for (unsigned int r=0; revaluation_points.size(), ExcInternalError()); + std::vector> velocities(current_solution_at_points.size(), Tensor<1,dim>()); + if (velocities.size()>1) + { + velocities[0][dim-1]=30.0; + velocities[1][dim-1]=-3.5; + } + return velocities; + } + + + + Tensor<1,dim> + compute_initial_deformation_on_boundary(const types::boundary_id /*boundary_indicator*/, + const Point &position) const override + { + const Tensor<1,dim> gravity = this->get_gravity_model().gravity_vector(position); + Tensor<1,dim> topography_direction; + if (gravity.norm() > 0.0) + topography_direction = -gravity / gravity.norm(); + + const double topography_amplitude = (position[0]>=0.5) ? (0.05 * (1.+std::cos(2.*numbers::PI*position[0]))) : 0.0; + return topography_amplitude * topography_direction; + } + }; + } +} + + +// explicit instantiation of the functions we implement in this file +namespace aspect +{ + namespace MeshDeformation + { + ASPECT_REGISTER_MESH_DEFORMATION_MODEL(TestExternalDeformation, + "external deformation", + "") + } +} diff --git a/tests/mesh_deformation_external_01.prm b/tests/mesh_deformation_external_01.prm new file mode 100644 index 00000000000..583e11c7536 --- /dev/null +++ b/tests/mesh_deformation_external_01.prm @@ -0,0 +1,96 @@ +# Test the ExternalToolInterface class for mesh deformation +# derived class written as a plugin, 2d box geometry, initial +# topography coming from the plugin (and not initial topography) +# +# MPI: 2 + +set Dimension = 2 +set Use years in output instead of seconds = false +set End time = 0.004 +set Maximum time step = 0.0005 +set Nonlinear solver scheme = no Advection, no Stokes +set Pressure normalization = surface +set Surface pressure = 0 + +subsection Geometry model + set Model name = box + + subsection Box + set X extent = 1 + set Y extent = 1 + end +end + +# Temperature effects are ignored +subsection Initial temperature model + set Model name = function + + subsection Function + set Function expression = 100+x + end +end + +subsection Boundary temperature model + set Fixed temperature boundary indicators = bottom, top + set List of model names = initial temperature +end + +# Free slip on all boundaries +subsection Boundary velocity model + set Tangential velocity boundary indicators = left, right, bottom, top +end + +subsection Mesh deformation + set Mesh deformation boundary indicators = top: external deformation + set Additional tangential mesh velocity boundary indicators = left, right +end + +# Vertical gravity +subsection Gravity model + set Model name = vertical + + subsection Vertical + set Magnitude = 1 + end +end + +# One material with unity properties +subsection Material model + set Model name = simple + + subsection Simple model + set Reference density = 1 + set Reference specific heat = 1 + set Reference temperature = 0 + set Thermal conductivity = 1 + set Thermal expansion coefficient = 1 + set Viscosity = 1 + end +end + +# We also have to specify that we want to use the Boussinesq +# approximation (assuming the density in the temperature +# equation to be constant, and incompressibility). +subsection Formulation + set Formulation = Boussinesq approximation +end + +# We here use a globally refined mesh without +# adaptive mesh refinement. +subsection Mesh refinement + set Initial global refinement = 4 + set Initial adaptive refinement = 0 + set Time steps between mesh refinement = 0 +end + +subsection Postprocess + set List of postprocessors = visualization, topography + + subsection Visualization + set Time between graphical output = 0 + set Output mesh velocity = true + set Output mesh displacement = true + set Output undeformed mesh = false + set Interpolate output = false + end +end diff --git a/tests/mesh_deformation_external_01/screen-output b/tests/mesh_deformation_external_01/screen-output new file mode 100644 index 00000000000..a62ec02d987 --- /dev/null +++ b/tests/mesh_deformation_external_01/screen-output @@ -0,0 +1,159 @@ +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- + +Loading shared library <./libmesh_deformation_external_01.debug.so> + +initialize() +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- +Number of active cells: 256 (on 5 levels) +Number of degrees of freedom: 3,556 (2,178+289+1,089) + +Number of mesh deformation degrees of freedom: 578 + Solving mesh displacement system... 6 iterations. +*** Timestep 0: t=0 seconds, dt=0 seconds +update() + setting points +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.3 1 0 0 0 100.3 +0.5 1 0 0 0 100.5 +rank 1: +0.72 1 0 0 0 100.72 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01/solution/solution-00000 + Topography min/max: 0 m, 0.1 m + +*** Timestep 1: t=0.0005 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.3 1 0 0 0 100.3 +0.5 1 0 0 0 100.5 +rank 1: +0.72 1 0 0 0 100.72 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01/solution/solution-00001 + Topography min/max: -0.00175 m, 0.1 m + +*** Timestep 2: t=0.001 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.3 1 0 0 0 100.3 +0.5 1 0 0 0 100.5 +rank 1: +0.72 1 0 0 0 100.72 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01/solution/solution-00002 + Topography min/max: -0.0035 m, 0.1 m + +*** Timestep 3: t=0.0015 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.3 1 0 0 0 100.3 +0.5 1 0 0 0 100.5 +rank 1: +0.72 1 0 0 0 100.72 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01/solution/solution-00003 + Topography min/max: -0.00525 m, 0.1 m + +*** Timestep 4: t=0.002 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.3 1 0 0 0 100.3 +0.5 1 0 0 0 100.5 +rank 1: +0.72 1 0 0 0 100.72 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01/solution/solution-00004 + Topography min/max: -0.007 m, 0.1 m + +*** Timestep 5: t=0.0025 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.3 1 0 0 0 100.3 +0.5 1 0 0 0 100.5 +rank 1: +0.72 1 0 0 0 100.72 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01/solution/solution-00005 + Topography min/max: -0.00875 m, 0.1 m + +*** Timestep 6: t=0.003 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.3 1 0 0 0 100.3 +0.5 1 0 0 0 100.5 +rank 1: +0.72 1 0 0 0 100.72 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01/solution/solution-00006 + Topography min/max: -0.0105 m, 0.1 m + +*** Timestep 7: t=0.0035 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.3 1 0 0 0 100.3 +0.5 1 0 0 0 100.5 +rank 1: +0.72 1 0 0 0 100.72 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01/solution/solution-00007 + Topography min/max: -0.01225 m, 0.105 m + +*** Timestep 8: t=0.004 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.3 1 0 0 0 100.3 +0.5 1 0 0 0 100.5 +rank 1: +0.72 1 0 0 0 100.72 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01/solution/solution-00008 + Topography min/max: -0.014 m, 0.12 m + +Termination requested by criterion: end time + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- diff --git a/tests/mesh_deformation_external_01/statistics b/tests/mesh_deformation_external_01/statistics new file mode 100644 index 00000000000..9ee80130b94 --- /dev/null +++ b/tests/mesh_deformation_external_01/statistics @@ -0,0 +1,19 @@ +# 1: Time step number +# 2: Time (seconds) +# 3: Time step size (seconds) +# 4: Number of mesh cells +# 5: Number of Stokes degrees of freedom +# 6: Number of temperature degrees of freedom +# 7: Number of nonlinear iterations +# 8: Visualization file name +# 9: Minimum topography (m) +# 10: Maximum topography (m) +0 0.000000000000e+00 0.000000000000e+00 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00000 0.00000000e+00 1.00000000e-01 +1 5.000000000000e-04 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00001 -1.75000000e-03 1.00000000e-01 +2 1.000000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00002 -3.50000000e-03 1.00000000e-01 +3 1.500000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00003 -5.25000000e-03 1.00000000e-01 +4 2.000000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00004 -7.00000000e-03 1.00000000e-01 +5 2.500000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00005 -8.75000000e-03 1.00000000e-01 +6 3.000000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00006 -1.05000000e-02 1.00000000e-01 +7 3.500000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00007 -1.22500000e-02 1.05000000e-01 +8 4.000000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00008 -1.40000000e-02 1.20000000e-01 diff --git a/tests/mesh_deformation_external_01_3d.cc b/tests/mesh_deformation_external_01_3d.cc new file mode 100644 index 00000000000..23ce928a720 --- /dev/null +++ b/tests/mesh_deformation_external_01_3d.cc @@ -0,0 +1,21 @@ +/* + Copyright (C) 2025 - 2025 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + +#include "mesh_deformation_external_01.cc" diff --git a/tests/mesh_deformation_external_01_3d.prm b/tests/mesh_deformation_external_01_3d.prm new file mode 100644 index 00000000000..fe4a8002c23 --- /dev/null +++ b/tests/mesh_deformation_external_01_3d.prm @@ -0,0 +1,97 @@ +# Test the ExternalToolInterface class for mesh deformation +# derived class written as a plugin, 3d box geometry, initial +# topography coming from the plugin (and not initial topography) +# +# MPI: 2 + +set Dimension = 3 +set Use years in output instead of seconds = false +set End time = 0.001 +set Maximum time step = 0.0005 +set Nonlinear solver scheme = no Advection, no Stokes +set Pressure normalization = surface +set Surface pressure = 0 + +subsection Geometry model + set Model name = box + + subsection Box + set X extent = 1 + set Y extent = 1 + set Z extent = 1 + end +end + +# Temperature effects are ignored +subsection Initial temperature model + set Model name = function + + subsection Function + set Function expression = 100*x+y + end +end + +subsection Boundary temperature model + set Fixed temperature boundary indicators = bottom, top + set List of model names = initial temperature +end + +# Free slip on all boundaries +subsection Boundary velocity model + set Tangential velocity boundary indicators = left, right, bottom, top +end + +subsection Mesh deformation + set Mesh deformation boundary indicators = top: external deformation + set Additional tangential mesh velocity boundary indicators = left, right +end + +# Vertical gravity +subsection Gravity model + set Model name = vertical + + subsection Vertical + set Magnitude = 1 + end +end + +# One material with unity properties +subsection Material model + set Model name = simple + + subsection Simple model + set Reference density = 1 + set Reference specific heat = 1 + set Reference temperature = 0 + set Thermal conductivity = 1 + set Thermal expansion coefficient = 1 + set Viscosity = 1 + end +end + +# We also have to specify that we want to use the Boussinesq +# approximation (assuming the density in the temperature +# equation to be constant, and incompressibility). +subsection Formulation + set Formulation = Boussinesq approximation +end + +# We here use a globally refined mesh without +# adaptive mesh refinement. +subsection Mesh refinement + set Initial global refinement = 2 + set Initial adaptive refinement = 0 + set Time steps between mesh refinement = 0 +end + +subsection Postprocess + set List of postprocessors = visualization, topography + + subsection Visualization + set Time between graphical output = 0 + set Output mesh velocity = true + set Output mesh displacement = true + set Output undeformed mesh = false + set Interpolate output = false + end +end diff --git a/tests/mesh_deformation_external_01_3d/screen-output b/tests/mesh_deformation_external_01_3d/screen-output new file mode 100644 index 00000000000..7effc939cde --- /dev/null +++ b/tests/mesh_deformation_external_01_3d/screen-output @@ -0,0 +1,75 @@ +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- + +Loading shared library <./libmesh_deformation_external_01_3d.debug.so> + +initialize() +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- +Number of active cells: 64 (on 3 levels) +Number of degrees of freedom: 3,041 (2,187+125+729) + +Number of mesh deformation degrees of freedom: 375 + Solving mesh displacement system... 6 iterations. +*** Timestep 0: t=0 seconds, dt=0 seconds +update() + setting points +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.3 0.7 1 0 0 0 -6.70092e-18 30.7 +0.5 0.5 1 0 0 0 -6.70092e-18 50.5 +0.7 0.2 1 0 0 0 -6.70092e-18 70.2 +0.7 0.91 1 0 0 0 -6.70092e-18 70.91 +rank 1: +0.9 0.9 1 0 0 0 -6.70092e-18 90.9 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01_3d/solution/solution-00000 + Topography min/max: 0 m, 0.1 m + +*** Timestep 1: t=0.0005 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.3 0.7 1 0 0 0 -6.70092e-18 30.7 +0.5 0.5 1 0 0 0 -6.70092e-18 50.5 +0.7 0.2 1 0 0 0 -6.70092e-18 70.2 +0.7 0.91 1 0 0 0 -6.70092e-18 70.91 +rank 1: +0.9 0.9 1 0 0 0 -6.70092e-18 90.9 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01_3d/solution/solution-00001 + Topography min/max: -0.00175 m, 0.1 m + +*** Timestep 2: t=0.001 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.3 0.7 1 0 0 0 -6.70092e-18 30.7 +0.5 0.5 1 0 0 0 -6.70092e-18 50.5 +0.7 0.2 1 0 0 0 -6.70092e-18 70.2 +0.7 0.91 1 0 0 0 -6.70092e-18 70.91 +rank 1: +0.9 0.9 1 0 0 0 -6.70092e-18 90.9 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01_3d/solution/solution-00002 + Topography min/max: -0.0035 m, 0.1 m + +Termination requested by criterion: end time + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- diff --git a/tests/mesh_deformation_external_01_3d/statistics b/tests/mesh_deformation_external_01_3d/statistics new file mode 100644 index 00000000000..5ad43574272 --- /dev/null +++ b/tests/mesh_deformation_external_01_3d/statistics @@ -0,0 +1,13 @@ +# 1: Time step number +# 2: Time (seconds) +# 3: Time step size (seconds) +# 4: Number of mesh cells +# 5: Number of Stokes degrees of freedom +# 6: Number of temperature degrees of freedom +# 7: Number of nonlinear iterations +# 8: Visualization file name +# 9: Minimum topography (m) +# 10: Maximum topography (m) +0 0.000000000000e+00 0.000000000000e+00 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00000 0.00000000e+00 1.00000000e-01 +1 5.000000000000e-04 5.000000000000e-04 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00001 -1.75000000e-03 1.00000000e-01 +2 1.000000000000e-03 5.000000000000e-04 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00002 -3.50000000e-03 1.00000000e-01 From b1f32adae5ec9392f54ea0928728ebba9cb61cc7 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Sat, 13 Sep 2025 17:38:23 -0400 Subject: [PATCH 06/10] address comments --- .../external_tool_interface.h | 17 ++-- .../external_tool_interface.cc | 93 +++++++++++-------- 2 files changed, 63 insertions(+), 47 deletions(-) diff --git a/include/aspect/mesh_deformation/external_tool_interface.h b/include/aspect/mesh_deformation/external_tool_interface.h index 83dceec22c9..4d95d0c19d2 100644 --- a/include/aspect/mesh_deformation/external_tool_interface.h +++ b/include/aspect/mesh_deformation/external_tool_interface.h @@ -166,12 +166,15 @@ namespace aspect evaluate_aspect_variables_at_points () const; /** - * Given velocities (typically vertical, but not always) at - * all of the points this process has previously set via - * set_evaluation_points(), compute a (global) finite element - * field vector that in the velocity components of surface - * nodes corresponds to an interpolation of the velocities - * provided. + * Interpolate from velocities given in the evaluation points + * to ASPECT velocities on the surface in form of a finite + * element field. + * + * The @p velocities, which are typically vertical, were previously + * computed by the external tool and belong to the current process. + * The Output is a (global) finite element field vector that in + * the velocity components of surface nodes corresponds to an + * interpolation of the velocities provided. * * The output of this function can then be used as input for a * function that implements the @@ -179,7 +182,7 @@ namespace aspect * base class. */ LinearAlgebra::Vector - interpolate_velocities_to_surface_points (const std::vector> &vertical_velocities) const; + interpolate_external_velocities_to_surface_support_points (const std::vector> &velocities) const; protected: diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc index c37b67e0414..e563fe9b5be 100644 --- a/source/mesh_deformation/external_tool_interface.cc +++ b/source/mesh_deformation/external_tool_interface.cc @@ -37,7 +37,9 @@ namespace aspect { - + /** + * Produce graphical output defined in points provided by the user. + */ template class PointDataOut : public dealii::DataOutInterface<0, spacedim> { @@ -72,14 +74,21 @@ namespace aspect DataComponentInterpretation::DataComponentInterpretation> &data_component_interpretations_ = {}) { - Assert( - data_component_names.size() == data_component_interpretations_.size(), - ExcMessage( - "When calling Particles::DataOut::build_patches with data component " - "names and interpretations you need to provide as many data component " - "names as interpretations. Provide the same name for components that " - "belong to a single vector or tensor.")); - + Assert(data_component_names.size() == data_component_interpretations_.size(), + ExcMessage( + "When calling PointDataOut::build_patches() with data component " + "names and interpretations you need to provide as many data component " + "names as interpretations. Provide the same name for components that " + "belong to a single vector or tensor.")); + + Assert(data.size() == 0 || locations.size(), + ExcMessage("You need to either provide no data or data for each point.")); + for (const auto &datum : data) + Assert(datum.size() == data_component_names.size(), + ExcMessage("The data provided in each point needs to have the same number " + "of components as names were provided.")); + + // Prepend the "id" to the data fields provided by the user: dataset_names.clear(); dataset_names.emplace_back("id"); dataset_names.insert(dataset_names.end(), @@ -104,26 +113,17 @@ namespace aspect patches[i].vertices[0] = locations[i]; patches[i].patch_index = i; - // We have one more data components than dataset_names (the particle id) + // Store id and properties given by the user: patches[i].data.reinit(n_data_components, 1); - - patches[i].data(0, 0) = i; // id - - if (n_data_components > 1) - { - for (unsigned int property_index = 0; - property_index < n_property_components; - ++property_index) - patches[i].data(property_index + 1, 0) = - data[i][property_index]; - } + patches[i].data(0, 0) = i; // store id + for (unsigned int property_index = 0; property_index < n_property_components; ++property_index) + patches[i].data(property_index + 1, 0) = data[i][property_index]; } } protected: /** - * Returns the patches built by the data_out class which was previously - * built using a particle handler + * Returns the patches previously built by the build_patches() function. */ virtual const std::vector> & get_patches() const override @@ -318,7 +318,7 @@ namespace aspect = compute_updated_velocities_at_points(aspect_surface_velocities); const LinearAlgebra::Vector v_interpolated - = interpolate_velocities_to_surface_points(external_surface_velocities); + = interpolate_external_velocities_to_surface_support_points(external_surface_velocities); // TODO: need ghost values of v_interpolated? @@ -363,10 +363,12 @@ namespace aspect // for later function calls describe the same number of points. this->evaluation_points = evaluation_points; - // Set up RemotePointEvaluation with a Mapping of the undeformed mesh: - remote_point_evaluator = std::make_unique>(); - // TODO: does this need to be a higher order mapping for spherical problems? + // Set up RemotePointEvaluation. The evaluation points are given in reference coordinates, + // so we need to use a simple mapping instead of the one stored in the Simulator. The latter + // would produce the deformed mesh. We currently always use a Q1 mapping when mesh deformation + // is enabled, so a Q1 mapping is the right choice. static MappingQ mapping(1); + remote_point_evaluator = std::make_unique>(); remote_point_evaluator->reinit(this->evaluation_points, this->get_triangulation(), mapping); // Finally, also ensure that upon mesh refinement, all of the @@ -401,13 +403,13 @@ namespace aspect // to put directly into deal.II... for (unsigned int c=0; c values = VectorTools::point_values<1>(*this->remote_point_evaluator, - this->get_dof_handler(), - this->get_solution(), - dealii::VectorTools::EvaluationFlags::avg, - c); + const std::vector values = VectorTools::point_values<1>(*this->remote_point_evaluator, + this->get_dof_handler(), + this->get_solution(), + dealii::VectorTools::EvaluationFlags::avg, + c); for (unsigned int i=0; i LinearAlgebra::Vector ExternalToolInterface:: - interpolate_velocities_to_surface_points (const std::vector> &velocities) const + interpolate_external_velocities_to_surface_support_points (const std::vector> &velocities) const { Assert (remote_point_evaluator != nullptr, ExcMessage("You can only call this function if you have previously " @@ -458,16 +460,24 @@ namespace aspect const DoFHandler &mesh_dof_handler = this->get_mesh_deformation_handler().get_mesh_deformation_dof_handler(); LinearAlgebra::Vector vector_with_surface_velocities(mesh_dof_handler.locally_owned_dofs(), this->get_mpi_communicator()); + + // The remote_point_evaluator will gives us the velocities in all evaluation points that are within one of our locally + // owned cells. The lambda defined below receives a list of points and their velocities for each cell. The coordinates + // are given in coordinates of the unit cell. + // For each support point of the velocity DoFHandler, we will set the velocity from the closest evaluation point. We + // do this by keeping track of 1/distance of the closest evaluation point checked so far. The initial value of 0.0 + // denotes an infinite distance. LinearAlgebra::Vector one_over_distance_vec(mesh_dof_handler.locally_owned_dofs(), this->get_mpi_communicator()); - const unsigned int dofs_per_cell = mesh_dof_handler.get_fe().dofs_per_cell; - std::vector cell_dof_indices (dofs_per_cell); + // Note: We assume that process_and_evaluate() does not call our lambda concurrently, otherwise we would have write + // conflicts when updating vector_with_surface_velocities and one_over_distance_vec. const auto eval_func = [&](const ArrayView< const Tensor<1,dim>> &values, const typename Utilities::MPI::RemotePointEvaluation::CellData &cell_data) { + std::vector cell_dof_indices (dofs_per_cell); for (const auto cell_index : cell_data.cell_indices()) { const auto cell_dofs = @@ -478,6 +488,9 @@ namespace aspect const auto local_values = cell_data.get_data_view(cell_index, values); cell_dofs->get_dof_indices(cell_dof_indices); + + // Note: This search is a nested loop with the inner part executed #evaluation_point_in_this_cell * #dofs_per_cell + // times. We could precompute this information as the point locations and do not change (outside of mesh refinement). const std::vector< Point< dim >> &support_points = mesh_dof_handler.get_fe().get_unit_support_points(); for (unsigned int i=0; i one_over_distance_vec(cell_dof_indices[j])) { + // The point i is closer to support point j than anything we have seen so far. Keep track + // of the distance and write the correct velocity component into the result: one_over_distance_vec(cell_dof_indices[j]) = one_over_distance; - const unsigned int c = mesh_dof_handler.get_fe().system_to_component_index(j).first; - vector_with_surface_velocities(cell_dof_indices[j]) = local_values[i][c]; + const unsigned int component = mesh_dof_handler.get_fe().system_to_component_index(j).first; + vector_with_surface_velocities(cell_dof_indices[j]) = local_values[i][component]; } } } @@ -496,8 +511,6 @@ namespace aspect }; this->remote_point_evaluator->template process_and_evaluate,1>(velocities, eval_func, /*sort_data*/ true); - - one_over_distance_vec.compress(VectorOperation::insert); vector_with_surface_velocities.compress(VectorOperation::insert); return vector_with_surface_velocities; From 6e7884e56c843cb2d7f698b2080e4647fc0d8305 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Sat, 13 Dec 2025 19:14:33 -0500 Subject: [PATCH 07/10] sequential version --- .../external_tool_interface.h | 16 + .../external_tool_interface.cc | 138 ++++-- tests/mesh_deformation_external_01.cc | 30 +- tests/mesh_deformation_external_01.prm | 6 +- .../screen-output | 194 ++++++-- tests/mesh_deformation_external_01/statistics | 18 +- tests/mesh_deformation_external_01_3d.prm | 6 +- .../screen-output | 446 +++++++++++++++++- .../statistics | 12 +- 9 files changed, 721 insertions(+), 145 deletions(-) diff --git a/include/aspect/mesh_deformation/external_tool_interface.h b/include/aspect/mesh_deformation/external_tool_interface.h index 4d95d0c19d2..656e39d946a 100644 --- a/include/aspect/mesh_deformation/external_tool_interface.h +++ b/include/aspect/mesh_deformation/external_tool_interface.h @@ -189,6 +189,22 @@ namespace aspect std::vector> evaluation_points; std::unique_ptr> remote_point_evaluator; + + /** + * A struct to map between DoF indices and evaluation points + */ + struct dof_to_eval_point_data + { + types::global_dof_index dof_index; + unsigned int evaluation_point_index; + unsigned int component; + }; + + /** + * A vector to store a map between Dof indices and evaluation points + */ + std::vector map_dof_to_eval_point; + }; } } diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc index e563fe9b5be..52d76877ff2 100644 --- a/source/mesh_deformation/external_tool_interface.cc +++ b/source/mesh_deformation/external_tool_interface.cc @@ -371,6 +371,91 @@ namespace aspect remote_point_evaluator = std::make_unique>(); remote_point_evaluator->reinit(this->evaluation_points, this->get_triangulation(), mapping); + + // Create a mapping from evaluation points to support points. Note that one evaluation point can map to + // multiple support points. + { + // For now, we only support a single MPI rank for ASPECT and the external mesh deformation class. Because deciding + // which evaluation point is closest to a support point requires somewhat complex MPI communication. + Assert(Utilities::MPI::n_mpi_processes(this->get_mpi_communicator()) == 1, + ExcNotImplemented("This function is not implemented for parallel computations.")); + + const DoFHandler &mesh_dof_handler = this->get_mesh_deformation_handler().get_mesh_deformation_dof_handler(); + + std::vector squared_distances(mesh_dof_handler.locally_owned_dofs().size(), std::numeric_limits::max()); + std::vector closest_evaluation_point_and_component(mesh_dof_handler.locally_owned_dofs().size(), + dof_to_eval_point_data {numbers::invalid_dof_index, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}); + + // TODO: do we need to support the case of more than one different mesh deformation plugin to be active? + const auto boundary_ids = this->get_mesh_deformation_boundary_indicators(); + + IndexSet boundary_dofs = DoFTools::extract_boundary_dofs(mesh_dof_handler, ComponentMask(dim, true), boundary_ids); + + const unsigned int dofs_per_cell = mesh_dof_handler.get_fe().dofs_per_cell; + std::vector local_dof_indices (dofs_per_cell); + + // The remote_point_evaluator will gives us the velocities in all evaluation points that are within one of our locally + // owned cells. The lambda defined below receives a list of points and their velocities for each cell. The coordinates + // are given in coordinates of the unit cell. + // For each support point of the velocity DoFHandler, we will try to find the closest evaluation point. We + // do this by keeping track of the squared distance of the closest evaluation point checked so far. + + + // Note: We assume that process_and_evaluate() does not call our lambda concurrently, otherwise we would have write + // conflicts when updating vector_with_surface_velocities and one_over_distance_vec. + + const auto eval_func = [&](const ArrayView &values, + const typename Utilities::MPI::RemotePointEvaluation::CellData &cell_data) + { + std::vector cell_dof_indices (dofs_per_cell); + for (const auto cell_index : cell_data.cell_indices()) + { + const auto cell_dofs = + cell_data.get_active_cell_iterator(cell_index)->as_dof_handler_iterator( + mesh_dof_handler); + cell_dofs->get_dof_indices(cell_dof_indices); + + const ArrayView> unit_points = cell_data.get_unit_points(cell_index); + const auto local_values = cell_data.get_data_view(cell_index, values); + + const std::vector< Point< dim >> &support_points = mesh_dof_handler.get_fe().get_unit_support_points(); + for (unsigned int i=0; i indices (evaluation_points.size()); + std::iota(indices.begin(), indices.end(), 0); + this->remote_point_evaluator->template process_and_evaluate(indices, eval_func, /*sort_data*/ true); + + map_dof_to_eval_point.clear(); + for (unsigned int i=0; iget_signals().pre_refinement_store_user_data @@ -434,7 +519,7 @@ namespace aspect static unsigned int output_no = 0; PointDataOut out; - // const auto &mapping = this->get_mapping(); + //const auto &mapping = this->get_mapping(); std::vector> real_evaluation_points(evaluation_points.size()); std::vector> data(evaluation_points.size(), std::vector(dim, 0.0)); for (unsigned int i=0; iget_mpi_communicator()); - // The remote_point_evaluator will gives us the velocities in all evaluation points that are within one of our locally - // owned cells. The lambda defined below receives a list of points and their velocities for each cell. The coordinates - // are given in coordinates of the unit cell. - // For each support point of the velocity DoFHandler, we will set the velocity from the closest evaluation point. We - // do this by keeping track of 1/distance of the closest evaluation point checked so far. The initial value of 0.0 - // denotes an infinite distance. - LinearAlgebra::Vector one_over_distance_vec(mesh_dof_handler.locally_owned_dofs(), - this->get_mpi_communicator()); - - const unsigned int dofs_per_cell = mesh_dof_handler.get_fe().dofs_per_cell; - - // Note: We assume that process_and_evaluate() does not call our lambda concurrently, otherwise we would have write - // conflicts when updating vector_with_surface_velocities and one_over_distance_vec. - const auto eval_func = [&](const ArrayView< const Tensor<1,dim>> &values, - const typename Utilities::MPI::RemotePointEvaluation::CellData &cell_data) - { - std::vector cell_dof_indices (dofs_per_cell); - for (const auto cell_index : cell_data.cell_indices()) - { - const auto cell_dofs = - cell_data.get_active_cell_iterator(cell_index)->as_dof_handler_iterator( - mesh_dof_handler); - - const ArrayView> unit_points = cell_data.get_unit_points(cell_index); - const auto local_values = cell_data.get_data_view(cell_index, values); - - cell_dofs->get_dof_indices(cell_dof_indices); - - // Note: This search is a nested loop with the inner part executed #evaluation_point_in_this_cell * #dofs_per_cell - // times. We could precompute this information as the point locations and do not change (outside of mesh refinement). - const std::vector< Point< dim >> &support_points = mesh_dof_handler.get_fe().get_unit_support_points(); - for (unsigned int i=0; i one_over_distance_vec(cell_dof_indices[j])) - { - // The point i is closer to support point j than anything we have seen so far. Keep track - // of the distance and write the correct velocity component into the result: - one_over_distance_vec(cell_dof_indices[j]) = one_over_distance; - const unsigned int component = mesh_dof_handler.get_fe().system_to_component_index(j).first; - vector_with_surface_velocities(cell_dof_indices[j]) = local_values[i][component]; - } - } - } - } - }; + for (const auto &entry : map_dof_to_eval_point) + vector_with_surface_velocities[entry.dof_index] = velocities[entry.evaluation_point_index][entry.component]; - this->remote_point_evaluator->template process_and_evaluate,1>(velocities, eval_func, /*sort_data*/ true); vector_with_surface_velocities.compress(VectorOperation::insert); return vector_with_surface_velocities; diff --git a/tests/mesh_deformation_external_01.cc b/tests/mesh_deformation_external_01.cc index ed803861046..79a8d92b4ba 100644 --- a/tests/mesh_deformation_external_01.cc +++ b/tests/mesh_deformation_external_01.cc @@ -60,28 +60,19 @@ namespace aspect { if constexpr (dim == 2) { - points.emplace_back(0.3, 1.0); - points.emplace_back(0.5, 1.0); + for (unsigned int x=0; x<11; ++x) + points.emplace_back(0.05+x*1.0/11, 1.0); } else { - points.emplace_back(0.3, 0.7, 1.0); - points.emplace_back(0.5, 0.5, 1.0); - points.emplace_back(0.7, 0.2, 1.0); - points.emplace_back(0.7, 0.91, 1.0); + for (unsigned int x=0; x<5; ++x) + for (unsigned int y=0; y<5; ++y) + points.emplace_back(0.1+x*1.0/5, 0.1+y*1.0/5, 1.0); } } else - { - if constexpr (dim == 2) - { - points.emplace_back(0.72, 1.0); - } - else - { - points.emplace_back(0.9, 0.9, 1.0); - } - } + Assert(false, ExcNotImplemented()); + this->set_evaluation_points (points); } } @@ -123,10 +114,11 @@ namespace aspect Assert(current_solution_at_points.size() == this->evaluation_points.size(), ExcInternalError()); std::vector> velocities(current_solution_at_points.size(), Tensor<1,dim>()); - if (velocities.size()>1) + if (velocities.size()>6) { - velocities[0][dim-1]=30.0; - velocities[1][dim-1]=-3.5; + velocities[1][dim-1]=30.0; + velocities[4][dim-1]=-5.0; + velocities[6][dim-1]=10.0; } return velocities; } diff --git a/tests/mesh_deformation_external_01.prm b/tests/mesh_deformation_external_01.prm index 583e11c7536..9cba69cd6bc 100644 --- a/tests/mesh_deformation_external_01.prm +++ b/tests/mesh_deformation_external_01.prm @@ -2,10 +2,10 @@ # derived class written as a plugin, 2d box geometry, initial # topography coming from the plugin (and not initial topography) # -# MPI: 2 +# MPI: 1 set Dimension = 2 -set Use years in output instead of seconds = false +set Use years instead of seconds = false set End time = 0.004 set Maximum time step = 0.0005 set Nonlinear solver scheme = no Advection, no Stokes @@ -78,7 +78,7 @@ end # We here use a globally refined mesh without # adaptive mesh refinement. subsection Mesh refinement - set Initial global refinement = 4 + set Initial global refinement = 3 set Initial adaptive refinement = 0 set Time steps between mesh refinement = 0 end diff --git a/tests/mesh_deformation_external_01/screen-output b/tests/mesh_deformation_external_01/screen-output index a62ec02d987..a0a7222240c 100644 --- a/tests/mesh_deformation_external_01/screen-output +++ b/tests/mesh_deformation_external_01/screen-output @@ -7,21 +7,65 @@ Loading shared library <./libmesh_deformation_external_01.debug.so> initialize() ----------------------------------------------------------------------------- ----------------------------------------------------------------------------- -Number of active cells: 256 (on 5 levels) -Number of degrees of freedom: 3,556 (2,178+289+1,089) +Number of active cells: 64 (on 4 levels) +Number of degrees of freedom: 948 (578+81+289) -Number of mesh deformation degrees of freedom: 578 +Number of mesh deformation degrees of freedom: 162 Solving mesh displacement system... 6 iterations. *** Timestep 0: t=0 seconds, dt=0 seconds update() setting points +map_dof_to_eval_point (dof, evaluation_point_index, component): +110 0 0 +111 0 1 +112 1 0 +113 1 1 +114 2 0 +115 2 1 +116 0 0 +117 0 1 +118 1 0 +119 1 1 +120 2 0 +121 2 1 +122 4 0 +123 4 1 +124 5 0 +125 5 1 +126 4 0 +127 4 1 +128 5 0 +129 5 1 +146 6 0 +147 6 1 +148 8 0 +149 8 1 +150 6 0 +151 6 1 +152 8 0 +153 8 1 +154 9 0 +155 9 1 +156 10 0 +157 10 1 +158 9 0 +159 9 1 +160 10 0 +161 10 1 compute_updated_velocities_at_points() Solution at evaluation points: rank 0: -0.3 1 0 0 0 100.3 -0.5 1 0 0 0 100.5 -rank 1: -0.72 1 0 0 0 100.72 +0.05 1 0 0 0 100.05 +0.140909 1 0 0 0 100.141 +0.231818 1 0 0 0 100.232 +0.322727 1 0 0 0 100.323 +0.413636 1 0 0 0 100.414 +0.504545 1 0 0 0 100.505 +0.595455 1 0 0 0 100.595 +0.686364 1 0 0 0 100.686 +0.777273 1 0 0 0 100.777 +0.868182 1 0 0 0 100.868 +0.959091 1 0 0 0 100.959 Solving mesh displacement system... 4 iterations. Postprocessing: @@ -33,120 +77,176 @@ update() compute_updated_velocities_at_points() Solution at evaluation points: rank 0: -0.3 1 0 0 0 100.3 -0.5 1 0 0 0 100.5 -rank 1: -0.72 1 0 0 0 100.72 +0.05 1 0 0 0 100.05 +0.140909 1 0 0 0 100.141 +0.231818 1 0 0 0 100.232 +0.322727 1 0 0 0 100.323 +0.413636 1 0 0 0 100.414 +0.504545 1 0 0 0 100.505 +0.595455 1 0 0 0 100.595 +0.686364 1 0 0 0 100.686 +0.777273 1 0 0 0 100.777 +0.868182 1 0 0 0 100.868 +0.959091 1 0 0 0 100.959 Solving mesh displacement system... 4 iterations. Postprocessing: Writing graphical output: output-mesh_deformation_external_01/solution/solution-00001 - Topography min/max: -0.00175 m, 0.1 m + Topography min/max: -0.0025 m, 0.1 m *** Timestep 2: t=0.001 seconds, dt=0.0005 seconds update() compute_updated_velocities_at_points() Solution at evaluation points: rank 0: -0.3 1 0 0 0 100.3 -0.5 1 0 0 0 100.5 -rank 1: -0.72 1 0 0 0 100.72 +0.05 1 0 0 0 100.05 +0.140909 1 0 0 0 100.141 +0.231818 1 0 0 0 100.232 +0.322727 1 0 0 0 100.323 +0.413636 1 0 0 0 100.414 +0.504545 1 0 0 0 100.505 +0.595455 1 0 0 0 100.595 +0.686364 1 0 0 0 100.686 +0.777273 1 0 0 0 100.777 +0.868182 1 0 0 0 100.868 +0.959091 1 0 0 0 100.959 Solving mesh displacement system... 4 iterations. Postprocessing: Writing graphical output: output-mesh_deformation_external_01/solution/solution-00002 - Topography min/max: -0.0035 m, 0.1 m + Topography min/max: -0.005 m, 0.1 m *** Timestep 3: t=0.0015 seconds, dt=0.0005 seconds update() compute_updated_velocities_at_points() Solution at evaluation points: rank 0: -0.3 1 0 0 0 100.3 -0.5 1 0 0 0 100.5 -rank 1: -0.72 1 0 0 0 100.72 +0.05 1 0 0 0 100.05 +0.140909 1 0 0 0 100.141 +0.231818 1 0 0 0 100.232 +0.322727 1 0 0 0 100.323 +0.413636 1 0 0 0 100.414 +0.504545 1 0 0 0 100.505 +0.595455 1 0 0 0 100.595 +0.686364 1 0 0 0 100.686 +0.777273 1 0 0 0 100.777 +0.868182 1 0 0 0 100.868 +0.959091 1 0 0 0 100.959 Solving mesh displacement system... 4 iterations. Postprocessing: Writing graphical output: output-mesh_deformation_external_01/solution/solution-00003 - Topography min/max: -0.00525 m, 0.1 m + Topography min/max: -0.0075 m, 0.1 m *** Timestep 4: t=0.002 seconds, dt=0.0005 seconds update() compute_updated_velocities_at_points() Solution at evaluation points: rank 0: -0.3 1 0 0 0 100.3 -0.5 1 0 0 0 100.5 -rank 1: -0.72 1 0 0 0 100.72 +0.05 1 0 0 0 100.05 +0.140909 1 0 0 0 100.141 +0.231818 1 0 0 0 100.232 +0.322727 1 0 0 0 100.323 +0.413636 1 0 0 0 100.414 +0.504545 1 0 0 0 100.505 +0.595455 1 0 0 0 100.595 +0.686364 1 0 0 0 100.686 +0.777273 1 0 0 0 100.777 +0.868182 1 0 0 0 100.868 +0.959091 1 0 0 0 100.959 Solving mesh displacement system... 4 iterations. Postprocessing: Writing graphical output: output-mesh_deformation_external_01/solution/solution-00004 - Topography min/max: -0.007 m, 0.1 m + Topography min/max: -0.01 m, 0.1 m *** Timestep 5: t=0.0025 seconds, dt=0.0005 seconds update() compute_updated_velocities_at_points() Solution at evaluation points: rank 0: -0.3 1 0 0 0 100.3 -0.5 1 0 0 0 100.5 -rank 1: -0.72 1 0 0 0 100.72 +0.05 1 0 0 0 100.05 +0.140909 1 0 0 0 100.141 +0.231818 1 0 0 0 100.232 +0.322727 1 0 0 0 100.323 +0.413636 1 0 0 0 100.414 +0.504545 1 0 0 0 100.505 +0.595455 1 0 0 0 100.595 +0.686364 1 0 0 0 100.686 +0.777273 1 0 0 0 100.777 +0.868182 1 0 0 0 100.868 +0.959091 1 0 0 0 100.959 Solving mesh displacement system... 4 iterations. Postprocessing: Writing graphical output: output-mesh_deformation_external_01/solution/solution-00005 - Topography min/max: -0.00875 m, 0.1 m + Topography min/max: -0.0125 m, 0.1 m *** Timestep 6: t=0.003 seconds, dt=0.0005 seconds update() compute_updated_velocities_at_points() Solution at evaluation points: rank 0: -0.3 1 0 0 0 100.3 -0.5 1 0 0 0 100.5 -rank 1: -0.72 1 0 0 0 100.72 +0.05 1 0 0 0 100.05 +0.140909 1 0 0 0 100.141 +0.231818 1 0 0 0 100.232 +0.322727 1 0 0 0 100.323 +0.413636 1 0 0 0 100.414 +0.504545 1 0 0 0 100.505 +0.595455 1 0 0 0 100.595 +0.686364 1 0 0 0 100.686 +0.777273 1 0 0 0 100.777 +0.868182 1 0 0 0 100.868 +0.959091 1 0 0 0 100.959 Solving mesh displacement system... 4 iterations. Postprocessing: Writing graphical output: output-mesh_deformation_external_01/solution/solution-00006 - Topography min/max: -0.0105 m, 0.1 m + Topography min/max: -0.015 m, 0.1 m *** Timestep 7: t=0.0035 seconds, dt=0.0005 seconds update() compute_updated_velocities_at_points() Solution at evaluation points: rank 0: -0.3 1 0 0 0 100.3 -0.5 1 0 0 0 100.5 -rank 1: -0.72 1 0 0 0 100.72 +0.05 1 0 0 0 100.05 +0.140909 1 0 0 0 100.141 +0.231818 1 0 0 0 100.232 +0.322727 1 0 0 0 100.323 +0.413636 1 0 0 0 100.414 +0.504545 1 0 0 0 100.505 +0.595455 1 0 0 0 100.595 +0.686364 1 0 0 0 100.686 +0.777273 1 0 0 0 100.777 +0.868182 1 0 0 0 100.868 +0.959091 1 0 0 0 100.959 Solving mesh displacement system... 4 iterations. Postprocessing: Writing graphical output: output-mesh_deformation_external_01/solution/solution-00007 - Topography min/max: -0.01225 m, 0.105 m + Topography min/max: -0.0175 m, 0.105 m *** Timestep 8: t=0.004 seconds, dt=0.0005 seconds update() compute_updated_velocities_at_points() Solution at evaluation points: rank 0: -0.3 1 0 0 0 100.3 -0.5 1 0 0 0 100.5 -rank 1: -0.72 1 0 0 0 100.72 +0.05 1 0 0 0 100.05 +0.140909 1 0 0 0 100.141 +0.231818 1 0 0 0 100.232 +0.322727 1 0 0 0 100.323 +0.413636 1 0 0 0 100.414 +0.504545 1 0 0 0 100.505 +0.595455 1 0 0 0 100.595 +0.686364 1 0 0 0 100.686 +0.777273 1 0 0 0 100.777 +0.868182 1 0 0 0 100.868 +0.959091 1 0 0 0 100.959 Solving mesh displacement system... 4 iterations. Postprocessing: Writing graphical output: output-mesh_deformation_external_01/solution/solution-00008 - Topography min/max: -0.014 m, 0.12 m + Topography min/max: -0.02 m, 0.12 m Termination requested by criterion: end time diff --git a/tests/mesh_deformation_external_01/statistics b/tests/mesh_deformation_external_01/statistics index 9ee80130b94..74e9abd6681 100644 --- a/tests/mesh_deformation_external_01/statistics +++ b/tests/mesh_deformation_external_01/statistics @@ -8,12 +8,12 @@ # 8: Visualization file name # 9: Minimum topography (m) # 10: Maximum topography (m) -0 0.000000000000e+00 0.000000000000e+00 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00000 0.00000000e+00 1.00000000e-01 -1 5.000000000000e-04 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00001 -1.75000000e-03 1.00000000e-01 -2 1.000000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00002 -3.50000000e-03 1.00000000e-01 -3 1.500000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00003 -5.25000000e-03 1.00000000e-01 -4 2.000000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00004 -7.00000000e-03 1.00000000e-01 -5 2.500000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00005 -8.75000000e-03 1.00000000e-01 -6 3.000000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00006 -1.05000000e-02 1.00000000e-01 -7 3.500000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00007 -1.22500000e-02 1.05000000e-01 -8 4.000000000000e-03 5.000000000000e-04 256 2467 1089 0 output-mesh_deformation_external_01/solution/solution-00008 -1.40000000e-02 1.20000000e-01 +0 0.000000000000e+00 0.000000000000e+00 64 659 289 0 output-mesh_deformation_external_01/solution/solution-00000 0.00000000e+00 1.00000000e-01 +1 5.000000000000e-04 5.000000000000e-04 64 659 289 0 output-mesh_deformation_external_01/solution/solution-00001 -2.50000000e-03 1.00000000e-01 +2 1.000000000000e-03 5.000000000000e-04 64 659 289 0 output-mesh_deformation_external_01/solution/solution-00002 -5.00000000e-03 1.00000000e-01 +3 1.500000000000e-03 5.000000000000e-04 64 659 289 0 output-mesh_deformation_external_01/solution/solution-00003 -7.50000000e-03 1.00000000e-01 +4 2.000000000000e-03 5.000000000000e-04 64 659 289 0 output-mesh_deformation_external_01/solution/solution-00004 -1.00000000e-02 1.00000000e-01 +5 2.500000000000e-03 5.000000000000e-04 64 659 289 0 output-mesh_deformation_external_01/solution/solution-00005 -1.25000000e-02 1.00000000e-01 +6 3.000000000000e-03 5.000000000000e-04 64 659 289 0 output-mesh_deformation_external_01/solution/solution-00006 -1.50000000e-02 1.00000000e-01 +7 3.500000000000e-03 5.000000000000e-04 64 659 289 0 output-mesh_deformation_external_01/solution/solution-00007 -1.75000000e-02 1.05000000e-01 +8 4.000000000000e-03 5.000000000000e-04 64 659 289 0 output-mesh_deformation_external_01/solution/solution-00008 -2.00000000e-02 1.20000000e-01 diff --git a/tests/mesh_deformation_external_01_3d.prm b/tests/mesh_deformation_external_01_3d.prm index fe4a8002c23..2de99d4e8df 100644 --- a/tests/mesh_deformation_external_01_3d.prm +++ b/tests/mesh_deformation_external_01_3d.prm @@ -2,11 +2,11 @@ # derived class written as a plugin, 3d box geometry, initial # topography coming from the plugin (and not initial topography) # -# MPI: 2 +# MPI: 1 set Dimension = 3 -set Use years in output instead of seconds = false -set End time = 0.001 +set Use years instead of seconds = false +set End time = 0.004 set Maximum time step = 0.0005 set Nonlinear solver scheme = no Advection, no Stokes set Pressure normalization = surface diff --git a/tests/mesh_deformation_external_01_3d/screen-output b/tests/mesh_deformation_external_01_3d/screen-output index 7effc939cde..9bc8e35d3e8 100644 --- a/tests/mesh_deformation_external_01_3d/screen-output +++ b/tests/mesh_deformation_external_01_3d/screen-output @@ -15,14 +15,184 @@ Number of mesh deformation degrees of freedom: 375 *** Timestep 0: t=0 seconds, dt=0 seconds update() setting points +map_dof_to_eval_point (dof, evaluation_point_index, component): +225 0 0 +226 0 1 +227 0 2 +228 5 0 +229 5 1 +230 5 2 +231 1 0 +232 1 1 +233 1 2 +234 6 0 +235 6 1 +236 6 2 +237 10 0 +238 10 1 +239 10 2 +240 11 0 +241 11 1 +242 11 2 +243 2 0 +244 2 1 +245 2 2 +246 7 0 +247 7 1 +248 7 2 +249 12 0 +250 12 1 +251 12 2 +252 0 0 +253 0 1 +254 0 2 +255 5 0 +256 5 1 +257 5 2 +258 1 0 +259 1 1 +260 1 2 +261 6 0 +262 6 1 +263 6 2 +264 10 0 +265 10 1 +266 10 2 +267 11 0 +268 11 1 +269 11 2 +270 2 0 +271 2 1 +272 2 2 +273 7 0 +274 7 1 +275 7 2 +276 12 0 +277 12 1 +278 12 2 +279 15 0 +280 15 1 +281 15 2 +282 16 0 +283 16 1 +284 16 2 +285 20 0 +286 20 1 +287 20 2 +288 21 0 +289 21 1 +290 21 2 +291 17 0 +292 17 1 +293 17 2 +294 22 0 +295 22 1 +296 22 2 +297 15 0 +298 15 1 +299 15 2 +300 16 0 +301 16 1 +302 16 2 +303 20 0 +304 20 1 +305 20 2 +306 21 0 +307 21 1 +308 21 2 +309 17 0 +310 17 1 +311 17 2 +312 22 0 +313 22 1 +314 22 2 +315 3 0 +316 3 1 +317 3 2 +318 8 0 +319 8 1 +320 8 2 +321 13 0 +322 13 1 +323 13 2 +324 4 0 +325 4 1 +326 4 2 +327 9 0 +328 9 1 +329 9 2 +330 14 0 +331 14 1 +332 14 2 +333 3 0 +334 3 1 +335 3 2 +336 8 0 +337 8 1 +338 8 2 +339 13 0 +340 13 1 +341 13 2 +342 4 0 +343 4 1 +344 4 2 +345 9 0 +346 9 1 +347 9 2 +348 14 0 +349 14 1 +350 14 2 +351 18 0 +352 18 1 +353 18 2 +354 23 0 +355 23 1 +356 23 2 +357 19 0 +358 19 1 +359 19 2 +360 24 0 +361 24 1 +362 24 2 +363 18 0 +364 18 1 +365 18 2 +366 23 0 +367 23 1 +368 23 2 +369 19 0 +370 19 1 +371 19 2 +372 24 0 +373 24 1 +374 24 2 compute_updated_velocities_at_points() Solution at evaluation points: rank 0: +0.1 0.1 1 0 0 0 -6.70092e-18 10.1 +0.1 0.3 1 0 0 0 -6.70092e-18 10.3 +0.1 0.5 1 0 0 0 -6.70092e-18 10.5 +0.1 0.7 1 0 0 0 -6.70092e-18 10.7 +0.1 0.9 1 0 0 0 -6.70092e-18 10.9 +0.3 0.1 1 0 0 0 -6.70092e-18 30.1 +0.3 0.3 1 0 0 0 -6.70092e-18 30.3 +0.3 0.5 1 0 0 0 -6.70092e-18 30.5 0.3 0.7 1 0 0 0 -6.70092e-18 30.7 +0.3 0.9 1 0 0 0 -6.70092e-18 30.9 +0.5 0.1 1 0 0 0 -6.70092e-18 50.1 +0.5 0.3 1 0 0 0 -6.70092e-18 50.3 0.5 0.5 1 0 0 0 -6.70092e-18 50.5 -0.7 0.2 1 0 0 0 -6.70092e-18 70.2 -0.7 0.91 1 0 0 0 -6.70092e-18 70.91 -rank 1: +0.5 0.7 1 0 0 0 -6.70092e-18 50.7 +0.5 0.9 1 0 0 0 -6.70092e-18 50.9 +0.7 0.1 1 0 0 0 -6.70092e-18 70.1 +0.7 0.3 1 0 0 0 -6.70092e-18 70.3 +0.7 0.5 1 0 0 0 -6.70092e-18 70.5 +0.7 0.7 1 0 0 0 -6.70092e-18 70.7 +0.7 0.9 1 0 0 0 -6.70092e-18 70.9 +0.9 0.1 1 0 0 0 -6.70092e-18 90.1 +0.9 0.3 1 0 0 0 -6.70092e-18 90.3 +0.9 0.5 1 0 0 0 -6.70092e-18 90.5 +0.9 0.7 1 0 0 0 -6.70092e-18 90.7 0.9 0.9 1 0 0 0 -6.70092e-18 90.9 Solving mesh displacement system... 4 iterations. @@ -35,34 +205,288 @@ update() compute_updated_velocities_at_points() Solution at evaluation points: rank 0: +0.1 0.1 1 0 0 0 -6.70092e-18 10.1 +0.1 0.3 1 0 0 0 -6.70092e-18 10.3 +0.1 0.5 1 0 0 0 -6.70092e-18 10.5 +0.1 0.7 1 0 0 0 -6.70092e-18 10.7 +0.1 0.9 1 0 0 0 -6.70092e-18 10.9 +0.3 0.1 1 0 0 0 -6.70092e-18 30.1 +0.3 0.3 1 0 0 0 -6.70092e-18 30.3 +0.3 0.5 1 0 0 0 -6.70092e-18 30.5 0.3 0.7 1 0 0 0 -6.70092e-18 30.7 +0.3 0.9 1 0 0 0 -6.70092e-18 30.9 +0.5 0.1 1 0 0 0 -6.70092e-18 50.1 +0.5 0.3 1 0 0 0 -6.70092e-18 50.3 0.5 0.5 1 0 0 0 -6.70092e-18 50.5 -0.7 0.2 1 0 0 0 -6.70092e-18 70.2 -0.7 0.91 1 0 0 0 -6.70092e-18 70.91 -rank 1: +0.5 0.7 1 0 0 0 -6.70092e-18 50.7 +0.5 0.9 1 0 0 0 -6.70092e-18 50.9 +0.7 0.1 1 0 0 0 -6.70092e-18 70.1 +0.7 0.3 1 0 0 0 -6.70092e-18 70.3 +0.7 0.5 1 0 0 0 -6.70092e-18 70.5 +0.7 0.7 1 0 0 0 -6.70092e-18 70.7 +0.7 0.9 1 0 0 0 -6.70092e-18 70.9 +0.9 0.1 1 0 0 0 -6.70092e-18 90.1 +0.9 0.3 1 0 0 0 -6.70092e-18 90.3 +0.9 0.5 1 0 0 0 -6.70092e-18 90.5 +0.9 0.7 1 0 0 0 -6.70092e-18 90.7 0.9 0.9 1 0 0 0 -6.70092e-18 90.9 Solving mesh displacement system... 4 iterations. Postprocessing: Writing graphical output: output-mesh_deformation_external_01_3d/solution/solution-00001 - Topography min/max: -0.00175 m, 0.1 m + Topography min/max: 0 m, 0.1 m *** Timestep 2: t=0.001 seconds, dt=0.0005 seconds update() compute_updated_velocities_at_points() Solution at evaluation points: rank 0: +0.1 0.1 1 0 0 0 -6.70092e-18 10.1 +0.1 0.3 1 0 0 0 -6.70092e-18 10.3 +0.1 0.5 1 0 0 0 -6.70092e-18 10.5 +0.1 0.7 1 0 0 0 -6.70092e-18 10.7 +0.1 0.9 1 0 0 0 -6.70092e-18 10.9 +0.3 0.1 1 0 0 0 -6.70092e-18 30.1 +0.3 0.3 1 0 0 0 -6.70092e-18 30.3 +0.3 0.5 1 0 0 0 -6.70092e-18 30.5 0.3 0.7 1 0 0 0 -6.70092e-18 30.7 +0.3 0.9 1 0 0 0 -6.70092e-18 30.9 +0.5 0.1 1 0 0 0 -6.70092e-18 50.1 +0.5 0.3 1 0 0 0 -6.70092e-18 50.3 0.5 0.5 1 0 0 0 -6.70092e-18 50.5 -0.7 0.2 1 0 0 0 -6.70092e-18 70.2 -0.7 0.91 1 0 0 0 -6.70092e-18 70.91 -rank 1: +0.5 0.7 1 0 0 0 -6.70092e-18 50.7 +0.5 0.9 1 0 0 0 -6.70092e-18 50.9 +0.7 0.1 1 0 0 0 -6.70092e-18 70.1 +0.7 0.3 1 0 0 0 -6.70092e-18 70.3 +0.7 0.5 1 0 0 0 -6.70092e-18 70.5 +0.7 0.7 1 0 0 0 -6.70092e-18 70.7 +0.7 0.9 1 0 0 0 -6.70092e-18 70.9 +0.9 0.1 1 0 0 0 -6.70092e-18 90.1 +0.9 0.3 1 0 0 0 -6.70092e-18 90.3 +0.9 0.5 1 0 0 0 -6.70092e-18 90.5 +0.9 0.7 1 0 0 0 -6.70092e-18 90.7 0.9 0.9 1 0 0 0 -6.70092e-18 90.9 Solving mesh displacement system... 4 iterations. Postprocessing: Writing graphical output: output-mesh_deformation_external_01_3d/solution/solution-00002 - Topography min/max: -0.0035 m, 0.1 m + Topography min/max: 0 m, 0.1 m + +*** Timestep 3: t=0.0015 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.1 0.1 1 0 0 0 -6.70092e-18 10.1 +0.1 0.3 1 0 0 0 -6.70092e-18 10.3 +0.1 0.5 1 0 0 0 -6.70092e-18 10.5 +0.1 0.7 1 0 0 0 -6.70092e-18 10.7 +0.1 0.9 1 0 0 0 -6.70092e-18 10.9 +0.3 0.1 1 0 0 0 -6.70092e-18 30.1 +0.3 0.3 1 0 0 0 -6.70092e-18 30.3 +0.3 0.5 1 0 0 0 -6.70092e-18 30.5 +0.3 0.7 1 0 0 0 -6.70092e-18 30.7 +0.3 0.9 1 0 0 0 -6.70092e-18 30.9 +0.5 0.1 1 0 0 0 -6.70092e-18 50.1 +0.5 0.3 1 0 0 0 -6.70092e-18 50.3 +0.5 0.5 1 0 0 0 -6.70092e-18 50.5 +0.5 0.7 1 0 0 0 -6.70092e-18 50.7 +0.5 0.9 1 0 0 0 -6.70092e-18 50.9 +0.7 0.1 1 0 0 0 -6.70092e-18 70.1 +0.7 0.3 1 0 0 0 -6.70092e-18 70.3 +0.7 0.5 1 0 0 0 -6.70092e-18 70.5 +0.7 0.7 1 0 0 0 -6.70092e-18 70.7 +0.7 0.9 1 0 0 0 -6.70092e-18 70.9 +0.9 0.1 1 0 0 0 -6.70092e-18 90.1 +0.9 0.3 1 0 0 0 -6.70092e-18 90.3 +0.9 0.5 1 0 0 0 -6.70092e-18 90.5 +0.9 0.7 1 0 0 0 -6.70092e-18 90.7 +0.9 0.9 1 0 0 0 -6.70092e-18 90.9 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01_3d/solution/solution-00003 + Topography min/max: 0 m, 0.1 m + +*** Timestep 4: t=0.002 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.1 0.1 1 0 0 0 -6.70092e-18 10.1 +0.1 0.3 1 0 0 0 -6.70092e-18 10.3 +0.1 0.5 1 0 0 0 -6.70092e-18 10.5 +0.1 0.7 1 0 0 0 -6.70092e-18 10.7 +0.1 0.9 1 0 0 0 -6.70092e-18 10.9 +0.3 0.1 1 0 0 0 -6.70092e-18 30.1 +0.3 0.3 1 0 0 0 -6.70092e-18 30.3 +0.3 0.5 1 0 0 0 -6.70092e-18 30.5 +0.3 0.7 1 0 0 0 -6.70092e-18 30.7 +0.3 0.9 1 0 0 0 -6.70092e-18 30.9 +0.5 0.1 1 0 0 0 -6.70092e-18 50.1 +0.5 0.3 1 0 0 0 -6.70092e-18 50.3 +0.5 0.5 1 0 0 0 -6.70092e-18 50.5 +0.5 0.7 1 0 0 0 -6.70092e-18 50.7 +0.5 0.9 1 0 0 0 -6.70092e-18 50.9 +0.7 0.1 1 0 0 0 -6.70092e-18 70.1 +0.7 0.3 1 0 0 0 -6.70092e-18 70.3 +0.7 0.5 1 0 0 0 -6.70092e-18 70.5 +0.7 0.7 1 0 0 0 -6.70092e-18 70.7 +0.7 0.9 1 0 0 0 -6.70092e-18 70.9 +0.9 0.1 1 0 0 0 -6.70092e-18 90.1 +0.9 0.3 1 0 0 0 -6.70092e-18 90.3 +0.9 0.5 1 0 0 0 -6.70092e-18 90.5 +0.9 0.7 1 0 0 0 -6.70092e-18 90.7 +0.9 0.9 1 0 0 0 -6.70092e-18 90.9 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01_3d/solution/solution-00004 + Topography min/max: 0 m, 0.1 m + +*** Timestep 5: t=0.0025 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.1 0.1 1 0 0 0 -6.70092e-18 10.1 +0.1 0.3 1 0 0 0 -6.70092e-18 10.3 +0.1 0.5 1 0 0 0 -6.70092e-18 10.5 +0.1 0.7 1 0 0 0 -6.70092e-18 10.7 +0.1 0.9 1 0 0 0 -6.70092e-18 10.9 +0.3 0.1 1 0 0 0 -6.70092e-18 30.1 +0.3 0.3 1 0 0 0 -6.70092e-18 30.3 +0.3 0.5 1 0 0 0 -6.70092e-18 30.5 +0.3 0.7 1 0 0 0 -6.70092e-18 30.7 +0.3 0.9 1 0 0 0 -6.70092e-18 30.9 +0.5 0.1 1 0 0 0 -6.70092e-18 50.1 +0.5 0.3 1 0 0 0 -6.70092e-18 50.3 +0.5 0.5 1 0 0 0 -6.70092e-18 50.5 +0.5 0.7 1 0 0 0 -6.70092e-18 50.7 +0.5 0.9 1 0 0 0 -6.70092e-18 50.9 +0.7 0.1 1 0 0 0 -6.70092e-18 70.1 +0.7 0.3 1 0 0 0 -6.70092e-18 70.3 +0.7 0.5 1 0 0 0 -6.70092e-18 70.5 +0.7 0.7 1 0 0 0 -6.70092e-18 70.7 +0.7 0.9 1 0 0 0 -6.70092e-18 70.9 +0.9 0.1 1 0 0 0 -6.70092e-18 90.1 +0.9 0.3 1 0 0 0 -6.70092e-18 90.3 +0.9 0.5 1 0 0 0 -6.70092e-18 90.5 +0.9 0.7 1 0 0 0 -6.70092e-18 90.7 +0.9 0.9 1 0 0 0 -6.70092e-18 90.9 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01_3d/solution/solution-00005 + Topography min/max: 0 m, 0.1 m + +*** Timestep 6: t=0.003 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.1 0.1 1 0 0 0 -6.70092e-18 10.1 +0.1 0.3 1 0 0 0 -6.70092e-18 10.3 +0.1 0.5 1 0 0 0 -6.70092e-18 10.5 +0.1 0.7 1 0 0 0 -6.70092e-18 10.7 +0.1 0.9 1 0 0 0 -6.70092e-18 10.9 +0.3 0.1 1 0 0 0 -6.70092e-18 30.1 +0.3 0.3 1 0 0 0 -6.70092e-18 30.3 +0.3 0.5 1 0 0 0 -6.70092e-18 30.5 +0.3 0.7 1 0 0 0 -6.70092e-18 30.7 +0.3 0.9 1 0 0 0 -6.70092e-18 30.9 +0.5 0.1 1 0 0 0 -6.70092e-18 50.1 +0.5 0.3 1 0 0 0 -6.70092e-18 50.3 +0.5 0.5 1 0 0 0 -6.70092e-18 50.5 +0.5 0.7 1 0 0 0 -6.70092e-18 50.7 +0.5 0.9 1 0 0 0 -6.70092e-18 50.9 +0.7 0.1 1 0 0 0 -6.70092e-18 70.1 +0.7 0.3 1 0 0 0 -6.70092e-18 70.3 +0.7 0.5 1 0 0 0 -6.70092e-18 70.5 +0.7 0.7 1 0 0 0 -6.70092e-18 70.7 +0.7 0.9 1 0 0 0 -6.70092e-18 70.9 +0.9 0.1 1 0 0 0 -6.70092e-18 90.1 +0.9 0.3 1 0 0 0 -6.70092e-18 90.3 +0.9 0.5 1 0 0 0 -6.70092e-18 90.5 +0.9 0.7 1 0 0 0 -6.70092e-18 90.7 +0.9 0.9 1 0 0 0 -6.70092e-18 90.9 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01_3d/solution/solution-00006 + Topography min/max: 0 m, 0.1 m + +*** Timestep 7: t=0.0035 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.1 0.1 1 0 0 0 -6.70092e-18 10.1 +0.1 0.3 1 0 0 0 -6.70092e-18 10.3 +0.1 0.5 1 0 0 0 -6.70092e-18 10.5 +0.1 0.7 1 0 0 0 -6.70092e-18 10.7 +0.1 0.9 1 0 0 0 -6.70092e-18 10.9 +0.3 0.1 1 0 0 0 -6.70092e-18 30.1 +0.3 0.3 1 0 0 0 -6.70092e-18 30.3 +0.3 0.5 1 0 0 0 -6.70092e-18 30.5 +0.3 0.7 1 0 0 0 -6.70092e-18 30.7 +0.3 0.9 1 0 0 0 -6.70092e-18 30.9 +0.5 0.1 1 0 0 0 -6.70092e-18 50.1 +0.5 0.3 1 0 0 0 -6.70092e-18 50.3 +0.5 0.5 1 0 0 0 -6.70092e-18 50.5 +0.5 0.7 1 0 0 0 -6.70092e-18 50.7 +0.5 0.9 1 0 0 0 -6.70092e-18 50.9 +0.7 0.1 1 0 0 0 -6.70092e-18 70.1 +0.7 0.3 1 0 0 0 -6.70092e-18 70.3 +0.7 0.5 1 0 0 0 -6.70092e-18 70.5 +0.7 0.7 1 0 0 0 -6.70092e-18 70.7 +0.7 0.9 1 0 0 0 -6.70092e-18 70.9 +0.9 0.1 1 0 0 0 -6.70092e-18 90.1 +0.9 0.3 1 0 0 0 -6.70092e-18 90.3 +0.9 0.5 1 0 0 0 -6.70092e-18 90.5 +0.9 0.7 1 0 0 0 -6.70092e-18 90.7 +0.9 0.9 1 0 0 0 -6.70092e-18 90.9 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01_3d/solution/solution-00007 + Topography min/max: 0 m, 0.105 m + +*** Timestep 8: t=0.004 seconds, dt=0.0005 seconds +update() +compute_updated_velocities_at_points() +Solution at evaluation points: +rank 0: +0.1 0.1 1 0 0 0 -6.70092e-18 10.1 +0.1 0.3 1 0 0 0 -6.70092e-18 10.3 +0.1 0.5 1 0 0 0 -6.70092e-18 10.5 +0.1 0.7 1 0 0 0 -6.70092e-18 10.7 +0.1 0.9 1 0 0 0 -6.70092e-18 10.9 +0.3 0.1 1 0 0 0 -6.70092e-18 30.1 +0.3 0.3 1 0 0 0 -6.70092e-18 30.3 +0.3 0.5 1 0 0 0 -6.70092e-18 30.5 +0.3 0.7 1 0 0 0 -6.70092e-18 30.7 +0.3 0.9 1 0 0 0 -6.70092e-18 30.9 +0.5 0.1 1 0 0 0 -6.70092e-18 50.1 +0.5 0.3 1 0 0 0 -6.70092e-18 50.3 +0.5 0.5 1 0 0 0 -6.70092e-18 50.5 +0.5 0.7 1 0 0 0 -6.70092e-18 50.7 +0.5 0.9 1 0 0 0 -6.70092e-18 50.9 +0.7 0.1 1 0 0 0 -6.70092e-18 70.1 +0.7 0.3 1 0 0 0 -6.70092e-18 70.3 +0.7 0.5 1 0 0 0 -6.70092e-18 70.5 +0.7 0.7 1 0 0 0 -6.70092e-18 70.7 +0.7 0.9 1 0 0 0 -6.70092e-18 70.9 +0.9 0.1 1 0 0 0 -6.70092e-18 90.1 +0.9 0.3 1 0 0 0 -6.70092e-18 90.3 +0.9 0.5 1 0 0 0 -6.70092e-18 90.5 +0.9 0.7 1 0 0 0 -6.70092e-18 90.7 +0.9 0.9 1 0 0 0 -6.70092e-18 90.9 + Solving mesh displacement system... 4 iterations. + + Postprocessing: + Writing graphical output: output-mesh_deformation_external_01_3d/solution/solution-00008 + Topography min/max: 0 m, 0.12 m Termination requested by criterion: end time diff --git a/tests/mesh_deformation_external_01_3d/statistics b/tests/mesh_deformation_external_01_3d/statistics index 5ad43574272..b51c1014899 100644 --- a/tests/mesh_deformation_external_01_3d/statistics +++ b/tests/mesh_deformation_external_01_3d/statistics @@ -8,6 +8,12 @@ # 8: Visualization file name # 9: Minimum topography (m) # 10: Maximum topography (m) -0 0.000000000000e+00 0.000000000000e+00 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00000 0.00000000e+00 1.00000000e-01 -1 5.000000000000e-04 5.000000000000e-04 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00001 -1.75000000e-03 1.00000000e-01 -2 1.000000000000e-03 5.000000000000e-04 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00002 -3.50000000e-03 1.00000000e-01 +0 0.000000000000e+00 0.000000000000e+00 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00000 0.00000000e+00 1.00000000e-01 +1 5.000000000000e-04 5.000000000000e-04 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00001 0.00000000e+00 1.00000000e-01 +2 1.000000000000e-03 5.000000000000e-04 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00002 0.00000000e+00 1.00000000e-01 +3 1.500000000000e-03 5.000000000000e-04 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00003 0.00000000e+00 1.00000000e-01 +4 2.000000000000e-03 5.000000000000e-04 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00004 0.00000000e+00 1.00000000e-01 +5 2.500000000000e-03 5.000000000000e-04 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00005 0.00000000e+00 1.00000000e-01 +6 3.000000000000e-03 5.000000000000e-04 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00006 0.00000000e+00 1.00000000e-01 +7 3.500000000000e-03 5.000000000000e-04 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00007 0.00000000e+00 1.05000000e-01 +8 4.000000000000e-03 5.000000000000e-04 64 2312 729 0 output-mesh_deformation_external_01_3d/solution/solution-00008 0.00000000e+00 1.20000000e-01 From 49d6495c63622259a41efdd2a98b19166ad86733 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Fri, 9 Jan 2026 09:51:34 -0500 Subject: [PATCH 08/10] comments --- .../external_tool_interface.h | 25 ++++++++++++++++--- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/include/aspect/mesh_deformation/external_tool_interface.h b/include/aspect/mesh_deformation/external_tool_interface.h index 656e39d946a..2e66148df32 100644 --- a/include/aspect/mesh_deformation/external_tool_interface.h +++ b/include/aspect/mesh_deformation/external_tool_interface.h @@ -106,6 +106,9 @@ namespace aspect class ExternalToolInterface : public Interface, public SimulatorAccess { public: + /** + * Main routine handling the mesh deformation + */ virtual void compute_velocity_constraints_on_boundary(const DoFHandler &mesh_deformation_dof_handler, @@ -184,10 +187,15 @@ namespace aspect LinearAlgebra::Vector interpolate_external_velocities_to_surface_support_points (const std::vector> &velocities) const; - protected: - + /** + * The list of evaluation points owned by the current process. These are the points where the + * external tool will receive the ASPECT solution values and return the updated velocities. + */ std::vector> evaluation_points; + /** + * deal.II RemotePointEvaluation object used to do point evaluation in the evaluation points. + */ std::unique_ptr> remote_point_evaluator; /** @@ -201,10 +209,19 @@ namespace aspect }; /** - * A vector to store a map between Dof indices and evaluation points + * A vector to store a map between Dof indices and evaluation points. + * + * The map will contain an entry for each DoF on the surface of the ASPECT mesh and contains + * the index of the evaluation point that is closest to the DoF. As each support point has + * several components (x,y,z velocity), the map contains one entry for each component. + * + * In a parallel computation, this map only contains entries for evaluation points owned by the + * current process. Note that the DoF indices are not necessarily locally owned. + * + * This map is used in interpolate_external_velocities_to_surface_support_points() to copy + * external velocities to each surface DoF from the closest evaluation point. */ std::vector map_dof_to_eval_point; - }; } } From fdff57ebd6ac5a9865ee98ea9fe3eb5bb03b3304 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Fri, 16 Jan 2026 11:12:39 -0500 Subject: [PATCH 09/10] remove debug output --- .../external_tool_interface.cc | 302 ------------------ tests/mesh_deformation_external_01.cc | 8 + 2 files changed, 8 insertions(+), 302 deletions(-) diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc index 52d76877ff2..8ab388f6e08 100644 --- a/source/mesh_deformation/external_tool_interface.cc +++ b/source/mesh_deformation/external_tool_interface.cc @@ -26,7 +26,6 @@ /* TODO: -- PointDataOut: move to deal.II - initial topography from external tool: for now compute_initial_deformation_on_boundary(), later refactor to provide FE vector */ @@ -36,273 +35,6 @@ namespace aspect namespace MeshDeformation { - - /** - * Produce graphical output defined in points provided by the user. - */ - template - class PointDataOut : public dealii::DataOutInterface<0, spacedim> - { - public: - /** - * Default constructor. - */ - PointDataOut() = default; - - /** - * Default destructor. - */ - ~PointDataOut() = default; - - - /** - * Build the patches for a given set of points and optionally data in each point. - * - * @param [in] locations The point locations. - * @param [in] data A vector of data values for each point. - * @param [in] data_component_names An optional vector of strings that - * describe the properties of each datum. - * @param [in] data_component_interpretations An optional vector that - * controls if the properties are interpreted as scalars, vectors, - * or tensors. Has to be of the same length as @p data_component_names. - */ - void - build_patches(const std::vector> &locations, - const std::vector> &data = {}, - const std::vector &data_component_names = {}, - const std::vector< - DataComponentInterpretation::DataComponentInterpretation> - &data_component_interpretations_ = {}) - { - Assert(data_component_names.size() == data_component_interpretations_.size(), - ExcMessage( - "When calling PointDataOut::build_patches() with data component " - "names and interpretations you need to provide as many data component " - "names as interpretations. Provide the same name for components that " - "belong to a single vector or tensor.")); - - Assert(data.size() == 0 || locations.size(), - ExcMessage("You need to either provide no data or data for each point.")); - for (const auto &datum : data) - Assert(datum.size() == data_component_names.size(), - ExcMessage("The data provided in each point needs to have the same number " - "of components as names were provided.")); - - // Prepend the "id" to the data fields provided by the user: - dataset_names.clear(); - dataset_names.emplace_back("id"); - dataset_names.insert(dataset_names.end(), - data_component_names.begin(), - data_component_names.end()); - - data_component_interpretations.clear(); - data_component_interpretations.emplace_back( - DataComponentInterpretation::component_is_scalar); - data_component_interpretations.insert( - data_component_interpretations.end(), - data_component_interpretations_.begin(), - data_component_interpretations_.end()); - - const unsigned int n_property_components = data_component_names.size(); - const unsigned int n_data_components = dataset_names.size(); - - patches.resize(locations.size()); - - for (unsigned int i = 0; i < locations.size(); ++i) - { - patches[i].vertices[0] = locations[i]; - patches[i].patch_index = i; - - // Store id and properties given by the user: - patches[i].data.reinit(n_data_components, 1); - patches[i].data(0, 0) = i; // store id - for (unsigned int property_index = 0; property_index < n_property_components; ++property_index) - patches[i].data(property_index + 1, 0) = data[i][property_index]; - } - } - - protected: - /** - * Returns the patches previously built by the build_patches() function. - */ - virtual const std::vector> & - get_patches() const override - { - return patches; - } - - /** - * Virtual function through which the names of data sets are obtained from - * this class - */ - virtual std::vector - get_dataset_names() const override - { - return dataset_names; - } - - - /** - * Overload of the respective DataOutInterface::get_nonscalar_data_ranges() - * function. See there for a more extensive documentation. - * This function is a reimplementation of the function - * DataOut_DoFData::get_nonscalar_data_ranges(). - */ - virtual std::vector< - std::tuple> - get_nonscalar_data_ranges() const override - { - std::vector< - std::tuple> - ranges; - - // Make sure the data structures were set up correctly. Since they - // can only be filled by build_patches() above, they should have - // been checked already. - Assert(dataset_names.size() == data_component_interpretations.size(), - ExcInternalError()); - - // collect the ranges of particle data - const unsigned int n_output_components = - data_component_interpretations.size(); - unsigned int output_component = 0; - for (unsigned int i = 0; i < n_output_components; /* i is updated below */) - // see what kind of data we have here. note that for the purpose of the - // current function all we care about is vector data - switch (data_component_interpretations[i]) - { - case DataComponentInterpretation::component_is_scalar: - { - // Just move component forward by one - ++i; - ++output_component; - - break; - } - case DataComponentInterpretation::component_is_part_of_vector: - { - // ensure that there is a continuous number of next space_dim - // components that all deal with vectors - Assert( - i + spacedim <= n_output_components, - Exceptions::DataOutImplementation::ExcInvalidVectorDeclaration( - i, dataset_names[i])); - for (unsigned int dd = 1; dd < spacedim; ++dd) - Assert( - data_component_interpretations[i + dd] == - DataComponentInterpretation::component_is_part_of_vector, - Exceptions::DataOutImplementation:: - ExcInvalidVectorDeclaration(i, dataset_names[i])); - - // all seems right, so figure out whether there is a common - // name to these components. if not, leave the name empty and - // let the output format writer decide what to do here - std::string name = dataset_names[i]; - for (unsigned int dd = 1; dd < spacedim; ++dd) - if (name != dataset_names[i + dd]) - { - name = ""; - break; - } - - // Finally add a corresponding range. - // - // This sort of logic is also explained in some detail in - // DataOut::build_one_patch(). - ranges.emplace_back(std::forward_as_tuple( - output_component, - output_component + spacedim - 1, - name, - DataComponentInterpretation::component_is_part_of_vector)); - - // increase the 'component' counter by the appropriate amount, - // same for 'i', since we have already dealt with all these - // components - output_component += spacedim; - i += spacedim; - - break; - } - - case DataComponentInterpretation::component_is_part_of_tensor: - { - const unsigned int size = spacedim * spacedim; - // ensure that there is a continuous number of next - // spacedim*spacedim components that all deal with tensors - Assert( - i + size <= n_output_components, - Exceptions::DataOutImplementation::ExcInvalidTensorDeclaration( - i, dataset_names[i])); - for (unsigned int dd = 1; dd < size; ++dd) - Assert( - data_component_interpretations[i + dd] == - DataComponentInterpretation::component_is_part_of_tensor, - Exceptions::DataOutImplementation:: - ExcInvalidTensorDeclaration(i, dataset_names[i])); - - // all seems right, so figure out whether there is a common - // name to these components. if not, leave the name empty and - // let the output format writer decide what to do here - std::string name = dataset_names[i]; - for (unsigned int dd = 1; dd < size; ++dd) - if (name != dataset_names[i + dd]) - { - name = ""; - break; - } - - // Finally add a corresponding range. - ranges.emplace_back(std::forward_as_tuple( - output_component, - output_component + size - 1, - name, - DataComponentInterpretation::component_is_part_of_tensor)); - - // increase the 'component' counter by the appropriate amount, - // same for 'i', since we have already dealt with all these - // components - output_component += size; - i += size; - break; - } - - default: - Assert(false, ExcNotImplemented()); - } - return ranges; - } - - private: - /** - * This is a vector of patches that is created each time build_patches() is - * called. These patches are used in the output routines of the base - * classes. - */ - std::vector> patches; - - /** - * A vector of field names for all data components stored in patches. - */ - std::vector dataset_names; - - /** - * A vector that for each of the data components of the - * current data set indicates whether they are scalar fields, parts of a - * vector-field, or any of the other supported kinds of data. - */ - std::vector - data_component_interpretations; - }; - - - - template void ExternalToolInterface:: @@ -446,16 +178,8 @@ namespace aspect if (closest_evaluation_point_and_component[i].dof_index != numbers::invalid_dof_index) map_dof_to_eval_point.push_back(closest_evaluation_point_and_component[i]); } - - // print all information: - std::cout << "map_dof_to_eval_point (dof, evaluation_point_index, component): " << std::endl; - for (const auto &dof_to_eval_point : map_dof_to_eval_point) - { - std::cout << dof_to_eval_point.dof_index << " " << dof_to_eval_point.evaluation_point_index << " " << dof_to_eval_point.component << std::endl; - } } - // Finally, also ensure that upon mesh refinement, all of the // information set herein is invalidated: this->get_signals().pre_refinement_store_user_data @@ -514,32 +238,6 @@ namespace aspect "refinement step.")); AssertDimension(velocities.size(), evaluation_points.size()); - { - // Visualize the evaluation points and their velocities - static unsigned int output_no = 0; - - PointDataOut out; - //const auto &mapping = this->get_mapping(); - std::vector> real_evaluation_points(evaluation_points.size()); - std::vector> data(evaluation_points.size(), std::vector(dim, 0.0)); - for (unsigned int i=0; i data_component_names(dim, "velocity"); - const std::vector data_component_interpretations(dim, DataComponentInterpretation::component_is_part_of_vector); - - out.build_patches(real_evaluation_points, data, data_component_names, data_component_interpretations); - - out.write_vtu_with_pvtu_record(this->get_output_directory(), "surf_points", output_no, this->get_mpi_communicator(), 4, 0); - - ++output_no; - } - - // Create the output vector. const DoFHandler &mesh_dof_handler = this->get_mesh_deformation_handler().get_mesh_deformation_dof_handler(); diff --git a/tests/mesh_deformation_external_01.cc b/tests/mesh_deformation_external_01.cc index 79a8d92b4ba..08912d7daf0 100644 --- a/tests/mesh_deformation_external_01.cc +++ b/tests/mesh_deformation_external_01.cc @@ -74,6 +74,14 @@ namespace aspect Assert(false, ExcNotImplemented()); this->set_evaluation_points (points); + + // print all information: + this->get_pcout() << "map_dof_to_eval_point (dof, evaluation_point_index, component): " << std::endl; + for (const auto &dof_to_eval_point : this->map_dof_to_eval_point) + { + this->get_pcout() << "\t" << dof_to_eval_point.dof_index << " " << dof_to_eval_point.evaluation_point_index << " " << dof_to_eval_point.component << std::endl; + } + } } From 3694cfed6f7b65660fc2eb3fe058bbf20ca62d19 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Sun, 18 Jan 2026 18:55:42 -0500 Subject: [PATCH 10/10] address comments --- .../aspect/mesh_deformation/external_tool_interface.h | 4 ++-- source/mesh_deformation/external_tool_interface.cc | 10 +++++----- tests/mesh_deformation_external_01.cc | 1 + 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/include/aspect/mesh_deformation/external_tool_interface.h b/include/aspect/mesh_deformation/external_tool_interface.h index 2e66148df32..9dc98d64b40 100644 --- a/include/aspect/mesh_deformation/external_tool_interface.h +++ b/include/aspect/mesh_deformation/external_tool_interface.h @@ -201,7 +201,7 @@ namespace aspect /** * A struct to map between DoF indices and evaluation points */ - struct dof_to_eval_point_data + struct DofToEvalPointData { types::global_dof_index dof_index; unsigned int evaluation_point_index; @@ -221,7 +221,7 @@ namespace aspect * This map is used in interpolate_external_velocities_to_surface_support_points() to copy * external velocities to each surface DoF from the closest evaluation point. */ - std::vector map_dof_to_eval_point; + std::vector map_dof_to_eval_point; }; } } diff --git a/source/mesh_deformation/external_tool_interface.cc b/source/mesh_deformation/external_tool_interface.cc index 8ab388f6e08..eec97508982 100644 --- a/source/mesh_deformation/external_tool_interface.cc +++ b/source/mesh_deformation/external_tool_interface.cc @@ -115,13 +115,13 @@ namespace aspect const DoFHandler &mesh_dof_handler = this->get_mesh_deformation_handler().get_mesh_deformation_dof_handler(); std::vector squared_distances(mesh_dof_handler.locally_owned_dofs().size(), std::numeric_limits::max()); - std::vector closest_evaluation_point_and_component(mesh_dof_handler.locally_owned_dofs().size(), - dof_to_eval_point_data {numbers::invalid_dof_index, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}); + std::vector closest_evaluation_point_and_component(mesh_dof_handler.locally_owned_dofs().size(), + DofToEvalPointData {numbers::invalid_dof_index, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int}); // TODO: do we need to support the case of more than one different mesh deformation plugin to be active? const auto boundary_ids = this->get_mesh_deformation_boundary_indicators(); - IndexSet boundary_dofs = DoFTools::extract_boundary_dofs(mesh_dof_handler, ComponentMask(dim, true), boundary_ids); + const IndexSet boundary_dofs = DoFTools::extract_boundary_dofs(mesh_dof_handler, ComponentMask(dim, true), boundary_ids); const unsigned int dofs_per_cell = mesh_dof_handler.get_fe().dofs_per_cell; std::vector local_dof_indices (dofs_per_cell); @@ -134,7 +134,7 @@ namespace aspect // Note: We assume that process_and_evaluate() does not call our lambda concurrently, otherwise we would have write - // conflicts when updating vector_with_surface_velocities and one_over_distance_vec. + // conflicts when updating closest_evaluation_point_and_component and squared_distances. const auto eval_func = [&](const ArrayView &values, const typename Utilities::MPI::RemotePointEvaluation::CellData &cell_data) @@ -161,7 +161,7 @@ namespace aspect squared_distances[cell_dof_indices[j]] = distance_sq; const unsigned int component = mesh_dof_handler.get_fe().system_to_component_index(j).first; closest_evaluation_point_and_component[cell_dof_indices[j]] = - dof_to_eval_point_data {cell_dof_indices[j], local_values[i], component}; + DofToEvalPointData {cell_dof_indices[j], local_values[i], component}; } } } diff --git a/tests/mesh_deformation_external_01.cc b/tests/mesh_deformation_external_01.cc index 08912d7daf0..177136f4182 100644 --- a/tests/mesh_deformation_external_01.cc +++ b/tests/mesh_deformation_external_01.cc @@ -120,6 +120,7 @@ namespace aspect } } + // Generate some velocities: Assert(current_solution_at_points.size() == this->evaluation_points.size(), ExcInternalError()); std::vector> velocities(current_solution_at_points.size(), Tensor<1,dim>()); if (velocities.size()>6)