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
55 namespace 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 
146 namespace
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 
382 namespace 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
dolfin
Definition: adapt.h:29
dolfin::Point
Definition: Point.h:40
dolfin::info
void info(std::string msg,...)
Print message.
Definition: log.cpp:72
dolfin::dolfin_error
void dolfin_error(std::string location, std::string task, std::string reason,...)
Definition: log.cpp:129
Generated by   doxygen 1.8.11