Add TetGen output capture and mesh quality warning system

src/vasp/preprocessing/vmtkmeshgeneratorfsi.py

@@ -24,6 +24,8 @@
 from vmtk import vtkvmtk, vmtkscripts, pypes
 import sys
 import numpy as np
+import os
+import re
 
 
 class vmtkMeshGeneratorFsi(pypes.pypeScript):
@@ -134,6 +136,67 @@ def __init__(self):
             ['RemeshedSurface', 'remeshedsurface', 'vtkPolyData', 1, '', 'the output surface', 'vmtksurfacewriter'],
         ])
 
+    def _run_and_analyze_tetgen(self, tetgen):
+        """
+        Execute TetGen and analyze the printed output for mesh quality.
+        """
+        old_stdout = os.dup(1)
+        old_stderr = os.dup(2)
+        r_stdout, w_stdout = os.pipe()
+        r_stderr, w_stderr = os.pipe()
+        os.dup2(w_stdout, 1)
+        os.dup2(w_stderr, 2)
+        os.close(w_stdout)
+        os.close(w_stderr)
+
+        try:
+            tetgen.Execute()
+            os.close(1)
+            os.close(2)
+            output = os.read(r_stdout, 100000).decode()
+            _ = os.read(r_stderr, 100000)  # discard stderr output for now
+        finally:
+            os.dup2(old_stdout, 1)
+            os.dup2(old_stderr, 2)
+            os.close(old_stdout)
+            os.close(old_stderr)
+
+        if output:
+            self._parse_tetgen_output(output)
+
+    def _parse_tetgen_output(self, text):
+        """
+        Parse TetGen output and print warnings for low-quality elements.
+        """
+        def extract(pattern):
+            match = re.search(pattern, text)
+            return float(match.group(1)) if match else None
+
+        smallest_volume = extract(r"Smallest volume:\s+([\d.e+-]+)")
+        largest_volume = extract(r"Largest volume:\s+([\d.e+-]+)")
+        _ = extract(r"Smallest aspect ratio:\s+([\d.e+-]+)")  # not currently used
+        largest_aspect = extract(r"Largest aspect ratio:\s+([\d.e+-]+)")
+        smallest_dihedral = extract(r"Smallest dihedral:\s+([\d.e+-]+)")
+        largest_dihedral = extract(r"Largest dihedral:\s+([\d.e+-]+)")
+        smallest_facangle = extract(r"Smallest facangle:\s+([\d.e+-]+)")
+        largest_facangle = extract(r"Largest facangle:\s+([\d.e+-]+)")
+
+        self.PrintLog("TetGen mesh quality analysis:")
+        if largest_aspect and largest_aspect > 10:
+            self.PrintLog(f"  Warning: High aspect ratio (max {largest_aspect:.2f})")
+        if smallest_dihedral and smallest_dihedral < 10:
+            self.PrintLog(f"  Warning: Small dihedral angle (min {smallest_dihedral:.2f}°)")
+        if largest_dihedral and largest_dihedral > 160:
+            self.PrintLog(f"  Warning: Large dihedral angle (max {largest_dihedral:.2f}°)")
+        if smallest_facangle and smallest_facangle < 20:
+            self.PrintLog(f"  Warning: Small face angle (min {smallest_facangle:.2f}°)")
+        if largest_facangle and largest_facangle > 140:
+            self.PrintLog(f"  Warning: Large face angle (max {largest_facangle:.2f}°)")
+        if smallest_volume and largest_volume:
+            ratio = smallest_volume / largest_volume
+            if ratio < 0.01:
+                self.PrintLog(f"  Warning: Smallest volume is <1% of largest volume (ratio {ratio:.4f})")
+
     def Execute(self):
 
         if self.Surface is None:
@@ -368,7 +431,9 @@ def Execute(self):
             tetgen.OutputSurfaceElements = 1
             tetgen.OutputVolumeElements = 1
             tetgen.RegionAttrib = 0
-            tetgen.Execute()
+            tetgen.Verbose = 1
+
+            self._run_and_analyze_tetgen(tetgen)
 
             if tetgen.Mesh.GetNumberOfCells() == 0 and surfaceToMesh.Mesh.GetNumberOfCells() > 0:
                 self.PrintError('Running TetGen failed. Try to re-mesh.')
@@ -453,8 +518,9 @@ def VisitNeighbors(i, cellEntityId):
             tetgen.RemoveSliver = 1
             tetgen.OutputSurfaceElements = 1
             tetgen.OutputVolumeElements = 1
+            tetgen.Verbose = 1
 
-            tetgen.Execute()
+            self._run_and_analyze_tetgen(tetgen)
 
             self.Mesh = tetgen.Mesh