diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8ce1d1fa6e6..1fc8abbf0fa 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -278,6 +278,28 @@ if(ASPECT_WITH_FASTSCAPE)
endif()
+# Fastscapelib (C++)
+set(ASPECT_WITH_FASTSCAPELIB OFF CACHE BOOL "Whether the user wants to compile ASPECT with the landscape evolution C++ code Fastscape(lib), or not.")
+if(ASPECT_WITH_FASTSCAPELIB)
+ find_package(fastscapelib CONFIG)
+ if(${fastscapelib_FOUND})
+ set(FS_WITH_HEALPIX OFF CACHE BOOL "Turn off Fastscapelib support of Healpix grid" FORCE)
+ message(STATUS "Using ASPECT_WITH_FASTSCAPELIB = '${ASPECT_WITH_FASTSCAPELIB}'")
+ message(STATUS " fastscapelib_INCLUDE_DIR: ${fastscapelib_INCLUDE_DIRS}")
+ message(STATUS " fastscapelib_VERSION: ${fastscapelib_VERSION}")
+ # fastscapelib dependencies
+ find_package(xtensor REQUIRED)
+ message(STATUS " xtensor_INCLUDE_DIR: ${xtensor_INCLUDE_DIRS}")
+ message(STATUS " xtensor_VERSION: ${xtensor_VERSION}")
+ else()
+ message(STATUS "Fastscapelib not found. Disabling ASPECT_WITH_FASTSCAPELIB. You can specify a hint to your installation directory with fastscapelib_DIR.")
+ set(ASPECT_WITH_FASTSCAPELIB OFF CACHE BOOL "" FORCE)
+ endif()
+else()
+ message(STATUS "Using ASPECT_WITH_FASTSCAPELIB = 'OFF'")
+endif()
+
+
# NetCDF (c including parallel)
set(ASPECT_WITH_NETCDF ON CACHE BOOL "Check if the user wants to compile ASPECT with the NetCDF libraries.")
@@ -408,6 +430,7 @@ if(ASPECT_WITH_WORLD_BUILDER)
endif()
message(STATUS "")
+#include "datatypes.h"
#
# Other stuff about external tools and how we interact with the system:
@@ -789,6 +812,16 @@ if (FASTSCAPE_LIB_NAME)
endif()
endif()
+# Fastscapelib
+if (${fastscapelib_FOUND})
+ message(STATUS "Linking ASPECT against Fastscapelib (header-only)")
+ foreach(_T ${TARGET_EXECUTABLES})
+ target_include_directories(${_T} PRIVATE ${fastscapelib_INCLUDE_DIRS})
+ target_include_directories(${_T} PRIVATE ${xtensor_INCLUDE_DIRS})
+ target_include_directories(${_T} PRIVATE ${xtl_INCLUDE_DIRS})
+ endforeach()
+endif()
+
# NetCDF
if(${NETCDF_FOUND})
message(STATUS "Linking ASPECT against NetCDF")
diff --git a/WITH_FASTSCAPELIB.md b/WITH_FASTSCAPELIB.md
new file mode 100644
index 00000000000..5d15216f0bf
--- /dev/null
+++ b/WITH_FASTSCAPELIB.md
@@ -0,0 +1,117 @@
+# Using Aspect with Fastscapelib (C++)
+
+## Install Fastscapelib using conda
+
+You can use conda to install fastscapelib and all of its dependencies (we assume
+that you are familiar with managing conda environments):
+
+```
+conda install fastscapelib -c conda-forge
+```
+
+## Install Fastscapelib from source
+
+You can follow the steps below to install Fastscapelib and its dependencies from
+source.
+
+### Set the installation prefix path for Fastscapelib and its dependencies
+
+Choose a common location where you want to download the source files of
+Fastscapelib and its dependencies:
+
+```
+$ export FS_SOURCE_PATH=/path/to/fastscapelib/source
+```
+
+Choose a common location where you want to install Fastscapelib and its
+dependencies:
+
+```
+export CMAKE_INSTALL_PREFIX=/path/to/fastscapelib/install/dir
+```
+
+### Install xtl and xtensor
+
+Download the `xtl` git repository and checkout the last stable version:
+
+```
+$ cd $FS_SOURCE_PATH
+$ git clone https://github.com/xtensor-stack/xtl
+$ cd xtl
+$ git checkout 0.7.7
+```
+
+Build and install `xtl`:
+
+```
+$ cmake -S. -Bbuild
+$ cmake --install build
+```
+
+Download the `xtensor` git repository and checkout the last stable version:
+
+```
+$ cd $FS_SOURCE_PATH
+$ git clone https://github.com/xtensor-stack/xtensor
+$ cd xtensor
+$ git checkout 0.25.0
+```
+
+Build and install `xtensor`:
+
+```
+$ cmake -S. -Bbuild
+$ cmake --install build
+```
+
+### Install Fastscapelib
+
+Download the `fastscapelib` git repository and checkout the last stable version:
+
+```
+$ cd $FS_SOURCE_PATH
+$ git clone https://github.com/fastscape-lem/fastscapelib
+$ cd fastscapelib
+$ git checkout v0.2.2
+```
+
+IMPORTANT: edit the `cmake/fastscapelibConfig.cmake` file with the editor of
+your choice and add the following line just below `@PACKAGE_INIT@` (this will be
+required every time you checkout another version or branch and re-install
+Fastscapelib, but this won't be required for Fastscapelib version >0.2.2):
+
+```
+include(CMakeFindDependencyMacro)
+```
+
+Build and install `fastscapelib`:
+
+```
+$ cmake -S. -Bbuild
+$ cmake --install build
+```
+
+## Compile Aspect with Fastscapelib enabled
+
+If you have installed Fastscapelib using conda, don't forget to activate the
+environment where you have installed Fastscapelib. Once activated, configuring
+`aspect` with CMake should then find Fastscapelib automatically.
+
+If you have installed Fastscapelib from source like described above, you can
+help CMake find the Fastscapelib installation path using:
+
+```
+export CMAKE_PREFIX_PATH=/path/to/fastscapelib/install/dir
+```
+
+where the given path corresponds to the path that you set during installation
+via the `CMAKE_INSTALL_PREFIX` environment variable (see instructions above).
+
+Configuring `aspect` should then find Fastscapelib. Check the output of the
+following command run from aspect's source root directory (note: you might need
+to set additional arguments or run extra commands in order to find the other
+dependencies of `aspect`, we skipped it here):
+
+```
+cmake -S. -Bbuild -DASPECT_WITH_FASTSCAPELIB=ON
+```
diff --git a/contrib/fastscape/fastscape_rewrite_VTK_with_time.py b/contrib/fastscape/fastscape_rewrite_VTK_with_time.py
index ffde75c5bb6..f39d656de85 100644
--- a/contrib/fastscape/fastscape_rewrite_VTK_with_time.py
+++ b/contrib/fastscape/fastscape_rewrite_VTK_with_time.py
@@ -33,7 +33,7 @@ def get_times_pvd(filename):
times[i] = float(temp_str)
- return times
+ return times
#%% Get file paths (absolute, not relative paths!)
@@ -54,73 +54,73 @@ def absoluteFilePaths(directory):
# trace generated using paraview version 5.8.0
#
-# To ensure correct image size when batch processing, please search
+# To ensure correct image size when batch processing, please search
# for and uncomment the line `# renderView*.ViewSize = [*,*]`
for File in FileNames:
#### disable automatic camera reset on 'Show'
paraview.simple._DisableFirstRenderCameraReset()
-
+
# create a new 'Legacy VTK Reader'
topography0000 = LegacyVTKReader(FileNames=File)
-
+
# get animation scene
animationScene1 = GetAnimationScene()
-
+
# get the time-keeper
timeKeeper1 = GetTimeKeeper()
-
+
# update animation scene based on data timesteps
animationScene1.UpdateAnimationUsingDataTimeSteps()
-
+
# get active view
renderView1 = GetActiveViewOrCreate('RenderView')
# uncomment following to set a specific view size
# renderView1.ViewSize = [882, 554]
-
+
# get layout
layout1 = GetLayout()
-
+
# show data in view
topography0000Display = Show(topography0000, renderView1, 'StructuredGridRepresentation')
-
+
# trace defaults for the display properties.
topography0000Display.Representation = 'Surface'
-
+
# reset view to fit data
renderView1.ResetCamera()
-
+
# get the material library
materialLibrary1 = GetMaterialLibrary()
-
+
# show color bar/color legend
topography0000Display.SetScalarBarVisibility(renderView1, True)
-
+
# update the view to ensure updated data information
renderView1.Update()
-
+
# get color transfer function/color map for 'H'
hLUT = GetColorTransferFunction('H')
-
+
# get opacity transfer function/opacity map for 'H'
hPWF = GetOpacityTransferFunction('H')
-
+
# save data
SaveData(File[:-4] + '.vts', proxy=topography0000, ChooseArraysToWrite=1,
PointDataArrays=['HHHHH','basement','catchment','drainage_area','erosion_rate','topography','total_erosion'])
-
+
#### saving camera placements for all active views
-
+
# current camera placement for renderView1
renderView1.CameraPosition = [225625.0, 25625.0, 878497.0779980461]
renderView1.CameraFocalPoint = [225625.0, 25625.0, 1135.7277145385742]
renderView1.CameraParallelScale = 227077.82689023562
-
+
#### uncomment the following to render all views
# RenderAllViews()
# alternatively, if you want to write images, you can use SaveScreenshot(...).
-
+
Delete(topography0000)
del topography0000
@@ -140,17 +140,9 @@ def absoluteFilePaths(directory):
for n, File in enumerate(FileNames):
dataset_attributes = {"timestep": str(times[n]), "group": "", "group": "", "file": os.path.basename(File)[:-4]+'.vts'}
-
+
ET.SubElement(doc, "DataSet", attrib=dataset_attributes).text="\n"
tree = ET.ElementTree(root)
-tree.write('./VTK/topography.pvd', xml_declaration=True, encoding='utf-8', method="xml")
-
-
-
-
-
-
-
-
+tree.write('./VTK/topography.pvd', xml_declaration=True, encoding='utf-8', method="xml")
diff --git a/include/aspect/compat.h b/include/aspect/compat.h
index 00ae9ef9bb5..4a157671b32 100644
--- a/include/aspect/compat.h
+++ b/include/aspect/compat.h
@@ -213,4 +213,220 @@ namespace aspect
#endif
+
+// deal.II version 9.7 introduces a new class VectorFunctionFromTensorFunctionObject
+// that we would like to use also for earlier versions
+#if !DEAL_II_VERSION_GTE(9,7,0)
+
+namespace aspect
+{
+ using namespace dealii;
+
+ /**
+ * This class is built as a means of translating the Tensor<1,dim,
+ * RangeNumberType> values produced by function objects that
+ * for a given point return a tensor,
+ * and returning them as a multiple component version of the same thing as a
+ * Vector for use in, for example, the VectorTools::interpolate or the many
+ * other functions taking Function objects. It allows the user to place the
+ * desired components into an n_components long vector starting at
+ * the selected_component location in that vector and have all other
+ * components be 0.
+ *
+ * For example: Say you created a function object that returns the gravity
+ * (here, a radially inward pointing vector of magnitude 9.81):
+ * @code
+ * const auto gravity
+ * = [](const Point &p) -> Tensor<1,dim> { return -9.81 * (p /
+ * p.norm()); }
+ * @endcode
+ * and you want to interpolate this onto your mesh using the
+ * VectorTools::interpolate function, with a finite element for the
+ * DoFHandler object has 3 copies of a finite element with dim
+ * components, for a total of 3*dim components. To interpolate onto that
+ * DoFHandler, you need an object of type Function that has 3*dim vector
+ * components. Creating such an object from the existing gravity
+ * object is done using this piece of code:
+ * @code
+ * VectorFunctionFromTensorFunctionObject
+ * gravity_vector_function(gravity, 0, 3*dim);
+ * @endcode
+ * where the dim components of the `gravity` function are placed
+ * into the first dim components of the function object.
+ *
+ * @ingroup functions
+ */
+ template
+ class VectorFunctionFromTensorFunctionObject
+ : public Function
+ {
+ public:
+ /**
+ * Given a function object that takes a Point and returns a
+ * Tensor<1,dim, RangeNumberType> value, convert this into an object
+ * that matches the Function@ interface.
+ *
+ * By default, create a Vector object of the same size as
+ * tensor_function returns, i.e., with dim components.
+ *
+ * @param tensor_function_object The TensorFunction that will form `dim` components of
+ * the resulting Vector Function object.
+ * @param n_components The total number of vector components of the
+ * resulting TensorFunction object. This clearly has to be at least `dim`.
+ * @param selected_component The first component that should be filled by
+ * the first argument. This should be such that the entire tensor_function
+ * fits inside the n_component length return vector.
+ */
+ explicit VectorFunctionFromTensorFunctionObject(
+ const std::function(const Point &)>
+ &tensor_function_object,
+ const unsigned int selected_component = 0,
+ const unsigned int n_components = dim);
+
+ /**
+ * This destructor is defined as virtual so as to coincide with all other
+ * aspects of class.
+ */
+ virtual ~VectorFunctionFromTensorFunctionObject() override = default;
+
+ /**
+ * Return a single component of a vector-valued function at a given point.
+ */
+ virtual RangeNumberType
+ value(const Point &p, const unsigned int component = 0) const override;
+
+ /**
+ * Return all components of a vector-valued function at a given point.
+ *
+ * values shall have the right size beforehand, i.e. #n_components.
+ */
+ virtual void
+ vector_value(const Point &p,
+ Vector &values) const override;
+
+ /**
+ * Return all components of a vector-valued function at a list of points.
+ *
+ * value_list shall be the same size as points and each
+ * element of the vector will be passed to vector_value() to evaluate the
+ * function
+ */
+ virtual void
+ vector_value_list(
+ const std::vector> &points,
+ std::vector> &value_list) const override;
+
+ private:
+ /**
+ * The TensorFunction object which we call when this class's vector_value()
+ * or vector_value_list() functions are called.
+ */
+ const std::function(const Point &)>
+ tensor_function_object;
+
+ /**
+ * The first vector component whose value is to be filled by the given
+ * TensorFunction. The values will be placed in components
+ * selected_component to selected_component+dim-1 for a
+ * TensorFunction<1,dim, RangeNumberType> object.
+ */
+ const unsigned int selected_component;
+ };
+
+
+ template
+ VectorFunctionFromTensorFunctionObject::
+ VectorFunctionFromTensorFunctionObject(
+ const std::function(const Point &)>
+ &tensor_function_object,
+ const unsigned int selected_component,
+ const unsigned int n_components)
+ : Function(n_components)
+ , tensor_function_object(tensor_function_object)
+ , selected_component(selected_component)
+ {
+ // Verify that the Tensor<1,dim,RangeNumberType> will fit in the given length
+ // selected_components and not hang over the end of the vector.
+ AssertIndexRange(selected_component + dim - 1, this->n_components);
+ }
+
+
+
+ template
+ inline RangeNumberType
+ VectorFunctionFromTensorFunctionObject::value(
+ const Point &p,
+ const unsigned int component) const
+ {
+ AssertIndexRange(component, this->n_components);
+
+ // if the requested component is out of the range selected, then we can
+ // return early
+ if ((component < selected_component) ||
+ (component >= selected_component + dim))
+ return 0;
+
+ // otherwise retrieve the values from the tensor_function to be
+ // placed at the selected_component to
+ // selected_component + dim - 1 elements of the Vector
+ // values and pick the correct one
+ const Tensor<1, dim, RangeNumberType> tensor_value =
+ tensor_function_object(p);
+
+ return tensor_value[component - selected_component];
+ }
+
+
+ template
+ inline void
+ VectorFunctionFromTensorFunctionObject::vector_value(
+ const Point &p,
+ Vector &values) const
+ {
+ Assert(values.size() == this->n_components,
+ ExcDimensionMismatch(values.size(), this->n_components));
+
+ // Retrieve the values from the tensor_function to be placed at
+ // the selected_component to
+ // selected_component + dim - 1 elements of the Vector
+ // values.
+ const Tensor<1, dim, RangeNumberType> tensor_value =
+ tensor_function_object(p);
+
+ // First we make all elements of values = 0
+ values = 0;
+
+ // Second we adjust the desired components to take on the values in
+ // tensor_value.
+ for (unsigned int i = 0; i < dim; ++i)
+ values(i + selected_component) = tensor_value[i];
+ }
+
+
+ /**
+ * Member function vector_value_list is the interface for giving a
+ * list of points (vector>) of which to evaluate
+ * using the vector_value member function. Again, this function is
+ * written so as to not replicate the function definition but passes each
+ * point on to vector_value to be evaluated.
+ */
+ template
+ void
+ VectorFunctionFromTensorFunctionObject::vector_value_list(
+ const std::vector> &points,
+ std::vector> &value_list) const
+ {
+ Assert(value_list.size() == points.size(),
+ ExcDimensionMismatch(value_list.size(), points.size()));
+
+ const unsigned int n_points = points.size();
+
+ for (unsigned int p = 0; p < n_points; ++p)
+ VectorFunctionFromTensorFunctionObject::vector_value(
+ points[p], value_list[p]);
+ }
+}
+
+#endif
+
#endif
diff --git a/include/aspect/config.h.in b/include/aspect/config.h.in
index 19f699d4c41..16231850be8 100644
--- a/include/aspect/config.h.in
+++ b/include/aspect/config.h.in
@@ -32,6 +32,7 @@
#cmakedefine ASPECT_WITH_WORLD_BUILDER
#cmakedefine ASPECT_USE_FP_EXCEPTIONS
#cmakedefine ASPECT_WITH_FASTSCAPE
+#cmakedefine ASPECT_WITH_FASTSCAPELIB
#cmakedefine ASPECT_HAVE_FASTSCAPE_NAMED_VTK
#define ASPECT_MAX_NUM_PARTICLE_SYSTEMS @ASPECT_MAX_NUM_PARTICLE_SYSTEMS@
diff --git a/include/aspect/mesh_deformation/fastscapecc.h b/include/aspect/mesh_deformation/fastscapecc.h
new file mode 100644
index 00000000000..9d96bf9981c
--- /dev/null
+++ b/include/aspect/mesh_deformation/fastscapecc.h
@@ -0,0 +1,350 @@
+/*
+ Copyright (C) 2023 - 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
+ .
+*/
+
+#ifndef _aspect_mesh_deformation_fastscapecc_h
+#define _aspect_mesh_deformation_fastscapecc_h
+
+#include
+
+#ifdef ASPECT_WITH_FASTSCAPELIB
+
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+
+
+#include
+
+#include
+#include
+
+namespace aspect
+{
+ using namespace dealii;
+
+ namespace MeshDeformation
+ {
+ /**
+ * A comparator for dealii::Point to use in std::map
+ * ------------------------------------------------------
+ * std::map requires a way to compare keys (i.e., define a "less than" relation)
+ * dealii::Point does not provide operator< by default,
+ * so we must define one ourselves.
+ *
+ * This comparator compares two Point objects lexicographically:
+ * First compares x, then y (and z if dim == 3).
+ *
+ * This is required to allow using std::map, T>
+ * to uniquely index surface vertices on unstructured spherical meshes.
+ */
+ template
+ struct PointComparator
+ {
+ bool operator()(const dealii::Point &a, const dealii::Point &b) const
+ {
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ if (a[d] < b[d]) return true;
+ if (a[d] > b[d]) return false;
+ }
+ return false;
+ }
+ };
+
+ template
+ void output_fastscape_visualization(
+ const unsigned int timestep_index,
+ const dealii::DoFHandler &dof_handler,
+ const xt::xarray &elevation,
+ const xt::xarray &uplift_rate,
+ const xt::xarray &erosion,
+ const xt::xarray &drainage_area,
+ const xt::xarray &sediment_flux,
+ const std::string &output_directory,
+ const MPI_Comm mpi_communicator,
+ const dealii::ConditionalOStream &pcout);
+
+
+ /**
+ * A plugin that utilizes the landscape evolution code Fastscapelib (C++)
+ * to deform the ASPECT boundary through erosion.
+ */
+ template
+ class FastScapecc : public Interface, public SimulatorAccess
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ FastScapecc();
+
+ /**
+ * Initialize variables for FastScape.
+ */
+ void initialize () override;
+
+ /**
+ * A function that creates constraints for the velocity of certain mesh
+ * vertices (namely, the surface vertices) for a specific boundary.
+ * The calling class will respect
+ * these constraints when computing the new vertex positions.
+ */
+ void
+ compute_velocity_constraints_on_boundary(const DoFHandler &mesh_deformation_dof_handler,
+ AffineConstraints &mesh_velocity_constraints,
+ const std::set &boundary_id) const override;
+
+
+ double interpolate_surface_velocity(const Point &p,
+ const std::vector &V) const;
+ /**
+ * Returns whether or not the plugin requires surface stabilization.
+ */
+ bool needs_surface_stabilization () const override;
+
+ /**
+ * Declare parameters for the FastScape plugin.
+ */
+ static
+ void declare_parameters (ParameterHandler &prm);
+
+ /**
+ * Parse parameters for the FastScape plugin.
+ */
+ void parse_parameters (ParameterHandler &prm) override;
+
+
+ private:
+
+ // friend class boost::serialization::access;
+
+ /**
+ * For a given location, typically the location of a vertex of ASPECT's
+ * volume mesh, return the DoF index that corresponds to this vertex
+ * location within the surface_dof_handler. This index can then also
+ * be used to index into vectors defined on surface_dof_handler.
+ */
+ unsigned int vertex_index(const Point &p) const;
+
+ /*
+ * A map, used only for the spherical geometry, that maps the location
+ * of a vertex to an index between zero and the number of vertices
+ * in the surface mesh.
+ *
+ * TODO: This map is only used for the spherical shell geometry. We
+ * should prevent this object from being used altogether in all
+ * other cases, and this could be facilitated by wrapping it in
+ * std_cxx17::optional<>. Only when we encounter the spherical
+ * geometry do we actually set anything in the optional.
+ */
+ std::map, unsigned int, PointComparator> spherical_vertex_index_map;
+
+ /**
+ * The surface mesh is a distributed Triangulation to support large-scale
+ * parallel simulations (e.g., global models).
+ */
+ using SurfaceMeshType = parallel::distributed::Triangulation;
+
+ SurfaceMeshType surface_mesh;
+ DoFHandler surface_mesh_dof_handler;
+ mutable AffineConstraints surface_constraints;
+
+ FE_Q surface_fe;
+
+ /**
+ * Initialize the surface mesh based on the given geometry model.
+ */
+ void init_surface_mesh(const GeometryModel::Interface &);
+
+ /**
+ * Pointers to Fastscapelib objects
+ */
+ using GridAdapterType = typename fastscapelib::dealii_grid;
+ using FlowGraphType = typename fastscapelib::flow_graph;
+ std::unique_ptr grid;
+ std::unique_ptr flow_graph;
+ std::unique_ptr> spl_eroder;
+ // mutable std::shared_ptr> surface_cache;
+
+
+ void project_surface_solution(const std::set &boundary_ids,
+ dealii::LinearAlgebra::distributed::Vector &surface_vertical_velocity,
+ dealii::LinearAlgebra::distributed::Vector &surface_elevation) const;
+
+
+ /**
+ * Project the Stokes velocity solution onto the
+ * free surface. Called by make_constraints()
+ */
+ // void project_velocity_onto_boundary (const DoFHandler &free_surface_dof_handler,
+ // const IndexSet &mesh_locally_owned,
+ // const IndexSet &mesh_locally_relevant,
+ // LinearAlgebra::Vector &output) const;
+
+
+ /**
+ * Fill velocity data table to be interpolated back onto the ASPECT mesh.
+ */
+ // Table fill_data_table(std::vector &values,
+ // TableIndices &size_idx,
+ // const int &array_size) const;
+
+ /**
+ * Run-time parameters
+ * @{
+ */
+ /**
+ * Suggestion for the number of FastScape steps to run for every ASPECT timestep,
+ * where the FastScape time step is determined by ASPECT's time step length divided by
+ * this parameter.
+ */
+ unsigned int fastscape_steps_per_aspect_step;
+
+ /**
+ * Maximum time step allowed for FastScape. If the suggested time step exceeds this
+ * limit it is repeatedly divided by 2 until the final time step is smaller than this parameter.
+ */
+ double maximum_fastscape_timestep;
+
+ /**
+ * Total number of Fastscapelib grid nodes.
+ */
+ unsigned int n_grid_nodes;
+
+ /**
+ * How many levels FastScape should be refined above the maximum ASPECT surface resolution.
+ */
+ unsigned int additional_refinement_levels;
+
+ /**
+ * Maximum expected refinement level at ASPECT's surface.
+ * This and resolution_difference are required to properly transfer node data from
+ * ASPECT to FastScape.
+ */
+ unsigned int surface_refinement_level;
+
+ /**
+ * Difference in refinement levels expected at the ASPECT surface,
+ * where this would be set to 2 if 3 refinement levels are set at the surface.
+ * This and surface_resolution are required to properly transfer node data from
+ * ASPECT to FastScape.
+ *
+ * TODO: Should this be kept this way, or make it so the input is the expected levels
+ * of refinement at the surface, and we can subtract one within the code? Also,
+ * it would be good to find a way to check these are correct, because they are a
+ * common source of errors.
+ */
+ int surface_refinement_difference;
+
+ /**
+ * Magnitude (m) of the initial noise applied to FastScape.
+ * Applied as either a + or - value to the topography
+ * such that the total difference can be up to 2*noise_h.
+ */
+ double noise_h;
+
+ // FastScape erosional parameters //
+ /**
+ * Drainage area exponent for the stream power law (m variable in FastScape surface equation).
+ */
+ double m;
+
+ /**
+ * Slope exponent for the steam power law (n variable in FastScape surface equation).
+ */
+ double n;
+
+ /**
+ * Bedrock river incision rate for the stream power law
+ * (meters^(1-2m)/yr, kf variable in FastScape surface equation).
+ */
+ double kff;
+
+ /**
+ * Sediment river incision rate for the stream power law (meters^(1-2m)/yr).
+ * When set to -1 this is identical to the bedrock value.
+ * (kf variable in FastScape surface equation applied to sediment).
+ */
+ double kfsed;
+
+ /**
+ * Bedrock transport coefficient for hillslope diffusion (m^2/yr, kd in FastScape surface equation).
+ */
+ double kdd;
+
+ /**
+ * Bedrock transport coefficient for hillslope diffusion (m^2/yr). When set to -1 this is
+ * identical to the bedrock value.
+ * (kd in FastScape surface equation applied to sediment).
+ */
+ double kdsed;
+ /**
+ * @}
+ */
+
+
+ // Grid extent in each direction [min, max]
+ // std::array, 2> grid_extent_surface;
+
+ std::array,dim> box_grid_extent;
+
+
+ // Grid resolution (spacing between points)
+ double dx;
+ double dy;
+
+ // Number of grid cells in x and y directions
+ unsigned int nx;
+ unsigned int ny;
+
+ mutable bool elevation_initialized = false;
+ mutable xt::xarray elevation;
+
+ /**
+ * The number of cells in each coordinate direction.
+ */
+ std::array repetitions;
+ std::array, dim - 1> grid_extent_surface;
+
+ // A way to ouput internal fastscape iteration for visualisation
+ bool output_internal_fastscape_steps;
+
+
+ };
+ }
+}
+
+#endif // #ifdef ASPECT_WITH_FASTSCAPELIB
+
+#endif
diff --git a/include/aspect/mesh_deformation/fastscapecc_adapter.h b/include/aspect/mesh_deformation/fastscapecc_adapter.h
new file mode 100644
index 00000000000..cf18b867db7
--- /dev/null
+++ b/include/aspect/mesh_deformation/fastscapecc_adapter.h
@@ -0,0 +1,394 @@
+/*
+ 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_fastscapecc_adapter_h
+#define _aspect_mesh_deformation_fastscapecc_adapter_h
+
+#include
+
+#ifdef ASPECT_WITH_FASTSCAPELIB
+
+#include
+#include
+
+// Include the Xtensor configuration file so that we can query the
+// Xtensor version number below to address the fact that between
+// 0.25.0 and 0.26.0, all Xtensor include files were moved around,
+// see https://github.com/xtensor-stack/xtensor/pull/2829 and the
+// complaints therein.
+#if __has_include() // 0.26.0 and later
+# include
+#else // 0.25.0 and earlier
+# include
+#endif
+
+#if XTENSOR_VERSION_MAJOR ==0 && XTENSOR_VERSION_MINOR <= 25
+# include
+# include
+# include
+# include
+#else
+# include
+# include
+# include
+# include
+#endif
+
+#include
+
+// Up until FastScape++ 0.2.2, there was a file xtensor_utils.hpp that declared
+// alias types that are used elsewhere. After that, the file went away, see
+// https://github.com/fastscape-lem/fastscapelib/pull/155
+// but we still need the contents and so have to declare them ourselves.
+#include
+#if FASTSCAPELIB_VERSION_MAJOR*10000 + FASTSCAPELIB_VERSION_MINOR*100 + FASTSCAPELIB_VERSION_PATCH <= 202
+# include
+#else
+
+namespace fastscapelib
+{
+ /**
+ * Used to get the actual xtensor container type from a given selector.
+ *
+ * @tparam S The xtensor selector type.
+ * @tparam T The container value type.
+ * @tparam N The number of dimensions (only for static dimension containers)
+ */
+ template
+ struct xt_container
+ {
+ };
+
+ template
+ struct xt_container
+ {
+ using tensor_type = xt::xtensor;
+ using array_type = xt::xarray;
+ };
+
+ /**
+ * Alias for the selected (static dimension) xtensor container type.
+ *
+ * @tparam S The xtensor selector type.
+ * @tparam T The container value type.
+ * @tparam N The fixed number of dimensions.
+ */
+ template
+ using xt_tensor_t = typename xt_container::tensor_type;
+
+ /**
+ * Alias for the selected (dynamic dimension) xtensor container type.
+ *
+ * @tparam S The xtensor selector type.
+ * @tparam T The container value type.
+ */
+ template
+ using xt_array_t = typename xt_container::array_type;
+
+} // namespace fastscapelib
+#endif
+
+
+namespace fastscapelib
+{
+ template
+ class dealii_grid;
+
+ /**
+ * Dealii grid specialized types
+ */
+ template
+ struct grid_inner_types>
+ {
+ static constexpr bool is_structured = false;
+ static constexpr bool is_uniform = false;
+
+ using grid_data_type = double;
+
+ using xt_selector = S;
+ static constexpr std::size_t xt_ndims = 1;
+
+ static constexpr uint8_t n_neighbors_max = 8;
+ using neighbors_cache_type = neighbors_no_cache;
+ };
+
+ /**
+ * @brief Grid adapter class for using a DEAL.II surface or boundary mesh with
+ * Fastscapelib
+ *
+ * This adapter should work with various kinds of DEAL.II meshes (e.g., box,
+ * sphere, etc.). The only requirement is that the mesh must be uniform (i.e.,
+ * all cells have an equal area) at a given refinement level.
+ *
+ * Note: a ``dealii_grid`` node cooresponds to a DEAL.II mesh cell.
+ *
+ * @tparam T The DEAL.II mesh type.
+ * @tparam S The container selector for data array members.
+ */
+ template
+ class dealii_grid : public grid>
+ {
+ public:
+ using self_type = dealii_grid;
+ using base_type = grid;
+ using inner_types = grid_inner_types;
+
+ using grid_data_type = typename base_type::grid_data_type;
+
+ using xt_selector = typename base_type::xt_selector;
+ using container_type = xt_tensor_t;
+ using neighbors_type = typename base_type::neighbors_type;
+ using neighbors_indices_type = typename base_type::neighbors_indices_type;
+ using neighbors_distances_type = typename base_type::neighbors_distances_type;
+
+ using nodes_status_type = typename base_type::nodes_status_type;
+ using nodes_status_array_type = xt_tensor_t;
+
+ using size_type = typename base_type::size_type;
+ using shape_type = typename base_type::shape_type;
+
+ dealii_grid(T &triangulation, double cell_area);
+
+ void set_nodes_status(const nodes_status_array_type &nodes_status);
+
+ protected:
+ T &m_triangulation;
+
+ shape_type m_shape;
+ size_type m_size;
+ double m_node_area;
+
+ nodes_status_type m_nodes_status;
+
+ using neighbors_distances_impl_type = typename base_type::neighbors_distances_impl_type;
+ using neighbors_indices_impl_type = typename base_type::neighbors_indices_impl_type;
+
+ std::vector m_neighbors_count;
+ std::vector m_neighbors_indices;
+ std::vector m_neighbors_distances;
+
+ void compute_connectivity();
+
+ inline container_type nodes_areas_impl() const;
+ inline grid_data_type nodes_areas_impl(const size_type &idx) const noexcept;
+
+ inline size_type neighbors_count_impl(const size_type &idx) const;
+
+ void neighbors_indices_impl(neighbors_indices_impl_type &neighbors,
+ const size_type &idx) const;
+
+ inline const neighbors_distances_impl_type &neighbors_distances_impl(
+ const size_type &idx) const;
+
+ static constexpr std::size_t dimension_impl() noexcept;
+
+ friend class grid;
+ };
+
+
+ /**
+ * Creates a new adapter from an existing DEAL.II surface mesh.
+ *
+ * @param triangulation The DEAL.II triangulation object.
+ * @param cell_area the uniform area of the mesh cells.
+ */
+ template
+ dealii_grid::dealii_grid(T &triangulation, double cell_area)
+ : base_type(0)
+ , m_triangulation(triangulation)
+ , m_node_area(cell_area)
+ {
+
+ m_size = triangulation.n_global_active_cells();
+ m_shape = {{ static_cast(m_size) }};
+
+ // pre-compute explicit grid connectivity
+ compute_connectivity();
+
+ // Set all grid nodes as "active"
+ // - Fastscapelib doesn't support yet labelling grid nodes as ghost nodes
+ // - we assume that all the erosion processes that we apply here depend on flow routing
+ // (i.e., they rely on a `fastscapelib::flow_graph` object built on top of this grid,
+ // such as `fastscapelib::spl_eroder`)
+ // - It is better (more explicit) to set base level nodes via the flow graph
+ // object (i.e., `fastscapelib::flow_graph::set_base_levels`)
+ // - Masks (e.g., oceans) can be defined via the flow graph too
+ // (i.e., `fastscapelib::flow_graph::set_mask`)
+ // - TODO: periodic boundary conditions (if any) should be handled in `compute_connectivity`
+ // (there is no need to label nodes as periodic boundaries since the grid connectivity
+ // is computed explicitly here)
+ // m_nodes_status = node_status::core;
+ m_nodes_status = xt::xarray::from_shape({m_size});
+ std::fill(m_nodes_status.begin(), m_nodes_status.end(), node_status::core);
+
+
+ }
+
+ template
+ void dealii_grid::set_nodes_status(const nodes_status_array_type &nodes_status)
+ {
+ if (!xt::same_shape(nodes_status.shape(), m_shape))
+ {
+ throw std::invalid_argument(
+ "invalid shape for nodes_status array (expects shape [N] where N is the total number of nodes)");
+ }
+ m_nodes_status = nodes_status;
+ }
+
+// template
+// void dealii_grid::compute_connectivity()
+// {
+// m_neighbors_count.resize(m_size);
+// m_neighbors_indices.resize(m_size);
+// m_neighbors_distances.resize(m_size);
+// unsigned int counter = 0;
+// std::unordered_map global_to_local_index;
+// unsigned int local_index = 0;
+// // for (const auto &cell : m_triangulation.get_dof_handler().active_cell_iterators())
+// for (const auto &cell : m_triangulation.active_cell_iterators())
+// {
+// m_neighbors_count[counter] = cell->n_faces();
+// auto center = cell.center();
+// for (auto face_i = 0; face_i < m_neighbors_count[counter]; ++face_i)
+// {
+// unsigned int global_neighbor_index = 0;
+// for (const auto &neighbor_cell : cell.neighbor(face_i))
+// {
+// global_neighbor_index = neighbor_cell.global_active_cell_index();
+// m_neighbors_distances[counter][face_i] = center.distance(neighbor_cell.center());
+// // TODO: I assume we only have one neighbor per face, since we start with a uniform grid, we will only have one, but cwe could generatlize this
+// }
+// // TODO: add assert to check if neighbor_index is correclty set
+// global_to_local_index.emplace(global_neighbor_index, counter);
+// }
+
+// counter++;
+// }
+// counter = 0;
+// // for (const auto &cell : this->get_dof_handler().active_cell_iterators())
+// for (const auto &cell : m_triangulation.active_cell_iterators())
+// {
+// for (auto k = 0; k <= m_neighbors_count[counter]; k++)
+// {
+// auto global_nb_index = m_neighbors_indices[counter][k];
+// m_neighbors_indices[counter][k] = global_to_local_index.at(global_nb_index);
+// }
+// counter++;
+// }
+// }
+
+
+ template
+ void fastscapelib::dealii_grid::compute_connectivity()
+ {
+ // Resize internal storage for neighbor counts, indices, and distances
+ m_neighbors_count.resize(m_size);
+ m_neighbors_indices.resize(m_size);
+ m_neighbors_distances.resize(m_size);
+
+ // Step 1: Map global active cell indices to local indices
+ // This is necessary to translate global indices into local array indices
+ std::unordered_map global_to_local_index;
+ unsigned int local_index = 0;
+ for (const auto &cell : m_triangulation.active_cell_iterators())
+ {
+ global_to_local_index[cell->global_active_cell_index()] = local_index++;
+ }
+
+ // Step 2: Loop through each active cell and determine valid neighbors
+ local_index = 0;
+ for (const auto &cell : m_triangulation.active_cell_iterators())
+ {
+ const auto ¢er = cell->center(); // Get cell center for distance computation
+ unsigned int valid_neighbors = 0;
+
+ for (unsigned int face_i = 0; face_i < cell->n_faces(); ++face_i)
+ {
+ // Only consider internal faces with valid neighbor
+ // if (!cell->at_boundary(face_i) && cell->neighbor_index_is_valid(face_i))
+ if (!cell->at_boundary(face_i) && cell->neighbor(face_i).state() == dealii::IteratorState::valid)
+ {
+ const auto &neighbor = cell->neighbor(face_i);
+ const unsigned int global_nb_index = neighbor->global_active_cell_index();
+
+ auto it = global_to_local_index.find(global_nb_index);
+ AssertThrow(it != global_to_local_index.end(),
+ aspect::ExcMessage("Neighbor index not found in global_to_local_index map."));
+
+
+ // Fill neighbor index and distance arrays
+ m_neighbors_indices[local_index][valid_neighbors] = it->second;
+ m_neighbors_distances[local_index][valid_neighbors] = center.distance(neighbor->center());
+
+ ++valid_neighbors;
+ }
+ }
+
+ // Save how many neighbors this node has
+ m_neighbors_count[local_index] = valid_neighbors;
+
+ ++local_index;
+ }
+ }
+
+
+ template
+ inline auto dealii_grid::nodes_areas_impl() const -> container_type
+ {
+ return xt::broadcast(m_node_area, m_shape);
+ }
+
+ template
+ inline auto dealii_grid::nodes_areas_impl(const size_type & /*idx*/) const noexcept
+ -> grid_data_type
+ {
+ return m_node_area;
+ }
+
+ template
+ inline auto dealii_grid::neighbors_count_impl(const size_type &idx) const -> size_type
+ {
+ return m_neighbors_count[idx];
+ }
+
+ template
+ void dealii_grid::neighbors_indices_impl(neighbors_indices_impl_type &neighbors,
+ const size_type &idx) const
+ {
+ const auto &size = m_neighbors_count[idx];
+
+ for (size_type i = 0; i < size; i++)
+ {
+ neighbors[i] = m_neighbors_indices[idx][i];
+ }
+ }
+
+ template
+ auto dealii_grid::neighbors_distances_impl(const size_type &idx) const
+ -> const neighbors_distances_impl_type &
+ {
+ return m_neighbors_distances[idx];
+ }
+}
+
+#endif // #ifdef ASPECT_WITH_FASTSCAPELIB
+
+
+#endif // _aspect_mesh_deformation_fastscapecc_adapter_h
diff --git a/source/mesh_deformation/fastscapecc.cc b/source/mesh_deformation/fastscapecc.cc
new file mode 100644
index 00000000000..b0cd79e4e80
--- /dev/null
+++ b/source/mesh_deformation/fastscapecc.cc
@@ -0,0 +1,1146 @@
+/*
+ Copyright (C) 2011 - 2022 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
+
+#ifdef ASPECT_WITH_FASTSCAPELIB
+
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include // For Utilities::MPI::this_mpi_process
+#include // For Utilities::create_directory
+#include // C++17 or later
+
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+
+
+
+
+namespace aspect
+{
+ namespace MeshDeformation
+ {
+ template
+ FastScapecc::FastScapecc()
+ :
+ surface_mesh(MPI_COMM_WORLD),
+ surface_fe(1) // Use a Q1 element for the surface mesh
+ // , surface_mesh_dof_handler(surface_mesh) // Link DoFHandler to surface_mesh
+ {}
+
+ struct ValuesAtSurfaceVertex
+ {
+ double topography;
+ double surface_uplift_rate;
+
+ template
+ void serialize(Archive &ar, const unsigned int /*version*/)
+ {
+ ar & topography;
+ ar & surface_uplift_rate;
+ }
+ };
+
+ template
+ void FastScapecc::initialize()
+ {
+ const auto &geom_model = this->get_geometry_model();
+ init_surface_mesh(geom_model);
+
+ //Fastscapelib operates on grid nodes, each of which conceptually represents a cell center
+ //and has associated area and neighbor connections.
+ // n_grid_nodes = surface_mesh.n_active_cells();
+ n_grid_nodes = surface_mesh.n_used_vertices();
+ }
+
+ template
+ void FastScapecc::init_surface_mesh(const GeometryModel::Interface &geom_model)
+ {
+ // First create the surface triangulation objects. This is necessarily different
+ // between geometries.
+ if (const auto *box = dynamic_cast *>(&geom_model))
+ {
+ this->get_pcout() << "Box geometry detected. Initializing FastScape for Box geometry..." << std::endl;
+
+ // Get origin and extent of the ASPECT domain
+ const auto origin = box->get_origin();
+ const auto extent = box->get_extents();
+
+ for (unsigned int i = 0; i < dim; ++i)
+ {
+ box_grid_extent[i].first = origin[i];
+ box_grid_extent[i].second = origin[i] + extent[i];
+ }
+
+ // Build surface mesh corners from grid extent
+ Point p1, p2;
+ for (unsigned int i = 0; i < dim - 1; ++i)
+ {
+ p1[i] = box_grid_extent[i].first;
+ p2[i] = box_grid_extent[i].second;
+ }
+
+ // Extract surface repetitions from full domain repetitions
+ std::vector surface_repetitions(repetitions.begin(), repetitions.begin() + (dim - 1));
+
+ // Apply surface refinement factor
+ const unsigned int total_refinement = surface_refinement_level + additional_refinement_levels;
+ const unsigned int refinement_factor = static_cast(std::pow(2.0, total_refinement));
+
+ for (unsigned int i = 0; i < dim - 1; ++i)
+ surface_repetitions[i] *= refinement_factor;
+
+ // Store grid dimensions for indexing and spacing
+ nx = surface_repetitions[0]+1;
+ ny = surface_repetitions[1]+1;
+ dx = (box_grid_extent[0].second - box_grid_extent[0].first) / static_cast(nx-1);
+ dy = (box_grid_extent[1].second - box_grid_extent[1].first) / static_cast(ny-1);
+
+ // Create refined surface mesh
+ GridGenerator::subdivided_hyper_rectangle(surface_mesh,
+ surface_repetitions,
+ p1,
+ p2,
+ /* colorize */ true);
+ Assert (surface_mesh.n_used_vertices() == nx*ny, ExcInternalError());
+ }
+ else if (const auto *spherical_shell = dynamic_cast *>(&geom_model))
+ {
+ this->get_pcout() << "Spherical Shell geometry detected. Initializing FastScape for Spherical Shell geometry..." << std::endl;
+
+ const Point center;
+ GridGenerator::hyper_sphere(surface_mesh, center, spherical_shell->outer_radius());
+ surface_mesh.refine_global(3); // Adjust as needed
+
+ // Create a unique index for each surface vertex
+ unsigned int counter = 0;
+ for (const auto &cell : surface_mesh.active_cell_iterators())
+ for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; ++v)
+ {
+ const Point &vertex = cell->vertex(v);
+ if (spherical_vertex_index_map.find(vertex) == spherical_vertex_index_map.end())
+ {
+ spherical_vertex_index_map[vertex] = counter++;
+ }
+ }
+
+ n_grid_nodes = spherical_vertex_index_map.size();
+ }
+ else
+ AssertThrow(false, ExcMessage("FastScapecc plugin only supports Box or Spherical Shell geometries."));
+
+ // Having create the mesh, now set up the DoFHandlers and constraints objects
+ surface_mesh_dof_handler.reinit(surface_mesh);
+ surface_mesh_dof_handler.distribute_dofs(surface_fe);
+
+ surface_constraints.clear();
+ DoFTools::make_hanging_node_constraints(surface_mesh_dof_handler, surface_constraints);
+ surface_constraints.close();
+ }
+
+
+
+ template
+ void FastScapecc::project_surface_solution(const std::set & /*boundary_ids*/,
+ dealii::LinearAlgebra::distributed::Vector &surface_vertical_velocity,
+ dealii::LinearAlgebra::distributed::Vector &surface_elevation ) const
+ {
+ TimerOutput::Scope timer_section(this->get_computing_timer(), "Project surface solution");
+
+ const auto *spherical_model = dynamic_cast *>(&this->get_geometry_model());
+ const auto *box_model = dynamic_cast *>(&this->get_geometry_model());
+
+ AssertThrow(spherical_model || box_model,
+ ExcMessage("FastScapecc only supports Box or SphericalShell geometries."));
+
+ // Initialize the surface solution vector on the FastScape surface mesh
+ const IndexSet locally_relevant_dofs_surface
+ = DoFTools::extract_locally_relevant_dofs(surface_mesh_dof_handler);
+
+ surface_vertical_velocity.reinit(surface_mesh_dof_handler.locally_owned_dofs(),
+ locally_relevant_dofs_surface,
+ this->get_mpi_communicator());
+
+ surface_elevation.reinit(surface_mesh_dof_handler.locally_owned_dofs(),
+ locally_relevant_dofs_surface,
+ this->get_mpi_communicator());
+ const types::boundary_id top_boundary =
+ this->get_geometry_model().translate_symbolic_boundary_name_to_id("top");
+ // std::cout<<"here it works 0a"<get_dof_handler().active_cell_iterators())
+ if (volume_cell->is_locally_owned())
+ {
+ for (unsigned int face = 0; face < volume_cell->n_faces(); ++face)
+ {
+ if (!volume_cell->face(face)->at_boundary() ||
+ volume_cell->face(face)->boundary_id() != top_boundary)
+ continue;
+
+ for (unsigned int v = 0; v < volume_cell->face(face)->n_vertices(); ++v)
+ {
+ // std::cout<<"here it works 0a"< pos = volume_cell->face(face)->vertex(v);
+ // std::cout<<"here it works 0b"< velocity;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ // Query velocity Dofs using introspection
+ const unsigned int component_index = this->introspection().component_indices.velocities[d];
+ velocity[d] = this->get_solution()[volume_cell->face(face)->vertex_dof_index(v, component_index)];
+ }
+ // std::cout<<"here it works 0c"< gravity = this->get_gravity_model().gravity_vector(pos);
+ Tensor<1, dim> gravity_dir = gravity;
+
+ const double gravity_norm = gravity.norm();
+ if (gravity_norm > 0.0)
+ gravity_dir /= gravity_norm;
+ // std::cout<<"here it works 0d"<get_geometry_model().height_above_reference_surface(pos);
+ const unsigned int index = this->vertex_index(pos);
+ // std::cout<<"here it works 0e"<get_pcout() << "Projected surface solution onto FastScape mesh." << std::endl;
+
+ }
+
+ template
+ unsigned int FastScapecc::vertex_index(const Point &p) const
+ {
+ if (const auto *box_model = dynamic_cast *>(&this->get_geometry_model()))
+ {
+ if constexpr (dim == 3)
+ {
+ const unsigned int i =std::min(static_cast((p[0] - box_grid_extent[0].first) / dx + 0.5), nx - 1);
+ const unsigned int j =std::min(static_cast((p[1] - box_grid_extent[1].first) / dy + 0.5), ny - 1);
+ const unsigned int index = j * nx + i;
+
+ // std::cout << "[DEBUG vertex_index] p: (" << x << ", " << y << "), "
+ // << "box_grid_extent x0: " << x0 << ", y0: " << y0 << ", "
+ // << "dx: " << dx << ", dy: " << dy << ", "
+ // << "i: " << i << ", j: " << j << ", nx: " << nx << ", "
+ // << "index: " << index << std::endl;
+
+ return index;
+ }
+ else if constexpr (dim == 2)
+ {
+ const double x = p[0];
+ const double x0 = box_grid_extent[0].first;
+ const unsigned int i = static_cast((x - x0) / dx + 0.5);
+
+ std::cout << "[DEBUG vertex_index 2D] p: " << x
+ << ", box_grid_extent x0: " << x0
+ << ", dx: " << dx << ", i: " << i << std::endl;
+
+ return i;
+ }
+ }
+ else if (const auto *spherical_model = dynamic_cast *>(&this->get_geometry_model()))
+ {
+ const auto it = spherical_vertex_index_map.find(p);
+ AssertThrow(it != spherical_vertex_index_map.end(),
+ ExcMessage("Point not found in spherical_vertex_index_map."));
+ return it->second;
+ }
+ else
+ {
+ AssertThrow(false, ExcMessage("Unsupported geometry in vertex_index()."));
+ return numbers::invalid_unsigned_int;
+ }
+ }
+
+
+ // Make visualisation its own function for later
+ // template
+ // void output_fastscape_visualization(
+ // const unsigned int timestep_index,
+ // const dealii::DoFHandler &dof_handler,
+ // const xt::xarray &elevation,
+ // const xt::xarray &uplift_rate,
+ // const xt::xarray &erosion,
+ // const xt::xarray &drainage_area,
+ // const xt::xarray &sediment_flux,
+ // const std::string &output_directory,
+ // const MPI_Comm mpi_communicator,
+ // const dealii::ConditionalOStream &pcout)
+ // {
+ // if (!std::filesystem::exists(output_directory))
+ // std::filesystem::create_directories(output_directory);
+
+ // // Flatten in case they're 2D arrays
+ // const auto elev_flat = xt::flatten(elevation);
+ // const auto uplift_flat = xt::flatten(uplift_rate);
+ // const auto erosion_flat = xt::flatten(erosion);
+ // const auto drainage_flat = xt::flatten(drainage_area);
+ // const auto sediment_flat = xt::flatten(sediment_flux);
+
+ // AssertThrow(elev_flat.size() == uplift_flat.size() &&
+ // elev_flat.size() == erosion_flat.size() &&
+ // elev_flat.size() == drainage_flat.size() &&
+ // elev_flat.size() == sediment_flat.size(),
+ // dealii::ExcMessage("All FastScape input arrays must have same number of elements."));
+
+ // const unsigned int n_grid_nodes = elev_flat.size();
+
+ // // === Prepare deal.II vectors ===
+ // dealii::Vector elevation_output(n_grid_nodes);
+ // dealii::Vector uplift_rate_output(n_grid_nodes);
+ // dealii::Vector erosion_output(n_grid_nodes);
+ // dealii::Vector drainage_area_output(n_grid_nodes);
+ // dealii::Vector sediment_flux_output(n_grid_nodes);
+
+ // for (unsigned int j = 0; j < n_grid_nodes; ++j)
+ // {
+ // elevation_output[j] = elev_flat(j);
+ // uplift_rate_output[j] = uplift_flat(j);
+ // erosion_output[j] = erosion_flat(j);
+ // drainage_area_output[j] = drainage_flat(j);
+ // sediment_flux_output[j] = sediment_flat(j);
+ // }
+
+ // // === Build VTU file ===
+ // dealii::DataOut data_out;
+ // data_out.attach_dof_handler(dof_handler);
+ // data_out.add_data_vector(elevation_output, "Elevation");
+ // data_out.add_data_vector(uplift_rate_output, "UpliftRate");
+ // data_out.add_data_vector(erosion_output, "Erosion");
+ // data_out.add_data_vector(drainage_area_output, "DrainageArea");
+ // data_out.add_data_vector(sediment_flux_output, "SedimentFlux");
+ // data_out.build_patches();
+
+ // const unsigned int mpi_rank = dealii::Utilities::MPI::this_mpi_process(mpi_communicator);
+ // const unsigned int mpi_size = dealii::Utilities::MPI::n_mpi_processes(mpi_communicator);
+
+ // const std::string basename = "fastscape_timestep_" + dealii::Utilities::int_to_string(timestep_index, 4);
+ // const std::string filename_base = output_directory + "/" + basename;
+ // const std::string vtu_filename = filename_base + "." + dealii::Utilities::int_to_string(mpi_rank, 4) + ".vtu";
+
+ // {
+ // std::ofstream output(vtu_filename);
+ // data_out.write_vtu(output);
+ // }
+
+ // if (mpi_rank == 0)
+ // {
+ // // Write .pvtu file
+ // std::vector vtu_filenames;
+ // for (unsigned int r = 0; r < mpi_size; ++r)
+ // vtu_filenames.push_back(basename + "." + dealii::Utilities::int_to_string(r, 4) + ".vtu");
+
+ // const std::string pvtu_filename = filename_base + ".pvtu";
+ // {
+ // std::ofstream pvtu(pvtu_filename);
+ // data_out.write_pvtu_record(pvtu, vtu_filenames);
+ // }
+
+ // // Append or create .pvd file
+ // const std::string pvd_filename = output_directory + "/fastscape_timeseries.pvd";
+ // std::vector existing_entries;
+
+ // if (std::filesystem::exists(pvd_filename))
+ // {
+ // std::ifstream in(pvd_filename);
+ // std::string line;
+ // while (std::getline(in, line))
+ // if (line.find("";
+ // existing_entries.push_back(new_entry.str());
+
+ // std::ofstream pvd_file(pvd_filename);
+ // pvd_file << "\n"
+ // << "\n"
+ // << "