Model dimensions are not correct

[AlexanderJuestel]

Describe the bug
I am creating a simple model with the extent [0, 1000, 0, 1000, -600, 0]. When looking at the PyVista plots, I noticed that the model is smaller than the extent, see screenshot.

The behavior becomes more apparent, when you want to calculate smaller submodels with a higher resolution and merge them later on. See second screenshot and the gaps between the models.

To Reproduce
Provide detailed steps to reproduce the behavior:

Building a model with the shown extent

Expected behavior
The meshes stretch all the way to the provided extent.

Screenshots

[javoha]

Hi Alex,
pretty sure I know why this is happening. It looks like the Pyvista plotter is not given the extent but instead matches what is displayed. The Dual contouring meshes that gempy v3 create seem not to cover the full area to the extent. If you set show_lith=True the correct extend should be displayed.
If this did not happen in the pervious gempy version, I would be pretty optimistic that you get meshes that cover the complete extent if you instead use marching cubes. I would need to try to confirm though.

Edit: Yes, in fact in my case the mc meshes work. Let me know if you need the code for that.

[AlexanderJuestel]

Yes, I can confirm that the lith_block covers the entire area.

If you have the code for the marching cubes at hand, I gladly take it.

[javoha]

Of course! This gets all the mc_meshes from the lithology block of a gempy model. Just note that it does not include faults. You would need to get them the same way from the corresponding fault blocks,

# get mc meshes
mc_vertices=[]
mc_edges=[]
block = geo_model.solutions.raw_arrays.lith_block.reshape(geo_model.grid.regular_grid.resolution)
for i in np.unique(block)[:-1]:
    verts, faces, _, _ = measure.marching_cubes(block, i,
                                                spacing=(geo_model.grid.regular_grid.dx,
                                                        geo_model.grid.regular_grid.dy,
                                                        geo_model.grid.regular_grid.dz))
    mc_vertices.append(verts)
    mc_edges.append(faces)

And the plotting part, the loop here hacked to loop over the number of elements (without faults)

# plotter = pv.Plotter(off_screen=not show_plotter)
plotter = pv.Plotter()

colors = ['#4285f4', '#ea4335', '#fbbc05', '#34a853', '#673ab7']
for i in range(3):
    plotter.add_mesh(
        pv.PolyData(mc_vertices[i],
                    np.insert(mc_edges[i], 0, 3, axis=1).ravel()),
                    color=colors[i])

plotter.show_bounds()
plotter.show()

[AlexanderJuestel]

I am trying to make my way around that now. The idea is to create multiple models and combine the meshes afterwards. Actually, since we are extracting the meshes using marching cubes from the lith blocks, we could combine the lith_blocks first and then extract the meshes from the lith_blocks. However, some formations may not be present in one submodel, therefore, we have to filter the interface points and orientations for each submodel and automatically create the mapping object and fault_group. We must also ensure that all elements that are present in the big model in one block are also present in the small blocks. This is i.e. the case for a fault that crops out in one block and is located in another block at a greater depth.....

This sounds more like hierarchical modeling again.....