Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Changes to include a bounding box in computing fault composition. #264

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
20 changes: 20 additions & 0 deletions include/world_builder/features/fault.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -179,6 +192,13 @@ namespace WorldBuilder
double maximum_total_fault_length;
double maximum_fault_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_fault_cartesian;

};
} // namespace Features
Expand Down
19 changes: 19 additions & 0 deletions include/world_builder/features/subducting_plate.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand Down
67 changes: 64 additions & 3 deletions source/features/fault.cc
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
#include "world_builder/types/point.h"
#include "world_builder/types/unsigned_int.h"
#include "world_builder/world.h"
#include <algorithm>

using namespace std;

Expand Down Expand Up @@ -325,6 +326,59 @@ namespace WorldBuilder
total_fault_length[i] = local_total_fault_length;
maximum_total_fault_length = std::max(maximum_total_fault_length, local_total_fault_length);
}


std::vector<double> x_list(original_number_of_coordinates,0.0);
std::vector<double> y_list(original_number_of_coordinates,0.0);

for (size_t j=0; j<original_number_of_coordinates; ++j)
{
x_list[j] = coordinates[j][0];
y_list[j] = coordinates[j][1];
}


min_along_x = *std::min_element(x_list.begin(), x_list.end());
max_along_x = *std::max_element(x_list.begin(), x_list.end());
min_along_y = *std::min_element(y_list.begin(), y_list.end());
max_along_y = *std::max_element(y_list.begin(), y_list.end());

min_lat_cos_inv = 1 / std::cos(min_along_y);
max_lat_cos_inv = 1 / std::cos(max_along_y);

buffer_around_fault_cartesian = (maximum_fault_thickness + maximum_total_fault_length);

}


BoundingBox<2>
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>, 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>, 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;
}


Expand All @@ -348,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
Expand Down Expand Up @@ -468,6 +524,7 @@ namespace WorldBuilder
return temperature;
}


double
Fault::composition(const Point<3> &position,
const NaturalCoordinate &position_in_natural_coordinates,
Expand All @@ -479,7 +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)
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
Expand Down Expand Up @@ -613,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
Expand Down
69 changes: 66 additions & 3 deletions source/features/subducting_plate.cc
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
#include "world_builder/types/point.h"
#include "world_builder/types/unsigned_int.h"
#include "world_builder/world.h"
#include <algorithm>

using namespace std;

Expand Down Expand Up @@ -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<double> x_list(original_number_of_coordinates,0.0);
std::vector<double> y_list(original_number_of_coordinates,0.0);

for (size_t j=0; j<original_number_of_coordinates; ++j)
{
x_list[j] = coordinates[j][0];
y_list[j] = coordinates[j][1];
}

// Here, we compute the spherical bounding box using the two extreme points of the box containing 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.

min_along_x = *std::min_element(x_list.begin(), x_list.end());
max_along_x = *std::max_element(x_list.begin(), x_list.end());
min_along_y = *std::min_element(y_list.begin(), y_list.end());
max_along_y = *std::max_element(y_list.begin(), y_list.end());

min_lat_cos_inv = 1 / std::cos(min_along_y);
max_lat_cos_inv = 1 / std::cos(max_along_y);

buffer_around_slab_cartesian = (maximum_slab_thickness + maximum_total_slab_length);
}

BoundingBox<2>
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>, 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>, 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;
}


Expand All @@ -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()
Expand Down Expand Up @@ -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 =
Expand Down Expand Up @@ -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 =
Expand Down