# Masked Grid for Two Sides of a Fault

In this example, I demonstrate how to use a surface mesh of a fault in the subsurface to create a data mask on a modeling grid. This is a particularly useful exercise for scenarios where you may want to perform some sort of modeling in a different manner due to geological differences on the two sides of the fault - but still have a single modeling grid.

Let’s get to it!

```# sphinx_gallery_thumbnail_number = 4
import numpy as np
import pooch
import pyvista as pv
```
```url = "https://raw.githubusercontent.com/pyvista/vtk-data/master/Data/opal_mound_fault.vtk"
file_path = pooch.retrieve(url=url, known_hash=None)
fault
```
PolyDataInformation
N Cells3834
N Points2029
N Strips0
X Bounds3.361e+05, 3.389e+05
Y Bounds4.258e+06, 4.264e+06
Z Bounds-2.606e+03, 2.906e+03
N Arrays0

Create the modelling grid if you don’t already have one

```grid = pv.UniformGrid()
# Bottom south-west corner
grid.origin = (329700, 4252600, -2700)
# Cell sizes
grid.spacing = (500, 500, 500)
# Number of cells in each direction
grid.dimensions = (30, 35, 10)
grid
```
UniformGridInformation
N Cells8874
N Points10500
X Bounds3.297e+05, 3.442e+05
Y Bounds4.253e+06, 4.270e+06
Z Bounds-2.700e+03, 1.800e+03
Dimensions30, 35, 10
Spacing5.000e+02, 5.000e+02, 5.000e+02
N Arrays0

Take a quick preview to see where the fault is inside of the grid

```p = pv.Plotter()
p.show()
``` You may notice that the modeling grid’s extent is far greater than that of the fault – not to worry! PyVista’s clip_surface filter and the utility I’m going to share below handles this quite well by interpolating the fault’s plane outward.

This is a reusable utility for performing the mask:

```def mask_mesh_by_surface(mesh, surface):
grid = mesh.copy()
# Split the mesh by the fault
grid["pids"] = np.arange(grid.n_points)
grid["cids"] = np.arange(grid.n_cells)
a = grid.clip_surface(surface, invert=False, compute_distance=True)
b = grid.clip_surface(surface, invert=True, compute_distance=True)
# Use implicit distance to get point mask
lpids = np.argwhere(grid["implicit_distance"] >= 0)
gpids = np.argwhere(grid["implicit_distance"] < 0)
return grid
```

Let’s run it and take a look at the result!

```masked = mask_mesh_by_surface(grid, fault)

p = pv.Plotter()
p.show()
``` And here is how you might use that mask to do some sort of fancy modeling. In my example, I’m going to use a rather sophisticated distance calculation:

```ids = np.argwhere(masked["point_mask"] == 1).ravel()
pts = grid.points[ids]
len(pts)
```
```5659
```

Compute distance from TNE corner

```compute = lambda a, b: np.sqrt(np.sum((b - a) ** 2, axis=1))
dist = compute(pts, np.repeat([masked.bounds[1::2]], pts.shape, axis=0))
```

Add those results back to the source grid

```masked["cool_math"] = np.zeros(grid.n_points)  # Need to preallocate

# Do some different math for the other side
...
```
```Ellipsis
```

Display!

```masked.plot(scalars="cool_math")
``` Visualize one side of the masked grid

```a = masked.threshold(1.5, scalars="cell_mask", invert=True)

p = pv.Plotter() 