PeriDEM 0.3.0
PeriDEM -- Peridynamics-based high-fidelity model for granular media
Loading...
Searching...
No Matches
geomObjects.cpp
Go to the documentation of this file.
1/*
2 * -------------------------------------------
3 * Copyright (c) 2021 - 2026 Prashant K. Jha
4 * -------------------------------------------
5 * PeriDEM https://github.com/prashjha/PeriDEM
6 *
7 * Distributed under the Boost Software License, Version 1.0. (See accompanying
8 * file LICENSE)
9 */
10
11#include "geomObjects.h"
12#include "geomUtilFunctions.h"
13#include "util/function.h"
14#include "util/vecMethods.h"
15#include "util/io.h"
16#include <algorithm>
17#include <cmath>
18#include <set>
19#include <stdexcept>
20#include <vector>
21
22namespace {
23 std::string printErrMsg(const std::string &geom_type,
24 const std::vector<double> &params,
25 const std::vector<size_t> &num_params_needed) {
26
27 std::ostringstream oss;
28
29 oss << "Error: Number of parameters needed to create geometry = "
30 << geom_type << " are "
31 << util::io::printStr(num_params_needed, 0)
32 << ". But the number of parameters provided are "
33 << params.size()
34 << " and the parameters are "
35 << util::io::printStr(params, 0)
36 << ". Exiting.\n";
37
38 return oss.str();
39 }
40};
41
42//
43// Plane
44//
45namespace geom {
47 WallContactHit &hit) const {
48 hit = WallContactHit();
49 const double ln = d_n.length();
50 if (!(ln > 1.e-16))
51 return false;
52 const util::Point n = d_n / ln;
53 const double gap = (x - d_p) * n;
54 hit.active = true;
55 hit.signed_gap = gap;
56 hit.outward_n = n;
57 hit.closest = x - gap * n;
58 return true;
59 }
60} // Plane
61
62//
63// Line
64//
65namespace geom {
66
68 WallContactHit &hit) const {
69 hit = WallContactHit();
70 if (d_vertices.size() < 2)
71 return false;
72 const util::Point &a = d_vertices[0];
73 const util::Point &b = d_vertices[1];
74 const util::Point t = b - a;
75 const double t2 = t.lengthSq();
76 if (!(t2 > 1.e-24))
77 return false;
78
79 // Outward = left of directed segment in xy (CCW rotate).
80 util::Point n(-t.d_y, t.d_x, 0.);
81 const double nl = n.length();
82 if (!(nl > 1.e-16))
83 return false;
84 n = n / nl;
85
86 double s = ((x - a) * t) / t2;
87 if (s < 0.)
88 s = 0.;
89 else if (s > 1.)
90 s = 1.;
91 hit.closest = a + s * t;
92 const util::Point dx = x - hit.closest;
93 if (s > 1.e-12 && s < 1. - 1.e-12) {
94 hit.signed_gap = dx * n;
95 hit.outward_n = n;
96 } else {
97 const double dist = dx.length();
98 if (!(dist > 1.e-16)) {
99 hit.signed_gap = 0.;
100 hit.outward_n = n;
101 } else {
102 hit.outward_n = dx / dist;
103 hit.signed_gap = (hit.outward_n * n >= 0.) ? dist : -dist;
104 if (hit.outward_n * n < 0.)
105 hit.outward_n = n;
106 }
107 }
108 hit.active = true;
109 return true;
110 }
111
112 double Line::volume() const {
113 return d_L;
114 }
115
117 return d_x;
118 }
119
120 std::pair<util::Point, util::Point> Line::box() const {
121
122 return box(0.);
123 }
124
125 std::pair<util::Point, util::Point> Line::box(const
126 double &tol) const {
127
128 return {d_vertices[0] - tol, d_vertices[1] + tol};
129 }
130
131 double Line::inscribedRadius() const {
132
133 return d_r;
134 }
135
136 double Line::boundingRadius() const {
137
138 return d_r;
139 }
140
141 bool Line::isInside(const util::Point &x) const {
142
143 auto da = (d_vertices[1] - d_vertices[0]) / d_L;
144 auto db = x - d_vertices[0];
145 double dot = db * da;
146
147 if (util::isLess(dot, 0.) or util::isGreater(dot, d_L))
148 return false;
149
150
151 auto dx = db - dot * da;
152
153 return util::isLess(dx.length(), 1.0E-10);
154 }
155
156 bool Line::isOutside(const util::Point &x) const {
157
158 return !isInside(x);
159 }
160
162 const double &tol) const {
163
164 auto da = (d_vertices[1] - d_vertices[0]) / d_L;
165 auto db = x - d_vertices[0];
166 double dot = db * da;
167
168 if (util::isLess(dot, 0.) or util::isGreater(dot, d_L))
169 return false;
170
171 auto dx = db - dot * da;
172
173 return util::isLess(dx.length(), tol);
174 }
175
177 const double &tol, const bool
178 &within) const {
179
180 // check if particle is inside the object
181 if (!isNear(x, within ? 0. : tol))
182 return false;
183
184 auto da = (d_vertices[1] - d_vertices[0]) / d_L;
185 auto db = x - d_vertices[0];
186 double dot = db * da;
187
188 if (util::isLess(dot, 0.) or util::isGreater(dot, tol)
189 or util::isGreater(dot, d_L) or util::isLess(dot, d_L - tol))
190 return false;
191
192 auto dx = db - dot * da;
193
194 return util::isLess(dx.length(), tol);
195 }
196
197 bool Line::doesIntersect(const util::Point &x) const {
198
199 return isNearBoundary(x, 1.0E-8, false);
200 }
201
203 const std::pair<util::Point, util::Point> &box) const {
204
205 return false;
206 }
207
209 const std::pair<util::Point, util::Point> &box) const {
210
211 return true;
212 }
213
215 const std::pair<util::Point, util::Point> &box,
216 const double &tol) const {
217
218 return true;
219 }
220
222 const std::pair<util::Point, util::Point> &box) const {
223
224 return false;
225 }
226
227 std::string Line::printStr(int nt, int lvl) const {
228
229 auto tabS = util::io::getTabS(nt);
230
231 std::ostringstream oss;
232
233 oss << tabS << "------- Line --------" << std::endl << std::endl;
234 oss << tabS << "Name = " << d_name << std::endl;
235 oss << tabS << "Length = " << d_L << std::endl;
236 oss << tabS << "Point 1 = " << d_vertices[0].printStr(0, lvl) << std::endl;
237 oss << tabS << "Point 2 = " << d_vertices[1].printStr(0, lvl) << std::endl;
238 oss << std::endl;
239
240 if (lvl > 0)
241 oss << tabS << "Bounding box: "
242 << util::io::printBoxStr(box(0.), nt + 1);
243
244 if (lvl == 0)
245 oss << std::endl;
246
247 return oss.str();
248 }
249} //Line
250
251//
252// Triangle
253//
254namespace geom {
255 double Triangle::volume() const {
256 // 2D area of the stored triangle (works for equilateral and vertex-defined).
257 return std::abs(geom::triangleArea(d_vertices[0], d_vertices[1], d_vertices[2]));
258 }
259
261 return d_x;
262 }
263
264 std::pair<util::Point, util::Point> Triangle::box() const {
265
266 return box(0.);
267 }
268
269 std::pair<util::Point, util::Point> Triangle::box(const
270 double &tol) const {
271
272 util::Point p1 = d_vertices[0], p2 = d_vertices[0];
273 for (const auto &v : d_vertices) {
274 p1.d_x = std::min(p1.d_x, v.d_x);
275 p1.d_y = std::min(p1.d_y, v.d_y);
276 p1.d_z = std::min(p1.d_z, v.d_z);
277 p2.d_x = std::max(p2.d_x, v.d_x);
278 p2.d_y = std::max(p2.d_y, v.d_y);
279 p2.d_z = std::max(p2.d_z, v.d_z);
280 }
281 return {p1 - tol, p2 + tol};
282 }
283
285
286 return d_r * std::sin(M_PI / 3);
287 }
288
290
291 return d_r;
292 }
293
294 bool Triangle::isInside(const util::Point &x) const {
295
296 if ((x - d_x).length() > d_r)
297 return false;
298
299 if ((x - d_x).length() < this->inscribedRadius())
300 return true;
301
302 double a = this->volume();
303 double a1 =
304 std::abs(geom::triangleArea(x, d_vertices[1], d_vertices[2]));
305 double a2 =
306 std::abs(geom::triangleArea(d_vertices[0], x, d_vertices[2]));
307 double a3 =
308 std::abs(geom::triangleArea(d_vertices[0], d_vertices[1], x));
309
310 return (a1 + a2 + a3 <= a * (1.0 + 1.0e-10));
311 }
312
313 bool Triangle::isOutside(const util::Point &x) const {
314 return !isInside(x);
315 }
316
318 const double &tol) const {
319
320 // get a bigger box containing this object
321 auto bbox = box(tol);
322
323 return geom::isPointInsideBox(x, 2, bbox);
324 }
325
327 const double &tol, const bool
328 &within) const {
329
330 // check if particle is inside the object
331 if (!isNear(x, within ? 0. : tol))
332 return false;
333
334 double a = this->volume();
335 double l = 0.5 * std::sqrt(a);
336
337 double a1 =
338 std::abs(geom::triangleArea(x, d_vertices[1], d_vertices[2]));
339 if (a1 < tol * l)
340 return true;
341
342 double a2 =
343 std::abs(geom::triangleArea(d_vertices[0], x, d_vertices[2]));
344 if (a2 < tol * l)
345 return true;
346
347 double a3 =
348 std::abs(geom::triangleArea(d_vertices[0], d_vertices[1], x));
349 if (a3 < tol * l)
350 return true;
351
352 return false;
353 }
354
356
357 return isNearBoundary(x, 1.0E-8, false);
358 }
359
361 const std::pair<util::Point, util::Point> &box) const {
362
363 for (auto p: geom::getCornerPoints(2, box))
364 if (!this->isInside(p))
365 return false;
366
367 return true;
368 }
369
371 const std::pair<util::Point, util::Point> &box) const {
372
373 bool intersect = false;
374 for (auto p: geom::getCornerPoints(2, box))
375 if (!intersect)
376 intersect = this->isInside(p);
377
378 return !intersect;
379 }
380
382 const std::pair<util::Point, util::Point> &box,
383 const double &tol) const {
384
385 return geom::areBoxesNear(this->box(), box, tol, 2);
386 }
387
389 const std::pair<util::Point, util::Point> &box) const {
390
391 // need to check all four corner points
392 for (auto p: geom::getCornerPoints(2, box))
393 if (this->isInside(p))
394 return true;
395
396 return false;
397 }
398
399 std::string Triangle::printStr(int nt, int lvl) const {
400
401 auto tabS = util::io::getTabS(nt);
402
403 std::ostringstream oss;
404
405 oss << tabS << "------- Triangle --------" << std::endl << std::endl;
406 oss << tabS << "Name = " << d_name << std::endl;
407 oss << tabS << "Center = " << d_x.printStr(0, lvl) << std::endl;
408 oss << tabS << "Radius = " << d_r << std::endl;
409 oss << tabS << "Vertices = " << util::io::printStr(d_vertices)
410 << std::endl;
411 oss << std::endl;
412
413 if (lvl == 0)
414 oss << std::endl;
415
416 return oss.str();
417 }
418}// Triangle
419
420//
421// Square
422//
423namespace geom {
424 double Square::volume() const {
425 return std::pow(d_L, 2);
426 }
427
429 return d_x;
430 }
431
432 std::pair<util::Point, util::Point> Square::box() const {
433
434 return box(0.);
435 }
436 std::pair<util::Point, util::Point> Square::box(const
437 double &tol) const {
438
439 return {util::Point(d_vertices[0].d_x - tol, d_vertices[0].d_y - tol,
440 0.),
441 util::Point(d_vertices[2].d_x + tol, d_vertices[2].d_y + tol,
442 0.)};
443 }
444
445 double Square::inscribedRadius() const {
446
447 return 0.5*d_L;
448 }
449
450 double Square::boundingRadius() const {
451
452 return d_r;
453 }
454
455 bool Square::isInside(const util::Point &x) const {
457 }
458
459 bool Square::isOutside(const util::Point &x) const {
460 return !isInside(x);
461 }
462
464 const double &tol) const {
465
466 // get a bigger box containing this object
467 auto bbox = box(tol);
468
469 return geom::isPointInsideBox(x, 2, bbox);
470 }
471
473 const double &tol, const bool
474 &within) const {
475
476 // check if particle is inside the object
477 if (!isNear(x, within ? 0. : tol))
478 return false;
479
480 bool near_x_edge = util::isLess(std::abs(x.d_x - d_vertices[0].d_x), tol) or
481 util::isLess(std::abs(x.d_x - d_vertices[2].d_x), tol);
482
483 bool near_y_edge = util::isLess(std::abs(x.d_y - d_vertices[0].d_y), tol) or
484 util::isLess(std::abs(x.d_y - d_vertices[2].d_y), tol);
485
486 return near_x_edge || near_y_edge;
487 }
488
489 bool Square::doesIntersect(const util::Point &x) const {
490
491 return isNearBoundary(x, 1.0E-8, false);
492 }
493
495 const std::pair<util::Point, util::Point> &box) const {
496
497 for (auto p : geom::getCornerPoints(2, box))
498 if(!this->isInside(p))
499 return false;
500
501 return true;
502 }
503
505 const std::pair<util::Point, util::Point> &box) const {
506
507 bool intersect = false;
508 for (auto p : geom::getCornerPoints(2, box))
509 if (!intersect)
510 intersect = this->isInside(p);
511
512 return !intersect;
513 }
514
516 const std::pair<util::Point, util::Point> &box, const double &tol) const {
517
518 return geom::areBoxesNear(this->box(), box, tol, 2);
519 }
520
522 const std::pair<util::Point, util::Point> &box) const {
523
524 // need to check all four corner points
525 for (auto p : geom::getCornerPoints(2, box))
526 if (this->isInside(p))
527 return true;
528
529 return false;
530 }
531
532 std::string Square::printStr(int nt, int lvl) const {
533
534 auto tabS = util::io::getTabS(nt);
535
536 std::ostringstream oss;
537
538 oss << tabS << "------- Rectangle --------" << std::endl << std::endl;
539 oss << tabS << "Name = " << d_name << std::endl;
540 oss << tabS << "Length = " << d_L << std::endl;
541 oss << tabS << "Bounding radius = " << d_r << std::endl;
542 oss << tabS << "Center = " << d_x.printStr(0, lvl) << std::endl;
543 oss << tabS << "Vertices = " << util::io::printStr(d_vertices, 0) << std::endl;
544 oss << std::endl;
545
546 if (lvl > 0)
547 oss << tabS << "Bounding box: " << util::io::printBoxStr(box(0.), nt + 1);
548
549 if (lvl == 0)
550 oss << std::endl;
551
552 return oss.str();
553 }
554}// Square
555
556//
557// Rectangle
558//
559namespace geom {
561 WallContactHit &hit) const {
562 hit = WallContactHit();
563 // Axis-aligned SDF from current bounding box (walls are typically AA).
564 const auto bb = box();
565 const double x0 = bb.first.d_x;
566 const double y0 = bb.first.d_y;
567 const double x1 = bb.second.d_x;
568 const double y1 = bb.second.d_y;
569 if (!(x1 > x0) || !(y1 > y0))
570 return false;
571
572 const double dx = std::max(x0 - x.d_x, x.d_x - x1);
573 const double dy = std::max(y0 - x.d_y, x.d_y - y1);
574 const bool outside = dx > 0. || dy > 0.;
575
576 if (outside) {
577 // Closest point on AABB, gap = exterior distance.
578 util::Point c = x;
579 if (x.d_x < x0)
580 c.d_x = x0;
581 else if (x.d_x > x1)
582 c.d_x = x1;
583 if (x.d_y < y0)
584 c.d_y = y0;
585 else if (x.d_y > y1)
586 c.d_y = y1;
587 c.d_z = x.d_z;
588 const util::Point d = x - c;
589 const double dist = d.length();
590 if (!(dist > 1.e-16)) {
591 // On a face: pick dominant exterior axis.
592 if (dx >= dy) {
593 hit.outward_n = util::Point(x.d_x < x0 ? -1. : 1., 0., 0.);
594 } else {
595 hit.outward_n = util::Point(0., x.d_y < y0 ? -1. : 1., 0.);
596 }
597 hit.signed_gap = 0.;
598 hit.closest = c;
599 } else {
600 hit.outward_n = d / dist;
601 hit.signed_gap = dist;
602 hit.closest = c;
603 }
604 } else {
605 // Inside: negative distance to nearest face.
606 const double dl = x.d_x - x0;
607 const double dr = x1 - x.d_x;
608 const double db = x.d_y - y0;
609 const double dt = y1 - x.d_y;
610 const double m = std::min(std::min(dl, dr), std::min(db, dt));
611 hit.signed_gap = -m;
612 if (m == dl) {
613 hit.outward_n = util::Point(-1., 0., 0.);
614 hit.closest = util::Point(x0, x.d_y, x.d_z);
615 } else if (m == dr) {
616 hit.outward_n = util::Point(1., 0., 0.);
617 hit.closest = util::Point(x1, x.d_y, x.d_z);
618 } else if (m == db) {
619 hit.outward_n = util::Point(0., -1., 0.);
620 hit.closest = util::Point(x.d_x, y0, x.d_z);
621 } else {
622 hit.outward_n = util::Point(0., 1., 0.);
623 hit.closest = util::Point(x.d_x, y1, x.d_z);
624 }
625 }
626 hit.active = true;
627 return true;
628 }
629
630 double Rectangle::volume() const {
631 return d_Lx * d_Ly;
632 }
633
635 return d_x;
636 }
637
638 std::pair<util::Point, util::Point> Rectangle::box() const {
639
640 return box(0.);
641 }
642
643 std::pair<util::Point, util::Point> Rectangle::box(const
644 double &tol) const {
645
646 return {util::Point(d_vertices[0].d_x - tol, d_vertices[0].d_y - tol,
647 0.),
648 util::Point(d_vertices[2].d_x + tol, d_vertices[2].d_y + tol,
649 0.)};
650 }
651
653
654 return util::isLess(d_Lx, d_Ly) ? d_Lx : d_Ly;
655 }
656
658
659 return d_r;
660 }
661
662 bool Rectangle::isInside(const util::Point &x) const {
664 }
665
666 bool Rectangle::isOutside(const util::Point &x) const {
667 return !isInside(x);
668 }
669
671 const double &tol) const {
672
673 // get a bigger box containing this object
674 auto bbox = box(tol);
675
676 return geom::isPointInsideBox(x, 2, bbox);
677 }
678
680 const double &tol, const bool
681 &within) const {
682
683 // check if particle is inside the object
684 if (!isNear(x, within ? 0. : tol))
685 return false;
686
687 bool near_x_edge = util::isLess(std::abs(x.d_x - d_vertices[0].d_x), tol) or
688 util::isLess(std::abs(x.d_x - d_vertices[2].d_x), tol);
689
690 bool near_y_edge = util::isLess(std::abs(x.d_y - d_vertices[0].d_y), tol) or
691 util::isLess(std::abs(x.d_y - d_vertices[2].d_y), tol);
692
693 return near_x_edge || near_y_edge;
694 }
695
697
698 return isNearBoundary(x, 1.0E-8, false);
699 }
700
702 const std::pair<util::Point, util::Point> &box) const {
703
704 for (auto p: geom::getCornerPoints(2, box))
705 if (!this->isInside(p))
706 return false;
707
708 return true;
709 }
710
712 const std::pair<util::Point, util::Point> &box) const {
713
714 bool intersect = false;
715 for (auto p: geom::getCornerPoints(2, box))
716 if (!intersect)
717 intersect = this->isInside(p);
718
719 return !intersect;
720 }
721
723 const std::pair<util::Point, util::Point> &box,
724 const double &tol) const {
725
726 return geom::areBoxesNear(this->box(), box, tol, 2);
727 }
728
730 const std::pair<util::Point, util::Point> &box) const {
731
732 // need to check all four corner points
733 for (auto p: geom::getCornerPoints(2, box))
734 if (this->isInside(p))
735 return true;
736
737 return false;
738 }
739
740 std::string Rectangle::printStr(int nt, int lvl) const {
741
742 auto tabS = util::io::getTabS(nt);
743
744 std::ostringstream oss;
745
746 oss << tabS << "------- Rectangle --------" << std::endl << std::endl;
747 oss << tabS << "Name = " << d_name << std::endl;
748 oss << tabS << "Lengths (Lx, Ly) = (" << d_Lx << ", " << d_Ly << ")" << std::endl;
749 oss << tabS << "Bounding circle radius = " << d_r << std::endl;
750 oss << tabS << "Vertices = " << util::io::printStr(d_vertices, 0) << std::endl;
751 oss << std::endl;
752
753 if (lvl > 0)
754 oss << tabS << "Bounding box: "
755 << util::io::printBoxStr(box(0.), nt + 1);
756
757 if (lvl == 0)
758 oss << std::endl;
759
760 return oss.str();
761 }
762}// Rectangle
763
764//
765// Hexagon
766//
767namespace geom {
768 double Hexagon::volume() const {
769 // https://en.wikipedia.org/wiki/Hexagon
770 double r_small = this->inscribedRadius();
771 return 2. * std::sqrt(3.) * r_small * r_small;
772 }
773
775 return d_x;
776 }
777
778 std::pair<util::Point, util::Point> Hexagon::box() const {
779
780 return box(0.);
781 }
782
783 std::pair<util::Point, util::Point> Hexagon::box(const
784 double &tol) const {
785
786 auto p1 = d_x - util::Point(d_r + tol, d_r + tol, d_x[2] + tol);
787 auto p2 = d_x + util::Point(d_r + tol, d_r + tol, d_x[2] + tol);
788 return {p1, p2};
789 }
790
792
793 return d_r * 0.5 * std::sqrt(3.);
794 }
795
796 double Hexagon::boundingRadius() const {
797
798 return d_r;
799 }
800
801 bool Hexagon::isInside(const util::Point &x) const {
802
803 if ((x - d_x).length() > d_r)
804 return false;
805
806 if ((x - d_x).length() < inscribedRadius())
807 return true;
808
809 return false;
810 }
811
812 bool Hexagon::isOutside(const util::Point &x) const {
813 return !isInside(x);
814 }
815
817 const double &tol) const {
818
819 // get a bigger box containing this object
820 auto bbox = box(tol);
821
822 return geom::isPointInsideBox(x, 2, bbox);
823 }
824
826 const double &tol, const bool
827 &within) const {
828
829 if ((x - d_x).length() > d_r + tol)
830 return false;
831
832 if ((x - d_x).length() < inscribedRadius() - tol)
833 return false;
834
835 return true;
836 }
837
838 bool Hexagon::doesIntersect(const util::Point &x) const {
839
840 return isNearBoundary(x, 1.0E-8, false);
841 }
842
844 const std::pair<util::Point, util::Point> &box) const {
845
846 for (auto p: geom::getCornerPoints(2, box))
847 if (!this->isInside(p))
848 return false;
849
850 return true;
851 }
852
854 const std::pair<util::Point, util::Point> &box) const {
855
856 bool intersect = false;
857 for (auto p: geom::getCornerPoints(2, box))
858 if (!intersect)
859 intersect = this->isInside(p);
860
861 return !intersect;
862 }
863
865 const std::pair<util::Point, util::Point> &box,
866 const double &tol) const {
867
868 return geom::areBoxesNear(this->box(), box, tol, 2);
869 }
870
872 const std::pair<util::Point, util::Point> &box) const {
873
874 // need to check all four corner points
875 for (auto p: geom::getCornerPoints(2, box))
876 if (this->isInside(p))
877 return true;
878
879 return false;
880 }
881
882 std::string Hexagon::printStr(int nt, int lvl) const {
883
884 auto tabS = util::io::getTabS(nt);
885
886 std::ostringstream oss;
887
888 oss << tabS << "------- Hexagon --------" << std::endl << std::endl;
889 oss << tabS << "Name = " << d_name << std::endl;
890 oss << tabS << "Radius = " << d_r << std::endl;
891 oss << tabS << "Center = " << d_x.printStr(0, lvl) << std::endl;
892 oss << tabS << "Axis = " << d_a.printStr(0, lvl) << std::endl;
893 oss << tabS << "Vertices = " << util::io::printStr(d_vertices, lvl) <<
894 std::endl;
895 oss << std::endl;
896
897 if (lvl == 0)
898 oss << std::endl;
899
900 return oss.str();
901 }
902}// Hexagon
903
904//
905// Drum2D
906//
907namespace geom {
908 double Drum2D::volume() const {
909
910 return (2. * d_r * d_r - d_r * (d_r - 2. * d_w)) * std::sin(M_PI / 3.);
911 }
912
914 return d_x;
915 }
916
917 std::pair<util::Point, util::Point> Drum2D::box() const {
918
919 return box(0.);
920 }
921
922 std::pair<util::Point, util::Point> Drum2D::box(const
923 double &tol) const {
924
925 auto p1 = d_x - util::Point(d_r + tol, d_r + tol, d_x[2] + tol);
926 auto p2 = d_x + util::Point(d_r + tol, d_r + tol, d_x[2] + tol);
927 return {p1, p2};
928 }
929
930 double Drum2D::inscribedRadius() const {
931
932 return d_w;
933 }
934
935 double Drum2D::boundingRadius() const {
936
937 return d_r;
938 }
939
940 bool Drum2D::isInside(const util::Point &x) const {
941
942 if ((x - d_x).length() > d_r)
943 return false;
944
945 if ((x - d_x).length() < inscribedRadius())
946 return true;
947
948 // rotate axis to get orthogonal axis
949 auto ortho_axis = util::rotate(d_a, M_PI * 0.5, util::Point(0., 0., 1.));
950
951 //
952 // + v2
953 // /
954 // / x
955 // /
956 // /
957 // o----+v1
958 //
959 auto ox = x - d_x;
960 double angle_ox_ov1 = std::acos(std::abs(d_a.dot(ox)) / ox.length());
961 double max_length = d_w + angle_ox_ov1 * (d_r - d_w) / (M_PI / 3.);
962
963 return ox.length() <= max_length;
964 }
965
966 bool Drum2D::isOutside(const util::Point &x) const {
967 return !isInside(x);
968 }
969
971 const double &tol) const {
972
973 // get a bigger box containing this object
974 auto bbox = box(tol);
975
976 return geom::isPointInsideBox(x, 2, bbox);
977 }
978
980 const double &tol, const bool
981 &within) const {
982
983 if ((x - d_x).length() > d_r + tol)
984 return false;
985
986 if ((x - d_x).length() < this->inscribedRadius() - tol)
987 return false;
988
989 return true;
990 }
991
992 bool Drum2D::doesIntersect(const util::Point &x) const {
993
994 return isNearBoundary(x, 1.0E-8, false);
995 }
996
998 const std::pair<util::Point, util::Point> &box) const {
999
1000 for (auto p: geom::getCornerPoints(2, box))
1001 if (!this->isInside(p))
1002 return false;
1003
1004 return true;
1005 }
1006
1008 const std::pair<util::Point, util::Point> &box) const {
1009
1010 bool intersect = false;
1011 for (auto p: geom::getCornerPoints(2, box))
1012 if (!intersect)
1013 intersect = this->isInside(p);
1014
1015 return !intersect;
1016 }
1017
1019 const std::pair<util::Point, util::Point> &box,
1020 const double &tol) const {
1021
1022 return geom::areBoxesNear(this->box(), box, tol, 2);
1023 }
1024
1026 const std::pair<util::Point, util::Point> &box) const {
1027
1028 // need to check all four corner points
1029 for (auto p: geom::getCornerPoints(2, box))
1030 if (this->isInside(p))
1031 return true;
1032
1033 return false;
1034 }
1035
1036 std::string Drum2D::printStr(int nt, int lvl) const {
1037
1038 auto tabS = util::io::getTabS(nt);
1039
1040 std::ostringstream oss;
1041
1042 oss << tabS << "------- Drum2D --------" << std::endl << std::endl;
1043 oss << tabS << "Name = " << d_name << std::endl;
1044 oss << tabS << "Radius = " << d_r << std::endl;
1045 oss << tabS << "Neck half-width = " << d_w << std::endl;
1046 oss << tabS << "Center = " << d_x.printStr(0, lvl) << std::endl;
1047 oss << tabS << "Axis = " << d_a.printStr(0, lvl) << std::endl;
1048 oss << tabS << "Vertices = " << util::io::printStr(d_vertices, lvl) <<
1049 std::endl;
1050 oss << std::endl;
1051
1052 if (lvl == 0)
1053 oss << std::endl;
1054
1055 return oss.str();
1056 }
1057}// Drum2D
1058
1059//
1060// Cube
1061//
1062namespace geom {
1063 double Cube::volume() const {
1064 return std::pow(d_L, 3);
1065 }
1066
1068 return d_x;
1069 }
1070
1071 std::pair<util::Point, util::Point> Cube::box() const {
1072
1073 return box(0.);
1074 }
1075
1076 std::pair<util::Point, util::Point> Cube::box(const
1077 double &tol) const {
1078 return {util::Point(d_vertices[0].d_x - tol, d_vertices[0].d_y - tol,
1079 d_vertices[0].d_z - tol),
1080 util::Point(d_vertices[6].d_x + tol, d_vertices[6].d_y + tol,
1081 d_vertices[6].d_z + tol)};
1082 }
1083
1084 double Cube::inscribedRadius() const {
1085
1086 return d_L;
1087 }
1088
1089 double Cube::boundingRadius() const {
1090
1091 return d_r;
1092 }
1093
1094 bool Cube::isInside(const util::Point &x) const {
1096 }
1097
1098 bool Cube::isOutside(const util::Point &x) const {
1099 return !isInside(x);
1100 }
1101
1103 const double &tol) const {
1104
1105 // get a bigger box containing this object
1106 auto bbox = box(tol);
1107
1108 return geom::isPointInsideBox(x, 3, bbox);
1109 }
1110
1112 const double &tol, const bool
1113 &within) const {
1114
1115 // check if particle is within the tolerance distance
1116 if (!isNear(x, within ? 0. : tol))
1117 return false;
1118
1119 bool near_x_edge = util::isLess(std::abs(x.d_x - d_vertices[0].d_x), tol) or
1120 util::isLess(std::abs(x.d_x - d_vertices[6].d_x), tol);
1121
1122 bool near_y_edge = util::isLess(std::abs(x.d_y - d_vertices[0].d_y), tol) or
1123 util::isLess(std::abs(x.d_y - d_vertices[6].d_y), tol);
1124
1125 bool near_z_edge = util::isLess(std::abs(x.d_z - d_vertices[0].d_z), tol) or
1126 util::isLess(std::abs(x.d_z - d_vertices[6].d_z), tol);
1127
1128 return near_x_edge || near_y_edge || near_z_edge;
1129 }
1130
1131 bool Cube::doesIntersect(const util::Point &x) const {
1132
1133 return isNearBoundary(x, 1.0E-8, false);
1134 }
1135
1137 const std::pair<util::Point, util::Point> &box) const {
1138
1139 for (auto p: geom::getCornerPoints(3, box))
1140 if (!this->isInside(p))
1141 return false;
1142
1143 return true;
1144 }
1145
1147 const std::pair<util::Point, util::Point> &box) const {
1148
1149 bool intersect = false;
1150 for (auto p: geom::getCornerPoints(3, box))
1151 if (!intersect)
1152 intersect = this->isInside(p);
1153
1154 return !intersect;
1155 }
1156
1158 const std::pair<util::Point, util::Point> &bbox, const double &tol)
1159 const {
1160
1161 return geom::areBoxesNear(box(), bbox, tol, 3);
1162 }
1163
1165 const std::pair<util::Point, util::Point> &box) const {
1166
1167 // need to check all four corner points
1168 for (auto p: geom::getCornerPoints(3, box))
1169 if (this->isInside(p))
1170 return true;
1171
1172 return false;
1173 }
1174
1175 std::string Cube::printStr(int nt, int lvl) const {
1176
1177 auto tabS = util::io::getTabS(nt);
1178
1179 std::ostringstream oss;
1180
1181 oss << tabS << "------- Cube --------" << std::endl << std::endl;
1182 oss << tabS << "Name = " << d_name << std::endl;
1183 oss << tabS << "Length = " << d_L << std::endl;
1184 oss << tabS << "Bounding sphere radius = " << d_r << std::endl;
1185 oss << tabS << "Center = " << d_x.printStr(0, 0) << std::endl;
1186 oss << tabS << "Vertices = " << util::io::printStr(d_vertices, 0) << std::endl;
1187 oss << std::endl;
1188
1189 if (lvl > 0)
1190 oss << tabS << "Bounding box: "
1191 << util::io::printBoxStr(box(0.), nt + 1);
1192
1193 if (lvl == 0)
1194 oss << std::endl;
1195
1196 return oss.str();
1197 }
1198}// Cube
1199
1200//
1201// Cuboid
1202//
1203namespace geom {
1204 double Cuboid::volume() const {
1205 return d_Lx * d_Ly * d_Lz;
1206 }
1207
1209 return d_x;
1210 }
1211
1212 std::pair<util::Point, util::Point> Cuboid::box() const {
1213
1214 return box(0.);
1215 }
1216
1217 std::pair<util::Point, util::Point> Cuboid::box(const
1218 double &tol) const {
1219 return {util::Point(d_vertices[0].d_x - tol, d_vertices[0].d_y - tol,
1220 d_vertices[0].d_z - tol),
1221 util::Point(d_vertices[6].d_x + tol, d_vertices[6].d_y + tol,
1222 d_vertices[6].d_z + tol)};
1223 }
1224
1226
1227 auto l = util::isLess(d_Lx, d_Ly) ? d_Lx : d_Ly;
1228 return util::isLess(l, d_Lz) ? l : d_Lz;
1229 }
1230
1231 double Cuboid::boundingRadius() const {
1232
1233 return d_r;
1234 }
1235
1236 bool Cuboid::isInside(const util::Point &x) const {
1238 }
1239
1240 bool Cuboid::isOutside(const util::Point &x) const {
1241 return !isInside(x);
1242 }
1243
1245 const double &tol) const {
1246
1247 // get a bigger box containing this object
1248 auto bbox = box(tol);
1249
1250 return geom::isPointInsideBox(x, 3, bbox);
1251 }
1252
1254 const double &tol, const bool
1255 &within) const {
1256
1257 // check if particle is within the tolerance distance
1258 if (!isNear(x, within ? 0. : tol))
1259 return false;
1260
1261 bool near_x_edge = util::isLess(std::abs(x.d_x - d_vertices[0].d_x), tol) or
1262 util::isLess(std::abs(x.d_x - d_vertices[6].d_x), tol);
1263
1264 bool near_y_edge = util::isLess(std::abs(x.d_y - d_vertices[0].d_y), tol) or
1265 util::isLess(std::abs(x.d_y - d_vertices[6].d_y), tol);
1266
1267 bool near_z_edge = util::isLess(std::abs(x.d_z - d_vertices[0].d_z), tol) or
1268 util::isLess(std::abs(x.d_z - d_vertices[6].d_z), tol);
1269
1270 return near_x_edge || near_y_edge || near_z_edge;
1271 }
1272
1273 bool Cuboid::doesIntersect(const util::Point &x) const {
1274
1275 return isNearBoundary(x, 1.0E-8, false);
1276 }
1277
1279 const std::pair<util::Point, util::Point> &box) const {
1280
1281 for (auto p: geom::getCornerPoints(3, box))
1282 if (!this->isInside(p))
1283 return false;
1284
1285 return true;
1286 }
1287
1289 const std::pair<util::Point, util::Point> &box) const {
1290
1291 bool intersect = false;
1292 for (auto p: geom::getCornerPoints(3, box))
1293 if (!intersect)
1294 intersect = this->isInside(p);
1295
1296 return !intersect;
1297 }
1298
1300 const std::pair<util::Point, util::Point> &bbox, const double &tol)
1301 const {
1302
1303 return geom::areBoxesNear(box(), bbox, tol, 3);
1304 }
1305
1307 const std::pair<util::Point, util::Point> &box) const {
1308
1309 // need to check all four corner points
1310 for (auto p: geom::getCornerPoints(3, box))
1311 if (this->isInside(p))
1312 return true;
1313
1314 return false;
1315 }
1316
1317 std::string Cuboid::printStr(int nt, int lvl) const {
1318
1319 auto tabS = util::io::getTabS(nt);
1320
1321 std::ostringstream oss;
1322
1323 oss << tabS << "------- Cuboid --------" << std::endl << std::endl;
1324 oss << tabS << "Name = " << d_name << std::endl;
1325 oss << tabS << "Lengths (Lx, Ly, Lz) = "
1326 << util::io::printStr(std::vector<double>{d_Lx, d_Ly, d_Lz}, 0)
1327 << std::endl;
1328 oss << tabS << "Bounding sphere radius = " << d_r << std::endl;
1329 oss << tabS << "Center = " << d_x.printStr(0, 0) << std::endl;
1330 oss << tabS << "Vertices = " << util::io::printStr(d_vertices, 0) << std::endl;
1331 oss << std::endl;
1332
1333 if (lvl > 0)
1334 oss << tabS << "Bounding box: "
1335 << util::io::printBoxStr(box(0.), nt + 1);
1336
1337 if (lvl == 0)
1338 oss << std::endl;
1339
1340 return oss.str();
1341 }
1342}// Cuboid
1343
1344//
1345// Circle
1346//
1347namespace geom {
1348 double Circle::volume() const {
1349 return M_PI * d_r * d_r;
1350 }
1351
1353 return d_x;
1354 }
1355
1356 std::pair<util::Point, util::Point> Circle::box() const {
1357
1358 return box(0.);
1359 }
1360
1361 std::pair<util::Point, util::Point> Circle::box(const
1362 double &tol) const {
1363 double r = d_r + tol;
1364 return {
1365 util::Point(d_x.d_x - r, d_x.d_y - r, 0.),
1366 util::Point(d_x.d_x + r, d_x.d_y + r, 0.)
1367 };
1368 }
1369
1371
1372 return d_r;
1373 }
1374
1375 double Circle::boundingRadius() const {
1376
1377 return d_r;
1378 }
1379
1380 bool Circle::isInside(const util::Point &x) const {
1381
1382 return util::isLess(d_x.dist(x), d_r + 1.0E-12);
1383 }
1384
1385 bool Circle::isOutside(const util::Point &x) const {
1386 return !isInside(x);
1387 }
1388
1390 const double &tol) const {
1391
1392 // translate to origin
1393 auto x0 = x - d_x;
1394
1395 return util::isLess(x0.length(), d_r + tol);
1396 }
1397
1399 const double &tol, const bool
1400 &within) const {
1401
1402 // check if particle is within the tolerance distance
1403 if (!isNear(x, within ? 0. : tol))
1404 return false;
1405
1406 // check if it is close enough to circumference
1407 auto x0 = x - d_x;
1408
1409 return util::isLess(x0.length(), d_r + tol) ||
1410 util::isLess(x0.length(), d_r - tol);
1411 }
1412
1413 bool Circle::doesIntersect(const util::Point &x) const {
1414
1415 return isNearBoundary(x, 1.0E-8, false);
1416 }
1417
1419 const std::pair<util::Point, util::Point> &box) const {
1420
1421 for (auto p: geom::getCornerPoints(2, box))
1422 if (!this->isInside(p))
1423 return false;
1424
1425 return true;
1426 }
1427
1429 const std::pair<util::Point, util::Point> &box) const {
1430
1431 bool intersect = false;
1432 for (auto p: geom::getCornerPoints(2, box))
1433 if (!intersect)
1434 intersect = this->isInside(p);
1435
1436 return !intersect;
1437 }
1438
1440 const std::pair<util::Point, util::Point> &box,
1441 const double &tol) const {
1442
1443 if (this->isInside(box))
1444 return true;
1445
1446 // get corner points of box
1447 auto cp = geom::getCornerPoints(2, box);
1448
1449 for (auto p: cp) {
1450
1451 // check the distance of corner point with the center
1452 auto dx = p - d_x;
1453 if (util::isLess(dx.length(), d_r + tol))
1454 return true;
1455 }
1456
1457 // check center to center distance
1458 auto dxc = geom::getCenter(2, box) - d_x;
1459
1460 // check wrt inscribed circle
1461 auto r = geom::inscribedRadiusInBox(2, box);
1462 if (util::isLess(dxc.length(), d_r + r + tol))
1463 return true;
1464
1465 // check wrt circumscribed circle
1467 return util::isLess(dxc.length(), d_r + r + tol);
1468 }
1469
1471 const std::pair<util::Point, util::Point> &box) const {
1472
1473 // need to check all four corner points
1474 for (auto p: geom::getCornerPoints(2, box))
1475 if (this->isInside(p))
1476 return true;
1477
1478 return false;
1479 }
1480
1481 std::string Circle::printStr(int nt, int lvl) const {
1482
1483 auto tabS = util::io::getTabS(nt);
1484
1485 std::ostringstream oss;
1486
1487 oss << tabS << "------- Circle --------" << std::endl << std::endl;
1488 oss << tabS << "Name = " << d_name << std::endl;
1489 oss << tabS << "Center = " << d_x.printStr(0, lvl) << std::endl;
1490 oss << tabS << "Radius = " << d_r << std::endl;
1491
1492 if (lvl > 0)
1493 oss << tabS << "Bounding box: "
1494 << util::io::printBoxStr(box(0.), nt + 1);
1495
1496 if (lvl == 0)
1497 oss << std::endl;
1498
1499 return oss.str();
1500 }
1501}// Circle
1502
1503//
1504// Ellipse
1505//
1506namespace geom {
1507
1508 namespace {
1509 bool ellipseMetricInside(double u, double v, double a, double b) {
1510 if (a <= 0. || b <= 0.)
1511 return false;
1512 return util::isLess(u * u / (a * a) + v * v / (b * b), 1. + 1.0E-12);
1513 }
1514
1515 void ellipseLocal(const Ellipse &e, const util::Point &x, double &u, double &v) {
1516 const double dx = x.d_x - e.d_x.d_x;
1517 const double dy = x.d_y - e.d_x.d_y;
1518 const double c = std::cos(e.d_theta);
1519 const double s = std::sin(e.d_theta);
1520 u = c * dx + s * dy;
1521 v = -s * dx + c * dy;
1522 }
1523 } // namespace
1524
1525 double Ellipse::volume() const {
1526 return M_PI * d_a * d_b;
1527 }
1528
1530 return d_x;
1531 }
1532
1533 std::pair<util::Point, util::Point> Ellipse::box() const {
1534 return box(0.);
1535 }
1536
1537 std::pair<util::Point, util::Point> Ellipse::box(const double &tol) const {
1538
1539 const double at = d_a + tol;
1540 const double bt = d_b + tol;
1541 const double ex = std::hypot(at * std::cos(d_theta), bt * std::sin(d_theta));
1542 const double ey = std::hypot(at * std::sin(d_theta), bt * std::cos(d_theta));
1543
1544 return {util::Point(d_x.d_x - ex, d_x.d_y - ey, 0.),
1545 util::Point(d_x.d_x + ex, d_x.d_y + ey, 0.)};
1546 }
1547
1549 return (d_a < d_b) ? d_a : d_b;
1550 }
1551
1553 return (d_a > d_b) ? d_a : d_b;
1554 }
1555
1556 bool Ellipse::isInside(const util::Point &x) const {
1557 double u = 0., v = 0.;
1558 ellipseLocal(*this, x, u, v);
1559 return ellipseMetricInside(u, v, d_a, d_b);
1560 }
1561
1562 bool Ellipse::isOutside(const util::Point &x) const {
1563 return !isInside(x);
1564 }
1565
1566 bool Ellipse::isNear(const util::Point &x, const double &tol) const {
1567 double u = 0., v = 0.;
1568 ellipseLocal(*this, x, u, v);
1569 return ellipseMetricInside(u, v, d_a + tol, d_b + tol);
1570 }
1571
1572 bool Ellipse::isNearBoundary(const util::Point &x, const double &tol, const bool &within) const {
1573
1574 double u = 0., v = 0.;
1575 ellipseLocal(*this, x, u, v);
1576 const bool inOuter = ellipseMetricInside(u, v, d_a + tol, d_b + tol);
1577 const bool inInner = ellipseMetricInside(u, v, d_a - tol, d_b - tol);
1578 if (within)
1579 return inOuter && !inInner;
1580 return inOuter;
1581 }
1582
1584 return isNearBoundary(x, 1.0E-8, false);
1585 }
1586
1587 bool Ellipse::isInside(const std::pair<util::Point, util::Point> &box) const {
1588
1589 for (auto p : geom::getCornerPoints(2, box))
1590 if (!this->isInside(p))
1591 return false;
1592
1593 return true;
1594 }
1595
1596 bool Ellipse::isOutside(const std::pair<util::Point, util::Point> &box) const {
1597
1598 bool intersect = false;
1599 for (auto p : geom::getCornerPoints(2, box))
1600 if (!intersect)
1601 intersect = this->isInside(p);
1602
1603 return !intersect;
1604 }
1605
1606 bool Ellipse::isNear(const std::pair<util::Point, util::Point> &box, const double &tol) const {
1607
1608 if (this->isInside(box))
1609 return true;
1610
1611 for (auto p : geom::getCornerPoints(2, box)) {
1612 if (isNear(p, tol))
1613 return true;
1614 }
1615
1616 auto dxc = geom::getCenter(2, box) - d_x;
1617 auto r = geom::inscribedRadiusInBox(2, box);
1618 if (util::isLess(dxc.length(), boundingRadius() + r + tol))
1619 return true;
1620
1622 return util::isLess(dxc.length(), boundingRadius() + r + tol);
1623 }
1624
1625 bool Ellipse::doesIntersect(const std::pair<util::Point, util::Point> &box) const {
1626
1627 for (auto p : geom::getCornerPoints(2, box))
1628 if (this->isInside(p))
1629 return true;
1630
1631 return false;
1632 }
1633
1634 std::string Ellipse::printStr(int nt, int lvl) const {
1635
1636 auto tabS = util::io::getTabS(nt);
1637
1638 std::ostringstream oss;
1639
1640 oss << tabS << "------- Ellipse --------" << std::endl << std::endl;
1641 oss << tabS << "Name = " << d_name << std::endl;
1642 oss << tabS << "Center = " << d_x.printStr(0, lvl) << std::endl;
1643 oss << tabS << "Semi-axes a, b = " << d_a << ", " << d_b << std::endl;
1644 oss << tabS << "Theta (rad) = " << d_theta << std::endl;
1645
1646 if (lvl > 0)
1647 oss << tabS << "Bounding box: "
1648 << util::io::printBoxStr(box(0.), nt + 1);
1649
1650 if (lvl == 0)
1651 oss << std::endl;
1652
1653 return oss.str();
1654 }
1655
1656}// Ellipse
1657
1658//
1659// Cylinder
1660//
1661namespace geom {
1662 double Cylinder::volume() const {
1663 return M_PI * d_r * d_r * d_l;
1664 }
1665
1667 return d_x;
1668 }
1669
1670 std::pair<util::Point, util::Point> Cylinder::box() const {
1671
1672 return box(0.);
1673 }
1674
1675 std::pair<util::Point, util::Point>
1676 Cylinder::box(const double &tol) const {
1677
1678 if (d_xa.length() < 1.0E-10)
1679 return {util::Point(), util::Point()};
1680
1681 auto xb = d_xBegin - tol * d_xa;
1682 auto xt = d_xBegin + (d_l + tol) * d_xa;
1683
1684 double r = d_r + tol;
1685
1686 return {xb - r, xt + r};
1687 }
1688
1690
1691 auto box = this->box();
1692
1693 return 0.5 * (box.second - box.first).length();
1694 }
1695
1697
1698 return 0.5 * std::sqrt(d_l * d_l + 4. * d_r * d_r);
1699 }
1700
1701 bool Cylinder::isInside(const util::Point &x) const {
1702
1703 auto dx = x - d_xBegin;
1704
1705 if (dx.length() < 1.0E-10)
1706 return true;
1707
1708 double dx_dot_xa = dx * d_xa;
1709 if (util::isLess(dx_dot_xa, 0.) or
1710 util::isGreater(dx_dot_xa, d_l))
1711 return false;
1712 else {
1713
1714 // project dx onto cross-section plane of cylinder
1715 auto dx_project = dx - dx_dot_xa * d_xa;
1716
1717 return !util::isGreater(dx_project.length(), d_r + 1.0E-12);
1718 }
1719 }
1720
1721 bool Cylinder::isOutside(const util::Point &x) const {
1722 return !isInside(x);
1723 }
1724
1725 bool Cylinder::isNear(const util::Point &x, const double
1726 &tol) const {
1727
1728 auto dx = x - d_xBegin;
1729
1730 if (dx.length() < tol)
1731 return true;
1732
1733 double dx_dot_xa = dx * d_xa;
1734 if (util::isLess(dx_dot_xa, -tol) or
1735 util::isGreater(dx_dot_xa, d_l + tol))
1736 return false;
1737 else {
1738
1739 // project dx onto cross-section plane of cylinder
1740 auto dx_project = dx - dx_dot_xa * d_xa;
1741
1742 return !util::isGreater(dx_project.length(), d_r + tol);
1743 }
1744 }
1745
1747 const double &tol, const bool
1748 &within) const {
1749
1750 auto dx = x - d_xBegin;
1751
1752 if (dx.length() < tol)
1753 return true;
1754
1755 double dx_dot_xa = dx * d_xa;
1756 if (util::isLess(dx_dot_xa, -tol) or
1757 util::isGreater(dx_dot_xa, tol) or
1758 util::isGreater(dx_dot_xa, d_l + tol) or
1759 util::isLess(dx_dot_xa, d_l - tol))
1760 return false;
1761 else {
1762
1763 // project dx onto cross-section plane of cylinder
1764 auto dx_project = dx - dx_dot_xa * d_xa;
1765
1766 return !(util::isLess(dx_project.length(), d_r - tol) or
1767 util::isGreater(dx_project.length(), d_r + tol));
1768 }
1769 }
1770
1772
1773 return isNearBoundary(x, 1.0E-8, false);
1774 }
1775
1777 const std::pair<util::Point, util::Point> &box) const {
1778
1779 for (auto p: geom::getCornerPoints(3, box))
1780 if (!this->isInside(p))
1781 return false;
1782
1783 return true;
1784 }
1785
1787 const std::pair<util::Point, util::Point> &box) const {
1788
1789 bool intersect = false;
1790 for (auto p: geom::getCornerPoints(3, box))
1791 if (!intersect)
1792 intersect = this->isInside(p);
1793
1794 return !intersect;
1795 }
1796
1798 const std::pair<util::Point, util::Point> &box,
1799 const double &tol) const {
1800
1801 return geom::areBoxesNear(this->box(), box, tol, 3);
1802 }
1803
1805 const std::pair<util::Point, util::Point> &box) const {
1806
1807 // need to check all four corner points
1808 for (auto p: geom::getCornerPoints(3, box))
1809 if (this->isInside(p))
1810 return true;
1811
1812 return false;
1813 }
1814
1815
1816 std::string Cylinder::printStr(int nt, int lvl) const {
1817
1818 auto tabS = util::io::getTabS(nt);
1819
1820 std::ostringstream oss;
1821
1822 oss << tabS << "------- Cylinder --------" << std::endl << std::endl;
1823 oss << tabS << "Name = " << d_name << std::endl;
1824 oss << tabS << "Center = " << d_xBegin.printStr(0, lvl) << std::endl;
1825 oss << tabS << "Axis = " << d_xa.printStr(0, lvl) << std::endl;
1826 oss << tabS << "Radius = " << d_r << std::endl;
1827 oss << tabS << "Center = " << d_x.printStr(0, 0) << std::endl;
1828
1829 if (lvl > 0)
1830 oss << tabS << "Bounding box: "
1831 << util::io::printBoxStr(box(0.), nt + 1);
1832
1833 if (lvl == 0)
1834 oss << std::endl;
1835
1836 return oss.str();
1837 }
1838
1839}// Cylinder
1840
1841//
1842// Sphere
1843//
1844namespace geom {
1845 double Sphere::volume() const {
1846 return 4. * M_PI * d_r * d_r * d_r / 3.;
1847 }
1848
1850 return d_x;
1851 }
1852
1853 std::pair<util::Point, util::Point> Sphere::box() const {
1854
1855 return box(0.);
1856 }
1857
1858 std::pair<util::Point, util::Point> Sphere::box(const
1859 double &tol) const {
1860 double r = d_r + tol;
1861
1862 return {
1863 util::Point(d_x.d_x - r, d_x.d_y - r, d_x.d_z - r),
1864 util::Point(d_x.d_x + r, d_x.d_y + r, d_x.d_z + r)
1865 };
1866 }
1867
1869
1870 return d_r;
1871 }
1872
1873 double Sphere::boundingRadius() const {
1874
1875 return d_r;
1876 }
1877
1878 bool Sphere::isInside(const util::Point &x) const {
1879
1880 return util::isLess(d_x.dist(x), d_r + 1.0E-12);
1881 }
1882
1883 bool Sphere::isOutside(const util::Point &x) const {
1884 return !isInside(x);
1885 }
1886
1888 const double &tol) const {
1889
1890 // translate to origin
1891 auto x0 = x - d_x;
1892
1893 return util::isLess(x0.length(), d_r + tol);
1894 }
1895
1897 const double &tol, const bool
1898 &within) const {
1899
1900 // check if particle is within the tolerance distance
1901 if (!isNear(x, within ? 0. : tol))
1902 return false;
1903
1904 // check if it is close enough to circumference
1905 auto x0 = x - d_x;
1906
1907 return util::isLess(x0.length(), d_r + tol) ||
1908 util::isLess(x0.length(), d_r - tol);
1909 }
1910
1911 bool Sphere::doesIntersect(const util::Point &x) const {
1912
1913 return isNearBoundary(x, 1.0E-8, false);
1914 }
1915
1917 const std::pair<util::Point, util::Point> &box) const {
1918
1919 for (auto p: geom::getCornerPoints(3, box))
1920 if (!this->isInside(p))
1921 return false;
1922
1923 return true;
1924 }
1925
1927 const std::pair<util::Point, util::Point> &box) const {
1928
1929 bool intersect = false;
1930 for (auto p: geom::getCornerPoints(3, box))
1931 if (!intersect)
1932 intersect = this->isInside(p);
1933
1934 return !intersect;
1935 }
1936
1938 const std::pair<util::Point, util::Point> &box,
1939 const double &tol) const {
1940
1941 if (this->isInside(box))
1942 return true;
1943
1944 // get corner points of box
1945 auto cp = geom::getCornerPoints(3, box);
1946
1947 for (auto p: cp) {
1948
1949 // check the distance of corner point with the center
1950 auto dx = p - d_x;
1951 if (util::isLess(dx.length(), d_r + tol))
1952 return true;
1953 }
1954
1955 // check center to center distance
1956 auto dxc = geom::getCenter(3, box) - d_x;
1957
1958 // check wrt inscribed circle
1959 auto r = geom::inscribedRadiusInBox(3, box);
1960 if (util::isLess(dxc.length(), d_r + r + tol))
1961 return true;
1962
1963 // check wrt circumscribed circle
1965 return util::isLess(dxc.length(), d_r + r + tol);
1966 }
1967
1969 const std::pair<util::Point, util::Point> &box) const {
1970
1971 // need to check all four corner points
1972 for (auto p: geom::getCornerPoints(3, box))
1973 if (this->isInside(p))
1974 return true;
1975
1976 return false;
1977 }
1978
1979 std::string Sphere::printStr(int nt, int lvl) const {
1980
1981 auto tabS = util::io::getTabS(nt);
1982
1983 std::ostringstream oss;
1984
1985 oss << tabS << "------- Sphere --------" << std::endl << std::endl;
1986 oss << tabS << "Name = " << d_name << std::endl;
1987 oss << tabS << "Center = " << d_x.printStr(0, lvl) << std::endl;
1988 oss << tabS << "Radius = " << d_r << std::endl;
1989
1990 if (lvl > 0)
1991 oss << tabS << "Bounding box: "
1992 << util::io::printBoxStr(box(0.), nt + 1);
1993
1994 if (lvl == 0)
1995 oss << std::endl;
1996
1997 return oss.str();
1998 }
1999
2000}// Sphere
2001
2002//
2003// Ellipsoid
2004//
2005namespace geom {
2006
2007 namespace {
2008
2009 static void quatMul(double aw, double ax, double ay, double az,
2010 double bw, double bx, double by, double bz,
2011 double &cw, double &cx, double &cy, double &cz) {
2012 cw = aw * bw - ax * bx - ay * by - az * bz;
2013 cx = aw * bx + ax * bw + ay * bz - az * by;
2014 cy = aw * by - ax * bz + ay * bw + az * bx;
2015 cz = aw * bz + ax * by - ay * bx + az * bw;
2016 }
2017
2018 static void axisAngleToQuatSafe(const util::Point &axisIn, double theta,
2019 double &w, double &x, double &y, double &z) {
2020 if (std::abs(theta) < 1.0e-30) {
2021 w = 1.;
2022 x = y = z = 0.;
2023 return;
2024 }
2025 const double L = axisIn.length();
2026 if (L < 1.0e-30) {
2027 w = 1.;
2028 x = y = z = 0.;
2029 return;
2030 }
2031 const util::Point k = axisIn / L;
2032 const double half = 0.5 * theta;
2033 w = std::cos(half);
2034 const double s = std::sin(half);
2035 x = s * k.d_x;
2036 y = s * k.d_y;
2037 z = s * k.d_z;
2038 }
2039
2040 static void quatToAxisAngle(double qw, double qx, double qy, double qz,
2041 util::Point &axis, double &theta) {
2042 double n = std::sqrt(qw * qw + qx * qx + qy * qy + qz * qz);
2043 if (n < 1.0e-30) {
2044 axis = util::Point(0., 0., 1.);
2045 theta = 0.;
2046 return;
2047 }
2048 qw /= n;
2049 qx /= n;
2050 qy /= n;
2051 qz /= n;
2052 if (qw < 0.) {
2053 qw = -qw;
2054 qx = -qx;
2055 qy = -qy;
2056 qz = -qz;
2057 }
2058 qw = std::max(-1., std::min(1., qw));
2059 theta = 2. * std::acos(qw);
2060 const double sv = std::sqrt(qx * qx + qy * qy + qz * qz);
2061 if (sv < 1.0e-15) {
2062 axis = util::Point(0., 0., 1.);
2063 theta = 0.;
2064 return;
2065 }
2066 axis = util::Point(qx / sv, qy / sv, qz / sv);
2067 }
2068
2069 void ellipsoidBodyCoords(const Ellipsoid &e, const util::Point &p, double R[9], double &v0,
2070 double &v1, double &v2) {
2072 v0 = R[0] * p.d_x + R[3] * p.d_y + R[6] * p.d_z;
2073 v1 = R[1] * p.d_x + R[4] * p.d_y + R[7] * p.d_z;
2074 v2 = R[2] * p.d_x + R[5] * p.d_y + R[8] * p.d_z;
2075 }
2076
2077 double ellipsoidMetric(const Ellipsoid &e, const util::Point &x, double R[9], const double &ra,
2078 const double &rb, const double &rc) {
2079 const util::Point p{x.d_x - e.d_x.d_x, x.d_y - e.d_x.d_y, x.d_z - e.d_x.d_z};
2080 double v0, v1, v2;
2081 ellipsoidBodyCoords(e, p, R, v0, v1, v2);
2082 const double dx = v0 / ra;
2083 const double dy = v1 / rb;
2084 const double dz = v2 / rc;
2085 return dx * dx + dy * dy + dz * dz;
2086 }
2087
2088 } // namespace
2089
2090 void Ellipsoid::transform(const util::Point &translation, const double &scale, const double &angle,
2091 const util::Point &axis, const util::Point *rotationPoint) {
2092 const util::Point c0 = d_x;
2093 d_a *= scale;
2094 d_b *= scale;
2095 d_c *= scale;
2096
2097 double ow, ox, oy, oz;
2098 axisAngleToQuatSafe(d_axis, d_theta, ow, ox, oy, oz);
2099
2100 double rw, rx, ry, rz;
2101 axisAngleToQuatSafe(axis, angle, rw, rx, ry, rz);
2102
2103 double nw, nx, ny, nz;
2104 quatMul(rw, rx, ry, rz, ow, ox, oy, oz, nw, nx, ny, nz);
2105
2106 quatToAxisAngle(nw, nx, ny, nz, d_axis, d_theta);
2107 d_x = mapSimilarity(c0, c0, translation, scale, angle, axis, rotationPoint);
2108 }
2109
2110 double Ellipsoid::volume() const {
2111 return (4. / 3.) * M_PI * d_a * d_b * d_c;
2112 }
2113
2115 return d_x;
2116 }
2117
2118 std::pair<util::Point, util::Point> Ellipsoid::box() const {
2119 return box(0.);
2120 }
2121
2122 std::pair<util::Point, util::Point> Ellipsoid::box(const double &tol) const {
2123
2124 const double ra = d_a + tol;
2125 const double rb = d_b + tol;
2126 const double rc = d_c + tol;
2127
2128 double R[9];
2129 ellipsoidRotationMatrix(*this, R);
2130
2131 const double hx =
2132 std::sqrt((ra * R[0]) * (ra * R[0]) + (rb * R[3]) * (rb * R[3]) + (rc * R[6]) * (rc * R[6]));
2133 const double hy =
2134 std::sqrt((ra * R[1]) * (ra * R[1]) + (rb * R[4]) * (rb * R[4]) + (rc * R[7]) * (rc * R[7]));
2135 const double hz =
2136 std::sqrt((ra * R[2]) * (ra * R[2]) + (rb * R[5]) * (rb * R[5]) + (rc * R[8]) * (rc * R[8]));
2137
2138 return {util::Point(d_x.d_x - hx, d_x.d_y - hy, d_x.d_z - hz),
2139 util::Point(d_x.d_x + hx, d_x.d_y + hy, d_x.d_z + hz)};
2140 }
2141
2143 double m = d_a;
2144 if (d_b < m)
2145 m = d_b;
2146 if (d_c < m)
2147 m = d_c;
2148 return m;
2149 }
2150
2152 double m = d_a;
2153 if (d_b > m)
2154 m = d_b;
2155 if (d_c > m)
2156 m = d_c;
2157 return m;
2158 }
2159
2160 bool Ellipsoid::isInside(const util::Point &x) const {
2161 if (d_a <= 0. || d_b <= 0. || d_c <= 0.)
2162 return false;
2163 double R[9];
2164 return util::isLess(ellipsoidMetric(*this, x, R, d_a, d_b, d_c), 1. + 1.0E-12);
2165 }
2166
2167 bool Ellipsoid::isOutside(const util::Point &x) const {
2168 return !isInside(x);
2169 }
2170
2171 bool Ellipsoid::isNear(const util::Point &x, const double &tol) const {
2172 if (d_a + tol <= 0. || d_b + tol <= 0. || d_c + tol <= 0.)
2173 return false;
2174 double R[9];
2175 return util::isLess(ellipsoidMetric(*this, x, R, d_a + tol, d_b + tol, d_c + tol),
2176 1. + 1.0E-12);
2177 }
2178
2179 bool Ellipsoid::isNearBoundary(const util::Point &x, const double &tol, const bool &within) const {
2180
2181 const bool inOuter = isNear(x, tol);
2182 Ellipsoid shrunk = *this;
2183 shrunk.d_a = d_a - tol;
2184 shrunk.d_b = d_b - tol;
2185 shrunk.d_c = d_c - tol;
2186 if (shrunk.d_a <= 0. || shrunk.d_b <= 0. || shrunk.d_c <= 0.)
2187 return inOuter;
2188 const bool inInner = shrunk.isInside(x);
2189 if (within)
2190 return inOuter && !inInner;
2191 return inOuter;
2192 }
2193
2195 return isNearBoundary(x, 1.0E-8, false);
2196 }
2197
2198 bool Ellipsoid::isInside(const std::pair<util::Point, util::Point> &box) const {
2199
2200 for (auto p : geom::getCornerPoints(3, box))
2201 if (!this->isInside(p))
2202 return false;
2203
2204 return true;
2205 }
2206
2207 bool Ellipsoid::isOutside(const std::pair<util::Point, util::Point> &box) const {
2208
2209 bool intersect = false;
2210 for (auto p : geom::getCornerPoints(3, box))
2211 if (!intersect)
2212 intersect = this->isInside(p);
2213
2214 return !intersect;
2215 }
2216
2217 bool Ellipsoid::isNear(const std::pair<util::Point, util::Point> &box, const double &tol) const {
2218
2219 return geom::areBoxesNear(this->box(), box, tol, 3);
2220 }
2221
2222 bool Ellipsoid::doesIntersect(const std::pair<util::Point, util::Point> &box) const {
2223
2224 for (auto p : geom::getCornerPoints(3, box))
2225 if (this->isInside(p))
2226 return true;
2227
2228 return false;
2229 }
2230
2231 std::string Ellipsoid::printStr(int nt, int lvl) const {
2232
2233 auto tabS = util::io::getTabS(nt);
2234
2235 std::ostringstream oss;
2236
2237 oss << tabS << "------- Ellipsoid --------" << std::endl << std::endl;
2238 oss << tabS << "Name = " << d_name << std::endl;
2239 oss << tabS << "Center = " << d_x.printStr(0, lvl) << std::endl;
2240 oss << tabS << "Semi-axes a, b, c = " << d_a << ", " << d_b << ", " << d_c << std::endl;
2241 if (std::abs(d_theta) > 1.0e-14)
2242 oss << tabS << "Rotation axis (unit) = " << d_axis.printStr(0, lvl) << ", theta = " << d_theta
2243 << std::endl;
2244
2245 if (lvl > 0)
2246 oss << tabS << "Bounding box: "
2247 << util::io::printBoxStr(box(0.), nt + 1);
2248
2249 if (lvl == 0)
2250 oss << std::endl;
2251
2252 return oss.str();
2253 }
2254
2255}// Ellipsoid
2256
double inscribedRadius() const override
Computes the radius of biggest circle/sphere completely within the object.
bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override
cons
util::Point center() const override
Computes the center of object.
double d_r
Radius.
double boundingRadius() const override
Computes the radius of smallest circle/sphere such that object can be fit into it.
util::Point d_x
Center.
bool doesIntersect(const util::Point &x) const override
Checks if point lies exactly on the boundary.
bool isInside(const util::Point &x) const override
Checks if point is inside this object.
std::pair< util::Point, util::Point > box() const override
Computes the bounding box of object.
bool isNear(const util::Point &x, const double &tol) const override
Checks if point is within given distance of this object.
std::string printStr(int nt, int lvl) const override
Returns the string containing printable information about the object.
bool isOutside(const util::Point &x) const override
Checks if point is outside of this object.
double volume() const override
Computes the volume (area in 2d, length in 1d) of object.
util::Point center() const override
Computes the center of object.
util::Point d_x
Center.
std::pair< util::Point, util::Point > box() const override
Computes the bounding box of object.
bool isOutside(const util::Point &x) const override
Checks if point is outside of this object.
std::string printStr(int nt, int lvl) const override
Returns the string containing printable information about the object.
double d_r
Radius of bounding circle.
bool isInside(const util::Point &x) const override
Checks if point is inside this object.
bool doesIntersect(const util::Point &x) const override
Checks if point lies exactly on the boundary.
std::vector< util::Point > d_vertices
Vertices.
double inscribedRadius() const override
Computes the radius of biggest circle/sphere completely within the object.
bool isNear(const util::Point &x, const double &tol) const override
Checks if point is within given distance of this object.
double boundingRadius() const override
Computes the radius of smallest circle/sphere such that object can be fit into it.
bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override
cons
double d_L
Edge length of cube.
double volume() const override
Computes the volume (area in 2d, length in 1d) of object.
double boundingRadius() const override
Computes the radius of smallest circle/sphere such that object can be fit into it.
double inscribedRadius() const override
Computes the radius of biggest circle/sphere completely within the object.
double d_Ly
Edge length of cuboid in y-direction.
bool isNear(const util::Point &x, const double &tol) const override
Checks if point is within given distance of this object.
bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override
cons
std::vector< util::Point > d_vertices
Vertices.
util::Point center() const override
Computes the center of object.
double d_Lx
Edge length of cuboid in x-direction.
double volume() const override
Computes the volume (area in 2d, length in 1d) of object.
std::pair< util::Point, util::Point > box() const override
Computes the bounding box of object.
util::Point d_x
Center.
bool isInside(const util::Point &x) const override
Checks if point is inside this object.
double d_r
Radius of bounding circle.
std::string printStr(int nt, int lvl) const override
Returns the string containing printable information about the object.
double d_Lz
Edge length of cuboid in z-direction.
bool doesIntersect(const util::Point &x) const override
Checks if point lies exactly on the boundary.
bool isOutside(const util::Point &x) const override
Checks if point is outside of this object.
std::pair< util::Point, util::Point > box() const override
Computes the bounding box of object.
double boundingRadius() const override
Computes the radius of smallest circle/sphere such that object can be fit into it.
bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override
cons
util::Point d_xBegin
Center point of cross-section at the beginning.
double volume() const override
Computes the volume (area in 2d, length in 1d) of object.
util::Point d_xa
Axis of cylinder (unit vector)
bool isInside(const util::Point &x) const override
Checks if point is inside this object.
double inscribedRadius() const override
Computes the radius of biggest circle/sphere completely within the object.
double d_l
Length.
bool isOutside(const util::Point &x) const override
Checks if point is outside of this object.
bool isNear(const util::Point &x, const double &tol) const override
Checks if point is within given distance of this object.
std::string printStr(int nt, int lvl) const override
Returns the string containing printable information about the object.
util::Point d_x
Center.
util::Point center() const override
Computes the center of object.
double d_r
Radius.
bool doesIntersect(const util::Point &x) const override
Checks if point lies exactly on the boundary.
std::string printStr(int nt, int lvl) const override
Returns the string containing printable information about the object.
util::Point center() const override
Computes the center of object.
double inscribedRadius() const override
Computes the radius of biggest circle/sphere completely within the object.
double volume() const override
Computes the volume (area in 2d, length in 1d) of object.
bool doesIntersect(const util::Point &x) const override
Checks if point lies exactly on the boundary.
util::Point d_a
Axis: defined as the vector pointing from center to the first vertex v3 v2.
double boundingRadius() const override
Computes the radius of smallest circle/sphere such that object can be fit into it.
std::pair< util::Point, util::Point > box() const override
Computes the bounding box of object.
std::vector< util::Point > d_vertices
Vertices.
bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override
cons
double d_r
Distance between center and the farthest vertex.
bool isOutside(const util::Point &x) const override
Checks if point is outside of this object.
bool isNear(const util::Point &x, const double &tol) const override
Checks if point is within given distance of this object.
double d_w
Half width of neck.
bool isInside(const util::Point &x) const override
Checks if point is inside this object.
util::Point d_x
Center.
Filled ellipse in the plane z = center.d_z, semi-axes in the xy plane.
bool isOutside(const util::Point &x) const override
Checks if point is outside of this object.
double boundingRadius() const override
Computes the radius of smallest circle/sphere such that object can be fit into it.
double volume() const override
Computes the volume (area in 2d, length in 1d) of object.
double d_a
Semi-axis along local x before rotation (in-plane)
util::Point center() const override
Computes the center of object.
util::Point d_x
Center.
double d_b
Semi-axis along local y before rotation (in-plane)
double d_theta
Counter-clockwise rotation about +z through the center (radians)
bool isNear(const util::Point &x, const double &tol) const override
Checks if point is within given distance of this object.
std::string printStr(int nt, int lvl) const override
Returns the string containing printable information about the object.
bool isInside(const util::Point &x) const override
Checks if point is inside this object.
bool doesIntersect(const util::Point &x) const override
Checks if point lies exactly on the boundary.
double inscribedRadius() const override
Computes the radius of biggest circle/sphere completely within the object.
bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override
Checks if point is within given distance of this object.
std::pair< util::Point, util::Point > box() const override
Computes the bounding box of object.
Ellipsoid: center , semi-axes in a body frame rotated from world by axis–angle (Rodrigues)....
double inscribedRadius() const override
Computes the radius of biggest circle/sphere completely within the object.
double volume() const override
Computes the volume (area in 2d, length in 1d) of object.
double d_theta
Rotation angle (radians) about d_axis.
bool doesIntersect(const util::Point &x) const override
Checks if point lies exactly on the boundary.
util::Point center() const override
Computes the center of object.
bool isOutside(const util::Point &x) const override
Checks if point is outside of this object.
std::pair< util::Point, util::Point > box() const override
Computes the bounding box of object.
double boundingRadius() const override
Computes the radius of smallest circle/sphere such that object can be fit into it.
util::Point d_x
bool isInside(const util::Point &x) const override
Checks if point is inside this object.
util::Point d_axis
Unit rotation axis (axis–angle); default (0,0,1) when .
bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override
Checks if point is within given distance of this object.
std::string printStr(int nt, int lvl) const override
Returns the string containing printable information about the object.
void transform(const util::Point &translation, const double &scale, const double &angle, const util::Point &axis, const util::Point *rotationPoint) override
Similarity about pivot (default: old center d_x), then rigid displacement = translation: ....
bool isNear(const util::Point &x, const double &tol) const override
Checks if point is within given distance of this object.
const std::string d_name
name of object
bool isInside(const util::Point &x) const override
Checks if point is inside this object.
bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override
cons
double boundingRadius() const override
Computes the radius of smallest circle/sphere such that object can be fit into it.
util::Point d_x
Center.
std::pair< util::Point, util::Point > box() const override
Computes the bounding box of object.
double volume() const override
Computes the volume (area in 2d, length in 1d) of object.
util::Point center() const override
Computes the center of object.
bool doesIntersect(const util::Point &x) const override
Checks if point lies exactly on the boundary.
bool isNear(const util::Point &x, const double &tol) const override
Checks if point is within given distance of this object.
util::Point d_a
Axis: defined as the vector pointing from center to the first vertex.
std::string printStr(int nt, int lvl) const override
Returns the string containing printable information about the object.
double d_r
Distance between center and the farthest vertex of hexagon.
double inscribedRadius() const override
Computes the radius of biggest circle/sphere completely within the object.
bool isOutside(const util::Point &x) const override
Checks if point is outside of this object.
std::vector< util::Point > d_vertices
Vertices.
double inscribedRadius() const override
Computes the radius of biggest circle/sphere completely within the object.
std::pair< util::Point, util::Point > box() const override
Computes the bounding box of object.
std::vector< util::Point > d_vertices
Vertices.
double d_r
Radius of bounding circle.
bool doesIntersect(const util::Point &x) const override
Checks if point lies exactly on the boundary.
bool wallContactQuery(const util::Point &x, WallContactHit &hit) const override
util::Point center() const override
Computes the center of object.
double volume() const override
Computes the volume (area in 2d, length in 1d) of object.
bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override
cons
double boundingRadius() const override
Computes the radius of smallest circle/sphere such that object can be fit into it.
bool isOutside(const util::Point &x) const override
Checks if point is outside of this object.
double d_L
Length of line.
bool isNear(const util::Point &x, const double &tol) const override
Checks if point is within given distance of this object.
bool isInside(const util::Point &x) const override
Checks if point is inside this object.
std::string printStr(int nt, int lvl) const override
Returns the string containing printable information about the object.
util::Point d_x
Center.
util::Point d_p
a point on the plane
util::Point d_n
outward (free-space) normal, stored unit-length when possible
bool wallContactQuery(const util::Point &x, WallContactHit &hit) const override
Closest-point / signed-gap query for analytical wall contact. Default: unsupported (returns false,...
std::vector< util::Point > d_vertices
Vertices.
double boundingRadius() const override
Computes the radius of smallest circle/sphere such that object can be fit into it.
std::string printStr(int nt, int lvl) const override
Returns the string containing printable information about the object.
bool doesIntersect(const util::Point &x) const override
Checks if point lies exactly on the boundary.
double volume() const override
Computes the volume (area in 2d, length in 1d) of object.
bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override
cons
std::pair< util::Point, util::Point > box() const override
Computes the bounding box of object.
double d_r
Radius of bounding circle.
double inscribedRadius() const override
Computes the radius of biggest circle/sphere completely within the object.
bool wallContactQuery(const util::Point &x, WallContactHit &hit) const override
bool isInside(const util::Point &x) const override
Checks if point is inside this object.
bool isOutside(const util::Point &x) const override
Checks if point is outside of this object.
double d_Lx
Edge length of a rectangle in x-direction.
double d_Ly
Edge length of a rectangle in y-direction.
util::Point center() const override
Computes the center of object.
bool isNear(const util::Point &x, const double &tol) const override
Checks if point is within given distance of this object.
util::Point d_x
Center.
bool doesIntersect(const util::Point &x) const override
Checks if point lies exactly on the boundary.
bool isInside(const util::Point &x) const override
Checks if point is inside this object.
std::string printStr(int nt, int lvl) const override
Returns the string containing printable information about the object.
double boundingRadius() const override
Computes the radius of smallest circle/sphere such that object can be fit into it.
double d_r
Radius.
util::Point d_x
Center.
double inscribedRadius() const override
Computes the radius of biggest circle/sphere completely within the object.
util::Point center() const override
Computes the center of object.
std::pair< util::Point, util::Point > box() const override
Computes the bounding box of object.
bool isNear(const util::Point &x, const double &tol) const override
Checks if point is within given distance of this object.
bool isOutside(const util::Point &x) const override
Checks if point is outside of this object.
bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override
cons
double volume() const override
Computes the volume (area in 2d, length in 1d) of object.
util::Point center() const override
Computes the center of object.
std::string printStr(int nt, int lvl) const override
Returns the string containing printable information about the object.
bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override
cons
std::pair< util::Point, util::Point > box() const override
Computes the bounding box of object.
bool isNear(const util::Point &x, const double &tol) const override
Checks if point is within given distance of this object.
double boundingRadius() const override
Computes the radius of smallest circle/sphere such that object can be fit into it.
bool isInside(const util::Point &x) const override
Checks if point is inside this object.
bool isOutside(const util::Point &x) const override
Checks if point is outside of this object.
std::vector< util::Point > d_vertices
Vertices.
bool doesIntersect(const util::Point &x) const override
Checks if point lies exactly on the boundary.
double inscribedRadius() const override
Computes the radius of biggest circle/sphere completely within the object.
util::Point d_x
Center.
double volume() const override
Computes the volume (area in 2d, length in 1d) of object.
double d_L
Edge length of a square.
double d_r
Radius of bounding circle.
util::Point d_x
Center.
std::vector< util::Point > d_vertices
Vertices.
util::Point center() const override
Computes the center of object.
bool isNear(const util::Point &x, const double &tol) const override
Checks if point is within given distance of this object.
std::pair< util::Point, util::Point > box() const override
Computes the bounding box of object.
bool isOutside(const util::Point &x) const override
Checks if point is outside of this object.
double inscribedRadius() const override
Computes the radius of biggest circle/sphere completely within the object.
bool doesIntersect(const util::Point &x) const override
Checks if point lies exactly on the boundary.
bool isInside(const util::Point &x) const override
Checks if point is inside this object.
double boundingRadius() const override
Computes the radius of smallest circle/sphere such that object can be fit into it.
bool isNearBoundary(const util::Point &x, const double &tol, const bool &within) const override
cons
double volume() const override
Computes the volume (area in 2d, length in 1d) of object.
std::string printStr(int nt, int lvl) const override
Returns the string containing printable information about the object.
double d_r
Distance between center and the farthest vertex of triangle.
std::string printErrMsg(const std::string &geom_type, const std::vector< double > &params, const std::vector< size_t > &num_params_needed)
static void quatToAxisAngle(double qw, double qx, double qy, double qz, util::Point &axis, double &theta)
static void axisAngleToQuatSafe(const util::Point &axisIn, double theta, double &w, double &x, double &y, double &z)
void ellipseLocal(const Ellipse &e, const util::Point &x, double &u, double &v)
static void quatMul(double aw, double ax, double ay, double az, double bw, double bx, double by, double bz, double &cw, double &cx, double &cy, double &cz)
void ellipsoidBodyCoords(const Ellipsoid &e, const util::Point &p, double R[9], double &v0, double &v1, double &v2)
bool ellipseMetricInside(double u, double v, double a, double b)
double ellipsoidMetric(const Ellipsoid &e, const util::Point &x, double R[9], const double &ra, const double &rb, const double &rc)
util::Point getCenter(size_t dim, const std::pair< util::Point, util::Point > &box)
Returns center point.
bool isPointInsideRectangle(util::Point x, double x_min, double x_max, double y_min, double y_max)
Checks if point is inside a rectangle.
util::Point mapSimilarity(const util::Point &x, const util::Point &pivotDefault, const util::Point &t, double scale, double angle, const util::Point &axis, const util::Point *rotationPoint)
Definition geomObjects.h:55
bool isPointInsideBox(util::Point x, size_t dim, const std::pair< util::Point, util::Point > &box)
Returns true if point is inside box.
double circumscribedRadiusInBox(size_t dim, const std::pair< util::Point, util::Point > &box)
Computes the radius of smallest circle/sphere which can have the box inside.
bool areBoxesNear(const std::pair< util::Point, util::Point > &b1, const std::pair< util::Point, util::Point > &b2, const double &tol, size_t dim)
Checks if given two boxes are within given distance from each other.
std::vector< util::Point > getCornerPoints(size_t dim, const std::pair< util::Point, util::Point > &box)
Returns all corner points in the box.
double triangleArea(const util::Point &x1, const util::Point &x2, const util::Point &x3)
Compute area of triangle.
void ellipsoidRotationMatrix(const Ellipsoid &e, double R[9])
Row-major 3×3 rotation from ellipsoid axis–angle (identity if ).
double inscribedRadiusInBox(size_t dim, const std::pair< util::Point, util::Point > &box)
Computes the radius of biggest circle/sphere completely within the object.
bool isPointInsideCuboid(util::Point x, util::Point x_lbb, util::Point x_rtf)
Checks if point is inside a cuboid.
Collection of methods and database related to reading and writing.
std::string printBoxStr(const std::pair< util::Point, util::Point > &box, int nt=print_default_tab)
Returns formatted string for output.
Definition io.h:231
std::string getTabS(int nt)
Returns tab spaces of a given size.
Definition io.h:82
std::string printStr(const T &msg, int nt=print_default_tab)
Returns formatted string for output.
Definition io.h:96
bool isGreater(const double &a, const double &b)
Returns true if a > b.
Definition function.cpp:17
bool isLess(const double &a, const double &b)
Returns true if a < b.
Definition function.cpp:22
util::Point rotate(const util::Point &p, const double &theta, const util::Point &axis)
Returns the vector after rotating by desired angle.
Result of a wall-contact query against a geom (analytical walls).
util::Point outward_n
A structure to represent 3d vectors.
Definition point.h:30
double d_y
the y coordinate
Definition point.h:36
std::string printStr(int nt=0, int lvl=0) const
Returns the string containing printable information about the object.
Definition point.h:94
double dist(const Point &b) const
Computes the distance between a given point from this point.
Definition point.h:146
double d_z
the z coordinate
Definition point.h:39
double dot(const Point &b) const
Computes the dot product of this vector with another point.
Definition point.h:138
double length() const
Computes the Euclidean length of the vector.
Definition point.h:124
double lengthSq() const
Computes the Euclidean length of the vector.
Definition point.h:130
double d_x
the x coordinate
Definition point.h:33