Random uniform distribution deflected as rotation with tests

include/world_builder/features/continental_plate_models/grains/random_uniform_distribution_deflected.h

@@ -26,7 +26,6 @@
 
 namespace WorldBuilder
 {
-  class Parameters;
 
   namespace Features
   {
@@ -40,7 +39,7 @@ namespace WorldBuilder
          * what the returned temperature or grains of the temperature and grains
          * functions of this class will be.
          */
-        class RandomUniformDistributionDeflected final : public Interface
+        class RandomUniformDistributionDeflected final: public Interface
         {
           public:
             /**
@@ -106,6 +105,7 @@ namespace WorldBuilder
             std::vector<double> grain_sizes;
             std::vector<bool> normalize_grain_sizes;
             std::vector<double> deflections;
+            std::vector<std::array<std::array<double, 3>, 3>> basis_rotation_matrices;
         };
       } // namespace Grains
     }   // namespace ContinentalPlateModels

include/world_builder/features/fault_models/grains/random_uniform_distribution_deflected.h

@@ -104,6 +104,7 @@ namespace WorldBuilder
             std::vector<double> grain_sizes;
             std::vector<bool> normalize_grain_sizes;
             std::vector<double> deflections;
+            std::vector<std::array<std::array<double, 3>, 3>> basis_rotation_matrices;
         };
       } // namespace Grains
     }   // namespace FaultModels

include/world_builder/features/mantle_layer_models/grains/random_uniform_distribution_deflected.h

@@ -105,6 +105,7 @@ namespace WorldBuilder
             std::vector<double> grain_sizes;
             std::vector<bool> normalize_grain_sizes;
             std::vector<double> deflections;
+            std::vector<std::array<std::array<double, 3>, 3>> basis_rotation_matrices;
         };
       } // namespace Grains
     }   // namespace MantleLayerModels

include/world_builder/features/oceanic_plate_models/grains/random_uniform_distribution_deflected.h

@@ -34,7 +34,7 @@ namespace WorldBuilder
       namespace Grains
       {
         /**
-         * This class represents a continental plate and can implement
+         * This class represents an oceanic plate and can implement
          * submodules for temperature and grains. These submodules determine
          * what the returned temperature or grains of the temperature and grains
          * functions of this class will be.
@@ -105,6 +105,7 @@ namespace WorldBuilder
             std::vector<double> grain_sizes;
             std::vector<bool> normalize_grain_sizes;
             std::vector<double> deflections;
+            std::vector<std::array<std::array<double, 3>, 3>> basis_rotation_matrices;
         };
       } // namespace Grains
     }   // namespace OceanicPlateModels

include/world_builder/features/plume_models/grains/random_uniform_distribution_deflected.h

@@ -0,0 +1,112 @@
+/*
+  Copyright (C) 2018-2024 by the authors of the World Builder code.
+
+  This file is part of the World Builder.
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published
+   by the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program 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 Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef WORLD_BUILDER_FEATURES_PLUME_MODELS_GRAINS_RANDOM_UNIFORM_DISTRIBUTION_DEFLECTED_H
+#define WORLD_BUILDER_FEATURES_PLUME_MODELS_GRAINS_RANDOM_UNIFORM_DISTRIBUTION_DEFLECTED_H
+
+
+#include "world_builder/features/plume_models/grains/interface.h"
+#include "world_builder/objects/surface.h"
+
+namespace WorldBuilder
+{
+  namespace Features
+  {
+    namespace PlumeModels
+    {
+      namespace Grains
+      {
+        /**
+         * This class represents a plume and can implement
+         * submodules for temperature and grains. These submodules determine
+         * what the returned temperature or grains of the temperature and grains
+         * functions of this class will be.
+         */
+        class RandomUniformDistributionDeflected final : public Interface
+        {
+          public:
+            /**
+             * constructor
+             */
+            RandomUniformDistributionDeflected(WorldBuilder::World *world);
+
+            /**
+             * Destructor
+             */
+            ~RandomUniformDistributionDeflected() override final;
+
+            /**
+             * declare and read in the world builder file into the parameters
+             * class
+             */
+            static void declare_entries(Parameters &prm,
+                                        const std::string &parent_name = "");
+
+            /**
+             * declare and read in the world builder file into the parameters
+             * class
+             */
+            static void declare_grain_size_model_entries(
+              Parameters &prm, const std::string &parent_name,
+              const std::vector<std::string> &required_entries);
+
+            /**
+             * declare and read in the world builder file into the parameters
+             * class
+             */
+            static void
+            declare_fixed_size_model_entries(Parameters &prm,
+                                             const std::string &parent_name = "");
+
+            /**
+             * declare and read in the world builder file into the parameters
+             * class
+             */
+            void parse_entries(Parameters &prm) override final;
+
+            /**
+             * Returns a grains based on the given position, composition (e.g.
+             * olivine and/or enstatite)depth in the model, gravity and current grains.
+             */
+            WorldBuilder::grains
+            get_grains(const Point<3> &position,
+                       const Objects::NaturalCoordinate &position_in_natural_coordinates,
+                       const double depth,
+                       const unsigned int composition_number,
+                       WorldBuilder::grains grains,
+                       const double feature_min_depth,
+                       const double feature_max_depth) const override final;
+
+          private:
+            double min_depth;
+            double max_depth;
+            std::vector<unsigned int> grains;
+            std::vector<unsigned int> compositions;
+            std::string operation;
+            std::vector<double> grain_sizes;
+            std::vector<bool> normalize_grain_sizes;
+            std::vector<double> deflections;
+            std::vector<std::array<std::array<double, 3>, 3>> basis_rotation_matrices;
+        };
+      } // namespace Grains
+    }   // namespace PlumeModels
+  }     // namespace Features
+} // namespace WorldBuilder
+
+#endif

include/world_builder/features/subducting_plate_models/grains/random_uniform_distribution_deflected.h

@@ -103,6 +103,7 @@ namespace WorldBuilder
             std::vector<double> grain_sizes;
             std::vector<bool> normalize_grain_sizes;
             std::vector<double> deflections;
+            std::vector<std::array<std::array<double, 3>, 3>> basis_rotation_matrices;
 
         };
       } // namespace Grains

include/world_builder/utilities.h

@@ -487,6 +487,14 @@ namespace WorldBuilder
     std::vector<double>
     calculate_effective_trench_and_plate_ages(std::vector<double> ridge_parameters, double distance_along_plane);
 
+    /*
+     * Returns the result of the multiplication of two 3*3 matrix,
+     * used in applying the random uniform distribution rotation matrix
+     * to a given orientation (rotation matrix)
+     */
+    std::array<std::array<double,3>,3>
+    multiply_3x3_matrices(const std::array<std::array<double,3>,3> mat1, std::array<std::array<double,3>,3> const mat2);
+
   } // namespace Utilities
 } // namespace WorldBuilder
 

source/world_builder/features/continental_plate_models/grains/random_uniform_distribution_deflected.cc

@@ -19,8 +19,6 @@
 
 #include "world_builder/features/continental_plate_models/grains/random_uniform_distribution_deflected.h"
 
-#include <algorithm>
-
 #include "world_builder/nan.h"
 #include "world_builder/types/array.h"
 #include "world_builder/types/bool.h"
@@ -29,6 +27,7 @@
 #include "world_builder/types/one_of.h"
 #include "world_builder/types/unsigned_int.h"
 #include "world_builder/types/value_at_points.h"
+#include "world_builder/utilities.h"
 #include "world_builder/world.h"
 
 namespace WorldBuilder
@@ -89,7 +88,11 @@ namespace WorldBuilder
                             Types::Array(Types::Double(1),0),
                             "A list of the deflections of all of the grains in each composition between 0 and 1.");
 
+          prm.declare_entry("basis rotation matrices", Types::Array(Types::Array(Types::Array(Types::Double(0),3,3),3,3),0),
+                            "A list with the rotation matrices of the grains which are present there for each compositions.");
 
+          prm.declare_entry("basis Euler angles z-x-z", Types::Array(Types::Array(Types::Double(0),3,3),0),
+                            "A list with the z-x-z Euler angles of the grains which are present there for each compositions.");
 
         }
 
@@ -102,6 +105,31 @@ namespace WorldBuilder
           max_depth = max_depth_surface.maximum;
           compositions = prm.get_vector<unsigned int>("compositions");
 
+          const bool set_euler_angles = prm.check_entry("basis Euler angles z-x-z");
+          const bool set_rotation_matrices = prm.check_entry("basis rotation matrices");
+
+          WBAssertThrow(!(set_euler_angles == true && set_rotation_matrices == true),
+                        "Only Euler angles or Rotation matrices may be set, but both are set for " << prm.get_full_json_path());
+
+
+          WBAssertThrow(!(set_euler_angles == false && set_rotation_matrices == false),
+                        "Euler angles or Rotation matrices have to be set, but neither are set for " << prm.get_full_json_path());
+
+          if (set_euler_angles)
+            {
+              std::vector<std::array<double,3> > basis_euler_angles_vector = prm.get_vector<std::array<double,3> >("basis Euler angles z-x-z");
+              basis_rotation_matrices.resize(basis_euler_angles_vector.size());
+              for (size_t i = 0; i<basis_euler_angles_vector.size(); ++i)
+                {
+                  basis_rotation_matrices[i] = Utilities::euler_angles_to_rotation_matrix(basis_euler_angles_vector[i][0],basis_euler_angles_vector[i][1],basis_euler_angles_vector[i][2]);
+                }
+
+            }
+          else
+            {
+              basis_rotation_matrices = prm.get_vector<std::array<std::array<double,3>,3> >("basis rotation matrices");
+            }
+
           operation = prm.get<std::string>("orientation operation");
           grain_sizes = prm.get_vector<double>("grain sizes");
           normalize_grain_sizes = prm.get_vector<bool>("normalize grain sizes");
@@ -116,9 +144,11 @@ namespace WorldBuilder
           WBAssertThrow(compositions.size() == deflections.size(),
                         "There are not the same amount of compositions (" << compositions.size()
                         << ") and deflections (" << deflections.size() << ").");
+          WBAssertThrow(compositions.size() == basis_rotation_matrices.size(),
+                        "There are not the same amount of compositions (" << compositions.size()
+                        << ") and rotation_matrices (" << basis_rotation_matrices.size() << ").");
         }
 
-
         WorldBuilder::grains
         RandomUniformDistributionDeflected::get_grains(const Point<3> & /*position_in_cartesian_coordinates*/,
                                                        const Objects::NaturalCoordinate &position_in_natural_coordinates,
@@ -161,9 +191,10 @@ namespace WorldBuilder
                               const double two = dist(world->get_random_number_engine());
                               const double three = dist(world->get_random_number_engine());
 
+                              // the distribution is restricted by the deflection (between 0 and 1)
                               const double theta = 2.0 * Consts::PI * one * deflections[i]; // Rotation about the pole (Z)
                               const double phi = 2.0 * Consts::PI * two; // For direction of pole deflection.
-                              const double z = 2.0* three * deflections[i]; //For magnitude of pole deflection.
+                              const double z = 2.0 * three * deflections[i]; //For magnitude of pole deflection.
 
                               // Compute a vector V used for distributing points over the sphere
                               // via the reflection I - V Transpose(V).  This formulation of V
@@ -188,17 +219,33 @@ namespace WorldBuilder
                               // Construct the rotation matrix  ( V Transpose(V) - I ) R, which
                               // is equivalent to V S - R.
 
-                              it_rotation_matrices[0][0] = Vx * Sx - ct;
-                              it_rotation_matrices[0][1] = Vx * Sy - st;
-                              it_rotation_matrices[0][2] = Vx * Vz;
-
-                              it_rotation_matrices[1][0] = Vy * Sx + st;
-                              it_rotation_matrices[1][1] = Vy * Sy - ct;
-                              it_rotation_matrices[1][2] = Vy * Vz;
-
-                              it_rotation_matrices[2][0] = Vz * Sx;
-                              it_rotation_matrices[2][1] = Vz * Sy;
-                              it_rotation_matrices[2][2] = 1.0 - z;   // This equals Vz * Vz - 1.0
+                              std::array<std::array<double,3>,3> rotation_matrices;
+                              rotation_matrices[0][0] = (Vx * Sx - ct);
+                              rotation_matrices[0][1] = (Vx * Sy - st);
+                              rotation_matrices[0][2] = Vx * Vz;
+
+                              rotation_matrices[1][0] = (Vy * Sx + st);
+                              rotation_matrices[1][1] = (Vy * Sy - ct);
+                              rotation_matrices[1][2] = Vy * Vz;
+
+                              rotation_matrices[2][0] = Vz * Sx;
+                              rotation_matrices[2][1] = Vz * Sy;
+                              rotation_matrices[2][2] = 1.0 - z;   // This equals Vz * Vz - 1.0
+
+                              // Rotate the basis rotation matrix with the random uniform distribution rotation matrix
+                              // First get the transpose of the rotation matrix
+                              std::array<std::array<double, 3>, 3> rot_T= rotation_matrices;
+                              rot_T[0][1] = rotation_matrices[1][0];
+                              rot_T[1][0] = rotation_matrices[0][1];
+                              rot_T[1][2] = rotation_matrices[2][1];
+                              rot_T[2][1] = rotation_matrices[1][2];
+                              rot_T[0][2] = rotation_matrices[2][0];
+                              rot_T[2][0] = rotation_matrices[0][2];
+
+                              // Then U' = R * U * R^T
+                              std::array<std::array<double,3>,3> result1 = Utilities::multiply_3x3_matrices(rotation_matrices, basis_rotation_matrices[i]);
+
+                              it_rotation_matrices = result1;
                             }
 
                           double total_size = 0;
@@ -211,8 +258,15 @@ namespace WorldBuilder
                           if (normalize_grain_sizes[i])
                             {
                               const double one_over_total_size = 1/total_size;
-                              std::transform(grains_local.sizes.begin(), grains_local.sizes.end(), grains_local.sizes.begin(),
-                                             [one_over_total_size](double sizes) -> double { return sizes *one_over_total_size; });
+                              // std::transform is a c++17 feature, while current GWB is c++14
+                              // update this after switching to c+=17
+                              // std::transform(grains_local.sizes.begin(), grains_local.sizes.end(), grains_local.sizes.begin(),
+                              //                [one_over_total_size](double sizes) -> double { return sizes *one_over_total_size; });
+                              // Apply the transformation using a loop
+                              for (auto &&size : grains_local.sizes)
+                                {
+                                  size = size*one_over_total_size;
+                                }
                             }