precompute natural coodinates.

CHANGELOG.md

@@ -16,6 +16,7 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm
 ### Changed
 - The World Builder Visualizer will now use zlib compression for vtu files by default. If zlib is not available binary output will be used. \[Menno Fraters; 2021-06-26; [#282](github.com/GeodynamicWorldBuilder/WorldBuilder/pull/282)\]
 - The return argument type of the distance_point_from_curved_planes function has been converted from a map to a struct, requiring a change in the plugin interfaces for 'fault_models' and 'subducting_plate_models', but significantly increasing the speed of the function and all functions that access its returned values. \[Rene Gassmoeller; 2021-06-27; [#289](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/289)\]
+- The plugin systems (temperature, composition and grains) and the distance_point_from_curved_planes function now all pass a precomputed NaturalCoordinate, besides just the cartesian position. It turns out that this can make a significant performance differce. \[Issue found and solution suggested by Wolfgang Bangerth, implemented and tested by Menno Fraters; 2021-07-03; [#300](github.com/GeodynamicWorldBuilder/WorldBuilder/pull/300) and [#219](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/219)\]
 
 ## [0.4.0]
 ### Added

include/world_builder/features/continental_plate.h

@@ -85,6 +85,7 @@ namespace WorldBuilder
          * gravity and current temperature.
          */
         double temperature(const Point<3> &position,
+                           const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                            const double depth,
                            const double gravity,
                            double temperature) const override final;
@@ -96,6 +97,7 @@ namespace WorldBuilder
          * of that composition at this location and depth.
          */
         double composition(const Point<3> &position,
+                           const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                            const double depth,
                            const unsigned int composition_number,
                            double composition) const override final;
@@ -109,6 +111,7 @@ namespace WorldBuilder
 
         WorldBuilder::grains
         grains(const Point<3> &position,
+               const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                const double depth,
                const unsigned int composition_number,
                WorldBuilder::grains grains) const override final;

include/world_builder/features/fault.h

@@ -87,6 +87,7 @@ namespace WorldBuilder
          * gravity and current temperature.
          */
         double temperature(const Point<3> &position,
+                           const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                            const double depth,
                            const double gravity,
                            double temperature) const override final;
@@ -98,6 +99,7 @@ namespace WorldBuilder
          * of that composition at this location and depth.
          */
         double composition(const Point<3> &position,
+                           const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                            const double depth,
                            const unsigned int composition_number,
                            double composition_value) const override final;
@@ -111,6 +113,7 @@ namespace WorldBuilder
 
         WorldBuilder::grains
         grains(const Point<3> &position,
+               const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                const double depth,
                const unsigned int composition_number,
                WorldBuilder::grains grains) const override final;

include/world_builder/features/interface.h

@@ -87,6 +87,7 @@ namespace WorldBuilder
          */
         virtual
         double temperature(const Point<3> &position,
+                           const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                            const double depth,
                            const double gravity,
                            double temperature) const = 0;
@@ -96,6 +97,7 @@ namespace WorldBuilder
          */
         virtual
         double composition(const Point<3> &position,
+                           const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                            const double depth,
                            const unsigned int composition_number,
                            double value) const = 0;
@@ -106,6 +108,7 @@ namespace WorldBuilder
          */
         virtual
         WorldBuilder::grains grains(const Point<3> &position,
+                                    const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                                     const double depth,
                                     const unsigned int composition_number,
                                     WorldBuilder::grains value) const = 0;

include/world_builder/features/mantle_layer.h

@@ -86,6 +86,7 @@ namespace WorldBuilder
          * gravity and current temperature.
          */
         double temperature(const Point<3> &position,
+                           const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                            const double depth,
                            const double gravity,
                            double temperature) const override final;
@@ -98,6 +99,7 @@ namespace WorldBuilder
          * of that composition at this location and depth.
          */
         double composition(const Point<3> &position,
+                           const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                            const double depth,
                            const unsigned int composition_number,
                            double composition) const override final;
@@ -112,6 +114,7 @@ namespace WorldBuilder
 
         WorldBuilder::grains
         grains(const Point<3> &position,
+               const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                const double depth,
                const unsigned int composition_number,
                WorldBuilder::grains grains) const override final;

include/world_builder/features/oceanic_plate.h

@@ -86,6 +86,7 @@ namespace WorldBuilder
          * gravity and current temperature.
          */
         double temperature(const Point<3> &position,
+                           const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                            const double depth,
                            const double gravity,
                            double temperature) const override final;
@@ -98,6 +99,7 @@ namespace WorldBuilder
          * of that composition at this location and depth.
          */
         double composition(const Point<3> &position,
+                           const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                            const double depth,
                            const unsigned int composition_number,
                            double composition) const override final;
@@ -112,6 +114,7 @@ namespace WorldBuilder
 
         WorldBuilder::grains
         grains(const Point<3> &position,
+               const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                const double depth,
                const unsigned int composition_number,
                WorldBuilder::grains grains) const override final;

include/world_builder/features/subducting_plate.h

@@ -89,6 +89,7 @@ namespace WorldBuilder
          * gravity and current temperature.
          */
         double temperature(const Point<3> &position,
+                           const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                            const double depth,
                            const double gravity,
                            double temperature) const override final;
@@ -101,6 +102,7 @@ namespace WorldBuilder
          * of that composition at this location and depth.
          */
         double composition(const Point<3> &position,
+                           const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                            const double depth,
                            const unsigned int composition_number,
                            double composition_value) const override final;
@@ -114,6 +116,7 @@ namespace WorldBuilder
 
         WorldBuilder::grains
         grains(const Point<3> &position,
+               const  WorldBuilder::Utilities::NaturalCoordinate &natural_coordinate,
                const double depth,
                const unsigned int composition_number,
                WorldBuilder::grains grains) const override final;

include/world_builder/utilities.h

@@ -89,7 +89,7 @@ namespace WorldBuilder
          * Returns the coordinates in the given coordinate system, which may
          * not be Cartesian.
          */
-        const std::array<double,3> &get_coordinates();
+        const std::array<double,3> &get_coordinates() const;
 
         /**
          * The coordinate that represents the 'surface' directions in the
@@ -299,8 +299,10 @@ namespace WorldBuilder
     /**
      * Computes the distance of a point to a curved plane.
      * TODO: add more info on how this works/is implemented.
-     * \param point This is the cartesian point of which we want to know the
+     * \param check_point This is the cartesian point of which we want to know the
      * distance to the curved planes
+     * \param check_point_natural the check_point in the natural coordinates of the
+     * current coordinate system.
      * \param reference_point This is a 2d point in natural coordinates at the
      * surface which the curved planes dip towards. Natural coordinates are in
      * cartesian (x,y,z) in meters and in spherical radius in meters and longitude
@@ -342,7 +344,8 @@ namespace WorldBuilder
      * of the point from the plane and the distance of the point along the plane,
      * and the average angle of the closest segment/section.
      */
-    PointDistanceFromCurvedPlanes distance_point_from_curved_planes(const Point<3> &point,
+    PointDistanceFromCurvedPlanes distance_point_from_curved_planes(const Point<3> &check_point,
+                                                                    const NaturalCoordinate &check_point_natural,
                                                                     const Point<2> &reference_point,
                                                                     const std::vector<Point<2> > &point_list,
                                                                     const std::vector<std::vector<double> > &plane_segment_lengths,

source/features/continental_plate.cc

@@ -135,12 +135,11 @@ namespace WorldBuilder
 
     double
     ContinentalPlate::temperature(const Point<3> &position,
+                                  const  NaturalCoordinate &natural_coordinate,
                                   const double depth,
                                   const double gravity_norm,
                                   double temperature) const
     {
-      WorldBuilder::Utilities::NaturalCoordinate natural_coordinate = WorldBuilder::Utilities::NaturalCoordinate(position,
-                                                                      *(world->parameters.coordinate_system));
       if (depth <= max_depth && depth >= min_depth &&
           Utilities::polygon_contains_point(coordinates, Point<2>(natural_coordinate.get_surface_coordinates(),
                                                                   world->parameters.coordinate_system->natural_coordinate_system())))
@@ -167,12 +166,11 @@ namespace WorldBuilder
 
     double
     ContinentalPlate::composition(const Point<3> &position,
+                                  const  NaturalCoordinate &natural_coordinate,
                                   const double depth,
                                   const unsigned int composition_number,
                                   double composition) const
     {
-      WorldBuilder::Utilities::NaturalCoordinate natural_coordinate = WorldBuilder::Utilities::NaturalCoordinate(position,
-                                                                      *(world->parameters.coordinate_system));
 
       if (depth <= max_depth && depth >= min_depth &&
           Utilities::polygon_contains_point(coordinates, Point<2>(natural_coordinate.get_surface_coordinates(),
@@ -200,12 +198,11 @@ namespace WorldBuilder
 
     WorldBuilder::grains
     ContinentalPlate::grains(const Point<3> &position,
+                             const  NaturalCoordinate &natural_coordinate,
                              const double depth,
                              const unsigned int composition_number,
                              WorldBuilder::grains grains) const
     {
-      WorldBuilder::Utilities::NaturalCoordinate natural_coordinate = WorldBuilder::Utilities::NaturalCoordinate(position,
-                                                                      *(world->parameters.coordinate_system));
 
       if (depth <= max_depth && depth >= min_depth &&
           Utilities::polygon_contains_point(coordinates, Point<2>(natural_coordinate.get_surface_coordinates(),

source/features/fault.cc

@@ -329,13 +329,11 @@ namespace WorldBuilder
 
     double
     Fault::temperature(const Point<3> &position,
+                       const  NaturalCoordinate &natural_coordinate,
                        const double depth,
                        const double gravity_norm,
                        double temperature) const
     {
-      WorldBuilder::Utilities::NaturalCoordinate natural_coordinate = WorldBuilder::Utilities::NaturalCoordinate(position,
-                                                                      *(world->parameters.coordinate_system));
-
       // The depth variable is the distance from the surface to the position, the depth
       // coordinate is the distance from the bottom of the model to the position and
       // the starting radius is the distance from the bottom of the model to the surface.
@@ -356,6 +354,7 @@ namespace WorldBuilder
           // the fault to be centered around the line provided by the user.
           WorldBuilder::Utilities::PointDistanceFromCurvedPlanes distance_from_planes =
             WorldBuilder::Utilities::distance_point_from_curved_planes(position,
+                                                                       natural_coordinate,
                                                                        reference_point,
                                                                        coordinates,
                                                                        fault_segment_lengths,
@@ -470,12 +469,11 @@ namespace WorldBuilder
 
     double
     Fault::composition(const Point<3> &position,
+                       const  NaturalCoordinate &natural_coordinate,
                        const double depth,
                        const unsigned int composition_number,
                        double composition) const
     {
-      WorldBuilder::Utilities::NaturalCoordinate natural_coordinate = WorldBuilder::Utilities::NaturalCoordinate(position,
-                                                                      *(world->parameters.coordinate_system));
       // todo: explain
       const double starting_radius = natural_coordinate.get_depth_coordinate() + depth - starting_depth;
 
@@ -487,6 +485,7 @@ namespace WorldBuilder
           // the fault to be centered around the line provided by the user.
           WorldBuilder::Utilities::PointDistanceFromCurvedPlanes distance_from_planes =
             WorldBuilder::Utilities::distance_point_from_curved_planes(position,
+                                                                       natural_coordinate,
                                                                        reference_point,
                                                                        coordinates,
                                                                        fault_segment_lengths,
@@ -604,12 +603,11 @@ namespace WorldBuilder
 
     WorldBuilder::grains
     Fault::grains(const Point<3> &position,
+                  const  NaturalCoordinate &natural_coordinate,
                   const double depth,
                   const unsigned int composition_number,
                   WorldBuilder::grains grains) const
     {
-      WorldBuilder::Utilities::NaturalCoordinate natural_coordinate = WorldBuilder::Utilities::NaturalCoordinate(position,
-                                                                      *(world->parameters.coordinate_system));
       // todo: explain
       const double starting_radius = natural_coordinate.get_depth_coordinate() + depth - starting_depth;
 
@@ -621,6 +619,7 @@ namespace WorldBuilder
           // the fault to be centered around the line provided by the user.
           WorldBuilder::Utilities::PointDistanceFromCurvedPlanes distance_from_planes =
             WorldBuilder::Utilities::distance_point_from_curved_planes(position,
+                                                                       natural_coordinate,
                                                                        reference_point,
                                                                        coordinates,
                                                                        fault_segment_lengths,

source/features/mantle_layer.cc

@@ -128,13 +128,11 @@ namespace WorldBuilder
 
     double
     MantleLayer::temperature(const Point<3> &position,
+                             const  NaturalCoordinate &natural_coordinate,
                              const double depth,
                              const double gravity_norm,
                              double temperature) const
     {
-      WorldBuilder::Utilities::NaturalCoordinate natural_coordinate = WorldBuilder::Utilities::NaturalCoordinate(position,
-                                                                      *(world->parameters.coordinate_system));
-
       if (depth <= max_depth && depth >= min_depth &&
           WorldBuilder::Utilities::polygon_contains_point(coordinates, Point<2>(natural_coordinate.get_surface_coordinates(),
                                                                                 world->parameters.coordinate_system->natural_coordinate_system())))
@@ -161,13 +159,11 @@ namespace WorldBuilder
 
     double
     MantleLayer::composition(const Point<3> &position,
+                             const  NaturalCoordinate &natural_coordinate,
                              const double depth,
                              const unsigned int composition_number,
                              double composition) const
     {
-      WorldBuilder::Utilities::NaturalCoordinate natural_coordinate = WorldBuilder::Utilities::NaturalCoordinate(position,
-                                                                      *(world->parameters.coordinate_system));
-
       if (depth <= max_depth && depth >= min_depth &&
           WorldBuilder::Utilities::polygon_contains_point(coordinates, Point<2>(natural_coordinate.get_surface_coordinates(),
                                                                                 world->parameters.coordinate_system->natural_coordinate_system())))
@@ -196,13 +192,11 @@ namespace WorldBuilder
 
     WorldBuilder::grains
     MantleLayer::grains(const Point<3> &position,
+                        const  NaturalCoordinate &natural_coordinate,
                         const double depth,
                         const unsigned int composition_number,
                         WorldBuilder::grains grains) const
     {
-      WorldBuilder::Utilities::NaturalCoordinate natural_coordinate = WorldBuilder::Utilities::NaturalCoordinate(position,
-                                                                      *(world->parameters.coordinate_system));
-
       if (depth <= max_depth && depth >= min_depth &&
           WorldBuilder::Utilities::polygon_contains_point(coordinates, Point<2>(natural_coordinate.get_surface_coordinates(),
                                                                                 world->parameters.coordinate_system->natural_coordinate_system())))