-
Notifications
You must be signed in to change notification settings - Fork 192
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Add coordinate map SphericalToCartesianPfaffian
This map will be used in spherical shell elements in order to use a spherical harmonic basis.
- Loading branch information
Showing
5 changed files
with
292 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
133 changes: 133 additions & 0 deletions
133
src/Domain/CoordinateMaps/SphericalToCartesianPfaffian.cpp
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,133 @@ | ||
// Distributed under the MIT License. | ||
// See LICENSE.txt for details. | ||
|
||
#include "Domain/CoordinateMaps/SphericalToCartesianPfaffian.hpp" | ||
|
||
#include <cmath> | ||
#include <cstddef> | ||
#include <pup.h> | ||
|
||
#include "DataStructures/DataVector.hpp" | ||
#include "DataStructures/Tensor/Tensor.hpp" | ||
#include "Utilities/DereferenceWrapper.hpp" | ||
#include "Utilities/GenerateInstantiations.hpp" | ||
#include "Utilities/MakeWithValue.hpp" | ||
|
||
namespace domain::CoordinateMaps { | ||
SphericalToCartesianPfaffian::SphericalToCartesianPfaffian() = default; | ||
SphericalToCartesianPfaffian::SphericalToCartesianPfaffian( | ||
SphericalToCartesianPfaffian&&) = default; | ||
SphericalToCartesianPfaffian::SphericalToCartesianPfaffian( | ||
const SphericalToCartesianPfaffian&) = default; | ||
SphericalToCartesianPfaffian& SphericalToCartesianPfaffian::operator=( | ||
const SphericalToCartesianPfaffian&) = default; | ||
SphericalToCartesianPfaffian& SphericalToCartesianPfaffian::operator=( | ||
SphericalToCartesianPfaffian&&) = default; | ||
|
||
template <typename T> | ||
std::array<tt::remove_cvref_wrap_t<T>, 3> | ||
SphericalToCartesianPfaffian::operator()( | ||
const std::array<T, 3>& source_coords) const { | ||
const auto& [r, theta, phi] = source_coords; | ||
return { | ||
{r * sin(theta) * cos(phi), r * sin(theta) * sin(phi), r * cos(theta)}}; | ||
} | ||
|
||
// NOLINTNEXTLINE(readability-convert-member-functions-to-static) | ||
std::optional<std::array<double, 3>> SphericalToCartesianPfaffian::inverse( | ||
const std::array<double, 3>& target_coords) const { | ||
const auto& [x, y, z] = target_coords; | ||
if (UNLIKELY(y == 0.0 and x == 0.0)) { | ||
if (UNLIKELY(z == 0.0)) { | ||
return std::array{0.0, 0.5 * M_PI, 0.0}; | ||
} else { | ||
return std::array{std::abs(z), z > 0.0 ? 0.0 : M_PI, 0.0}; | ||
} | ||
} else { | ||
const double r = std::hypot(x, y, z); | ||
return std::array{r, acos(z / r), | ||
x > 0.0 ? atan2(y, x) : atan2(y, x) + 2.0 * M_PI}; | ||
} | ||
} | ||
|
||
template <typename T> | ||
tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> | ||
SphericalToCartesianPfaffian::jacobian( | ||
const std::array<T, 3>& source_coords) const { | ||
const auto& [r, theta, phi] = source_coords; | ||
using DataType = tt::remove_cvref_wrap_t<T>; | ||
tnsr::Ij<DataType, 3, Frame::NoFrame> jacobian_matrix{ | ||
make_with_value<DataType>(dereference_wrapper(r), 0.0)}; | ||
// Pfaffian basis means phi components are 1 / sin_theta times those of a | ||
// coord basis | ||
const auto& cos_theta = get<2, 0>(jacobian_matrix) = cos(theta); | ||
const auto& sin_theta = get<2, 1>(jacobian_matrix) = sin(theta); | ||
const auto& cos_phi = get<1, 2>(jacobian_matrix) = cos(phi); | ||
const auto& sin_phi = get<0, 2>(jacobian_matrix) = sin(phi); | ||
get<0, 0>(jacobian_matrix) = sin_theta * cos_phi; | ||
get<1, 0>(jacobian_matrix) = sin_theta * sin_phi; | ||
get<0, 1>(jacobian_matrix) = r * cos_theta * cos_phi; | ||
get<1, 1>(jacobian_matrix) = r * cos_theta * sin_phi; | ||
get<2, 1>(jacobian_matrix) *= -r; | ||
get<0, 2>(jacobian_matrix) *= -r; | ||
get<1, 2>(jacobian_matrix) *= r; | ||
// get<2, 2>(jacobian_matrix) is zero | ||
return jacobian_matrix; | ||
} | ||
|
||
template <typename T> | ||
tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> | ||
SphericalToCartesianPfaffian::inv_jacobian( | ||
const std::array<T, 3>& source_coords) const { | ||
const auto& [r, theta, phi] = source_coords; | ||
using DataType = tt::remove_cvref_wrap_t<T>; | ||
tnsr::Ij<DataType, 3, Frame::NoFrame> inv_jacobian_matrix{ | ||
make_with_value<DataType>(dereference_wrapper(r), 0.0)}; | ||
// Pfaffian basis means phi components are sin_theta times those of a coord | ||
// basis | ||
const auto& cos_theta = get<0, 2>(inv_jacobian_matrix) = cos(theta); | ||
const auto& sin_theta = get<1, 2>(inv_jacobian_matrix) = sin(theta); | ||
const auto& cos_phi = get<2, 1>(inv_jacobian_matrix) = cos(phi); | ||
const auto& sin_phi = get<2, 0>(inv_jacobian_matrix) = sin(phi); | ||
const auto& one_over_r = get<2, 2>(inv_jacobian_matrix) = 1.0/r; | ||
get<0, 0>(inv_jacobian_matrix) = sin_theta * cos_phi; | ||
get<0, 1>(inv_jacobian_matrix) = sin_theta * sin_phi; | ||
get<1, 0>(inv_jacobian_matrix) = cos_theta * cos_phi * one_over_r; | ||
get<1, 1>(inv_jacobian_matrix) = cos_theta * sin_phi * one_over_r; | ||
get<1, 2>(inv_jacobian_matrix) *= -one_over_r; | ||
get<2, 0>(inv_jacobian_matrix) *= -one_over_r; | ||
get<2, 1>(inv_jacobian_matrix) *= one_over_r; | ||
get<2, 2>(inv_jacobian_matrix) *= 0.0; | ||
return inv_jacobian_matrix; | ||
} | ||
|
||
void SphericalToCartesianPfaffian::pup(PUP::er& /*p*/) {} | ||
|
||
bool operator==(const SphericalToCartesianPfaffian& /*lhs*/, | ||
const SphericalToCartesianPfaffian& /*rhs*/) { | ||
return true; | ||
} | ||
|
||
bool operator!=(const SphericalToCartesianPfaffian& lhs, | ||
const SphericalToCartesianPfaffian& rhs) { | ||
return not(lhs == rhs); | ||
} | ||
|
||
#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) | ||
|
||
#define INSTANTIATE_DTYPE(_, data) \ | ||
template std::array<tt::remove_cvref_wrap_t<DTYPE(data)>, 3> \ | ||
SphericalToCartesianPfaffian::operator()( \ | ||
const std::array<DTYPE(data), 3>& source_coords) const; \ | ||
template tnsr::Ij<tt::remove_cvref_wrap_t<DTYPE(data)>, 3, Frame::NoFrame> \ | ||
SphericalToCartesianPfaffian::jacobian( \ | ||
const std::array<DTYPE(data), 3>& source_coords) const; \ | ||
template tnsr::Ij<tt::remove_cvref_wrap_t<DTYPE(data)>, 3, Frame::NoFrame> \ | ||
SphericalToCartesianPfaffian::inv_jacobian( \ | ||
const std::array<DTYPE(data), 3>& source_coords) const; | ||
|
||
GENERATE_INSTANTIATIONS(INSTANTIATE_DTYPE, | ||
(double, DataVector, | ||
std::reference_wrapper<const double>, | ||
std::reference_wrapper<const DataVector>)) | ||
} // namespace domain::CoordinateMaps |
87 changes: 87 additions & 0 deletions
87
src/Domain/CoordinateMaps/SphericalToCartesianPfaffian.hpp
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,87 @@ | ||
// Distributed under the MIT License. | ||
// See LICENSE.txt for details. | ||
|
||
#pragma once | ||
|
||
#include <array> | ||
#include <cstddef> | ||
#include <optional> | ||
|
||
#include "DataStructures/Tensor/TypeAliases.hpp" | ||
#include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp" | ||
|
||
/// \cond | ||
namespace PUP { | ||
class er; | ||
} // namespace PUP | ||
/// \endcond | ||
|
||
namespace domain::CoordinateMaps { | ||
|
||
/*! | ||
* \ingroup CoordinateMapsGroup | ||
* | ||
* \brief Pfaffian transformation from spherical to Cartesian coordinates. | ||
* | ||
* \note This map is designed to be used together with Spherepack! Spherepack | ||
* expects a Pfaffian transformation of the derivatives. | ||
* | ||
* \details This is a Pfaffian mapping from \f$(r,\theta,\phi) \rightarrow | ||
* (x,y,z) \f$. | ||
* | ||
* The formula for the mapping is... | ||
* \f{eqnarray*} | ||
* x &=& r \sin\theta \cos\phi \\ | ||
* y &=& r \sin\theta \sin\phi \\ | ||
* z &=& r \cos\theta | ||
* \f} | ||
* | ||
* The Pfaffian basis vectors | ||
* \f$ (e_{\hat r}, e_{\hat \theta}, e_{\hat \phi})\f$ | ||
* are related to the coordinate basis vectors | ||
* \f$ (e_r, e_{\theta}, e_{\phi})\f$ | ||
* by... | ||
* \f{eqnarray*} | ||
* e_{\hat r} &=& e_r \\ | ||
* e_{\hat \theta} &=& e_{\theta} \\ | ||
* e_{\hat \phi} &=& \frac{1}{\sin \theta} e_{\phi} | ||
* \f} | ||
*/ | ||
class SphericalToCartesianPfaffian { | ||
public: | ||
static constexpr size_t dim = 3; | ||
SphericalToCartesianPfaffian(); | ||
~SphericalToCartesianPfaffian() = default; | ||
SphericalToCartesianPfaffian(SphericalToCartesianPfaffian&&); | ||
SphericalToCartesianPfaffian(const SphericalToCartesianPfaffian&); | ||
SphericalToCartesianPfaffian& operator=(const SphericalToCartesianPfaffian&); | ||
SphericalToCartesianPfaffian& operator=(SphericalToCartesianPfaffian&&); | ||
|
||
template <typename T> | ||
std::array<tt::remove_cvref_wrap_t<T>, 3> operator()( | ||
const std::array<T, 3>& source_coords) const; | ||
|
||
// NOLINTNEXTLINE(readability-convert-member-functions-to-static) | ||
std::optional<std::array<double, 3>> inverse( | ||
const std::array<double, 3>& target_coords) const; | ||
|
||
template <typename T> | ||
tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> jacobian( | ||
const std::array<T, 3>& source_coords) const; | ||
|
||
template <typename T> | ||
tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> inv_jacobian( | ||
const std::array<T, 3>& source_coords) const; | ||
|
||
// NOLINTNEXTLINE(google-runtime-references) | ||
void pup(PUP::er& p); | ||
|
||
static constexpr bool is_identity() { return false; } | ||
}; | ||
|
||
bool operator==(const SphericalToCartesianPfaffian& lhs, | ||
const SphericalToCartesianPfaffian& rhs); | ||
|
||
bool operator!=(const SphericalToCartesianPfaffian& lhs, | ||
const SphericalToCartesianPfaffian& rhs); | ||
} // namespace domain::CoordinateMaps |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
69 changes: 69 additions & 0 deletions
69
tests/Unit/Domain/CoordinateMaps/Test_SphericalToCartesianPfaffian.cpp
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,69 @@ | ||
// Distributed under the MIT License. | ||
// See LICENSE.txt for details. | ||
|
||
#include "Framework/TestingFramework.hpp" | ||
|
||
#include <array> | ||
#include <optional> | ||
#include <pup.h> | ||
|
||
#include "DataStructures/DataVector.hpp" | ||
#include "DataStructures/Tensor/Tensor.hpp" | ||
#include "Domain/CoordinateMaps/SphericalToCartesianPfaffian.hpp" | ||
#include "Framework/TestHelpers.hpp" | ||
#include "Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp" | ||
#include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" | ||
#include "NumericalAlgorithms/Spectral/Mesh.hpp" | ||
|
||
namespace domain { | ||
namespace { | ||
void test_map_at_point(const CoordinateMaps::SphericalToCartesianPfaffian& map, | ||
const std::array<double, 3>& source_point, | ||
const std::array<double, 3>& target_point) { | ||
test_inverse_map(map, source_point); | ||
if (source_point != | ||
std::array{0.0, 0.5 * M_PI, 0.0}) { // inv jac singular at origin | ||
test_coordinate_map_argument_types(map, source_point); | ||
test_inv_jacobian(map, source_point); | ||
} | ||
CHECK_ITERABLE_APPROX(map(source_point), target_point); | ||
CHECK_ITERABLE_APPROX(map.inverse(target_point).value(), source_point); | ||
} | ||
|
||
void test_map(const CoordinateMaps::SphericalToCartesianPfaffian& map) { | ||
CHECK(not map.is_identity()); | ||
CHECK_FALSE(map != map); | ||
test_serialization(map); | ||
const std::array<double, 3> source_origin{{0.0, 0.5 * M_PI, 0.0}}; | ||
const std::array<double, 3> target_origin{{0.0, 0.0, 0.0}}; | ||
test_map_at_point(map, source_origin, target_origin); | ||
const std::array<double, 3> source_north_pole{{0.5, 0.0, 0.0}}; | ||
const std::array<double, 3> target_north_pole{{0.0, 0.0, 0.5}}; | ||
test_map_at_point(map, source_north_pole, target_north_pole); | ||
const std::array<double, 3> source_south_pole{{0.25, M_PI, 0.0}}; | ||
const std::array<double, 3> target_south_pole{{0.0, 0.0, -0.25}}; | ||
test_map_at_point(map, source_south_pole, target_south_pole); | ||
const Mesh<3> mesh{ | ||
{5, 3, 5}, | ||
{Spectral::Basis::Legendre, Spectral::Basis::SphericalHarmonic, | ||
Spectral::Basis::SphericalHarmonic}, | ||
{Spectral::Quadrature::GaussLobatto, Spectral::Quadrature::Gauss, | ||
Spectral::Quadrature::Equiangular}}; | ||
const auto xi = logical_coordinates(mesh); | ||
const std::array<DataVector, 3> source_coords{xi[0] + 2.0, xi[1], xi[2]}; | ||
test_inv_jacobian(map, source_coords); | ||
} | ||
|
||
void test() { | ||
const CoordinateMaps::SphericalToCartesianPfaffian original_map{}; | ||
test_map(original_map); | ||
const auto serialized_map = serialize_and_deserialize(original_map); | ||
test_map(serialized_map); | ||
} | ||
} // namespace | ||
|
||
SPECTRE_TEST_CASE("Unit.Domain.CoordinateMaps.SphericalToCartesianPfaffian", | ||
"[Domain][Unit]") { | ||
test(); | ||
} | ||
} // namespace domain |