Mantid
Loading...
Searching...
No Matches
Plane.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 +
12
13#include <limits>
14
15#ifdef ENABLE_OPENCASCADE
16// Opencascade defines _USE_MATH_DEFINES without checking whether it is already
17// used.
18// Undefine it here before we include the headers to avoid a warning
19#ifdef _MSC_VER
20#undef _USE_MATH_DEFINES
21#ifdef M_SQRT1_2
22#undef M_SQRT1_2
23#endif
24#endif
25
27GNU_DIAG_OFF("conversion")
28GNU_DIAG_OFF("cast-qual")
29#include <BRepAlgoAPI_Common.hxx>
30#include <BRepBuilderAPI_MakeFace.hxx>
31#include <BRepPrimAPI_MakeBox.hxx>
32#include <BRepPrimAPI_MakeHalfSpace.hxx>
33#include <gp_Pln.hxx>
34GNU_DIAG_ON("conversion")
35GNU_DIAG_ON("cast-qual")
36#endif
37
38namespace Mantid::Geometry {
39namespace {
40Kernel::Logger logger("Plane");
41}
43using Kernel::V3D;
44
46 : Quadratic(), m_normVec(1.0, 0.0, 0.0), m_distance(0)
50{
51 setBaseEqn();
52}
53
59{
60 return new Plane(*this);
61}
62
63std::unique_ptr<Plane> Plane::clone() const
68{
69 return std::unique_ptr<Plane>(doClone());
70}
71
72int Plane::setSurface(const std::string &Pstr)
82{
83 // Two types of plane string p[x-z] and p
84 std::string Line = Pstr;
85 std::string item;
86
87 if (!Mantid::Kernel::Strings::section(Line, item) || tolower(item[0]) != 'p')
88 return -1;
89 // Only 3 need to be declared
90 double surf[9] = {0.0, 0, 0, 0, 0};
91
92 if (item.size() == 1) // PROCESS BASIC PLANE:
93 {
94 int cnt;
95 for (cnt = 0; cnt < 9 && Mantid::Kernel::Strings::section(Line, surf[cnt]); cnt++)
96 ;
97
98 if (cnt != 4 && cnt != 9)
99 return -3;
100 if (cnt == 9) // V3D type
101 {
102 Kernel::V3D A = Kernel::V3D(surf[0], surf[1], surf[2]);
103 Kernel::V3D B = Kernel::V3D(surf[3], surf[4], surf[5]);
104 Kernel::V3D C = Kernel::V3D(surf[6], surf[7], surf[8]);
105 B -= A;
106 C -= A;
107 m_normVec = B * C;
110 } else // Norm Equation:
111 {
112 m_normVec = Kernel::V3D(surf[0], surf[1], surf[2]);
113 const double ll = m_normVec.normalize();
114 if (ll < Tolerance) // avoid
115 return -4;
116 m_distance = surf[3] / ll;
117 }
118 } else if (item.size() == 2) // PROCESS px type PLANE
119 {
120 const auto ptype = static_cast<int>(tolower(item[1]) - 'x');
121 if (ptype < 0 || ptype > 2) // Not x,y,z
122 return -5;
123 surf[ptype] = 1.0;
125 return -6; // Too short or no number
126 m_normVec = Kernel::V3D(surf[0], surf[1], surf[2]);
127 } else
128 return -3; // WRONG NAME
129
130 setBaseEqn();
131 return 0;
132}
133
141{
142 try {
143 m_normVec = normalize(N);
144 } catch (std::runtime_error &) {
145 throw std::invalid_argument("Attempt to create Plane with zero normal");
146 }
148 setBaseEqn();
149 return 0;
150}
151
157{
158 m_normVec.rotate(MA);
161}
162
169{
172}
173
174double Plane::distance(const Kernel::V3D &A) const
181{
182 return A.scalar_prod(m_normVec) - m_distance;
183}
184
185double Plane::dotProd(const Plane &A) const
190{
191 return m_normVec.scalar_prod(A.m_normVec);
192}
193
200{
201 return m_normVec.cross_prod(A.m_normVec);
202}
203
204int Plane::side(const Kernel::V3D &A) const
212{
213 const double Dp = m_normVec.scalar_prod(A) - m_distance;
214 if (Tolerance < std::abs(Dp))
215 return (Dp > 0) ? 1 : -1;
216 return 0;
217}
218
219bool Plane::onSurface(const Kernel::V3D &A) const
225{
226 return (side(A) == 0);
227}
228
229void Plane::print() const
234{
236 logger.debug() << "m_normVec == " << m_normVec << " : " << m_distance << '\n';
237}
238
239std::size_t Plane::planeType() const
247{
248 for (std::size_t i = 0; i < 3; i++)
249 if (fabs(m_normVec[i]) > (1.0 - Tolerance))
250 return i + 1;
251 return 0;
252}
253
258 BaseEqn[0] = 0.0; // A x^2
259 BaseEqn[1] = 0.0; // B y^2
260 BaseEqn[2] = 0.0; // C z^2
261 BaseEqn[3] = 0.0; // D xy
262 BaseEqn[4] = 0.0; // E xz
263 BaseEqn[5] = 0.0; // F yz
264 BaseEqn[6] = m_normVec[0]; // G x
265 BaseEqn[7] = m_normVec[1]; // H y
266 BaseEqn[8] = m_normVec[2]; // J z
267 BaseEqn[9] = -m_distance; // K const
268}
269
275void Plane::write(std::ostream &OX) const {
276 std::ostringstream cx;
278 cx.precision(Surface::Nprecision);
279 const std::size_t ptype = planeType();
280 if (!ptype)
281 cx << "p " << m_normVec[0] << " " << m_normVec[1] << " " << m_normVec[2] << " " << m_distance;
282 else if (m_normVec[ptype - 1] < 0)
283 cx << "p"
284 << "xyz"[ptype - 1] << " " << -m_distance;
285 else
286 cx << "p"
287 << "xyz"[ptype - 1] << " " << m_distance;
288
290}
291
299int Plane::LineIntersectionWithPlane(V3D startpt, V3D endpt, V3D &output) {
300 double const sprod = this->getNormal().scalar_prod(startpt - endpt);
301 if (sprod == 0)
302 return 0;
303 double const projection = m_normVec[0] * startpt[0] + m_normVec[1] * startpt[1] + m_normVec[2] * startpt[2];
304 double s1 = (projection - m_distance) / sprod;
305 if (s1 < 0 || s1 > 1)
306 return 0;
307 // The expressions below for resolving the point of intersection are
308 // resilient to the corner m_distance << sprod.
309 double const ratio = projection / sprod;
310 output[0] = ratio * endpt[0] + (1 - ratio) * startpt[0] - ((endpt[0] - startpt[0]) / sprod) * m_distance;
311 output[1] = ratio * endpt[1] + (1 - ratio) * startpt[1] - ((endpt[1] - startpt[1]) / sprod) * m_distance;
312 output[2] = ratio * endpt[2] + (1 - ratio) * startpt[2] - ((endpt[2] - startpt[2]) / sprod) * m_distance;
313 return 1;
314}
315
327void Plane::getBoundingBox(double &xmax, double &ymax, double &zmax, double &xmin, double &ymin, double &zmin) {
328 // to get the bounding box calculate the normal and the starting point
329 V3D vertex1(xmin, ymin, zmin);
330 V3D vertex2(xmax, ymin, zmin);
331 V3D vertex3(xmax, ymax, zmin);
332 V3D vertex4(xmin, ymax, zmin);
333 V3D vertex5(xmin, ymin, zmax);
334 V3D vertex6(xmax, ymin, zmax);
335 V3D vertex7(xmax, ymax, zmax);
336 V3D vertex8(xmin, ymax, zmax);
337 // find which points lie on which side of the plane
338 // find where the plane cuts the cube
339 //(xmin,ymin,zmin)--- (xmax,ymin,zmin) 1
340 //(xmax,ymin,zmin)--- (xmax,ymin,zmax) 2
341 //(xmax,ymin,zmax)--- (xmin,ymin,zmax) 3
342 //(xmin,ymin,zmax)--- (xmin,ymin,zmin) 4
343 //(xmin,ymax,zmin)--- (xmax,ymax,zmin) 5
344 //(xmax,ymax,zmin)--- (xmax,ymax,zmax) 6
345 //(xmax,ymax,zmax)--- (xmin,ymax,zmax) 7
346 //(xmin,ymax,zmax)--- (xmin,ymax,zmin) 8
347 //(xmin,ymin,zmin)--- (xmin,ymax,zmin) 9
348 //(xmax,ymin,zmin)--- (xmax,ymax,zmin) 10
349 //(xmax,ymin,zmax)--- (xmax,ymax,zmax) 11
350 //(xmin,ymin,zmax)--- (xmin,ymax,zmax) 12
351 std::vector<V3D> listOfPoints;
352 if (this->side(vertex1) <= 0)
353 listOfPoints.emplace_back(vertex1);
354 if (this->side(vertex2) <= 0)
355 listOfPoints.emplace_back(vertex2);
356 if (this->side(vertex3) <= 0)
357 listOfPoints.emplace_back(vertex3);
358 if (this->side(vertex4) <= 0)
359 listOfPoints.emplace_back(vertex4);
360 if (this->side(vertex5) <= 0)
361 listOfPoints.emplace_back(vertex5);
362 if (this->side(vertex6) <= 0)
363 listOfPoints.emplace_back(vertex6);
364 if (this->side(vertex7) <= 0)
365 listOfPoints.emplace_back(vertex7);
366 if (this->side(vertex8) <= 0)
367 listOfPoints.emplace_back(vertex8);
368 V3D edge1, edge2, edge3, edge4, edge5, edge6, edge7, edge8, edge9, edge10, edge11, edge12;
369 if (LineIntersectionWithPlane(vertex1, vertex2, edge1) == 1)
370 listOfPoints.emplace_back(edge1);
371 if (LineIntersectionWithPlane(vertex2, vertex3, edge2) == 1)
372 listOfPoints.emplace_back(edge2);
373 if (LineIntersectionWithPlane(vertex3, vertex4, edge3) == 1)
374 listOfPoints.emplace_back(edge3);
375 if (LineIntersectionWithPlane(vertex4, vertex1, edge4) == 1)
376 listOfPoints.emplace_back(edge4);
377 if (LineIntersectionWithPlane(vertex5, vertex6, edge5) == 1)
378 listOfPoints.emplace_back(edge5);
379 if (LineIntersectionWithPlane(vertex6, vertex7, edge6) == 1)
380 listOfPoints.emplace_back(edge6);
381 if (LineIntersectionWithPlane(vertex7, vertex8, edge7) == 1)
382 listOfPoints.emplace_back(edge7);
383 if (LineIntersectionWithPlane(vertex8, vertex5, edge8) == 1)
384 listOfPoints.emplace_back(edge8);
385 if (LineIntersectionWithPlane(vertex1, vertex5, edge9) == 1)
386 listOfPoints.emplace_back(edge9);
387 if (LineIntersectionWithPlane(vertex2, vertex6, edge10) == 1)
388 listOfPoints.emplace_back(edge10);
389 if (LineIntersectionWithPlane(vertex3, vertex7, edge11) == 1)
390 listOfPoints.emplace_back(edge11);
391 if (LineIntersectionWithPlane(vertex4, vertex8, edge12) == 1)
392 listOfPoints.emplace_back(edge12);
393 // now sort the vertices to find the mins and max
394 // std::cout<<listOfPoints.size()<<'\n';
395 if (!listOfPoints.empty()) {
396 xmin = ymin = zmin = std::numeric_limits<double>::max();
397 xmax = ymax = zmax = std::numeric_limits<double>::lowest();
398 for (std::vector<V3D>::const_iterator it = listOfPoints.begin(); it != listOfPoints.end(); ++it) {
399 // std::cout<<(*it)<<'\n';
400 if ((*it)[0] < xmin)
401 xmin = (*it)[0];
402 if ((*it)[1] < ymin)
403 ymin = (*it)[1];
404 if ((*it)[2] < zmin)
405 zmin = (*it)[2];
406 if ((*it)[0] > xmax)
407 xmax = (*it)[0];
408 if ((*it)[1] > ymax)
409 ymax = (*it)[1];
410 if ((*it)[2] > zmax)
411 zmax = (*it)[2];
412 }
413 }
414}
415
416#ifdef ENABLE_OPENCASCADE
417TopoDS_Shape Plane::createShape() {
418 // Get Plane normal and distance.
419 V3D normal = this->getNormal();
420 double norm2 = normal.norm2();
421 if (norm2 == 0.0) {
422 throw std::runtime_error("Cannot create a plane with zero normal");
423 }
424 double distance = this->getDistance();
425 // Find point closest to origin
426 double t = distance / norm2;
427 // Create Half Space
428 TopoDS_Face P = BRepBuilderAPI_MakeFace(gp_Pln(normal[0], normal[1], normal[2], -distance)).Face();
429
430 TopoDS_Shape Result =
431 BRepPrimAPI_MakeHalfSpace(P, gp_Pnt(normal[0] * (1 + t), normal[1] * (1 + t), normal[2] * (1 + t))).Solid();
432 return Result.Complemented();
433}
434#endif
435
436} // namespace Mantid::Geometry
#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.
Impliments a line.
Definition: Line.h:43
Holds a simple Plane.
Definition: Plane.h:35
void setBaseEqn() override
set up to be eqn based
Definition: Plane.cpp:257
void getBoundingBox(double &xmax, double &ymax, double &zmax, double &xmin, double &ymin, double &zmin) override
Returns the bounding box values for plane, double max is infinity and double min is -infinity A very ...
Definition: Plane.cpp:327
int side(const Kernel::V3D &) const override
Calcualates the side that the point is on.
Definition: Plane.cpp:204
double m_distance
Distance from origin.
Definition: Plane.h:38
std::unique_ptr< Plane > clone() const
Makes a clone (implicit virtual copy constructor)
Definition: Plane.cpp:63
bool onSurface(const Kernel::V3D &) const override
Calcuate the side that the point is on and returns success if it is on the surface.
Definition: Plane.cpp:219
double distance(const Kernel::V3D &) const override
distance from a point
Definition: Plane.cpp:174
double getDistance() const
Distance from origin.
Definition: Plane.h:65
int setSurface(const std::string &) override
processes a standard MCNPX plane string: There are three types :
Definition: Plane.cpp:72
int LineIntersectionWithPlane(Kernel::V3D startpt, Kernel::V3D endpt, Kernel::V3D &output)
Returns the point of intersection of line with the plane.
Definition: Plane.cpp:299
double dotProd(const Plane &) const
returns normal dot product
Definition: Plane.cpp:185
const Kernel::V3D & getNormal() const
Normal to plane (+ve surface)
Definition: Plane.h:66
void print() const override
Prints out the surface info and the Plane info.
Definition: Plane.cpp:229
int setPlane(const Kernel::V3D &, const Kernel::V3D &)
Given a point and a normal direction set the plane.
Definition: Plane.cpp:134
void write(std::ostream &) const override
Write in MCNPX form.
Definition: Plane.cpp:275
Plane * doClone() const override
Makes a clone (implicit virtual copy constructor)
Definition: Plane.cpp:54
Plane()
Constructor: sets plane in y-z plane and throught origin.
Definition: Plane.cpp:45
std::size_t planeType() const
are we alined on an axis
Definition: Plane.cpp:239
void rotate(const Kernel::Matrix< double > &) override
Rotate the plane about the origin by MA.
Definition: Plane.cpp:152
Kernel::V3D crossProd(const Plane &) const
returns normal cross product
Definition: Plane.cpp:194
Kernel::V3D m_normVec
Normal vector.
Definition: Plane.h:37
void displace(const Kernel::V3D &) override
Displace the plane by Point Sp.
Definition: Plane.cpp:163
Holds a basic quadratic surface.
Definition: Quadratic.h:29
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 writeHeader(std::ostream &) const
Writes out the start of an MCNPX surface description .
Definition: Surface.cpp:71
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
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
constexpr V3D cross_prod(const V3D &v) const noexcept
Cross product (this * argument)
Definition: V3D.h:278
void rotate(const Matrix< double > &) noexcept
Rotate a point by a matrix.
Definition: V3D.cpp:217
constexpr double norm2() const noexcept
Vector length squared.
Definition: V3D.h:265
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
int convert(const std::string &A, T &out)
Convert a string into a number.
Definition: Strings.cpp:665
constexpr double Tolerance
Standard tolerance value.
Definition: Tolerance.h:12
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
Definition: V3D.h:341