[WIP] Allow for Oblique Slabs/Faults Relative to Surface

include/world_builder/features/subducting_plate.h

@@ -217,6 +217,13 @@ namespace WorldBuilder
          */
         Point<2> reference_point;
 
+        /**
+         * A vector on the surface that indicates the direction of convergence
+         * of the subducting plate relative to the trench.
+         */
+        // Point<2> obliquity_vector;
+        std::vector<double> obliquity_vector;
+
         std::vector<std::vector<double> > slab_segment_lengths;
         std::vector<std::vector<Point<2> > > slab_segment_thickness;
         std::vector<std::vector<Point<2> > > slab_segment_top_truncation;

include/world_builder/utilities.h

@@ -420,7 +420,8 @@ namespace WorldBuilder
                                                                     const double start_radius,
                                                                     const std::unique_ptr<CoordinateSystems::Interface> &coordinate_system,
                                                                     const bool only_positive,
-                                                                    const Objects::BezierCurve &bezier_curve);
+                                                                    const Objects::BezierCurve &bezier_curve,
+                                                                    std::vector<double> obliquity_vector = {std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity()});
 
 
 

source/world_builder/features/subducting_plate.cc

@@ -106,7 +106,8 @@ namespace WorldBuilder
                         "The depth to which this feature is present");
       prm.declare_entry("dip point", Types::Point<2>(),
                         "The depth to which this feature is present");
-
+      prm.declare_entry("obliquity vector", Types::Array(Types::Double(std::numeric_limits<double>::infinity()),2),
+                        "A vector on the surface that indicates the direction of convergence of the subducting plate relative to the trench.");
       /*prm.declare_entry("segments", Types::Array(Types::Segment(0,Point<2>(0,0,invalid),Point<2>(0,0,invalid),Point<2>(0,0,invalid),
                                                                 Types::PluginSystem("", Features::SubductingPlateModels::Temperature::Interface::declare_entries, {"model"}),
                                                                 Types::PluginSystem("", Features::SubductingPlateModels::Composition::Interface::declare_entries, {"model"}),
@@ -172,12 +173,17 @@ namespace WorldBuilder
 
       reference_point = prm.get<Point<2> >("dip point");
 
+      obliquity_vector = prm.get_vector<double>("obliquity vector");
+
       if (coordinate_system == spherical)
         {
           // When spherical, input is in degrees, so change to radians for internal use.
           reference_point *= (Consts::PI/180.);
+          for (double &value : obliquity_vector)
+            value *= (Consts::PI/180.);
         }
 
+
       default_temperature_models.resize(0);
       default_composition_models.resize(0);
       default_grains_models.resize(0);
@@ -544,7 +550,8 @@ namespace WorldBuilder
                                                                        starting_radius,
                                                                        this->world->parameters.coordinate_system,
                                                                        false,
-                                                                       this->bezier_curve);
+                                                                       this->bezier_curve,
+                                                                       obliquity_vector);
 
           const double distance_from_plane = distance_from_planes.distance_from_plane;
           const double distance_along_plane = distance_from_planes.distance_along_plane;
@@ -845,7 +852,8 @@ namespace WorldBuilder
                                                                    starting_radius,
                                                                    this->world->parameters.coordinate_system,
                                                                    false,
-                                                                   this->bezier_curve);
+                                                                   this->bezier_curve,
+                                                                   obliquity_vector);
 
       Objects::PlaneDistances plane_distances(distance_from_planes.distance_from_plane, distance_from_planes.distance_along_plane);
       return plane_distances;

source/world_builder/utilities.cc

@@ -412,7 +412,8 @@ namespace WorldBuilder
                                       const double start_radius,
                                       const std::unique_ptr<CoordinateSystems::Interface> &coordinate_system,
                                       const bool only_positive,
-                                      const Objects::BezierCurve &bezier_curve)
+                                      const Objects::BezierCurve &bezier_curve,
+                                      std::vector<double> obliquity_vector)
     {
       // Initialize variables that this function will return
       double distance = std::numeric_limits<double>::infinity();
@@ -453,9 +454,182 @@ namespace WorldBuilder
       double fraction_CPL_P1P2 = std::numeric_limits<double>::infinity();
 
       // Get the closest point on the curve (the trench or the fault trace) to the check point.
-      const Objects::ClosestPointOnCurve closest_point_on_curve = bezier_curve.closest_point_on_curve_segment(check_point_surface_2d);
+      Objects::ClosestPointOnCurve closest_point_on_curve = bezier_curve.closest_point_on_curve_segment(check_point_surface_2d);
       Point<2> closest_point_on_line_2d = closest_point_on_curve.point;
 
+      // If the obliquity_vector is not infinite, this means that the user has specified an obliquity vector.
+      // This will require a potential modification of the `closest_point_on_line_2d` variable. We will take
+      // the check point and use it with the obliquity vector to parameterize a line and see where this line
+      // intersects the Bezier curve. If it does intersect the Bezier curve, then this intersection point will
+      // become the new `closest_point_on_line_2d`, otherwise it will be set to NaN.
+      // if (obliquity_vector_point.norm() > std::numeric_limits<double>::min())
+      if (std::isfinite(obliquity_vector[0]) && std::isfinite(obliquity_vector[1]))
+        {
+          Point<2> obliquity_vector_point(obliquity_vector[0], obliquity_vector[1], natural_coordinate_system);
+          // Normalize the vector
+          obliquity_vector_point = obliquity_vector_point / obliquity_vector_point.norm();
+
+          // If the check point is near the ends of the Bezier curve, there is a chance that the Bezier curve will not
+          // find a closest point. Check for this case by seeing if the line created by the checkpoint and the obliquity
+          // vector intersects the end segments of the Bezier curve.
+          Point<2> iterable_check_point_surface_2d = check_point_surface_2d;
+          Objects::ClosestPointOnCurve iterable_closest_point_on_curve = bezier_curve.closest_point_on_curve_segment(iterable_check_point_surface_2d);
+
+          if (std::isnan(closest_point_on_line_2d[0]))
+            {
+              // The closest point on the bezier curve has not been found. Extract the terminal points of the bezier curve
+              // and see which is closer to the checkpoint.
+              const Point<2> p1 = point_list[0];
+              const Point<2> p2 = point_list[point_list.size()-1];
+              const double dist_p1 = std::abs(check_point_surface_2d.distance(p1));
+              const double dist_p2 = std::abs(check_point_surface_2d.distance(p2));
+              const Point<2> closest_terminal_point = dist_p1 < dist_p2 ? p1 : p2;
+              double closest_distance_to_trench = dist_p1 < dist_p2 ? dist_p1 : dist_p2;
+
+              // Because the bezier curve geometry can be very complex, still track the second closest point. It's unlikely
+              // that the second closest point will end up closer to the check point than the first closest point during
+              // the iterations, but this may help with edge cases.
+              const Point<2> second_closest_terminal_point = dist_p1 < dist_p2 ? point_list[1] : point_list[point_list.size()-2];
+              for (unsigned int i = 0; i < 100; ++i)
+                {
+                  Point<2> check_point_surface_2d_temp_ob = iterable_check_point_surface_2d;
+
+                  if (!bool_cartesian)
+                    {
+                      const double normal = std::fabs(point_list[i_section_min_distance+static_cast<size_t>(std::round(fraction_CPL_P1P2))][0]-check_point_surface_2d[0]);
+                      const double plus   = std::fabs(point_list[i_section_min_distance+static_cast<size_t>(std::round(fraction_CPL_P1P2))][0]-(check_point_surface_2d[0]+2*Consts::PI));
+                      const double min    = std::fabs(point_list[i_section_min_distance+static_cast<size_t>(std::round(fraction_CPL_P1P2))][0]-(check_point_surface_2d[0]-2*Consts::PI));
+
+                      // find out whether the check point, checkpoint + 2pi or check point -2 pi is closest to the point list.
+                      if (plus < normal)
+                        {
+                          check_point_surface_2d_temp_ob[0]+= 2*Consts::PI;
+                        }
+                      else if (min < normal)
+                        {
+                          check_point_surface_2d_temp_ob[0]-= 2*Consts::PI;
+                        }
+                    }
+
+                  // Check whether the obliquity vector points toward the curve from the current check point.
+                  const Point<2> vector_checkpoint_to_curve = closest_terminal_point - check_point_surface_2d_temp_ob;
+                  const double line_factor = (obliquity_vector_point * vector_checkpoint_to_curve) >= 0.0 ? 1.0 : -1.0;
+
+                  // Parameterize a line through the checkpoint along the obliquity vector.
+                  Point<2> parameterized_line(iterable_check_point_surface_2d[0] + line_factor * closest_distance_to_trench * obliquity_vector_point[0],
+                                              iterable_check_point_surface_2d[1] + line_factor * closest_distance_to_trench * obliquity_vector_point[1],
+                                              natural_coordinate_system);
+
+                  const double new_distance1 = std::abs(parameterized_line.distance(closest_terminal_point));
+                  const double new_distance2 = std::abs(parameterized_line.distance(second_closest_terminal_point));
+                  const double new_distance_search = std::min(new_distance1, new_distance2);
+
+                  if (new_distance_search < closest_distance_to_trench)
+                    {
+                      // We are getting closer to the curve, update the closest point and distance.
+                      iterable_check_point_surface_2d = parameterized_line;
+                      closest_distance_to_trench = new_distance_search;
+
+                      if (!std::isnan(bezier_curve.closest_point_on_curve_segment(iterable_check_point_surface_2d).point[0]))
+                        {
+                          // We have found a point on the bezier curve. Set this as the checkpoints closest point on the curve.
+                          closest_point_on_curve = bezier_curve.closest_point_on_curve_segment(iterable_check_point_surface_2d);
+                          closest_point_on_line_2d = closest_point_on_curve.point;
+                          iterable_closest_point_on_curve = closest_point_on_curve;
+                          break;
+                        }
+                    }
+
+                  // We are getting farther from the curve, stop the search.
+                  else
+                    break;
+                }
+            }
+
+          // Now check if the bezier_curve has found a point on the curve that is closest to the checkpoint (without considering the obliquity vector).
+          // If so, we need to check whether this point is on the correct side of the curve, given the reference point.
+          // If the bezier_curve has not found a point on the curve that is closest to the checkpoint, check to see if the point will still
+          // intersect the trench given the obliquity_vector.
+          if (!std::isnan(closest_point_on_line_2d[0]))
+            {
+              Point<2> check_point_surface_2d_temp_ob = iterable_check_point_surface_2d;
+
+              if (!bool_cartesian)
+                {
+                  const double normal = std::fabs(point_list[i_section_min_distance+static_cast<size_t>(std::round(fraction_CPL_P1P2))][0]-check_point_surface_2d[0]);
+                  const double plus   = std::fabs(point_list[i_section_min_distance+static_cast<size_t>(std::round(fraction_CPL_P1P2))][0]-(check_point_surface_2d[0]+2*Consts::PI));
+                  const double min    = std::fabs(point_list[i_section_min_distance+static_cast<size_t>(std::round(fraction_CPL_P1P2))][0]-(check_point_surface_2d[0]-2*Consts::PI));
+
+                  // find out whether the check point, checkpoint + 2pi or check point -2 pi is closest to the point list.
+                  if (plus < normal)
+                    {
+                      check_point_surface_2d_temp_ob[0]+= 2*Consts::PI;
+                    }
+                  else if (min < normal)
+                    {
+                      check_point_surface_2d_temp_ob[0]-= 2*Consts::PI;
+                    }
+                }
+
+              // Check whether the obliquity vector points toward the curve from the current check point.
+              const Point<2> vector_checkpoint_to_curve = closest_point_on_line_2d - check_point_surface_2d_temp_ob;
+              const double line_factor = (obliquity_vector_point * vector_checkpoint_to_curve) >= 0.0 ? 1.0 : -1.0;
+
+              bool converged = false;
+              // Move the checkpoint along the obliquity vector and search for a converged
+              // intersection with the curve.
+              for (unsigned int i = 0; i < 100; ++i)
+                {
+                  const double old_dist = std::abs(iterable_closest_point_on_curve.distance);
+
+                  Point<2> parameterized_line(iterable_check_point_surface_2d[0] + line_factor * old_dist * obliquity_vector_point[0],
+                                              iterable_check_point_surface_2d[1] + line_factor * old_dist * obliquity_vector_point[1],
+                                              natural_coordinate_system);
+
+                  iterable_closest_point_on_curve = bezier_curve.closest_point_on_curve_segment(parameterized_line);
+                  iterable_check_point_surface_2d = parameterized_line;
+                  const double new_dist = std::abs(iterable_closest_point_on_curve.distance);
+
+                  const double convergence_tol = bool_cartesian ? 1.0 : 1e-6;
+
+                  // We are getting farther away, break the loop.
+                  if (new_dist > old_dist)
+                    break;
+
+                  // We converged to a point on the curve, break the loop, break the loop and update the check point.
+                  if (new_dist < convergence_tol)
+                    {
+                      converged = true;
+                      break;
+                    }
+                }
+
+              // If we converged, update the closest point on the bezier curve.
+              if (converged)
+                {
+                  closest_point_on_curve   = bezier_curve.closest_point_on_curve_segment(iterable_check_point_surface_2d);
+                  closest_point_on_line_2d = closest_point_on_curve.point;
+                }
+
+              // If we did not converge, assign the closest point as NaN.
+              else
+                {
+                  closest_point_on_line_2d[0] = NaN::DQNAN;
+                  closest_point_on_line_2d[1] = NaN::DQNAN;
+                }
+
+
+
+
+
+
+
+
+
+
+            }
+        }
+
       // We now need 3d points from this point on, so make them.
       // The order of a Cartesian coordinate is x,y,z and the order of
       // a spherical coordinate it radius, long, lat (in rad).
@@ -664,7 +838,7 @@ namespace WorldBuilder
 
               // check whether the check point and the reference point are on the same side, if not, change the side.
               const Point<2> AB_normal = closest_point_on_curve.normal*closest_point_on_line_2d.distance(reference_point);//*AB.norm();
-              const Point<2> local_reference_point = AB_normal*1.+closest_point_on_line_2d;
+              const Point<2> local_reference_point = AB_normal + closest_point_on_line_2d;
               const bool reference_normal_on_side_of_line =  (closest_point_on_line_2d-local_reference_point).norm_square() < (check_point_surface_2d_temp-local_reference_point).norm_square();
               const bool reference_point_on_side_of_line =  (point_list[point_list.size()-1][0] - point_list[0][0])*(reference_point[1] - point_list[0][1]) - (reference_point[0] - point_list[0][0])*(point_list[point_list.size()-1][1] - point_list[0][1]) < 0.;
               const double reference_on_side_of_line =  reference_normal_on_side_of_line == reference_point_on_side_of_line ? 1 : -1;

tests/gwb-dat/linear_cartesian_obliquity.dat

@@ -0,0 +1,11 @@
+# This is a comment in the data
+# file. 
+# Now define parameters:
+# dim = 3
+# compositions = 0
+# x y z d T C0
+3000e3 510e3 0e3 20e3
+3050e3 510e3 0e3 30e3
+3100e3 510e3 0e3 50e3
+3150e3 510e3 0e3 90e3
+3200e3 510e3 0e3 150e3