diff --git a/include/world_builder/features/subducting_plate_contours.h b/include/world_builder/features/subducting_plate_contours.h
new file mode 100644
index 000000000..6ea760d3a
--- /dev/null
+++ b/include/world_builder/features/subducting_plate_contours.h
@@ -0,0 +1,248 @@
+/*
+ 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_SUBDUCTING_PLATE_CONTOURS_H
+#define WORLD_BUILDER_FEATURES_SUBDUCTING_PLATE_CONTOURS_H
+
+
+#include "world_builder/features/subducting_plate_models/composition/interface.h"
+#include "world_builder/features/subducting_plate_models/grains/interface.h"
+#include "world_builder/features/subducting_plate_models/temperature/interface.h"
+#include "world_builder/features/subducting_plate_models/velocity/interface.h"
+#include "world_builder/objects/segment.h"
+#include "world_builder/bounding_box.h"
+#include "world_builder/objects/distance_from_surface.h"
+
+
+namespace WorldBuilder
+{
+ class Parameters;
+ class World;
+
+ namespace Features
+ {
+ using namespace FeatureUtilities;
+ namespace SubductingPlateModels
+ {
+ namespace Composition
+ {
+ class Interface;
+ } // namespace Composition
+ namespace Grains
+ {
+ class Interface;
+ } // namespace Grains
+ namespace Temperature
+ {
+ class Interface;
+ } // namespace Temperature
+ namespace Velocity
+ {
+ class Interface;
+ } // namespace Velocity
+ } // namespace SubductingPlateModels
+
+ /**
+ * This class represents a subducting plate and can implement submodules
+ * for temperature and composition. These submodules determine what
+ * the returned temperature or composition of the temperature and composition
+ * functions of this class will be.
+ */
+ class SubductingPlateContours final: public Interface
+ {
+ public:
+ /**
+ * constructor
+ */
+ SubductingPlateContours(WorldBuilder::World *world);
+
+ /**
+ * Destructor
+ */
+ ~SubductingPlateContours() 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 = "",
+ const std::vector &required_entries = {});
+
+ /**
+ * Produce a JSON snippet for the schema
+ */
+ static
+ void make_snippet(Parameters &prm);
+
+ /**
+ * declare and read in the world builder file into the parameters class
+ */
+ void parse_entries(Parameters &prm) override final;
+
+
+ /**
+ * Returns the bounding points for a BoundingBox object using two extreme points in all the surface
+ * coordinates and an additional buffer zone that accounts for the fault thickness and length. The first and second
+ * points correspond to the lower left and the upper right corners of the bounding box, respectively (see the
+ * documentation in include/bounding_box.h).
+ * For the spherical system, the buffer zone along the longitudal direction is calculated using the
+ * corresponding latitude points.
+ */
+ const BoundingBox<2> &get_surface_bounding_box () const;
+
+ /**
+ * Returns different values at a single point in one go stored in a vector of doubles.
+ *
+ * The properties input decides what each entry means, and the output is generated in the
+ * same order as the properties input. The properties input consists of
+ * a 3D array, where the first entry identifies the property and the last two entries
+ * provide extra information about that property.
+ *
+ * Temperature is identified by 1 and no extra information is needed. So temperature
+ * input usually looks like {1,0,0}. A temperature query prodoces one entry in the output
+ * vector.
+ *
+ * Composition is identified by 2. This produces one
+ * value in the output. The second entry identifies the composition number and the third
+ * number is not used. So a commposition query asking about composition 1 looks like this:
+ * {2,1,0}. A composition query prodoces one entry in the output vector.
+ *
+ * Grains are identified by 2. The second entry is the grain composition number and the third
+ * entry is the number of grains. A query about the grains, where it asks about composition 1
+ * (for example enstatite) and 500 grains, looks like this: {2,1,500}.
+ * A composition query prodoces n_grains*10 entries in the output vector. The first n_grains
+ * entries are the sizes of all the grains, and the other 9 entries are sets of rotation
+ * matrices. The rotation matrix entries are ordered [0][0],[0][1],[0][2],[1][0],[1][1],etc.
+ *
+ * The entries in output variable relates the index of the property to the index in the output.
+ */
+ void
+ properties(const Point<3> &position_in_cartesian_coordinates,
+ const Objects::NaturalCoordinate &position_in_natural_coordinates,
+ const double depth,
+ const std::vector> &properties,
+ const double gravity,
+ const std::vector &entry_in_output,
+ std::vector &output) const override final;
+
+ /**
+ * Returns a PlaneDistances object that has the distance from and along a subducting plate plane,
+ * calculated from the coordinates and the depth of the point.
+ */
+ // Objects::PlaneDistances
+ // distance_to_feature_plane(const Point<3> &position_in_cartesian_coordinates,
+ // const Objects::NaturalCoordinate &position_in_natural_coordinates,
+ // const double depth) const override;
+
+
+ private:
+ std::vector > default_temperature_models;
+ std::vector > default_composition_models;
+ std::vector > default_grains_models;
+ std::vector > default_velocity_models;
+
+ std::vector > default_segment_vector;
+
+ std::vector< std::vector > > sections_segment_vector;
+
+ // This vector stores segments to this coordinate/section.
+ //First used (raw) pointers to the segment relevant to this coordinate/section,
+ // but I do not trust it won't fail when memory is moved. So storing the all the data now.
+ std::vector > > segment_vector;
+
+ // todo: the memory of this can be greatly improved by
+ // or using a plugin system for the submodules, or
+ // putting the variables in a union. Although the memory
+ // used by this program is in real cases expected to be
+ // Insignificant compared to what a calling program may
+ // use, a smaller amount of memory used in here, could
+ // theoretically speed up the computation, because more
+ // relevant data could be stored in the cache. But this
+ // not a urgent problem, and would require testing.
+
+ /**
+ * This variable stores the depth at which the subducting
+ * plate starts. It makes this depth effectively the surface
+ * of the model for the slab.
+ */
+ double starting_depth;
+
+ /**
+ * The depth which below the subducting plate may no longer
+ * be present. This can not only help setting up models with
+ * less effort, but can also improve performance, because
+ * the algorithm doesn't have to search in locations below
+ * this depth.
+ */
+ double maximum_depth;
+
+ /**
+ * Stores the bounding points for a BoundingBox object using two extreme points in all the surface
+ * coordinates and an additional buffer zone that accounts for the fault thickness and length. The first and second
+ * points correspond to the lower left and the upper right corners of the bounding box, respectively (see the
+ * documentation in include/bounding_box.h).
+ * For the spherical system, the buffer zone along the longitudal direction is calculated using the
+ * corresponding latitude points.
+ */
+ BoundingBox<2> surface_bounding_box;
+
+ /**
+ * A point on the surface to which the subducting plates subduct.
+ */
+ Point<2> reference_point;
+
+ std::vector > slab_segment_lengths;
+ std::vector > slab_segment_thickness;
+ std::vector > slab_segment_top_truncation;
+ std::vector > > slab_segment_angles;
+ std::vector total_slab_length;
+ double maximum_total_slab_length;
+ double maximum_slab_thickness;
+ double slab_thickness = 100e3;
+ double top_truncation = 0;
+ double min_along_x;
+ double max_along_x;
+ double min_along_y;
+ double max_along_y;
+ double min_lat_cos_inv;
+ double max_lat_cos_inv;
+
+ double min_depth_local = 0;
+ double max_depth_local = 100e3;
+
+ std::vector>> regularized_points;
+ double buffer_around_slab_cartesian;
+ std::vector depths_of_contours = {0, 400e3};
+
+ unsigned int number_of_intervals = 5;
+ };
+ } // namespace Features
+} // namespace WorldBuilder
+
+#endif
diff --git a/include/world_builder/objects/bezier_curve.h b/include/world_builder/objects/bezier_curve.h
index 79401619b..c4b885c2b 100644
--- a/include/world_builder/objects/bezier_curve.h
+++ b/include/world_builder/objects/bezier_curve.h
@@ -66,24 +66,36 @@ namespace WorldBuilder
ClosestPointOnCurve closest_point_on_curve_segment(const Point<2> &p, const bool verbose = false) const;
/**
- * @brief
+ * @brief Returns a point that lies on the bezier curve at some interval x between coordinate i and coordinate
+ * i + 1. If x = 0, returns point i, if x = 1, returns point i + 1.
*
- * @param i
- * @param x
+ * @param i The index of the coordinate that defines the bezier curve.
+ * @param x The value used to determine additional points that lie on the bezier curve between coordinates i
+ * and i + 1
* @return Point<2>
*/
Point<2> operator()(const size_t i, const double x) const;
+ double get_bezier_arclength() const
+ {
+ return total_arclength;
+ }
+
+ std::vector get_bezier_segment_lengths() const
+ {
+ return lengths;
+ }
+
private:
std::vector > points;
- std::vector,2 > > control_points;
+ std::vector,2 > > control_points; // [ { {1,2}, {2,3} Point } array ] vector
std::vector lengths;
std::vector angles;
-
+ double total_arclength; // The arc length of the bezier curve from the first coordinate to the last coordinate
};
}
}
-#endif
+#endif
\ No newline at end of file
diff --git a/include/world_builder/objects/contour.h b/include/world_builder/objects/contour.h
new file mode 100644
index 000000000..4f915378a
--- /dev/null
+++ b/include/world_builder/objects/contour.h
@@ -0,0 +1,78 @@
+/*
+ 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_OBJECTS_CONTOUR_H
+#define WORLD_BUILDER_OBJECTS_CONTOUR_H
+
+
+#include "world_builder/types/plugin_system.h"
+
+
+namespace WorldBuilder
+{
+ class Parameters;
+
+ namespace Objects
+ {
+
+ /**
+ * This class represents an actual segment
+ */
+ template
+ class Contour
+ {
+ public:
+
+ /**
+ * A constructor for the clone and set_entry function
+ */
+ Contour(const double default_depth,
+ const WorldBuilder::Point<2> &default_thickness,
+ const WorldBuilder::Point<2> &default_top_truncation,
+ std::vector > temperature_systems,
+ std::vector > composition_systems,
+ std::vector > grains_systems,
+ std::vector > velocity_systems);
+
+ /**
+ * Copy constructor
+ */
+ Contour(Contour const &other);
+
+ /**
+ * Destructor
+ */
+ ~Contour();
+
+ double value_depth;
+ WorldBuilder::Point<2> value_thickness;
+ WorldBuilder::Point<2> value_top_truncation;
+ std::vector > temperature_systems;
+ std::vector > composition_systems;
+ std::vector > grains_systems;
+ std::vector > velocity_systems;
+
+ protected:
+ private:
+
+ };
+ } // namespace Objects
+} // namespace WorldBuilder
+
+#endif
diff --git a/include/world_builder/types/contour.h b/include/world_builder/types/contour.h
new file mode 100644
index 000000000..e434f4d33
--- /dev/null
+++ b/include/world_builder/types/contour.h
@@ -0,0 +1,102 @@
+/*
+ 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_TYPES_CONTOUR_H
+#define WORLD_BUILDER_TYPES_CONTOUR_H
+
+
+#include "world_builder/types/plugin_system.h"
+
+
+namespace WorldBuilder
+{
+ class Parameters;
+
+ namespace Types
+ {
+
+ /**
+ * This class represents a segment value with documentation
+ */
+ class Contour : public Interface
+ {
+ public:
+ /**
+ * A constructor
+ */
+ Contour(const double default_length,
+ const WorldBuilder::Point<2> &default_thickness,
+ const WorldBuilder::Point<2> &default_top_truncation,
+ const WorldBuilder::Point<2> &default_angle,
+ const Types::Interface &temperature_plugin_system,
+ const Types::Interface &composition_plugin_system,
+ const Types::Interface &grains_plugin_system_,
+ const Types::Interface &gvelocity_plugin_system_);
+
+ /**
+ * A constructor for the load_entry function
+ */
+ Contour(double default_length,
+ WorldBuilder::Point<2> default_thickness,
+ WorldBuilder::Point<2> default_angle,
+ std::string description);
+
+ /**
+ * Copy constructor
+ */
+ Contour(Contour const &other);
+
+
+ /**
+ * Destructor
+ */
+ ~Contour() override;
+
+ /**
+ * Todo
+ */
+ void write_schema(Parameters &prm,
+ const std::string &name,
+ const std::string &documentation) const override final;
+
+
+ double value_length;
+ double default_length;
+ WorldBuilder::Point<2> value_thickness;
+ WorldBuilder::Point<2> default_thickness;
+ WorldBuilder::Point<2> default_top_truncation;
+ WorldBuilder::Point<2> value_angle;
+ WorldBuilder::Point<2> default_angle;
+ std::unique_ptr temperature_plugin_system;
+ std::unique_ptr composition_plugin_system;
+ std::unique_ptr grains_plugin_system;
+ std::unique_ptr velocity_plugin_system;
+
+ protected:
+ Contour *clone_impl() const override final
+ {
+ return new Contour(*this);
+ };
+ private:
+
+ };
+ } // namespace Types
+} // namespace WorldBuilder
+
+#endif
diff --git a/include/world_builder/types/interface.h b/include/world_builder/types/interface.h
index 40cd8b586..8251cdc43 100644
--- a/include/world_builder/types/interface.h
+++ b/include/world_builder/types/interface.h
@@ -39,7 +39,7 @@ namespace WorldBuilder
enum class type
{
- None,Bool,String,Double,Int,UnsignedInt,Array,Object,List,Point2D,Point3D,CoordinateSystem,PluginSystem,Segment,ConstantLayer,ValueAtPoints,OneOf
+ None,Bool,String,Double,Int,UnsignedInt,Array,Object,List,Point2D,Point3D,CoordinateSystem,PluginSystem,Segment,ConstantLayer,ValueAtPoints,Contour,OneOf
};
class Interface
diff --git a/include/world_builder/utilities.h b/include/world_builder/utilities.h
index 68d389714..2615bb80f 100644
--- a/include/world_builder/utilities.h
+++ b/include/world_builder/utilities.h
@@ -487,6 +487,19 @@ namespace WorldBuilder
std::vector
calculate_effective_trench_and_plate_ages(std::vector ridge_parameters, double distance_along_plane);
+ /*
+ *
+ * @brief Calculate the distance of a point to a the closest point on a rectilinear plane
+ * @param p1 The first point defining the plane
+ * @param p2 The second point defining the plane
+ * @param p3 The third point defining the plane
+ * @param check_point The point to calculate the distance to
+ */
+ double
+ calculate_distance_from_point_to_rectilinear_plane(const Point<3> p1,
+ const Point<3> p2,
+ const Point<3> p3,
+ const Point<3> check_point);
/*
* Returns the result of the multiplication of two 3*3 matrix,
* used in applying the random uniform distribution rotation matrix
diff --git a/source/world_builder/features/interface.cc b/source/world_builder/features/interface.cc
index 0678ad063..2b771f503 100644
--- a/source/world_builder/features/interface.cc
+++ b/source/world_builder/features/interface.cc
@@ -114,7 +114,8 @@ namespace WorldBuilder
"If the tag is not provided or empty, it is set to the model name.");
prm.declare_entry("coordinates", Types::Array(Types::Point<2>(), 1),
"An array of 2d Points representing an array of coordinates where the feature is located.");
-
+ prm.declare_entry("geometry type", Types::String("contours"),
+ "The model geometry type");
prm.declare_entry("interpolation",Types::String("global"),
"What type of interpolation should be used to enforce the minimum points per "
"distance parameter. Options are 'global' and "
diff --git a/source/world_builder/features/subducting_plate_contours.cc b/source/world_builder/features/subducting_plate_contours.cc
new file mode 100644
index 000000000..30d3d5d10
--- /dev/null
+++ b/source/world_builder/features/subducting_plate_contours.cc
@@ -0,0 +1,356 @@
+/*
+ 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/subducting_plate_contours.h"
+
+
+#include "glm/glm.h"
+#include "world_builder/types/array.h"
+#include "world_builder/types/double.h"
+#include "world_builder/types/object.h"
+#include "world_builder/types/point.h"
+#include "world_builder/types/segment.h"
+#include "world_builder/types/contour.h"
+#include "world_builder/types/unsigned_int.h"
+#include "world_builder/world.h"
+#include
+#include
+
+/*
+The general form of the contours geometry type will probably look like this.
+
+{
+ "model": "subducting plate",
+ "name": "My SLAB",
+ "geometry type": "contours",
+
+ "contours": [[0, 0],
+ [1,2]
+ {"coordinate":[1, 1]},
+ {"coordinate":[0, 0], "thickness": 10e3},
+ [3,5]
+ .
+ .
+ }], "thickness": 10e3
+
+ OR
+
+ "contours": [{"coordinate":[0, 0], "thickness": 10e3},
+ {"coordinate":[1, 1], "thickness": 20e3},
+ .
+ .
+ .
+ }]
+
+ "depth contours": [{"depth":0.0, "contours":[{"coord":[0, 0],"thickness":10e3}, [1, 0], [1, 1], [0, 1]], "thickness":1000, "temperature models":[]},
+ {"depth":100.0, "contours":[[0, 0], [1, 0], [1, 1], [0, 1]], "thickness":1000}],
+ "temperature models":[{"model":"uniform", "temperature":600.0}],
+ "composition models":[{"model":"uniform", "compositions": [0], "fractions":[1.0]}]
+}
+
+This should result in a slab that has a temperature of 600.0 and a composition of 0.0 at all depths.
+This can then be extended to have variable temperature and composition models at different depths and
+different along strike points.
+
+The workflow will be to first fit the contour with a bezier curve, determine the arc-length of the curve,
+and then divide the arc-length into even intervals. All contours will be divided into the same number of
+intervals. Then the corresponding indices of the contours will be used to connect the contour and create
+a 3D surface.
+*/
+
+using namespace std;
+
+namespace WorldBuilder
+{
+ using namespace Utilities;
+
+ namespace Features
+ {
+ SubductingPlateContours::SubductingPlateContours(WorldBuilder::World *world_)
+ :
+ reference_point(0,0,cartesian)
+ {
+ this->world = world_;
+ this->name = "subducting plate contours";
+ }
+
+ SubductingPlateContours::~SubductingPlateContours()
+ = default;
+
+
+
+ void SubductingPlateContours::make_snippet(Parameters &prm)
+ {
+ using namespace rapidjson;
+ Document &declarations = prm.declarations;
+
+ const std::string path = prm.get_full_json_path();
+
+ /*
+ ideally:
+ {
+ "model": "fault",
+ "name": "${1:default name}",
+ "dip point":[0.0,0.0],
+ "coordinates": [[0.0,0.0]],
+ "segments": [],
+ "sections": [],
+ "temperature models":[{"model":"uniform", "temperature":600.0}],
+ "composition models":[{"model":"uniform", "compositions": [0], "fractions":[1.0]}]
+ }
+ */
+
+ Pointer((path + "/body").c_str()).Set(declarations,"object");
+ Pointer((path + "/body/model").c_str()).Set(declarations,"subducting plate");
+ Pointer((path + "/body/name").c_str()).Set(declarations,"${1:My Subducting Plate}");
+ Pointer((path + "/body/coordinates").c_str()).Create(declarations).SetArray();
+ Pointer((path + "/body/temperature models").c_str()).Create(declarations).SetArray();
+ Pointer((path + "/body/composition models").c_str()).Create(declarations).SetArray();
+ }
+
+
+
+ void
+ SubductingPlateContours::declare_entries(Parameters &prm,
+ const std::string &parent_name,
+ const std::vector &required_entries)
+ {
+ prm.declare_entry("geometry type", Types::String("contours"),
+ "The model geometry type");
+ prm.declare_entry("temperature models",
+ Types::PluginSystem("", Features::SubductingPlateModels::Temperature::Interface::declare_entries, {"model"}),
+ "A list of temperature models.");
+ prm.declare_entry("composition models",
+ Types::PluginSystem("", Features::SubductingPlateModels::Composition::Interface::declare_entries, {"model"}),
+ "A list of composition models.");
+ prm.declare_entry("grains models",
+ Types::PluginSystem("", Features::SubductingPlateModels::Grains::Interface::declare_entries, {"model"}),
+ "A list of grains models.");
+ prm.declare_entry("velocity models",
+ Types::PluginSystem("", Features::SubductingPlateModels::Velocity::Interface::declare_entries, {"model"}),
+ "A list of velocity models.");
+ }
+
+ void
+ SubductingPlateContours::parse_entries(Parameters &prm)
+ {
+ const CoordinateSystem coordinate_system = prm.coordinate_system->natural_coordinate_system();
+
+ this->name = prm.get("name");
+ std::string geometry_type = prm.get("geometry type");
+
+ std::string tag = prm.get("tag");
+ if (tag == "")
+ {
+ tag = "subducting plate";
+ }
+ this->tag_index = FeatureUtilities::add_vector_unique(this->world->feature_tags,tag);
+
+ regularized_points.resize(depths_of_contours.size());
+ slab_segment_lengths.resize(depths_of_contours.size());
+ slab_segment_thickness.resize(depths_of_contours.size());
+ slab_segment_top_truncation.resize(depths_of_contours.size());
+ slab_segment_angles.resize(depths_of_contours.size());
+
+ // In reality, the variable below should be calculated like this:
+ // std::vector bezier_sub_lengths = this->bezier_curve.get_bezier_segment_lengths();
+ // This is a VECTOR equal to N-1, where N is the number of coordinates which defines the feature.
+ // This is because for N points, you need N-1 bezier curves, and each bezier curve can be a different
+ // length.
+ const double bezier_sublengths = 2500e3;
+
+ // In reality, the variable below should be calculated like this:
+ // const double total_arclength = this->bezier_curve.get_bezier_arclength();
+ const double total_arclength = 5000e3;
+
+
+ std::vector> coordinates_1 = {Point<2>(1500e3, 0.0, cartesian), Point<2>(1500e3, 2500e3, cartesian), Point<2>(1500e3, 5000e3, cartesian)};
+ // std::vector> coordinates_2 = {Point<2>(500e3, 0.0, cartesian), Point<2>(500e3, 2500e3, cartesian), Point<2>(500e3, 5000e3, cartesian)};
+ std::vector> coordinates_2 = {Point<2>(1200e3, 0.0, cartesian), Point<2>(1200e3, 2500e3, cartesian), Point<2>(1200e3, 5000e3, cartesian)};
+ std::vector> coordinates_3 = {Point<2>(1200e3, 0.0, cartesian), Point<2>(1200e3, 2500e3, cartesian), Point<2>(1200e3, 5000e3, cartesian)};
+ // std::vector> coordinates_3 = {Point<2>(1000e3, 0.0, cartesian), Point<2>(1000e3, 2500e3, cartesian), Point<2>(1000e3, 5000e3, cartesian)};
+ Objects::BezierCurve bezier_curve;
+ for (unsigned int i = 0; i < depths_of_contours.size(); ++i)
+ {
+ // Somehow here I need to specify which coordinates I am using. This should be possible by using the
+ // index i which corresponds to the index of the depth value within the array of depth points.
+ if (i == 0)
+ {
+ // The first bezier curve
+ bezier_curve = Objects::BezierCurve(coordinates_1);
+ }
+
+ if (i == 1)
+ {
+ // The second bezier curve
+ bezier_curve = Objects::BezierCurve(coordinates_2);
+ }
+
+ if (i == 2)
+ {
+ // The second bezier curve
+ bezier_curve = Objects::BezierCurve(coordinates_3);
+ }
+ // The sum of all the previous bezier curves WITHIN the current contour
+ double distance_along_contour = 0;
+ double distance_along_current_bezier_curve = 0;
+
+ // The index of the bezier curve. Currently this doesnt do too much, but this will be used
+ // to index the length of the bezier curve as well.
+ unsigned int which_bezier_curve = 0;
+
+ for (unsigned int current_interval = 0; current_interval < number_of_intervals; ++current_interval)
+ {
+ // The increment in distance along the current bezier curve
+ const double distance_increment = total_arclength / static_cast(number_of_intervals);
+ distance_along_current_bezier_curve += distance_increment;
+
+ if (distance_along_current_bezier_curve <= bezier_sublengths)
+ {
+ const double t_value = distance_along_current_bezier_curve / bezier_sublengths;
+ regularized_points[i].emplace_back(bezier_curve(which_bezier_curve, t_value));
+ distance_along_contour += distance_increment;
+ }
+
+ else
+ {
+ ++which_bezier_curve;
+ const double remaining_distance = distance_along_current_bezier_curve -
+ bezier_sublengths;
+
+ const double t_value = remaining_distance / bezier_sublengths;
+ regularized_points[i].emplace_back(bezier_curve(which_bezier_curve, t_value));
+ distance_along_contour += distance_increment;
+ distance_along_current_bezier_curve = remaining_distance;
+ }
+ }
+ }
+
+ for (unsigned int i = 0; i < depths_of_contours.size(); ++i)
+ {
+ slab_segment_lengths[i].resize(number_of_intervals);
+ slab_segment_thickness[i].resize(number_of_intervals);
+ slab_segment_top_truncation[i].resize(number_of_intervals);
+ slab_segment_angles[i].resize(number_of_intervals, Point<2>(invalid));
+ for (unsigned int t = 0; t < number_of_intervals; ++t)
+ {
+ std::array P1 = {regularized_points[i][t][0], regularized_points[i][t][1], depths_of_contours[i]};
+ std::array P2;
+ if (i == depths_of_contours.size() - 1)
+ {
+ P2 = {regularized_points[i][t][0], regularized_points[i][t][1], depths_of_contours[i]};
+ }
+ else
+ {
+ P2 = {regularized_points[i + 1][t][0], regularized_points[i + 1][t][1], depths_of_contours[i + 1]};
+ }
+ // std::array P2 = {regularized_points[i + 1][t][0], regularized_points[i + 1][t][1], depths_of_contours[i + 1]};
+ std::array P1P2 = {P2[0] - P1[0], P2[1] - P1[1], P2[2] - P1[2]};
+
+ slab_segment_lengths[i][t] = std::sqrt(P1P2[0] * P1P2[0] + P1P2[1] * P1P2[1] + P1P2[2] * P1P2[2]);
+ slab_segment_thickness[i][t] = slab_thickness;
+
+ slab_segment_angles[i][t][0] = std::sin(P1P2[2] / slab_segment_lengths[i][t]);
+
+ if (i == depths_of_contours.size() - 1)
+ // Angle with the Z-axis
+ slab_segment_angles[i][t][0] = slab_segment_angles[0][t][0];
+
+ slab_segment_top_truncation[i][t] = top_truncation;
+ }
+ }
+
+
+ if (coordinate_system == spherical)
+ {
+ // When spherical, input is in degrees, so change to radians for internal use.
+ reference_point *= (Consts::PI/180.);
+ }
+
+ default_temperature_models.resize(0);
+ default_composition_models.resize(0);
+ default_grains_models.resize(0);
+ default_velocity_models.resize(0);
+ prm.get_shared_pointers("temperature models", default_temperature_models);
+ prm.get_shared_pointers("composition models", default_composition_models);
+ prm.get_shared_pointers("grains models", default_grains_models);
+ prm.get_shared_pointers("velocity models", default_velocity_models);
+ }
+
+
+
+ const BoundingBox<2> &
+ SubductingPlateContours::get_surface_bounding_box () const
+ {
+ return surface_bounding_box;
+ }
+
+
+ void
+ SubductingPlateContours::properties(const Point<3> &position_in_cartesian_coordinates,
+ const Objects::NaturalCoordinate &position_in_natural_coordinates,
+ const double depth,
+ const std::vector> &properties,
+ const double gravity_norm,
+ const std::vector &entry_in_output,
+ std::vector &output) const
+ {
+ // std::vector depths_of_contours = {0, 100e3};
+ std::vector absolute_distances_from_planes_to_check_point;
+ std::vector signed_distances_from_planes_to_check_point;
+ for (unsigned int i = 0; i < depths_of_contours.size() - 1; ++i)
+ {
+ std::vector> contour_i = regularized_points[i];
+ std::vector> contour_i_plus_1 = regularized_points[i + 1];
+ for (unsigned int j = 0; j < number_of_intervals - 1; ++j)
+ {
+ const Point<3> p1 = Point<3>(contour_i[j][0], contour_i[j][1], 500e3 - depths_of_contours[i], cartesian);
+ const Point<3> p2 = Point<3>(contour_i[j + 1][0], contour_i[j + 1][1], 500e3 - depths_of_contours[i], cartesian);
+ const Point<3> p3 = Point<3>(contour_i_plus_1[j][0], contour_i_plus_1[j][1], 500e3 - depths_of_contours[i + 1], cartesian);
+ const Point<3> p4 = Point<3>(contour_i_plus_1[j + 1][0], contour_i_plus_1[j + 1][1], 500e3 - depths_of_contours[i + 1], cartesian);
+ const double distance_1 = WorldBuilder::Utilities::calculate_distance_from_point_to_rectilinear_plane(p1, p2, p3, position_in_cartesian_coordinates);
+ const double distance_2 = WorldBuilder::Utilities::calculate_distance_from_point_to_rectilinear_plane(p4, p3, p2, position_in_cartesian_coordinates);
+ absolute_distances_from_planes_to_check_point.emplace_back(std::abs(distance_1));
+ absolute_distances_from_planes_to_check_point.emplace_back(std::abs(distance_2));
+ signed_distances_from_planes_to_check_point.emplace_back(distance_1);
+ signed_distances_from_planes_to_check_point.emplace_back(distance_2);
+ }
+ }
+
+ // Find the iterator to the minimum element
+ auto minIt = std::min_element(absolute_distances_from_planes_to_check_point.begin(), absolute_distances_from_planes_to_check_point.end());
+
+ // Calculate the index
+ const unsigned int minIndex = std::distance(absolute_distances_from_planes_to_check_point.begin(), minIt);
+ const double min_distance_to_plane = signed_distances_from_planes_to_check_point[minIndex];
+
+ if (min_distance_to_plane >= -slab_segment_thickness[0][0] &&
+ min_distance_to_plane <= 0 &&
+ depth <= depths_of_contours[depths_of_contours.size() - 1])
+ {
+ output[entry_in_output[0]] = 10.0;
+ }
+ }
+ /**
+ * Register plugin
+ */
+ WB_REGISTER_FEATURE(SubductingPlateContours, subducting plate contours)
+ } // namespace Features
+} // namespace WorldBuilder
+
diff --git a/source/world_builder/objects/bezier_curve.cc b/source/world_builder/objects/bezier_curve.cc
index 38c974970..1021b839a 100644
--- a/source/world_builder/objects/bezier_curve.cc
+++ b/source/world_builder/objects/bezier_curve.cc
@@ -46,6 +46,9 @@ namespace WorldBuilder
std::vector angle_constraints = angle_constraints_input;
angle_constraints.resize(n_points,NaN::DQNAN);
+ const unsigned int max_arclength_discretization = 10;
+
+
// if no angle is provided, compute the angle as the average angle between the previous and next point.
// The first angle points at the second point and the last angle points at the second to last point.
// The check points are set at a distance of 1/10th the line length from the point in the direction of the angle.
@@ -123,11 +126,28 @@ namespace WorldBuilder
const bool side_of_line_2 = (p1[0] - p2[0]) * (p3[1] - p1[1])
- (p1[1] - p2[1]) * (p3[0] - p1[0])
< 0;
- if (side_of_line_1 == side_of_line_2)
+
+ // Check to see if the angles are different. If the angles are the same, points p1, and p2
+ // are colinear, and therefore the control points will also be colinear with p1 and p2. This
+ // makes determining which 'side' the control points lie meaningless.
+ if ( (std::abs(std::abs(angles[0]) - std::abs(angles[1]))) > 1e-10)
+ {
+ if (side_of_line_1 == side_of_line_2)
+ {
+ // use a 180 degree rotated angle to create this control_point
+ control_points[0][1][0] = cos(angles[1]+Consts::PI)*length*fraction_of_length+p2[0];
+ control_points[0][1][1] = sin(angles[1]+Consts::PI)*length*fraction_of_length+p2[1];
+ }
+ }
+
+ // There is no closed-form analytic way to express the arc-length of a cubic bezier curve. We approximate
+ // the arc-length by dividing the curve into 20 points piecewise linearly connect them. We also store the
+ // length of the bezier curve within each of these intervals. We calculate the points that lie on the bezier
+ // curve using the operator function below.
+ for (unsigned int t_value = 1; t_value <= max_arclength_discretization; ++t_value)
{
- // use a 180 degree rotated angle to create this control_point
- control_points[0][1][0] = cos(angles[1]+Consts::PI)*length*fraction_of_length+p2[0];
- control_points[0][1][1] = sin(angles[1]+Consts::PI)*length*fraction_of_length+p2[1];
+ lengths[0] = (operator()(0, static_cast (t_value)/max_arclength_discretization) - operator()(0, static_cast (t_value - 1)/max_arclength_discretization)).norm();
+ total_arclength += lengths[0];
}
}
}
@@ -139,7 +159,6 @@ namespace WorldBuilder
const double length = (points[p_i]-points[p_i+1]).norm(); // can be squared
control_points[p_i][0][0] = cos(angles[p_i])*length*fraction_of_length+p1[0];
control_points[p_i][0][1] = sin(angles[p_i])*length*fraction_of_length+p1[1];
-
{
// Determine which side of the line the control points lie on
const bool side_of_line_1 = (p1[0] - p2[0]) * (control_points[p_i-1][1][1] - p1[1])
@@ -148,6 +167,10 @@ namespace WorldBuilder
const bool side_of_line_2 = (p1[0] - p2[0]) * (control_points[p_i][0][1] - p1[1])
- (p1[1] - p2[1]) * (control_points[p_i][0][0] - p1[0])
< 0;
+
+ // Check to see if the angles are different. If the angles are the same, points p1, p2, and p3
+ // are colinear, and therefore the control points will also be colinear with p1, p2 and p3. This
+ // makes determining which 'side' the control points lie meaningless.
if (side_of_line_1 == side_of_line_2)
{
// use a 180 degree rotated angle to create this control_point
@@ -156,8 +179,8 @@ namespace WorldBuilder
}
}
- control_points[p_i][1][0] = cos(angles[p_i+1])*length*fraction_of_length+points[p_i+1][0];
- control_points[p_i][1][1] = sin(angles[p_i+1])*length*fraction_of_length+points[p_i+1][1];
+ control_points[p_i][1][0] = cos(angles[p_i+1])*length*fraction_of_length+p2[0];
+ control_points[p_i][1][1] = sin(angles[p_i+1])*length*fraction_of_length+p2[1];
if (p_i+1 < n_points-1)
{
@@ -168,15 +191,37 @@ namespace WorldBuilder
const bool side_of_line_2 = (p1[0] - p2[0]) * (p3[1] - p1[1])
- (p1[1] - p2[1]) * (p3[0] - p1[0])
< 0;
- if (side_of_line_1 == side_of_line_2)
- {
- // use a 180 degree rotated angle to create this control_point
- control_points[p_i][1][0] = cos(angles[p_i+1]+Consts::PI)*length*fraction_of_length+p2[0];
- control_points[p_i][1][1] = sin(angles[p_i+1]+Consts::PI)*length*fraction_of_length+p2[1];
- }
+
+ // Check to see if the angles are different. If the angles are the same, points p1, p2, and p3
+ // are colinear, and therefore the control points will also be colinear with p1, p2 and p3. This
+ // makes determining which 'side' the control points lie meaningless.
+ if (std::abs(std::abs(angles[p_i]) - std::abs(angles[p_i + 1])) > 1e-10)
+ if (side_of_line_1 == side_of_line_2)
+ {
+ // use a 180 degree rotated angle to create this control_point
+ control_points[p_i][1][0] = cos(angles[p_i+1]+Consts::PI)*length*fraction_of_length+p2[0];
+ control_points[p_i][1][1] = sin(angles[p_i+1]+Consts::PI)*length*fraction_of_length+p2[1];
+ }
+ }
+
+ // There is no closed-form analytic way to express the arc-length of a cubic bezier curve. We approximate
+ // the arc-length by dividing the curve into 20 points piecewise linearly connect them. We also store the
+ // length of the bezier curve within each of these intervals. We calculate the points that lie on the bezier
+ // curve using the operator function below.
+ for (unsigned int t_value = 1; t_value <= max_arclength_discretization; ++t_value)
+ {
+ lengths[p_i] = (operator()(p_i, static_cast (t_value)/max_arclength_discretization) -
+ operator()(p_i, static_cast ((t_value - 1))/max_arclength_discretization)).norm();
+ total_arclength += lengths[p_i];
}
}
}
+
+ else
+ {
+ lengths[0] = (points[0]-points[1]).norm();
+ total_arclength = lengths[0];
+ }
}
@@ -188,7 +233,6 @@ namespace WorldBuilder
return (1-t)*(1-t)*(1-t)*points[i] + 3*(1-t)*(1-t)*t*control_points[i][0] + 3.*(1-t)*t*t*control_points[i][1]+t*t*t*points[i+1];
}
-
ClosestPointOnCurve
BezierCurve::closest_point_on_curve_segment(const Point<2> &check_point,
const bool verbose) const
@@ -205,7 +249,7 @@ namespace WorldBuilder
#endif
const Point<2> &p1 = points[cp_i];
const Point<2> &p2 = points[cp_i+1];
- //min_squared_distance = std::min(std::min(min_squared_distance,(check_point-p1).norm_square()),(check_point-p1).norm_square());
+ // min_squared_distance = std::min(std::min(min_squared_distance,(check_point-p1).norm_square()),(check_point-p1).norm_square());
// Getting an estimate for where the closest point is with a linear approximation
const Point<2> P1P2 = p2-p1;
@@ -560,5 +604,6 @@ namespace WorldBuilder
}
return closest_point_on_curve;
}
+
} // namespace Objects
-} // namespace WorldBuilder
+} // namespace WorldBuilder
\ No newline at end of file
diff --git a/source/world_builder/objects/contour.cc b/source/world_builder/objects/contour.cc
new file mode 100644
index 000000000..cce120ac7
--- /dev/null
+++ b/source/world_builder/objects/contour.cc
@@ -0,0 +1,132 @@
+/*
+ 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
+
+#include "world_builder/objects/contour.h"
+
+
+namespace WorldBuilder
+{
+
+ namespace Features
+ {
+ namespace FaultModels
+ {
+ namespace Composition
+ {
+ class Interface;
+ } // namespace Composition
+ namespace Grains
+ {
+ class Interface;
+ } // namespace Grains
+ namespace Temperature
+ {
+ class Interface;
+ } // namespace Temperature
+ namespace Velocity
+ {
+ class Interface;
+ } // namespace Velocity
+ } // namespace FaultModels
+ namespace SubductingPlateModels
+ {
+ namespace Composition
+ {
+ class Interface;
+ } // namespace Composition
+ namespace Grains
+ {
+ class Interface;
+ } // namespace Grains
+ namespace Temperature
+ {
+ class Interface;
+ } // namespace Temperature
+ namespace Velocity
+ {
+ class Interface;
+ } // namespace Velocity
+ } // namespace SubductingPlateModels
+ } // namespace Features
+
+
+ namespace Objects
+ {
+ // todo update function
+ template
+ Contour::Contour(const double default_depth_,
+ const WorldBuilder::Point<2> &default_thickness_,
+ const WorldBuilder::Point<2> &default_top_truncation_,
+ const std::vector > temperature_systems_,
+ const std::vector > composition_systems_,
+ const std::vector > grains_systems_,
+ const std::vector > velocity_systems_)
+ :
+ value_depth(default_depth_),
+ value_thickness(default_thickness_),
+ value_top_truncation(default_top_truncation_),
+ temperature_systems(std::move(temperature_systems_)),
+ composition_systems(std::move(composition_systems_)),
+ grains_systems(std::move(grains_systems_)),
+ velocity_systems(std::move(velocity_systems_))
+ {
+ }
+
+ template
+ Contour::Contour(Contour const &other)
+ :
+ value_depth(other.value_depth),
+ value_thickness(other.value_thickness),
+ value_top_truncation(other.value_top_truncation),
+ temperature_systems(other.temperature_systems),
+ composition_systems(other.composition_systems),
+ grains_systems(other.grains_systems),
+ velocity_systems(other.velocity_systems)
+ {
+ }
+
+ template
+ Contour::~Contour ()
+ = default;
+
+
+ /**
+ * Todo: Returns a vector of pointers to the Point<3> Type based on the provided name.
+ * Note that the variable with this name has to be loaded before this function is called.
+ */
+ template class
+ Contour;
+
+ /**
+ * Todo: Returns a vector of pointers to the Point<3> Type based on the provided name.
+ * Note that the variable with this name has to be loaded before this function is called.
+ */
+ template class
+ Contour;
+
+ /**
+ * Todo: Returns a vector of pointers to the Point<3> Type based on the provided name.
+ * Note that the variable with this name has to be loaded before this function is called.
+ */
+ //template class
+ //Contour;
+
+ } // namespace Objects
+} // namespace WorldBuilder
\ No newline at end of file
diff --git a/source/world_builder/parameters.cc b/source/world_builder/parameters.cc
index 5e1b9d382..0f4f6144d 100644
--- a/source/world_builder/parameters.cc
+++ b/source/world_builder/parameters.cc
@@ -36,6 +36,7 @@
#include "world_builder/features/plume_models/temperature/interface.h"
#include "world_builder/features/plume_models/velocity/interface.h"
#include "world_builder/features/subducting_plate.h"
+#include "world_builder/features/subducting_plate_contours.h"
#include "world_builder/features/subducting_plate_models/velocity/interface.h"
#include "world_builder/gravity_model/interface.h"
#include "world_builder/types/object.h"
@@ -1815,6 +1816,29 @@ namespace WorldBuilder
return true;
}
+ template<>
+ bool
+ Parameters::get_unique_pointers(const std::string &name, std::vector > &vector)
+ {
+ vector.resize(0);
+ const std::string strict_base = this->get_full_json_path();
+ if (Pointer((strict_base + "/" + name).c_str()).Get(parameters) != nullptr)
+ {
+ Value *array = Pointer((strict_base + "/" + name).c_str()).Get(parameters);
+
+ for (size_t i = 0; i < array->Size(); ++i )
+ {
+ vector.push_back(std::make_unique(&world));
+ }
+ }
+ else
+ {
+ return false;
+ }
+
+ return true;
+ }
+
template<>
bool
Parameters::get_unique_pointers(const std::string &name, std::vector > &vector)
diff --git a/source/world_builder/types/contour.cc b/source/world_builder/types/contour.cc
new file mode 100644
index 000000000..f1767cea6
--- /dev/null
+++ b/source/world_builder/types/contour.cc
@@ -0,0 +1,182 @@
+/*
+ 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/types/contour.h"
+
+
+namespace WorldBuilder
+{
+ namespace Features
+ {
+ namespace FaultModels
+ {
+ namespace Composition
+ {
+ class Interface;
+ } // namespace Composition
+ namespace Grains
+ {
+ class Interface;
+ } // namespace Grains
+ namespace Temperature
+ {
+ class Interface;
+ } // namespace Temperature
+ namespace Velocity
+ {
+ class Interface;
+ } // namespace Velocity
+ } // namespace FaultModels
+ namespace SubductingPlateModels
+ {
+ namespace Composition
+ {
+ class Interface;
+ } // namespace Composition
+ namespace Grains
+ {
+ class Interface;
+ } // namespace Grains
+ namespace Temperature
+ {
+ class Interface;
+ } // namespace Temperature
+ namespace Velocity
+ {
+ class Interface;
+ } // namespace Velocity
+ } // namespace SubductingPlateModels
+ } // namespace Features
+
+ namespace Types
+ {
+ Contour::Contour(const double default_length_,
+ const WorldBuilder::Point<2> &default_thickness_,
+ const WorldBuilder::Point<2> &default_top_truncation_,
+ const WorldBuilder::Point<2> &default_angle_,
+ const Types::Interface &temperature_plugin_system_,
+ const Types::Interface &composition_plugin_system_,
+ const Types::Interface &grains_plugin_system_,
+ const Types::Interface &velocity_plugin_system_)
+ :
+ value_length(default_length_),
+ default_length(default_length_),
+ value_thickness(default_thickness_),
+ default_thickness(default_thickness_),
+ default_top_truncation(default_top_truncation_),
+ value_angle(default_angle_),
+ default_angle(default_angle_),
+ temperature_plugin_system(temperature_plugin_system_.clone()),
+ composition_plugin_system(composition_plugin_system_.clone()),
+ grains_plugin_system(grains_plugin_system_.clone()),
+ velocity_plugin_system(velocity_plugin_system_.clone())
+ {
+ this->type_name = Types::type::Segment;
+ }
+
+
+ Contour::Contour(Contour const &other)
+ :
+ value_length(other.default_length),
+ default_length(other.default_length),
+ value_thickness(other.default_thickness),
+ default_thickness(other.default_thickness),
+ default_top_truncation(other.default_top_truncation),
+ value_angle(other.default_angle),
+ default_angle(other.default_angle),
+ temperature_plugin_system(other.temperature_plugin_system->clone()),
+ composition_plugin_system(other.composition_plugin_system->clone()),
+ grains_plugin_system(other.grains_plugin_system->clone()),
+ velocity_plugin_system(other.velocity_plugin_system->clone())
+ {
+ this->type_name = Types::type::Contour;
+ }
+
+
+ Contour::~Contour ()
+ = default;
+
+
+
+ void
+ Contour::write_schema(Parameters &prm,
+ const std::string &name,
+ const std::string &documentation) const
+ {
+ using namespace rapidjson;
+ prm.enter_subsection(name);
+ {
+ Document &declarations = prm.declarations;
+ std::string base = prm.get_full_json_path();
+
+ Pointer((base + "/type").c_str()).Set(declarations,"object");
+ Pointer((base + "/additionalProperties").c_str()).Set(declarations,false);
+ Pointer((base + "/description").c_str()).Set(declarations,documentation.c_str());
+ std::vector restricted_values = {"length", "thickness", "angle"};
+ for (unsigned int i = 0; i < restricted_values.size(); ++i)
+ {
+ if (!restricted_values[i].empty())
+ {
+ if (i == 0 && Pointer((base + "/required").c_str()).Get(declarations) == nullptr)
+ {
+ // The enum array doesn't exist yet, so we create it and fill it.
+ Pointer((base + "/required/0").c_str()).Create(declarations);
+ Pointer((base + "/required/0").c_str()).Set(declarations, restricted_values[i].c_str());
+ }
+ else
+ {
+ // The enum array already exist yet, so we add an element to the end.
+ Pointer((base + "/required/-").c_str()).Set(declarations, restricted_values[i].c_str());
+ }
+ }
+ }
+
+ prm.enter_subsection("properties");
+ {
+ base = prm.get_full_json_path();
+ Pointer((base + "/length/type").c_str()).Set(declarations,"number");
+
+ Pointer((base + "/thickness/type").c_str()).Set(declarations,"array");
+ Pointer((base + "/thickness/minItems").c_str()).Set(declarations,1);
+ Pointer((base + "/thickness/maxItems").c_str()).Set(declarations,2);
+ Pointer((base + "/thickness/items/type").c_str()).Set(declarations,"number");
+
+ Pointer((base + "/top truncation/type").c_str()).Set(declarations,"array");
+ Pointer((base + "/top truncation/minItems").c_str()).Set(declarations,1);
+ Pointer((base + "/top truncation/maxItems").c_str()).Set(declarations,2);
+ Pointer((base + "/top truncation/items/type").c_str()).Set(declarations,"number");
+
+ Pointer((base + "/angle/type").c_str()).Set(declarations,"array");
+ Pointer((base + "/angle/minItems").c_str()).Set(declarations,1);
+ Pointer((base + "/angle/maxItems").c_str()).Set(declarations,2);
+ Pointer((base + "/angle/items/type").c_str()).Set(declarations,"number");
+
+ temperature_plugin_system->write_schema(prm, "temperature models", "");
+ composition_plugin_system->write_schema(prm, "composition models", "");
+ grains_plugin_system->write_schema(prm, "grains models", "");
+ velocity_plugin_system->write_schema(prm, "velocity models", "");
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+
+ }
+ } // namespace Types
+} // namespace WorldBuilder
+
diff --git a/source/world_builder/utilities.cc b/source/world_builder/utilities.cc
index 6380322fb..736730d2c 100644
--- a/source/world_builder/utilities.cc
+++ b/source/world_builder/utilities.cc
@@ -1506,6 +1506,33 @@ namespace WorldBuilder
}
+ double
+ calculate_distance_from_point_to_rectilinear_plane(const Point<3> p1,
+ const Point<3> p2,
+ const Point<3> p3,
+ const Point<3> check_point)
+ {
+ // First Determine the Plane formed by p1, p2, p3
+ Point<3> P2P1 = p2 - p1;
+ Point<3> P3P1 = p3 - p1;
+ Point<3> normal_vector = cross_product(P2P1, P3P1);
+ Point<3> normal_vector_unit = normal_vector / normal_vector.norm();
+
+ // Determine the equation of the plane using the Normal vector and point p1
+ const double A = normal_vector_unit[0];
+ const double B = normal_vector_unit[1];
+ const double C = normal_vector_unit[2];
+ const double D = -A*p1[0] - B*p1[1] - C*p1[2];
+
+ // Calculate the denominator (squared magnitude of the normal vector)
+ const double denominator = A * A + B * B + C * C;
+
+ // Calculate the signed distance from the point to the plane
+ const double distance = (A * check_point[0] + B * check_point[1] + C * check_point[2] + D) / denominator;
+
+ return distance;
+ }
+
std::array,3>
multiply_3x3_matrices(const std::array,3> mat1, const std::array,3> mat2)
{
diff --git a/source/world_builder/world.cc b/source/world_builder/world.cc
index 5bfd422f4..2b0bba839 100644
--- a/source/world_builder/world.cc
+++ b/source/world_builder/world.cc
@@ -22,6 +22,7 @@
#include "world_builder/config.h"
#include "world_builder/features/subducting_plate.h"
+#include "world_builder/features/subducting_plate_contours.h"
#include "world_builder/gravity_model/interface.h"
#include "world_builder/nan.h"
#include "world_builder/point.h"