Skip to content

Commit

Permalink
Merge c5934bf into c3446f1
Browse files Browse the repository at this point in the history
  • Loading branch information
MFraters committed Nov 20, 2021
2 parents c3446f1 + c5934bf commit 66949d0
Show file tree
Hide file tree
Showing 7 changed files with 268 additions and 15 deletions.
111 changes: 111 additions & 0 deletions include/world_builder/coordinate_systems/invalid.h
@@ -0,0 +1,111 @@
/*
Copyright (C) 2018 - 2021 by the authors of the World Builder code.
This file is part of the World Builder.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <https://www.gnu.org/licenses/>.
*/

#ifndef WORLD_BUILDER_COORDINATE_SYSTEMS_INVALID_H
#define WORLD_BUILDER_COORDINATE_SYSTEMS_INVALID_H

#include "world_builder/coordinate_systems/interface.h"

#include "world_builder/utilities.h"


namespace WorldBuilder
{

namespace CoordinateSystems
{
/**
* This implements a Invalid geometry model.The Invalid geometry model
* doesn't do anything with the coordinates, but is needed to have a common
* interface for all the geometry models.
*/
class Invalid final : public Interface
{
public:
/**
* constructor
*/
Invalid(WorldBuilder::World *world);

/**
* Destructor
*/
~Invalid() override final;

/**
* declare and read in the world builder file into the parameters class
*/
static
void declare_entries(Parameters &prm, const std::string &parent_name = "");

/**
* declare and read in the world builder file into the parameters class
*/
void parse_entries(Parameters &prm) override final;


/**
* Returns what the natural coordinate system for this Coordinate System is.
*/
CoordinateSystem natural_coordinate_system() const override final;

/**
* Returns what method should be used to go down with an angle into
* the domain.
* \sa DepthMethod
*/
DepthMethod depth_method() const override final;

/**
* Takes the Invalid points (x,z or x,y,z) and returns standardized
* coordinates which are most 'natural' to the geometry model. For a box
* this will be (x,z) in 2d or (x,y,z) in 3d, and for a spheroid geometry
* model it will be (radius, longitude) in 2d and (radius, longitude,
* latitude) in 3d.
*/
std::array<double,3> cartesian_to_natural_coordinates(const std::array<double,3> &position) const override final;

/**
* Undoes the action of invalid_to_natural_coordinates, and turns the
* coordinate system which is most 'natural' to the geometry model into
* Invalid coordinates.
*/
std::array<double,3> natural_to_cartesian_coordinates(const std::array<double,3> &position) const override final;


/**
* Computes the distance between two points which are on the same depth.
* The input is two 3d points at that depth. It is implemented as the
*/
double distance_between_points_at_same_depth(const Point<3> &point_1, const Point<3> &point_2) const override final;

/**
* Returns the max model depth. This always returns infinity for Invalid
* models.
*/
virtual
double max_model_depth() const override final;

private:

};
} // namespace CoordinateSystems
} // namespace WorldBuilder

#endif
6 changes: 6 additions & 0 deletions include/world_builder/utilities.h
Expand Up @@ -104,6 +104,12 @@ namespace WorldBuilder
*/
double get_depth_coordinate() const;

/**
* Return a reference to the coordinate that represents the 'depth' direction
* in the chosen coordinate system.
*/
double &get_depth_coordinate();

/**
* get the coordinate system type of this coordinate.
*/
Expand Down
88 changes: 88 additions & 0 deletions source/coordinate_systems/invalid.cc
@@ -0,0 +1,88 @@
/*
Copyright (C) 2018-2021 by the authors of the World Builder code.
This file is part of the World Builder.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <https://www.gnu.org/licenses/>.
*/

#include "world_builder/coordinate_systems/invalid.h"
#include "world_builder/types/object.h"


namespace WorldBuilder
{
namespace CoordinateSystems
{
Invalid::Invalid(WorldBuilder::World *world_)
{
this->world = world_;
}

Invalid::~Invalid()
= default;

void
Invalid::parse_entries(Parameters & /*prm*/)
{}


CoordinateSystem
Invalid::natural_coordinate_system() const
{
return CoordinateSystem::invalid;
}


DepthMethod
Invalid::depth_method() const
{
return DepthMethod::none;
}


std::array<double,3>
Invalid::cartesian_to_natural_coordinates(const std::array<double,3> &position) const
{
return {NaN::DSNAN,NaN::DSNAN,NaN::DSNAN};
}


std::array<double,3>
Invalid::natural_to_cartesian_coordinates(const std::array<double,3> &position) const
{
return {NaN::DSNAN,NaN::DSNAN,NaN::DSNAN};
}


double
Invalid::distance_between_points_at_same_depth(const Point<3> &point_1, const Point<3> &point_2) const
{
return NaN::DSNAN;
}


double
Invalid::max_model_depth() const
{
return NaN::DSNAN;
}

/**
* Not registering plugin because this is only used for testing.
*/
//WB_REGISTER_COORDINATE_SYSTEM(Invalid, invalid)
} // namespace CoordinateSystems
} // namespace WorldBuilder

Expand Up @@ -120,8 +120,9 @@ namespace WorldBuilder
{
if (depth <= max_depth && depth >= min_depth)
{
WorldBuilder::Utilities::NaturalCoordinate position_in_natural_coordinates = WorldBuilder::Utilities::NaturalCoordinate(position,
*(world->parameters.coordinate_system));
WorldBuilder::Utilities::NaturalCoordinate position_in_natural_coordinates_at_min_depth = WorldBuilder::Utilities::NaturalCoordinate(position,
*(world->parameters.coordinate_system));
position_in_natural_coordinates_at_min_depth.get_depth_coordinate() += depth-min_depth;


double bottom_temperature_local = bottom_temperature;
Expand All @@ -140,7 +141,8 @@ namespace WorldBuilder

// first find if the coordinate is on this side of a ridge
unsigned int relevant_ridge = 0;
const Point<2> check_point(position_in_natural_coordinates.get_surface_coordinates(),position_in_natural_coordinates.get_coordinate_system());
const Point<2> check_point(position_in_natural_coordinates_at_min_depth.get_surface_coordinates(),
position_in_natural_coordinates_at_min_depth.get_coordinate_system());

// if there is only one ridge, there is no transform
if (mid_oceanic_ridges.size() > 1)
Expand Down Expand Up @@ -198,11 +200,14 @@ namespace WorldBuilder

Point<3> compare_point(coordinate_system);

compare_point[0] = coordinate_system == cartesian ? Pb[0] : position_in_natural_coordinates.get_depth_coordinate();
compare_point[0] = coordinate_system == cartesian ? Pb[0] : position_in_natural_coordinates_at_min_depth.get_depth_coordinate();
compare_point[1] = coordinate_system == cartesian ? Pb[1] : Pb[0];
compare_point[2] = coordinate_system == cartesian ? position_in_natural_coordinates.get_depth_coordinate() : Pb[1];
compare_point[2] = coordinate_system == cartesian ? position_in_natural_coordinates_at_min_depth.get_depth_coordinate() : Pb[1];

distance_ridge = std::min(distance_ridge,this->world->parameters.coordinate_system->distance_between_points_at_same_depth(Point<3>(position_in_natural_coordinates.get_coordinates(),position_in_natural_coordinates.get_coordinate_system()),compare_point));
distance_ridge = std::min(distance_ridge,
this->world->parameters.coordinate_system->distance_between_points_at_same_depth(Point<3>(position_in_natural_coordinates_at_min_depth.get_coordinates(),
position_in_natural_coordinates_at_min_depth.get_coordinate_system()),
compare_point));

}

Expand Down
14 changes: 7 additions & 7 deletions source/features/oceanic_plate_models/temperature/plate_model.cc
Expand Up @@ -117,9 +117,9 @@ namespace WorldBuilder
{
if (depth <= max_depth && depth >= min_depth)
{
WorldBuilder::Utilities::NaturalCoordinate position_in_natural_coordinates = WorldBuilder::Utilities::NaturalCoordinate(position,
*(world->parameters.coordinate_system));

WorldBuilder::Utilities::NaturalCoordinate position_in_natural_coordinates_at_min_depth = WorldBuilder::Utilities::NaturalCoordinate(position,
*(world->parameters.coordinate_system));
position_in_natural_coordinates_at_min_depth.get_depth_coordinate() += depth-min_depth;

double bottom_temperature_local = bottom_temperature;

Expand All @@ -139,7 +139,7 @@ namespace WorldBuilder

// first find if the coordinate is on this side of a ridge
unsigned int relevant_ridge = 0;
const Point<2> check_point(position_in_natural_coordinates.get_surface_coordinates(),position_in_natural_coordinates.get_coordinate_system());
const Point<2> check_point(position_in_natural_coordinates_at_min_depth.get_surface_coordinates(),position_in_natural_coordinates_at_min_depth.get_coordinate_system());

// if there is only one ridge, there is no transform
if (mid_oceanic_ridges.size() > 1)
Expand Down Expand Up @@ -197,11 +197,11 @@ namespace WorldBuilder

Point<3> compare_point(coordinate_system);

compare_point[0] = coordinate_system == cartesian ? Pb[0] : position_in_natural_coordinates.get_depth_coordinate();
compare_point[0] = coordinate_system == cartesian ? Pb[0] : position_in_natural_coordinates_at_min_depth.get_depth_coordinate();
compare_point[1] = coordinate_system == cartesian ? Pb[1] : Pb[0];
compare_point[2] = coordinate_system == cartesian ? position_in_natural_coordinates.get_depth_coordinate() : Pb[1];
compare_point[2] = coordinate_system == cartesian ? position_in_natural_coordinates_at_min_depth.get_depth_coordinate() : Pb[1];

distance_ridge = std::min(distance_ridge,this->world->parameters.coordinate_system->distance_between_points_at_same_depth(Point<3>(position_in_natural_coordinates.get_coordinates(),position_in_natural_coordinates.get_coordinate_system()),compare_point));
distance_ridge = std::min(distance_ridge,this->world->parameters.coordinate_system->distance_between_points_at_same_depth(Point<3>(position_in_natural_coordinates_at_min_depth.get_coordinates(),position_in_natural_coordinates_at_min_depth.get_coordinate_system()),compare_point));

}

Expand Down
21 changes: 19 additions & 2 deletions source/utilities.cc
Expand Up @@ -270,7 +270,7 @@ namespace WorldBuilder
break;

default:
WBAssert (false, "Coordinate system not implemented.");
WBAssertThrow (false, "Coordinate system not implemented.");
}

return coordinate;
Expand All @@ -295,12 +295,29 @@ namespace WorldBuilder
return coordinates[0];

default:
WBAssert (false, "Coordinate system not implemented.");
WBAssertThrow (false, "Coordinate system not implemented.");
}

return 0;
}

double &NaturalCoordinate::get_depth_coordinate()
{
switch (coordinate_system)
{
case CoordinateSystem::cartesian:
return coordinates[2];

case CoordinateSystem::spherical:
return coordinates[0];

default:
WBAssertThrow (false, "Coordinate system not implemented.");
}

return coordinates[2];
}


std::array<double,3>
cartesian_to_spherical_coordinates(const Point<3> &position)
Expand Down
26 changes: 26 additions & 0 deletions tests/unit_tests/unit_test_world_builder.cc
Expand Up @@ -24,6 +24,7 @@
#include "world_builder/config.h"
#include "world_builder/coordinate_system.h"
#include "world_builder/coordinate_systems/cartesian.h"
#include "world_builder/coordinate_systems/invalid.h"
#include "world_builder/coordinate_systems/interface.h"
#include "world_builder/features/continental_plate.h"
#include "world_builder/features/interface.h"
Expand Down Expand Up @@ -613,6 +614,31 @@ TEST_CASE("WorldBuilder Utilities: Natural Coordinate")
CHECK(nsp1_surface_array[1] == Approx(0.9302740141));
CHECK(nsp1.get_depth_coordinate() == Approx(std::sqrt(1.0 * 1.0 + 2.0 * 2.0 + 3.0 * 3.0)));

// Invalid tests
std::string file_name = WorldBuilder::Data::WORLD_BUILDER_SOURCE_DIR + "/tests/data/subducting_plate_different_angles_cartesian.wb";
WorldBuilder::World world(file_name);
Parameters prm(world);
std::unique_ptr<CoordinateSystems::Interface> invalid(new CoordinateSystems::Invalid(nullptr));
invalid->parse_entries(prm);
Utilities::NaturalCoordinate ivp1(Point<3>(1,2,3,CoordinateSystem::invalid),*invalid);
std::array<double,3> ivp1_array = ivp1.get_coordinates();
CHECK(std::isnan(ivp1_array[0]));
CHECK(std::isnan(ivp1_array[1]));
CHECK(std::isnan(ivp1_array[2]));
CHECK_THROWS(ivp1.get_surface_coordinates());
CHECK_THROWS(ivp1.get_depth_coordinate());
CHECK(std::isnan(invalid->distance_between_points_at_same_depth(Point<3>(1,2,3,CoordinateSystem::invalid),
Point<3>(1,2,3,CoordinateSystem::invalid))));
CHECK(invalid->depth_method() == DepthMethod::none);

std::array<double,3> iv_array = invalid->natural_to_cartesian_coordinates({1,2,3});
CHECK(std::isnan(iv_array[0]));
CHECK(std::isnan(iv_array[1]));
CHECK(std::isnan(iv_array[2]));

CHECK(std::isnan(invalid->max_model_depth()));


}

TEST_CASE("WorldBuilder Utilities: Coordinate systems transformations")
Expand Down

0 comments on commit 66949d0

Please sign in to comment.