Exporting faults issue

[s-pearsongrant]

I have a 3D geological model and I am exporting faults for use in a reservoir model. If I only have faults offsetting formations it works fine, but if I have faults offsetting everything the fault surfaces become really random! I am using

rocktypes = [
geo_data.input_transform.apply_inverse(meshes.vertices)
for meshes in geo_data.solutions.dc_meshes
].

Am I doing something wrong, or is there more going on when the faults offset each other (hence why they have a double line in the attached cross-section)?
geological_model.html
geological_model_old.html

Many thanks,
Sophie

[NilsChudalla]

Hi @s-pearsongrant ,
that looks like a cool model! Would you mind sharing the structural frame with the fault matrix?

[s-pearsongrant]

@NilsChudalla
Thank you, I'm excited to get to grips with GemPy! Here is the structural frame - ohaaki and fault_h should truncate against other faults.
StructuralFrame(.txt

[NilsChudalla]

Have you tried only setting those faults to "True" in the fault relations matrix, if they actually are in proximity to each other and have an effect on each other within the model extent? I have the feeling that the scalarfield of the faults is off because it tries to consider all layers. Maybe could you also share a screen shot of the scalarfields of the faults in X direction, if the problem still persists?

[s-pearsongrant]

I tried to just set those faults to True, but the ones below in the matrix set to all False including the formation. One of the faults should be truncated against 2 faults, the other should be truncated against one fault (which also truncates the other one). They don't intersect each other. I don't know if I am setting it wrong?

These figures show the fault surfaces and I've mapped onto a grid for my modelling software.
Only offset formation

Offset all

The geological model with faults colour-coded

[NilsChudalla]

That second screenshot looks a lot different than the surfaces though... Interesting. Would you mind sharing the python script of the geomodel (only the Gempy part) including in the input? I would like to try a few different settings.

[s-pearsongrant]

I plotted the data points as well as the surfaces, and realised the surface reconstructions are so crazy because the faults don't actually truncate against other faults, they bend to follow the whole length of the other fault as well. Is there a way to make them truncate, or do I manually need to delete any that run alongside the fault they truncate against? This is causing issues when I try to map the faults onto my grid and make them continuous as well.

Many thanks,
Sophie

[NilsChudalla]

I think the refinement, which is setup in the model initialisation or before computation, is the required screw to turn. I might be able to give it a try, if you share the gempy Modeling steps with me. You can either upload a zip file here or invite me to a private repo if you are not comfortable putting it out in the open.

[s-pearsongrant]

Thank you very much! I have shared a private repo with you. This is my first time making a GemPy model so any help gratefully received.

[NilsChudalla]

I am having a go at it and I am experiencing some odd behaviour, trying to control the faults. I am trying to understand it, but I doubt that this on you, rather than Gempy :D

The doubling however is simply an artefact of the mesh extraction and Gempy being an implicit modeling software (example). Finite faults are tricky in many modeling approaches. However, this behaviour is consistent, so I do not understand, how behaviour from your screen shot "Offset all" could be explained.

Since your workflow seems to work as intendet for the "non-finite" faults (not affecting each other) and I do not see any urgency for the addition of this behaviour I would propose a intermediate solution: you could utilize the voxels from fault masking to clean up the cells in your final output. You can retrieve them in the order of the structural frame here:

fault_mask = geo_model.solutions.raw_arrays.mask_matrix[i]
# if you have to you can reshape it according to the models resolution
fault_mask_reshaped = fault_mask.reshape((geo_model.grid.regular_grid.resolution))

If possible in your simulation software, I would try adding the faults iteratively, starting with the minor ones, merging into the bigger ones (e.g. fault 4 and 5 before fault 1). Then, I would apply the mask of fault 1 to remove the values overlapping.

Let me know, if this workflow works for you.

Furthermore, I am really intrigued by your usecase. I would like to encourage you to contribute it to our examples. Especially, if the other software you are using is open source as well. Many other geoscientists may highly appreciate this!

[s-pearsongrant]

I can't get the mask to work because where faults are interacting the mask is filling all the area in between. I've tried deleting points within a certain distance of the previous fault but I have to put a wide tolerance to catch some points and then the faults that are truncating end up too short. Is there something I am missing or am I asking too much? If this works we are hoping to write up what we are doing and then would be happy to share it as a tutorial. Many thanks, Sophie

…

________________________________ From: NilsChudalla ***@***.***> Sent: Friday, October 24, 2025 3:09 AM To: gempy-project/gempy ***@***.***> Cc: Sophie Pearson-Grant ***@***.***>; Mention ***@***.***> Subject: Re: [gempy-project/gempy] Exporting faults issue (Discussion #1066) CAUTION: This email originated from outside of the organization. Do not click links or open attachments unless you recognize the sender and know the content is safe: I am having a go at it and I am experiencing some odd behaviour, trying to control the faults. I am trying to understand it, but I doubt that this on you, rather than Gempy :D The doubling however is simply an artefact of the mesh extraction and Gempy being an implicit modeling software (example<https://docs.gempy.org/tutorials/ch1_fundamentals/ch1_5_fault_relations.html>). Finite faults are tricky in many modeling approaches. However, this behaviour is consistent, so I do not understand, how behaviour from your screen shot "Offset all" could be explained. Since your workflow seems to work as intendet for the "non-finite" faults (not affecting each other) and I do not see any urgency for the addition of this behaviour I would propose a intermediate solution: you could utilize the voxels from fault masking to clean up the cells in your final output. You can retrieve them in the order of the structural frame here: fault_mask = geo_model.solutions.raw_arrays.mask_matrix[i] # if you have to you can reshape it according to the models resolution fault_mask_reshaped = fault_mask.reshape((geo_model.grid.regular_grid.resolution)) If possible in your simulation software, I would try adding the faults iteratively, starting with the minor ones, merging into the bigger ones (e.g. fault 4 and 5 before fault 1). Then, I would apply the mask of fault 1 to remove the values overlapping. Let me know, if this workflow works for you. Furthermore, I am really intrigued by your usecase. I would like to encourage you to contribute it to our examples. Especially, if the other software you are using is open source as well. Many other geoscientists may highly appreciate this! — Reply to this email directly, view it on GitHub<#1066 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/BG6J5T7YVHMRFH3TWOVCPLT3ZDOSDAVCNFSM6AAAAACJ25BFN6VHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTINZWGMYDAOI>. You are receiving this because you were mentioned.Message ID: ***@***.***> Sophie Pearson-Grant Geothermal Modeller +64 4 570 4712 1 Fairway Drive, Avalon 5011, New Zealand Earth Sciences New Zealand [Earth Sciences New Zealand]<https://earthsciences.nz> The Institute of Geological and Nuclear Sciences Limited and the National Institute of Water and Atmospheric Research Limited joined to become the New Zealand Institute for Earth Science Limited. We are known as Earth Sciences New Zealand. For more information on the Earth Sciences transition click here<https://www.gns.cri.nz/assets/share/ESNZ-2-Stakeholder-information-for-website.pdf>. Notice: This email may contain information which is confidential and/or subject to copyright, and may not be used, published or redistributed without the prior written consent of Earth Sciences New Zealand. If you are not the intended recipient, please immediately notify the sender and delete the email and any attachments. Any opinion or views expressed in this email are those of the individual sender and may not represent those of Earth Sciences New Zealand. For information about how we process data and monitor communications please see our privacy policy<https://www.gns.cri.nz/contact-us/privacy-policy>.

[NilsChudalla]

Just a small update: I am still looking into it and this hasn't gotten dust yet ;) I am looking into alternate ways of surface extraction.

[NilsChudalla]

In case, you have not received a notification, I did push some additions by using marching cubes to extract surfaces, rather than the default octree/dual contouring.

[NilsChudalla]

So, it took some time, but I got some progress. I had the issue that the model seemed to be quite unstable and I had to use pytorch to compute (otherwise computation was impossible and I got a singular matrix, which suprised me, given that we had more complex models already). This meant, I could not use marching cubes in the current state of gempy, therefore I wrote a small pull request. With marching cubes, I was able to compute the model with the faults offsetting in the way you asked and without these weird blobs that dual contouring produced before (seen in your previous screenshot).

You can adapt those changes in your Gempy Version and run the marching cubes function (see #1073). I did not know how you set your voxels, so built something around pyvista kdtree, that actually allows you controll the width of the faults.

It is not perfect, because I buffer around sections of the fault block, but you can reduce the amount of artefacts from mesh extraction. These artefacts are unavoidable in the current state (the surfaces extracted from the scalarfield will always be continuous). Without these buffers, but simply controling the order in which the voxels are assigned, you also get okay results.

Is this useable for your workflow?

[s-pearsongrant]

That is amazing, thank you! I will try it out asap. We noticed the singular issue, were thinking it must be because of the data we used - we tried to only put in what we know. We're planning to have more of a look at the data resolution soon. Unless you have thoughts what it might be? Thank you so so much, Sophie

…

________________________________ From: NilsChudalla ***@***.***> Sent: Friday, November 7, 2025 4:44:50 AM To: gempy-project/gempy ***@***.***> Cc: Sophie Pearson-Grant ***@***.***>; Mention ***@***.***> Subject: Re: [gempy-project/gempy] Exporting faults issue (Discussion #1066) CAUTION: This email originated from outside of the organization. Do not click links or open attachments unless you recognize the sender and know the content is safe: So, it took some time, but I got some progress. I had the issue that the model seemed to be quite unstable and I had to use pytorch to compute (otherwise computation was impossible and I got a singular matrix, which suprised me, given that we had more complex models already). This meant, I could not use marching cubes in the current state of gempy, therefore I wrote a small pull request. With marching cubes, I was able to compute the model with the faults offsetting in the way you asked and without these weird blobs that dual contouring produced before (seen in your previous screenshot). grafik.png (view on web)<https://github.com/user-attachments/assets/0ba41c71-aa0d-4b0f-944a-826759d9e8a3> You can adapt those changes in your Gempy Version and run the marching cubes function (see #1073<#1073>). I did not know how you set your voxels, so built something around pyvista kdtree<https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.KDTree.html>, that actually allows you controll the width of the faults. grafik.png (view on web)<https://github.com/user-attachments/assets/d090bf70-1506-4003-91a1-4b5f68907599> It is not perfect, because I buffer around sections of the fault block, but you can reduce the amount of artefacts from mesh extraction. These artefacts are unavoidable in the current state (the surfaces extracted from the scalarfield will always be continuous). Without these buffers, but simply controling the order in which the voxels are assigned, you also get okay results. grafik.png (view on web)<https://github.com/user-attachments/assets/e3294858-af8f-4fd8-a2e1-f7154d0f9b43> Is this useable for your workflow? — Reply to this email directly, view it on GitHub<#1066 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/BG6J5T6WW4B66OP2S5JWFUT33NUHFAVCNFSM6AAAAACJ25BFN6VHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTIOBZGMYDSNY>. You are receiving this because you were mentioned.Message ID: ***@***.***> Sophie Pearson-Grant Geothermal Modeller +64 4 570 4712 1 Fairway Drive, Avalon 5011, New Zealand Earth Sciences New Zealand [Earth Sciences New Zealand]<https://earthsciences.nz> The Institute of Geological and Nuclear Sciences Limited and the National Institute of Water and Atmospheric Research Limited joined to become the New Zealand Institute for Earth Science Limited. We are known as Earth Sciences New Zealand. For more information on the Earth Sciences transition click here<https://www.gns.cri.nz/assets/share/ESNZ-2-Stakeholder-information-for-website.pdf>. Notice: This email may contain information which is confidential and/or subject to copyright, and may not be used, published or redistributed without the prior written consent of Earth Sciences New Zealand. If you are not the intended recipient, please immediately notify the sender and delete the email and any attachments. Any opinion or views expressed in this email are those of the individual sender and may not represent those of Earth Sciences New Zealand. For information about how we process data and monitor communications please see our privacy policy<https://www.gns.cri.nz/contact-us/privacy-policy>.

[NilsChudalla]

I was playing around some more with the data as well. I tried applying some random noise in various strengths, doubling the amount of orientations (screenshot), and disabling the transformation (geo_model.update_transform(gp.data.GlobalAnisotropy.NONE)). Sadly, no success for the computation on the numpy backend.

I also noticed some offset for fault 4. The scalarfield seems decent though. I am a bit confused and am trying to have a look at it with a colleague.

And as a general FYI by the side (also for passerbys): when applying marching cubes, grid has to be adjusted and reset. The code before and after computation should look something like this

# Change the grid to only be the dense grid
dense_grid: RegularGrid = RegularGrid(
    extent=geo_model.grid.extent,
    resolution=np.array([40, 40, 40])
)

geo_model.grid.dense_grid = dense_grid
gp.set_active_grid(
    grid=geo_model.grid,
    grid_type=[geo_model.grid.GridTypes.DENSE],
    reset=True
)

geo_model.interpolation_options = gp.data.InterpolationOptions.init_dense_grid_options()

gp.compute_model(geo_model)

# # use marching cubes to extract the mesh
from gempy.modules.mesh_extranction import marching_cubes
marching_cubes.set_meshes_with_marching_cubes(geo_model)

Verts are stored per element and can be extracted (already transformed, unlike dc meshes) for example like this

# Extract vertices of faults from marching cubes solution
fault_vertices = []
for igroup in geo_model.structural_frame.structural_groups:
    if igroup.is_fault:
        for ielem in igroup.elements:
            fault_vertices.append(ielem.vertices)