Make the pinch and distorted element removal consistent

README.md

@@ -161,7 +161,9 @@ options:
   --translate TRANSLATE TRANSLATE
                         Translate the (x, y) coordinates by this amount
   --mapaxes             Use the MAPAXES coordinates for an Eclipse file
-  --pinch               Remove pinched elements
+  --pinch               Remove pinched elements (coincident corners within
+                        --pinch-tol). Off by default; faithful conversion is
+                        produced without this flag.
   --pinch-tol PINCH_TOL
                         Tolerance for coincident corners when removing pinched
                         elements (default: 1e-3)
@@ -202,6 +204,10 @@ options:
   --strict-jacobians    Treat any non-positive element Jacobian as a fatal
                         error and exit non-zero. By default such elements only
                         produce a warning. Useful for CI / scripted workflows.
+  --remove-distorted    Remove elements with non-positive Jacobians (degenerate
+                        or inverted) from the output mesh, reporting a count of
+                        those removed. By default such elements are kept and
+                        only a warning is printed.
 ```
 
 ### Lateral refinement (Eclipse only)
@@ -311,32 +317,43 @@ A few practical notes:
 
 ### Element Jacobian check
 
-After the mesh is built, em2ex automatically computes the Jacobian at all 8 corners of every HEX8 element and prints a one-line summary:
+After conversion, em2ex evaluates the Jacobian at all 8 corners of every HEX8 element and prints a one-line summary:
 
 ```
 Element Jacobian check: 1000000 / 1000000 elements OK
 ```
 
-A non-positive Jacobian (negative = inverted element, zero = degenerate) almost always means one of two things:
+A non-positive Jacobian (zero = degenerate, negative = inverted) almost always means one of two things:
 
 - **The input file uses a `z increases upward` convention** rather than Eclipse's default `z increases downward`. The cells come out "upside down" and need `--flip` to invert the z values.
-- **There's something genuinely wrong with the input geometry** — corner ordering inconsistent within a cell, self-intersecting cells that survived pinch removal, etc.
+- **There are a few genuinely distorted cells** — for example near faults that are so skewed the element is inverted, or near-pinched cells that slipped through the tolerance filter.
 
-The relevance is downstream: most finite-element solvers (MOOSE, libMesh, MFEM, deal.II) reject elements with non-positive Jacobian at the assembly stage. Catching the problem at conversion time is cheaper than chasing it through a failed simulation.
+When all elements are non-positive the output is expanded with a hint:
 
-When the check fails, em2ex prints an expanded report with element IDs and centroid locations:
+```
+Element Jacobian check: 27 negative, 0 zero, 0 OK (out of 27)
+  All elements have non-positive Jacobians. Possible causes:
+    - z-up coordinate convention: try --flip to invert the z-axis
+    - Orientation-reversing coordinate system (e.g. MAPAXES handedness)
+  Use --remove-distorted to remove these elements and proceed anyway.
+```
+
+When only some elements are bad, the report includes element IDs and centroid locations:
 
 ```
-Element Jacobian check: 1 negative, 0 zero, 0 OK (out of 1)
-  Examples of negative-Jacobian elements (showing up to 5 of 1):
-    element 1: centroid (0.5, 0.5, 0.5), min Jacobian = -1.000e+00
+Element Jacobian check: 3 negative, 0 zero, 24 OK (out of 27)
+  Examples of negative-Jacobian elements (showing up to 5 of 3):
+    element 7: centroid (0.5, 0.5, 0.5), min Jacobian = -1.000e-03
 ```
 
-By default the conversion continues despite warnings (the Exodus file is still written; the user is informed). Two flags adjust this:
+By default the output file is still written regardless of warnings. Three flags adjust this:
 
-- `--strict-jacobians` upgrades any non-positive Jacobian to a fatal error (exit code 1). Useful in CI / scripted workflows where a bad mesh should stop the pipeline rather than propagate downstream.
+- `--remove-distorted` removes elements with non-positive Jacobians before writing, printing a count of those dropped. Useful when a small number of distorted cells near faults would otherwise cause solver failures.
+- `--strict-jacobians` upgrades any non-positive Jacobian to a fatal error (exit code 1). Useful in CI / scripted workflows where a bad mesh should stop the pipeline.
 - `--no-check-jacobians` skips the check entirely (a small time saving on very large grids, but disables a useful safety net).
 
+**Relationship to `--pinch`.** `--remove-distorted` and `--pinch` are complementary rather than interchangeable. `--pinch` detects cells where any two corners are within `--pinch-tol` of each other — it catches *near*-coincident corners (e.g. a 0.5 m thick cell in a grid measured in metres) even when the Jacobian is still technically positive. `--remove-distorted` catches cells whose Jacobian has already reached zero or gone negative, which only happens once corners are *exactly* coincident or the cell has become inverted. In practice, running both is the safest option for grids with thin reservoir layers near faults: `--pinch` handles the near-zero cells that `--remove-distorted` would miss, and `--remove-distorted` catches any remaining inverted cells.
+
 `em2ex` attempts to guess the reservoir model format from the file extension (see supported formats below). If the reservoir model has a non-standard file extension, the user can force
 `em2ex` to read the correct format using the `--filetype` commandline option.
 

em2ex.py

@@ -117,7 +117,7 @@ def get_parser():
     parser.add_argument('--flip', dest = 'flip_z', action = 'store_true', help = 'Flip the sign of the Z coordinates')
     parser.add_argument('--translate', nargs = 2, dest = 'translate', type = float, help = 'Translate the (x, y) coordinates by this amount')
     parser.add_argument('--mapaxes', dest = 'use_mapaxes', action = 'store_true', help = 'Use the MAPAXES coordinates for an Eclipse file')
-    parser.add_argument('--pinch', default = True, dest = 'no_pinch', action = 'store_true', help = 'Remove pinched elements')
+    parser.add_argument('--pinch', default = False, dest = 'no_pinch', action = 'store_true', help = 'Remove pinched elements (coincident corners within --pinch-tol). Off by default; faithful conversion is produced without this flag.')
     parser.add_argument('--pinch-tol', default = 1e-3, dest = 'pinch_tol', type = float, help = 'Tolerance for coincident corners when removing pinched elements (default: 1e-3)')
     parser.add_argument('--refine-xy', nargs = 2, dest = 'refine_xy', type = _positive_int,
         metavar = ('RX', 'RY'),
@@ -142,6 +142,8 @@ def get_parser():
         help = 'Skip the per-element Jacobian sanity check. By default em2ex computes the Jacobian at all 8 corners of every HEX8 element and warns if any are non-positive (degenerate or inverted), since those elements would be rejected by most FEM solvers.')
     parser.add_argument('--strict-jacobians', dest = 'strict_jacobians', action = 'store_true',
         help = 'Treat any non-positive element Jacobian as a fatal error and exit non-zero. By default such elements only produce a warning. Useful for CI / scripted workflows.')
+    parser.add_argument('--remove-distorted', dest = 'remove_distorted', action = 'store_true',
+        help = 'Remove elements with non-positive Jacobians (degenerate or inverted) from the output mesh, reporting a count of those removed. By default such elements are kept and only a warning is printed.')
     return parser
 
 def main():

readers/eclipse.py

@@ -367,6 +367,14 @@ def parseEclipse(f, args):
             xvec = xvec / np.sqrt(xvec[0]**2 + xvec[1]**2)
             yvec = yvec / np.sqrt(yvec[0]**2 + yvec[1]**2)
 
+            # If the MAPAXES transform is orientation-reversing (det < 0), negate xvec so
+            # the resulting mesh has positive Jacobians instead of all-negative ones.
+            det = xvec[0] * yvec[1] - xvec[1] * yvec[0]
+            if det < 0:
+                print("MAPAXES transform is orientation-reversing (det={:.6g}); "
+                      "negating x-axis direction to restore positive Jacobians".format(det))
+                xvec = -xvec
+
             # Transform top and bottom pillar x,y: world = origin + local_x*xhat + local_y*yhat
             for xi, yi in [(0, 1), (3, 4)]:
                 xdata = coord[:,:,xi].copy()
@@ -404,12 +412,40 @@ def parseEclipse(f, args):
         # All elements are active
         active_elements = np.ones((nz, ny, nx), dtype = int)
 
-    # If pinched or distorted elements are removed, these elements will also be set to
-    # be inactive.
-    # Check for pinched out elements (where the vertical distance between nodes is less than tol)
+    # Check for pinched elements (coincident corners within pinch_tol).
+    # Always detect so a count can be reported; only remove when --pinch is passed.
+    distorted = distortedElem(elemcornx, elemcorny, elemcornz, args.pinch_tol).reshape(nz, ny, nx)
+    n_pinched = int((distorted & (active_elements > 0)).sum())
     if args.no_pinch:
-        distorted = distortedElem(elemcornx, elemcorny, elemcornz, args.pinch_tol).reshape(nz, ny, nx)
+        if n_pinched > 0:
+            print('{} pinched element(s) removed'.format(n_pinched))
         active_elements[distorted] = 0
+    elif n_pinched > 0:
+        print('Note: {} element(s) have coincident corners (pinched) and will be '
+              'invalid in FEM solvers. Use --pinch to remove them.'.format(n_pinched))
+
+    # Optionally remove elements with non-positive Jacobian (degenerate or inverted).
+    # Off by default (faithful conversion); enabled with --remove-distorted.
+    if getattr(args, 'remove_distorted', False):
+        from readers.reader_utils import nonPositiveJacobianElems
+        # When z is flipped, the top/bottom corner swap applied later to elemNodes
+        # (corners [4,5,6,7,0,1,2,3]) restores positive Jacobians. Apply the same
+        # permutation here so the check reflects the final assembled element orientation.
+        cx, cy, cz = elemcornx, elemcorny, elemcornz
+        if args.flip_z:
+            perm = [4, 5, 6, 7, 0, 1, 2, 3]
+            cx, cy, cz = cx[:, perm], cy[:, perm], cz[:, perm]
+        bad_jac = nonPositiveJacobianElems(cx, cy, cz).reshape(nz, ny, nx)
+        newly_removed = bad_jac & (active_elements > 0)
+        n_removed = int(newly_removed.sum())
+        if n_removed > 0:
+            n_active = int((active_elements > 0).sum())
+            if n_removed == n_active:
+                print('Warning: --remove-distorted would remove all {} active element(s). '
+                      'Skipping removal — check grid orientation (try --flip).'.format(n_active))
+            else:
+                print('{} element(s) with zero or negative element Jacobian removed'.format(n_removed))
+                active_elements[bad_jac] = 0
 
     # The number of active elements is
     num_active_elements = np.count_nonzero(active_elements)

readers/reader_utils.py

@@ -154,6 +154,12 @@ def checkElementJacobians(model, strict=False):
     print('Element Jacobian check: {} negative, {} zero, {} OK (out of {})'.format(
         num_neg, num_zero, num_ok, model.numElems))
 
+    if num_ok == 0:
+        print('  All elements have non-positive Jacobians. Possible causes:')
+        print('    - z-up coordinate convention: try --flip to invert the z-axis')
+        print('    - Orientation-reversing coordinate system (e.g. MAPAXES handedness)')
+        print('  Use --remove-distorted to remove these elements and proceed anyway.')
+
     # Map elemNodes row -> Exodus element ID. elemIds[k, j, i] holds the
     # 1-based Exodus ID (0 for inactive); the non-zero values in (k, j, i)
     # flat order align with elemNodes' row ordering.
@@ -182,6 +188,23 @@ def checkElementJacobians(model, strict=False):
     return False
 
 
+def nonPositiveJacobianElems(elemcornx, elemcorny, elemcornz):
+    ''' Return a boolean array (numelems,) that is True for elements whose
+    minimum per-corner Jacobian is <= 0 (degenerate or inverted).
+
+    Inputs are per-element corner coordinate arrays of shape (numelems, 8),
+    ordered in the HEX8 element corner layout (matching elemCornerCoords output).
+    '''
+    P = np.stack([elemcornx, elemcorny, elemcornz], axis=-1)  # (N, 8, 3)
+    jacobians = np.empty((P.shape[0], 8))
+    for c, ((xi_a, xi_b), (eta_a, eta_b), (zeta_a, zeta_b)) in enumerate(_HEX8_JAC_EDGES):
+        e_xi   = P[:, xi_a]   - P[:, xi_b]
+        e_eta  = P[:, eta_a]  - P[:, eta_b]
+        e_zeta = P[:, zeta_a] - P[:, zeta_b]
+        jacobians[:, c] = np.einsum('ij,ij->i', e_xi, np.cross(e_eta, e_zeta))
+    return jacobians.min(axis=1) <= 0
+
+
 def nonZeroValues(arr):
     ''' Utility to determine first non-zero values in the top plane of an array '''
 

run_tests.py

@@ -66,6 +66,13 @@ def test_em2ex(key, use_official_api, exodiff):
 
             assert tests[key]['expected_error'] in str(excinfo.value)
 
+    # If the test type is output, check that expected text appears in stdout
+    elif tests[key]['type'] == 'output':
+        if 'expected_output' not in tests[key].keys():
+            pytest.skip(key + ': Skipped as expected_output not specified')
+        else:
+            check_output(key)
+
     else:
         # Skip unknown test type
         pytest.skip(key + ': Skipped as unknown test type')
@@ -125,3 +132,24 @@ def expected_error(key):
     raise Em2exException(result.stdout + result.stderr)
 
     return
+
+def check_output(key):
+    ''' Check that expected text appears in em2ex stdout on a successful run '''
+
+    filepath = tests[key]['filepath']
+    filename = tests[key]['filename']
+    testfilename = os.path.join(filepath, filename)
+
+    command = ['./em2ex.py', '-f']
+    if 'cli_args' in tests[key].keys():
+        command.extend(tests[key]['cli_args'].split())
+    command.append(testfilename)
+
+    result = subprocess.run(command, capture_output=True, text=True)
+    combined = result.stdout + result.stderr
+
+    assert tests[key]['expected_output'] in combined, \
+        '{}: expected output not found.\nExpected: {}\nGot: {}'.format(
+            key, tests[key]['expected_output'], combined)
+
+    return

test/eclipse/tests

@@ -38,6 +38,8 @@ simple_cube_translated:
 simple_cube_mapaxes:
   filename: simple_cube_mapaxes.grdecl
   type: exodiff
+  # The MAPAXES spec (1.0 0.0 0.0 0.0 0.0 1.0) is orientation-reversing (det=-1).
+  # The reader detects this and negates xvec so the output has positive Jacobians.
   cli_args: --mapaxes
   gold: simple_cube_mapaxes.e
 
@@ -151,15 +153,22 @@ simple_cube_quoted_unit:
   type: exodiff
   gold: simple_cube.e
 
-# Z-up convention file (pillar tops at z=1, bots at z=0). em2ex's default
-# corner ordering treats this as inverted; --strict-jacobians should error
-# out rather than write a mesh
+# Z-up convention file (pillar tops at z=1, bots at z=0). em2ex produces a
+# faithful conversion but the Jacobian check sees all-negative elements;
+# --strict-jacobians upgrades that to a fatal error.
 strict_jacobians_inverted:
   filename: simple_cube_zup.grdecl
   type: exception
   cli_args: --strict-jacobians
   expected_error: exiting due to invalid Jacobians
 
+# Same z-up file without --strict-jacobians: conversion succeeds but the
+# all-negative hint (suggesting --flip) should be printed.
+all_negative_jacobians_hint:
+  filename: simple_cube_zup.grdecl
+  type: output
+  expected_output: All elements have non-positive Jacobians
+
 # Driven entirely by a YAML config file (extract + refine settings). The
 # output must match the same gold as the CLI-driven simple_cube_extract_refine
 # test, proving the config loader applies values correctly.