From 7840ac68770ceac08b4a85b5db14b81ab48f44a7 Mon Sep 17 00:00:00 2001 From: asaxena Date: Fri, 11 Jun 2021 14:06:39 -0500 Subject: [PATCH 1/3] Changes to include a bounding box in computing fault composition. PR comments, not addressed completely. PR comments addressed, fixed the previous bugs. Added a larger buffer zone. --- include/world_builder/features/fault.h | 2 + include/world_builder/utilities.h | 18 ++++++++- source/features/fault.cc | 10 ++++- source/utilities.cc | 54 ++++++++++++++++++++++++++ 4 files changed, 82 insertions(+), 2 deletions(-) diff --git a/include/world_builder/features/fault.h b/include/world_builder/features/fault.h index e1ddd3356..aad6780b2 100644 --- a/include/world_builder/features/fault.h +++ b/include/world_builder/features/fault.h @@ -179,6 +179,8 @@ namespace WorldBuilder double maximum_total_fault_length; double maximum_fault_thickness; + std::vector > bounding_box_coordinates; + }; } // namespace Features diff --git a/include/world_builder/utilities.h b/include/world_builder/utilities.h index 8f0a286ce..bd2d5ca3f 100644 --- a/include/world_builder/utilities.h +++ b/include/world_builder/utilities.h @@ -65,8 +65,24 @@ namespace WorldBuilder double signed_distance_to_polygon(const std::vector > &point_list_, const Point<2> &point_); + /** + * Takes the list of the fault coordinates and computes a bounding box + * based on the minimum and maximum x,y values in the list, and an + * additional buffer zone computed using the maximum fault length. + * Currenly, only implemented for the spherical coordinates. + */ + std::vector > + get_bounding_box (const Point<2> &point, + const std::vector > &point_list, + std::size_t original_number_of_coordinates, + const double buffer_around_fault); - + /** + * Takes the list of coordinates and the maximum fault extent to define the + * bounding box, and returns true if the point lies inside the box and false otherwise. + */ + bool bounding_box_contains_point (const std::vector > &point_list, + const Point<2> &point); /* * A class that represents a point in a chosen coordinate system. */ diff --git a/source/features/fault.cc b/source/features/fault.cc index 3bed5bbd2..660bb9ba6 100644 --- a/source/features/fault.cc +++ b/source/features/fault.cc @@ -26,6 +26,7 @@ #include "world_builder/types/point.h" #include "world_builder/types/unsigned_int.h" #include "world_builder/world.h" +#include using namespace std; @@ -325,6 +326,10 @@ namespace WorldBuilder total_fault_length[i] = local_total_fault_length; maximum_total_fault_length = std::max(maximum_total_fault_length, local_total_fault_length); } + + const double buffer_zone = maximum_fault_thickness + maximum_fault_thickness; + bounding_box_coordinates = WorldBuilder::Utilities::get_bounding_box ( + reference_point, coordinates, original_number_of_coordinates, buffer_zone); } @@ -468,6 +473,7 @@ namespace WorldBuilder return temperature; } + double Fault::composition(const Point<3> &position, const NaturalCoordinate &position_in_natural_coordinates, @@ -479,7 +485,9 @@ namespace WorldBuilder const double starting_radius = position_in_natural_coordinates.get_depth_coordinate() + depth - starting_depth; // todo: explain and check -starting_depth - if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness) + if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness + && WorldBuilder::Utilities::bounding_box_contains_point (bounding_box_coordinates, Point<2>(natural_coordinate.get_surface_coordinates(), + world->parameters.coordinate_system->natural_coordinate_system())) ) { // todo: explain // This function only returns positive values, because we want diff --git a/source/utilities.cc b/source/utilities.cc index 0bddb67cf..636dbc59d 100644 --- a/source/utilities.cc +++ b/source/utilities.cc @@ -207,6 +207,60 @@ namespace WorldBuilder } + std::vector > + get_bounding_box (const Point<2> &point, + const std::vector > &point_list, + std::size_t original_number_of_coordinates, + const double buffer_around_fault) + { + std::vector x_list(original_number_of_coordinates,0.0); + std::vector y_list(original_number_of_coordinates,0.0); + + for (size_t j=0; j > bounding_box (4, point); + + // The bounding box currently is hard-coded to include buffer zone based on the spherical geometry. + const double bounding_buffer_zone = 2.* (buffer_around_fault) / (const_pi * 2. * 6371000.) ; + + // The bounding box corners defined from (x0, y0) in clockwise direction + bounding_box[0][0] = *std::min_element(x_list.begin(), x_list.end()) - bounding_buffer_zone; + bounding_box[0][1] = *std::min_element(y_list.begin(), y_list.end()) - bounding_buffer_zone; + + bounding_box[1][0] = *std::min_element(x_list.begin(), x_list.end()) - bounding_buffer_zone; + bounding_box[1][1] = *std::max_element(y_list.begin(), y_list.end()) + bounding_buffer_zone; + + bounding_box[2][0] = *std::max_element(x_list.begin(), x_list.end()) + bounding_buffer_zone; + bounding_box[2][1] = *std::max_element(y_list.begin(), y_list.end()) + bounding_buffer_zone; + + bounding_box[3][0] = *std::max_element(x_list.begin(), x_list.end()) + bounding_buffer_zone; + bounding_box[3][1] = *std::min_element(y_list.begin(), y_list.end()) - bounding_buffer_zone; + + return bounding_box; + } + + + bool + bounding_box_contains_point (const std::vector > &point_list, + const Point<2> &point) + { + for (size_t i=0; i point_list[0][0] && point[1] < point_list[1][1] && point[1] > point_list[0][1]) + return true; + else + return false; + } + return true; + } + + NaturalCoordinate::NaturalCoordinate(const std::array &position, const CoordinateSystems::Interface &coordinate_system_) { From f27590fe45db4ce21fca2746f65fe7be5aedc837 Mon Sep 17 00:00:00 2001 From: asaxena Date: Wed, 16 Jun 2021 10:50:36 -0500 Subject: [PATCH 2/3] Added bounding box in subducting plate. Added comments about bounding box. Include the computation of buffer zone in the composition phase. Changes to use two points in calculating bounding box. Addressed suggestions from the PR. Use the bounding box class to compute the spherical bounding box. Rewritten the bounding box with more accurate computation and bounding box class. Moving bounding box computation around. Pre definition of most variables in parse stage. Addressed recent suggestions on the PR for better performance. --- include/world_builder/features/fault.h | 22 +++++- .../world_builder/features/subducting_plate.h | 19 +++++ include/world_builder/utilities.h | 18 +---- source/features/fault.cc | 69 ++++++++++++++++--- source/features/subducting_plate.cc | 69 ++++++++++++++++++- source/utilities.cc | 54 --------------- 6 files changed, 167 insertions(+), 84 deletions(-) diff --git a/include/world_builder/features/fault.h b/include/world_builder/features/fault.h index aad6780b2..e26137254 100644 --- a/include/world_builder/features/fault.h +++ b/include/world_builder/features/fault.h @@ -25,6 +25,7 @@ #include "world_builder/features/fault_models/grains/interface.h" #include "world_builder/features/fault_models/temperature/interface.h" #include "world_builder/types/segment.h" +#include "world_builder/bounding_box.h" namespace WorldBuilder @@ -104,6 +105,18 @@ namespace WorldBuilder const unsigned int composition_number, double composition_value) const override final; + /** + * Computes 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 + * correponding latitude points. + */ + + BoundingBox<2> get_bounding_box (const WorldBuilder::Utilities::NaturalCoordinate &position_in_natural_coordinates, + const double depth) const; + /** * Returns a grains (rotation matrix and grain size) based on the * given position, depth in the model, the composition (e.g. representing @@ -179,8 +192,13 @@ namespace WorldBuilder double maximum_total_fault_length; double maximum_fault_thickness; - std::vector > bounding_box_coordinates; - + 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 buffer_around_fault_cartesian; }; } // namespace Features diff --git a/include/world_builder/features/subducting_plate.h b/include/world_builder/features/subducting_plate.h index f06097b77..4889cc979 100644 --- a/include/world_builder/features/subducting_plate.h +++ b/include/world_builder/features/subducting_plate.h @@ -25,6 +25,7 @@ #include "world_builder/features/subducting_plate_models/grains/interface.h" #include "world_builder/features/subducting_plate_models/temperature/interface.h" #include "world_builder/types/segment.h" +#include "world_builder/bounding_box.h" namespace WorldBuilder @@ -83,6 +84,17 @@ namespace WorldBuilder void parse_entries(Parameters &prm) override final; + /** + * Computes 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 + * correponding latitude points. + */ + BoundingBox<2> get_bounding_box (const WorldBuilder::Utilities::NaturalCoordinate &position_in_natural_coordinates, + const double depth) const; + /** * Returns a temperature based on the given position, depth in the model, @@ -182,6 +194,13 @@ namespace WorldBuilder double maximum_total_slab_length; double maximum_slab_thickness; + 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 buffer_around_slab_cartesian; }; } // namespace Features } // namespace WorldBuilder diff --git a/include/world_builder/utilities.h b/include/world_builder/utilities.h index bd2d5ca3f..8f0a286ce 100644 --- a/include/world_builder/utilities.h +++ b/include/world_builder/utilities.h @@ -65,24 +65,8 @@ namespace WorldBuilder double signed_distance_to_polygon(const std::vector > &point_list_, const Point<2> &point_); - /** - * Takes the list of the fault coordinates and computes a bounding box - * based on the minimum and maximum x,y values in the list, and an - * additional buffer zone computed using the maximum fault length. - * Currenly, only implemented for the spherical coordinates. - */ - std::vector > - get_bounding_box (const Point<2> &point, - const std::vector > &point_list, - std::size_t original_number_of_coordinates, - const double buffer_around_fault); - /** - * Takes the list of coordinates and the maximum fault extent to define the - * bounding box, and returns true if the point lies inside the box and false otherwise. - */ - bool bounding_box_contains_point (const std::vector > &point_list, - const Point<2> &point); + /* * A class that represents a point in a chosen coordinate system. */ diff --git a/source/features/fault.cc b/source/features/fault.cc index 660bb9ba6..d315c883f 100644 --- a/source/features/fault.cc +++ b/source/features/fault.cc @@ -327,9 +327,58 @@ namespace WorldBuilder maximum_total_fault_length = std::max(maximum_total_fault_length, local_total_fault_length); } - const double buffer_zone = maximum_fault_thickness + maximum_fault_thickness; - bounding_box_coordinates = WorldBuilder::Utilities::get_bounding_box ( - reference_point, coordinates, original_number_of_coordinates, buffer_zone); + + std::vector x_list(original_number_of_coordinates,0.0); + std::vector y_list(original_number_of_coordinates,0.0); + + for (size_t j=0; j + Fault::get_bounding_box (const NaturalCoordinate &position_in_natural_coordinates, + const double depth) const + { + BoundingBox<2> surface_bounding_box; + const double starting_radius_inv = 1 / (position_in_natural_coordinates.get_depth_coordinate() + depth - starting_depth); + if (world->parameters.coordinate_system->natural_coordinate_system() == CoordinateSystem::spherical) + { + std::pair, Point<2> > &spherical_bounding_box = surface_bounding_box.get_boundary_points(); + + const double buffer_around_fault_spherical = 2 * const_pi * buffer_around_fault_cartesian * starting_radius_inv; + + spherical_bounding_box.first = {(min_along_x - buffer_around_fault_spherical * min_lat_cos_inv) , + (min_along_y - buffer_around_fault_spherical), spherical + } ; + + spherical_bounding_box.second = {(max_along_x + buffer_around_fault_spherical * max_lat_cos_inv) , + (max_along_y + buffer_around_fault_spherical), spherical + }; + } + else if (world->parameters.coordinate_system->natural_coordinate_system() == CoordinateSystem::cartesian) + { + std::pair, Point<2> > &bounding_box = surface_bounding_box.get_boundary_points(); + bounding_box.first = {min_along_x, min_along_y, cartesian}; + bounding_box.second = {max_along_x, max_along_x, cartesian}; + surface_bounding_box.extend(buffer_around_fault_cartesian); + } + return surface_bounding_box; } @@ -353,7 +402,9 @@ namespace WorldBuilder ); // todo: explain and check -starting_depth - if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness) + if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness && + get_bounding_box(position_in_natural_coordinates, depth).point_inside(Point<2>(position_in_natural_coordinates.get_surface_coordinates(), + world->parameters.coordinate_system->natural_coordinate_system()))) { // todo: explain // This function only returns positive values, because we want @@ -485,9 +536,9 @@ namespace WorldBuilder const double starting_radius = position_in_natural_coordinates.get_depth_coordinate() + depth - starting_depth; // todo: explain and check -starting_depth - if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness - && WorldBuilder::Utilities::bounding_box_contains_point (bounding_box_coordinates, Point<2>(natural_coordinate.get_surface_coordinates(), - world->parameters.coordinate_system->natural_coordinate_system())) ) + if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness && + get_bounding_box(position_in_natural_coordinates, depth).point_inside(Point<2>(position_in_natural_coordinates.get_surface_coordinates(), + world->parameters.coordinate_system->natural_coordinate_system()))) { // todo: explain // This function only returns positive values, because we want @@ -621,7 +672,9 @@ namespace WorldBuilder const double starting_radius = position_in_natural_coordinates.get_depth_coordinate() + depth - starting_depth; // todo: explain and check -starting_depth - if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness) + if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_fault_length + maximum_fault_thickness && + get_bounding_box(position_in_natural_coordinates, depth).point_inside(Point<2>(position_in_natural_coordinates.get_surface_coordinates(), + world->parameters.coordinate_system->natural_coordinate_system()))) { // todo: explain // This function only returns positive values, because we want diff --git a/source/features/subducting_plate.cc b/source/features/subducting_plate.cc index db79146ee..bc42bce89 100644 --- a/source/features/subducting_plate.cc +++ b/source/features/subducting_plate.cc @@ -26,6 +26,7 @@ #include "world_builder/types/point.h" #include "world_builder/types/unsigned_int.h" #include "world_builder/world.h" +#include using namespace std; @@ -329,6 +330,62 @@ namespace WorldBuilder total_slab_length[i] = local_total_slab_length; maximum_total_slab_length = std::max(maximum_total_slab_length, local_total_slab_length); } + + std::vector x_list(original_number_of_coordinates,0.0); + std::vector y_list(original_number_of_coordinates,0.0); + + for (size_t j=0; j + SubductingPlate::get_bounding_box (const NaturalCoordinate &position_in_natural_coordinates, + const double depth) const + { + BoundingBox<2> surface_bounding_box; + const double starting_radius_inv = 1 / (position_in_natural_coordinates.get_depth_coordinate() + depth - starting_depth); + if (world->parameters.coordinate_system->natural_coordinate_system() == CoordinateSystem::spherical) + { + std::pair, Point<2> > &spherical_bounding_box = surface_bounding_box.get_boundary_points(); + + const double buffer_around_fault_spherical = 2 * const_pi * buffer_around_slab_cartesian * starting_radius_inv; + + spherical_bounding_box.first = {(min_along_x - buffer_around_fault_spherical * min_lat_cos_inv) , + (min_along_y - buffer_around_fault_spherical), spherical + } ; + + spherical_bounding_box.second = {(max_along_x + buffer_around_fault_spherical * max_lat_cos_inv) , + (max_along_y + buffer_around_fault_spherical), spherical + }; + } + else if (world->parameters.coordinate_system->natural_coordinate_system() == CoordinateSystem::cartesian) + { + std::pair, Point<2> > &bounding_box = surface_bounding_box.get_boundary_points(); + bounding_box.first = {min_along_x, min_along_y, cartesian}; + bounding_box.second = {max_along_x, max_along_x, cartesian}; + surface_bounding_box.extend(buffer_around_slab_cartesian); + } + return surface_bounding_box; } @@ -352,7 +409,9 @@ namespace WorldBuilder ); // todo: explain and check -starting_depth - if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_slab_length + maximum_slab_thickness) + if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_slab_length + maximum_slab_thickness && + get_bounding_box(position_in_natural_coordinates, depth).point_inside(Point<2>(position_in_natural_coordinates.get_surface_coordinates(), + world->parameters.coordinate_system->natural_coordinate_system()))) { /*WBAssert(coordinates.size() == slab_segment_lengths.size(), "Internal error: The size of coordinates (" << coordinates.size() @@ -489,7 +548,9 @@ namespace WorldBuilder const double starting_radius = position_in_natural_coordinates.get_depth_coordinate() + depth - starting_depth; // todo: explain and check -starting_depth - if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_slab_length + maximum_slab_thickness) + if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_slab_length + maximum_slab_thickness && + get_bounding_box(position_in_natural_coordinates, depth).point_inside(Point<2>(position_in_natural_coordinates.get_surface_coordinates(), + world->parameters.coordinate_system->natural_coordinate_system()))) { // todo: explain WorldBuilder::Utilities::PointDistanceFromCurvedPlanes distance_from_planes = @@ -620,7 +681,9 @@ namespace WorldBuilder const double starting_radius = position_in_natural_coordinates.get_depth_coordinate() + depth - starting_depth; // todo: explain and check -starting_depth - if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_slab_length + maximum_slab_thickness) + if (depth <= maximum_depth && depth >= starting_depth && depth <= maximum_total_slab_length + maximum_slab_thickness && + get_bounding_box(position_in_natural_coordinates, depth).point_inside(Point<2>(position_in_natural_coordinates.get_surface_coordinates(), + world->parameters.coordinate_system->natural_coordinate_system()))) { // todo: explain WorldBuilder::Utilities::PointDistanceFromCurvedPlanes distance_from_planes = diff --git a/source/utilities.cc b/source/utilities.cc index 636dbc59d..0bddb67cf 100644 --- a/source/utilities.cc +++ b/source/utilities.cc @@ -207,60 +207,6 @@ namespace WorldBuilder } - std::vector > - get_bounding_box (const Point<2> &point, - const std::vector > &point_list, - std::size_t original_number_of_coordinates, - const double buffer_around_fault) - { - std::vector x_list(original_number_of_coordinates,0.0); - std::vector y_list(original_number_of_coordinates,0.0); - - for (size_t j=0; j > bounding_box (4, point); - - // The bounding box currently is hard-coded to include buffer zone based on the spherical geometry. - const double bounding_buffer_zone = 2.* (buffer_around_fault) / (const_pi * 2. * 6371000.) ; - - // The bounding box corners defined from (x0, y0) in clockwise direction - bounding_box[0][0] = *std::min_element(x_list.begin(), x_list.end()) - bounding_buffer_zone; - bounding_box[0][1] = *std::min_element(y_list.begin(), y_list.end()) - bounding_buffer_zone; - - bounding_box[1][0] = *std::min_element(x_list.begin(), x_list.end()) - bounding_buffer_zone; - bounding_box[1][1] = *std::max_element(y_list.begin(), y_list.end()) + bounding_buffer_zone; - - bounding_box[2][0] = *std::max_element(x_list.begin(), x_list.end()) + bounding_buffer_zone; - bounding_box[2][1] = *std::max_element(y_list.begin(), y_list.end()) + bounding_buffer_zone; - - bounding_box[3][0] = *std::max_element(x_list.begin(), x_list.end()) + bounding_buffer_zone; - bounding_box[3][1] = *std::min_element(y_list.begin(), y_list.end()) - bounding_buffer_zone; - - return bounding_box; - } - - - bool - bounding_box_contains_point (const std::vector > &point_list, - const Point<2> &point) - { - for (size_t i=0; i point_list[0][0] && point[1] < point_list[1][1] && point[1] > point_list[0][1]) - return true; - else - return false; - } - return true; - } - - NaturalCoordinate::NaturalCoordinate(const std::array &position, const CoordinateSystems::Interface &coordinate_system_) { From 54e3d2ada2c1df421b5ebc17a7a2478f210f79b3 Mon Sep 17 00:00:00 2001 From: asaxena Date: Thu, 15 Jul 2021 21:12:32 -0500 Subject: [PATCH 3/3] Added the change log entry. Fixed the heading for changelog. Fixed the link in the changelog. --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b2ddcbb9e..a5cecb8e1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,7 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm - The World Builder Visualizer will now use zlib compression for vtu files by default. If zlib is not available binary output will be used. \[Menno Fraters; 2021-06-26; [#282](github.com/GeodynamicWorldBuilder/WorldBuilder/pull/282)\] - The return argument type of the distance_point_from_curved_planes function has been converted from a map to a struct, requiring a change in the plugin interfaces for 'fault_models' and 'subducting_plate_models', but significantly increasing the speed of the function and all functions that access its returned values. \[Rene Gassmoeller; 2021-06-27; [#289](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/289)\] - The plugin systems (temperature, composition and grains) and the distance_point_from_curved_planes function now all pass a precomputed NaturalCoordinate, besides just the cartesian position. It turns out that this can make a significant performance differce. \[Issue found and solution suggested by Wolfgang Bangerth, implemented and tested by Menno Fraters; 2021-07-03; [#300](github.com/GeodynamicWorldBuilder/WorldBuilder/pull/300) and [#219](github.com/GeodynamicWorldBuilder/WorldBuilder/issues/219)\] +- Introduces a bounding box for searching the fault and the subducting plate model properties (temperature, composition, grains) in the lateral direction. This reduces the computation time by constraining the number of points for which the expensive distance_point_from_curved_planes function is called. \[Arushi Saxena; 2021-07-15; [#219](https://github.com/GeodynamicWorldBuilder/WorldBuilder/issues/219)\] ## [0.4.0] ### Added