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
132
133
134
135
136
137
138
139
|
// This file is part of VecGeom and is distributed under the
// conditions in the file LICENSE.txt in the top directory.
// For the full list of authors see CONTRIBUTORS.txt and `git log`.
/// Unit test for Paralleliped.
/// @file test/unit_tests/TestParaboloid.cpp
/// @author Raman Sehgal
#undef NDEBUG
#include "VecGeom/base/FpeEnable.h"
#include "VecGeom/base/Global.h"
#include "VecGeom/base/Vector3D.h"
#include "VecGeom/volumes/Box.h"
#include "VecGeom/volumes/Paraboloid.h"
#include "ApproxEqual.h"
#include <cmath>
#include <iomanip>
#define PI 3.14159265358979323846
template <class Paraboloid_t, class Vec_t = vecgeom::Vector3D<vecgeom::Precision>>
bool TestParaboloid()
{
Paraboloid_t p1("testPara", 6., 9., 10.);
// std::cout<< p1.GetK1() << std::endl;
// assert(false);
std::cout << std::setprecision(15);
Vec_t norm(0., 0., 0.);
double Dist = 0;
Vec_t pzero(0., 0., 0.), pbigx(100, 0., 0.), pbigmx(-100., 0., 0.), pbigy(0., 100., 0.), pbigmy(0., -100., 0.),
pbigz(0., 0., 100.), pbigmz(0., 0., -100);
Vec_t dirx(1., 0., 0.), dirmx(-1., 0., 0.), diry(0., 1., 0.), dirmy(0., -1., 0.), dirz(0., 0., 1), dirmz(0., 0., -1.);
Vec_t zSurfPt(0., 0., 10.), mzSurfPt(0., 0., -10.);
// Inside
assert(p1.Inside(pzero) == vecgeom::EInside::kInside);
assert(p1.Inside(pbigx) == vecgeom::EInside::kOutside);
assert(p1.Inside(pbigmx) == vecgeom::EInside::kOutside);
assert(p1.Inside(pbigy) == vecgeom::EInside::kOutside);
assert(p1.Inside(pbigmy) == vecgeom::EInside::kOutside);
assert(p1.Inside(pbigz) == vecgeom::EInside::kOutside);
assert(p1.Inside(pbigmz) == vecgeom::EInside::kOutside);
assert(p1.Inside(zSurfPt) == vecgeom::EInside::kSurface);
assert(p1.Inside(mzSurfPt) == vecgeom::EInside::kSurface);
Vec_t pInsideZ(0., 0., 10. - 2 * vecgeom::kTolerance);
assert(p1.Inside(pInsideZ) == vecgeom::EInside::kInside);
Vec_t pOutsideZ(0., 0., 10. + 2 * vecgeom::kTolerance);
assert(p1.Inside(pOutsideZ) == vecgeom::EInside::kOutside);
Vec_t pWithinZTolerance(0., 0., 10. + 0.25 * vecgeom::kTolerance);
assert(p1.Inside(pWithinZTolerance) == vecgeom::EInside::kSurface);
Dist = p1.DistanceToOut(pzero, dirz);
assert(Dist == 10.);
Dist = p1.DistanceToOut(pzero, dirmz);
assert(Dist == 10.);
Dist = p1.DistanceToOut(zSurfPt, dirz);
assert(Dist == 0.);
Dist = p1.DistanceToOut(mzSurfPt, dirmz);
assert(Dist == 0.);
Dist = p1.DistanceToIn(pbigz, dirmz);
assert(Dist == 90.);
Dist = p1.DistanceToIn(pbigmz, dirz);
assert(Dist == 90.);
Dist = p1.DistanceToIn(zSurfPt, dirmz);
assert(Dist == 0.);
Dist = p1.DistanceToIn(mzSurfPt, dirz);
assert(Dist == 0.);
Paraboloid_t p2("testPara", 0., 8., 10.);
Dist = p2.DistanceToOut(Vec_t(0., 0., -10), dirmz);
assert(Dist == 0.);
Dist = p2.DistanceToOut(Vec_t(0., 0., -10), dirz);
assert(Dist == 20.);
Dist = p2.DistanceToIn(Vec_t(0., 0., -10), dirz);
assert(Dist == 0.);
Vec_t tmpDir = Vec_t(1., 2., 2.).Unit();
Dist = p2.DistanceToOut(pzero, tmpDir);
Vec_t tmpPt = pzero + Dist * tmpDir;
Dist = p2.DistanceToOut(tmpPt, tmpDir);
assert(Dist == 0.);
Dist = p2.DistanceToIn(tmpPt, -tmpDir);
assert(Dist == 0.);
Vec_t normal(0., 0., 0.);
Vec_t pTopZ(0., 0., 10), pBottomZ(0., 0., -10.);
bool valid = p1.Normal(pTopZ, normal);
assert(valid && ApproxEqual(normal, Vec_t(0, 0, 1.)));
valid = p1.Normal(pBottomZ, normal);
assert(valid && ApproxEqual(normal, Vec_t(0, 0, -1.)));
pTopZ.Set(3., 4, 10.);
valid = p1.Normal(pTopZ, normal);
assert(valid && ApproxEqual(normal, Vec_t(0, 0, 1.)));
pBottomZ.Set(3., 4, -10.);
valid = p1.Normal(pBottomZ, normal);
assert(valid && ApproxEqual(normal, Vec_t(0, 0, -1.)));
pTopZ.Set(9., 0., 10.);
valid = p1.Normal(pTopZ, normal);
assert(valid);
// Testing noraml for outside point,
// Logic to calculate normal for this case is not yet written.
// This test is just to check the "valid", which should be false
pTopZ.Set(9., 0., 5.);
valid = p1.Normal(pTopZ, normal);
assert(!valid);
Vec_t pOutZ(0., 0., 11.);
valid = p1.Normal(pOutZ, normal);
assert(!valid && ApproxEqual(normal, Vec_t(0, 0, 1.)));
pOutZ.Set(0., 0., -11.);
valid = p1.Normal(pOutZ, normal);
assert(!valid && ApproxEqual(normal, Vec_t(0, 0, -1.)));
// Check Extent and cached BBox
Vec_t minExtent, maxExtent;
Vec_t minBBox, maxBBox;
p1.Extent(minExtent, maxExtent);
p1.GetUnplacedVolume()->GetBBox(minBBox, maxBBox);
assert(ApproxEqual(minExtent, minBBox));
assert(ApproxEqual(maxExtent, maxBBox));
return true;
}
int main(int argc, char *argv[])
{
assert(TestParaboloid<vecgeom::SimpleParaboloid>());
std::cout << "VecGeomParaboloid passed\n";
return 0;
}
|
