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..4d95d0c19d2
--- /dev/null
+++ b/include/aspect/mesh_deformation/external_tool_interface.h
@@ -0,0 +1,196 @@
+/*
+ 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;
+
+ /**
+ * 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;
+
+
+ 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;
+
+ /**
+ * 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
+ * compute_velocity_constraints_on_boundary() function of the
+ * base class.
+ */
+ LinearAlgebra::Vector
+ interpolate_external_velocities_to_surface_support_points (const std::vector> &velocities) const;
+
+ protected:
+
+ std::vector> evaluation_points;
+
+ 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..e563fe9b5be
--- /dev/null
+++ b/source/mesh_deformation/external_tool_interface.cc
@@ -0,0 +1,531 @@
+/*
+ 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
+#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
+ {
+
+
+ /**
+ * 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::
+ 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_external_velocities_to_surface_support_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
+ // ensure we also know about the locally relevant DoFs'
+ // constraints:
+ // 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 (const types::global_dof_index index : constrained_dofs)
+ {
+ if (mesh_velocity_constraints.can_store_line(index))
+ if (mesh_velocity_constraints.is_constrained(index)==false)
+ {
+#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
+ }
+ }
+
+ // TODO: make consistent?
+ }
+
+
+ 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;
+
+ // 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
+ // 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."));
+
+ 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
+ LinearAlgebra::Vector
+ ExternalToolInterface::
+ 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 "
+ "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());
+
+ {
+ // 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();
+ 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;
+
+ // 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];
+ }
+ }
+ }
+ }
+ };
+
+ 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;
+ }
+ }
+
+
+
+ namespace MeshDeformation
+ {
+#define INSTANTIATE(dim) \
+ template class ExternalToolInterface;
+
+ ASPECT_INSTANTIATE(INSTANTIATE)
+
+#undef INSTANTIATE
+ }
+}
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