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

new surface bounding box optimalizations. #408

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
13 changes: 11 additions & 2 deletions include/world_builder/features/fault.h
Original file line number Diff line number Diff line change
Expand Up @@ -113,8 +113,7 @@ namespace WorldBuilder
* For the spherical system, the buffer zone along the longitudal direction is calculated using the
* corresponding latitude points.
*/
BoundingBox<2> get_bounding_box (const Objects::NaturalCoordinate &position_in_natural_coordinates,
const double depth) const;
const BoundingBox<2> &get_surface_bounding_box () const;

/**
* Returns a grains (rotation matrix and grain size) based on the
Expand Down Expand Up @@ -178,6 +177,16 @@ namespace WorldBuilder
*/
double maximum_depth;

/**
* Computee 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 fault dips.
*/
Expand Down
15 changes: 12 additions & 3 deletions include/world_builder/features/subducting_plate.h
Original file line number Diff line number Diff line change
Expand Up @@ -85,15 +85,14 @@ namespace WorldBuilder


/**
* Computes the bounding points for a BoundingBox object using two extreme points in all the surface
* 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.
*/
BoundingBox<2> get_bounding_box (const Objects::NaturalCoordinate &position_in_natural_coordinates,
const double depth) const;
const BoundingBox<2> &get_surface_bounding_box () const;


/**
Expand Down Expand Up @@ -181,6 +180,16 @@ namespace WorldBuilder
*/
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.
*/
Expand Down
30 changes: 15 additions & 15 deletions source/features/fault.cc
Original file line number Diff line number Diff line change
Expand Up @@ -359,18 +359,12 @@ namespace WorldBuilder
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 Objects::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);
// Compute the surface bounding box
buffer_around_fault_cartesian = (maximum_fault_thickness + maximum_total_fault_length);
if (world->parameters.coordinate_system->natural_coordinate_system() == CoordinateSystem::spherical)
{
const double starting_radius_inv = 1 / (world->parameters.coordinate_system->max_model_depth());
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;
Expand All @@ -390,6 +384,12 @@ namespace WorldBuilder
bounding_box.second = {max_along_x, max_along_y, cartesian};
surface_bounding_box.extend(buffer_around_fault_cartesian);
}
}


const BoundingBox<2> &
Fault::get_surface_bounding_box () const
{
return surface_bounding_box;
}

Expand All @@ -415,8 +415,8 @@ namespace WorldBuilder

// todo: explain and check -starting_depth
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())))
get_surface_bounding_box().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 @@ -550,8 +550,8 @@ namespace WorldBuilder

// todo: explain and check -starting_depth
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())))
get_surface_bounding_box().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 @@ -687,8 +687,8 @@ namespace WorldBuilder

// todo: explain and check -starting_depth
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())))
get_surface_bounding_box().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
43 changes: 22 additions & 21 deletions source/features/subducting_plate.cc
Original file line number Diff line number Diff line change
Expand Up @@ -340,7 +340,7 @@ namespace WorldBuilder
}

// 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
// coordinates and an additional buffer zone that accounts for the slab 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
Expand Down Expand Up @@ -370,28 +370,22 @@ namespace WorldBuilder
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 Objects::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);
// Compute the surface bounding box
buffer_around_slab_cartesian = (maximum_slab_thickness + maximum_total_slab_length);
if (world->parameters.coordinate_system->natural_coordinate_system() == CoordinateSystem::spherical)
{
const double starting_radius_inv = 1 / (world->parameters.coordinate_system->max_model_depth());
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;
const double buffer_around_slab_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.first = {(min_along_x - buffer_around_slab_spherical * min_lat_cos_inv) ,
(min_along_y - buffer_around_slab_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
spherical_bounding_box.second = {(max_along_x + buffer_around_slab_spherical * max_lat_cos_inv) ,
(max_along_y + buffer_around_slab_spherical), spherical
};
}
else if (world->parameters.coordinate_system->natural_coordinate_system() == CoordinateSystem::cartesian)
Expand All @@ -401,6 +395,13 @@ namespace WorldBuilder
bounding_box.second = {max_along_x, max_along_y, cartesian};
surface_bounding_box.extend(buffer_around_slab_cartesian);
}
}



const BoundingBox<2> &
SubductingPlate::get_surface_bounding_box () const
{
return surface_bounding_box;
}

Expand All @@ -426,8 +427,8 @@ namespace WorldBuilder

// todo: explain and check -starting_depth
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())))
get_surface_bounding_box().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 @@ -565,8 +566,8 @@ namespace WorldBuilder

// todo: explain and check -starting_depth
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())))
get_surface_bounding_box().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 @@ -698,8 +699,8 @@ namespace WorldBuilder

// todo: explain and check -starting_depth
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())))
get_surface_bounding_box().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