3D Mesh creation

[tabdivider]

Problem description

I have a topography pointcloud and need to translate it into a 3D mesh with a given lower boundary, e.g. 200 m asl. Is there a simple way in pyGIMLi meshtools to import the xyz points and define a 3D mesh with topography to make a traveltime simulation similar to the Refraction in 3D example?

Your environment

Operating system: Windows 10/64b Python version: 3.7.3 pyGIMLi version: 0+unknown (??)/ 1.0.12 Way of installation: Conda package

[tabdivider]

Hello

I guess, I found a solution that work for me. I used BERT to generate the 3D mesh with topography and then I load the mesh into pyGIMLi

fname = 'mesh\mesh'
full_mesh = pg.meshtools.mesh.readTetgen(fname)

and subsequently separating the fine inner mesh from the course mesh around the model space by filtering the mesh using the region markers

mesh = pg.Mesh()
mesh.createMeshByMarker(full_mesh, 2)

I hope that the mesh is just fine and no additional adjustments or settings are needed (?). I have some first results made on the mesh and it seems to be allright. It would be great if the 3D mesh generation utility from BERT could be impoted into pyGIMLi as well.

I was also struggling a bit with the setup of the gradient initial model to force the raypath travel through the underground, as the useGradient function works apparently for 2D only. I fixed that using code parts from #186 .

vel_gradient = [] for node in mesh.nodes(): vel_gradient.append( a + (-b * (abs(node.z()) ** 2)) ) vel_gradient = pg.meshtools.nodeDataToCellData(mesh,np.array(vel_gradient))

[carsten-forty2]

Hello David,

sorry for the late reply.

Good you found a clever working solution.

The mesh for Refraction problems does not need to be to fine like meshes for ERT as you can play with the number of secondary nodes to increase numerical accuracy.

The 3D meshtools for pygimli-1.1 should by able to emulate the behavior of the BERT scripts, however I don't yet found the time to combine them into a nice single method.

When you are running conda you can also try the precompiled release candidates for pygimli 1.1 listed here:

https://anaconda.org/gimli/pygimli/files

with 

conda install -c gimli/label/test pygimli

I recommend using conda environments for these so you can easily switch versions.

Best regards,

Carsten

davidkusnirak <notifications@github.com> hat am 3. Februar 2020 17:09 geschrieben:

Hello I guess, I found a solution that work for me. I used BERT to generate the 3D mesh with topography and then I load the mesh into pyGIMLi fname = 'mesh\mesh' full_mesh = pg.meshtools.mesh.readTetgen(fname)

and subsequently separating the fine inner mesh from the course mesh around the model space by filtering the mesh using the region markers mesh = pg.Mesh() mesh.createMeshByMarker(full_mesh, 2)

I hope that the mesh is just fine and no additional adjustments or settings are needed (?). I have some first results made on the mesh and it seems to be allright. It would be great if the 3D mesh generation utility from BERT could be impoted into pyGIMLi as well. I was also struggling a bit with the setup of the gradient initial model to force the raypath travel through the underground, as the useGradient function works apparently for 2D only. I fixed that using code parts from #186 . vel_gradient = [] for node in mesh.nodes(): vel_gradient.append( a + (-b * (abs(node.z()) ** 2)) ) vel_gradient = pg.meshtools.nodeDataToCellData(mesh,np.array(vel_gradient)) —You are receiving this because you are subscribed to this thread.Reply to this email directly, view it on GitHub, or unsubscribe.

[tabdivider]

Hello Carsten,

thanks for your update. I do already use the pre-release of the 1.1, mainly for ploting and I have to say that I love the 3D Viewer. Really looking forward for the stable release of the 1.1!

Cheers, David

[tabdivider]

Hello again, I'm running the 3D inversion under 1.0.12.1, however the forward calculation seems to be the bottleneck at the moment as it is utilizing just single core, as discussed in #189. Is there a way how to use the new Dijkstra forward computation from v1.1? I haven't found the function in the 1.1rc6, but could be that I'm not looking in the right way.. Could you give a hint where to look for the new fowrard function or how to adopt it in 1.0.12.1?

[halbmy]

The forward computation should be the pretty much the same in pg 1.0 and 1.1.

You can set the number of CPUs by e.g. fop.setThreadCount(16) where fop is the forward operator. Note sure whether there is a good CPU autodetect or whether it works well on Windows. Just try.

[tabdivider]

Seems that setting

ra = Refraction('inputSGTfile.sgt',debug=True)
ra.fop.setThreadCount(16)  

had no effect on the perfomance..still running on a single core only. Happening during this stage min/max(dweight)=

[halbmy]

Please have a closer look again. The inverse solver (coming directly after the min/max messages) is still not parallelized and uses only 1 core. But when it comes to the forward computation and Jacobian matrix it should use the number of cores. If you cannot observe that, we will have a deeper look.

[tabdivider]

I will continue on that and let you know.. BTW..I was deploying the code on a Xeon based workstation, where the code used all the CPU power available at least for the Jacobian calculation. Now I switched on an AMD Ryzen maschine (due to the more memory and cores available) and seems that everyhing is running on a single core. Sure, that the architecture could also have an impact too..

[tabdivider]

I fixed the hardware problem...the problem was caused by another process that was running in parallel on the same python environment, so it was just a local issue, not related to the pyGIMLi.

I got an reasonable tomographic output, but I'd like to improve the model, by inserting the know cavity into the mesh and rerun the simulation. I'm fighting with the gmsh and the input and output formats, but apparently I'm not getting any closer. What I have at the moment is the original mesh (MainMesh.bms) and a known cavity mesh (cavity.bms). My goal is to properly merge the two meshes, where the known cavity is really an empty space. After many trials and errors, seems I'm in a dead end..

cavity.bms.zip

MainMesh.bms.zip

[carsten-forty2]

with 1.1-rc you can try the following (just a rough untested idea):

main = pg.load('MainMesh.bms')
cavity = pg.load('cavity.bms')

mainBoundary = main.createSubMesh(main.boundaries(main.boundaryMarkers()!=0))
#pg.show(mainBoundary)
cavityBoundary = cavity.createSubMesh(cavity.boundaries(cavity.boundaryMarkers()!=0))
#pg.show(cavityBoundary)
mesh = pg.meshtools.createMesh([mainBoundary, cavityBoundary])
#pg.show(mesh)

The api is supposed to working like this, if it doesn't work I will look in more detail or fix it this way.

[tabdivider]

Hi Carsten,

I'm happy that you see a straightforward solution! I got an error in the

mesh = pg.meshtools.createMesh([mainBoundary, cavityBoundary])

Opening C:\Users\DK\AppData\Local\Temp\tmpodk1vrdc.poly.
PLC Error:  A segment and a facet intersect at point (95.6495,16.157,242.464).
  Segment: [14308,17627] #-1 (-1)
  Facet:   [6986,7430,7596] #103
A self-intersection was detected. Program stopped.
Hint: use -d option to detect all self-intersections.

I also tried to use the -d switch but no success

mesh = pg.meshtools.createMesh([mainBoundary, cavityBoundary],switches='-d')

[frodo4fingers]

hmm... without any checking...

replacing the last line to

mesh = pg.meshtools.merge2Meshes(mainBoundary, cavityBoundary)

looks at least fine

[tabdivider]

that just joined the boundaries (no cells are created)

though, when I use the new merged boundary file to create the mesh

mesh_bound = pg.meshtools.merge2Meshes(mainBoundary, cavityBoundary)
pg.show(mesh_bound)
mesh = pg.meshtools.createMesh(mesh_bound,quality = 33)
pg.show(mesh)

I'm getting the same error

[florian-wagner]

This looks to me like a Gmsh issue rather than a pygimli one. The mesh with the cavity should be properly designed in Gmsh. No mesh modification should be needed afterwards. Or you define the geometry in Gmsh and follow @carsten-forty2 advice to fill the geometry with cells in pygimli afterwards. But either way, the geometry must be unambiguously defined in Gmsh before and I suppose there lies the problem. If you post your .geo file, I can have a look.

[tabdivider]

you are right...I changed the geometry of the cylinder in Gmsh, so it is little bit shorter than the extents of the mainMesh and now the the suggestion from @carsten-forty2 works just perfect! Thank you for the support. I will now design the real, complex shape of the cavity and re-use this mesh procedure.

Btw..I have now two markers classes in the mesh (1&2), where 1 is the cavity and 2 is the subsurface. My idea is now to exclude all the cells marked with 1 and create a 'swiss cheese' mesh with a hole inside. I hope this will work just fine. however, just out of curiosty..is there a way how to use the complete mesh and run the inversion just on cells with marker2? Is there any benefit of this approach (apart from slightly higher number of cells)?

[carsten-forty2]

ah .. I forgot the cavity.

just add

cavityBoundary.addHoleMarker([x,y,z])

with x,y,z some point inside the cavity

[carsten-forty2]

You can't avoid the cells with marker == 1 to be used in the forward calculation of the inversion. But in principle you can define the region with marker == 1 as a fixed region high slowness to be "not seen" by the inversion. However, I'm not sure if the current state of the TravelTimeManager is able to do so out of the box.