diff --git a/include/world_builder/features/subducting_plate.h b/include/world_builder/features/subducting_plate.h index c519738ed..e3cf54c9f 100644 --- a/include/world_builder/features/subducting_plate.h +++ b/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 obliquity_vector; + std::vector > slab_segment_lengths; std::vector > > slab_segment_thickness; std::vector > > slab_segment_top_truncation; diff --git a/include/world_builder/utilities.h b/include/world_builder/utilities.h index 862ef3895..c23eb904e 100644 --- a/include/world_builder/utilities.h +++ b/include/world_builder/utilities.h @@ -420,7 +420,8 @@ namespace WorldBuilder const double start_radius, const std::unique_ptr &coordinate_system, const bool only_positive, - const Objects::BezierCurve &bezier_curve); + const Objects::BezierCurve &bezier_curve, + std::vector obliquity_vector = {std::numeric_limits::infinity(), std::numeric_limits::infinity()}); diff --git a/source/world_builder/features/subducting_plate.cc b/source/world_builder/features/subducting_plate.cc index b8da9327e..0ab071487 100644 --- a/source/world_builder/features/subducting_plate.cc +++ b/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::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 >("dip point"); + obliquity_vector = prm.get_vector("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; diff --git a/source/world_builder/utilities.cc b/source/world_builder/utilities.cc index a31881e4b..2c2750cca 100644 --- a/source/world_builder/utilities.cc +++ b/source/world_builder/utilities.cc @@ -412,7 +412,8 @@ namespace WorldBuilder const double start_radius, const std::unique_ptr &coordinate_system, const bool only_positive, - const Objects::BezierCurve &bezier_curve) + const Objects::BezierCurve &bezier_curve, + std::vector obliquity_vector) { // Initialize variables that this function will return double distance = std::numeric_limits::infinity(); @@ -453,9 +454,182 @@ namespace WorldBuilder double fraction_CPL_P1P2 = std::numeric_limits::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::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(std::round(fraction_CPL_P1P2))][0]-check_point_surface_2d[0]); + const double plus = std::fabs(point_list[i_section_min_distance+static_cast(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(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(std::round(fraction_CPL_P1P2))][0]-check_point_surface_2d[0]); + const double plus = std::fabs(point_list[i_section_min_distance+static_cast(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(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; diff --git a/tests/gwb-dat/linear_cartesian_obliquity.dat b/tests/gwb-dat/linear_cartesian_obliquity.dat new file mode 100644 index 000000000..9977ff1a9 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/tests/gwb-dat/linear_cartesian_obliquity.wb b/tests/gwb-dat/linear_cartesian_obliquity.wb new file mode 100644 index 000000000..b6ea4878e --- /dev/null +++ b/tests/gwb-dat/linear_cartesian_obliquity.wb @@ -0,0 +1,58 @@ +// Composition 0 is slab +// Composition 1 is weakzone +// Composition 2 is overriding +// Composition 3 is fracture_zone + +{ + "version": "1.1", + "gravity model":{"model":"uniform", "magnitude":9.81}, + "surface temperature":273, "potential mantle temperature":1573, + "thermal expansion coefficient":3.0e-5, "thermal diffusivity":1.0e-6, + "features": + [ + {"model": "mantle layer", "name": "mantle", "coordinates": [[-500e3, 10000e3], [-500e3, -10000e3], [16000e3, -10000e3], [16000e3, 10000e3]], + "min depth": 0, "max depth":10000e3, + "temperature models": + [{"model":"uniform", "min depth": 0, "max depth": 10000e3, "temperature":1573}] + }, + + {"model": "oceanic plate", "name": "Subducting Plate", + "coordinates": [[3000e3, -10000e3], [3000e3, 500e3], + [2500e3, 1000e3], [2500e3, 10000e3], + [500e3, 10000e3], [500e3, -10000e3]], + "min depth": 0, + "max depth": 150e3, + "composition models": + [{"model": "uniform", "compositions": [0], "min depth":0, "max depth": 150e3, "operation":"replace"}], + + "temperature models": [{"model": "plate model", "ridge coordinates": [[[500e3,-10000e3],[500e3,10000e3]]], + "spreading velocity": 0.03, "min depth": 0, "max depth": 300e3, "bottom temperature": 1573}] + }, + + + {"model":"subducting plate", "name":"Slab", + + "coordinates":[[2500e3, 1000e3],[3000e3, 500e3]], + "dip point":[1e10,750e3],"max depth":1000e3, "obliquity vector":[1, 0], + "segments":[{"length":300e3,"thickness":[150e3],"top truncation":[-50e3],"angle":[0,70]}], + + "composition models":[ + {"model":"uniform", "compositions":[1], "min distance slab top":0, "max distance slab top":20e3, "operation":"replace"}, + {"model":"uniform", "compositions":[0], "min distance slab top":20e3, "max distance slab top":150e3, "operation":"replace"}], + + "temperature models":[{"model":"mass conserving", + "reference model name": "plate model", + "density":3300, + "thermal conductivity":3.3, + "adiabatic heating":false, + "subducting velocity":0.03, + "spreading velocity":0.03, + "ridge coordinates":[[[500e3,-10000e3],[500e3,10000e3]]], + "coupling depth":80e3, + "forearc cooling factor":1.0, + "taper distance":10e3, + "min distance slab top":-200e3, + "max distance slab top":300e3}] + } + ] +} diff --git a/tests/gwb-dat/linear_cartesian_obliquity/screen-output.log b/tests/gwb-dat/linear_cartesian_obliquity/screen-output.log new file mode 100644 index 000000000..42d834478 --- /dev/null +++ b/tests/gwb-dat/linear_cartesian_obliquity/screen-output.log @@ -0,0 +1,6 @@ +# x y z d g T vx vy vz tag +3000e3 510e3 0e3 20e3 550.128 0 0 0 2 +3050e3 510e3 0e3 30e3 573.728 0 0 0 2 +3100e3 510e3 0e3 50e3 559.492 0 0 0 2 +3150e3 510e3 0e3 90e3 567.813 0 0 0 2 +3200e3 510e3 0e3 150e3 504.47 0 0 0 2