Mantid
Loading...
Searching...
No Matches
Cylinder.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
11
12#include <limits>
13
14#ifdef ENABLE_OPENCASCADE
15// Opencascade defines _USE_MATH_DEFINES without checking whether it is already
16// used.
17// Undefine it here before we include the headers to avoid a warning
18#ifdef _MSC_VER
19#undef _USE_MATH_DEFINES
20#ifdef M_SQRT1_2
21#undef M_SQRT1_2
22#endif
23#endif
24
26GNU_DIAG_OFF("conversion")
27GNU_DIAG_OFF("cast-qual")
28#include <BRepPrimAPI_MakeCylinder.hxx>
29GNU_DIAG_ON("conversion")
30GNU_DIAG_ON("cast-qual")
31#endif
32
33namespace Mantid::Geometry {
34namespace {
35Kernel::Logger logger("Cylinder");
36}
38using Kernel::V3D;
39
41 : Quadratic(), m_centre(), m_normal(1, 0, 0), m_normVec(0), m_radius(0.0)
46{
47 // Called after it has been sized by Quadratic
49}
50
56{
57 return new Cylinder(*this);
58}
59
60std::unique_ptr<Cylinder> Cylinder::clone() const
65{
66 return std::unique_ptr<Cylinder>(doClone());
67}
68
69int Cylinder::setSurface(const std::string &Pstr)
79{
80 enum { errDesc = -1, errAxis = -2, errCent = -3, errRadius = -4 };
81
82 std::string Line = Pstr;
83 std::string item;
84 if (!Mantid::Kernel::Strings::section(Line, item) || tolower(item[0]) != 'c' || item.length() < 2 ||
85 item.length() > 3)
86 return errDesc;
87
88 // Cylinders on X/Y/Z axis
89 const std::size_t itemPt((item[1] == '/' && item.length() == 3) ? 2 : 1);
90 const auto ptype = static_cast<std::size_t>(tolower(item[itemPt]) - 'x');
91 if (ptype >= 3)
92 return errAxis;
93 std::vector<double> norm(3, 0.0);
94 std::vector<double> cent(3, 0.0);
95 norm[ptype] = 1.0;
96
97 if (itemPt != 1) {
98 // get the other two coordinates
99 std::size_t index((!ptype) ? 1 : 0);
100 while (index < 3 && Mantid::Kernel::Strings::section(Line, cent[index])) {
101 index++;
102 if (index == ptype)
103 index++;
104 }
105 if (index != 3)
106 return errCent;
107 }
108 // Now get radius
109 double R;
110 if (!Mantid::Kernel::Strings::section(Line, R) || R <= 0.0)
111 return errRadius;
112
113 m_centre = Kernel::V3D(cent[0], cent[1], cent[2]);
114 m_normal = Kernel::V3D(norm[0], norm[1], norm[2]);
115 m_normVec = ptype + 1;
117 setBaseEqn();
118 return 0;
119}
120
121int Cylinder::side(const Kernel::V3D &Pt) const
130{
131 if (m_normVec) // m_normVec =1-3 (point to exclude == m_normVec-1)
132 {
133 if (m_radius > 0.0) {
134 double x = Pt[m_normVec % 3] - m_centre[m_normVec % 3];
135 x *= x;
136 double y = Pt[(m_normVec + 1) % 3] - m_centre[(m_normVec + 1) % 3];
137 y *= y;
138 double displace = x + y - m_radius * m_radius;
139 if (fabs(displace * m_oneoverradius) < Tolerance)
140 return 0;
141 return (displace > 0.0) ? 1 : -1;
142 } else {
143 return -1;
144 }
145 }
146 return Quadratic::side(Pt);
147}
148
149bool Cylinder::onSurface(const Kernel::V3D &Pt) const
154{
155 if (m_normVec > 0) // m_normVec =1-3 (point to exclude == m_normVec-1)
156 {
157 double x = Pt[m_normVec % 3] - m_centre[m_normVec % 3];
158 x *= x;
159 double y = Pt[(m_normVec + 1) % 3] - m_centre[(m_normVec + 1) % 3];
160 y *= y;
161 return (std::abs((x + y) - m_radius * m_radius) <= Tolerance);
162 }
163 return Quadratic::onSurface(Pt);
164}
165
172{
173 m_normVec = 0;
174 for (std::size_t i = 0; i < 3; i++) {
175 if (fabs(m_normal[i]) > (1.0 - Tolerance)) {
176 m_normVec = i + 1;
177 return;
178 }
179 }
180}
181
188{
189 m_centre.rotate(MA);
190 m_normal.rotate(MA);
192 setNormVec();
194}
195
201{
202 if (m_normVec > 0) {
203 m_centre[m_normVec % 3] += Pt[m_normVec % 3];
204 m_centre[(m_normVec + 1) % 3] += Pt[(m_normVec + 1) % 3];
205 } else
206 m_centre += Pt;
208}
209
215{
216 m_centre = A;
217 setBaseEqn();
218}
219
226{
227 m_normal = normalize(A);
228 setBaseEqn();
229 setNormVec();
230}
231
237{
238 const double CdotN(m_centre.scalar_prod(m_normal));
239 BaseEqn[0] = 1.0 - m_normal[0] * m_normal[0]; // A x^2
240 BaseEqn[1] = 1.0 - m_normal[1] * m_normal[1]; // B y^2
241 BaseEqn[2] = 1.0 - m_normal[2] * m_normal[2]; // C z^2
242 BaseEqn[3] = -2 * m_normal[0] * m_normal[1]; // D xy
243 BaseEqn[4] = -2 * m_normal[0] * m_normal[2]; // E xz
244 BaseEqn[5] = -2 * m_normal[1] * m_normal[2]; // F yz
245 BaseEqn[6] = 2.0 * (m_normal[0] * CdotN - m_centre[0]); // G x
246 BaseEqn[7] = 2.0 * (m_normal[1] * CdotN - m_centre[1]); // H y
247 BaseEqn[8] = 2.0 * (m_normal[2] * CdotN - m_centre[2]); // J z
248 BaseEqn[9] = m_centre.scalar_prod(m_centre) - CdotN * CdotN - m_radius * m_radius; // K const
249}
250
251double Cylinder::distance(const Kernel::V3D &A) const
261{
262 // First find the normal going to the point
263 const Kernel::V3D Amov = A - m_centre;
264 double lambda = Amov.scalar_prod(m_normal);
265 const Kernel::V3D Ccut = m_normal * lambda;
266 // The distance is from the centre line to the
267 return fabs(Ccut.distance(Amov) - m_radius);
268}
269
270void Cylinder::write(std::ostream &OX) const
275{
276 // -3 -2 -1 0 1 2 3
277 const char Tailends[] = "zyx xyz";
278 const int Ndir = m_normal.masterDir(Tolerance);
279 if (Ndir == 0) {
280 // general surface
282 return;
283 }
284
285 const int Cdir = m_centre.masterDir(Tolerance);
286 std::ostringstream cx;
287
288 writeHeader(cx);
289 cx.precision(Surface::Nprecision);
290 // Name and transform
291
292 if (Cdir * Cdir == Ndir * Ndir || m_centre.nullVector(Tolerance)) {
293 cx << "c";
294 cx << Tailends[Ndir + 3] << " "; // set x,y,z based on Ndir
295 cx << m_radius;
296 } else {
297 cx << " c/";
298 cx << Tailends[Ndir + 3] << " "; // set x,y,z based on Ndir
299
300 if (Ndir == 1 || Ndir == -1)
301 cx << m_centre[1] << " " << m_centre[2] << " ";
302 else if (Ndir == 2 || Ndir == -2)
303 cx << m_centre[0] << " " << m_centre[2] << " ";
304 else
305 cx << m_centre[0] << " " << m_centre[1] << " ";
306
307 cx << m_radius;
308 }
309
311}
312
313double Cylinder::lineIntersect(const Kernel::V3D &Pt, const Kernel::V3D &uVec) const
322{
323 (void)Pt; // Avoid compiler warning
324 (void)uVec; // Avoid compiler warning
325 return -1;
326}
327
328void Cylinder::print() const
332{
334 logger.debug() << "Axis ==" << m_normal << " ";
335 logger.debug() << "Centre == " << m_centre << " ";
336 logger.debug() << "Radius == " << m_radius << '\n';
337}
338
339void Cylinder::getBoundingBox(double &xmax, double &ymax, double &zmax, double &xmin, double &ymin, double &zmin) {
355 std::vector<V3D> listOfPoints;
356 double txmax, tymax, tzmax, txmin, tymin, tzmin;
357 txmax = xmax;
358 tymax = ymax;
359 tzmax = zmax;
360 txmin = xmin;
361 tymin = ymin;
362 tzmin = zmin;
363 V3D xminPoint, xmaxPoint, yminPoint, ymaxPoint, zminPoint, zmaxPoint;
364 // xmin and max plane
365 if (m_normal[0] != 0) {
366 xminPoint = m_centre + m_normal * ((xmin - m_centre[0]) / m_normal[0]);
367 xmaxPoint = m_centre + m_normal * ((xmax - m_centre[0]) / m_normal[0]);
368 listOfPoints.emplace_back(xminPoint);
369 listOfPoints.emplace_back(xmaxPoint);
370 }
371
372 if (m_normal[1] != 0) {
373 // ymin plane
374 yminPoint = m_centre + m_normal * ((ymin - m_centre[1]) / m_normal[1]);
375 // ymax plane
376 ymaxPoint = m_centre + m_normal * ((ymax - m_centre[1]) / m_normal[1]);
377 listOfPoints.emplace_back(yminPoint);
378 listOfPoints.emplace_back(ymaxPoint);
379 }
380 if (m_normal[2] != 0) {
381 // zmin plane
382 zminPoint = m_centre + m_normal * ((zmin - m_centre[2]) / m_normal[2]);
383 // zmax plane
384 zmaxPoint = m_centre + m_normal * ((zmax - m_centre[2]) / m_normal[2]);
385 listOfPoints.emplace_back(zminPoint);
386 listOfPoints.emplace_back(zmaxPoint);
387 }
388 if (!listOfPoints.empty()) {
389 xmin = ymin = zmin = std::numeric_limits<double>::max();
390 xmax = ymax = zmax = std::numeric_limits<double>::lowest();
391 for (std::vector<V3D>::const_iterator it = listOfPoints.begin(); it != listOfPoints.end(); ++it) {
392 // std::cout<<(*it)<<'\n';
393 if ((*it)[0] < xmin)
394 xmin = (*it)[0];
395 if ((*it)[1] < ymin)
396 ymin = (*it)[1];
397 if ((*it)[2] < zmin)
398 zmin = (*it)[2];
399 if ((*it)[0] > xmax)
400 xmax = (*it)[0];
401 if ((*it)[1] > ymax)
402 ymax = (*it)[1];
403 if ((*it)[2] > zmax)
404 zmax = (*it)[2];
405 }
406 xmax += m_radius;
407 ymax += m_radius;
408 zmax += m_radius;
409 xmin -= m_radius;
410 ymin -= m_radius;
411 zmin -= m_radius;
412 if (xmax > txmax)
413 xmax = txmax;
414 if (xmin < txmin)
415 xmin = txmin;
416 if (ymax > tymax)
417 ymax = tymax;
418 if (ymin < tymin)
419 ymin = tymin;
420 if (zmax > tzmax)
421 zmax = tzmax;
422 if (zmin < tzmin)
423 zmin = tzmin;
424 }
425}
426
427#ifdef ENABLE_OPENCASCADE
428TopoDS_Shape Cylinder::createShape() {
429 gp_Pnt center;
430 center.SetX(m_centre[0] - m_normal[0] * 500.0);
431 center.SetY(m_centre[1] - m_normal[1] * 500.0);
432 center.SetZ(m_centre[2] - m_normal[2] * 500.0);
433 gp_Ax2 gpA(center, gp_Dir(m_normal[0], m_normal[1], m_normal[2]));
434 return BRepPrimAPI_MakeCylinder(gpA, m_radius, 1000.0, 2.0 * M_PI).Solid();
435}
436#endif
437} // namespace Mantid::Geometry
const std::vector< double > * lambda
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define fabs(x)
Definition: Matrix.cpp:22
#define GNU_DIAG_ON(x)
#define GNU_DIAG_OFF(x)
This is a collection of macros for turning compiler warnings off in a controlled manner.
Holds a cylinder as a vector form.
Definition: Cylinder.h:32
bool onSurface(const Kernel::V3D &) const override
Calculate if the point PT on the cylinder.
Definition: Cylinder.cpp:149
void setNormVec()
check to obtain orientation
Definition: Cylinder.cpp:166
int side(const Kernel::V3D &) const override
Calculate if the point PT within the middle of the cylinder.
Definition: Cylinder.cpp:121
void getBoundingBox(double &xmax, double &ymax, double &zmax, double &xmin, double &ymin, double &zmin) override
bounding box for the surface
Definition: Cylinder.cpp:339
std::unique_ptr< Cylinder > clone() const
Makes a clone (implicit virtual copy constructor)
Definition: Cylinder.cpp:60
Cylinder * doClone() const override
Makes a clone (implicit virtual copy constructor)
Definition: Cylinder.cpp:51
void print() const override
Debug routine to print out basic information.
Definition: Cylinder.cpp:328
double m_radius
Radius of cylinder.
Definition: Cylinder.h:40
void setRadiusInternal(const double &r)
Definition: Cylinder.h:46
void displace(const Kernel::V3D &) override
Apply a displacement Pt.
Definition: Cylinder.cpp:196
void setCentre(const Kernel::V3D &)
Sets the centre Kernel::V3D.
Definition: Cylinder.cpp:210
int setSurface(const std::string &) override
Processes a standard MCNPX cone string Recall that cones can only be specified on an axis Valid input...
Definition: Cylinder.cpp:69
void setNorm(const Kernel::V3D &)
Sets the centre line unit vector A does not need to be a unit vector.
Definition: Cylinder.cpp:220
Kernel::V3D m_normal
Direction of centre line.
Definition: Cylinder.h:37
double distance(const Kernel::V3D &) const override
Calculates the distance of point A from the surface of the cylinder.
Definition: Cylinder.cpp:251
std::size_t m_normVec
Normal vector is x,y or z :: (1-3) (0 if general)
Definition: Cylinder.h:38
void rotate(const Kernel::Matrix< double > &) override
Apply a rotation to the cylinder and re-check the status of the main axis.
Definition: Cylinder.cpp:182
virtual double lineIntersect(const Kernel::V3D &, const Kernel::V3D &) const
Given a track starting from Pt and traveling along uVec determine the intersection point (distance)
Definition: Cylinder.cpp:313
Kernel::V3D m_centre
Kernel::V3D for centre.
Definition: Cylinder.h:36
Cylinder()
Standard Constructor creats a cylinder (radius 0) along the x axis.
Definition: Cylinder.cpp:40
void setBaseEqn() override
Sets an equation of type (cylinder)
Definition: Cylinder.cpp:232
void write(std::ostream &) const override
Write out the cylinder for MCNPX.
Definition: Cylinder.cpp:270
Impliments a line.
Definition: Line.h:43
Holds a basic quadratic surface.
Definition: Quadratic.h:29
bool onSurface(const Kernel::V3D &) const override
is point valid on surface
Definition: Quadratic.cpp:233
std::vector< double > BaseEqn
Base equation (as a 10 point vector)
Definition: Quadratic.h:35
void displace(const Kernel::V3D &) override
Apply a general displacement to the surface.
Definition: Quadratic.cpp:244
void rotate(const Kernel::Matrix< double > &) override
Rotate the surface by matrix MX.
Definition: Quadratic.cpp:258
void print() const override
Print out the genreal equation for debugging.
Definition: Quadratic.cpp:298
void write(std::ostream &) const override
Writes out an MCNPX surface description Note : Swap since base equation is swapped in gq output (mcnp...
Definition: Quadratic.cpp:310
int side(const Kernel::V3D &) const override
Determine if the the Point is true to the surface or on the other side.
Definition: Quadratic.cpp:68
static const int Nprecision
Precision of the output.
Definition: Surface.h:41
Numerical Matrix class.
Definition: Matrix.h:42
Class for 3D vectors.
Definition: V3D.h:34
double distance(const V3D &v) const noexcept
Calculates the distance between two vectors.
Definition: V3D.h:287
constexpr double scalar_prod(const V3D &v) const noexcept
Calculates the cross product.
Definition: V3D.h:274
double normalize()
Make a normalized vector (return norm value)
Definition: V3D.cpp:130
void rotate(const Matrix< double > &) noexcept
Rotate a point by a matrix.
Definition: V3D.cpp:217
int section(std::string &A, T &out)
Convert and cut a string.
Definition: Strings.cpp:573
MANTID_KERNEL_DLL void writeMCNPX(const std::string &Line, std::ostream &OX)
Write file in standard MCNPX input form.
Definition: Strings.cpp:421
constexpr double Tolerance
Standard tolerance value.
Definition: Tolerance.h:12
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
Definition: V3D.h:341