diff --git a/CHANGELOG.md b/CHANGELOG.md index 6402487d6..d94bbd8b3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,7 +11,7 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm ## [Unreleased] ### Added - +- Added a `random uniform distribution deflected` grains model for all features that allows an initial texture computed from a random uniform distribution of rotation matrices applied to a given orientation specified as a set of Euler angles or a rotation matrix. \[Yijun Wang; 2024-06-06; [#713](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/713)\] ### Changed - In the "mass conserving" model, change the name of the entry "plate velocity" to "spreading velocity" \[Haoyuan Li; 2024-03-11; [#694](https://github.com/GeodynamicWorldBuilder/WorldBuilder/pull/694)\] diff --git a/include/world_builder/features/continental_plate_models/grains/random_uniform_distribution_deflected.h b/include/world_builder/features/continental_plate_models/grains/random_uniform_distribution_deflected.h index 4e8cf7173..54c389f25 100644 --- a/include/world_builder/features/continental_plate_models/grains/random_uniform_distribution_deflected.h +++ b/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 grain_sizes; std::vector normalize_grain_sizes; std::vector deflections; + std::vector, 3>> basis_rotation_matrices; }; } // namespace Grains } // namespace ContinentalPlateModels diff --git a/include/world_builder/features/fault_models/grains/random_uniform_distribution_deflected.h b/include/world_builder/features/fault_models/grains/random_uniform_distribution_deflected.h index 82d2c683f..b8ee7b15c 100644 --- a/include/world_builder/features/fault_models/grains/random_uniform_distribution_deflected.h +++ b/include/world_builder/features/fault_models/grains/random_uniform_distribution_deflected.h @@ -104,6 +104,7 @@ namespace WorldBuilder std::vector grain_sizes; std::vector normalize_grain_sizes; std::vector deflections; + std::vector, 3>> basis_rotation_matrices; }; } // namespace Grains } // namespace FaultModels diff --git a/include/world_builder/features/mantle_layer_models/grains/random_uniform_distribution_deflected.h b/include/world_builder/features/mantle_layer_models/grains/random_uniform_distribution_deflected.h index 4d0bd807b..5b28b4bce 100644 --- a/include/world_builder/features/mantle_layer_models/grains/random_uniform_distribution_deflected.h +++ b/include/world_builder/features/mantle_layer_models/grains/random_uniform_distribution_deflected.h @@ -105,6 +105,7 @@ namespace WorldBuilder std::vector grain_sizes; std::vector normalize_grain_sizes; std::vector deflections; + std::vector, 3>> basis_rotation_matrices; }; } // namespace Grains } // namespace MantleLayerModels diff --git a/include/world_builder/features/oceanic_plate_models/grains/random_uniform_distribution_deflected.h b/include/world_builder/features/oceanic_plate_models/grains/random_uniform_distribution_deflected.h index 4a8094581..f9b9c408a 100644 --- a/include/world_builder/features/oceanic_plate_models/grains/random_uniform_distribution_deflected.h +++ b/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 grain_sizes; std::vector normalize_grain_sizes; std::vector deflections; + std::vector, 3>> basis_rotation_matrices; }; } // namespace Grains } // namespace OceanicPlateModels diff --git a/include/world_builder/features/plume_models/grains/random_uniform_distribution_deflected.h b/include/world_builder/features/plume_models/grains/random_uniform_distribution_deflected.h new file mode 100644 index 000000000..da139e048 --- /dev/null +++ b/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 . +*/ + +#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 &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 grains; + std::vector compositions; + std::string operation; + std::vector grain_sizes; + std::vector normalize_grain_sizes; + std::vector deflections; + std::vector, 3>> basis_rotation_matrices; + }; + } // namespace Grains + } // namespace PlumeModels + } // namespace Features +} // namespace WorldBuilder + +#endif diff --git a/include/world_builder/features/subducting_plate_models/grains/random_uniform_distribution_deflected.h b/include/world_builder/features/subducting_plate_models/grains/random_uniform_distribution_deflected.h index ddaada0be..273764c13 100644 --- a/include/world_builder/features/subducting_plate_models/grains/random_uniform_distribution_deflected.h +++ b/include/world_builder/features/subducting_plate_models/grains/random_uniform_distribution_deflected.h @@ -103,6 +103,7 @@ namespace WorldBuilder std::vector grain_sizes; std::vector normalize_grain_sizes; std::vector deflections; + std::vector, 3>> basis_rotation_matrices; }; } // namespace Grains diff --git a/include/world_builder/utilities.h b/include/world_builder/utilities.h index fff53a6bc..0da073ea2 100644 --- a/include/world_builder/utilities.h +++ b/include/world_builder/utilities.h @@ -487,6 +487,14 @@ namespace WorldBuilder std::vector calculate_effective_trench_and_plate_ages(std::vector 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,3> + multiply_3x3_matrices(const std::array,3> mat1, std::array,3> const mat2); + } // namespace Utilities } // namespace WorldBuilder diff --git a/source/world_builder/features/continental_plate_models/grains/random_uniform_distribution_deflected.cc b/source/world_builder/features/continental_plate_models/grains/random_uniform_distribution_deflected.cc index 02150f29e..d1b077f17 100644 --- a/source/world_builder/features/continental_plate_models/grains/random_uniform_distribution_deflected.cc +++ b/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 - #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("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 > basis_euler_angles_vector = prm.get_vector >("basis Euler angles z-x-z"); + basis_rotation_matrices.resize(basis_euler_angles_vector.size()); + for (size_t i = 0; i,3> >("basis rotation matrices"); + } + operation = prm.get("orientation operation"); grain_sizes = prm.get_vector("grain sizes"); normalize_grain_sizes = prm.get_vector("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,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, 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,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; + } } diff --git a/source/world_builder/features/fault_models/grains/random_uniform_distribution_deflected.cc b/source/world_builder/features/fault_models/grains/random_uniform_distribution_deflected.cc index 1f970ee20..d2d10249d 100644 --- a/source/world_builder/features/fault_models/grains/random_uniform_distribution_deflected.cc +++ b/source/world_builder/features/fault_models/grains/random_uniform_distribution_deflected.cc @@ -19,8 +19,6 @@ #include "world_builder/features/fault_models/grains/random_uniform_distribution_deflected.h" -#include - #include "world_builder/nan.h" #include "world_builder/types/array.h" #include "world_builder/types/bool.h" @@ -88,6 +86,12 @@ 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."); + } @@ -99,6 +103,31 @@ namespace WorldBuilder max_depth = prm.get("max distance fault center"); compositions = prm.get_vector("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 > basis_euler_angles_vector = prm.get_vector >("basis Euler angles z-x-z"); + basis_rotation_matrices.resize(basis_euler_angles_vector.size()); + for (size_t i = 0; i,3> >("basis rotation matrices"); + } + operation = prm.get("orientation operation"); grain_sizes = prm.get_vector("grain sizes"); normalize_grain_sizes = prm.get_vector("normalize grain sizes"); @@ -113,9 +142,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 double /*depth*/, @@ -182,17 +213,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,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, 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,3> result1 = multiply_3x3_matrices(rotation_matrices, basis_rotation_matrices[i]); + + it_rotation_matrices = result1; } double total_size = 0; @@ -205,8 +252,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; + } } return grains_local; diff --git a/source/world_builder/features/mantle_layer_models/grains/random_uniform_distribution_deflected.cc b/source/world_builder/features/mantle_layer_models/grains/random_uniform_distribution_deflected.cc index d3cf85005..deb0dfa72 100644 --- a/source/world_builder/features/mantle_layer_models/grains/random_uniform_distribution_deflected.cc +++ b/source/world_builder/features/mantle_layer_models/grains/random_uniform_distribution_deflected.cc @@ -19,8 +19,6 @@ #include "world_builder/features/mantle_layer_models/grains/random_uniform_distribution_deflected.h" -#include - #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("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 > basis_euler_angles_vector = prm.get_vector >("basis Euler angles z-x-z"); + basis_rotation_matrices.resize(basis_euler_angles_vector.size()); + for (size_t i = 0; i,3> >("basis rotation matrices"); + } + operation = prm.get("orientation operation"); grain_sizes = prm.get_vector("grain sizes"); normalize_grain_sizes = prm.get_vector("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, @@ -164,7 +194,7 @@ namespace WorldBuilder // 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 @@ -189,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,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, 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,3> result1 = multiply_3x3_matrices(rotation_matrices, basis_rotation_matrices[i]); + + it_rotation_matrices = result1; } double total_size = 0; @@ -212,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; + } } return grains_local; diff --git a/source/world_builder/features/oceanic_plate_models/grains/random_uniform_distribution_deflected.cc b/source/world_builder/features/oceanic_plate_models/grains/random_uniform_distribution_deflected.cc index ffac67363..e95a4cf56 100644 --- a/source/world_builder/features/oceanic_plate_models/grains/random_uniform_distribution_deflected.cc +++ b/source/world_builder/features/oceanic_plate_models/grains/random_uniform_distribution_deflected.cc @@ -19,8 +19,6 @@ #include "world_builder/features/oceanic_plate_models/grains/random_uniform_distribution_deflected.h" -#include - #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("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 > basis_euler_angles_vector = prm.get_vector >("basis Euler angles z-x-z"); + basis_rotation_matrices.resize(basis_euler_angles_vector.size()); + for (size_t i = 0; i,3> >("basis rotation matrices"); + } + operation = prm.get("orientation operation"); grain_sizes = prm.get_vector("grain sizes"); normalize_grain_sizes = prm.get_vector("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,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, 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,3> result1 = 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; + } } return grains_local; diff --git a/source/world_builder/features/plume_models/grains/random_uniform_distribution_deflected.cc b/source/world_builder/features/plume_models/grains/random_uniform_distribution_deflected.cc new file mode 100644 index 000000000..c1b98b6db --- /dev/null +++ b/source/world_builder/features/plume_models/grains/random_uniform_distribution_deflected.cc @@ -0,0 +1,282 @@ +/* + 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 . +*/ + +#include "world_builder/features/plume_models/grains/random_uniform_distribution_deflected.h" + +#include "world_builder/nan.h" +#include "world_builder/types/array.h" +#include "world_builder/types/double.h" +#include "world_builder/types/object.h" +#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/types/bool.h" +#include "world_builder/world.h" + +namespace WorldBuilder +{ + + using namespace Utilities; + + namespace Features + { + namespace PlumeModels + { + namespace Grains + { + RandomUniformDistributionDeflected::RandomUniformDistributionDeflected(WorldBuilder::World *world_) + : + min_depth(NaN::DSNAN), + max_depth(NaN::DSNAN) + + { + this->world = world_; + this->name = "random uniform distribution deflected"; + } + + RandomUniformDistributionDeflected::~RandomUniformDistributionDeflected() + = default; + + void + RandomUniformDistributionDeflected::declare_entries(Parameters &prm, const std::string & /*unused*/) + { + // Document plugin and require entries if needed. + // Add compositions to the required parameters. + prm.declare_entry("", Types::Object({"compositions"}), + "Random uniform distribution grains model. The size of the grains can be independently set " + "to a single value or to a random distribution."); + + // Declare entries of this plugin + prm.declare_entry("min depth", Types::Double(0), + "The depth in meters from which the grains of this feature are present."); + prm.declare_entry("max depth", Types::Double(std::numeric_limits::max()), + "The depth in meters to which the grains of this feature are present."); + + prm.declare_entry("compositions", Types::Array(Types::UnsignedInt(),0), + "A list with the integer labels of the composition which are present there."); + + prm.declare_entry("orientation operation", Types::String("replace", std::vector {"replace", "multiply"}), + "Whether the value should replace any value previously defined at this location (replace) or " + "add the value to the previously define value (add, not implemented). Replacing implies that all values not " + "explicitly defined are set to zero."); + + prm.declare_entry("grain sizes", + Types::Array(Types::Double(1),0), + "A list of the size of all of the grains in each composition. If set to <0, the size will be randomized between 0 and 1."); + + prm.declare_entry("normalize grain sizes", + Types::Array(Types::Bool(true),0), + "A list of whether the sizes of the grains should be normalized or not. If normalized, the total of the grains of a composition will be equal to 1."); + + prm.declare_entry("deflections", + 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."); + + + } + + void + RandomUniformDistributionDeflected::parse_entries(Parameters &prm) + { + min_depth = prm.get("min depth"); + max_depth = prm.get("max depth"); + compositions = prm.get_vector("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 > basis_euler_angles_vector = prm.get_vector >("basis Euler angles z-x-z"); + basis_rotation_matrices.resize(basis_euler_angles_vector.size()); + for (size_t i = 0; i,3> >("basis rotation matrices"); + } + + operation = prm.get("orientation operation"); + grain_sizes = prm.get_vector("grain sizes"); + normalize_grain_sizes = prm.get_vector("normalize grain sizes"); + deflections = prm.get_vector("deflections"); + + + WBAssertThrow(compositions.size() == grain_sizes.size(), + "There are not the same amount of compositions (" << compositions.size() + << ") and grain_sizes (" << grain_sizes.size() << ")."); + WBAssertThrow(compositions.size() == normalize_grain_sizes.size(), + "There are not the same amount of compositions (" << compositions.size() + << ") and normalize_grain_sizes (" << normalize_grain_sizes.size() << ")."); + 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*/, + const double depth, + const unsigned int composition_number, + WorldBuilder::grains grains_, + const double /*feature_min_depth*/, + const double /*feature_max_depth*/) const + { + WorldBuilder::grains grains_local = grains_; + if (depth <= max_depth && depth >= min_depth) + { + for (unsigned int i =0; i < compositions.size(); ++i) + { + if (compositions[i] == composition_number) + { + std::uniform_real_distribution<> dist(0.0,1.0); + for (auto &&it_rotation_matrices : grains_local.rotation_matrices) + { + // set a uniform random a_cosine_matrix per grain + // This function is based on an article in Graphic Gems III, written by James Arvo, Cornell University (p 116-120). + // The original code can be found on http://www.realtimerendering.com/resources/GraphicsGems/gemsiii/rand_rotation.c + // and is licenend accourding to this website with the following licence: + // + // "The Graphics Gems code is copyright-protected. In other words, you cannot claim the text of the code as your own and + // resell it. Using the code is permitted in any program, product, or library, non-commercial or commercial. Giving credit + // is not required, though is a nice gesture. The code comes as-is, and if there are any flaws or problems with any Gems + // code, nobody involved with Gems - authors, editors, publishers, or webmasters - are to be held responsible. Basically, + // don't be a jerk, and remember that anything free comes with no guarantee."" + // + // The book saids in the preface the following: "As in the first two volumes, all of the C and C++ code in this book is in + // the public domain, and is yours to study, modify, and use." + + // first generate three random numbers between 0 and 1 and multiply them with 2 PI or 2 for z. Note that these are not the same as phi_1, theta and phi_2. + const double one = dist(world->get_random_number_engine()); + 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. + + // Compute a vector V used for distributing points over the sphere + // via the reflection I - V Transpose(V). This formulation of V + // will guarantee that if x[1] and x[2] are uniformly distributed, + // the reflected points will be uniform on the sphere. Note that V + // has length sqrt(2) to eliminate the 2 in the Householder matrix. + + const double r = std::sqrt( z ); + const double Vx = std::sin( phi ) * r; + const double Vy = std::cos( phi ) * r; + const double Vz = std::sqrt( 2.F - z ); + + // Compute the row vector S = Transpose(V) * R, where R is a simple + // rotation by theta about the z-axis. No need to compute Sz since + // it's just Vz. + + const double st = std::sin( theta ); + const double ct = std::cos( theta ); + const double Sx = Vx * ct - Vy * st; + const double Sy = Vx * st + Vy * ct; + + // Construct the rotation matrix ( V Transpose(V) - I ) R, which + // is equivalent to V S - R. + + std::array,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, 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,3> result1 = multiply_3x3_matrices(rotation_matrices, basis_rotation_matrices[i]); + + it_rotation_matrices = result1; + } + + double total_size = 0; + for (auto &&it_sizes : grains_local.sizes) + { + it_sizes = grain_sizes[i] < 0 ? dist(world->get_random_number_engine()) : grain_sizes[i]; + total_size += it_sizes; + } + + if (normalize_grain_sizes[i]) + { + const double one_over_total_size = 1/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; + } + } + + + return grains_local; + } + } + } + return grains_local; + } + WB_REGISTER_FEATURE_PLUME_GRAINS_MODEL(RandomUniformDistributionDeflected, random uniform distribution deflected) + } // namespace Grains + } // namespace PlumeModels + } // namespace Features +} // namespace WorldBuilder + + diff --git a/source/world_builder/features/subducting_plate_models/grains/random_uniform_distribution_deflected.cc b/source/world_builder/features/subducting_plate_models/grains/random_uniform_distribution_deflected.cc index bf9ac3a3e..dfbf267af 100644 --- a/source/world_builder/features/subducting_plate_models/grains/random_uniform_distribution_deflected.cc +++ b/source/world_builder/features/subducting_plate_models/grains/random_uniform_distribution_deflected.cc @@ -19,8 +19,6 @@ #include "world_builder/features/subducting_plate_models/grains/random_uniform_distribution_deflected.h" -#include - #include "world_builder/nan.h" #include "world_builder/types/array.h" #include "world_builder/types/bool.h" @@ -88,6 +86,12 @@ 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."); + } @@ -99,6 +103,31 @@ namespace WorldBuilder max_depth = prm.get("max distance slab top"); compositions = prm.get_vector("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 > basis_euler_angles_vector = prm.get_vector >("basis Euler angles z-x-z"); + basis_rotation_matrices.resize(basis_euler_angles_vector.size()); + for (size_t i = 0; i,3> >("basis rotation matrices"); + } + operation = prm.get("orientation operation"); grain_sizes = prm.get_vector("grain sizes"); normalize_grain_sizes = prm.get_vector("normalize grain sizes"); @@ -113,9 +142,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 double /*depth*/, @@ -182,17 +213,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,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, 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,3> result1 = multiply_3x3_matrices(rotation_matrices, basis_rotation_matrices[i]); + + it_rotation_matrices = result1; } double total_size = 0; @@ -205,8 +252,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; + } } return grains_local; diff --git a/source/world_builder/utilities.cc b/source/world_builder/utilities.cc index eadcaf66d..177f13c17 100644 --- a/source/world_builder/utilities.cc +++ b/source/world_builder/utilities.cc @@ -1505,6 +1505,25 @@ namespace WorldBuilder return result; } + + std::array,3> + multiply_3x3_matrices(const std::array,3> mat1, const std::array,3> mat2) + { + std::array,3> result; + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + result[i][j] = 0; + for (int k = 0; k < 3; k++) + { + result[i][j] += mat1[i][k] * mat2[k][j]; + } + } + } + + return result; + } } // namespace Utilities } // namespace WorldBuilder diff --git a/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1.dat b/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1.dat new file mode 100644 index 000000000..99dd6303d --- /dev/null +++ b/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1.dat @@ -0,0 +1,14 @@ +# This is a comment in the data +# file. +# Now define parameters: +# dim = 3 +# compositions = 6 +# grain compositions = 2 +# number of grains = 1 +# x y z d c1 c2 c3 c4 c5 c6 +5 5 0 15 +1 5 0 5 +7 5 0 5 +5 3 0 0.1 +4.5 0 0 15 +5 6 0 0.1 \ No newline at end of file diff --git a/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1.wb b/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1.wb new file mode 100644 index 000000000..0cd87ce47 --- /dev/null +++ b/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1.wb @@ -0,0 +1,80 @@ +{ + "version": "1.0", + "coordinate system":{"model":"cartesian"}, + "random number seed":1, + "features": + [ + { + "model":"mantle layer", "name":"Mantle", "min depth":10, "max depth": 20, + "coordinates":[[0,0],[0,10],[10,10],[10,0]], + "temperature models":[{"model":"uniform", "temperature":293}], + "composition models":[{"model":"uniform", "compositions":[0]}], + "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis Euler angles z-x-z": [[0,0,0]]}, + {"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis rotation matrices": [[[-0.5,0,0],[0,1,0],[0,0,-1]]]}] + }, + + { + "model":"continental plate", "name":"CP", "min depth":0, "max depth": 10, + "coordinates":[[0,0],[0,10],[5,10],[5,0]], + "temperature models":[{"model":"uniform", "temperature":293}], + "composition models":[{"model":"uniform", "compositions":[1]}], + "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis Euler angles z-x-z": [[0,30,0]]}, + {"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis rotation matrices": [[[-0.5,0,0],[0,1,0],[0,0,-1]]]}] + }, + + { + "model":"oceanic plate", "name":"OP", "min depth":0, "max depth": 10, + "coordinates":[[5,0],[5,10],[10,10],[10,0]], + "temperature models":[{"model":"uniform", "temperature":293}], + "composition models":[{"model":"uniform", "compositions":[2]}], + "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis Euler angles z-x-z": [[0,0,45]]}, + {"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis rotation matrices": [[[-0.5,0,0],[0,1,0],[0,0,-1]]]}] + }, + + { + "model":"fault", "name":"Fault", "min depth":0, "max depth": 10, + "coordinates":[[0,3],[10,3]], "dip point":[0,0], + "segments":[{"length":10, "thickness":[1], "angle":[45]}], + "temperature models":[{"model":"uniform", "temperature":293}], + "composition models":[{"model":"uniform", "compositions":[3]}], + "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis Euler angles z-x-z": [[60,60,60]]}, + {"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis rotation matrices": [[[-0.5,0,0],[0,1,0],[0,0,-1]]]}] + }, + + { + "model":"plume", + "name":"Plume", + "coordinates":[[3, 0],[3.7, 0],[5.5, 0],[4.5, 0],[4, 0]], + "cross section depths":[0, 1, 2.5, 5, 7.5], + "semi-major axis":[0.5, 0.7, 0.3, 0.8, 1], + "eccentricity":[0, 0, 0, 0, 0], + "rotation angles":[0, 0, 0, 0, 0], + "temperature models":[{"model":"uniform", "min depth": 10, "max depth":20, "temperature":1800}], + "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.5], "normalize grain sizes": [true], "basis Euler angles z-x-z": [[0,90,90]]}, + {"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis rotation matrices": [[[-0.5,0,0],[0,1,0],[0,0,-1]]]}] + }, + + { + "model":"subducting plate", "name":"SP", "min depth":0, "max depth": 10, + "coordinates":[[0,6],[10,6]], "dip point":[0,0], + "segments":[{"length":10, "thickness":[1], "angle":[45]}], + "temperature models":[{"model":"uniform", "temperature":150}], + "composition models":[{"model":"uniform", "compositions":[5]}], + "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis Euler angles z-x-z": [[30,30,60]]}, + {"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis rotation matrices": [[[-0.5,0,0],[0,1,0],[0,0,-1]]]}] + } + + ] +} \ No newline at end of file diff --git a/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1/screen-output.log b/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1/screen-output.log new file mode 100644 index 000000000..4c039fb16 --- /dev/null +++ b/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1/screen-output.log @@ -0,0 +1,7 @@ +# x y z d g T c0 c1 c2 c3 c4 c5 gs0-0 gm0-0[0:0] gm0-0[0:1] gm0-0[0:2] gm0-0[1:0] gm0-0[1:1] gm0-0[1:2] gm0-0[2:0] gm0-0[2:1] gm0-0[2:2] gs1-0 gm1-0[0:0] gm1-0[0:1] gm1-0[0:2] gm1-0[1:0] gm1-0[1:1] gm1-0[1:2] gm1-0[2:0] gm1-0[2:1] gm1-0[2:2] tag +5 5 0 15 293 1 0 0 0 0 0 1 0.499857 -0.0018567 -0.023828 -0.000554464 -0.999508 0.0313571 -0.0119373 -0.0313217 -0.999224 0 0 0 0 0 0 0 0 0 0 0 +1 5 0 5 293 0 1 0 0 0 0 1 0.499473 -0.00264771 -0.0458424 -0.00146095 -0.999978 -0.005915 -0.0229129 0.00604271 -0.998931 0 0 0 0 0 0 0 0 0 0 1 +7 5 0 5 293 0 0 1 0 0 0 1 0.499912 -0.00455493 0.0182169 -0.00222068 -0.99997 -0.00626999 0.00912247 0.00618798 -0.999814 0 0 0 0 0 0 0 0 0 0 2 +5 3 0 0.1 293 0 0 0 1 0 0 1 0.999529 -0.0039371 -0.0290701 -0.00342152 -0.825622 0.0162258 -0.0291353 -0.0161085 -0.825165 0 1 0 0 0 1 0 0 0 1 3 +4.5 0 0 15 1800 1 0 0 0 0 0 1 0.499757 -0.000725128 0.0311578 0.000175418 -0.999404 -0.0345134 0.0155821 0.0345076 -0.998918 0 0 0 0 0 0 0 0 0 0 4 +5 6 0 0.1 150 0 0 0 0 0 1 1 0.999996 -0.00262365 0.000402509 -0.0025966 -0.824507 -0.0517783 0.000550654 0.0517769 -0.824511 0 1 0 0 0 1 0 0 0 1 5 diff --git a/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1000.dat b/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1000.dat new file mode 100644 index 000000000..99dd6303d --- /dev/null +++ b/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1000.dat @@ -0,0 +1,14 @@ +# This is a comment in the data +# file. +# Now define parameters: +# dim = 3 +# compositions = 6 +# grain compositions = 2 +# number of grains = 1 +# x y z d c1 c2 c3 c4 c5 c6 +5 5 0 15 +1 5 0 5 +7 5 0 5 +5 3 0 0.1 +4.5 0 0 15 +5 6 0 0.1 \ No newline at end of file diff --git a/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1000.wb b/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1000.wb new file mode 100644 index 000000000..479bce675 --- /dev/null +++ b/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1000.wb @@ -0,0 +1,80 @@ +{ + "version": "1.0", + "coordinate system":{"model":"cartesian"}, + "random number seed":1000, + "features": + [ + { + "model":"mantle layer", "name":"Mantle", "min depth":10, "max depth": 20, + "coordinates":[[0,0],[0,10],[10,10],[10,0]], + "temperature models":[{"model":"uniform", "temperature":293}], + "composition models":[{"model":"uniform", "compositions":[0]}], + "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis Euler angles z-x-z": [[0,0,0]]}, + {"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis rotation matrices": [[[-0.5,0,0],[0,1,0],[0,0,-1]]]}] + }, + + { + "model":"continental plate", "name":"CP", "min depth":0, "max depth": 10, + "coordinates":[[0,0],[0,10],[5,10],[5,0]], + "temperature models":[{"model":"uniform", "temperature":293}], + "composition models":[{"model":"uniform", "compositions":[1]}], + "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis Euler angles z-x-z": [[0,30,0]]}, + {"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis rotation matrices": [[[-0.5,0,0],[0,1,0],[0,0,-1]]]}] + }, + + { + "model":"oceanic plate", "name":"OP", "min depth":0, "max depth": 10, + "coordinates":[[5,0],[5,10],[10,10],[10,0]], + "temperature models":[{"model":"uniform", "temperature":293}], + "composition models":[{"model":"uniform", "compositions":[2]}], + "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis Euler angles z-x-z": [[0,0,45]]}, + {"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis rotation matrices": [[[-0.5,0,0],[0,1,0],[0,0,-1]]]}] + }, + + { + "model":"fault", "name":"Fault", "min depth":0, "max depth": 10, + "coordinates":[[0,3],[10,3]], "dip point":[0,0], + "segments":[{"length":10, "thickness":[1], "angle":[45]}], + "temperature models":[{"model":"uniform", "temperature":293}], + "composition models":[{"model":"uniform", "compositions":[3]}], + "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis Euler angles z-x-z": [[60,60,60]]}, + {"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis rotation matrices": [[[-0.5,0,0],[0,1,0],[0,0,-1]]]}] + }, + + { + "model":"plume", + "name":"Plume", + "coordinates":[[3, 0],[3.7, 0],[5.5, 0],[4.5, 0],[4, 0]], + "cross section depths":[0, 1, 2.5, 5, 7.5], + "semi-major axis":[0.5, 0.7, 0.3, 0.8, 1], + "eccentricity":[0, 0, 0, 0, 0], + "rotation angles":[0, 0, 0, 0, 0], + "temperature models":[{"model":"uniform", "min depth": 10, "max depth":20, "temperature":1800}], + "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.5], "normalize grain sizes": [true], "basis Euler angles z-x-z": [[0,90,90]]}, + {"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis rotation matrices": [[[-0.5,0,0],[0,1,0],[0,0,-1]]]}] + }, + + { + "model":"subducting plate", "name":"SP", "min depth":0, "max depth": 10, + "coordinates":[[0,6],[10,6]], "dip point":[0,0], + "segments":[{"length":10, "thickness":[1], "angle":[45]}], + "temperature models":[{"model":"uniform", "temperature":150}], + "composition models":[{"model":"uniform", "compositions":[5]}], + "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis Euler angles z-x-z": [[30,30,60]]}, + {"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], + "deflections": [0.001], "normalize grain sizes": [true], "basis rotation matrices": [[[-0.5,0,0],[0,1,0],[0,0,-1]]]}] + } + + ] +} \ No newline at end of file diff --git a/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1000/screen-output.log b/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1000/screen-output.log new file mode 100644 index 000000000..7956bb5c0 --- /dev/null +++ b/tests/gwb-dat/random_uniform_distribution_deflected_asrotation_seed1000/screen-output.log @@ -0,0 +1,7 @@ +# x y z d g T c0 c1 c2 c3 c4 c5 gs0-0 gm0-0[0:0] gm0-0[0:1] gm0-0[0:2] gm0-0[1:0] gm0-0[1:1] gm0-0[1:2] gm0-0[2:0] gm0-0[2:1] gm0-0[2:2] gs1-0 gm1-0[0:0] gm1-0[0:1] gm1-0[0:2] gm1-0[1:0] gm1-0[1:1] gm1-0[1:2] gm1-0[2:0] gm1-0[2:1] gm1-0[2:2] tag +5 5 0 15 293 1 0 0 0 0 0 1 0.499991 -0.00528206 -0.00316564 -0.00269775 -0.999302 -0.0369657 -0.00148409 0.0369821 -0.999312 0 0 0 0 0 0 0 0 0 0 0 +1 5 0 5 293 0 1 0 0 0 0 1 0.499984 -0.000140805 -0.00798551 -0.000308283 -0.998222 -0.0596069 -0.00398146 0.0596099 -0.99819 0 0 0 0 0 0 0 0 0 0 1 +7 5 0 5 293 0 0 1 0 0 0 1 0.499796 -0.00317139 0.0283707 -0.00140799 -0.999917 -0.0125586 0.0142041 0.0124735 -0.999519 0 0 0 0 0 0 0 0 0 0 2 +5 3 0 0.1 293 0 0 0 1 0 0 1 0.999858 -0.00458949 -0.0154063 -0.00475835 -0.825733 -0.00998195 -0.015355 0.0100607 -0.825615 0 1 0 0 0 1 0 0 0 1 3 +4.5 0 0 15 1800 1 0 0 0 0 0 1 0.499936 -0.00455961 -0.0152841 -0.00258492 -0.999181 -0.0401272 -0.00754431 0.0402011 -0.999078 0 0 0 0 0 0 0 0 0 0 4 +5 6 0 0.1 150 0 0 0 0 0 1 1 0.999991 -0.00293585 0.00292161 -0.0029854 -0.8257 0.0156096 0.00287097 -0.015619 -0.8257 0 1 0 0 0 1 0 0 0 1 5 diff --git a/tests/gwb-dat/random_uniform_texture_all_features_seed1.dat b/tests/gwb-dat/random_uniform_texture_all_features_seed1.dat deleted file mode 100644 index d0c093e2e..000000000 --- a/tests/gwb-dat/random_uniform_texture_all_features_seed1.dat +++ /dev/null @@ -1,20 +0,0 @@ -# This is a comment in the data -# file. -# Now define parameters: -# dim = 3 -# compositions = 5 -# grain compositions = 1 -# number of grains = 1 -# x y z d c1 c2 c3 c4 -300e3 500e3 0 50e3 -1400e3 500e3 0 50e3 -1600e3 500e3 0 50e3 -1800e3 500e3 0 50e3 -500e3 500e3 0 300e3 -2000e3 500e3 0 400e3 -1500e3 500e3 0 50e3 -1400e3 500e3 0 200e3 -950e3 500e3 0 50e3 -1050e3 500e3 0 50e3 -950e3 500e3 0 210e3 -1050e3 500e3 0 210e3 diff --git a/tests/gwb-dat/random_uniform_texture_all_features_seed1.wb b/tests/gwb-dat/random_uniform_texture_all_features_seed1.wb deleted file mode 100644 index 76d0b5a69..000000000 --- a/tests/gwb-dat/random_uniform_texture_all_features_seed1.wb +++ /dev/null @@ -1,54 +0,0 @@ -{ - "version": "1.0", - "coordinate system":{"model":"cartesian"}, - "random number seed":1, - "features": - [ - { - "model":"continental plate", "name":"Overriding Plate", "max depth":100e3, - "coordinates":[[0,0],[0,1000e3],[1500e3,1000e3],[1500e3,0]], - "temperature models":[{"model":"uniform", "temperature":293}], - "composition models":[{"model":"uniform", "compositions":[0]}], - "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], - "normalize grain sizes":[true], "deflections":[0.1]}] - }, - - { - "model":"oceanic plate", "name":"Subducting Plate", "max depth":100e3, - "coordinates":[[1500e3,0],[1500e3,1000e3],[2500e3,1000e3],[2500e3,0]], - "temperature models":[{"model":"uniform", "temperature":293}], - "composition models":[{"model":"uniform", "compositions":[1]}], - "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], - "normalize grain sizes":[true], "deflections":[0.2]}] - }, - - { - "model":"mantle layer", "name":"Mantle", "min depth":100e3, - "coordinates":[[0,0],[0,1000e3],[2500e3,1000e3],[2500e3,0]], - "temperature models":[{"model":"uniform", "temperature":293}], - "composition models":[{"model":"uniform", "compositions":[2]}], - "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], - "normalize grain sizes":[true], "deflections":[0.3]}] - }, - - { - "model":"subducting plate", "name":"Slab", "min depth":0e3, - "coordinates":[[1500e3,0],[1500e3,500e3],[1500e3,1000e3]], "dip point":[10, 10], - "segments":[{"length":500e3, "thickness":[100e3], "angle":[50]}], - "temperature models":[{"model":"uniform", "temperature":10}], - "composition models":[{"model":"uniform", "compositions":[3]}], - "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], - "normalize grain sizes":[true], "deflections":[0.4]}] - }, - - { - "model":"fault", "name":"Fault", "min depth":0e3, - "coordinates":[[1000e3,0],[1000e3,500e3],[1000e3,1000e3]], "dip point":[10, 10], - "segments":[{"length":200e3, "thickness":[200e3], "angle":[90]}], - "temperature models":[{"model":"uniform", "temperature":10}], - "composition models":[{"model":"uniform", "compositions":[4]}], - "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], - "normalize grain sizes":[true], "deflections":[0.5]}] - } - ] -} \ No newline at end of file diff --git a/tests/gwb-dat/random_uniform_texture_all_features_seed1/screen-output.log b/tests/gwb-dat/random_uniform_texture_all_features_seed1/screen-output.log deleted file mode 100644 index e797cc6a9..000000000 --- a/tests/gwb-dat/random_uniform_texture_all_features_seed1/screen-output.log +++ /dev/null @@ -1,13 +0,0 @@ -# x y z d g T c0 c1 c2 c3 c4 gs0-0 gm0-0[0:0] gm0-0[0:1] gm0-0[0:2] gm0-0[1:0] gm0-0[1:1] gm0-0[1:2] gm0-0[2:0] gm0-0[2:1] gm0-0[2:2] tag -300e3 500e3 0 50e3 293 1 0 0 0 0 1 -0.800914 -0.591593 -0.0924877 0.566088 -0.798439 0.205035 -0.195143 0.111859 0.974375 0 -1400e3 500e3 0 50e3 293 1 0 0 0 0 1 -0.955406 -0.180561 0.233658 0.103564 -0.945899 -0.307488 0.276538 -0.269578 0.922418 0 -1600e3 500e3 0 50e3 293 0 1 0 0 0 1 -0.298222 -0.867278 -0.398613 0.862937 -0.423456 0.275724 -0.407925 -0.261751 0.874691 1 -1800e3 500e3 0 50e3 293 0 1 0 0 0 1 -0.684602 -0.394635 0.612848 0.550037 -0.831389 0.0790752 0.478309 0.391224 0.786234 1 -500e3 500e3 0 300e3 293 0 0 1 0 0 1 -0.422225 -0.823232 0.379494 0.267122 -0.513044 -0.81574 0.86624 -0.243055 0.436523 2 -2000e3 500e3 0 400e3 293 0 0 1 0 0 1 -0.191783 -0.93082 -0.311117 0.966116 -0.234847 0.107082 -0.172739 -0.280039 0.94432 2 -1500e3 500e3 0 50e3 10 0 0 0 1 0 1 -0.46567 -0.645243 0.60565 0.87519 -0.437189 0.207144 0.131125 0.62652 0.768296 3 -1400e3 500e3 0 200e3 10 0 0 0 1 0 1 -0.9776 -0.197078 0.0738833 0.155499 -0.912865 -0.377488 0.14184 -0.357544 0.923062 3 -950e3 500e3 0 50e3 10 0 0 0 0 1 1 0.407186 -0.895089 -0.181701 0.733709 0.202085 0.648716 -0.54394 -0.397464 0.739021 4 -1050e3 500e3 0 50e3 10 0 0 0 0 1 1 -0.728705 -0.338593 -0.595268 -0.202176 -0.724117 0.659379 -0.654305 0.600841 0.459212 4 -950e3 500e3 0 210e3 293 0 0 1 0 0 1 -0.815387 -0.468934 -0.339477 0.54922 -0.441202 -0.709717 0.183032 -0.765142 0.617298 2 -1050e3 500e3 0 210e3 293 0 0 1 0 0 1 -0.328436 -0.805474 -0.493297 0.917498 -0.148023 -0.36917 0.224338 -0.573848 0.787637 2 diff --git a/tests/gwb-dat/random_uniform_texture_all_features_seed1000.dat b/tests/gwb-dat/random_uniform_texture_all_features_seed1000.dat deleted file mode 100644 index d0c093e2e..000000000 --- a/tests/gwb-dat/random_uniform_texture_all_features_seed1000.dat +++ /dev/null @@ -1,20 +0,0 @@ -# This is a comment in the data -# file. -# Now define parameters: -# dim = 3 -# compositions = 5 -# grain compositions = 1 -# number of grains = 1 -# x y z d c1 c2 c3 c4 -300e3 500e3 0 50e3 -1400e3 500e3 0 50e3 -1600e3 500e3 0 50e3 -1800e3 500e3 0 50e3 -500e3 500e3 0 300e3 -2000e3 500e3 0 400e3 -1500e3 500e3 0 50e3 -1400e3 500e3 0 200e3 -950e3 500e3 0 50e3 -1050e3 500e3 0 50e3 -950e3 500e3 0 210e3 -1050e3 500e3 0 210e3 diff --git a/tests/gwb-dat/random_uniform_texture_all_features_seed1000.wb b/tests/gwb-dat/random_uniform_texture_all_features_seed1000.wb deleted file mode 100644 index 3a2cebd39..000000000 --- a/tests/gwb-dat/random_uniform_texture_all_features_seed1000.wb +++ /dev/null @@ -1,54 +0,0 @@ -{ - "version": "1.0", - "coordinate system":{"model":"cartesian"}, - "random number seed":1000, - "features": - [ - { - "model":"continental plate", "name":"Overriding Plate", "max depth":100e3, - "coordinates":[[0,0],[0,1000e3],[1500e3,1000e3],[1500e3,0]], - "temperature models":[{"model":"uniform", "temperature":293}], - "composition models":[{"model":"uniform", "compositions":[0]}], - "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], - "normalize grain sizes":[true], "deflections":[0.1]}] - }, - - { - "model":"oceanic plate", "name":"Subducting Plate", "max depth":100e3, - "coordinates":[[1500e3,0],[1500e3,1000e3],[2500e3,1000e3],[2500e3,0]], - "temperature models":[{"model":"uniform", "temperature":293}], - "composition models":[{"model":"uniform", "compositions":[1]}], - "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], - "normalize grain sizes":[true], "deflections":[0.2]}] - }, - - { - "model":"mantle layer", "name":"Mantle", "min depth":100e3, - "coordinates":[[0,0],[0,1000e3],[2500e3,1000e3],[2500e3,0]], - "temperature models":[{"model":"uniform", "temperature":293}], - "composition models":[{"model":"uniform", "compositions":[2]}], - "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], - "normalize grain sizes":[true], "deflections":[0.3]}] - }, - - { - "model":"subducting plate", "name":"Slab", "min depth":0e3, - "coordinates":[[1500e3,0],[1500e3,500e3],[1500e3,1000e3]], "dip point":[10, 10], - "segments":[{"length":500e3, "thickness":[100e3], "angle":[50]}], - "temperature models":[{"model":"uniform", "temperature":10}], - "composition models":[{"model":"uniform", "compositions":[3]}], - "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], - "normalize grain sizes":[true], "deflections":[0.4]}] - }, - - { - "model":"fault", "name":"Fault", "min depth":0e3, - "coordinates":[[1000e3,0],[1000e3,500e3],[1000e3,1000e3]], "dip point":[10, 10], - "segments":[{"length":200e3, "thickness":[200e3], "angle":[90]}], - "temperature models":[{"model":"uniform", "temperature":10}], - "composition models":[{"model":"uniform", "compositions":[4]}], - "grains models": [{"model":"random uniform distribution deflected", "compositions":[0], "grain sizes":[-1], - "normalize grain sizes":[true], "deflections":[0.5]}] - } - ] -} \ No newline at end of file diff --git a/tests/gwb-dat/random_uniform_texture_all_features_seed1000/screen-output.log b/tests/gwb-dat/random_uniform_texture_all_features_seed1000/screen-output.log deleted file mode 100644 index 3ecc83849..000000000 --- a/tests/gwb-dat/random_uniform_texture_all_features_seed1000/screen-output.log +++ /dev/null @@ -1,13 +0,0 @@ -# x y z d g T c0 c1 c2 c3 c4 gs0-0 gm0-0[0:0] gm0-0[0:1] gm0-0[0:2] gm0-0[1:0] gm0-0[1:1] gm0-0[1:2] gm0-0[2:0] gm0-0[2:1] gm0-0[2:2] tag -300e3 500e3 0 50e3 293 1 0 0 0 0 1 -0.990972 -0.127666 -0.0409292 0.129618 -0.990341 -0.0492332 -0.0342484 -0.0540939 0.997948 0 -1400e3 500e3 0 50e3 293 1 0 0 0 0 1 -0.863294 -0.503742 0.0311121 0.479282 -0.798937 0.363301 -0.158154 0.328547 0.931152 0 -1600e3 500e3 0 50e3 293 0 1 0 0 0 1 -0.55876 -0.762818 -0.325417 0.720801 -0.252629 -0.645465 0.410162 -0.595221 0.690999 1 -1800e3 500e3 0 50e3 293 0 1 0 0 0 1 -0.994372 -0.0277444 0.102248 0.0963064 -0.638921 0.76322 0.0441535 0.768772 0.637997 1 -500e3 500e3 0 300e3 293 0 0 1 0 0 1 -0.88545 -0.216658 -0.411142 0.0819602 -0.943615 0.32074 -0.45745 0.250302 0.853281 2 -2000e3 500e3 0 400e3 293 0 0 1 0 0 1 -0.505877 -0.721112 -0.473377 0.730473 -0.65 0.209545 -0.4588 -0.239785 0.855573 2 -1500e3 500e3 0 50e3 10 0 0 0 1 0 1 0.277217 -0.320014 -0.905948 0.59042 0.800607 -0.102137 0.757994 -0.506575 0.410885 3 -1400e3 500e3 0 200e3 10 0 0 0 1 0 1 0.784316 -0.497925 -0.370026 0.616522 0.559361 0.55409 -0.0689168 -0.662711 0.745697 3 -950e3 500e3 0 50e3 10 0 0 0 0 1 1 -0.459173 -0.634613 -0.621632 -0.305909 -0.544008 0.781329 -0.834014 0.548928 0.05566 4 -1050e3 500e3 0 50e3 10 0 0 0 0 1 1 0.615449 -0.728911 0.299852 0.271936 0.553452 0.787237 -0.739779 -0.402964 0.538839 4 -950e3 500e3 0 210e3 293 0 0 1 0 0 1 -0.577637 -0.64995 -0.493863 0.410363 -0.754211 0.51261 -0.705647 0.0934395 0.702375 2 -1050e3 500e3 0 210e3 293 0 0 1 0 0 1 0.248301 -0.967981 0.0368649 0.636139 0.134242 -0.759807 0.73053 0.212112 0.649103 2 diff --git a/tests/unit_tests/utilities.cc b/tests/unit_tests/utilities.cc new file mode 100644 index 000000000..a2ef6a924 --- /dev/null +++ b/tests/unit_tests/utilities.cc @@ -0,0 +1,96 @@ +/* + Copyright (C) 2021 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 . + */ + +#define DOCTEST_CONFIG_SUPER_FAST_ASSERTS + +#include "doctest/doctest.h" + +#include "world_builder/utilities.h" + +using namespace WorldBuilder; +using doctest::Approx; + + +TEST_CASE("multiply 3x3 matrices") +{ + std::array,3> mat1 = {{{{1,2,3}},{{4,5,6}},{{7,8,9}}}}; + std::array,3> mat2 = {{{{9,8,7}},{{6,5,4}},{{3,2,1}}}}; + std::array,3> mat3 = {{{{0.1,10,-7}},{{4,-0.5,9}},{{-9,0.3,-5}}}}; + std::array,3> result1 = {{{{30,24,18}},{{84,69,54}},{{138,114,90}}}}; //mat1*mat2 + std::array,3> result2 = {{{{90,114,138}},{{54,69,84}},{{18,24,30}}}}; //mat2*mat1 + std::array,3> result3 = {{{{-18.9,9.9,-4}},{{-33.6,39.3,-13}},{{-48.3,68.7,-22}}}}; //mat1*mat3 + std::array,3> result4 = {{{{-8.9,-5.8,-2.7}},{{65,77.5,90}},{{-42.8,-56.5,-70.2}}}}; //mat3*mat1 + std::array,3> result5 = {{{{-30.1,88.1,-26}},{{-15.4,58.7,-17}},{{-0.7,29.3,-8}}}}; //mat2*mat3 + std::array,3> result6 = {{{{39.9,36.8,33.7}},{{60,47.5,35}},{{-94.2,-80.5,-66.8}}}}; //mat3*mat2 + + std::array,3> result7 = Utilities::multiply_3x3_matrices(mat1, mat2); + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + CHECK(result1[i][j] == Approx(result7[i][j])); + } + } + + std::array,3> result8 = Utilities::multiply_3x3_matrices(mat2, mat1); + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + CHECK(result2[i][j] == Approx(result8[i][j])); + } + } + + std::array,3> result9 = Utilities::multiply_3x3_matrices(mat1, mat3); + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + CHECK(result3[i][j] == Approx(result9[i][j])); + } + } + + std::array,3> result10 = Utilities::multiply_3x3_matrices(mat3, mat1); + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + CHECK(result4[i][j] == Approx(result10[i][j])); + } + } + + std::array,3> result11 = Utilities::multiply_3x3_matrices(mat2, mat3); + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + CHECK(result5[i][j] == Approx(result11[i][j])); + } + } + + std::array,3> result12 = Utilities::multiply_3x3_matrices(mat3, mat2); + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + CHECK(result6[i][j] == Approx(result12[i][j])); + } + } + +} \ No newline at end of file