@@ -180,7 +180,7 @@ angleBetweenVecs(const double vecA[3], const double vecB[3],
180
180
}
181
181
182
182
static double
183
- triArea_old (const double coords[][3 ], const GMGW_int v0, const GMGW_int v1,
183
+ triArea (const double coords[][3 ], const GMGW_int v0, const GMGW_int v1,
184
184
const GMGW_int v2)
185
185
{
186
186
double e0 = distance (coords[v0], coords[v1]);
@@ -201,7 +201,12 @@ triArea_old(const double coords[][3], const GMGW_int v0, const GMGW_int v1,
201
201
}
202
202
assert (e0 >= e1 );
203
203
assert (e1 >= e2 );
204
- assert (e2 >= (e0 -e1 )); // astonishingly numerically possible for 180° inner edge angle (jww)
204
+ if (e2 - (e0 -e1 ) < 0 ) {
205
+ // These three edge lengths can't be a real triangle, because e0
206
+ // is larger than the sum of e1 and e2. Treat this triangle as
207
+ // exactly linear, and return zero area.
208
+ return 0 ;
209
+ }
205
210
// The parentheses enforce a numerically stable order of operations.
206
211
double Area = 0.25
207
212
* sqrt (
@@ -211,7 +216,7 @@ triArea_old(const double coords[][3], const GMGW_int v0, const GMGW_int v1,
211
216
}
212
217
213
218
static double // replaces previous edge-length based version (jww)
214
- triArea (const double coords[][3 ], const GMGW_int v0, const GMGW_int v1,
219
+ triArea_alternative (const double coords[][3 ], const GMGW_int v0, const GMGW_int v1,
215
220
const GMGW_int v2)
216
221
{
217
222
// the following is unneeded computational overshoot - only for the assertion below
0 commit comments