Description
Problem description
This issue is an extension of Gimli Issue #589, specifically related to creating a 2D quadrilateral mesh with a geometry surface that has an upward concave terrain.
Problem Summary:
I am trying to create an inversion mesh using quadrilateral grids and append a triangular boundary using the mt.appendTriangleBoundary method. When the terrain is upward concave, the method fails to append the triangular boundary. However, when the terrain is downward concave, the method works as expected, and the triangular boundary is successfully added.
Your environment
OS : Windows (10 10.0.22631 SP0 Multiprocessor Free)
CPU(s) : 16
Machine : AMD64
Architecture : 64bit
RAM : 23.6 GiB
Environment : Jupyter
Python 3.11.10 | packaged by conda-forge | (main, Oct 16 2024, 01:17:14)
[MSC v.1941 64 bit (AMD64)]
pygimli : 1.5.3
pgcore : 1.5.0
numpy : 1.26.4
matplotlib : 3.9.2
scipy : 1.14.1
tqdm : 4.67.0
IPython : 8.29.0
pyvista : 0.44.1
Steps to reproduce
Case 1: Upward Concave Terrain (Fails)
import numpy as np
import pygimli as pg
from pygimli.physics import ert # the module
import pygimli.meshtools as mt
topo = [[x, 1.0 + np.cos(2 * np.pi * 1/30 * x + 2 * np.pi)] for x in range(31)]
topo = np.array(topo)
depth = 10
yDevide = (1.0 - np.logspace(np.log10(1.0), np.log10(1 + depth), 31))[::-1]
grid = pg.createGrid(x=topo[:, 0], y=yDevide, worldBoundaryMarker=True)
pg.show(grid, boundaryMarkers=True)
mx = pg.x(grid)
meshtopo = np.interp(mx, topo[:, 0], topo[:, 1])
for i, n in enumerate(grid.nodes()):
n.setPos([n.x(), n.y() + meshtopo[i]])
mesh = mt.appendTriangleBoundary(grid, xbound=10, ybound=10, marker=2, isSubSurface=True)
ax, _ = pg.show(mesh, boundaryMarkers=True, showMesh=True, markers=True)Expected behavior
Case 2: Downward Concave Terrain (Succeeds)
topo = [[x, 1.0 + np.cos(2 * np.pi * 1/30 * x + 1 * np.pi)] for x in range(31)]
topo = np.array(topo)
depth = 10
yDevide = (1.0 - np.logspace(np.log10(1.0), np.log10(1 + depth), 31))[::-1]
grid = pg.createGrid(x=topo[:, 0], y=yDevide, worldBoundaryMarker=True)
pg.show(grid, boundaryMarkers=True)
mx = pg.x(grid)
meshtopo = np.interp(mx, topo[:, 0], topo[:, 1])
for i, n in enumerate(grid.nodes()):
n.setPos([n.x(), n.y() + meshtopo[i]])
mesh = mt.appendTriangleBoundary(grid, xbound=10, ybound=10, marker=2, isSubSurface=True)
ax, _ = pg.show(mesh, boundaryMarkers=True, showMesh=True, markers=True)