GCC Code Coverage Report


Directory: ./
File: lib/geogram/delaunay/delaunay_tetgen.cpp
Date: 2026-09-07 02:28:19
Exec Total Coverage
Lines: 88 133 66.2%
Functions: 5 7 71.4%
Branches: 51 132 38.6%

Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #include <geogram/delaunay/delaunay_tetgen.h>
41 #include <geogram/mesh/mesh.h>
42 #include <geogram/basic/logger.h>
43 #include <geogram/basic/command_line.h>
44 #include <geogram/bibliography/bibliography.h>
45
46 #include <geogram/mesh/mesh_io.h>
47
48 #ifdef GEOGRAM_WITH_TETGEN
49
50 namespace GEO {
51
52 1 DelaunayTetgen::DelaunayTetgen(coord_index_t dimension) :
53 1 Delaunay(3)
54 {
55
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if(dimension != 3) {
56 throw InvalidDimension(dimension, "DelaunayTetgen", "3");
57 }
58
59
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 geo_cite("DBLP:journals/toms/Si15");
60
61 1 tetgen_in_.firstnumber = 0;
62 1 tetgen_out_.firstnumber = 0;
63 1 }
64
65 4 DelaunayTetgen::~DelaunayTetgen() {
66 // We do not want that the destructor of tetgen_in_
67 // deallocates the points, since they are managed
68 // by a vector<>.
69 2 tetgen_in_.pointlist = nullptr;
70 2 tetgen_in_.numberofpoints = 0;
71 8 }
72
73 index_t DelaunayTetgen::region(index_t t) const {
74 geo_debug_assert(keep_regions_);
75 geo_debug_assert(t < nb_cells());
76 return index_t(tetgen_out_.tetrahedronattributelist[t]);
77 }
78
79 1 bool DelaunayTetgen::supports_constraints() const {
80 1 return true;
81 }
82
83 1 void DelaunayTetgen::set_vertices(
84 index_t nb_vertices, const double* vertices
85 ) {
86
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 if(constraints_ != nullptr) {
87 1 set_vertices_constrained(nb_vertices, vertices);
88 } else {
89 set_vertices_unconstrained(nb_vertices, vertices);
90 }
91 1 }
92
93 void DelaunayTetgen::set_vertices_unconstrained(
94 index_t nb_vertices, const double* vertices
95 ) {
96 // Q: quiet
97 // n: output tet neighbors
98 // V: verbose
99 if(CmdLine::get_arg_bool("dbg:tetgen")) {
100 tetgen_args_.parse_commandline(const_cast<char*>("Vn"));
101 } else {
102 tetgen_args_.parse_commandline(const_cast<char*>("Qn"));
103 }
104
105 Delaunay::set_vertices(nb_vertices, vertices);
106 tetgen_out_.deinitialize();
107 tetgen_in_.initialize();
108 tetgen_in_.numberofpoints = int(nb_vertices);
109 tetgen_in_.pointlist = const_cast<double*>(vertices);
110 try {
111 GEO_3rdParty::tetrahedralize(
112 &tetgen_args_, &tetgen_in_, &tetgen_out_
113 );
114 }
115 catch(int error_code) {
116 Logger::err("DelaunayTetgen")
117 << "Encountered a problem..." << std::endl;
118 throw Delaunay::InvalidInput(error_code);
119 }
120 set_arrays(
121 index_t(tetgen_out_.numberoftetrahedra),
122 reinterpret_cast<index_t*>(tetgen_out_.tetrahedronlist),
123 reinterpret_cast<index_t*>(tetgen_out_.neighborlist)
124 );
125 }
126
127
128 1 void DelaunayTetgen::set_vertices_constrained(
129 index_t nb_vertices, const double* vertices
130 ) {
131 index_t nb_borders = 0;
132
2/2
✓ Branch 0 taken 12000 times.
✓ Branch 1 taken 1 times.
12001 for(index_t c=0; c<constraints_->facet_corners.nb(); ++c) {
133
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12000 times.
12000 if(constraints_->facet_corners.adjacent_facet(c) == NO_FACET) {
134 ++nb_borders;
135 }
136 }
137
138
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if(nb_borders != 0) {
139 Logger::warn("DelaunayTetgen")
140 << "Constraints have " << nb_borders
141 << " edge(s) on the border"
142 << std::endl;
143 }
144
145 1 tetgen_out_.deinitialize();
146
147 // Q: quiet
148 // p: input data is surfacic
149 // n: output tet neighbors
150 // q: desired quality
151 // O0: do not optimize mesh
152 // V: verbose
153 // YY: prohibit steiner points on boundaries
154 // (first Y for exterior boundary, second Y for the
155 // other ones).
156 // AA: generate region tags for each shell.
157 // M: do not merge coplanar facets
158
159 std::string cmdline;
160
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 if(refine_) {
161
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
2 if(CmdLine::get_arg_bool("dbg:tetgen")) {
162 cmdline = "Vpnq" + String::to_string(quality_) + "YYAA";
163 } else {
164
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
3 cmdline = "Qpnq" + String::to_string(quality_) + "YYAA";
165 }
166 } else {
167 if(CmdLine::get_arg_bool("dbg:tetgen")) {
168 cmdline = "VpnO0YYAA";
169 } else {
170 cmdline = "QpnO0YYAA";
171 }
172 }
173
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 tetgen_args_.parse_commandline(const_cast<char*>(cmdline.c_str()));
174
175 1 tetgen_in_.deinitialize();
176 1 tetgen_in_.initialize();
177 1 tetgen_in_.firstnumber = 0 ;
178
179 //
180 // Copy vertices
181 //
182
183 1 tetgen_in_.numberofpoints = int(
184
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 constraints_->vertices.nb()+nb_vertices
185 );
186
3/6
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
1 tetgen_in_.pointlist = new double[3*tetgen_in_.numberofpoints];
187
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 if(constraints_->vertices.nb() != 0) {
188 1 Memory::copy(
189 tetgen_in_.pointlist, constraints_->vertices.point_ptr(0),
190 1 constraints_->vertices.nb()*3*sizeof(double)
191 );
192 }
193
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if(nb_vertices != 0) {
194 Memory::copy(
195 &tetgen_in_.pointlist[3*constraints_->vertices.nb()],
196 vertices, nb_vertices*3*sizeof(double)
197 );
198 }
199
200 // Edges constraints
201 // (no need to copy, we make tetgen_in_
202 // point to the edges of the input
203 // constraints mesh)
204
205
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if(constraints_->edges.nb() != 0) {
206 tetgen_in_.numberofedges = int(
207 constraints_->edges.nb()
208 );
209 tetgen_in_.edgelist = const_cast<int*>(
210 (const int*)constraints_->edges.vertex_index_ptr(0)
211 );
212 }
213
214 // Copy facet constraints
215 //
216 // All the polygons are allocated in one go, in a contiguous array.
217
218 GEO_3rdParty::tetgenio::polygon* polygons =
219
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 new GEO_3rdParty::tetgenio::polygon[constraints_->facets.nb()];
220 1 tetgen_in_.numberoffacets = int(constraints_->facets.nb()) ;
221 1 tetgen_in_.facetlist =
222
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
1 new GEO_3rdParty::tetgenio::facet[tetgen_in_.numberoffacets];
223
2/2
✓ Branch 0 taken 4000 times.
✓ Branch 1 taken 1 times.
4001 for(index_t f=0; f<constraints_->facets.nb(); ++f) {
224
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4000 times.
4000 GEO_3rdParty::tetgenio::facet& F = tetgen_in_.facetlist[f];
225 GEO_3rdParty::tetgenio::init(&F);
226 4000 F.numberofpolygons = 1;
227
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4000 times.
4000 F.polygonlist = &polygons[f];
228 GEO_3rdParty::tetgenio::polygon& P = F.polygonlist[0];
229 GEO_3rdParty::tetgenio::init(&P) ;
230
1/2
✓ Branch 0 taken 4000 times.
✗ Branch 1 not taken.
4000 P.numberofvertices = int(constraints_->facets.nb_vertices(f));
231 4000 P.vertexlist = reinterpret_cast<int*>(
232 const_cast<Mesh*>(constraints_)->facet_corners.vertex_index_ptr(
233 constraints_->facets.corners_begin(f)
234 )
235 );
236 4000 F.numberofholes = 0 ;
237 4000 F.holelist = nullptr ;
238 }
239
240 bool there_was_an_error = false;
241 int error_code = 0;
242
243 try {
244
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 GEO_3rdParty::tetrahedralize(
245 &tetgen_args_, &tetgen_in_, &tetgen_out_
246 );
247 } catch(...) {
248 Logger::err("DelaunayTetgen")
249 << "Encountered a problem..."
250 << std::endl;
251 there_was_an_error = true;
252 /*
253 Logger::err("DelaunayTetgen")
254 << "Encountered a problem, relaunching in diagnose mode..."
255 << std::endl;
256 cmdline += "d";
257 tetgen_args_.parse_commandline(const_cast<char*>(cmdline.c_str()));
258 try {
259 GEO_3rdParty::tetrahedralize(
260 &tetgen_args_, &tetgen_in_, &tetgen_out_
261 );
262 } catch(...) {
263 }
264 */
265 }
266
267
268 // Deallocate the datastructures used by tetgen,
269 // and disconnect them from tetgen,
270 // so that tetgen does not try to deallocate them.
271
272 // Pointlist was allocated in local array
273 1 tetgen_in_.numberofpoints = 0;
274
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 delete[] tetgen_in_.pointlist;
275 1 tetgen_in_.pointlist = nullptr;
276
277 // Edges were shared with constraint mesh
278 // (no need to deallocate)
279 1 tetgen_in_.numberofedges = 0;
280 1 tetgen_in_.edgelist = nullptr;
281
282 // Facets structures were allocated in local
283 // array, and vertices indices were shared
284 // with constraint mesh
285
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 delete[] tetgen_in_.facetlist;
286 1 tetgen_in_.facetlist = nullptr;
287 1 tetgen_in_.numberoffacets = 0;
288 1 delete[] polygons;
289
290
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
1 if(there_was_an_error) {
291 InvalidInput error_report(error_code);
292 for(auto f : tetgen_in_.isectfaces) {
293 // Note: tetgen reports facets with
294 // 1-based indexing ([1...nf]) !!
295 error_report.invalid_facets.push_back(index_t(f)-1);
296 }
297 Logger::err("DelaunayTetgen")
298 << "Found "
299 << error_report.invalid_facets.size()
300 << " facets with intersections." << std::endl;
301 throw(error_report);
302 }
303
304 1 index_t nb_tets = index_t(tetgen_out_.numberoftetrahedra);
305
306
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 if(!keep_regions_) {
307 // Determine which regions are incident to
308 // the 'exterior' (neighbor = -1 or tet is adjacent to
309 // a tet in region 0).
310 // The region Id of tet t is determined by:
311 // tetgen_out_.tetrahedronattributelist[t]
312
313 std::set<double> good_regions;
314 10536 for(
315 index_t t = 0;
316
2/2
✓ Branch 0 taken 10535 times.
✓ Branch 1 taken 1 times.
10536 t < index_t(tetgen_out_.numberoftetrahedra); ++t
317 ) {
318
2/2
✓ Branch 0 taken 39757 times.
✓ Branch 1 taken 6730 times.
46487 for(index_t f=0; f<4; ++f) {
319 39757 signed_index_t n = (tetgen_out_.neighborlist[t*4+f]);
320
2/2
✓ Branch 0 taken 35952 times.
✓ Branch 1 taken 3805 times.
39757 if(
321 35952 n == -1 ||
322
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 35952 times.
35952 tetgen_out_.tetrahedronattributelist[n] == 0.0
323 ) {
324 good_regions.insert(
325
1/2
✓ Branch 1 taken 3805 times.
✗ Branch 2 not taken.
3805 tetgen_out_.tetrahedronattributelist[t]
326 );
327 break;
328 }
329 }
330 }
331
332 // Remove the tets that are not in good_region.
333 vector<index_t> old2new(
334 index_t(tetgen_out_.numberoftetrahedra),NO_INDEX
335 );
336 nb_tets = 0;
337 10536 for(
338 index_t t = 0;
339
2/2
✓ Branch 0 taken 10535 times.
✓ Branch 1 taken 1 times.
10536 t < index_t(tetgen_out_.numberoftetrahedra); ++t
340 ) {
341
1/2
✓ Branch 0 taken 10535 times.
✗ Branch 1 not taken.
10535 if(
342 good_regions.find(
343 10535 tetgen_out_.tetrahedronattributelist[t]
344 ) != good_regions.end()
345 ) {
346
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10535 times.
10535 if(t != nb_tets) {
347 Memory::copy(
348 &tetgen_out_.tetrahedronlist[nb_tets * 4],
349 &tetgen_out_.tetrahedronlist[t * 4],
350 4 * sizeof(signed_index_t)
351 );
352 Memory::copy(
353 &tetgen_out_.neighborlist[nb_tets * 4],
354 &tetgen_out_.neighborlist[t * 4],
355 4 * sizeof(signed_index_t)
356 );
357 }
358 10535 old2new[t] = nb_tets;
359 10535 ++nb_tets;
360 }
361 }
362
2/2
✓ Branch 0 taken 42140 times.
✓ Branch 1 taken 1 times.
42141 for(index_t i = 0; i < 4 * nb_tets; ++i) {
363 42140 signed_index_t t = tetgen_out_.neighborlist[i];
364
2/2
✓ Branch 0 taken 38140 times.
✓ Branch 1 taken 4000 times.
42140 if(t != -1) {
365 38140 t = signed_index_t(old2new[t]);
366 }
367 42140 tetgen_out_.neighborlist[i] = t;
368 }
369 }
370
371 // Link tetgen's output to Delaunay class data structures.
372
373 1 Delaunay::set_vertices(
374
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 index_t(tetgen_out_.numberofpoints), tetgen_out_.pointlist
375 );
376
377 1 set_arrays(
378 nb_tets,
379 1 reinterpret_cast<index_t*>(tetgen_out_.tetrahedronlist),
380
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 reinterpret_cast<index_t*>(tetgen_out_.neighborlist)
381 );
382 1 }
383
384
385 }
386
387 #else
388
389 // Declare a dummy variable so that
390 // MSVC does not complain that it
391 // generated an empty object file.
392 extern int dummy_delaunay_tetgen_compiled;
393 int dummy_delaunay_tetgen_compiled = 1;
394
395 #endif
396