FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
cgal.cpp
1// FEAT3: Finite Element Analysis Toolbox, Version 3
2// Copyright (C) 2010 by Stefan Turek & the FEAT group
3// FEAT3 is released under the GNU General Public License version 3,
4// see the file 'copyright.txt' in the top level directory for details.
5
6#include <kernel/base_header.hpp>
7#include <kernel/util/assertion.hpp>
8#include <iostream>
9#include <fstream>
10#include <algorithm>
11
12#ifdef FEAT_HAVE_CGAL
13
14#define CGAL_HEADER_ONLY 1
15
16#ifdef FEAT_COMPILER_MICROSOFT
17#pragma warning(disable: 4244)
18#endif // FEAT_COMPILER_MICROSOFT
19
20FEAT_DISABLE_WARNINGS
21#define CGAL_NO_DEPRECATION_WARNINGS 1
22#include <CGAL/Simple_cartesian.h>
23#include <CGAL/AABB_tree.h>
24#include <CGAL/AABB_traits.h> // deprecated in CGAL 6.x
25#include <CGAL/Surface_mesh.h>
26#include <CGAL/AABB_face_graph_triangle_primitive.h>
27#include <CGAL/algorithm.h>
28#include <CGAL/Side_of_triangle_mesh.h>
29#include <CGAL/Aff_transformation_3.h>
30#include <CGAL/Polygon_mesh_processing/transform.h>
31#include <CGAL/Polygon_mesh_processing/compute_normal.h>
32#include <CGAL/Polygon_mesh_processing/detect_features.h>
33#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
34FEAT_RESTORE_WARNINGS
35
36// check if cgal threading support is enabled
37#ifdef FEAT_HAVE_OMP
38#ifndef CGAL_HAS_THREADS
39static_assert(false, "No cgal multithreading support");
40#endif
41
42#ifndef BOOST_HAS_THREADS
43static_assert(false, "Boost has no threads");
44#endif
45#endif
46
47
48#ifdef FEAT_COMPILER_INTEL
49#pragma warning disable 3280
50#pragma warning disable 1224
51#endif //FEAT_COMPILER_INTEL
52
53#ifdef FEAT_COMPILER_MICROSOFT
54#pragma warning(default: 4244)
55#endif // FEAT_COMPILER_MICROSOFT
56
57#undef FEAT_EICKT
58
59#include <kernel/geometry/cgal.hpp>
60
61#ifdef FEAT_HAVE_HALFMATH
62#include <kernel/util/half.hpp>
63#endif
64
65namespace FEAT::Geometry
66{
67 template<typename DT_>
68 struct CGALTypeWrapper
69 {
70 typedef DT_ CGALDT_;
71 typedef typename CGAL::Simple_cartesian<DT_> K;
72 typedef typename K::Point_3 Point_;
73 typedef typename K::Vector_3 Vector_;
74 typedef typename K::Segment_3 Segment_;
75 typedef typename CGAL::Surface_mesh<Point_> Polyhedron_;
76 typedef typename CGAL::AABB_face_graph_triangle_primitive<Polyhedron_> Primitive_;
77 typedef typename CGAL::AABB_traits<K, Primitive_> Traits_;
78 typedef typename CGAL::AABB_tree<Traits_> Tree_;
79 typedef typename CGAL::Side_of_triangle_mesh<Polyhedron_, K> Point_inside_;
80 typedef typename CGAL::Aff_transformation_3<K> Transformation_;
81 typedef typename CGAL::AABB_traits<K, Primitive_>::Point_and_primitive_id Point_and_primitive_id_;
82 typedef std::optional< typename Tree_::template Intersection_and_primitive_id<Segment_>::Type > Segment_intersection_;
83 };
84
85 template<typename DT_>
86 struct CGALWrapperData
87 {
88 typedef CGALTypeWrapper<DT_> TW_;
89 typename TW_::Polyhedron_ * _polyhedron;
90 typename TW_::Tree_ * _tree;
91 typename TW_::Point_inside_ * _inside_tester;
92
93 CGALWrapperData() :
94 _polyhedron(nullptr),
95 _tree(nullptr),
96 _inside_tester(nullptr)
97 {
98 }
99
100 CGALWrapperData(const CGALWrapperData&) = delete;
101 CGALWrapperData& operator=(const CGALWrapperData&) = delete;
102
103 void swap(CGALWrapperData& other)
104 {
105 std::swap(this->_polyhedron, other._polyhedron);
106 std::swap(this->_tree, other._tree);
107 std::swap(this->_inside_tester, other._inside_tester);
108 }
109
110 CGALWrapperData(CGALWrapperData&& other) noexcept
111 : _polyhedron(nullptr), _tree(nullptr), _inside_tester(nullptr)
112 {
113 this->swap(other);
114 }
115
116 CGALWrapperData& operator=(CGALWrapperData&& other) noexcept
117 {
118 this->swap(other);
119 return *this;
120 }
121
122
123 ~CGALWrapperData()
124 {
125 if(_inside_tester)
126 delete _inside_tester;
127 if(_tree)
128 delete _tree;
129 if(_polyhedron)
130 delete _polyhedron;
131 }
132 };
133
134 #ifdef FEAT_HAVE_QUADMATH
135 template<>
136 struct CGALTypeWrapper<__float128>
137 {
138 typedef double CGALDT_;
139 typedef typename CGAL::Simple_cartesian<double> K;
140 typedef typename K::Point_3 Point_;
141 typedef typename K::Vector_3 Vector_;
142 typedef typename K::Segment_3 Segment_;
143 typedef typename CGAL::Surface_mesh<Point_> Polyhedron_;
144 typedef typename CGAL::AABB_face_graph_triangle_primitive<Polyhedron_> Primitive_;
145 typedef typename CGAL::AABB_traits<K, Primitive_> Traits_;
146 typedef typename CGAL::AABB_tree<Traits_> Tree_;
147 typedef typename CGAL::Side_of_triangle_mesh<Polyhedron_, K> Point_inside_;
148 typedef typename CGAL::Aff_transformation_3<K> Transformation_;
149 typedef typename CGAL::AABB_traits<K, Primitive_>::Point_and_primitive_id Point_and_primitive_id_;
150 typedef std::optional< typename Tree_::template Intersection_and_primitive_id<Segment_>::Type > Segment_intersection_;
151 };
152 #endif
153
154 template<typename DT_>
155 using PointTypeAlias = typename FEAT::Geometry::template CGALWrapper<DT_>::PointType;
156
157 template<typename DT_>
158 using TransformMatrixAlias = typename FEAT::Geometry::template CGALWrapper<DT_>::TransformMatrix;
159
160 template<typename DT_>
161 CGALWrapper<DT_>::CGALWrapper(std::istream & file, CGALFileMode file_mode) :
162 _cgal_data(nullptr),
163 _expensive_intersection(false)
164 {
165 _parse_mesh(file, file_mode);
166 }
167
168 template<typename DT_>
169 CGALWrapper<DT_>::CGALWrapper(const String & filename, CGALFileMode file_mode) :
170 _cgal_data(nullptr),
171 _expensive_intersection(false)
172 {
173 std::ifstream file(filename.c_str());
174 XASSERTM(file.is_open(), "CGALWrapper: file read error: " + filename + " !");
175 _parse_mesh(file, file_mode);
176 file.close();
177 }
178
179 template<typename DT_>
180 void CGALWrapper<DT_>::set_expensive_intersection(bool expensive_test)
181 {
182 this->_expensive_intersection = expensive_test;
183 }
184
185 template <typename DT_>
186 static CGALWrapperData<DT_> *from_topology(
187 const std::vector<typename CGALWrapper<DT_>::PointType> &vertices,
188 const std::vector<std::array<Index, 3>> &faces)
189 {
190 using VertexType = typename CGALTypeWrapper<DT_>::Point_;
191 using SurfaceMesh = typename CGALTypeWrapper<DT_>::Polyhedron_;
192 using VertexIndex = typename SurfaceMesh::Vertex_index;
193 using FaceIndex = typename SurfaceMesh::Face_index;
194
195 auto* result = new CGALWrapperData<DT_>;
196 result->_polyhedron = new SurfaceMesh;
197
198 // Add vertices to surface mesh
199 std::vector<VertexIndex> vs(vertices.size());
200 for(Index i(0); i < vertices.size(); i++)
201 {
202 VertexType vertex(vertices[i][0], vertices[i][1], vertices[i][2]);
203 vs[i] = result->_polyhedron->add_vertex(vertex);
204 }
205
206 // Add faces to surface mesh
207 for(const std::array<Index, 3>& face : faces)
208 {
209 FaceIndex i = result->_polyhedron->add_face(vs[face[0]], vs[face[1]], vs[face[2]]);
210
211 XASSERTM(i != SurfaceMesh::null_face(), "Failed to add face to surface mesh!");
212 }
213
214 return result;
215 }
216
217 template<typename DT_>
218 CGALWrapper<DT_>::CGALWrapper(
219 const std::vector<typename CGALWrapper<DT_>::PointType> &vertices,
220 const std::vector<std::array<Index, 3>> &faces) :
221 _cgal_data(from_topology<DT_>(vertices, faces)),
222 _expensive_intersection(false)
223 {
224 _init_wrapper();
225 }
226
227 template<typename DT_>
228 CGALWrapper<DT_>::CGALWrapper::~CGALWrapper()
229 {
230 if(_cgal_data)
231 delete (CGALWrapperData<DT_>*)_cgal_data;
232 }
233
234
235 template<typename DT_>
236 bool CGALWrapper<DT_>::point_inside(DT_ x, DT_ y, DT_ z) const
237 {
238 typedef typename CGALTypeWrapper<DT_>::CGALDT_ IDT_;
239 typename CGALTypeWrapper<DT_>::Point_ query{IDT_(x), IDT_(y), IDT_(z)};
240
241 // Determine the side and return true if inside!
242 CGALWrapperData<DT_> * cd = (CGALWrapperData<DT_>*)_cgal_data;
243 return (*(cd->_inside_tester))(query) == CGAL::ON_BOUNDED_SIDE;
244 }
245
246 template<typename DT_>
247 bool CGALWrapper<DT_>::point_not_outside(DT_ x, DT_ y, DT_ z) const
248 {
249 typedef typename CGALTypeWrapper<DT_>::CGALDT_ IDT_;
250 typename CGALTypeWrapper<DT_>::Point_ query{IDT_(x), IDT_(y), IDT_(z)};
251
252 // Determine the side and return true if inside!
253 CGALWrapperData<DT_> * cd = (CGALWrapperData<DT_>*)_cgal_data;
254 CGAL::Bounded_side test_result = (*(cd->_inside_tester))(query);
255 return test_result != CGAL::ON_UNBOUNDED_SIDE;
256 }
257
258 template<typename DT_>
259 DT_ CGALWrapper<DT_>::squared_distance(DT_ x, DT_ y, DT_ z) const
260 {
261 typedef typename CGALTypeWrapper<DT_>::CGALDT_ IDT_;
262 typename CGALTypeWrapper<DT_>::Point_ query{IDT_(x), IDT_(y), IDT_(z)};
263
264 CGALWrapperData<DT_> * cd = (CGALWrapperData<DT_>*)_cgal_data;
265 return DT_(cd->_tree->squared_distance(query));
266 }
267
268 template<typename DT_>
269 DT_ CGALWrapper<DT_>::squared_distance_to_feature(const CGALFeature& f, const PointType& p) const
270 {
271 return (p - closest_point_on_feature(f, p)).norm_euclid_sqr();
272 }
273
274
275 template<typename DT_>
276 typename CGALWrapper<DT_>::PointType CGALWrapper<DT_>::closest_point(const PointType& point) const
277 {
278 return closest_point(point[0], point[1], point[2]);
279 }
280
281 template<typename DT_>
282 typename CGALWrapper<DT_>::PointType CGALWrapper<DT_>::closest_point(DT_ x, DT_ y, DT_ z) const
283 {
284 typedef typename CGALTypeWrapper<DT_>::CGALDT_ IDT_;
285 typename CGALTypeWrapper<DT_>::Point_ query{IDT_(x), IDT_(y), IDT_(z)};
286 CGALWrapperData<DT_> * cd = (CGALWrapperData<DT_>*)_cgal_data;
287 typename CGALTypeWrapper<DT_>::Point_ n_point = cd->_tree->closest_point(query);
288 return PointType{{DT_(n_point[0]), DT_(n_point[1]), DT_(n_point[2])}};
289 }
290
291 template<typename DT_>
292 typename CGALWrapper<DT_>::PointType CGALWrapper<DT_>::closest_point(const PointType& point, PointType& primitive_grad) const
293 {
294 typedef typename CGALTypeWrapper<DT_>::CGALDT_ IDT_;
295 typename CGALTypeWrapper<DT_>::Point_ query{IDT_(point[0]), IDT_(point[1]), IDT_(point[2])};
296 CGALWrapperData<DT_> * cd = (CGALWrapperData<DT_>*)_cgal_data;
297 const typename CGALTypeWrapper<DT_>::Point_and_primitive_id_& cl_p_query = cd->_tree->closest_point_and_primitive(query);
298 typename CGALTypeWrapper<DT_>::Polyhedron_::Face_index f = cl_p_query.second;
299 auto v = CGAL::Polygon_mesh_processing::compute_face_normal(f,*(cd->_polyhedron));
300 for(int i = 0; i < 3; ++i)
301 primitive_grad[i] = DT_(v[i]);
302 return PointType{{DT_(cl_p_query.first[0]), DT_(cl_p_query.first[1]), DT_(cl_p_query.first[2])}};
303 }
304
305 template<typename DT_>
306 typename CGALWrapper<DT_>::PointType
307 CGALWrapper<DT_>::closest_point_on_feature(const CGALFeature& f, const PointType& query) const
308 {
309 using PolyhedronType = typename CGALTypeWrapper<DT_>::Polyhedron_;
310 using VertexIndex = typename PolyhedronType::Vertex_index;
311
312 CGALWrapperData<DT_>* cd = (CGALWrapperData<DT_>*)_cgal_data;
313 PolyhedronType* poly = cd->_polyhedron;
314
315 const auto to_feat_point = [](const auto& p)
316 {
317 return PointType{p.x(), p.y(), p.z()};
318 };
319
320 PointType closest = to_feat_point(poly->point(VertexIndex(f[0])));
321 DT_ distance_to_closest = (query - closest).norm_euclid_sqr();
322
323 const auto closest_point_on_segment = [&](PointType start, PointType end, PointType q)
324 {
325 PointType segment = end - start;
326 PointType segment_dir = segment.normalize();
327
328 // Orthogonal projection of (query - start) onto the segment
329 DT_ coefficient = Math::clamp(dot(segment_dir, q - start), DT_(0.0), DT_(1.0));
330
331 return start + coefficient * segment_dir;
332 };
333
334 for(std::size_t i(0); i < f.size() - 1; i++)
335 {
336 PointType segment_start = to_feat_point(poly->point(VertexIndex(f[i])));
337 PointType segment_end = to_feat_point(poly->point(VertexIndex(f[i + 1])));
338
339 PointType closest_on_segment = closest_point_on_segment(segment_start, segment_end, query);
340
341 DT_ distance_to_segment = (query - closest_on_segment).norm_euclid_sqr();
342 if(distance_to_segment < distance_to_closest)
343 {
344 closest = closest_on_segment;
345 distance_to_closest = distance_to_segment;
346 }
347 }
348
349 return closest;
350 }
351
352 template<typename DT_>
353 typename CGALWrapper<DT_>::PointType CGALWrapper<DT_>::point(std::uint32_t idx) const
354 {
355 using PolyhedronType = typename CGALTypeWrapper<DT_>::Polyhedron_;
356 using VertexIndex = typename PolyhedronType::Vertex_index;
357
358 const auto to_feat_point = [](const auto& p)
359 {
360 return PointType{p.x(), p.y(), p.z()};
361 };
362
363 CGALWrapperData<DT_>* cd = (CGALWrapperData<DT_>*)_cgal_data;
364
365 return to_feat_point(cd->_polyhedron->point(VertexIndex(idx)));
366 }
367
368 template<typename DT_>
369 std::vector<typename CGALWrapper<DT_>::PointType> CGALWrapper<DT_>::points() const
370 {
371 std::vector<PointType> result(get_num_entities(0));
372 for(std::uint32_t i(0); i < result.size(); i++)
373 {
374 result[i] = point(i);
375 }
376 return result;
377 }
378
379 template<typename DT_>
380 void CGALWrapper<DT_>::_parse_mesh(std::istream & file, CGALFileMode file_mode)
381 {
382 delete (CGALWrapperData<DT_>*)_cgal_data;
383
384 _cgal_data = new CGALWrapperData<DT_>;
385 CGALWrapperData<DT_> * cd = (CGALWrapperData<DT_>*)_cgal_data;
386
387 cd->_polyhedron = new typename CGALTypeWrapper<DT_>::Polyhedron_;
388
389 bool status(false);
390 switch(file_mode)
391 {
392 case CGALFileMode::fm_off:
393 status = CGAL::IO::read_OFF(file, *(cd->_polyhedron));
394 XASSERTM(status == true, "CGAL::IO::read_OFF: read/parse error !");
395 break;
396 case CGALFileMode::fm_obj:
397 status = CGAL::IO::read_OBJ(file, *(cd->_polyhedron));
398 XASSERTM(status == true, "CGAL::IO::read_OBJ: read/parse error !");
399 break;
400 default:
401 XASSERTM(false, "CGAL FileMode not supported!");
402 }
403
404 _init_wrapper();
405
406 }
407
408 template<typename DT_>
409 CGALFeatureNetwork CGALWrapper<DT_>::detect_features(DT_ critical_angle)
410 {
411 using CGALDT = typename CGALTypeWrapper<DT_>::CGALDT_;
412 using PolyhedronType = typename CGALTypeWrapper<DT_>::Polyhedron_;
413 using VertexIndex = typename PolyhedronType::Vertex_index;
414 using EdgeIndex = typename PolyhedronType::Edge_index;
415 using HalfedgeIndex = typename PolyhedronType::Halfedge_index;
416 using EdgeIsSharpMapType = typename PolyhedronType::template Property_map<EdgeIndex, bool>;
417
418 CGALFeatureNetwork result;
419
420 auto* cd = static_cast<CGALWrapperData<DT_>*>(_cgal_data);
421 PolyhedronType* poly = cd->_polyhedron;
422
424 // Detect sharp edges
426
427 std::pair<EdgeIsSharpMapType, bool> map = poly->template add_property_map<EdgeIndex, bool>("e:is_sharp");
428
429 if(map.second)
430 {
431 // Property was newly created. Fill it with proper data.
432 CGAL::Polygon_mesh_processing::detect_sharp_edges(*poly, CGALDT(critical_angle), map.first);
433 }
434
436 // Find feature starting points
438
439 std::vector<VertexIndex> seeds;
440
441 for(VertexIndex v : poly->vertices())
442 {
443 Index incident_feature_edges(0);
444
445 for(HalfedgeIndex he : poly->halfedges_around_target(poly->halfedge(v)))
446 {
447 EdgeIndex edge = poly->edge(he);
448 if(map.first[edge])
449 {
450 incident_feature_edges++;
451 }
452 }
453
454 // At least one incident feature that is not just a feature passing through this vertex
455 if(incident_feature_edges != 0 && incident_feature_edges != 2)
456 {
457 seeds.push_back(v);
458 }
459 }
460
462 // Flood fill for feature lines
464
465 std::vector<bool> visited(poly->number_of_edges(), false);
466 const auto flood_feature = [&](VertexIndex v, HalfedgeIndex h)
467 {
468 CGALFeature feature;
469 feature.push_back(static_cast<std::uint32_t>(v));
470 VertexIndex current = poly->source(h);
471 feature.push_back(static_cast<std::uint32_t>(current));
472
473 visited[poly->edge(h)] = true;
474
475 while(current != v && std::find(seeds.begin(), seeds.end(), current) == seeds.end())
476 {
477 for(HalfedgeIndex he : poly->halfedges_around_target(poly->halfedge(current)))
478 {
479 EdgeIndex e = poly->edge(he);
480 if(!map.first[e] || visited[e])
481 {
482 // Edge is either not sharp or has already been visited.
483 // Either way it is not the way forward
484 continue;
485 }
486
487 current = poly->source(he);
488 feature.push_back(static_cast<std::uint32_t>(current));
489 visited[e] = true;
490 }
491 }
492
493 return feature;
494 };
495
497 // Fill in features
499
500 for(VertexIndex seed : seeds)
501 {
502 for(HalfedgeIndex he : poly->halfedges_around_target(poly->halfedge(seed)))
503 {
504 EdgeIndex e = poly->edge(he);
505 if(map.first[e] && !visited[e])
506 {
507 // Edge is sharp and unvisited. It must be the start of a new feature.
508 result.push_back(flood_feature(seed, he));
509 }
510 }
511 }
512
513 // On isolated circular features, all vertices have exactly two incident sharp edges.
514 // They thus contain to seed vertices and have to be handled specially.
515
516 for(EdgeIndex e : poly->edges())
517 {
518 if(map.first[e] && !visited[e])
519 {
520 HalfedgeIndex he = poly->halfedge(e);
521 VertexIndex start = poly->target(he);
522
523 result.push_back(flood_feature(start, he));
524 }
525 }
526
527 return result;
528 }
529
530 template<typename DT_>
531 void CGALWrapper<DT_>::transform(const TransformMatrix& scale_rot, const PointType& translation, DT_ scale)
532 {
533 typedef typename CGALTypeWrapper<DT_>::CGALDT_ IDT_;
534 CGALWrapperData<DT_> * cd = (CGALWrapperData<DT_>*)_cgal_data;
535 //delete tree and _inside_tester
536 _delete_tree();
537 //build affine transformation
538 typename CGALTypeWrapper<DT_>::Transformation_ trafo_mat{(IDT_)(scale_rot[0][0]), (IDT_)(scale_rot[0][1]), (IDT_)(scale_rot[0][2]), (IDT_)(translation[0]),
539 (IDT_)(scale_rot[1][0]), (IDT_)(scale_rot[1][1]), (IDT_)(scale_rot[1][2]), (IDT_)(translation[1]),
540 (IDT_)(scale_rot[2][0]), (IDT_)(scale_rot[2][1]), (IDT_)(scale_rot[2][2]), (IDT_)(translation[2]),
541 IDT_(scale)};
542 //apply affine transformation on polyhedron
543 CGAL::Polygon_mesh_processing::transform(trafo_mat, *(cd->_polyhedron));
544 //create new trees
545 _init_wrapper();
546
547
548 }
549
550 template<typename DT_>
551 void CGALWrapper<DT_>::displace_vertices(const std::vector<PointType>& offsets)
552 {
553 using Point = typename CGALTypeWrapper<DT_>::Point_;
554 using Vector = typename CGALTypeWrapper<DT_>::Vector_;
555
556 ASSERT(get_num_vertices() == offsets.size());
557
558 auto* cd = static_cast<CGALWrapperData<DT_>*>(_cgal_data);
559
560 _delete_tree();
561 std::size_t i(0);
562 for(auto vertex_idx : cd->_polyhedron->vertices())
563 {
564 Point& vertex = cd->_polyhedron->point(vertex_idx);
565 const auto& offset = offsets[i++];
566 vertex += Vector(offset[0], offset[1], offset[2]);
567 }
568 _init_wrapper();
569 }
570
571 template<typename DT_>
572 Index CGALWrapper<DT_>::get_num_vertices() const
573 {
574 auto* cd = static_cast<CGALWrapperData<DT_>*>(_cgal_data);
575 return cd->_polyhedron->number_of_vertices();
576 }
577
578 template<typename DT_>
579 Index CGALWrapper<DT_>::get_num_entities(int dim) const
580 {
581 ASSERT(0<= dim && dim <= 2);
582
583 auto* cd = static_cast<CGALWrapperData<DT_>*>(_cgal_data);
584 switch(dim)
585 {
586 case 0: return cd->_polyhedron->number_of_vertices();
587 case 1: return cd->_polyhedron->number_of_edges();
588 case 2: return cd->_polyhedron->number_of_faces();
589 default: XABORTM("Unsupported dimension!");
590 }
591
592 return 0;
593 }
594
595 template<typename DT_>
596 std::vector<typename CGALWrapper<DT_>::PointType> CGALWrapper<DT_>::outer_normals_at_faces() const
597 {
598 using PolyhedronType = typename CGALTypeWrapper<DT_>::Polyhedron_;
599
600 CGALWrapperData<DT_>* cd = (CGALWrapperData<DT_>*)_cgal_data;
601 PolyhedronType* poly = cd->_polyhedron;
602
603 std::vector<PointType> normals(this->get_num_entities(2));
604 for(const auto facet : poly->faces())
605 {
606 auto outer_normal = CGAL::Polygon_mesh_processing::compute_face_normal(facet, *poly);
607 normals.at(facet) = PointType{outer_normal.x(), outer_normal.y(), outer_normal.z()};
608 }
609
610 return normals;
611 }
612
613 template<typename DT_>
614 Index CGALVerticesAroundFaceAdjactor<DT_>::get_num_nodes_domain() const
615 {
616 return Index(_wrapper->get_num_entities(2));
617 }
618
619 template<typename DT_>
620 Index CGALVerticesAroundFaceAdjactor<DT_>::get_num_nodes_image() const
621 {
622 return Index(_wrapper->get_num_entities(0));
623 }
624
625 template<typename DT_>
626 CGALValueIteratorWrapper<Index> CGALVerticesAroundFaceAdjactor<DT_>::image_begin(Index face) const
627 {
628 using PolyhedronType = typename CGALTypeWrapper<DT_>::Polyhedron_;
629 using FaceIndex = typename PolyhedronType::Face_index;
630
631 ASSERT(face < _wrapper->get_num_entities(2));
632
633 auto* cd = static_cast<CGALWrapperData<DT_>*>(_wrapper->_cgal_data);
634
635 auto range = cd->_polyhedron->vertices_around_face(cd->_polyhedron->halfedge(FaceIndex(typename PolyhedronType::size_type(face))));
636 return CGALValueIteratorWrapper<Index>([it = range.begin(), end = range.end()]() mutable {
637 if(it != end)
638 {
639 auto val = *it;
640 it++;
641 return std::optional<Index>(val);
642 }
643 else
644 {
645 return std::optional<Index>{};
646 }
647 });
648 }
649
650 template<typename DT_>
651 CGALValueIteratorWrapper<Index> CGALVerticesAroundFaceAdjactor<DT_>::image_end(Index /*face*/) const
652 {
653 return CGALValueIteratorWrapper<Index>([]() { return std::nullopt; });
654 }
655
656 template<typename DT_>
657 CGALVerticesAroundFaceAdjactor<DT_> CGALWrapper<DT_>::vertices_around_face() const
658 {
659 return CGALVerticesAroundFaceAdjactor<DT_>(this);
660 }
661
662 template<typename DT_>
663 std::size_t CGALWrapper<DT_>::bytes() const
664 {
665 auto* cd = static_cast<CGALWrapperData<DT_>*>(_cgal_data);
666
667 // CGALs surface mesh uses std::uint32_t to store it indices
668 constexpr std::size_t index_size = sizeof(std::uint32_t);
669
670 // CGAL stores an halfedge-index, the vertex coordinates, and a removed flag for each vertex
671 constexpr std::size_t bytes_per_vertex = index_size + sizeof(typename CGALTypeWrapper<DT_>::Point_) + sizeof(bool);
672
673 // CGAL stores 4 indices per halfedge
674 constexpr std::size_t bytes_per_halfedge = 4 * index_size;
675
676 // CGAL stores just the removed flag for edges
677 constexpr std::size_t bytes_per_edge = sizeof(bool);
678
679 // CGAL stores 1 index per face and the removed flag
680 constexpr std::size_t bytes_per_face = index_size + sizeof(bool);
681
682 std::size_t result = 0;
683
684 if(cd->_polyhedron)
685 {
686 result += cd->_polyhedron->number_of_vertices() * bytes_per_vertex +
687 cd->_polyhedron->number_of_halfedges() * bytes_per_halfedge +
688 cd->_polyhedron->number_of_edges() * bytes_per_edge +
689 cd->_polyhedron->number_of_faces() * bytes_per_face;
690 }
691
692 if(cd->_tree != nullptr)
693 {
694 // See https://doc.cgal.org/latest/AABB_tree/index.html, section 4.2 Memory
695 constexpr std::size_t bytes_per_primitive = 150;
696
697 result += cd->_tree->size() * bytes_per_primitive;
698 }
699
700 // Inside tester stores no meaningful own data
701
702 return result;
703 }
704
705 template<typename DT_>
706 void CGALWrapper<DT_>::write_off(std::ostream& stream) const
707 {
708 CGALWrapperData<DT_> * cd = (CGALWrapperData<DT_>*)_cgal_data;
709
710 bool status(false);
711 status = CGAL::IO::write_OFF(stream, *(cd->_polyhedron));
712 if(!status)
713 throw FileError("Failed to write off file to stream");
714 }
715
716 template<typename DT_>
717 void CGALWrapper<DT_>::_delete_tree()
718 {
719 CGALWrapperData<DT_> * cd = (CGALWrapperData<DT_>*)_cgal_data;
720 XASSERTM(cd->_polyhedron != nullptr, "ERROR: Polyhedron is not initialized!");
721 delete cd->_inside_tester;
722 cd->_inside_tester = nullptr;
723 delete cd->_tree;
724 cd->_tree = nullptr;
725 }
726
727 template<typename DT_>
728 void CGALWrapper<DT_>::_init_wrapper()
729 {
730 CGALWrapperData<DT_> * cd = (CGALWrapperData<DT_>*)_cgal_data;
731 XASSERTM(cd->_polyhedron != nullptr, "ERROR: Polyhedron is not initialized!");
732 XASSERTM(cd->_tree == nullptr && cd->_inside_tester == nullptr, "ERROR: Tree or Inside Tester are already initialized");
733
734 // Construct AABB tree with a KdTree
735 cd->_tree = new typename CGALTypeWrapper<DT_>::Tree_(faces(*(cd->_polyhedron)).first, faces(*(cd->_polyhedron)).second, *(cd->_polyhedron));
736 cd->_tree->accelerate_distance_queries();
737 // Initialize the point-in-polyhedron tester
738 cd->_inside_tester = new typename CGALTypeWrapper<DT_>::Point_inside_(*(cd->_tree));
739 }
740
741 template<typename DT_>
742 CGALWrapper<DT_>::CGALWrapper(CGALWrapper<DT_>&& other) noexcept
743 : _cgal_data(nullptr),
744 _expensive_intersection(other._expensive_intersection)
745 {
746 std::swap(this->_cgal_data, other._cgal_data);
747 }
748
749 template<typename DT_>
750 CGALWrapper<DT_>& CGALWrapper<DT_>::operator=(CGALWrapper<DT_>&& other) noexcept
751 {
752 this->_expensive_intersection = other._expensive_intersection;
753 std::swap(this->_cgal_data, other._cgal_data);
754 return *this;
755 }
756
757 template<typename DT_>
758 bool FEAT::Geometry::CGALWrapper<DT_>::intersects_line(const FEAT::Geometry::CGALWrapper<DT_>::PointType& a, const FEAT::Geometry::CGALWrapper<DT_>::PointType& b) const
759 {
760 typename CGALTypeWrapper<DT_>::Point_ t(a[0], a[1], a[2]), z(b[0], b[1], b[2]);
761 typename CGALTypeWrapper<DT_>::Segment_ seg(t,z);
762 return static_cast<CGALWrapperData<DT_>*>(this->_cgal_data)->_tree->do_intersect(seg);
763 }
764
765 namespace Intern
766 {
767 template<typename DT_>
768 static inline void split_hexaeder_surface(typename CGALTypeWrapper<DT_>::Polyhedron_& surface, const std::array<typename FEAT::Geometry::CGALWrapper<DT_>::PointType, 8>& points)
769 {
770 typename CGALTypeWrapper<DT_>::Polyhedron_::Vertex_index vtx_desc[8];
771 for(std::size_t k = 0; k < std::size_t (8); ++k)
772 {
773 vtx_desc[k] = surface.add_vertex(typename FEAT::Geometry::CGALTypeWrapper<DT_>::Point_(points[k][0], points[k][1], points[k][2]));
774 }
775
776 constexpr std::size_t faces[12][3] =
777 {
778 {0, 1, 2}, {1, 3, 2}, //bottom
779 {6, 5, 4}, {7, 5, 6}, //top
780 {0, 2, 4}, {2, 6, 4}, //left
781 {3, 1, 5}, {3, 5, 7}, //right
782 {4, 1, 0}, {5, 1, 4}, //front
783 {2, 3, 6}, {3, 7, 6} //back
784 };
785
786 for(std::size_t f = 0; f < std::size_t(12); ++f)
787 {
788 surface.add_face(vtx_desc[faces[f][0]], vtx_desc[faces[f][1]], vtx_desc[faces[f][2]]);
789 }
790 }
791 }
792
793 template<typename DT_>
794 bool FEAT::Geometry::CGALWrapper<DT_>::intersects_hexaedron(const std::array<FEAT::Geometry::CGALWrapper<DT_>::PointType, 8>& points) const
795 {
796 typename CGALTypeWrapper<DT_>::Polyhedron_ hexaeder;
797 typedef typename CGALTypeWrapper<DT_>::Point_ Point;
798 Intern::template split_hexaeder_surface<DT_>(hexaeder, points);
799 const auto* tree = static_cast<CGALWrapperData<DT_>*>(this->_cgal_data)->_tree;
800
801 // first of all, check if the first surface point is contained in our mesh cell, if yes, we are done
802 // if not we know, that the surface mesh is not wholly contained in our mesh cell
803 {
804 const auto* poly = static_cast<CGALWrapperData<DT_>*>(this->_cgal_data)->_polyhedron;
805 typename CGALTypeWrapper<DT_>::Point_inside_ inside(hexaeder);
806 auto res = inside(poly->point(*std::begin(poly->vertices())));
807 if(res == CGAL::ON_BOUNDED_SIDE || res == CGAL::ON_BOUNDARY)
808 {
809 return true;
810 }
811 }
812
813 // first, create boundingbox
814 auto boxA = hexaeder.point(*std::begin(hexaeder.vertices())).bbox();
815 for(auto v : hexaeder.vertices())
816 {
817 boxA += hexaeder.point(v).bbox();
818 }
819
820 // check if we overlap with the surface bb
821 {
822 auto boxB = tree->bbox();
823 if(!CGAL::do_overlap(boxA, boxB))
824 {
825 return false;
826 }
827 }
828
829 // first a few easy tests, check if any vertex of our mesh (and the midpoint) is inside the surface mesh
830 {
831 typename FEAT::Geometry::CGALWrapper<DT_>::PointType mid_point(DT_(0));
832 for(std::size_t k = 0; k < 8u; ++k)
833 {
834 mid_point += points[k];
835 }
836 mid_point *= DT_(1)/DT_(8);
837
838 auto& inside_tester = *static_cast<CGALWrapperData<DT_>*>(this->_cgal_data)->_inside_tester;
839 if(inside_tester(Point(mid_point[0], mid_point[1], mid_point[2])) != CGAL::ON_UNBOUNDED_SIDE)
840 return true;
841
842 for(auto v : hexaeder.vertices())
843 {
844 auto p = hexaeder.point(v);
845 if(inside_tester(p) != CGAL::ON_UNBOUNDED_SIDE)
846 return true;
847 }
848 }
849
850 if(!this->_expensive_intersection)
851 {
852 return tree->do_intersect(boxA);
853 }
854 else
855 {
856 // we construct a list of all primitives (i.e. surface elements) intersected with our boundinbox
857 std::list<typename FEAT::Geometry::CGALTypeWrapper<DT_>::Tree_::Primitive_id> primitives;
858
859 // gather all intersections
860 tree->all_intersected_primitives(boxA, std::back_inserter(primitives));
861 const auto* poly = static_cast<CGALWrapperData<DT_>*>(this->_cgal_data)->_polyhedron;
862
863 // iterate over all primitives (i.e. faces) that possibly intersect our boundingbox
864 for(const auto& sf: primitives)
865 {
866 std::array<Point, 3> triang;
867 {
868 std::size_t k = 0;
869 for(auto v : poly->vertices_around_face(poly->halfedge(sf)))
870 {
871 triang[k++] = poly->point(v);
872 }
873 }
874 CGAL::Triangle_3<typename CGALTypeWrapper<DT_>::K> triangle_s{triang[0], triang[1], triang[2]};
875
876 // now test if we actually intersect, testing only intersection is enough, since we know the geometries do not
877 // contain each other
878 for(auto f : hexaeder.faces())
879 {
880 {
881 std::size_t k = 0;
882 for(auto v : hexaeder.vertices_around_face(hexaeder.halfedge(f)))
883 {
884 triang[k++] = hexaeder.point(v);
885 }
886 }
887 CGAL::Triangle_3<typename CGALTypeWrapper<DT_>::K> triangle_a(triang[0], triang[1], triang[2]);
888
889 if(CGAL::do_intersect(triangle_s, triangle_a))
890 return true;
891 }
892 }
893 }
894
895
896 // no intersection
897 return false;
898 }
899
900} // namespace Geometry
901
902
903// explicitly instantiate templates for all sensible datatype
904
905template class FEAT::Geometry::CGALWrapper<double>;
906template class FEAT::Geometry::CGALVerticesAroundFaceAdjactor<double>;
907
908template class FEAT::Geometry::CGALWrapper<float>;
909template class FEAT::Geometry::CGALVerticesAroundFaceAdjactor<float>;
910
911#ifdef FEAT_HAVE_HALFMATH
912// TODO: Figure out missing operators and ambigous overloads to support this properly
913//template class FEAT::Geometry::CGALWrapper<FEAT::Half>;
914//template class FEAT::Geometry::CGALVerticesAroundFaceAdjactor<FEAT::Half>;
915#endif
916
917#ifdef FEAT_HAVE_QUADMATH
918template class FEAT::Geometry::CGALWrapper<__float128>;
919template class FEAT::Geometry::CGALVerticesAroundFaceAdjactor<__float128>;
920#endif
921
922#elif !defined(DOXYGEN)
923
924// Dummy class instance to silence ipo linker optimization warnings about empty libgeometry
925class ipo_foobar_geometry_cgal
926{
927public:
928 int i;
929 ipo_foobar_geometry_cgal() :
930 i(0)
931 {
932 (void)i;
933 }
934} ipo_barfoo_geometry_cgal;
935
936#endif // FEAT_HAVE_CGAL
ASSERT
#define ASSERT(expr)
Debug-Assertion macro definition.
Definition: assertion.hpp:229
XABORTM
#define XABORTM(msg)
Abortion macro definition with custom message.
Definition: assertion.hpp:192
XASSERTM
#define XASSERTM(expr, msg)
Assertion macro definition with custom message.
Definition: assertion.hpp:263
base_header.hpp
FEAT Kernel base header.
FEAT::Geometry::CGALVerticesAroundFaceAdjactor
Adjactor for vertices around a face.
Definition: cgal.hpp:139
FEAT::Geometry::CGALWrapper
Wrapper for the CGAL Library.
Definition: cgal.hpp:173
FEAT::Geometry::CGALWrapper::intersects_line
bool intersects_line(const PointType &a, const PointType &b) const
tests whether the cgal mesh intersects with the line segment a->b
FEAT::Geometry::CGALWrapper::CGALWrapper
CGALWrapper(const CGALWrapper &)=delete
rule of five
FEAT::Tiny::Vector
Tiny Vector class template.
Definition: tiny_algebra.hpp:56
FEAT::Geometry
Geometry namespace.
Definition: adaptive_export_vtu.hpp:17
FEAT::Geometry::CGALFeatureNetwork
std::vector< CGALFeature > CGALFeatureNetwork
A FeatureNetwork is a list of features.
Definition: cgal.hpp:47
FEAT::Geometry::CGALFeature
std::vector< std::uint32_t > CGALFeature
A feature is an edge-path on a surface mesh, stored as a list of vertex indices.
Definition: cgal.hpp:40
FEAT::Geometry::PermutationStrategy::other
@ other
generic/other permutation strategy
FEAT::Tiny::dot
CUDA_HOST_DEVICE T_ dot(const T_ &a, const T_ &b)
Computes the dot-product of two scalars.
Definition: tiny_algebra.hpp:3001
FEAT::Index
std::uint64_t Index
Index data type.
Definition: base_header.hpp:122