4138821718

Build Type

push

github

Committed by GitHub

Commit Message

Merge pull request #4714 from nilsdeppe/fix_jac_diag_missing_include

Run Details

63323 of 66009 relevant lines covered (95.93%)

403761.71 hits per line

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

77.08

/src/Domain/CoordinateMaps/FocallyLiftedMapHelpers.cpp

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

#include "Domain/CoordinateMaps/FocallyLiftedMapHelpers.hpp"

#include <cmath>
#include <gsl/gsl_poly.h>
#include <limits>

#include "DataStructures/DataVector.hpp"
#include "NumericalAlgorithms/RootFinding/QuadraticEquation.hpp"
#include "Utilities/ConstantExpressions.hpp"
#include "Utilities/DereferenceWrapper.hpp"
#include "Utilities/EqualWithinRoundoff.hpp"
#include "Utilities/ErrorHandling/Assert.hpp"
#include "Utilities/GenerateInstantiations.hpp"
#include "Utilities/MakeWithValue.hpp"

namespace domain::CoordinateMaps::FocallyLiftedMapHelpers {

template <typename T>
void scale_factor(const gsl::not_null<tt::remove_cvref_wrap_t<T>*>& result,
                  const std::array<T, 3>& src_point,
                  const std::array<double, 3>& proj_center,
                  const std::array<double, 3>& sphere_center, double radius,
                  const bool src_is_between_proj_and_target) {
  using ReturnType = tt::remove_cvref_wrap_t<T>;
  // quadratic equation is
  // a scale_factor^2 + b scale_factor + c = 0
  //
  // For this case
  // a = |x_0^i-P^i|^2
  // b = 2(x_0^i-P^i)(P^j-C^j)\delta_{ij}
  // c = |P^i-C^i|^2 - R^2
  //
  const ReturnType a = square(src_point[0] - proj_center[0]) +
                       square(src_point[1] - proj_center[1]) +
                       square(src_point[2] - proj_center[2]);
  const ReturnType b =
      2.0 *
      ((src_point[0] - proj_center[0]) * (proj_center[0] - sphere_center[0]) +
       (src_point[1] - proj_center[1]) * (proj_center[1] - sphere_center[1]) +
       (src_point[2] - proj_center[2]) * (proj_center[2] - sphere_center[2]));
  const double c = square(sphere_center[0] - proj_center[0]) +
                   square(sphere_center[1] - proj_center[1]) +
                   square(sphere_center[2] - proj_center[2]) - square(radius);
  if (src_is_between_proj_and_target) {
    // Here we assume that src_point is between proj_center and
    // target_point.  There are three cases: 1) src and proj are both
    // inside the sphere, 2) src and proj are both outside the sphere,
    // and 3) proj is outside the sphere and src is inside the sphere.
    // Note that a root of zero corresponds to target_point = proj_center,
    // and a root of unity corresponds to target_point = src_point.
    // To cover all 3 cases, we choose the smallest root that is
    // greater than or equal to unity. This means for case 2) we are
    // choosing the point closest to src.
    *result = smallest_root_greater_than_value_within_roundoff(
        a, b, make_with_value<ReturnType>(a, c), 1.0);
  } else {
    // Here we assume that target_point is between proj_center and
    // src_point. There are three cases: 1) proj is inside the sphere
    // and src is outside the sphere 2) src is inside the sphere and proj
    // is outside the sphere, and 3) the sphere is between src and proj.
    // Note that a root of zero corresponds to target_point = proj_center,
    // and a root of unity corresponds to target_point = src_point.
    // To cover all 3 cases, we choose the largest root that is less than
    // or equal to unity, and we require that this root is positive.
    // This means that for 3) we are choosing the point closest to src.
    *result = largest_root_between_values_within_roundoff(
        a, b, make_with_value<ReturnType>(a, c), 0.0, 1.0);
  }
}

std::optional<double> try_scale_factor(
    const std::array<double, 3>& src_point,
    const std::array<double, 3>& proj_center,
    const std::array<double, 3>& sphere_center, double radius,
    const bool pick_larger_root, const bool pick_root_greater_than_one,
    const bool solve_for_root_minus_one) {
  double x0 = std::numeric_limits<double>::signaling_NaN();
  double x1 = std::numeric_limits<double>::signaling_NaN();
  int num_real_roots = 0;
  if (solve_for_root_minus_one) {
    // We solve the quadratic for (scale_factor-1) instead of
    // scale_factor to avoid roundoff problems when scale_factor is very
    // nearly equal to unity.  Note that scale_factor==1 will occur in
    // the inverse map when src_point (which is x^i) is
    // on the sphere.
    //
    // Roundoff is not a problem when forward-mapping and solving only
    // for lambda.

    // If the original quadratic equation is
    //
    // A scale_factor^2 + B scale_factor + C = 0,
    //
    // and if x = scale_factor-1, then the quadratic equation for x is
    //
    // a x^2 + b x + c = 0
    //
    // where
    //
    // a = A
    // b = 2A + B
    // c = A + B + C
    //
    // For this case
    // A = |x_0^i-P^i|^2
    // B = 2(x_0^i-P^i)(P^j-C^j)\delta_{ij}
    // C = |P^i-C^i|^2 - R^2
    //
    // Note that A + B + C is |x_0^i-P^i+P^i-C^i|^2 - R^2
    // and 2A+B is 2(x_0^i-P^i)(x_0^j-P^j+P^j-C^j), so
    //
    // a = |x_0^i-P^i|^2
    // b = 2(x_0^i-P^i)(x_0^j-C^j)\delta_{ij}
    // c = |x_0^i-C^i|^2 - R^2

    const double a = square(src_point[0] - proj_center[0]) +
                     square(src_point[1] - proj_center[1]) +
                     square(src_point[2] - proj_center[2]);
    const double b =
        2.0 *
        ((src_point[0] - proj_center[0]) * (src_point[0] - sphere_center[0]) +
         (src_point[1] - proj_center[1]) * (src_point[1] - sphere_center[1]) +
         (src_point[2] - proj_center[2]) * (src_point[2] - sphere_center[2]));
    const double c = square(sphere_center[0] - src_point[0]) +
                     square(sphere_center[1] - src_point[1]) +
                     square(sphere_center[2] - src_point[2]) - square(radius);
    num_real_roots = gsl_poly_solve_quadratic(a, b, c, &x0, &x1);
  } else {
    const double a = square(src_point[0] - proj_center[0]) +
                     square(src_point[1] - proj_center[1]) +
                     square(src_point[2] - proj_center[2]);
    const double b =
        2.0 *
        ((src_point[0] - proj_center[0]) * (proj_center[0] - sphere_center[0]) +
         (src_point[1] - proj_center[1]) * (proj_center[1] - sphere_center[1]) +
         (src_point[2] - proj_center[2]) * (proj_center[2] - sphere_center[2]));
    const double c = square(proj_center[0] - sphere_center[0]) +
                     square(proj_center[1] - sphere_center[1]) +
                     square(proj_center[2] - sphere_center[2]) - square(radius);
    num_real_roots = gsl_poly_solve_quadratic(a, b, c, &x0, &x1);
  }
  if (num_real_roots == 2) {
    if(solve_for_root_minus_one) {
      // We solved for scale_factor-1 above, so add 1 to get scale_factor.
      x0 += 1.0;
      x1 += 1.0;
    }
    if (UNLIKELY(equal_within_roundoff(x0, 1.0))) {
      x0 = 1.0;
    }
    if (UNLIKELY(equal_within_roundoff(x1, 1.0))) {
      x1 = 1.0;
    }
    if (pick_root_greater_than_one) {
      // For the inverse map, we want the a scale_factor s such that
      // s >= 1. Note that gsl_poly_solve_quadratic returns x0 < x1.
      // have three cases:
      //  a) x0 < x1 < 1         ->   no root
      //  b) x0 < 1 <= x1        ->   Choose x1
      //  c) 1 <= x0 < x1        ->   choose based on pick_larger_root
      if (x0 >= 1.0 and not pick_larger_root) {
        return x0;
      } else if (x1 >= 1.0) {
        return x1;
      } else {
        return std::optional<double>{};
      }
    } else {
      // For the inverse map, we want a scale_factor s such that 0 < s <= 1.
      // Note that gsl_poly_solve_quadratic returns x0 < x1.
      // So we have six cases:
      //  a) x0 < x1 <= 0        ->   no root
      //  b) x0 <= 0 < x1 <= 1   ->   Choose x1
      //  c) x0 <= 0 and x1 > 1  ->   no root
      //  d) 0 < x0 < x1 <= 1      ->   Choose according to pick_larger_root
      //  e) 0 < x0 <= 1 < x1      ->   Choose x0
      //  f) 1 < x0 < x1           ->   no root
      if (x0 <= 0.0) {
        if (x1 > 0.0 and x1 <= 1.0) {
          return x1;  // b)
        } else {
          return std::optional<double>{};  // a) and c)
        }
      } else if (x0 <= 1.0) {
        if (x1 > 1.0) {
          return x0;  // e)
        } else {
          return pick_larger_root ? x1 : x0;  // d)
        }
      } else {
        return std::optional<double>{};  // f)
      }
    }
  } else if (num_real_roots == 1) {
    if(solve_for_root_minus_one) {
      // We solved for scale_factor-1 above, so add 1 to get scale_factor.
      x0 += 1.0;
    }
    if (equal_within_roundoff(x0, 1.0)) {
      x0 = 1.0;
    }
    if (pick_root_greater_than_one) {
      if (x0 < 1.0) {
        return std::optional<double>{};
      } else {
        return x0;
      }
    } else {
      if (x0 <= 0.0 or x0 > 1.0) {
        return std::optional<double>{};
      } else {
        return x0;
      }
    }
  } else {
    return std::optional<double>{};
  }
}

template <typename T>
void d_scale_factor_d_src_point(
    const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>& result,
    const std::array<T, 3>& intersection_point,
    const std::array<double, 3>& proj_center,
    const std::array<double, 3>& sphere_center, const T& lambda) {
  using ReturnType = tt::remove_cvref_wrap_t<T>;
  ASSERT(not(intersection_point[0] == proj_center[0] and
             intersection_point[1] == proj_center[1] and
             intersection_point[2] == proj_center[2]),
         "d_scale_factor_d_src_point: trying to divide by zero.  If this"
         "happens, then the map is singular.");
  const ReturnType lambda_squared_over_denominator =
      square(lambda) / (square(intersection_point[0] - proj_center[0]) +
                        square(intersection_point[1] - proj_center[1]) +
                        square(intersection_point[2] - proj_center[2]) +
                        ((intersection_point[0] - proj_center[0]) *
                             (proj_center[0] - sphere_center[0]) +
                         (intersection_point[1] - proj_center[1]) *
                             (proj_center[1] - sphere_center[1]) +
                         (intersection_point[2] - proj_center[2]) *
                             (proj_center[2] - sphere_center[2])));
  for (size_t i = 0; i < 3; ++i) {
    gsl::at(*result, i) =
        lambda_squared_over_denominator *
        (gsl::at(sphere_center, i) - gsl::at(intersection_point, i));
  }
}

// Explicit instantiations
#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)

#define INSTANTIATE(_, data)                                              \
  template void scale_factor<DTYPE(data)>(                                \
      const gsl::not_null<tt::remove_cvref_wrap_t<DTYPE(data)>*>& result, \
      const std::array<DTYPE(data), 3>& src_point,                        \
      const std::array<double, 3>& proj_center,                           \
      const std::array<double, 3>& sphere_center, double radius,          \
      bool src_is_between_proj_and_target);                               \
  template void d_scale_factor_d_src_point<DTYPE(data)>(                  \
      const gsl::not_null<                                                \
          std::array<tt::remove_cvref_wrap_t<DTYPE(data)>, 3>*>& result,  \
      const std::array<DTYPE(data), 3>& intersection_point,               \
      const std::array<double, 3>& proj_center,                           \
      const std::array<double, 3>& sphere_center, const DTYPE(data) & lambda);

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

}  // namespace domain::CoordinateMaps::FocallyLiftedMapHelpers

1	// Distributed under the MIT License.
2	// See LICENSE.txt for details.
3
4	#include "Domain/CoordinateMaps/FocallyLiftedMapHelpers.hpp"
5
6	#include <cmath>
7	#include <gsl/gsl_poly.h>
8	#include <limits>
9
10	#include "DataStructures/DataVector.hpp"
11	#include "NumericalAlgorithms/RootFinding/QuadraticEquation.hpp"
12	#include "Utilities/ConstantExpressions.hpp"
13	#include "Utilities/DereferenceWrapper.hpp"
14	#include "Utilities/EqualWithinRoundoff.hpp"
15	#include "Utilities/ErrorHandling/Assert.hpp"
16	#include "Utilities/GenerateInstantiations.hpp"
17	#include "Utilities/MakeWithValue.hpp"
18
19	namespace domain::CoordinateMaps::FocallyLiftedMapHelpers {
20
21	template <typename T>
22	void scale_factor(const gsl::not_null<tt::remove_cvref_wrap_t<T>*>& result,	3,444✔
23	const std::array<T, 3>& src_point,
24	const std::array<double, 3>& proj_center,
25	const std::array<double, 3>& sphere_center, double radius,
26	const bool src_is_between_proj_and_target) {
27	using ReturnType = tt::remove_cvref_wrap_t<T>;
28	// quadratic equation is
29	// a scale_factor^2 + b scale_factor + c = 0
30	//
31	// For this case
32	// a = \|x_0^i-P^i\|^2
33	// b = 2(x_0^i-P^i)(P^j-C^j)\delta_{ij}
34	// c = \|P^i-C^i\|^2 - R^2
35	//
36	const ReturnType a = square(src_point[0] - proj_center[0]) +	4,956✔
37	square(src_point[1] - proj_center[1]) +	3,948✔
38	square(src_point[2] - proj_center[2]);	3,444✔
39	const ReturnType b =	3,948✔
40	2.0 *	3,444✔
41	((src_point[0] - proj_center[0]) * (proj_center[0] - sphere_center[0]) +	3,444✔
42	(src_point[1] - proj_center[1]) * (proj_center[1] - sphere_center[1]) +	3,948✔
43	(src_point[2] - proj_center[2]) * (proj_center[2] - sphere_center[2]));	3,444✔
44	const double c = square(sphere_center[0] - proj_center[0]) +	3,444✔
45	square(sphere_center[1] - proj_center[1]) +	3,444✔
46	square(sphere_center[2] - proj_center[2]) - square(radius);	6,888✔
47	if (src_is_between_proj_and_target) {	3,444✔
48	// Here we assume that src_point is between proj_center and
49	// target_point. There are three cases: 1) src and proj are both
50	// inside the sphere, 2) src and proj are both outside the sphere,
51	// and 3) proj is outside the sphere and src is inside the sphere.
52	// Note that a root of zero corresponds to target_point = proj_center,
53	// and a root of unity corresponds to target_point = src_point.
54	// To cover all 3 cases, we choose the smallest root that is
55	// greater than or equal to unity. This means for case 2) we are
56	// choosing the point closest to src.
57	*result = smallest_root_greater_than_value_within_roundoff(	1,476✔
58	a, b, make_with_value<ReturnType>(a, c), 1.0);	2,520✔
59	} else {
60	// Here we assume that target_point is between proj_center and
61	// src_point. There are three cases: 1) proj is inside the sphere
62	// and src is outside the sphere 2) src is inside the sphere and proj
63	// is outside the sphere, and 3) the sphere is between src and proj.
64	// Note that a root of zero corresponds to target_point = proj_center,
65	// and a root of unity corresponds to target_point = src_point.
66	// To cover all 3 cases, we choose the largest root that is less than
67	// or equal to unity, and we require that this root is positive.
68	// This means that for 3) we are choosing the point closest to src.
69	*result = largest_root_between_values_within_roundoff(	1,968✔
70	a, b, make_with_value<ReturnType>(a, c), 0.0, 1.0);	3,360✔
71	}
72	}	3,444✔
73
74	std::optional<double> try_scale_factor(	144✔
75	const std::array<double, 3>& src_point,
76	const std::array<double, 3>& proj_center,
77	const std::array<double, 3>& sphere_center, double radius,
78	const bool pick_larger_root, const bool pick_root_greater_than_one,
79	const bool solve_for_root_minus_one) {
80	double x0 = std::numeric_limits<double>::signaling_NaN();	144✔
81	double x1 = std::numeric_limits<double>::signaling_NaN();	144✔
82	int num_real_roots = 0;	144✔
83	if (solve_for_root_minus_one) {	144✔
84	// We solve the quadratic for (scale_factor-1) instead of
85	// scale_factor to avoid roundoff problems when scale_factor is very
86	// nearly equal to unity. Note that scale_factor==1 will occur in
87	// the inverse map when src_point (which is x^i) is
88	// on the sphere.
89	//
90	// Roundoff is not a problem when forward-mapping and solving only
91	// for lambda.
92
93	// If the original quadratic equation is
94	//
95	// A scale_factor^2 + B scale_factor + C = 0,
96	//
97	// and if x = scale_factor-1, then the quadratic equation for x is
98	//
99	// a x^2 + b x + c = 0
100	//
101	// where
102	//
103	// a = A
104	// b = 2A + B
105	// c = A + B + C
106	//
107	// For this case
108	// A = \|x_0^i-P^i\|^2
109	// B = 2(x_0^i-P^i)(P^j-C^j)\delta_{ij}
110	// C = \|P^i-C^i\|^2 - R^2
111	//
112	// Note that A + B + C is \|x_0^i-P^i+P^i-C^i\|^2 - R^2
113	// and 2A+B is 2(x_0^i-P^i)(x_0^j-P^j+P^j-C^j), so
114	//
115	// a = \|x_0^i-P^i\|^2
116	// b = 2(x_0^i-P^i)(x_0^j-C^j)\delta_{ij}
117	// c = \|x_0^i-C^i\|^2 - R^2
118
119	const double a = square(src_point[0] - proj_center[0]) +	108✔
120	square(src_point[1] - proj_center[1]) +	108✔
121	square(src_point[2] - proj_center[2]);	108✔
122	const double b =
123	2.0 *
124	((src_point[0] - proj_center[0]) * (src_point[0] - sphere_center[0]) +	108✔
125	(src_point[1] - proj_center[1]) * (src_point[1] - sphere_center[1]) +	108✔
126	(src_point[2] - proj_center[2]) * (src_point[2] - sphere_center[2]));	108✔
127	const double c = square(sphere_center[0] - src_point[0]) +	108✔
128	square(sphere_center[1] - src_point[1]) +	108✔
129	square(sphere_center[2] - src_point[2]) - square(radius);	216✔
130	num_real_roots = gsl_poly_solve_quadratic(a, b, c, &x0, &x1);	108✔
131	} else {
132	const double a = square(src_point[0] - proj_center[0]) +	36✔
133	square(src_point[1] - proj_center[1]) +	36✔
134	square(src_point[2] - proj_center[2]);	36✔
135	const double b =
136	2.0 *
137	((src_point[0] - proj_center[0]) * (proj_center[0] - sphere_center[0]) +	36✔
138	(src_point[1] - proj_center[1]) * (proj_center[1] - sphere_center[1]) +	36✔
139	(src_point[2] - proj_center[2]) * (proj_center[2] - sphere_center[2]));	36✔
140	const double c = square(proj_center[0] - sphere_center[0]) +	36✔
141	square(proj_center[1] - sphere_center[1]) +	36✔
142	square(proj_center[2] - sphere_center[2]) - square(radius);	72✔
143	num_real_roots = gsl_poly_solve_quadratic(a, b, c, &x0, &x1);	36✔
144	}
145	if (num_real_roots == 2) {	144✔
146	if(solve_for_root_minus_one) {	144✔
147	// We solved for scale_factor-1 above, so add 1 to get scale_factor.
148	x0 += 1.0;	108✔
149	x1 += 1.0;	108✔
150	}
151	if (UNLIKELY(equal_within_roundoff(x0, 1.0))) {	288✔
152	x0 = 1.0;	×
153	}
154	if (UNLIKELY(equal_within_roundoff(x1, 1.0))) {	288✔
155	x1 = 1.0;	42✔
156	}
157	if (pick_root_greater_than_one) {	144✔
158	// For the inverse map, we want the a scale_factor s such that
159	// s >= 1. Note that gsl_poly_solve_quadratic returns x0 < x1.
160	// have three cases:
161	// a) x0 < x1 < 1 -> no root
162	// b) x0 < 1 <= x1 -> Choose x1
163	// c) 1 <= x0 < x1 -> choose based on pick_larger_root
164	if (x0 >= 1.0 and not pick_larger_root) {	60✔
165	return x0;	×
166	} else if (x1 >= 1.0) {	60✔
167	return x1;	60✔
168	} else {
169	return std::optional<double>{};	×
170	}
171	} else {
172	// For the inverse map, we want a scale_factor s such that 0 < s <= 1.
173	// Note that gsl_poly_solve_quadratic returns x0 < x1.
174	// So we have six cases:
175	// a) x0 < x1 <= 0 -> no root
176	// b) x0 <= 0 < x1 <= 1 -> Choose x1
177	// c) x0 <= 0 and x1 > 1 -> no root
178	// d) 0 < x0 < x1 <= 1 -> Choose according to pick_larger_root
179	// e) 0 < x0 <= 1 < x1 -> Choose x0
180	// f) 1 < x0 < x1 -> no root
181	if (x0 <= 0.0) {	84✔
182	if (x1 > 0.0 and x1 <= 1.0) {	84✔
183	return x1; // b)	84✔
184	} else {
185	return std::optional<double>{}; // a) and c)	×
186	}
187	} else if (x0 <= 1.0) {	×
188	if (x1 > 1.0) {	×
189	return x0; // e)	×
190	} else {
191	return pick_larger_root ? x1 : x0; // d)	×
192	}
193	} else {
194	return std::optional<double>{}; // f)	×
195	}
196	}
197	} else if (num_real_roots == 1) {	×
198	if(solve_for_root_minus_one) {	×
199	// We solved for scale_factor-1 above, so add 1 to get scale_factor.
200	x0 += 1.0;	×
201	}
202	if (equal_within_roundoff(x0, 1.0)) {	×
203	x0 = 1.0;	×
204	}
205	if (pick_root_greater_than_one) {	×
206	if (x0 < 1.0) {	×
207	return std::optional<double>{};	×
208	} else {
209	return x0;	×
210	}
211	} else {
212	if (x0 <= 0.0 or x0 > 1.0) {	×
213	return std::optional<double>{};	×
214	} else {
215	return x0;	×
216	}
217	}
218	} else {
219	return std::optional<double>{};	×
220	}
221	}
222
223	template <typename T>
224	void d_scale_factor_d_src_point(	1,748✔
225	const gsl::not_null<std::array<tt::remove_cvref_wrap_t<T>, 3>*>& result,
226	const std::array<T, 3>& intersection_point,
227	const std::array<double, 3>& proj_center,
228	const std::array<double, 3>& sphere_center, const T& lambda) {
229	using ReturnType = tt::remove_cvref_wrap_t<T>;
230	ASSERT(not(intersection_point[0] == proj_center[0] and	1,748✔
231	intersection_point[1] == proj_center[1] and
232	intersection_point[2] == proj_center[2]),
233	"d_scale_factor_d_src_point: trying to divide by zero. If this"
234	"happens, then the map is singular.");
235	const ReturnType lambda_squared_over_denominator =	2,204✔
236	square(lambda) / (square(intersection_point[0] - proj_center[0]) +	2,660✔
237	square(intersection_point[1] - proj_center[1]) +	2,204✔
238	square(intersection_point[2] - proj_center[2]) +	2,204✔
239	((intersection_point[0] - proj_center[0]) *	1,748✔
240	(proj_center[0] - sphere_center[0]) +	2,204✔
241	(intersection_point[1] - proj_center[1]) *	1,748✔
242	(proj_center[1] - sphere_center[1]) +	2,204✔
243	(intersection_point[2] - proj_center[2]) *	1,748✔
244	(proj_center[2] - sphere_center[2])));	1,748✔
245	for (size_t i = 0; i < 3; ++i) {	6,992✔
246	gsl::at(*result, i) =	10,488✔
247	lambda_squared_over_denominator *	5,244✔
248	(gsl::at(sphere_center, i) - gsl::at(intersection_point, i));	11,856✔
249	}
250	}	1,748✔
251
252	// Explicit instantiations
253	#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
254
255	#define INSTANTIATE(_, data) \
256	template void scale_factor<DTYPE(data)>( \
257	const gsl::not_null<tt::remove_cvref_wrap_t<DTYPE(data)>*>& result, \
258	const std::array<DTYPE(data), 3>& src_point, \
259	const std::array<double, 3>& proj_center, \
260	const std::array<double, 3>& sphere_center, double radius, \
261	bool src_is_between_proj_and_target); \
262	template void d_scale_factor_d_src_point<DTYPE(data)>( \
263	const gsl::not_null< \
264	std::array<tt::remove_cvref_wrap_t<DTYPE(data)>, 3>*>& result, \
265	const std::array<DTYPE(data), 3>& intersection_point, \
266	const std::array<double, 3>& proj_center, \
267	const std::array<double, 3>& sphere_center, const DTYPE(data) & lambda);
268
269	GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector,
270	std::reference_wrapper<const double>,
271	std::reference_wrapper<const DataVector>))
272	#undef INSTANTIATE
273	#undef DTYPE
274
275	} // namespace domain::CoordinateMaps::FocallyLiftedMapHelpers

sxs-collaboration / spectre / 4138821718

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