Description
centroid_3DPolygon() in ComputationalGeometry.hpp calls GEOS_ERROR when a polygon has zero area. This aborts the run on corner-point meshes that contain pinch-outs. A pinch-out is ordinary reservoir geology, not a defect in the mesh.
Regression
Before #3977 (01a4f15b2e, "Change to local Newell's normal") the branch returned zero:
if( area > areaTolerance ) { /* normalize */ }
else if( area < -areaTolerance ) { GEOS_ERROR( "Negative area found" ); }
else { return 0.0; }
#3977 replaced that return 0.0 with GEOS_ERROR( "Null area found" ). The return 0.0 looks deliberate: a zero-area face gives zero transmissibility, which is the physically correct outcome, because no flow crosses a face of zero area. #3977 was about Newell normals, so the change to the tolerance branch appears incidental.
Concrete case
A corner-point mesh of 1761 hexahedra, converted from a standard Eclipse-style deck. Eight consecutive cells each have one vertical edge of zero length, while the other three measure 3.9 to 4.7 m. One lateral face of each cell collapses to a triangle.
The cells are valid and active: non-zero volume, and porosity between 0.05 and 0.14. The mesh holds 8 degenerate faces among several thousand.
On current develop this mesh cannot be run at all. The first degenerate face aborts initialization.
Suggested fix
Change GEOS_ERROR to GEOS_WARNING in the null-area branch and keep return 0.0. This restores the pre-#3977 behavior and keeps the diagnostic that #3977 added. The negative-area branch should stay a GEOS_ERROR, because a negative area does signal a real defect.
A stricter alternative would be an input flag that lets the user choose between abort and warn, with warn as the default for corner-point meshes.
Impact
With the warning, two such decks run to completion and reproduce our reference results bit for bit over 500 simulated years. Each reports 16 warnings, from the 8 degenerate faces.
Environment
GEOS develop at 24a8410, gcc 13.3.0, VTK 9.4.2, serial run.
Description
centroid_3DPolygon()inComputationalGeometry.hppcallsGEOS_ERRORwhen a polygon has zero area. This aborts the run on corner-point meshes that contain pinch-outs. A pinch-out is ordinary reservoir geology, not a defect in the mesh.Regression
Before #3977 (
01a4f15b2e, "Change to local Newell's normal") the branch returned zero:#3977 replaced that
return 0.0withGEOS_ERROR( "Null area found" ). Thereturn 0.0looks deliberate: a zero-area face gives zero transmissibility, which is the physically correct outcome, because no flow crosses a face of zero area. #3977 was about Newell normals, so the change to the tolerance branch appears incidental.Concrete case
A corner-point mesh of 1761 hexahedra, converted from a standard Eclipse-style deck. Eight consecutive cells each have one vertical edge of zero length, while the other three measure 3.9 to 4.7 m. One lateral face of each cell collapses to a triangle.
The cells are valid and active: non-zero volume, and porosity between 0.05 and 0.14. The mesh holds 8 degenerate faces among several thousand.
On current
developthis mesh cannot be run at all. The first degenerate face aborts initialization.Suggested fix
Change
GEOS_ERRORtoGEOS_WARNINGin the null-area branch and keepreturn 0.0. This restores the pre-#3977 behavior and keeps the diagnostic that #3977 added. The negative-area branch should stay aGEOS_ERROR, because a negative area does signal a real defect.A stricter alternative would be an input flag that lets the user choose between abort and warn, with warn as the default for corner-point meshes.
Impact
With the warning, two such decks run to completion and reproduce our reference results bit for bit over 500 simulated years. Each reports 16 warnings, from the 8 degenerate faces.
Environment
GEOS
developat 24a8410, gcc 13.3.0, VTK 9.4.2, serial run.