DOLFIN
DOLFIN C++ interface
CGALExactArithmetic.h
1// Copyright (C) 2016-2017 Benjamin Kehlet, August Johansson, and Anders Logg
2//
3// This file is part of DOLFIN.
4//
5// DOLFIN is free software: you can redistribute it and/or modify
6// it under the terms of the GNU Lesser General Public License as published by
7// the Free Software Foundation, either version 3 of the License, or
8// (at your option) any later version.
9//
10// DOLFIN is distributed in the hope that it will be useful,
11// but WITHOUT ANY WARRANTY; without even the implied warranty of
12// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13// GNU Lesser General Public License for more details.
14//
15// You should have received a copy of the GNU Lesser General Public License
16// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
17//
18// First added: 2016-05-03
19// Last changed: 2017-10-26
20//
21// Developer note:
22//
23// This file contains reference implementations of collision detection
24// algorithms using exact arithmetic with CGAL. It is not included in
25// a normal build but is used as a reference for verification and
26// debugging of the inexact DOLFIN collision detection algorithms.
27// To enable, set the option DOLFIN_ENABLE_GEOMETRY_DEBUGGING when
28// configuring DOLFIN
29
30#ifndef __CGAL_EXACT_ARITHMETIC_H
31#define __CGAL_EXACT_ARITHMETIC_H
32
33#ifndef DOLFIN_ENABLE_GEOMETRY_DEBUGGING
34
35// Comparison macro just bypasses CGAL and test when not enabled
36#define CHECK_CGAL(RESULT_DOLFIN, RESULT_CGAL) RESULT_DOLFIN
37
38#define CGAL_INTERSECTION_CHECK(RESULT_DOLFIN, RESULT_CGAL) RESULT_DOLFIN
39
40#else
41
42#define CGAL_CHECK_TOLERANCE 1e-10
43
44#include "Point.h"
45#include "predicates.h"
46#include <dolfin/log/log.h>
47#include <dolfin/log/LogStream.h>
48#include <dolfin/math/basic.h>
49#include <vector>
50#include <algorithm>
51#include <sstream>
52#include <iomanip>
53
54// Check that results from DOLFIN and CGAL match
55namespace dolfin
56{
57 //---------------------------------------------------------------------------
58 // Functions to compare results between DOLFIN and CGAL
59 //---------------------------------------------------------------------------
60 inline bool
61 check_cgal(bool result_dolfin,
62 bool result_cgal,
63 const std::string& function)
64 {
65 if (result_dolfin != result_cgal)
66 {
67 // Convert results to strings
68 std::stringstream s_dolfin;
69 std::stringstream s_cgal;
70 s_dolfin << result_dolfin;
71 s_cgal << result_cgal;
72
73 // Issue error
74 dolfin_error("CGALExactArithmetic.h",
75 "verify geometric predicate with exact types",
76 "Error in predicate %s\n DOLFIN: %s\n CGAL: %s",
77 function.c_str(), s_dolfin.str().c_str(), s_cgal.str().c_str());
78 }
79
80 return result_dolfin;
81 }
82 //-----------------------------------------------------------------------------
83 inline
84 std::vector<Point>
85 cgal_intersection_check(const std::vector<Point>& dolfin_result,
86 const std::vector<Point>& cgal_result,
87 const std::string& function)
88 {
89 if (dolfin_result.size() != cgal_result.size())
90 {
91 dolfin_error("CGALExactArithmetic.h",
92 "verify intersection",
93 "Intersection function %s and CGAL give different size of point sets (%d vs %d)",
94 function.c_str(), dolfin_result.size(), cgal_result.size());
95 }
96
97 for (const Point& p1 : dolfin_result)
98 {
99 bool found = false;
100 for (const Point& p2 : cgal_result)
101 {
102 if ( (p1-p2).norm() < 1e-15 )
103 {
104 found = true;
105 break;
106 }
107 }
108
109 if (!found)
110 dolfin_error("CGALExactArithmetic.h",
111 "verify intersection construction result",
112 "Error in intersection function %s\nPoint (%f, %f, %f) in dolfin result not found",
113 function.c_str(), p1[0], p1[1], p1[2]);
114 }
115 return dolfin_result;
116 }
117} // end namespace dolfin
118//-----------------------------------------------------------------------------
119// Comparison macro that calls comparison function
120#define CHECK_CGAL(RESULT_DOLFIN, RESULT_CGAL) \
121 check_cgal(RESULT_DOLFIN, RESULT_CGAL, __FUNCTION__)
122
123#define CGAL_INTERSECTION_CHECK(RESULT_DOLFIN, RESULT_CGAL) \
124 cgal_intersection_check(RESULT_DOLFIN, RESULT_CGAL, __FUNCTION__)
125
126// CGAL includes
127#define CGAL_HEADER_ONLY
128#include <CGAL/Cartesian.h>
129#include <CGAL/Quotient.h>
130#include <CGAL/MP_Float.h>
131#include <CGAL/Point_2.h>
132#include <CGAL/Triangle_2.h>
133#include <CGAL/Segment_2.h>
134#include <CGAL/Point_3.h>
135#include <CGAL/Triangle_3.h>
136#include <CGAL/Segment_3.h>
137#include <CGAL/Tetrahedron_3.h>
138#include <CGAL/Polyhedron_3.h>
139#include <CGAL/convex_hull_3.h>
140#include <CGAL/intersections.h>
141#include <CGAL/intersection_of_Polyhedra_3.h>
142#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
143#include <CGAL/Triangulation_2.h>
144#include <CGAL/Nef_polyhedron_3.h>
145
146namespace
147{
148 // CGAL typedefs
149 /* typedef CGAL::Quotient<CGAL::MP_Float> ExactNumber; */
150 /* typedef CGAL::Cartesian<ExactNumber> ExactKernel; */
151 typedef CGAL::Exact_predicates_exact_constructions_kernel ExactKernel;
152 typedef ExactKernel::FT ExactNumber;
153
154 typedef ExactKernel::Point_2 Point_2;
155 typedef ExactKernel::Triangle_2 Triangle_2;
156 typedef ExactKernel::Segment_2 Segment_2;
157 typedef ExactKernel::Intersect_2 Intersect_2;
158 typedef ExactKernel::Point_3 Point_3;
159 typedef ExactKernel::Vector_3 Vector_3;
160 typedef ExactKernel::Triangle_3 Triangle_3;
161 typedef ExactKernel::Segment_3 Segment_3;
162 typedef ExactKernel::Tetrahedron_3 Tetrahedron_3;
163 typedef ExactKernel::Intersect_3 Intersect_3;
164 typedef CGAL::Nef_polyhedron_3<ExactKernel> Nef_polyhedron_3;
165 typedef CGAL::Polyhedron_3<ExactKernel> Polyhedron_3;
166 typedef CGAL::Triangulation_2<ExactKernel> Triangulation_2;
167
168 //---------------------------------------------------------------------------
169 // CGAL utility functions
170 //---------------------------------------------------------------------------
171 inline Point_2 convert_to_cgal_2d(double a, double b)
172 {
173 return Point_2(a, b);
174 }
175 //-----------------------------------------------------------------------------
176 inline Point_3 convert_to_cgal_3d(double a, double b, double c)
177 {
178 return Point_3(a, b, c);
179 }
180 //-----------------------------------------------------------------------------
181 inline Point_2 convert_to_cgal_2d(const dolfin::Point& p)
182 {
183 return Point_2(p[0], p[1]);
184 }
185 //-----------------------------------------------------------------------------
186 inline Point_3 convert_to_cgal_3d(const dolfin::Point& p)
187 {
188 return Point_3(p[0], p[1], p[2]);
189 }
190 //-----------------------------------------------------------------------------
191 inline Segment_2 convert_to_cgal_2d(const dolfin::Point& a,
192 const dolfin::Point& b)
193 {
194 return Segment_2(convert_to_cgal_2d(a), convert_to_cgal_2d(b));
195 }
196 //-----------------------------------------------------------------------------
197 inline Segment_3 convert_to_cgal_3d(const dolfin::Point& a,
198 const dolfin::Point& b)
199 {
200 return Segment_3(convert_to_cgal_3d(a), convert_to_cgal_3d(b));
201 }
202 //-----------------------------------------------------------------------------
203 inline Triangle_2 convert_to_cgal_2d(const dolfin::Point& a,
204 const dolfin::Point& b,
205 const dolfin::Point& c)
206 {
207 return Triangle_2(convert_to_cgal_2d(a),
208 convert_to_cgal_2d(b),
209 convert_to_cgal_2d(c));
210 }
211 //-----------------------------------------------------------------------------
212 inline Triangle_3 convert_to_cgal_3d(const dolfin::Point& a,
213 const dolfin::Point& b,
214 const dolfin::Point& c)
215 {
216 return Triangle_3(convert_to_cgal_3d(a),
217 convert_to_cgal_3d(b),
218 convert_to_cgal_3d(c));
219 }
220 //-----------------------------------------------------------------------------
221 inline Tetrahedron_3 convert_to_cgal_3d(const dolfin::Point& a,
222 const dolfin::Point& b,
223 const dolfin::Point& c,
224 const dolfin::Point& d)
225 {
226 return Tetrahedron_3(convert_to_cgal_3d(a),
227 convert_to_cgal_3d(b),
228 convert_to_cgal_3d(c),
229 convert_to_cgal_3d(d));
230 }
231 //-----------------------------------------------------------------------------
232 inline bool is_degenerate_2d(const dolfin::Point& a,
233 const dolfin::Point& b)
234 {
235 const Segment_2 s(convert_to_cgal_2d(a), convert_to_cgal_2d(b));
236 return s.is_degenerate();
237 }
238 //-----------------------------------------------------------------------------
239 inline bool is_degenerate_3d(const dolfin::Point& a,
240 const dolfin::Point& b)
241 {
242 const Segment_3 s(convert_to_cgal_3d(a), convert_to_cgal_3d(b));
243 return s.is_degenerate();
244 }
245 //-----------------------------------------------------------------------------
246 inline bool is_degenerate_2d(const dolfin::Point& a,
247 const dolfin::Point& b,
248 const dolfin::Point& c)
249 {
250 const Triangle_2 t(convert_to_cgal_2d(a),
251 convert_to_cgal_2d(b),
252 convert_to_cgal_2d(c));
253 return t.is_degenerate();
254 }
255 //-----------------------------------------------------------------------------
256 inline bool is_degenerate_3d(const dolfin::Point& a,
257 const dolfin::Point& b,
258 const dolfin::Point& c)
259 {
260 const Triangle_3 t(convert_to_cgal_3d(a),
261 convert_to_cgal_3d(b),
262 convert_to_cgal_3d(c));
263 return t.is_degenerate();
264 }
265 //-----------------------------------------------------------------------------
266 inline bool is_degenerate_3d(const dolfin::Point& a,
267 const dolfin::Point& b,
268 const dolfin::Point& c,
269 const dolfin::Point& d)
270 {
271 const Tetrahedron_3 t(convert_to_cgal_3d(a),
272 convert_to_cgal_3d(b),
273 convert_to_cgal_3d(c),
274 convert_to_cgal_3d(d));
275 return t.is_degenerate();
276 }
277 //-----------------------------------------------------------------------------
278 inline dolfin::Point convert_from_cgal(const Point_2& p)
279 {
280 return dolfin::Point(CGAL::to_double(p.x()),CGAL::to_double(p.y()));
281 }
282 //-----------------------------------------------------------------------------
283 inline dolfin::Point convert_from_cgal(const Point_3& p)
284 {
285 return dolfin::Point(CGAL::to_double(p.x()),
286 CGAL::to_double(p.y()),
287 CGAL::to_double(p.z()));
288 }
289 //-----------------------------------------------------------------------------
290 inline std::vector<dolfin::Point> convert_from_cgal(const Segment_2& s)
291 {
292 const std::vector<dolfin::Point> triangulation =
293 {{ dolfin::Point(CGAL::to_double(s.vertex(0)[0]),
294 CGAL::to_double(s.vertex(0)[1])),
295 dolfin::Point(CGAL::to_double(s.vertex(1)[0]),
296 CGAL::to_double(s.vertex(1)[1]))
297 }};
298 return triangulation;
299 }
300 //-----------------------------------------------------------------------------
301 inline std::vector<dolfin::Point> convert_from_cgal(const Segment_3& s)
302 {
303 const std::vector<dolfin::Point> triangulation =
304 {{ dolfin::Point(CGAL::to_double(s.vertex(0)[0]),
305 CGAL::to_double(s.vertex(0)[1]),
306 CGAL::to_double(s.vertex(0)[2])),
307 dolfin::Point(CGAL::to_double(s.vertex(1)[0]),
308 CGAL::to_double(s.vertex(1)[1]),
309 CGAL::to_double(s.vertex(1)[2]))
310 }};
311 return triangulation;
312 }
313 //-----------------------------------------------------------------------------
314 inline std::vector<dolfin::Point> convert_from_cgal(const Triangle_2& t)
315 {
316 const std::vector<dolfin::Point> triangulation =
317 {{ dolfin::Point(CGAL::to_double(t.vertex(0)[0]),
318 CGAL::to_double(t.vertex(0)[1])),
319 dolfin::Point(CGAL::to_double(t.vertex(2)[0]),
320 CGAL::to_double(t.vertex(2)[1])),
321 dolfin::Point(CGAL::to_double(t.vertex(1)[0]),
322 CGAL::to_double(t.vertex(1)[1]))
323 }};
324 return triangulation;
325 }
326 //-----------------------------------------------------------------------------
327 inline std::vector<dolfin::Point> convert_from_cgal(const Triangle_3& t)
328 {
329 const std::vector<dolfin::Point> triangulation =
330 {{ dolfin::Point(CGAL::to_double(t.vertex(0)[0]),
331 CGAL::to_double(t.vertex(0)[1]),
332 CGAL::to_double(t.vertex(0)[2])),
333 dolfin::Point(CGAL::to_double(t.vertex(2)[0]),
334 CGAL::to_double(t.vertex(2)[1]),
335 CGAL::to_double(t.vertex(2)[2])),
336 dolfin::Point(CGAL::to_double(t.vertex(1)[0]),
337 CGAL::to_double(t.vertex(1)[1]),
338 CGAL::to_double(t.vertex(1)[2]))
339 }};
340 return triangulation;
341 }
342 //-----------------------------------------------------------------------------
343 inline
344 std::vector<std::vector<dolfin::Point>>
345 triangulate_polygon_2d(const std::vector<dolfin::Point>& points)
346 {
347 // Convert points
348 std::vector<Point_2> pcgal(points.size());
349 for (std::size_t i = 0; i < points.size(); ++i)
350 pcgal[i] = convert_to_cgal_2d(points[i]);
351
352 // Triangulate
353 Triangulation_2 tcgal;
354 tcgal.insert(pcgal.begin(), pcgal.end());
355
356 // Convert back
357 std::vector<std::vector<dolfin::Point>> t;
358 for (Triangulation_2::Finite_faces_iterator fit = tcgal.finite_faces_begin();
359 fit != tcgal.finite_faces_end(); ++fit)
360 {
361 t.push_back({{ convert_from_cgal(tcgal.triangle(fit)[0]),
362 convert_from_cgal(tcgal.triangle(fit)[1]),
363 convert_from_cgal(tcgal.triangle(fit)[2]) }});
364 }
365
366 return t;
367 }
368 //-----------------------------------------------------------------------------
369 inline
370 std::vector<std::vector<dolfin::Point>>
371 triangulate_polygon_3d(const std::vector<dolfin::Point>& points)
372 {
373 // FIXME
374 dolfin::dolfin_error("CGALExactArithmetic.h",
375 "triangulate_polygon_3d",
376 "Not implemented");
377 return std::vector<std::vector<dolfin::Point>>();
378 }
379 //-----------------------------------------------------------------------------
380}
381
382namespace dolfin
383{
384 //---------------------------------------------------------------------------
385 // Reference implementations of DOLFIN collision detection predicates
386 // using CGAL exact arithmetic
387 // ---------------------------------------------------------------------------
388 inline bool cgal_collides_segment_point_2d(const Point& q0,
389 const Point& q1,
390 const Point& p,
391 bool only_interior=false)
392 {
393 const Point_2 q0_ = convert_to_cgal_2d(q0);
394 const Point_2 q1_ = convert_to_cgal_2d(q1);
395 const Point_2 p_ = convert_to_cgal_2d(p);
396
397 const bool intersects = CGAL::do_intersect(Segment_2(q0_, q1_), p_);
398 return only_interior ? intersects && p_ != q0_ && p_ != q1_ : intersects;
399 }
400 //-----------------------------------------------------------------------------
401 inline bool cgal_collides_segment_point_3d(const Point& q0,
402 const Point& q1,
403 const Point& p,
404 bool only_interior=false)
405 {
406 const Point_3 q0_ = convert_to_cgal_3d(q0);
407 const Point_3 q1_ = convert_to_cgal_3d(q1);
408 const Point_3 p_ = convert_to_cgal_3d(p);
409
410 const Segment_3 segment(q0_, q1_);
411 const bool intersects = segment.has_on(p_);
412 return only_interior ? intersects && p_ != q0_ && p_ != q1_ : intersects;
413 }
414 //-----------------------------------------------------------------------------
415 inline bool cgal_collides_segment_segment_2d(const Point& p0,
416 const Point& p1,
417 const Point& q0,
418 const Point& q1)
419 {
420 return CGAL::do_intersect(convert_to_cgal_2d(p0, p1),
421 convert_to_cgal_2d(q0, q1));
422 }
423 //-----------------------------------------------------------------------------
424 inline bool cgal_collides_segment_segment_3d(const Point& p0,
425 const Point& p1,
426 const Point& q0,
427 const Point& q1)
428 {
429 return CGAL::do_intersect(convert_to_cgal_3d(p0, p1),
430 convert_to_cgal_3d(q0, q1));
431 }
432 //-----------------------------------------------------------------------------
433 inline bool cgal_collides_triangle_point_2d(const Point& p0,
434 const Point& p1,
435 const Point& p2,
436 const Point &point)
437 {
438 return CGAL::do_intersect(convert_to_cgal_2d(p0, p1, p2),
439 convert_to_cgal_2d(point));
440 }
441 //-----------------------------------------------------------------------------
442 inline bool cgal_collides_triangle_point_3d(const Point& p0,
443 const Point& p1,
444 const Point& p2,
445 const Point &point)
446 {
447 const Triangle_3 tri = convert_to_cgal_3d(p0, p1, p2);
448 return tri.has_on(convert_to_cgal_3d(point));
449 }
450 //-----------------------------------------------------------------------------
451 inline bool cgal_collides_triangle_segment_2d(const Point& p0,
452 const Point& p1,
453 const Point& p2,
454 const Point& q0,
455 const Point& q1)
456 {
457 return CGAL::do_intersect(convert_to_cgal_2d(p0, p1, p2),
458 convert_to_cgal_2d(q0, q1));
459 }
460 //-----------------------------------------------------------------------------
461 inline bool cgal_collides_triangle_segment_3d(const Point& p0,
462 const Point& p1,
463 const Point& p2,
464 const Point& q0,
465 const Point& q1)
466 {
467 return CGAL::do_intersect(convert_to_cgal_3d(p0, p1, p2),
468 convert_to_cgal_3d(q0, q1));
469 }
470 //-----------------------------------------------------------------------------
471 inline bool cgal_collides_triangle_triangle_2d(const Point& p0,
472 const Point& p1,
473 const Point& p2,
474 const Point& q0,
475 const Point& q1,
476 const Point& q2)
477 {
478 return CGAL::do_intersect(convert_to_cgal_2d(p0, p1, p2),
479 convert_to_cgal_2d(q0, q1, q2));
480 }
481 //-----------------------------------------------------------------------------
482 inline bool cgal_collides_triangle_triangle_3d(const Point& p0,
483 const Point& p1,
484 const Point& p2,
485 const Point& q0,
486 const Point& q1,
487 const Point& q2)
488 {
489 return CGAL::do_intersect(convert_to_cgal_3d(p0, p1, p2),
490 convert_to_cgal_3d(q0, q1, q2));
491 }
492 //-----------------------------------------------------------------------------
493 inline bool cgal_collides_tetrahedron_point_3d(const Point& p0,
494 const Point& p1,
495 const Point& p2,
496 const Point& p3,
497 const Point& q0)
498 {
499 const Tetrahedron_3 tet = convert_to_cgal_3d(p0, p1, p2, p3);
500 return !tet.has_on_unbounded_side(convert_to_cgal_3d(q0));
501 }
502 //-----------------------------------------------------------------------------
503 inline bool cgal_collides_tetrahedron_segment_3d(const Point& p0,
504 const Point& p1,
505 const Point& p2,
506 const Point& p3,
507 const Point& q0,
508 const Point& q1)
509 {
510 if (cgal_collides_tetrahedron_point_3d(p0, p1, p2, p3, q0) or
511 cgal_collides_tetrahedron_point_3d(p0, p1, p2, p3, q1))
512 return true;
513
514 if (cgal_collides_triangle_segment_3d(p0, p1, p2, q0, q1) or
515 cgal_collides_triangle_segment_3d(p0, p2, p3, q0, q1) or
516 cgal_collides_triangle_segment_3d(p0, p3, p1, q0, q1) or
517 cgal_collides_triangle_segment_3d(p1, p3, p2, q0, q1))
518 return true;
519
520 return false;
521 }
522 //-----------------------------------------------------------------------------
523 inline bool cgal_collides_tetrahedron_triangle_3d(const Point& p0,
524 const Point& p1,
525 const Point& p2,
526 const Point& p3,
527 const Point& q0,
528 const Point& q1,
529 const Point& q2)
530 {
531 return CGAL::do_intersect(convert_to_cgal_3d(p0, p1, p2, p3),
532 convert_to_cgal_3d(q0, q1, q2));
533 }
534 //-----------------------------------------------------------------------------
535 inline bool cgal_collides_tetrahedron_tetrahedron_3d(const Point& p0,
536 const Point& p1,
537 const Point& p2,
538 const Point& p3,
539 const Point& q0,
540 const Point& q1,
541 const Point& q2,
542 const Point& q3)
543 {
544 // Check volume collisions
545 if (cgal_collides_tetrahedron_point_3d(p0, p1, p2, p3, q0)) return true;
546 if (cgal_collides_tetrahedron_point_3d(p0, p1, p2, p3, q1)) return true;
547 if (cgal_collides_tetrahedron_point_3d(p0, p1, p2, p3, q2)) return true;
548 if (cgal_collides_tetrahedron_point_3d(p0, p1, p2, p3, q3)) return true;
549 if (cgal_collides_tetrahedron_point_3d(q0, q1, q2, q3, p0)) return true;
550 if (cgal_collides_tetrahedron_point_3d(q0, q1, q2, q3, p1)) return true;
551 if (cgal_collides_tetrahedron_point_3d(q0, q1, q2, q3, p2)) return true;
552 if (cgal_collides_tetrahedron_point_3d(q0, q1, q2, q3, p3)) return true;
553
554 Polyhedron_3 tet_a;
555 tet_a.make_tetrahedron(convert_to_cgal_3d(p0),
556 convert_to_cgal_3d(p1),
557 convert_to_cgal_3d(p2),
558 convert_to_cgal_3d(p3));
559
560 Polyhedron_3 tet_b;
561 tet_b.make_tetrahedron(convert_to_cgal_3d(q0),
562 convert_to_cgal_3d(q1),
563 convert_to_cgal_3d(q2),
564 convert_to_cgal_3d(q3));
565
566 // Check for polyhedron intersection (recall that a polyhedron is
567 // only its vertices, edges and faces)
568 std::size_t cnt = 0;
569 CGAL::Counting_output_iterator out(&cnt);
570 CGAL::intersection_Polyhedron_3_Polyhedron_3<Polyhedron_3>(tet_a,
571 tet_b,
572 out);
573 // The tetrahedra does not intersect if cnt == 0
574 return cnt != 0;
575 }
576 //----------------------------------------------------------------------------
577 // Reference implementations of DOLFIN intersection triangulation
578 // functions using CGAL with exact arithmetic
579 // ---------------------------------------------------------------------------
580 inline
581 std::vector<Point> cgal_intersection_segment_segment_2d(const Point& p0,
582 const Point& p1,
583 const Point& q0,
584 const Point& q1)
585 {
586 dolfin_assert(!is_degenerate_2d(p0, p1));
587 dolfin_assert(!is_degenerate_2d(q0, q1));
588
589 const auto I0 = convert_to_cgal_2d(p0, p1);
590 const auto I1 = convert_to_cgal_2d(q0, q1);
591
592 if (const auto ii = CGAL::intersection(I0, I1))
593 {
594 if (const Point_2* p = boost::get<Point_2>(&*ii))
595 {
596 return std::vector<Point>{convert_from_cgal(*p)};
597 }
598 else if (const Segment_2* s = boost::get<Segment_2>(&*ii))
599 {
600 return convert_from_cgal(*s);
601 }
602 else
603 {
604 dolfin_error("CGALExactArithmetic.h",
605 "cgal_intersection_segment_segment_2d",
606 "Unexpected behavior");
607 }
608 }
609
610 return std::vector<Point>();
611 }
612 //-----------------------------------------------------------------------------
613 inline
614 std::vector<Point> cgal_intersection_segment_segment_3d(const Point& p0,
615 const Point& p1,
616 const Point& q0,
617 const Point& q1)
618 {
619 dolfin_assert(!is_degenerate_3d(p0, p1));
620 dolfin_assert(!is_degenerate_3d(q0, q1));
621
622 const auto I0 = convert_to_cgal_3d(p0, p1);
623 const auto I1 = convert_to_cgal_3d(q0, q1);
624
625 if (const auto ii = CGAL::intersection(I0, I1))
626 {
627 if (const Point_3* p = boost::get<Point_3>(&*ii))
628 {
629 return std::vector<Point>{convert_from_cgal(*p)};
630 }
631 else if (const Segment_3* s = boost::get<Segment_3>(&*ii))
632 {
633 return convert_from_cgal(*s);
634 }
635 else
636 {
637 dolfin_error("CGALExactArithmetic.h",
638 "cgal_intersection_segment_segment_3d",
639 "Unexpected behavior");
640 }
641 }
642
643 return std::vector<Point>();
644 }
645 //-----------------------------------------------------------------------------
646 inline
647 std::vector<Point> cgal_triangulate_segment_segment_2d(const Point& p0,
648 const Point& p1,
649 const Point& q0,
650 const Point& q1)
651 {
652 return cgal_intersection_segment_segment_2d(p0, p1, q0, q1);
653 }
654 //-----------------------------------------------------------------------------
655 inline
656 std::vector<Point> cgal_triangulate_segment_segment_3d(const Point& p0,
657 const Point& p1,
658 const Point& q0,
659 const Point& q1)
660 {
661 return cgal_intersection_segment_segment_3d(p0, p1, q0, q1);
662 }
663 //-----------------------------------------------------------------------------
664 inline
665 std::vector<Point> cgal_intersection_triangle_segment_2d(const Point& p0,
666 const Point& p1,
667 const Point& p2,
668 const Point& q0,
669 const Point& q1)
670 {
671 dolfin_assert(!is_degenerate_2d(p0, p1, p2));
672 dolfin_assert(!is_degenerate_2d(q0, q1));
673
674 const auto T = convert_to_cgal_2d(p0, p1, p2);
675 const auto I = convert_to_cgal_2d(q0, q1);
676
677 if (const auto ii = CGAL::intersection(T, I))
678 {
679 if (const Point_2* p = boost::get<Point_2>(&*ii))
680 {
681 return std::vector<Point>{convert_from_cgal(*p)};
682 }
683 else if (const Segment_2* s = boost::get<Segment_2>(&*ii))
684 {
685 return convert_from_cgal(*s);
686 }
687 else
688 {
689 dolfin_error("CGALExactArithmetic.h",
690 "cgal_intersection_triangle_segment_2d",
691 "Unexpected behavior");
692 }
693 }
694
695 return std::vector<Point>();
696 }
697 //-----------------------------------------------------------------------------
698 inline
699 std::vector<Point> cgal_intersection_triangle_segment_3d(const Point& p0,
700 const Point& p1,
701 const Point& p2,
702 const Point& q0,
703 const Point& q1)
704 {
705 dolfin_assert(!is_degenerate_3d(p0, p1, p2));
706 dolfin_assert(!is_degenerate_3d(q0, q1));
707
708 const auto T = convert_to_cgal_3d(p0, p1, p2);
709 const auto I = convert_to_cgal_3d(q0, q1);
710
711 if (const auto ii = CGAL::intersection(T, I))
712 {
713 if (const Point_3* p = boost::get<Point_3>(&*ii))
714 return std::vector<Point>{convert_from_cgal(*p)};
715 else if (const Segment_3* s = boost::get<Segment_3>(&*ii))
716 return convert_from_cgal(*s);
717 else
718 {
719 dolfin_error("CGALExactArithmetic.h",
720 "cgal_intersection_triangle_segment_3d",
721 "Unexpected behavior");
722 }
723 }
724
725 return std::vector<Point>();
726 }
727 //-----------------------------------------------------------------------------
728 inline
729 std::vector<Point> cgal_triangulate_triangle_segment_2d(const Point& p0,
730 const Point& p1,
731 const Point& p2,
732 const Point& q0,
733 const Point& q1)
734 {
735 return cgal_intersection_triangle_segment_2d(p0, p1, p2, q0, q1);
736 }
737 //-----------------------------------------------------------------------------
738 inline
739 std::vector<Point> cgal_triangulate_triangle_segment_3d(const Point& p0,
740 const Point& p1,
741 const Point& p2,
742 const Point& q0,
743 const Point& q1)
744 {
745 return cgal_intersection_triangle_segment_3d(p0, p1, p2, q0, q1);
746 }
747 //-----------------------------------------------------------------------------
748 inline
749 std::vector<Point> cgal_intersection_triangle_triangle_2d(const Point& p0,
750 const Point& p1,
751 const Point& p2,
752 const Point& q0,
753 const Point& q1,
754 const Point& q2)
755 {
756 dolfin_assert(!is_degenerate_2d(p0, p1, p2));
757 dolfin_assert(!is_degenerate_2d(q0, q1, q2));
758
759 const Triangle_2 T0 = convert_to_cgal_2d(p0, p1, p2);
760 const Triangle_2 T1 = convert_to_cgal_2d(q0, q1, q2);
761 std::vector<Point> intersection;
762
763 if (const auto ii = CGAL::intersection(T0, T1))
764 {
765 if (const Point_2* p = boost::get<Point_2>(&*ii))
766 {
767 intersection.push_back(convert_from_cgal(*p));
768 }
769 else if (const Segment_2* s = boost::get<Segment_2>(&*ii))
770 {
771 intersection = convert_from_cgal(*s);
772 }
773 else if (const Triangle_2* t = boost::get<Triangle_2>(&*ii))
774 {
775 intersection = convert_from_cgal(*t);;
776 }
777 else if (const std::vector<Point_2>* cgal_points = boost::get<std::vector<Point_2>>(&*ii))
778 {
779 for (Point_2 p : *cgal_points)
780 {
781 intersection.push_back(convert_from_cgal(p));
782 }
783 }
784 else
785 {
786 dolfin_error("CGALExactArithmetic.h",
787 "cgal_intersection_triangle_triangle_2d",
788 "Unexpected behavior");
789 }
790
791 // NB: the parsing can return triangulation of size 0, for example
792 // if it detected a triangle but it was found to be flat.
793 /* if (triangulation.size() == 0) */
794 /* dolfin_error("CGALExactArithmetic.h", */
795 /* "find intersection of two triangles in cgal_intersection_triangle_triangle function", */
796 /* "no intersection found"); */
797 }
798
799 return intersection;
800 }
801 //-----------------------------------------------------------------------------
802 inline
803 std::vector<Point> cgal_intersection_triangle_triangle_3d(const Point& p0,
804 const Point& p1,
805 const Point& p2,
806 const Point& q0,
807 const Point& q1,
808 const Point& q2)
809 {
810 dolfin_assert(!is_degenerate_3d(p0, p1, p2));
811 dolfin_assert(!is_degenerate_3d(q0, q1, q2));
812
813 const Triangle_3 T0 = convert_to_cgal_3d(p0, p1, p2);
814 const Triangle_3 T1 = convert_to_cgal_3d(q0, q1, q2);
815 std::vector<Point> intersection;
816
817 if (const auto ii = CGAL::intersection(T0, T1))
818 {
819 if (const Point_3* p = boost::get<Point_3>(&*ii))
820 {
821 intersection.push_back(convert_from_cgal(*p));
822 }
823 else if (const Segment_3* s = boost::get<Segment_3>(&*ii))
824 {
825 intersection = convert_from_cgal(*s);
826 }
827 else if (const Triangle_3* t = boost::get<Triangle_3>(&*ii))
828 {
829 intersection = convert_from_cgal(*t);;
830 }
831 else if (const std::vector<Point_3>* cgal_points = boost::get<std::vector<Point_3>>(&*ii))
832 {
833 for (Point_3 p : *cgal_points)
834 {
835 intersection.push_back(convert_from_cgal(p));
836 }
837 }
838 else
839 {
840 dolfin_error("CGALExactArithmetic.h",
841 "cgal_intersection_triangle_triangle_3d",
842 "Unexpected behavior");
843 }
844 }
845
846 return intersection;
847 }
848 //----------------------------------------------------------------------------
849 inline
850 std::vector<std::vector<Point>>
851 cgal_triangulate_triangle_triangle_2d(const Point& p0,
852 const Point& p1,
853 const Point& p2,
854 const Point& q0,
855 const Point& q1,
856 const Point& q2)
857 {
858 dolfin_assert(!is_degenerate_2d(p0, p1, p2));
859 dolfin_assert(!is_degenerate_2d(q0, q1, q2));
860
861 const std::vector<Point> intersection
862 = cgal_intersection_triangle_triangle_2d(p0, p1, p2, q0, q1, q2);
863
864 if (intersection.size() < 4)
865 {
866 return std::vector<std::vector<Point>>{intersection};
867 }
868 else
869 {
870 dolfin_assert(intersection.size() == 4 ||
871 intersection.size() == 5 ||
872 intersection.size() == 6);
873 return triangulate_polygon_2d(intersection);
874 }
875 }
876 //-----------------------------------------------------------------------------
877 inline
878 std::vector<std::vector<Point>>
879 cgal_triangulate_triangle_triangle_3d(const Point& p0,
880 const Point& p1,
881 const Point& p2,
882 const Point& q0,
883 const Point& q1,
884 const Point& q2)
885 {
886 dolfin_assert(!is_degenerate_3d(p0, p1, p2));
887 dolfin_assert(!is_degenerate_3d(q0, q1, q2));
888
889 const std::vector<Point> intersection
890 = cgal_intersection_triangle_triangle_3d(p0, p1, p2, q0, q1, q2);
891
892 if (intersection.size() < 4)
893 {
894 return std::vector<std::vector<Point>>{intersection};
895 }
896 else
897 {
898 dolfin_assert(intersection.size() == 4 ||
899 intersection.size() == 5 ||
900 intersection.size() == 6);
901 return triangulate_polygon_3d(intersection);
902 }
903 }
904 //-----------------------------------------------------------------------------
905 inline
906 std::vector<Point>
907 cgal_intersection_tetrahedron_triangle(const Point& p0,
908 const Point& p1,
909 const Point& p2,
910 const Point& p3,
911 const Point& q0,
912 const Point& q1,
913 const Point& q2)
914 {
915 dolfin_assert(!is_degenerate_3d(p0, p1, p2, p3));
916 dolfin_assert(!is_degenerate_3d(q0, q1, q2));
917
918 // const Tetrahedron_3 tet = convert_from_cgal(p0, p1, p2, p3);
919 // const Triangle_3 tri = convert_from_cgal(q0, q1, q2);
920
921 Polyhedron_3 tet;
922 tet.make_tetrahedron(convert_to_cgal_3d(p0),
923 convert_to_cgal_3d(p1),
924 convert_to_cgal_3d(p2),
925 convert_to_cgal_3d(p3));
926 Polyhedron_3 tri;
927 tri.make_triangle(convert_to_cgal_3d(q0),
928 convert_to_cgal_3d(q1),
929 convert_to_cgal_3d(q2));
930
931 std::list<std::vector<Point_3> > triangulation;
932 CGAL::intersection_Polyhedron_3_Polyhedron_3(tet,
933 tri,
934 std::back_inserter(triangulation));
935
936 // FIXME: do we need to add interior point checks? Maybe
937 // Polyhedron_3 is only top dim 2?
938
939 // Shouldn't get here
940 dolfin_error("CGALExactArithmetic.h",
941 "cgal_intersection_tetrahedron_triangle",
942 "Not implemented");
943
944 return std::vector<Point>();
945 }
946 //-----------------------------------------------------------------------------
947 inline std::vector<std::vector<Point>>
948 cgal_triangulate_tetrahedron_triangle(const Point& p0,
949 const Point& p1,
950 const Point& p2,
951 const Point& p3,
952 const Point& q0,
953 const Point& q1,
954 const Point& q2)
955 {
956 dolfin_assert(!is_degenerate_3d(p0, p1, p2, p3));
957 dolfin_assert(!is_degenerate_3d(q0, q1, q2));
958
959 std::vector<Point> intersection =
960 cgal_intersection_tetrahedron_triangle(p0, p1, p2, p3, q0, q1, q2);
961
962 // Shouldn't get here
963 dolfin_error("CGALExactArithmetic.h",
964 "cgal_triangulation_tetrahedron_triangle",
965 "Not implemented");
966
967 return std::vector<std::vector<Point>>();
968 }
969 //-----------------------------------------------------------------------------
970 inline std::vector<Point>
971 cgal_intersection_tetrahedron_tetrahedron_3d(const Point& p0,
972 const Point& p1,
973 const Point& p2,
974 const Point& p3,
975 const Point& q0,
976 const Point& q1,
977 const Point& q2,
978 const Point& q3)
979 {
980 dolfin_assert(!is_degenerate_3d(p0, p1, p2, p3));
981 dolfin_assert(!is_degenerate_3d(q0, q1, q2, q3));
982
983 Polyhedron_3 tet_a;
984 tet_a.make_tetrahedron(convert_to_cgal_3d(p0),
985 convert_to_cgal_3d(p1),
986 convert_to_cgal_3d(p2),
987 convert_to_cgal_3d(p3));
988 Polyhedron_3 tet_b;
989 tet_b.make_tetrahedron(convert_to_cgal_3d(q0),
990 convert_to_cgal_3d(q1),
991 convert_to_cgal_3d(q2),
992 convert_to_cgal_3d(q3));
993
994 const Nef_polyhedron_3 tet_a_nef(tet_a);
995 const Nef_polyhedron_3 tet_b_nef(tet_b);
996
997 const Nef_polyhedron_3 intersection_nef = tet_a_nef*tet_b_nef;
998
999 std::vector<Point> res;
1000
1001 for (auto vit = intersection_nef.vertices_begin();
1002 vit != intersection_nef.vertices_end(); ++vit)
1003 {
1004 res.push_back(Point(CGAL::to_double(vit->point().x()),
1005 CGAL::to_double(vit->point().y()),
1006 CGAL::to_double(vit->point().z())));
1007 }
1008
1009 return res;
1010 }
1011 //----------------------------------------------------------------------------
1012 // Reference implementations of DOLFIN is_degenerate
1013 //-----------------------------------------------------------------------------
1014 inline bool cgal_is_degenerate_2d(const std::vector<Point>& s)
1015 {
1016 if (s.size() < 2 or s.size() > 3)
1017 {
1018 info("Degenerate 2D simplex with %d vertices.", s.size());
1019 return true;
1020 }
1021
1022 switch (s.size())
1023 {
1024 case 2: return is_degenerate_2d(s[0], s[1]);
1025 case 3: return is_degenerate_2d(s[0], s[1], s[2]);
1026 }
1027
1028 // Shouldn't get here
1029 dolfin_error("CGALExactArithmetic.h",
1030 "call cgal_is_degenerate_2d",
1031 "Only implemented for simplices of tdim 0, 1 and 2, not tdim = %d", s.size() - 1);
1032
1033 return true;
1034 }
1035 //-----------------------------------------------------------------------------
1036 inline bool cgal_is_degenerate_3d(const std::vector<Point>& s)
1037 {
1038 if (s.size() < 2 or s.size() > 4)
1039 {
1040 info("Degenerate 3D simplex with %d vertices.", s.size());
1041 return true;
1042 }
1043
1044 switch (s.size())
1045 {
1046 case 2: return is_degenerate_3d(s[0], s[1]);
1047 case 3: return is_degenerate_3d(s[0], s[1], s[2]);
1048 case 4: return is_degenerate_3d(s[0], s[1], s[2], s[3]);
1049 }
1050
1051 // Shouldn't get here
1052 dolfin_error("CGALExactArithmetic.h",
1053 "call cgal_is_degenerate_3d",
1054 "Only implemented for simplices of tdim 0, 1, 2 and 3, not tdim = %d", s.size() - 1);
1055
1056 return true;
1057 }
1058 //-----------------------------------------------------------------------------
1059 // Computes the volume of the convex hull of the given points
1060 inline double cgal_polyhedron_volume(const std::vector<Point>& ch)
1061 {
1062 if (ch.size() < 4)
1063 return 0;
1064
1065 std::vector<Point_3> exact_points;
1066 exact_points.reserve(ch.size());
1067 for (const Point& p : ch)
1068 {
1069 exact_points.push_back(Point_3(p.x(), p.y(), p.z()));
1070 }
1071
1072 // Compute the convex hull as a cgal polyhedron_3
1073 Polyhedron_3 p;
1074 CGAL::convex_hull_3(exact_points.begin(), exact_points.end(), p);
1075
1076 ExactNumber volume = .0;
1077 for (Polyhedron_3::Facet_const_iterator it = p.facets_begin();
1078 it != p.facets_end(); it++)
1079 {
1080 const Polyhedron_3::Halfedge_const_handle h = it->halfedge();
1081 const Vector_3 V0 = h->vertex()->point()-CGAL::ORIGIN;
1082 const Vector_3 V1 = h->next()->vertex()->point()-CGAL::ORIGIN;
1083 const Vector_3 V2 = h->next()->next()->vertex()->point()-CGAL::ORIGIN;
1084
1085 volume += V0*CGAL::cross_product(V1, V2);
1086 }
1087
1088 return std::abs(CGAL::to_double(volume/6.0));
1089 }
1090 //-----------------------------------------------------------------------------
1091 inline double cgal_tet_volume(const std::vector<Point>& ch)
1092 {
1093 dolfin_assert(ch.size() == 3);
1094 return CGAL::to_double(CGAL::volume(Point_3(ch[0].x(), ch[0].y(), ch[0].z()),
1095 Point_3(ch[1].x(), ch[1].y(), ch[1].z()),
1096 Point_3(ch[2].x(), ch[2].y(), ch[2].z()),
1097 Point_3(ch[3].x(), ch[3].y(), ch[3].z())));
1098 }
1099 //-----------------------------------------------------------------------------
1100 inline bool cgal_tet_is_degenerate(const std::vector<Point>& t)
1101 {
1102 Tetrahedron_3 tet(Point_3(t[0].x(), t[0].y(), t[0].z()),
1103 Point_3(t[1].x(), t[1].y(), t[1].z()),
1104 Point_3(t[2].x(), t[2].y(), t[2].z()),
1105 Point_3(t[3].x(), t[3].y(), t[3].z()));
1106
1107 return tet.is_degenerate();
1108 }
1109
1110 //-----------------------------------------------------------------------------
1111 inline bool cgal_triangulation_has_degenerate(std::vector<std::vector<Point>> triangulation)
1112 {
1113 for (const std::vector<Point>& t : triangulation)
1114 {
1115 if (cgal_tet_is_degenerate(t))
1116 return true;
1117 }
1118
1119 return false;
1120 }
1121
1122 //-----------------------------------------------------------------------------
1123 inline bool cgal_triangulation_overlap(std::vector<std::vector<Point>> triangulation)
1124 {
1125 std::vector<Tetrahedron_3> tets;
1126
1127 for (const std::vector<Point>& t : triangulation)
1128 {
1129 Tetrahedron_3 tet(Point_3(t[0].x(), t[0].y(), t[0].z()),
1130 Point_3(t[1].x(), t[1].y(), t[1].z()),
1131 Point_3(t[2].x(), t[2].y(), t[2].z()),
1132 Point_3(t[3].x(), t[3].y(), t[3].z()));
1133
1134 for (const Tetrahedron_3& t0 : tets)
1135 {
1136 for (int i = 0; i < 4; i++)
1137 {
1138 if (t0.has_on_bounded_side(tet[i]) || tet.has_on_bounded_side(t0[i]))
1139 return true;
1140 }
1141 }
1142
1143 tets.push_back(tet);
1144 }
1145
1146 return false;
1147 }
1148
1149 //-----------------------------------------------------------------------------
1150
1151}
1152#endif
1153
1154#endif
Definition: Point.h:41
Definition: adapt.h:30
void dolfin_error(std::string location, std::string task, std::string reason,...)
Definition: log.cpp:129
double norm(const GenericVector &x, std::string norm_type="l2")
Definition: solve.cpp:173
void info(std::string msg,...)
Print message.
Definition: log.cpp:72