GCC Code Coverage Report

Directory: ./
File: mesh/mesh_geometry.cpp
Date: 2026-09-27 03:24:14
Exec Total Coverage
Lines: 120 241 49.8%
Functions: 11 18 61.1%
Branches: 138 478 28.9%
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/mesh/mesh_geometry.h>
41 #include <geogram/delaunay/LFS.h>
42 #include <geogram/voronoi/CVT.h>
43 #include <geogram/basic/attributes.h>
44 #include <geogram/basic/geometry.h>
45 #include <geogram/basic/logger.h>
46
47 namespace {
48
49 using namespace GEO;
50
51 /**
52 * \brief Computes a sizing field using local feature size
53 * \details The sizing field is stored into the vertex weights
54 * of the mesh
55 * \param[in] M the mesh
56 * \param[in] LFS the local feature size
57 * \param[in] gradation power to be applied to the sizing field
58 */
59 1 void compute_sizing_field_lfs(
60 Mesh& M, const LocalFeatureSize& LFS, double gradation
61 ) {
62 // Avoid LFS points that are too close to the surface
63
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 double min_distance2 = 0.1*surface_average_edge_length(M);
64 1 min_distance2 = min_distance2 * min_distance2;
65
66
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 Attribute<double> weight(M.vertices.attributes(),"weight");
67
2/2
✓ Branch 4 taken 31762 times.
✓ Branch 5 taken 1 times.
 31763 for(index_t v: M.vertices) {
68
2/4
✓ Branch 1 taken 31762 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 31762 times.
✗ Branch 5 not taken.
 31762 double lfs2 = LFS.squared_lfs(M.vertices.point_ptr(v));
69 31762 lfs2 = std::max(lfs2, min_distance2);
70 31762 double w = pow(lfs2, -2.0 * gradation);
71
1/2
✓ Branch 1 taken 31762 times.
✗ Branch 2 not taken.
 31762 weight[v] = w;
72 }
73 1 }
74 }
75
76 /****************************************************************************/
77
78 namespace GEO {
79
80 namespace Geom {
81
82 275218 vec3 mesh_facet_normal(const Mesh& M, index_t f) {
83 275218 vec3 result(0.0, 0.0, 0.0);
84
7/12
✓ Branch 1 taken 275218 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 275218 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 275218 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 275218 times.
✗ Branch 11 not taken.
✓ Branch 16 taken 550436 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 275218 times.
✓ Branch 19 taken 275218 times.
 550436 for(auto [p1, p2, p3]: M.facets.triangle_points(f)) {
85
1/2
✓ Branch 5 taken 275218 times.
✗ Branch 6 not taken.
 275218 result += cross(p2 - p1, p3 - p1);
86 }
87 275218 return result;
88 }
89
90 ✗ double mesh_unsigned_normal_angle(
91 const Mesh& M, index_t f1, index_t f2
92 ) {
93 ✗ vec3 n1 = mesh_facet_normal(M,f1);
94 ✗ vec3 n2 = mesh_facet_normal(M,f2);
95 ✗ double l = length(n1)*length(n2);
96 ✗ double cos_angle = (l > 1e-30) ? dot(n1, n2)/l : 1.0;
97 // Numerical precision problem may occur, and generate
98 // normalized dot products that are outside the valid
99 // range of acos.
100 ✗ geo_clamp(cos_angle, -1.0, 1.0);
101 ✗ return acos(cos_angle);
102 }
103
104 ✗ double mesh_normal_angle(const Mesh& M, index_t c) {
105 ✗ geo_debug_assert(M.facets.are_simplices());
106 ✗ index_t f1 = c/3;
107 ✗ index_t f2 = M.facet_corners.adjacent_facet(c);
108 ✗ geo_debug_assert(f2 != NO_FACET);
109 ✗ vec3 n1 = mesh_facet_normal(M,f1);
110 ✗ vec3 n2 = mesh_facet_normal(M,f2);
111 ✗ double l = length(n1)*length(n2);
112 ✗ double sign = 1.0;
113 ✗ if(dot(cross(n1,n2),mesh_corner_vector(M,c)) > 0.0) {
114 ✗ sign = -1.0;
115 }
116 ✗ double cos_angle = (l > 1e-30) ? dot(n1, n2)/l : 1.0;
117 // Numerical precision problem may occur, and generate
118 // normalized dot products that are outside the valid
119 // range of acos.
120 ✗ geo_clamp(cos_angle, -1.0, 1.0);
121 ✗ return sign*acos(cos_angle);
122 }
123
124 64 double mesh_area(const Mesh& M, index_t dim) {
125 64 double result = 0.0;
126
2/2
✓ Branch 4 taken 547633 times.
✓ Branch 5 taken 64 times.
 547697 for(index_t f: M.facets) {
127
1/2
✓ Branch 1 taken 547633 times.
✗ Branch 2 not taken.
 547633 result += mesh_facet_area(M, f, dim);
128 }
129 64 return result;
130 }
131
132 7 double GEOGRAM_API mesh_enclosed_volume(const Mesh& M) {
133 7 double result = 0.0;
134
2/2
✓ Branch 5 taken 69733 times.
✓ Branch 6 taken 7 times.
 69740 for(index_t f: M.facets) {
135
7/12
✓ Branch 1 taken 69733 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 69733 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 69733 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 69733 times.
✗ Branch 11 not taken.
✓ Branch 16 taken 139466 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 69733 times.
✓ Branch 19 taken 69733 times.
 139466 for(auto [ p1, p2, p3 ] : M.facets.triangle_points(f)) {
136
1/2
✓ Branch 3 taken 69733 times.
✗ Branch 4 not taken.
 69733 result += dot(p1,cross(p2,p3)) / 6.0;
137 }
138 }
139 7 return ::fabs(result);
140 }
141
142 }
143
144 6 void compute_normals(Mesh& M) {
145
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 if(M.vertices.dimension() < 6) {
146 6 M.vertices.set_dimension(6);
147 } else {
148 ✗ for(index_t i: M.vertices) {
149 ✗ Geom::mesh_vertex_normal_ref(M, i) = vec3(0.0, 0.0, 0.0);
150 }
151 }
152
2/2
✓ Branch 5 taken 69827 times.
✓ Branch 6 taken 6 times.
 69833 for(index_t f: M.facets) {
153
1/2
✓ Branch 1 taken 69827 times.
✗ Branch 2 not taken.
 69827 vec3 N = Geom::mesh_facet_normal(M, f);
154
3/4
✓ Branch 1 taken 69827 times.
✗ Branch 2 not taken.
✓ Branch 8 taken 209481 times.
✓ Branch 9 taken 69827 times.
 279308 for(index_t corner: M.facets.corners(f)) {
155
1/2
✓ Branch 1 taken 209481 times.
✗ Branch 2 not taken.
 209481 index_t v = M.facet_corners.vertex(corner);
156
1/2
✓ Branch 1 taken 209481 times.
✗ Branch 2 not taken.
 209481 Geom::mesh_vertex_normal_ref(M, v) += N;
157 }
158 }
159
2/2
✓ Branch 4 taken 35534 times.
✓ Branch 5 taken 6 times.
 35540 for(index_t i: M.vertices) {
160
1/2
✓ Branch 1 taken 35534 times.
✗ Branch 2 not taken.
 35534 Geom::mesh_vertex_normal_ref(M, i) = normalize(
161
2/4
✓ Branch 1 taken 35534 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35534 times.
✗ Branch 5 not taken.
 35534 Geom::mesh_vertex_normal(M, i)
162 );
163 }
164 6 }
165
166 6 void simple_Laplacian_smooth(Mesh& M, index_t nb_iter, bool normals_only) {
167
2/8
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 6 geo_assert(M.vertices.dimension() >= 6);
168
1/2
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 12 std::vector<vec3> p(M.vertices.nb());
169
1/2
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 6 std::vector<double> c(M.vertices.nb());
170
171
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 6 times.
 24 for(index_t k = 0; k < nb_iter; k++) {
172
1/2
✓ Branch 3 taken 18 times.
✗ Branch 4 not taken.
 18 p.assign(M.vertices.nb(), vec3(0.0, 0.0, 0.0));
173
1/2
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
 18 c.assign(M.vertices.nb(), 0);
174
2/2
✓ Branch 5 taken 209481 times.
✓ Branch 6 taken 18 times.
 209499 for(index_t f: M.facets) {
175
1/2
✓ Branch 1 taken 209481 times.
✗ Branch 2 not taken.
 209481 index_t b = M.facets.corners_begin(f);
176
1/2
✓ Branch 1 taken 209481 times.
✗ Branch 2 not taken.
 209481 index_t e = M.facets.corners_end(f);
177
2/2
✓ Branch 0 taken 628443 times.
✓ Branch 1 taken 209481 times.
 837924 for(index_t c1 = b; c1 != e; c1++) {
178
2/2
✓ Branch 0 taken 418962 times.
✓ Branch 1 taken 209481 times.
 628443 index_t c2 = (c1 == e - 1) ? b : c1 + 1;
179
1/2
✓ Branch 1 taken 628443 times.
✗ Branch 2 not taken.
 628443 index_t v1 = M.facet_corners.vertex(c1);
180
1/2
✓ Branch 1 taken 628443 times.
✗ Branch 2 not taken.
 628443 index_t v2 = M.facet_corners.vertex(c2);
181
2/2
✓ Branch 0 taken 315144 times.
✓ Branch 1 taken 313299 times.
 628443 if(v1 < v2) {
182 315144 double a = 1.0;
183 315144 c[v1] += a;
184 315144 c[v2] += a;
185
1/2
✓ Branch 0 taken 315144 times.
✗ Branch 1 not taken.
 315144 if(normals_only) {
186
1/2
✓ Branch 1 taken 315144 times.
✗ Branch 2 not taken.
 315144 p[v1] += a * Geom::mesh_vertex_normal(M, v2);
187
1/2
✓ Branch 1 taken 315144 times.
✗ Branch 2 not taken.
 315144 p[v2] += a * Geom::mesh_vertex_normal(M, v1);
188 } else {
189 ✗ p[v1] += a * M.vertices.point(v2);
190 ✗ p[v2] += a * M.vertices.point(v1);
191 }
192 }
193 }
194 }
195
2/2
✓ Branch 5 taken 106602 times.
✓ Branch 6 taken 18 times.
 106620 for(index_t v: M.vertices) {
196
1/2
✓ Branch 0 taken 106602 times.
✗ Branch 1 not taken.
 106602 if(normals_only) {
197
1/2
✓ Branch 2 taken 106602 times.
✗ Branch 3 not taken.
 106602 double l = length(p[v]);
198
1/2
✓ Branch 0 taken 106602 times.
✗ Branch 1 not taken.
 106602 if(l > 1e-30) {
199
1/2
✓ Branch 3 taken 106602 times.
✗ Branch 4 not taken.
 106602 Geom::mesh_vertex_normal_ref(M,v) = (1.0 / l) * p[v];
200 }
201 } else {
202 ✗ M.vertices.point(v) = 1.0 / c[v] * p[v];
203 }
204 }
205 }
206
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 6 if(!normals_only) {
207 ✗ compute_normals(M);
208 }
209 6 }
210
211 84 void get_bbox(const Mesh& M, double* xyzmin, double* xyzmax) {
212
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 84 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 84 geo_assert(M.vertices.dimension() >= 3);
213
2/2
✓ Branch 0 taken 252 times.
✓ Branch 1 taken 84 times.
 336 for(index_t c = 0; c < 3; c++) {
214 252 xyzmin[c] = Numeric::max_float64();
215 252 xyzmax[c] = Numeric::min_float64();
216 }
217
6/10
✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 84 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 84 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 307852 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 307852 times.
✓ Branch 15 taken 84 times.
 307936 for(const vec3& p: M.vertices.points()) {
218
2/2
✓ Branch 0 taken 923556 times.
✓ Branch 1 taken 307852 times.
 1231408 for(index_t c = 0; c < 3; c++) {
219
1/2
✓ Branch 1 taken 923556 times.
✗ Branch 2 not taken.
 923556 xyzmin[c] = std::min(xyzmin[c], p[c]);
220
1/2
✓ Branch 1 taken 923556 times.
✗ Branch 2 not taken.
 923556 xyzmax[c] = std::max(xyzmax[c], p[c]);
221 }
222 }
223 84 }
224
225 82 double bbox_diagonal(const Mesh& M) {
226
2/8
✓ Branch 1 taken 82 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 82 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 82 geo_assert(M.vertices.dimension() >= 3);
227 double xyzmin[3];
228 double xyzmax[3];
229
1/2
✓ Branch 1 taken 82 times.
✗ Branch 2 not taken.
 82 get_bbox(M, xyzmin, xyzmax);
230 82 return ::sqrt(
231 82 geo_sqr(xyzmax[0] - xyzmin[0]) +
232 82 geo_sqr(xyzmax[1] - xyzmin[1]) +
233 82 geo_sqr(xyzmax[2] - xyzmin[2])
234 164 );
235 }
236
237 6 void set_anisotropy(Mesh& M, double s) {
238
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 6 times.
 6 if(M.vertices.dimension() < 6) {
239 ✗ compute_normals(M);
240 }
241
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 6 if(s == 0.0) {
242 ✗ unset_anisotropy(M);
243 ✗ return;
244 }
245 6 s *= bbox_diagonal(M);
246
2/2
✓ Branch 4 taken 35534 times.
✓ Branch 5 taken 6 times.
 35540 for(index_t i: M.vertices) {
247
1/2
✓ Branch 1 taken 35534 times.
✗ Branch 2 not taken.
 35534 Geom::mesh_vertex_normal_ref(M, i) =
248
2/4
✓ Branch 1 taken 35534 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35534 times.
✗ Branch 5 not taken.
 71068 s * normalize(Geom::mesh_vertex_normal(M, i));
249 }
250 }
251
252 ✗ void unset_anisotropy(Mesh& M) {
253 ✗ if(M.vertices.dimension() < 6) {
254 ✗ return;
255 }
256 ✗ for(index_t i: M.vertices) {
257 ✗ Geom::mesh_vertex_normal_ref(M, i) = normalize(
258 ✗ Geom::mesh_vertex_normal(M, i)
259 );
260 }
261 }
262
263 1 void compute_sizing_field(
264 Mesh& M, double gradation, index_t nb_lfs_samples
265 ) {
266
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if(nb_lfs_samples != 0) {
267
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
 2 Logger::out("LFS") << "Sampling surface" << std::endl;
268
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 CentroidalVoronoiTesselation CVT(&M, 3);
269
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 CVT.compute_initial_sampling(nb_lfs_samples);
270
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
 1 Logger::out("LFS") << "Optimizing sampling (Lloyd)" << std::endl;
271
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 CVT.Lloyd_iterations(5);
272
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
 1 Logger::out("LFS") << "Optimizing sampling (Newton)" << std::endl;
273
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 CVT.Newton_iterations(10);
274
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
 1 Logger::out("LFS") << "Computing medial axis" << std::endl;
275
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 1 LocalFeatureSize LFS(CVT.nb_points(), CVT.embedding(0));
276
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
 1 Logger::out("LFS") << "Computing sizing field" << std::endl;
277
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 compute_sizing_field_lfs(M, LFS, gradation);
278 1 } else {
279 ✗ if(M.vertices.dimension() == 3) {
280 ✗ LocalFeatureSize LFS(M.vertices.nb(), M.vertices.point_ptr(0));
281 ✗ compute_sizing_field_lfs(M, LFS, gradation);
282 ✗ } else {
283 ✗ std::vector<double> pts;
284 ✗ pts.reserve(M.vertices.nb() * 3);
285 ✗ for(index_t v: M.vertices) {
286 ✗ pts.push_back(M.vertices.point_ptr(v)[0]);
287 ✗ pts.push_back(M.vertices.point_ptr(v)[1]);
288 ✗ pts.push_back(M.vertices.point_ptr(v)[2]);
289 }
290 ✗ LocalFeatureSize LFS(M.vertices.nb(), pts.data());
291 ✗ compute_sizing_field_lfs(M, LFS, gradation);
292 ✗ }
293 }
294 1 }
295
296 ✗ void normalize_embedding_area(Mesh& M) {
297 ✗ if(M.vertices.dimension() == 3) {
298 ✗ if(M.vertices.attributes().is_defined("weight")) {
299 ✗ M.vertices.attributes().delete_attribute_store("weight");
300 }
301 ✗ return;
302 }
303 ✗ Attribute<double> weight(M.vertices.attributes(), "weight");
304 ✗ std::vector<double> area3d(M.vertices.nb(), 0.0);
305 ✗ std::vector<double> areaNd(M.vertices.nb(), 0.0);
306 ✗ for(index_t f: M.facets) {
307 ✗ double A3d = Geom::mesh_facet_area(M, f, 3);
308 ✗ double ANd = Geom::mesh_facet_area(M, f);
309 ✗ for(index_t c: M.facets.corners(f)) {
310 ✗ index_t v = M.facet_corners.vertex(c);
311 ✗ area3d[v] += A3d;
312 ✗ areaNd[v] += ANd;
313 }
314 }
315 ✗ for(index_t v: M.vertices) {
316 ✗ double A3d = area3d[v];
317 ✗ double ANd = areaNd[v];
318 ✗ ANd = std::max(ANd, 1e-6);
319 ✗ double w = ::pow(A3d / ANd, 2.0);
320 ✗ weight[v] = w;
321 }
322 ✗ }
323
324 ✗ double mesh_cell_volume(
325 const Mesh& M, index_t c
326 ) {
327 ✗ geo_assert(M.vertices.dimension() >= 3);
328
329 // Connectors are virtual cells, they
330 // do not have a geometry.
331 ✗ if(M.cells.type(c) == MESH_CONNECTOR) {
332 ✗ return 0.0;
333 }
334
335 // Easy case: tetrahedra.
336 ✗ if(M.cells.type(c) == MESH_TET) {
337 ✗ const vec3& p0 = M.cells.point(c,0);
338 ✗ const vec3& p1 = M.cells.point(c,1);
339 ✗ const vec3& p2 = M.cells.point(c,2);
340 ✗ const vec3& p3 = M.cells.point(c,3);
341 ✗ return ::fabs(Geom::tetra_signed_volume(p0,p1,p2,p3));
342 }
343
344 // Arbitrary cells are decomposed into tetrahedra, with one
345 // vertex in the center of the cell, and one vertex in the
346 // center of each face. This ensures that two adjacent cells
347 // are not overlapping (if simply triangulating the faces,
348 // there would be tiny overlaps / tiny gaps that would introduce
349 // errors in the total volume of the mesh).
350 //
351 // Note that the center point may fall outside the cell in some
352 // degenerate configurations. It is not a problem since we compute
353 // signed tetrahedra volumes, in such a way that overlapping volumes
354 // will cancel-out in such a configuration.
355 // Therefore, we could take an arbitrary point as the enter point,
356 // including the origin, but taking the center probably makes
357 // computations more stable, by cancelling the translations.
358
359 ✗ double result = 0.0;
360 ✗ vec3 center{0.0, 0.0, 0.0};
361 ✗ index_t nbcv = M.cells.nb_vertices(c);
362 ✗ for(const vec3& p: M.cells.points(c)) {
363 ✗ center += p;
364 }
365 ✗ center /= double(nbcv);
366 ✗ for(index_t lf=0; lf<M.cells.nb_facets(c); ++lf) {
367 ✗ index_t nbcfv = M.cells.facet_nb_vertices(c,lf);
368 ✗ if(nbcfv == 3) {
369 ✗ const vec3& p1 = M.vertices.point(M.cells.facet_vertex(c,lf,0));
370 ✗ const vec3& p2 = M.vertices.point(M.cells.facet_vertex(c,lf,1));
371 ✗ const vec3& p3 = M.vertices.point(M.cells.facet_vertex(c,lf,2));
372 ✗ result += Geom::tetra_signed_volume(center, p1, p2, p3);
373 } else {
374 ✗ vec3 facet_center {0.0, 0.0, 0.0};
375 ✗ for(index_t lfv=0; lfv<nbcfv; ++lfv) {
376 facet_center +=
377 ✗ M.vertices.point(M.cells.facet_vertex(c,lf,lfv));
378 }
379 ✗ facet_center /= double(nbcfv);
380 ✗ for(index_t lfv1=0; lfv1<nbcfv; ++lfv1) {
381 ✗ index_t lfv2 = (lfv1 + 1) % nbcfv;
382 const vec3& p1 =
383 ✗ M.vertices.point(M.cells.facet_vertex(c,lf,lfv1));
384 const vec3& p2 =
385 ✗ M.vertices.point(M.cells.facet_vertex(c,lf,lfv2));
386 ✗ result +=
387 ✗ Geom::tetra_signed_volume(center, facet_center, p1, p2);
388 }
389 }
390 }
391 ✗ return ::fabs(result);
392 }
393
394
395 ✗ double mesh_cells_volume(const Mesh& M) {
396 ✗ double result = 0.0;
397 ✗ for(index_t c: M.cells) {
398 ✗ result += mesh_cell_volume(M,c);
399 }
400 ✗ return result;
401 }
402
403 ✗ vec3 GEOGRAM_API mesh_cell_facet_normal(
404 const Mesh& M, index_t c, index_t lf
405 ) {
406 ✗ geo_debug_assert(M.vertices.dimension() >= 3);
407 ✗ geo_debug_assert(c < M.cells.nb());
408 ✗ geo_debug_assert(lf < M.cells.nb_facets(c));
409
410 ✗ index_t v1 = M.cells.facet_vertex(c,lf,0);
411 ✗ index_t v2 = M.cells.facet_vertex(c,lf,1);
412 ✗ index_t v3 = M.cells.facet_vertex(c,lf,2);
413
414 ✗ const vec3& p1 = M.vertices.point(v1);
415 ✗ const vec3& p2 = M.vertices.point(v2);
416 ✗ const vec3& p3 = M.vertices.point(v3);
417
418 ✗ return cross(p2 - p1, p3 - p1);
419 }
420
421
422
423 3 double surface_average_edge_length(const Mesh& M) {
424 3 double result = 0.0;
425 3 index_t count = 0;
426
2/2
✓ Branch 5 taken 82200 times.
✓ Branch 6 taken 3 times.
 82203 for(index_t f: M.facets) {
427
3/4
✓ Branch 1 taken 82200 times.
✗ Branch 2 not taken.
✓ Branch 7 taken 246600 times.
✓ Branch 8 taken 82200 times.
 328800 for(index_t c1: M.facets.corners(f)) {
428
1/2
✓ Branch 1 taken 246600 times.
✗ Branch 2 not taken.
 246600 index_t c2 = M.facets.next_corner_around_facet(f,c1);
429
1/2
✓ Branch 1 taken 246600 times.
✗ Branch 2 not taken.
 246600 index_t v1 = M.facet_corners.vertex(c1);
430
1/2
✓ Branch 1 taken 246600 times.
✗ Branch 2 not taken.
 246600 index_t v2 = M.facet_corners.vertex(c2);
431
2/4
✓ Branch 1 taken 246600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 246600 times.
✗ Branch 5 not taken.
 246600 result += Geom::distance(
432 M.vertices.point_ptr(v1),
433 M.vertices.point_ptr(v2),
434
2/4
✓ Branch 1 taken 246600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 246600 times.
✗ Branch 5 not taken.
 246600 coord_index_t(M.vertices.dimension())
435 );
436 246600 ++count;
437 }
438 }
439
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 if(count != 0) {
440 3 result /= double(count);
441 }
442 3 return result;
443 }
444 }
445