Merge pull request #410 from gsb76/gab_add_equiangular_rectangle

Add EquiangularMap

src/Domain/CoordinateMaps/CMakeLists.txt

@@ -5,6 +5,7 @@ set(LIBRARY CoordinateMaps)
 
 set(LIBRARY_SOURCES
     AffineMap.cpp
+    Equiangular.cpp
     Identity.cpp
     Rotation.cpp
     Wedge2D.cpp

src/Domain/CoordinateMaps/Equiangular.cpp

@@ -0,0 +1,149 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Domain/CoordinateMaps/Equiangular.hpp"
+
+#include <pup.h>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Utilities/DereferenceWrapper.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+
+namespace CoordinateMaps {
+
+Equiangular::Equiangular(const double A, const double B, const double a,
+                         const double b) noexcept
+    : A_(A),
+      B_(B),
+      a_(a),
+      b_(b),
+      length_of_domain_over_m_pi_4_((B - A) / M_PI_4),
+      length_of_range_(b - a),
+      m_pi_4_over_length_of_domain_(1.0 / length_of_domain_over_m_pi_4_),
+      one_over_length_of_range_(1.0 / length_of_range_),
+      linear_jacobian_times_m_pi_4_(length_of_range_ /
+                                    length_of_domain_over_m_pi_4_),
+      linear_inverse_jacobian_over_m_pi_4_(length_of_domain_over_m_pi_4_ /
+                                           length_of_range_) {}
+
+template <typename T>
+std::array<std::decay_t<tt::remove_reference_wrapper_t<T>>, 1> Equiangular::
+operator()(const std::array<T, 1>& xi) const noexcept {
+  return {{0.5 * (a_ + b_ +
+                  length_of_range_ * tan(m_pi_4_over_length_of_domain_ *
+                                         (-B_ - A_ + 2.0 * xi[0])))}};
+}
+
+template <typename T>
+std::array<std::decay_t<tt::remove_reference_wrapper_t<T>>, 1>
+Equiangular::inverse(const std::array<T, 1>& x) const noexcept {
+  return {
+      {0.5 * (A_ + B_ +
+              length_of_domain_over_m_pi_4_ *
+                  atan(one_over_length_of_range_ * (-a_ - b_ + 2.0 * x[0])))}};
+}
+
+template <typename T>
+Tensor<std::decay_t<tt::remove_reference_wrapper_t<T>>,
+       tmpl::integral_list<std::int32_t, 2, 1>,
+       index_list<SpatialIndex<1, UpLo::Up, Frame::NoFrame>,
+                  SpatialIndex<1, UpLo::Lo, Frame::NoFrame>>>
+Equiangular::jacobian(const std::array<T, 1>& xi) const noexcept {
+  const std::decay_t<tt::remove_reference_wrapper_t<T>> tan_variable =
+      tan(m_pi_4_over_length_of_domain_ * (-B_ - A_ + 2.0 * xi[0]));
+  return Tensor<std::decay_t<tt::remove_reference_wrapper_t<T>>,
+                tmpl::integral_list<std::int32_t, 2, 1>,
+                index_list<SpatialIndex<1, UpLo::Up, Frame::NoFrame>,
+                           SpatialIndex<1, UpLo::Lo, Frame::NoFrame>>>{
+      linear_jacobian_times_m_pi_4_ * (1.0 + square(tan_variable))};
+}
+
+template <typename T>
+Tensor<std::decay_t<tt::remove_reference_wrapper_t<T>>,
+       tmpl::integral_list<std::int32_t, 2, 1>,
+       index_list<SpatialIndex<1, UpLo::Up, Frame::NoFrame>,
+                  SpatialIndex<1, UpLo::Lo, Frame::NoFrame>>>
+Equiangular::inv_jacobian(const std::array<T, 1>& xi) const noexcept {
+  const std::decay_t<tt::remove_reference_wrapper_t<T>> tan_variable =
+      tan(m_pi_4_over_length_of_domain_ * (-B_ - A_ + 2.0 * xi[0]));
+  return Tensor<std::decay_t<tt::remove_reference_wrapper_t<T>>,
+                tmpl::integral_list<std::int32_t, 2, 1>,
+                index_list<SpatialIndex<1, UpLo::Up, Frame::NoFrame>,
+                           SpatialIndex<1, UpLo::Lo, Frame::NoFrame>>>{
+      linear_inverse_jacobian_over_m_pi_4_ / (1.0 + square(tan_variable))};
+}
+
+void Equiangular::pup(PUP::er& p) noexcept {
+  p | A_;
+  p | B_;
+  p | a_;
+  p | b_;
+  p | length_of_domain_over_m_pi_4_;
+  p | length_of_range_;
+  p | m_pi_4_over_length_of_domain_;
+  p | one_over_length_of_range_;
+  p | linear_jacobian_times_m_pi_4_;
+  p | linear_inverse_jacobian_over_m_pi_4_;
+}
+
+bool operator==(const CoordinateMaps::Equiangular& lhs,
+                const CoordinateMaps::Equiangular& rhs) noexcept {
+  return lhs.A_ == rhs.A_ and lhs.B_ == rhs.B_ and lhs.a_ == rhs.a_ and
+         lhs.b_ == rhs.b_ and
+         lhs.length_of_domain_over_m_pi_4_ ==
+             rhs.length_of_domain_over_m_pi_4_ and
+         lhs.length_of_range_ == rhs.length_of_range_ and
+         lhs.m_pi_4_over_length_of_domain_ ==
+             rhs.m_pi_4_over_length_of_domain_ and
+         lhs.one_over_length_of_range_ == rhs.one_over_length_of_range_ and
+         lhs.linear_jacobian_times_m_pi_4_ ==
+             rhs.linear_jacobian_times_m_pi_4_ and
+         lhs.linear_inverse_jacobian_over_m_pi_4_ ==
+             rhs.linear_inverse_jacobian_over_m_pi_4_;
+}
+
+// Explicit instantiations
+/// \cond
+#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
+
+#define INSTANTIATE(_, data)                                                  \
+  template std::array<DTYPE(data), 1> Equiangular::operator()(                \
+      const std::array<std::reference_wrapper<const DTYPE(data)>, 1>& /*xi*/) \
+      const noexcept;                                                         \
+  template std::array<DTYPE(data), 1> Equiangular::operator()(                \
+      const std::array<DTYPE(data), 1>& /*xi*/) const noexcept;               \
+  template std::array<DTYPE(data), 1> Equiangular::inverse(                   \
+      const std::array<std::reference_wrapper<const DTYPE(data)>, 1>& /*xi*/) \
+      const noexcept;                                                         \
+  template std::array<DTYPE(data), 1> Equiangular::inverse(                   \
+      const std::array<DTYPE(data), 1>& /*xi*/) const noexcept;               \
+  template Tensor<DTYPE(data), tmpl::integral_list<std::int32_t, 2, 1>,       \
+                  index_list<SpatialIndex<1, UpLo::Up, Frame::NoFrame>,       \
+                             SpatialIndex<1, UpLo::Lo, Frame::NoFrame>>>      \
+  Equiangular::jacobian(                                                      \
+      const std::array<std::reference_wrapper<const DTYPE(data)>, 1>& /*xi*/) \
+      const noexcept;                                                         \
+  template Tensor<DTYPE(data), tmpl::integral_list<std::int32_t, 2, 1>,       \
+                  index_list<SpatialIndex<1, UpLo::Up, Frame::NoFrame>,       \
+                             SpatialIndex<1, UpLo::Lo, Frame::NoFrame>>>      \
+  Equiangular::jacobian(const std::array<DTYPE(data), 1>& /*xi*/)             \
+      const noexcept;                                                         \
+  template Tensor<DTYPE(data), tmpl::integral_list<std::int32_t, 2, 1>,       \
+                  index_list<SpatialIndex<1, UpLo::Up, Frame::NoFrame>,       \
+                             SpatialIndex<1, UpLo::Lo, Frame::NoFrame>>>      \
+  Equiangular::inv_jacobian(                                                  \
+      const std::array<std::reference_wrapper<const DTYPE(data)>, 1>& /*xi*/) \
+      const noexcept;                                                         \
+  template Tensor<DTYPE(data), tmpl::integral_list<std::int32_t, 2, 1>,       \
+                  index_list<SpatialIndex<1, UpLo::Up, Frame::NoFrame>,       \
+                             SpatialIndex<1, UpLo::Lo, Frame::NoFrame>>>      \
+  Equiangular::inv_jacobian(const std::array<DTYPE(data), 1>& /*xi*/)         \
+      const noexcept;
+
+GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector))
+
+#undef DTYPE
+#undef INSTANTIATE
+/// \endcond
+}  // namespace CoordinateMaps

src/Domain/CoordinateMaps/Equiangular.hpp

@@ -0,0 +1,109 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+/// \file
+/// Defines the class Equiangular.
+
+#pragma once
+
+#include <array>
+
+#include "DataStructures/Tensor/Tensor.hpp"
+
+namespace PUP {
+class er;
+}  // namespace PUP
+
+namespace CoordinateMaps {
+
+/*!
+ * \ingroup CoordinateMapsGroup
+ * \brief Non-linear map from \f$\xi \in [A, B]\rightarrow x \in [a, b]\f$.
+ *
+ * The formula for the mapping is:
+ * \f{align}
+ * x &= \frac{a}{2} \left(1-\mathrm{tan}\left(
+ *      \frac{\pi(2\xi-B-A)}{4(B-A)}\right)\right) +
+ *      \frac{b}{2} \left(1+\mathrm{tan}\left(
+ *      \frac{\pi(2\xi-B-A)}{4(B-A)}\right)\right)\\
+ * \xi &= \frac{A}{2} \left(1-\frac{4}{\pi}\mathrm{arctan}\left(
+ *        \frac{2x-a-b}{b-a}\right)\right)+
+ *        \frac{B}{2} \left(1+\frac{4}{\pi}\mathrm{arctan}\left(
+ *        \frac{2x-a-b}{b-a}\right)\right)
+ * \f}
+ *
+ * \note The intermediate step in which a tangent map is applied can be more
+ * clearly understood if we define the coordinates:
+ * \f{align}
+ * \xi_{logical} &:= \frac{2\xi-B-A}{B-A} \in [-1, 1]\\
+ * \Xi &:= \mathrm{tan}\left(\frac{\pi\xi_{logical}}{4}\right) \in [-1, 1]
+ * \f}
+ *
+ * This map is intended to be used with `Wedge2D` and `Wedge3D` when equiangular
+ * coordinates are chosen for those maps. For more information on this choice
+ * of coordinates, see the documentation for Wedge3D.
+ */
+class Equiangular {
+ public:
+  static constexpr size_t dim = 1;
+
+  Equiangular(double A, double B, double a, double b) noexcept;
+
+  Equiangular() = default;
+  ~Equiangular() = default;
+  Equiangular(const Equiangular&) = default;
+  Equiangular(Equiangular&&) noexcept = default;
+  Equiangular& operator=(const Equiangular&) = default;
+  Equiangular& operator=(Equiangular&&) = default;
+
+  template <typename T>
+  std::array<std::decay_t<tt::remove_reference_wrapper_t<T>>, 1> operator()(
+      const std::array<T, 1>& xi) const noexcept;
+
+  template <typename T>
+  std::array<std::decay_t<tt::remove_reference_wrapper_t<T>>, 1> inverse(
+      const std::array<T, 1>& x) const noexcept;
+
+  template <typename T>
+  Tensor<std::decay_t<tt::remove_reference_wrapper_t<T>>,
+         tmpl::integral_list<std::int32_t, 2, 1>,
+         index_list<SpatialIndex<1, UpLo::Up, Frame::NoFrame>,
+                    SpatialIndex<1, UpLo::Lo, Frame::NoFrame>>>
+  inv_jacobian(const std::array<T, 1>& xi) const noexcept;
+
+  template <typename T>
+  Tensor<std::decay_t<tt::remove_reference_wrapper_t<T>>,
+         tmpl::integral_list<std::int32_t, 2, 1>,
+         index_list<SpatialIndex<1, UpLo::Up, Frame::NoFrame>,
+                    SpatialIndex<1, UpLo::Lo, Frame::NoFrame>>>
+  jacobian(const std::array<T, 1>& xi) const noexcept;
+
+  // clang-tidy: google-runtime-references
+  void pup(PUP::er& p) noexcept;  // NOLINT
+
+ private:
+  friend bool operator==(const Equiangular& lhs,
+                         const Equiangular& rhs) noexcept;
+
+  double A_{-1.0};
+  double B_{1.0};
+  double a_{-1.0};
+  double b_{1.0};
+  double length_of_domain_over_m_pi_4_{(B_ - A_) / M_PI_4};  // 4(B-A)/\pi
+  double length_of_range_{2.0};                             // b-a
+  double m_pi_4_over_length_of_domain_{M_PI_4 / (B_ - A_)};
+  double one_over_length_of_range_{0.5};
+  // The jacobian for the affine map with the same parameters.
+  double linear_jacobian_times_m_pi_4_{length_of_range_ /
+                                       length_of_domain_over_m_pi_4_};
+  // The inverse jacobian for the affine map with the same parameters.
+  double linear_inverse_jacobian_over_m_pi_4_{length_of_domain_over_m_pi_4_ /
+                                              length_of_range_};
+};
+
+inline bool operator!=(const CoordinateMaps::Equiangular& lhs,
+                       const CoordinateMaps::Equiangular& rhs) noexcept {
+  return not(lhs == rhs);
+}
+
+}  // namespace CoordinateMaps

tests/Unit/Domain/CMakeLists.txt

@@ -4,6 +4,7 @@
 set(DOMAIN_TESTS
     Domain/CoordinateMaps/Test_AffineMap.cpp
     Domain/CoordinateMaps/Test_CoordinateMap.cpp
+    Domain/CoordinateMaps/Test_Equiangular.cpp
     Domain/CoordinateMaps/Test_Identity.cpp
     Domain/CoordinateMaps/Test_ProductMaps.cpp
     Domain/CoordinateMaps/Test_Rotation.cpp

tests/Unit/Domain/CoordinateMaps/Test_Equiangular.cpp

@@ -0,0 +1,71 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include <catch.hpp>
+
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Domain/CoordinateMaps/Equiangular.hpp"
+#include "Utilities/StdHelpers.hpp"
+#include "tests/Unit/Domain/CoordinateMaps/TestMapHelpers.hpp"
+#include "tests/Unit/TestHelpers.hpp"
+
+SPECTRE_TEST_CASE("Unit.Domain.CoordinateMaps.Equiangular", "[Domain][Unit]") {
+  const double xA = -1.0;
+  const double xB = 2.0;
+  const double xa = -3.0;
+  const double xb = 2.0;
+
+  CoordinateMaps::Equiangular equiangular_map(xA, xB, xa, xb);
+
+  const double xi = 0.5 * (xA + xB);
+  const double x = xb * (xi - xA) / (xB - xA) + xa * (xB - xi) / (xB - xA);
+  const double r1 = 0.172899453;
+  const double r2 = -0.234256219;
+
+  const std::array<double, 1> point_A{{xA}};
+  const std::array<double, 1> point_B{{xB}};
+  const std::array<double, 1> point_a{{xa}};
+  const std::array<double, 1> point_b{{xb}};
+  const std::array<double, 1> point_xi{{xi}};
+  const std::array<double, 1> point_x{{x}};
+  const std::array<double, 1> point_r1{{r1}};
+  const std::array<double, 1> point_r2{{r2}};
+
+  CHECK(equiangular_map(point_A)[0] == approx(point_a[0]));
+  CHECK(equiangular_map(point_B)[0] == approx(point_b[0]));
+  CHECK(equiangular_map(point_xi)[0] == approx(point_x[0]));
+  CHECK(equiangular_map(point_r1)[0] ==
+        approx(-0.5 + 2.5 * tan(M_PI_4 * (2.0 * r1 - 1.0) / 3.0)));
+  CHECK(equiangular_map(point_r2)[0] ==
+        approx(-0.5 + 2.5 * tan(M_PI_4 * (2.0 * r2 - 1.0) / 3.0)));
+
+  CHECK(equiangular_map.inverse(point_a)[0] == approx(point_A[0]));
+  CHECK(equiangular_map.inverse(point_b)[0] == approx(point_B[0]));
+  CHECK(equiangular_map.inverse(point_x)[0] == approx(point_xi[0]));
+  CHECK(equiangular_map.inverse(point_r1)[0] ==
+        approx(0.5 + 3.0 * atan(0.4 * r1 + 0.2) / M_PI_2));
+  CHECK(equiangular_map.inverse(point_r2)[0] ==
+        approx(0.5 + 3.0 * atan(0.4 * r2 + 0.2) / M_PI_2));
+
+  CHECK((get<0, 0>(equiangular_map.inv_jacobian(point_A))) *
+            (get<0, 0>(equiangular_map.jacobian(point_A))) ==
+        approx(1.0));
+  CHECK((get<0, 0>(equiangular_map.inv_jacobian(point_B))) *
+            (get<0, 0>(equiangular_map.jacobian(point_B))) ==
+        approx(1.0));
+  CHECK((get<0, 0>(equiangular_map.inv_jacobian(point_xi))) *
+            (get<0, 0>(equiangular_map.jacobian(point_xi))) ==
+        approx(1.0));
+
+  test_jacobian(equiangular_map, point_xi);
+  test_inv_jacobian(equiangular_map, point_xi);
+  test_inverse_map(equiangular_map, point_xi);
+
+  // Check inequivalence operator
+  CHECK_FALSE(equiangular_map != equiangular_map);
+  test_serialization(equiangular_map);
+
+  test_coordinate_map_argument_types(equiangular_map, point_xi);
+  test_coordinate_map_argument_types(equiangular_map, point_r1);
+  test_coordinate_map_argument_types(equiangular_map, point_r2);
+}