petro.resistivityArchie interpolation returns 'corrupt mesh' error

[Coastal0]

I'm certain this was working yesterday, and now this morning - it does not!!!!

Why would it say my mesh is now corrupt?!

I've tried re-creating the mesh without additional arguments and still get the error...

resBulk = petro.resistivityArchie(rFluid=(1. / dataInterpArray), porosity=0.3, m=1.8, mesh=datamesh_pointcloud,meshI = meshERT,fill=True)

meshERT is created from a set of sensors located on the topography.

e.g.

meshERT = pg.meshtools.createParaMesh(ertScheme, quality=33, paraMaxCellSize=100, smooth=[1, 3])

Resulting in the error below;

resBulk = petro.resistivityArchie(rFluid=(1. / dataInterpArray), porosity=0.3, m=1.8, mesh=datamesh_pointcloud,meshI = meshERT,fill=True)
Traceback (most recent call last):

  File "<ipython-input-9-42e74e8f87e9>", line 1, in <module>
    resBulk = petro.resistivityArchie(rFluid=(1. / dataInterpArray), porosity=0.3, m=1.8, mesh=datamesh_pointcloud,meshI = meshERT,fill=True)

  File "F:\WinPython-64bit-3.6.2.0Qt5\python-3.6.2.amd64\Lib\pygimli\physics\petro\resistivity.py", line 97, in resistivityArchie
    pg.interpolate(mesh, r, meshI.cellCenters(), rI)

  File "F:\WinPython-64bit-3.6.2.0Qt5\python-3.6.2.amd64\Lib\pygimli\meshtools\mapping.py", line 405, in interpolate
    return pg.core._pygimli_.interpolate(*args, **kwargs)

RuntimeError: C:/msys64/home/Guenther.T/py36/gimli/gimli/src/mesh.cpp: 524		GIMLI::Cell* GIMLI::Mesh::findCell(const RVector3&, size_t&, bool) const  no cells for this node. This is a corrupt mesh

Winpython-64bit-3.6.2.0Qt5, Spyder 3.2.4, pyGimli-1.0.3-cp36-cp36m-win_amd64.whl

[carsten-forty2]

If you don't updated GIMLI .. I assume your input has been changed .. such spooky behavior seems a little unlikely.

The error message says one of the meshes have free nodes (nodes that are not related to the mesh elements).
This usual happens if the input of the mesh creation is invalid. e.g., duplicated input nodes.
Please check ertScheme for duplicated sensor position and take a view to booth meshes, datamesh_pointcloud and meshERT. If the error persist please provide all data that we can reproduce the problem.

[Coastal0]

You're quite right, something must have changed. I just can't see what exactly...

There doesn't seem to be any duplicate sensors in the ertScheme.
(On that note - how do I inspect the ertScheme.sensorPositions() (or indeed, any RVector3) data?)

All of the files are below;
ConcaveBoundary.dat.txt
Pointcloud.pnt.txt

Alexs_Main_Workbook - Copy.py.txt

Edit;

What is the best way to check if free nodes exist?!

I can't understand why free nodes should exist on the parameter mesh, as it was generated with the sensor list.
Everything should have been triangulated!!! 😖

It's quite obviously something to do with the interpolation,
For example, the same error message occurs with this command;
resBulk = pg.interpolate(datamesh_pointcloud_border,rFluid,meshERT.cellCenter())

but I still can't work out why there would be free nodes

[carsten-forty2]

Maybe there are duplicated points in your point cloud? You can count through you mesh nodes e.g.

for n in mesh.nodes():
    if  len(n.cellSet()) == 0:
        print(n, n.pos(), " have no cells!")

btw. It would be easier for us to check stuff, if you provide either all meshes or data directly. e.g. to check the command resBulk = pg.interpolate(datamesh_pointcloud_border,rFluid,meshERT.cellCenter()) we need datamesh_pointcloud_border, rFluid and meshERT and/or the condensed script that create them. You are familiar to you own scripts .. for us, it takes always some time to puzzle it together.

[carsten-forty2]

In your concaveBoundary.dat the first and last positions are equal ..

As this seems to happen regularly .. I think we will drop a Warning about in mt.createPolygon().

If you want to close the polygon set isClosed=True to the arguments.
https://www.pygimli.org/pygimliapi/_generated/pygimli.meshtools.html#pygimli.meshtools.createPolygon

[Coastal0]

Good point - sorry about that.

My scripts probably quite poorly written too. It's all a bit intertwined at the moment.
The zip file below has the datamesh_pointcloud_border, rFluid, and meshERT saved in Ascii format from within pygimli.
Hopefully it's easy to just load them in...?

pyGimli_files.zip

I removed those equal points in the hull - I think it may have solved some instability issues I was having when making meshes, (:yay:), but hasn't immediately solved the other issues yet.

[carsten-forty2]

why you don't use mesh.save(filename) to create a copyable file .. its a little easier so we don't need a complicated back import

[Coastal0]

Sorry. I'd used mesh.saveAscii thinking it would be easier... lol.

pyGimli_files_new.zip

[carsten-forty2]

import pygimli as pg

m1 = pg.load('datamesh_pointcloud_border.bms')
print(m1) #Mesh: Nodes: 16519 Cells: 0 Boundaries: 0
m2 = pg.load('meshERT.bms')
print(m2) #Mesh: Nodes: 4607 Cells: 9067 Boundaries: 13673

m1 have no cells and hence can't be used for interpolation.

[Coastal0]

Thank you.

Meshing on that double node was causing my iPython to crash 40% of the time too. Seems stable so far!

print(pointcloud_data_mesh)
Mesh: Nodes: 8553 Cells: 16795 Boundaries: 25347

print(meshERT)
Mesh: Nodes: 4607 Cells: 9067 Boundaries: 13673

print(resBulk)
<class 'pygimli.core._pygimli_.RVector'> 4607 [0.0,...,0.0]