This may be a corrupt mesh error

[Abby123455r]

Problem description

i made a 3d geometry on freecad. and then mesh it on gmsh..after that i import it in pygimli for simulation. but when i try to invert the simulated data it gives me following error:
GIMLI::Cell* GIMLI::Mesh::findCell(const GIMLI::RVector3&, size_t&, bool) const no cells or boundaries for this node. This may be a corrupt mesh

Your environment

OS : Windows
CPU(s) : 32
Environment : Jupyter

Operating system: Windows,
Python version: 3.10,
pyGIMLi version: 1.5.0
Way of installation: Conda package

Steps to reproduce

import pygimli.meshtools as mt
import pygimli as pg
import pybert as pb
import gmsh
from pygimli.physics import ert
gmsh.initialize()
msh_nm = 'cone.msh'
gmsh.open(msh_nm)
mesh_mod = mt.readGmsh(msh_nm, verbose=True)

Reading bucket cone.msh...

Nodes: 3295
Entries: 18299
Points: 32
Lines: 0
Triangles: 2385
Quads: 0
Tetrahedra: 15882

Creating mesh object...

Dimension: 3-D
Boundary types: 2 (-2, -1)
Regions: 1 (1,)
Marked nodes: 32 (99,)

Done.

Mesh: Nodes: 3295 Cells: 15882 Boundaries: 33154

pg.show(mesh_mod, markers=True, showMesh=True)
bert_file = '01_1.udf'
BERT = pg.load(bert_file)

mgr = ert.ERTManager('01_1.udf')
# Define resistivity mapping (rhomap)
rhomap = [
    [1, 300.0]]
data_mod = ert.simulate(mesh_mod , res=rhomap  scheme=BERT, noiseLevel=1,
                    noiseAbs=1e-6, seed=1337 )
data_mod.save('practice.dat' )
mgr = ert.ERTManager('practice.dat')
inv = mgr.invert(lam=5, verbose=True)
np.testing.assert_approx_equal(mgr.inv.chi2(), 0.7, significant=1)

Paste your script output here.

15/05/24 - 12:12:42 - pyGIMLi - INFO - Found 3 regions.
15/05/24 - 12:12:42 - pyGIMLi - INFO - (ERTModelling) Region with smallest marker (0) set to background.
15/05/24 - 12:12:42 - pyGIMLi - INFO - Found 3 regions.
15/05/24 - 12:12:42 - pyGIMLi - INFO - (ERTModelling) Region with smallest marker (0) set to background.
15/05/24 - 12:12:42 - pyGIMLi - INFO - Creating forward mesh from region infos.
15/05/24 - 12:12:42 - pyGIMLi - INFO - Creating refined mesh (H2) to solve forward task.
15/05/24 - 12:12:42 - pyGIMLi - INFO - Mesh for forward task: Mesh: Nodes: 2093 Cells: 3824 Boundaries: 3042
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[112], line 1
----> 1 inv = mgr.invert(lam=5, verbose=True)
      2 np.testing.assert_approx_equal(mgr.inv.chi2(), 0.7, significant=1)

File C:\WINDOWS\system32\D\pg\Lib\site-packages\pygimli\frameworks\methodManager.py:825, in MeshMethodManager.invert(self, data, mesh, startModel, **kwargs)
    822 dataVals = self._ensureData(self.fop.data)
    823 errorVals = self._ensureError(self.fop.data, dataVals)
--> 825 if self.fop.mesh() is None:
    826     pg.critical('Please provide a mesh')
    828 kwargs['startModel'] = startModel

File C:\WINDOWS\system32\D\pg\Lib\site-packages\pygimli\frameworks\modelling.py:709, in MeshModelling.mesh(self)
    706 self._applyRegionProperties()
    708 if self._regionManagerInUse and self._regionChanged is True:
--> 709     self.createFwdMesh_()
    711 return super().mesh()

File C:\WINDOWS\system32\D\pg\Lib\site-packages\pygimli\frameworks\modelling.py:695, in MeshModelling.createFwdMesh_(self)
    685     pg.verbose("\tRegion: {0}, Parameter: {1}, PD: {2},"
    686                " Single: {3}, Background: {4}, Fixed: {5}".format(
    687                    iId,
   (...)
    691                    self.regionManager().region(iId).isBackground(),
    692                    self.regionManager().region(iId).fixValue()))
    694 m = self.createRefinedFwdMesh(m)
--> 695 self.setMeshPost(m)
    697 self._regionChanged = False
    698 super().setMesh(m, ignoreRegionManager=True)

File C:\WINDOWS\system32\D\pg\Lib\site-packages\pygimli\physics\ert\ertModelling.py:211, in ERTModelling.setMeshPost(self, mesh)
    209 def setMeshPost(self, mesh):
    210     """Set mesh (at a later stage)."""
--> 211     self._core.setMesh(mesh, ignoreRegionManager=True)

RuntimeError: ./core/src/mesh.cpp:928		GIMLI::Cell* GIMLI::Mesh::findCell(const GIMLI::RVector3&, size_t&, bool) const  no cells or boundaries for this node. This may be a corrupt mesh

If possible, please add one or more labels to your issue, e.g. if you expect that your issue is rather a question than a problem with the code, please add the label "question".

[andieie]

It doesn't look like you are passing a mesh to the inversion from your code so it does one for you. I haven't properly understood if you would like to use mesh_mod for your inversion. It seems you have to create another mesh for the inversion.

[Abby123455r]

i want to use same mesh for both modeling and inversion.is it not possible to use same mesh?

[halbmy]

Of couse you can, but using the same mesh for synthetic and inverse modelling is considered a so-called inverse crime as you already provide information that you normally don't have. At least you should reset all cell markers by mesh['marker']=0.

[Abby123455r]

thank you it is working.but it is showing large chi2 and negitive dphi:
inv.iter 0 ... chi² = 12550.59

inv.iter 1 ... chi² = 3388.10 (dPhi = 72.10%) lam: 20.0

inv.iter 2 ... chi² = 2880.36 (dPhi = 15.17%) lam: 20.0

inv.iter 3 ... chi² = 2713.54 (dPhi = 6.22%) lam: 20.0

inv.iter 4 ... chi² = 3033.46 (dPhi = -11.19%) lam: 20.0
################################################################################

Abort criterion reached: dPhi = -11.19 (< 2.0%)

[halbmy]

Hard to give any recommendations without knowing details. I guess the mesh is just inappropriate with this low node number, maybe not even containing the electrodes.

[Abby123455r]

import pygimli.meshtools as mt
import pygimli as pg
import pybert as pb
import gmsh
from pygimli.physics import ert
gmsh.initialize()
msh_nm = 'bucket cone.msh'
mesh_mod = mt.readGmsh(msh_nm, verbose=True)

output:Reading bucket cone.msh...
Nodes: 8455
Entries: 49328
Points: 32
Lines: 0
Triangles: 5586
Quads: 0
Tetrahedra: 43710

Creating mesh object...

Dimension: 3-D
Boundary types: 2 (-2, -1)
Regions: 1 (1,)
Marked nodes: 32 (99,)

Done.
Mesh: Nodes: 8455 Cells: 43710 Boundaries: 90213.

pg.show(mesh_mod, markers=True, showMesh=True)
bert_file = '01_1.udf'

# Import data using PyGimli
BERT = pg.load(bert_file)

# Display information about the imported data
print(BERT)
Data: Sensors: 32 data: 675, nonzero entries: ['a', 'b', 'm', 'n', 'valid']
mgr = ert.ERTManager('01_1.udf')
# Define resistivity mapping (rhomap)
rhomap = [
    [1, 300.0]]
data_mod = ert . simulate (mesh_mod , res=rhomap , scheme=BERT, noiseLevel=1,
                    noiseAbs=1e-6, seed=1337 )
data_mod . save ( 'practice.dat' )
data = pg.load("practice.dat")
print(data)
mgr = ert.ERTManager('practice.dat')
 mesh_mod['marker']=0
mgr.invert(verbose=True,mesh=mesh_mod,lam=10)

#out put of inversion:

23/05/24 - 12:20:52 - pyGIMLi - INFO - Found 1 regions.
23/05/24 - 12:20:52 - pyGIMLi - INFO - Creating forward mesh from region infos.
23/05/24 - 12:20:56 - pyGIMLi - INFO - Creating refined mesh (H2) to solve forward task.
23/05/24 - 12:21:07 - pyGIMLi - INFO - Mesh for forward task: Mesh: Nodes: 63412 Cells: 349680 Boundaries: 360852
23/05/24 - 12:21:15 - pyGIMLi - INFO - Use median(data values)=255.563433298874
23/05/24 - 12:21:15 - pyGIMLi - INFO - Created startmodel from forward operator: 43710, min/max=255.563433/255.563433
23/05/24 - 12:21:15 - pyGIMLi - INFO - Starting inversion.
fop: <pygimli.physics.ert.ertModelling.ERTModelling object at 0x00000193DFA5D2C0>
Data transformation: <pgcore._pygimli_.RTransLogLU object at 0x00000193FAFE8F40>
Model transformation: <pgcore._pygimli_.RTransLog object at 0x00000193DFA5D310>
min/max (data): 148/572
min/max (error): 1%/1%
min/max (start model): 256/256
--------------------------------------------------------------------------------
inv.iter 0 ... chi² =  256.64
--------------------------------------------------------------------------------
inv.iter 1 ... chi² =  243.72 (dPhi = 5.04%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 2 ... chi² =  236.36 (dPhi = 3.02%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 3 ... chi² =  230.17 (dPhi = 2.62%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² =  223.60 (dPhi = 2.85%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 5 ... chi² =  212.04 (dPhi = 5.17%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 6 ... chi² =  206.51 (dPhi = 2.61%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 7 ... chi² =  201.27 (dPhi = 2.53%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 8 ... chi² =  196.75 (dPhi = 2.25%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 9 ... chi² =  192.43 (dPhi = 2.20%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 10 ... chi² =  188.30 (dPhi = 2.15%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 11 ... chi² =  183.95 (dPhi = 2.31%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 12 ... chi² =  180.17 (dPhi = 2.05%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 13 ... chi² =  176.55 (dPhi = 2.01%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 14 ... chi² =  172.73 (dPhi = 2.16%) lam: 10.0
--------------------------------------------------------------------------------
inv.iter 15 ... chi² =  169.40 (dPhi = 1.93%) lam: 10.0
################################################################################
#                Abort criterion reached: dPhi = 1.93 (< 2.0%)                 #
################################################################################
43710 [255.56343329887392,...,255.69008917555843]

The value of chi2 is too high

[halbmy]

Nobody can really help without seeing the mesh (maybe include a figure) and without seeing the data.

First point for high chi-square is having a look at the data. So please compute geometric factors numerically and show the apparent resistivity.

If the data look good, I suggest using total field computation with quadratic shape functions to ensure the quality of the forward operator by something like:

mgr = ERTManager(data, sr=False)
mgr.setMesh(mesh)
mgr.fop._core.createRefinedForwardMesh(True, True)

At any rate, we clearly need an example of ERT on a closed 3D geometry and probably some options to the ERT manager controlling the refinement.