search

gimli-org / gimli

Geophysical Inversion and Modeling Library :earth_africa:
https://www.pygimli.org
Other
384 stars 137 forks source link

Region-wise regularization with changing water level #764

Open HenryWHR opened 3 months ago

HenryWHR commented 3 months ago

Problem description

Based on the example in 'Region-wise regularization', I want to fix the base mesh when I have different water levels. I tried with the example, but it seems that, if I have different water levels, the base mesh can change. How can I fix the base grid so that I can do direct difference analyses when the water level is different?

Your environment

        OS : Windows
        CPU(s) : 14
        Machine : AMD64
        Architecture : 64bit
        RAM : 31.5 GiB
        Environment : IPython

Python 3.11.9 | packaged by conda-forge | (main, Apr 19 2024, 18:27:10) [MSC v.1938 64 bit (AMD64)]

       pygimli : 1.5.1
        pgcore : 1.5.0
         numpy : 1.26.4
        matplotlib : 3.8.4
         scipy : 1.13.0
          tqdm : 4.66.4
       IPython : 8.20.0
       pyvista : 0.43.8

Steps to reproduce

import pygimli as pg
from pygimli.physics import ert
import pygimli.meshtools as mt

data = ert.load('lake.ohm')
plc = mt.createParaMeshPLC(data, paraDepth=20, boundary=1)

plc.createEdge(plc.node(11), plc.node(99), marker=-1)

# lower water level
# plc.createEdge(plc.node(10), plc.node(100), marker=-1)

plc.addRegionMarker([50, -0.1], marker=3)
ax, _ = pg.show(plc, markers=True)
mesh = mt.createMesh(plc, quality=34.5)
base = mesh.createSubMesh(mesh.cells(mesh.cellMarkers() == 2))
ax, cb = pg.show(base,showMesh=True)
print(base)

...
HenryWHR commented 3 months ago

By using merge2Meshes, It seems that I can fix the base mesh with varying water levels. If there are better options to achieve this, please let me know! Thank you!

One question, there are fewer Boundaries for the water body after the merge of two meshes. The Cells are the same. What does this mean?

from pygimli.physics import ert
import pygimli.meshtools as mt
import numpy as np

data = ert.load('lake.ohm')
plc = mt.createParaMeshPLC(data, paraDepth=20, boundary=1)

mesh_base = mt.createMesh(plc, quality=34.5)
mesh_base_inv = mesh_base.createSubMesh(mesh_base.cells(mesh_base.cellMarkers() == 2))
sensors = np.array(data.sensors())
sensors = np.delete(sensors, 1, axis=1)

waterbottom1 = sensors[1:47,:]
waterpoly1 = mt.createPolygon(waterbottom1, isClosed=True, marker=3)
mesh_water1 = mt.createMesh(waterpoly1, quality=34.5)
mesh_merge1 = mt.merge2Meshes(mesh_base,mesh_water1)
mesh_merge1_inv = mesh_merge1.createSubMesh(mesh_merge1.cells(mesh_merge1.cellMarkers() == 2))
mesh_merge1_water = mesh_merge1.createSubMesh(mesh_merge1.cells(mesh_merge1.cellMarkers() == 3))

waterbottom2 = sensors[3:45,:]
waterpoly2 = mt.createPolygon(waterbottom2, isClosed=True, marker=3)
mesh_water2 = mt.createMesh(waterpoly2, quality=34.5)
mesh_merge2 = mt.merge2Meshes(mesh_base,mesh_water2)
mesh_merge2_inv = mesh_merge2.createSubMesh(mesh_merge2.cells(mesh_merge2.cellMarkers() == 2))
mesh_merge2_water = mesh_merge2.createSubMesh(mesh_merge2.cells(mesh_merge2.cellMarkers() == 3))

print(mesh_base_inv)
print(mesh_merge1_inv)
print(mesh_merge2_inv)
print(mesh_water1)
print(mesh_merge1_water)
print(mesh_water2)
print(mesh_merge2_water)

Mesh: Nodes: 672 Cells: 1202 Boundaries: 140 Mesh: Nodes: 672 Cells: 1202 Boundaries: 140 Mesh: Nodes: 672 Cells: 1202 Boundaries: 140 Mesh: Nodes: 163 Cells: 185 Boundaries: 347 Mesh: Nodes: 163 Cells: 185 Boundaries: 139 Mesh: Nodes: 126 Cells: 129 Boundaries: 254 Mesh: Nodes: 126 Cells: 129 Boundaries: 121