4190347449

Build # 4190347449

Build Type

push

github

Committed by GitHub

Commit Message

Merge pull request #4729 from knelli2/single_horizon_measurement

Run Details

2 of 2 new or added lines in 1 file covered. (100.0%)

63512 of 66219 relevant lines covered (95.91%)

330033.47 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

94.44

/src/Domain/CoordinateMaps/FocallyLiftedEndcap.cpp

// Distributed under the MIT License.
// See LICENSE.txt for details.

#include "Domain/CoordinateMaps/FocallyLiftedEndcap.hpp"

#include <cmath>
#include <limits>
#include <pup.h>

#include "DataStructures/DataVector.hpp"
#include "DataStructures/Tensor/Tensor.hpp"
#include "Domain/CoordinateMaps/CylindricalEndcapHelpers.hpp"
#include "Domain/CoordinateMaps/FocallyLiftedMapHelpers.hpp"
#include "Parallel/PupStlCpp11.hpp"
#include "Utilities/ConstantExpressions.hpp"
#include "Utilities/ContainerHelpers.hpp"
#include "Utilities/DereferenceWrapper.hpp"
#include "Utilities/EqualWithinRoundoff.hpp"
#include "Utilities/GenerateInstantiations.hpp"
#include "Utilities/MakeWithValue.hpp"

namespace domain::CoordinateMaps::FocallyLiftedInnerMaps {

Endcap::Endcap(const std::array<double, 3>& center, const double radius,
               const double z_plane)
    : center_(center),
      radius_([&radius]() {
        // The equal_within_roundoff below has an implicit scale of 1,
        // so the ASSERT may trigger in the case where we really
        // want an entire domain that is very small.
        ASSERT(not equal_within_roundoff(radius, 0.0),
               "Cannot have zero radius");
        return radius;
      }()),
      theta_max_([&center, &radius, &z_plane]() {
        const double cos_theta_max = (z_plane - center[2]) / radius;
        ASSERT(abs(cos_theta_max) < 1.0,
               "Plane must intersect sphere, and at more than one point");
        return acos(cos_theta_max);
      }()) {}

template <typename T>
void Endcap::forward_map(
    const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
        target_coords,
    const std::array<T, 3>& source_coords) const {
  using ReturnType = tt::remove_cvref_wrap_t<T>;
  const ReturnType& xbar = source_coords[0];
  const ReturnType& ybar = source_coords[1];
  ReturnType& x = (*target_coords)[0];
  ReturnType& y = (*target_coords)[1];
  ReturnType& z = (*target_coords)[2];
  // Use z and y as temporary storage to avoid allocations,
  // before setting them to their actual values.
  z = sqrt(square(xbar) + square(ybar));
  y = cylindrical_endcap_helpers::sin_ax_over_x(z, theta_max_) * radius_;
  x = y * xbar + center_[0];
  y = y * ybar + center_[1];
  z = radius_ * cos(z * theta_max_) + center_[2];
}

template <typename T>
void Endcap::jacobian(const gsl::not_null<tnsr::Ij<tt::remove_cvref_wrap_t<T>,
                                                   3, Frame::NoFrame>*>
                          jacobian_out,
                      const std::array<T, 3>& source_coords) const {
  using ReturnType = tt::remove_cvref_wrap_t<T>;
  const ReturnType& xbar = source_coords[0];
  const ReturnType& ybar = source_coords[1];

  if constexpr (not std::is_same<ReturnType, double>::value) {
    destructive_resize_components(
        jacobian_out, static_cast<ReturnType>(source_coords[0]).size());
  }
  // Most of the jacobian components are zero.
  for (auto& jac_component : *jacobian_out) {
    jac_component = 0.0;
  }

  // Use parts of Jacobian as temp storage to reduce allocations.
  get<1, 1>(*jacobian_out) = sqrt(square(xbar) + square(ybar));
  get<1, 0>(*jacobian_out) =
      radius_ * cylindrical_endcap_helpers::sin_ax_over_x(
                    get<1, 1>(*jacobian_out), theta_max_);
  // 1/rhobar d/d(rhobar) [(sin (rhobar theta_max)/rhobar)]
  get<0, 1>(*jacobian_out) =
      radius_ * cylindrical_endcap_helpers::one_over_x_d_sin_ax_over_x(
                    get<1, 1>(*jacobian_out), theta_max_);

  // dy/dybar
  get<1, 1>(*jacobian_out) =
      get<0, 1>(*jacobian_out) * square(ybar) + get<1, 0>(*jacobian_out);
  // dx/dxbar
  get<0, 0>(*jacobian_out) =
      get<0, 1>(*jacobian_out) * square(xbar) + get<1, 0>(*jacobian_out);

  // Still a temporary variable here.
  get<1, 0>(*jacobian_out) *= -theta_max_;

  // dz/dxbar
  get<2, 0>(*jacobian_out) = get<1, 0>(*jacobian_out) * xbar;
  // dz/dybar
  get<2, 1>(*jacobian_out) = get<1, 0>(*jacobian_out) * ybar;

  // dx/dybar
  get<0, 1>(*jacobian_out) *= xbar * ybar;
  // dy/dxbar
  get<1, 0>(*jacobian_out) = get<0, 1>(*jacobian_out);
}

template <typename T>
void Endcap::inv_jacobian(
    const gsl::not_null<
        tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>*>
        inv_jacobian_out,
    const std::array<T, 3>& source_coords) const {
  using ReturnType = tt::remove_cvref_wrap_t<T>;
  const ReturnType& xbar = source_coords[0];
  const ReturnType& ybar = source_coords[1];

  if constexpr (not std::is_same<ReturnType, double>::value) {
    destructive_resize_components(
        inv_jacobian_out, static_cast<ReturnType>(source_coords[0]).size());
  }
  // Most of the inverse jacobian components are zero.
  for (auto& jac_component : *inv_jacobian_out) {
    jac_component = 0.0;
  }

  // Use parts of Jacobian as temp storage to reduce allocations.
  // We comment temporary quantities to make the code more understandable.

  // rhobar, Eq. 7 in the documentation.
  get<1, 0>(*inv_jacobian_out) = sqrt(square(xbar) + square(ybar));

  // q = sin(rhobar theta)/rhobar as defined in the documentation, Eq. 17.
  get<0, 1>(*inv_jacobian_out) = cylindrical_endcap_helpers::sin_ax_over_x(
      get<1, 0>(*inv_jacobian_out), theta_max_);

  // 1/rhobar dq/d(rhobar)
  get<1, 1>(*inv_jacobian_out) =
      cylindrical_endcap_helpers::one_over_x_d_sin_ax_over_x(
          get<1, 0>(*inv_jacobian_out), theta_max_);

  // Right-hand side of Eq. 23, without the factor of
  // xbar ybar/rq or the minus sign.
  get<1, 0>(*inv_jacobian_out) =
      get<1, 1>(*inv_jacobian_out) /
      (get<0, 1>(*inv_jacobian_out) +
       square(get<1, 0>(*inv_jacobian_out)) * get<1, 1>(*inv_jacobian_out));

  // 1/(r q), i.e. first term in Eq. 22.
  get<0, 0>(*inv_jacobian_out) = 1.0 / (get<0, 1>(*inv_jacobian_out) * radius_);

  // Right-hand side of Eq. 23, without the factor of
  // xbar ybar or the minus sign.
  get<1, 0>(*inv_jacobian_out) *= get<0, 0>(*inv_jacobian_out);

  // dybar/dy
  get<1, 1>(*inv_jacobian_out) = get<0, 0>(*inv_jacobian_out) -
                                 square(ybar) * get<1, 0>(*inv_jacobian_out);
  // dxbar/dx
  get<0, 0>(*inv_jacobian_out) -= square(xbar) * get<1, 0>(*inv_jacobian_out);
  // dybar/dx
  get<1, 0>(*inv_jacobian_out) *= -xbar * ybar;
  // dxbar/dy
  get<0, 1>(*inv_jacobian_out) = get<1, 0>(*inv_jacobian_out);
}

std::optional<std::array<double, 3>> Endcap::inverse(
    const std::array<double, 3>& target_coords, const double sigma_in) const {
  const double x = target_coords[0] - center_[0];
  const double y = target_coords[1] - center_[1];
  const double z = target_coords[2] - center_[2];

  // Are we in the range of the map?
  // The equal_within_roundoff below has an implicit scale of 1,
  // so the inverse may fail if radius_ is very small on purpose,
  // e.g. if we really want a tiny tiny domain for some reason.
  const double r = sqrt(square(x) + square(y) + square(z));
  if (not equal_within_roundoff(r, radius_)) {
    return std::optional<std::array<double, 3>>{};
  }

  // Compute zbar and check if we are in the range of the map.
  const double zbar = 2.0 * sigma_in - 1.0;
  if (abs(zbar) > 1.0 and not equal_within_roundoff(abs(zbar), 1.0)) {
    return std::optional<std::array<double, 3>>{};
  }

  const double rho = sqrt(square(x) + square(y));
  if (UNLIKELY(rho == 0.0)) {
    // If x and y are zero, so are xbar and ybar,
    // so we are done.
    return std::array<double, 3>{{0.0, 0.0, zbar}};
  }

  // Note: theta_max_ cannot be zero for a nonsingular map.
  const double rhobar = atan2(rho, z) / theta_max_;

  // Check if we are outside the range of the map.
  if (rhobar > 1.0 and not equal_within_roundoff(rhobar, 1.0)) {
    return std::optional<std::array<double, 3>>{};
  }

  const double xbar = x * rhobar / rho;
  const double ybar = y * rhobar / rho;
  return std::array<double, 3>{{xbar, ybar, zbar}};
}

template <typename T>
void Endcap::sigma(const gsl::not_null<tt::remove_cvref_wrap_t<T>*> sigma_out,
                   const std::array<T, 3>& source_coords) const {
  *sigma_out = 0.5 * (source_coords[2] + 1.0);
}

template <typename T>
void Endcap::deriv_sigma(
    const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
        deriv_sigma_out,
    const std::array<T, 3>& source_coords) const {
  using ReturnType = tt::remove_cvref_wrap_t<T>;
  if constexpr (not std::is_same<ReturnType, double>::value) {
    destructive_resize_components(
        deriv_sigma_out, static_cast<ReturnType>(source_coords[0]).size());
  }
  (*deriv_sigma_out)[0] = 0.0;
  (*deriv_sigma_out)[1] = 0.0;
  (*deriv_sigma_out)[2] = 0.5;
}

template <typename T>
void Endcap::dxbar_dsigma(
    const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
        dxbar_dsigma_out,
    const std::array<T, 3>& source_coords) const {
  using ReturnType = tt::remove_cvref_wrap_t<T>;
  if constexpr (not std::is_same<ReturnType, double>::value) {
    destructive_resize_components(
        dxbar_dsigma_out, static_cast<ReturnType>(source_coords[0]).size());
  }
  (*dxbar_dsigma_out)[0] = 0.0;
  (*dxbar_dsigma_out)[1] = 0.0;
  (*dxbar_dsigma_out)[2] = 2.0;
}

std::optional<double> Endcap::lambda_tilde(
    const std::array<double, 3>& parent_mapped_target_coords,
    const std::array<double, 3>& projection_point,
    const bool source_is_between_focus_and_target) const {
  // Try to find lambda_tilde going from target_coords to sphere.
  // If the target surface is outside the sphere (that is,
  // source_is_between_focus_and_target is true), then lambda_tilde should be
  // positive and less than or equal to unity. If the target surface is inside
  // the sphere, then lambda_tilde should be greater than or equal
  // to unity. If there are two such roots, we choose based on where the points
  // are.
  const bool choose_larger_root =
      parent_mapped_target_coords[2] > projection_point[2];
  return FocallyLiftedMapHelpers::try_scale_factor(
      parent_mapped_target_coords, projection_point, center_, radius_,
      choose_larger_root, not source_is_between_focus_and_target);
}

template <typename T>
void Endcap::deriv_lambda_tilde(
    const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
        deriv_lambda_tilde_out,
    const std::array<T, 3>& target_coords, const T& lambda_tilde,
    const std::array<double, 3>& projection_point) const {
  FocallyLiftedMapHelpers::d_scale_factor_d_src_point(
      deriv_lambda_tilde_out, target_coords, projection_point, center_,
      lambda_tilde);
}

void Endcap::pup(PUP::er& p) {
  p | center_;
  p | radius_;
  p | theta_max_;
}

bool operator==(const Endcap& lhs, const Endcap& rhs) {
  return lhs.center_ == rhs.center_ and lhs.radius_ == rhs.radius_ and
         lhs.theta_max_ == rhs.theta_max_;
}

bool operator!=(const Endcap& lhs, const Endcap& rhs) {
  return not(lhs == rhs);
}
// Explicit instantiations
#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)

#define INSTANTIATE(_, data)                                                  \
  template void Endcap::forward_map(                                          \
      const gsl::not_null<                                                    \
          std::array<tt::remove_cvref_wrap_t<DTYPE(data)>, 3>*>               \
          target_coords,                                                      \
      const std::array<DTYPE(data), 3>& source_coords) const;                 \
  template void Endcap::jacobian(                                             \
      const gsl::not_null<                                                    \
          tnsr::Ij<tt::remove_cvref_wrap_t<DTYPE(data)>, 3, Frame::NoFrame>*> \
          jacobian_out,                                                       \
      const std::array<DTYPE(data), 3>& source_coords) const;                 \
  template void Endcap::inv_jacobian(                                         \
      const gsl::not_null<                                                    \
          tnsr::Ij<tt::remove_cvref_wrap_t<DTYPE(data)>, 3, Frame::NoFrame>*> \
          inv_jacobian_out,                                                   \
      const std::array<DTYPE(data), 3>& source_coords) const;                 \
  template void Endcap::sigma(                                                \
      const gsl::not_null<tt::remove_cvref_wrap_t<DTYPE(data)>*> sigma_out,   \
      const std::array<DTYPE(data), 3>& source_coords) const;                 \
  template void Endcap::deriv_sigma(                                          \
      const gsl::not_null<                                                    \
          std::array<tt::remove_cvref_wrap_t<DTYPE(data)>, 3>*>               \
          deriv_sigma_out,                                                    \
      const std::array<DTYPE(data), 3>& source_coords) const;                 \
  template void Endcap::dxbar_dsigma(                                         \
      const gsl::not_null<                                                    \
          std::array<tt::remove_cvref_wrap_t<DTYPE(data)>, 3>*>               \
          dxbar_dsigma_out,                                                   \
      const std::array<DTYPE(data), 3>& source_coords) const;                 \
  template void Endcap::deriv_lambda_tilde(                                   \
      const gsl::not_null<                                                    \
          std::array<tt::remove_cvref_wrap_t<DTYPE(data)>, 3>*>               \
          deriv_lambda_tilde_out,                                             \
      const std::array<DTYPE(data), 3>& target_coords,                        \
      const DTYPE(data) & lambda_tilde,                                       \
      const std::array<double, 3>& projection_point) const;

GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector,
                                      std::reference_wrapper<const double>,
                                      std::reference_wrapper<const DataVector>))

#undef INSTANTIATE
#undef DTYPE

}  // namespace domain::CoordinateMaps::FocallyLiftedInnerMaps

1	// Distributed under the MIT License.
2	// See LICENSE.txt for details.
3
4	#include "Domain/CoordinateMaps/FocallyLiftedEndcap.hpp"
5
6	#include <cmath>
7	#include <limits>
8	#include <pup.h>
9
10	#include "DataStructures/DataVector.hpp"
11	#include "DataStructures/Tensor/Tensor.hpp"
12	#include "Domain/CoordinateMaps/CylindricalEndcapHelpers.hpp"
13	#include "Domain/CoordinateMaps/FocallyLiftedMapHelpers.hpp"
14	#include "Parallel/PupStlCpp11.hpp"
15	#include "Utilities/ConstantExpressions.hpp"
16	#include "Utilities/ContainerHelpers.hpp"
17	#include "Utilities/DereferenceWrapper.hpp"
18	#include "Utilities/EqualWithinRoundoff.hpp"
19	#include "Utilities/GenerateInstantiations.hpp"
20	#include "Utilities/MakeWithValue.hpp"
21
22	namespace domain::CoordinateMaps::FocallyLiftedInnerMaps {
23
24	Endcap::Endcap(const std::array<double, 3>& center, const double radius,	2✔
25	const double z_plane)	2✔
26	: center_(center),
27	radius_([&radius]() {	6✔
28	// The equal_within_roundoff below has an implicit scale of 1,
29	// so the ASSERT may trigger in the case where we really
30	// want an entire domain that is very small.
31	ASSERT(not equal_within_roundoff(radius, 0.0),	4✔
32	"Cannot have zero radius");
33	return radius;	2✔
34	}()),	×
35	theta_max_([&center, &radius, &z_plane]() {	4✔
36	const double cos_theta_max = (z_plane - center[2]) / radius;	2✔
37	ASSERT(abs(cos_theta_max) < 1.0,	2✔
38	"Plane must intersect sphere, and at more than one point");
39	return acos(cos_theta_max);	2✔
40	}()) {}	2✔
41
42	template <typename T>
43	void Endcap::forward_map(	984✔
44	const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
45	target_coords,
46	const std::array<T, 3>& source_coords) const {
47	using ReturnType = tt::remove_cvref_wrap_t<T>;
48	const ReturnType& xbar = source_coords[0];	984✔
49	const ReturnType& ybar = source_coords[1];	984✔
50	ReturnType& x = (*target_coords)[0];	984✔
51	ReturnType& y = (*target_coords)[1];	984✔
52	ReturnType& z = (*target_coords)[2];	984✔
53	// Use z and y as temporary storage to avoid allocations,
54	// before setting them to their actual values.
55	z = sqrt(square(xbar) + square(ybar));	1,128✔
56	y = cylindrical_endcap_helpers::sin_ax_over_x(z, theta_max_) * radius_;	984✔
57	x = y * xbar + center_[0];	984✔
58	y = y * ybar + center_[1];	984✔
59	z = radius_ * cos(z * theta_max_) + center_[2];	984✔
60	}	984✔
61
62	template <typename T>
63	void Endcap::jacobian(const gsl::not_null<tnsr::Ij<tt::remove_cvref_wrap_t<T>,	184✔
64	3, Frame::NoFrame>*>
65	jacobian_out,
66	const std::array<T, 3>& source_coords) const {
67	using ReturnType = tt::remove_cvref_wrap_t<T>;
68	const ReturnType& xbar = source_coords[0];	184✔
69	const ReturnType& ybar = source_coords[1];	184✔
70
71	if constexpr (not std::is_same<ReturnType, double>::value) {
72	destructive_resize_components(	96✔
73	jacobian_out, static_cast<ReturnType>(source_coords[0]).size());	48✔
74	}
75	// Most of the jacobian components are zero.
76	for (auto& jac_component : *jacobian_out) {	1,840✔
77	jac_component = 0.0;	1,656✔
78	}
79
80	// Use parts of Jacobian as temp storage to reduce allocations.
81	get<1, 1>(*jacobian_out) = sqrt(square(xbar) + square(ybar));	416✔
82	get<1, 0>(*jacobian_out) =	232✔
83	radius_ * cylindrical_endcap_helpers::sin_ax_over_x(	232✔
84	get<1, 1>(*jacobian_out), theta_max_);	368✔
85	// 1/rhobar d/d(rhobar) [(sin (rhobar theta_max)/rhobar)]
86	get<0, 1>(*jacobian_out) =	232✔
87	radius_ * cylindrical_endcap_helpers::one_over_x_d_sin_ax_over_x(	232✔
88	get<1, 1>(*jacobian_out), theta_max_);	368✔
89
90	// dy/dybar
91	get<1, 1>(*jacobian_out) =	232✔
92	get<0, 1>(jacobian_out) square(ybar) + get<1, 0>(*jacobian_out);	736✔
93	// dx/dxbar
94	get<0, 0>(*jacobian_out) =	232✔
95	get<0, 1>(jacobian_out) square(xbar) + get<1, 0>(*jacobian_out);	736✔
96
97	// Still a temporary variable here.
98	get<1, 0>(jacobian_out) = -theta_max_;	232✔
99
100	// dz/dxbar
101	get<2, 0>(jacobian_out) = get<1, 0>(jacobian_out) * xbar;	416✔
102	// dz/dybar
103	get<2, 1>(jacobian_out) = get<1, 0>(jacobian_out) * ybar;	416✔
104
105	// dx/dybar
106	get<0, 1>(jacobian_out) = xbar * ybar;	232✔
107	// dy/dxbar
108	get<1, 0>(jacobian_out) = get<0, 1>(jacobian_out);	416✔
109	}	184✔
110
111	template <typename T>
112	void Endcap::inv_jacobian(	184✔
113	const gsl::not_null<
114	tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>*>
115	inv_jacobian_out,
116	const std::array<T, 3>& source_coords) const {
117	using ReturnType = tt::remove_cvref_wrap_t<T>;
118	const ReturnType& xbar = source_coords[0];	184✔
119	const ReturnType& ybar = source_coords[1];	184✔
120
121	if constexpr (not std::is_same<ReturnType, double>::value) {
122	destructive_resize_components(	96✔
123	inv_jacobian_out, static_cast<ReturnType>(source_coords[0]).size());	48✔
124	}
125	// Most of the inverse jacobian components are zero.
126	for (auto& jac_component : *inv_jacobian_out) {	1,840✔
127	jac_component = 0.0;	1,656✔
128	}
129
130	// Use parts of Jacobian as temp storage to reduce allocations.
131	// We comment temporary quantities to make the code more understandable.
132
133	// rhobar, Eq. 7 in the documentation.
134	get<1, 0>(*inv_jacobian_out) = sqrt(square(xbar) + square(ybar));	416✔
135
136	// q = sin(rhobar theta)/rhobar as defined in the documentation, Eq. 17.
137	get<0, 1>(*inv_jacobian_out) = cylindrical_endcap_helpers::sin_ax_over_x(	232✔
138	get<1, 0>(*inv_jacobian_out), theta_max_);	368✔
139
140	// 1/rhobar dq/d(rhobar)
141	get<1, 1>(*inv_jacobian_out) =	232✔
142	cylindrical_endcap_helpers::one_over_x_d_sin_ax_over_x(	136✔
143	get<1, 0>(*inv_jacobian_out), theta_max_);	368✔
144
145	// Right-hand side of Eq. 23, without the factor of
146	// xbar ybar/rq or the minus sign.
147	get<1, 0>(*inv_jacobian_out) =	232✔
148	get<1, 1>(*inv_jacobian_out) /	232✔
149	(get<0, 1>(*inv_jacobian_out) +	232✔
150	square(get<1, 0>(inv_jacobian_out)) get<1, 1>(*inv_jacobian_out));	736✔
151
152	// 1/(r q), i.e. first term in Eq. 22.
153	get<0, 0>(inv_jacobian_out) = 1.0 / (get<0, 1>(inv_jacobian_out) * radius_);	416✔
154
155	// Right-hand side of Eq. 23, without the factor of
156	// xbar ybar or the minus sign.
157	get<1, 0>(inv_jacobian_out) = get<0, 0>(*inv_jacobian_out);	416✔
158
159	// dybar/dy
160	get<1, 1>(inv_jacobian_out) = get<0, 0>(inv_jacobian_out) -	688✔
161	square(ybar) * get<1, 0>(*inv_jacobian_out);	416✔
162	// dxbar/dx
163	get<0, 0>(inv_jacobian_out) -= square(xbar) get<1, 0>(*inv_jacobian_out);	464✔
164	// dybar/dx
165	get<1, 0>(inv_jacobian_out) = -xbar * ybar;	232✔
166	// dxbar/dy
167	get<0, 1>(inv_jacobian_out) = get<1, 0>(inv_jacobian_out);	416✔
168	}	184✔
169
170	std::optional<std::array<double, 3>> Endcap::inverse(	24✔
171	const std::array<double, 3>& target_coords, const double sigma_in) const {
172	const double x = target_coords[0] - center_[0];	24✔
173	const double y = target_coords[1] - center_[1];	24✔
174	const double z = target_coords[2] - center_[2];	24✔
175
176	// Are we in the range of the map?
177	// The equal_within_roundoff below has an implicit scale of 1,
178	// so the inverse may fail if radius_ is very small on purpose,
179	// e.g. if we really want a tiny tiny domain for some reason.
180	const double r = sqrt(square(x) + square(y) + square(z));	48✔
181	if (not equal_within_roundoff(r, radius_)) {	48✔
182	return std::optional<std::array<double, 3>>{};	×
183	}
184
185	// Compute zbar and check if we are in the range of the map.
186	const double zbar = 2.0 * sigma_in - 1.0;	24✔
187	if (abs(zbar) > 1.0 and not equal_within_roundoff(abs(zbar), 1.0)) {	24✔
188	return std::optional<std::array<double, 3>>{};	×
189	}
190
191	const double rho = sqrt(square(x) + square(y));	24✔
192	if (UNLIKELY(rho == 0.0)) {	24✔
193	// If x and y are zero, so are xbar and ybar,
194	// so we are done.
195	return std::array<double, 3>{{0.0, 0.0, zbar}};	×
196	}
197
198	// Note: theta_max_ cannot be zero for a nonsingular map.
199	const double rhobar = atan2(rho, z) / theta_max_;	24✔
200
201	// Check if we are outside the range of the map.
202	if (rhobar > 1.0 and not equal_within_roundoff(rhobar, 1.0)) {	26✔
203	return std::optional<std::array<double, 3>>{};	×
204	}
205
206	const double xbar = x * rhobar / rho;	24✔
207	const double ybar = y * rhobar / rho;	24✔
208	return std::array<double, 3>{{xbar, ybar, zbar}};	48✔
209	}
210
211	template <typename T>
212	void Endcap::sigma(const gsl::not_null<tt::remove_cvref_wrap_t<T>*> sigma_out,	984✔
213	const std::array<T, 3>& source_coords) const {
214	sigma_out = 0.5 (source_coords[2] + 1.0);	1,128✔
215	}	984✔
216
217	template <typename T>
218	void Endcap::deriv_sigma(	184✔
219	const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
220	deriv_sigma_out,
221	const std::array<T, 3>& source_coords) const {
222	using ReturnType = tt::remove_cvref_wrap_t<T>;
223	if constexpr (not std::is_same<ReturnType, double>::value) {
224	destructive_resize_components(	96✔
225	deriv_sigma_out, static_cast<ReturnType>(source_coords[0]).size());	48✔
226	}
227	(*deriv_sigma_out)[0] = 0.0;	184✔
228	(*deriv_sigma_out)[1] = 0.0;	184✔
229	(*deriv_sigma_out)[2] = 0.5;	184✔
230	}	184✔
231
232	template <typename T>
233	void Endcap::dxbar_dsigma(	184✔
234	const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
235	dxbar_dsigma_out,
236	const std::array<T, 3>& source_coords) const {
237	using ReturnType = tt::remove_cvref_wrap_t<T>;
238	if constexpr (not std::is_same<ReturnType, double>::value) {
239	destructive_resize_components(	96✔
240	dxbar_dsigma_out, static_cast<ReturnType>(source_coords[0]).size());	48✔
241	}
242	(*dxbar_dsigma_out)[0] = 0.0;	184✔
243	(*dxbar_dsigma_out)[1] = 0.0;	184✔
244	(*dxbar_dsigma_out)[2] = 2.0;	184✔
245	}	184✔
246
247	std::optional<double> Endcap::lambda_tilde(	24✔
248	const std::array<double, 3>& parent_mapped_target_coords,
249	const std::array<double, 3>& projection_point,
250	const bool source_is_between_focus_and_target) const {
251	// Try to find lambda_tilde going from target_coords to sphere.
252	// If the target surface is outside the sphere (that is,
253	// source_is_between_focus_and_target is true), then lambda_tilde should be
254	// positive and less than or equal to unity. If the target surface is inside
255	// the sphere, then lambda_tilde should be greater than or equal
256	// to unity. If there are two such roots, we choose based on where the points
257	// are.
258	const bool choose_larger_root =
259	parent_mapped_target_coords[2] > projection_point[2];	24✔
260	return FocallyLiftedMapHelpers::try_scale_factor(
261	parent_mapped_target_coords, projection_point, center_, radius_,	24✔
262	choose_larger_root, not source_is_between_focus_and_target);	24✔
263	}
264
265	template <typename T>
266	void Endcap::deriv_lambda_tilde(	184✔
267	const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>
268	deriv_lambda_tilde_out,
269	const std::array<T, 3>& target_coords, const T& lambda_tilde,
270	const std::array<double, 3>& projection_point) const {
271	FocallyLiftedMapHelpers::d_scale_factor_d_src_point(	184✔
272	deriv_lambda_tilde_out, target_coords, projection_point, center_,	184✔
273	lambda_tilde);
274	}	184✔
275
276	void Endcap::pup(PUP::er& p) {	18✔
277	p \| center_;	18✔
278	p \| radius_;	18✔
279	p \| theta_max_;	18✔
280	}	18✔
281
282	bool operator==(const Endcap& lhs, const Endcap& rhs) {	8✔
283	return lhs.center_ == rhs.center_ and lhs.radius_ == rhs.radius_ and	16✔
284	lhs.theta_max_ == rhs.theta_max_;	16✔
285	}
286
287	bool operator!=(const Endcap& lhs, const Endcap& rhs) {	×
288	return not(lhs == rhs);	×
289	}
290	// Explicit instantiations
291	#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
292
293	#define INSTANTIATE(_, data) \
294	template void Endcap::forward_map( \
295	const gsl::not_null< \
296	std::array<tt::remove_cvref_wrap_t<DTYPE(data)>, 3>*> \
297	target_coords, \
298	const std::array<DTYPE(data), 3>& source_coords) const; \
299	template void Endcap::jacobian( \
300	const gsl::not_null< \
301	tnsr::Ij<tt::remove_cvref_wrap_t<DTYPE(data)>, 3, Frame::NoFrame>*> \
302	jacobian_out, \
303	const std::array<DTYPE(data), 3>& source_coords) const; \
304	template void Endcap::inv_jacobian( \
305	const gsl::not_null< \
306	tnsr::Ij<tt::remove_cvref_wrap_t<DTYPE(data)>, 3, Frame::NoFrame>*> \
307	inv_jacobian_out, \
308	const std::array<DTYPE(data), 3>& source_coords) const; \
309	template void Endcap::sigma( \
310	const gsl::not_null<tt::remove_cvref_wrap_t<DTYPE(data)>*> sigma_out, \
311	const std::array<DTYPE(data), 3>& source_coords) const; \
312	template void Endcap::deriv_sigma( \
313	const gsl::not_null< \
314	std::array<tt::remove_cvref_wrap_t<DTYPE(data)>, 3>*> \
315	deriv_sigma_out, \
316	const std::array<DTYPE(data), 3>& source_coords) const; \
317	template void Endcap::dxbar_dsigma( \
318	const gsl::not_null< \
319	std::array<tt::remove_cvref_wrap_t<DTYPE(data)>, 3>*> \
320	dxbar_dsigma_out, \
321	const std::array<DTYPE(data), 3>& source_coords) const; \
322	template void Endcap::deriv_lambda_tilde( \
323	const gsl::not_null< \
324	std::array<tt::remove_cvref_wrap_t<DTYPE(data)>, 3>*> \
325	deriv_lambda_tilde_out, \
326	const std::array<DTYPE(data), 3>& target_coords, \
327	const DTYPE(data) & lambda_tilde, \
328	const std::array<double, 3>& projection_point) const;
329
330	GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector,
331	std::reference_wrapper<const double>,
332	std::reference_wrapper<const DataVector>))
333
334	#undef INSTANTIATE
335	#undef DTYPE
336
337	} // namespace domain::CoordinateMaps::FocallyLiftedInnerMaps

sxs-collaboration / spectre / 4190347449

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous