forked from celeritas-project/celeritas
-
Notifications
You must be signed in to change notification settings - Fork 0
/
QuadricCylConverter.hh
131 lines (116 loc) · 3.76 KB
/
QuadricCylConverter.hh
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
//----------------------------------*-C++-*----------------------------------//
// Copyright 2023-2024 UT-Battelle, LLC, and other Celeritas developers.
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: (Apache-2.0 OR MIT)
//---------------------------------------------------------------------------//
//! \file orange/surf/detail/QuadricCylConverter.hh
//---------------------------------------------------------------------------//
#pragma once
#include <optional>
#include "corecel/math/SoftEqual.hh"
#include "../CylAligned.hh"
#include "../SimpleQuadric.hh"
namespace celeritas
{
namespace detail
{
//---------------------------------------------------------------------------//
/*!
* Try to convert a simple quadric to a cylinder.
*
* The simple quadric must have already undergone normalization (one
* second-order term approximately zero, the others positive).
*/
class QuadricCylConverter
{
public:
// Construct with tolerance
inline QuadricCylConverter(real_type tol);
// Try converting to a cylinder with this orientation
template<Axis T>
std::optional<CylAligned<T>>
operator()(AxisTag<T>, SimpleQuadric const& sq) const;
private:
SoftEqual<> soft_equal_;
};
//---------------------------------------------------------------------------//
// INLINE DEFINITIONS
//---------------------------------------------------------------------------//
/*!
* Construct with tolerance.
*/
QuadricCylConverter::QuadricCylConverter(real_type tol) : soft_equal_{tol} {}
//---------------------------------------------------------------------------//
/*!
* Try converting to a cylinder with this orientation.
*
* Cone:
* \verbatim
(x - x_0)^2 + (y - y_0)^2 - r^2 = 0
* \endverbatim
* Expanded:
* \verbatim
x^2 + y^2
- 2x_0 x - 2y_0 y
+ (x_0^2 + y_0^2 - r^2) = 0
* \endverbatim
* SQ:
* \verbatim
a x^2 + (a + \epsilon) y^2 + e x + f y + h = 0
* \endverbatim
* Normalized (\c epsilon = 0)
* \verbatim
x^2 + y^2 + e/b x + f/b y + g/b z + h/b = 0
* \endverbatim
* Match terms:
* \verbatim
-2 x_0 = e/b --> x_0 = e / (-2 * b)
-2 y_0 = f/b --> y_0 = f / (-2 * b)
(x_0^2 + y_0^2 - r^2) = h/b
-> r^2 = x_0^2 + y_0^2 - h/b
* \endverbatim
*/
template<Axis T>
std::optional<CylAligned<T>>
QuadricCylConverter::operator()(AxisTag<T>, SimpleQuadric const& sq) const
{
// Other coordinate system
constexpr auto U = CylAligned<T>::u_axis();
constexpr auto V = CylAligned<T>::v_axis();
auto second = sq.second();
if (!soft_equal_(0, second[to_int(T)]))
{
// Not the zero component we're looking for
return {};
}
if (!soft_equal_(0, sq.first()[to_int(T)]))
{
// Not an axis-aligned cylinder
return {};
}
if (!soft_equal_(second[to_int(U)], second[to_int(V)]))
{
// Not a *circular* cylinder
return {};
}
CELER_ASSERT(second[to_int(U)] > 0 && second[to_int(V)] > 0);
// Normalize so U, V second-order coefficients are 1
auto const inv_norm = 2 / (second[to_int(U)] + second[to_int(V)]);
// Calculate origin from first-order coefficients
Real3 origin{0, 0, 0};
origin[to_int(U)] = real_type{-0.5} * inv_norm * sq.first()[to_int(U)];
origin[to_int(V)] = real_type{-0.5} * inv_norm * sq.first()[to_int(V)];
real_type radius_sq = ipow<2>(origin[to_int(U)])
+ ipow<2>(origin[to_int(V)]) - sq.zeroth() * inv_norm;
if (radius_sq <= 0)
{
// No real solution
return {};
}
// Clear potential signed zeros before returning
origin += real_type{0};
return CylAligned<T>::from_radius_sq(origin, radius_sq);
}
//---------------------------------------------------------------------------//
} // namespace detail
} // namespace celeritas