Merge pull request #5559 from bangerth/initial-topography-cleanups

Better document and check initial topography models.

include/aspect/geometry_model/initial_topography_model/interface.h

@@ -67,10 +67,17 @@ namespace aspect
         virtual void initialize ();
 
         /**
-         * Return the value of the elevation at the given point.
+         * Return the value of the elevation at the given surface point.
+         *
+         * Note that different geometry models use different conventions for
+         * how they describe surface points. In general, the models use
+         * their own "natural" coordinate system. For example, box-type
+         * geometry models will generally provide points as x-y coordinates
+         * on the surface, whereas spherical-type geometry models will generally
+         * provide surface points in spherical coordinates.
          */
         virtual
-        double value (const Point<dim-1> &p) const = 0;
+        double value (const Point<dim-1> &surface_point) const = 0;
 
         /**
          * Return the maximum value of the elevation.
@@ -97,7 +104,6 @@ namespace aspect
         virtual
         void
         parse_parameters (ParameterHandler &prm);
-
     };
 
 

source/geometry_model/initial_topography_model/ascii_data.cc

@@ -58,25 +58,38 @@ namespace aspect
     double
     AsciiData<dim>::value (const Point<dim-1> &surface_point) const
     {
-      // Because the get_data_component function of AsciiDataBoundary
-      // expects a dim-dimensional cartesian point, we have to
-      // add a coordinate here, and, for spherical geometries,
-      // change to cartesian coordinates.
+      // In a first step, create a global 'dim'-dimensional point that we can pass to the
+      // function expression as input -- because the function is a dim-dimensional
+      // function.
+      //
+      // Different geometry models pass the surface point in in different ways.
+      // In the following, we will first normalize the input to a dim-dimensional
+      // point with a dummy vertical/radial coordinate that, one would hope,
+      // the AsciiDataBoundary class will then simply ignore.
       Point<dim> global_point;
       if (Plugins::plugin_type_matches<const GeometryModel::Box<dim>> (this->get_geometry_model()) ||
           Plugins::plugin_type_matches<const GeometryModel::TwoMergedBoxes<dim>> (this->get_geometry_model()))
         {
-          // No need to set the vertical coordinate correctly,
-          // because it will be thrown away in get_data_component anyway
           for (unsigned int d=0; d<dim-1; ++d)
             global_point[d] = surface_point[d];
+
+          // Now for the vertical component:
+          global_point[dim-1] = 0;
         }
       else if (Plugins::plugin_type_matches<const GeometryModel::Sphere<dim>> (this->get_geometry_model()) ||
                Plugins::plugin_type_matches<const GeometryModel::SphericalShell<dim>> (this->get_geometry_model()) ||
                Plugins::plugin_type_matches<const GeometryModel::Chunk<dim>> (this->get_geometry_model()))
         {
-          // No need to set the radial coordinate correctly,
-          // because it will be thrown away in get_data_component anyway
+          // AsciiDataBoundary always expects to get the input
+          // parameters for its functions in Cartesian
+          // coordinates. The get_data_component function then changes
+          // the coordinate system, or more precisely it asks the
+          // geometry model to convert the point into its natural
+          // coordinate system, before doing the table lookup.
+          //
+          // This is of course all not very efficient (TODO: Think
+          // about a better scheme), but the first step then needs to
+          // be to convert what we have into Cartesian coordinates...
           std::array<double, dim> point;
           point[0] = 6371000.0;
           for (unsigned int d=0; d<dim-1; ++d)
@@ -87,20 +100,24 @@ namespace aspect
       else
         AssertThrow(false, ExcNotImplemented());
 
-      const double topo = Utilities::AsciiDataBoundary<dim>::get_data_component(surface_boundary_id,
-                                                                                global_point,
-                                                                                0);
+      const double topo = this->Utilities::AsciiDataBoundary<dim>::get_data_component(surface_boundary_id,
+                                                                                      global_point,
+                                                                                      0);
 
       return topo;
     }
 
+
+
     template <int dim>
     Tensor<1,dim-1>
     AsciiData<dim>::vector_gradient(const Point<dim> &point) const
     {
       return Utilities::AsciiDataBoundary<dim>::vector_gradient(surface_boundary_id, point,0);
     }
 
+
+
     template <int dim>
     double
     AsciiData<dim>::max_topography () const
@@ -109,6 +126,7 @@ namespace aspect
     }
 
 
+
     template <int dim>
     void
     AsciiData<dim>::declare_parameters (ParameterHandler &prm)

source/geometry_model/initial_topography_model/function.cc

@@ -41,32 +41,55 @@ namespace aspect
       coordinate_system(Utilities::Coordinates::CoordinateSystem::cartesian)
     {}
 
+
+
     template <int dim>
     double
     Function<dim>::
     value (const Point<dim-1> &surface_point) const
     {
+      // In a first step, create a global 'dim'-dimensional point that we can pass to the
+      // function expression as input -- because the function is a dim-dimensional
+      // function.
+      //
+      // Different geometry models pass the surface point in in different ways.
+      // In the following, we will first normalize the input to a dim-dimensional
+      // point with a dummy vertical/radial coordinate that, one would hope,
+      // the function expression will then simply ignore.
       Point<dim> global_point;
+
       if (Plugins::plugin_type_matches<GeometryModel::Box<dim>>(this->get_geometry_model()) ||
           Plugins::plugin_type_matches<const GeometryModel::TwoMergedBoxes<dim>> (this->get_geometry_model()))
         {
-          // No need to set the vertical coordinate correctly,
-          // because it will be thrown away in get_data_component anyway
           for (unsigned int d=0; d<dim-1; ++d)
             global_point[d] = surface_point[d];
+
+          // Now for the vertical component:
+          global_point[dim-1] = 0;
+
+          // The point as it is would have to be translated into a different
+          // coordinate system if that was requested in the input file.
+          // This is not currently implemented.
+          Assert (coordinate_system == Utilities::Coordinates::CoordinateSystem::cartesian,
+                  ExcNotImplemented());
         }
       else if (Plugins::plugin_type_matches<GeometryModel::Sphere<dim>>(this->get_geometry_model()) ||
                Plugins::plugin_type_matches<GeometryModel::SphericalShell<dim>>(this->get_geometry_model()) ||
                Plugins::plugin_type_matches<GeometryModel::Chunk<dim>>(this->get_geometry_model()) )
         {
-          // No need to set the radial coordinate correctly,
-          // because it will be thrown away in get_data_component anyway
           std::array<double, dim> point;
           point[0] = 6371000.0;
           for (unsigned int d=0; d<dim-1; ++d)
             point[d+1] = surface_point[d];
 
           global_point = Utilities::Coordinates::spherical_to_cartesian_coordinates<dim>(point);
+
+          // The point as it is would have to be translated into a different
+          // coordinate system (or, perhaps, better just left in the spherical
+          // coordinates we received) if that was requested in the input file.
+          // This is not currently implemented.
+          Assert (coordinate_system == Utilities::Coordinates::CoordinateSystem::cartesian,
+                  ExcNotImplemented());
         }
       else
         AssertThrow(false, ExcNotImplemented());

source/geometry_model/initial_topography_model/prm_polygon.cc

@@ -32,7 +32,7 @@ namespace aspect
     PrmPolygon<dim>::
     value (const Point<dim-1> &p) const
     {
-      Point<2> p1 = (dim == 2 ?Point<2>(p[0],0) : Point<2>(p[0],p[1]));
+      const Point<2> p1 = (dim == 2 ? Point<2>(p[0],0) : Point<2>(p[0],p[1]));
 
       /**
        * We go through the loop in the reverse order, because we