How to delete the points and cells above irregular surface?

[handechao]

Hi guys! Here is my question.
I have a segy model with an irregular surface. The Vp above the surface is 0 m/s. Using 'domain=rectangle' is not worked and mesh above the irregular surface, which is not desired, will also be generated.
Is there any function that can delete those points and cells above the surface? Or, whether it can skip the region in which Vp is 0? Or, Should I provide the coordinates of the irregular surface to define the geometry?

[krober10nd]

So the domain must be a signed distance function. You can create this function by thresholding the velocity where Vp > 0 m/s and setting these areas equal to 1 and elsewhere equal to -1. Then you can use a method such as the fast marching method (from skfmm) to turn this field into a discrete signed distance field. Finally, I would create a Scipy.RegularGriddedInterpolant from that signed distance field that takes a vector of points and returns a vector of signed distances.

Now once I have the signed distance function, I'd pass it to generate_mesh which excepts either a :class:Domain or a regular function. Since it's not a Domain, you'll have to pass an extent bbox as well.

Keep in mind if you use get_sizing_function_from_segy then you'll have to make sure the bbox and your velocity model are orientated in z, x, y axes order and this needs to be consistent when you define your signed distance function below.

So in summary here's some pseudo code that is untested since I do not have the Vp field. I will say, I have done this before and it worked well. Happy meshing.

   import numpy as np
   import skfmm
   
   bbox = () # your bounds in the same length units as your velocity model
   
   # create a signed distance field based on velocity Vp
 phi = np.ones_like(Vp)
 phi[Vp < 0] = -1
 # call fast marching method
 d = skfmm.distance(phi, [dx, dy, dz])
 
  # create a regular grid interpolator
  # often I smooth the discrete signed distance but try without that at first.
  # this can make the boundary have higher quality elements.
  # d = sp.ndimage.filters.gaussian_filter(d, sigma, mode="constant")
   
 nx, ny, nz = Vp.shape
 z = np.linspace(bbox[0], bbox[1], nz)
 x = np.linspace(bbox[2], bbox[3], nx)
 y = np.linspace(bbox[4], bbox[5], ny)
 X, Y, Z = np.meshgrid(x, y, z)
 interpolant = RegularGridInterpolator((x, y, z), d, bounds_error=False, fill_value=None)
 
  # now make your signed_distance_function
 def my_signed_distance_function(points):
     return interpolant(points)
   
   
 # make your edge length function per the README for the other examples if you'd like
 # ef = get_sizing_function_from_segy(...)
   
  # call the mesh generator
 points, cells = generate_mesh(
     domain=my_signed_distance_function, bbox=bbox, edge_length=ef
 )

[krober10nd]

@handechao did this help you out? Please let me know if you need for me to clarify anything further.

[handechao]

I have tried this code, but it still didn't work. And I found that, in my code, if domain_pad=0, the value of mesh coordinates will be huge. It seems like non-convergence.
Maybe I made some mistakes. So I attach a velocity file (binary file, float32) with an irregular surface for testing.

Here is my code.

# create a signed distance field based on velocity vp
# I read the segy file using _read_segy function in this repository
phi = np.ones_like(vp)
phi[vp==0] = -1
# call fast marching method
d = skfmm.distance(phi, [dz, dx])  # d.shape=(nz,nx) is it correct?

# create a regular grid interpolator  
z = np.linspace(bbox[0], bbox[1], nz)
x = np.linspace(bbox[2], bbox[3], nx)
X, Z = np.meshgrid(x, z)
interpolant = RegularGridInterpolator((z, x), d, bounds_error=False, fill_value=None) # d.shape=(nz,nx) so i pride (z,x) as the grid
 
# now make your signed_distance_function
def my_signed_distance_function(points):
    return interpolant(points)

# Construct mesh sizing object from velocity model
ef = get_sizing_function_from_segy(
    fname,
    bbox,
    hmin=18,
    hmax=80,
    wl=3,
    freq=25,
    grade=0.15,
    dt=0.001,
    domain_pad=200, # setting pad = 0 will lead to non-convergence mentioned above.
    space_order=2,
)

# points, cells = generate_mesh(domain=rectangle, edge_length=ef)
points, cells = generate_mesh(domain=my_signed_distance_function, bbox=bbox, edge_length=ef)

vp_nz1000dz10_nx1668dx15.zip

[krober10nd]

still a work in progress but I added some functions to the open PR (misc fixes) that enable you to read in your velocity model.

I'll come back to this in a couple days.

   import meshio
   # import matplotlib.pyplot as plt
   import numpy as np
   import skfmm
   from scipy.interpolate import RegularGridInterpolator
   from scipy.ndimage import filters
   
   import SeismicMesh as sm
   
   filename = "vp_nz1000dz10_nx1668dx15.bin"
   
   dz = 10  # m
   dx = 15  # m
   nz = 1000  # grid points
   nx = 1668  # grid poins
   hmin = 100  # m
   bbox = (  # m
       -10000.0,
       0.0,
       0.0,
       25020.0,
   )  # your bounds in the same length units as your velocity model
   sigma = 100  # m^2 for spatial filtering discrete signed distance function
   
   WL = 10  # 10 cells per wavelength
   FREQ = 2  # for a peak frequency of 2 Hz
   
   vp, nz, nx, ny = sm.read_velocity_model(
       filename,
       nz=1000,
       nx=1668,
       ny=1,
       axes_order_sort="F",
       axes_order=(0, 1, 2),
       byte_order="little",
       dtype="float32",
   )
   vp = np.squeeze(vp)
   print(f"THE MINIMUM VELOCITY IS {np.amin(vp)}")
   print(f"THE MAXIMUM VELOCITY IS {np.amax(vp)}")
   
   # create a signed distance field based on velocity Vp
   phi = np.ones_like(vp)
   phi[vp > 0] = -1
   # call fast marching method
   d = skfmm.distance(phi, [dz, dx])
   
   # create a regular grid interpolator
   # often I smooth the discrete signed distance but try without that at first.
   # this can make the boundary have higher quality elements.
   # d = filters.gaussian_filter(d, sigma, mode="constant")
   
   z = np.linspace(bbox[0], bbox[1], nz)
   x = np.linspace(bbox[2], bbox[3], nx)
   
   # Z, X = np.meshgrid(z, x)
   # fig, ax = plt.subplots()
   # ax.pcolor(X, Z, d.T, vmin=-0.01, vmax=0.01)
   # ax.axis("equal")
   # plt.title("Discrete signed distance function")
   # plt.savefig("discrete_signed_distance.png")
   
   interpolant = RegularGridInterpolator((z, x), d, bounds_error=False, fill_value=None)
   
   
   def my_signed_distance_function(points):
       return interpolant(points)
   
   
   # make your edge length function per the README for the other examples if you'd like
   ef = sm.get_sizing_function_from_segy(
       None,
       bbox,
       hmin=hmin,
       velocity_data=vp,
       nz=nz,
       nx=nx,
       wl=WL,
       freq=FREQ,
       axes_order_sort="F",
       axes_order=(0, 1, 2),
       byte_order="little",
       dtype="float32",
       units="m/s",
   )
   
   sm.plot_sizing_function(ef)
   
   # call the mesh generator
   points, cells = sm.generate_mesh(
       domain=my_signed_distance_function,
       bbox=bbox,
       edge_length=ef,
       mesh_improvement=False,
   )
   
   meshio.write_points_cells(
       "irregular_mesh.vtk",
       points[:, [1, 0]],  # FOR VISUALIZATION IN PARAVIEW
       [("triangle", cells)],
       file_format="vtk",
   )

[handechao]

Thank you for kindly help.
The code in my last reply looks similar to your pseudo code since I changed only 2 lines for testing.

I see that a new branch 'misc' is added. I tried pip install git+https://github.com/krober10nd/SeismicMesh.git@misc to rebuilding SeismicMesh, but it showed AttributeError: module 'SeismicMesh' has no attribute 'read_velocity_model'. Maybe I'm wrong.
So I copy the _read_velocity_model and related functions into my code to use it. And I added segy header in vp file since I can't use the newest get_sizing_function_from_segy.

I used the spatial filter on discrete signed distance and set a smaller sigma=5. The result is good with only a few poor quality mesh along the boundary. (picture below)

If I run the code without the spatial filter, here is the mesh picture.

If I run the code with the spatial filter sigma=100, here is the mesh picture.

If I run the code with domain_pad=200, sigma=5, it seems didn't work. Here is the mesh picture.

SO is there any suggestions for setting proper parameters, like sigma, domain_pad, mesh_improvement?

segy vp file is attached.

Here is my code.

import meshio
import matplotlib.pyplot as plt
import numpy as np
import skfmm
import warnings
import segyio
from scipy.interpolate import RegularGridInterpolator
from scipy.ndimage import filters

import SeismicMesh as sm

filename = "vp_nz1000dz10_nx1668dx15.segy"

####################################################
# read vp from file copy from sizing folder
def read_velocity_model(
    filename,
    nz=None,
    nx=None,
    ny=None,
    byte_order=None,
    axes_order=None,
    axes_order_sort=None,
    dtype=None,
):
    """Read a velocity model"""
    if filename.endswith(".segy"):
        return read_segy(filename)
    else:
        return read_bin(
            filename, nz, nx, ny, byte_order, axes_order, axes_order_sort, dtype
        )

def read_bin(filename, nz, nx, ny, byte_order, axes_order, axes_order_sort, dtype):
    """Read a velocity model from a binary"""
    if (nz is None) or (nx is None) or (ny is None):
        raise ValueError(
            "Please specify the number of grid points in each dimension (e.g., `nz`, `nx`, `ny`)..."
        )
    axes = [nz, nx, ny]
    ix = np.argsort(axes_order)
    axes = [axes[o] for o in ix]
    with open(filename, "r") as file:
        print(f"Reading binary file: {filename}")
        if byte_order == "big":
            vp = np.fromfile(file, dtype=np.dtype(dtype).newbyteorder(">"))
        elif byte_order == "little":
            vp = np.fromfile(file, dtype=np.dtype(dtype).newbyteorder("<"))
        else:
            raise ValueError("Please specify byte_order as either: little or big.")
        vp = vp.reshape(*axes, order=axes_order_sort)

        vp = np.flipud(vp.transpose((*axes_order,)))

        return vp, nz, nx, ny  # z, x and then y

def read_segy(filename):
    """Read a velocity model from a SEG-y file"""

    with segyio.open(filename, ignore_geometry=True) as f:
        nz, nx = len(f.samples), len(f.trace)
        vp = np.zeros(shape=(nz, nx))
        for index, trace in enumerate(f.trace):
            vp[:, index] = trace
        if np.amin(vp) < 1000.0:
            warnings.warn(
                "Velocity appear to be in km/s. Maybe pass `units` km-s key pair?"
            )
        return np.flipud(vp), nz, nx, 0
####################################################
# 
####################################################
dz = 10  # m
dx = 15  # m
nz = 1000  # grid points
nx = 1668  # grid poins
hmin = 40  # m
bbox = (  # m
    -10000.0,
    0.0,
    0.0,
    25020.0,
) 
sigma = 5  # small sigma result in high quality

WL = 3  # 3 cells per wavelength
FREQ = 30  # for a peak frequency of 30 Hz

vp, nz, nx, ny = read_velocity_model(
    filename,
)
vp = np.squeeze(vp)
print(f"THE MINIMUM VELOCITY IS {np.amin(vp)}")
print(f"THE MAXIMUM VELOCITY IS {np.amax(vp)}")

# create a signed distance field based on velocity Vp
phi = np.ones_like(vp)
phi[vp > 0] = -1
# call fast marching method
d = skfmm.distance(phi, [dz, dx])

# create a regular grid interpolator
# often I smooth the discrete signed distance but try without that at first.
# this can make the boundary have higher quality elements.
d = filters.gaussian_filter(d, sigma, mode="constant")

z = np.linspace(bbox[0], bbox[1], nz)
x = np.linspace(bbox[2], bbox[3], nx)

interpolant = RegularGridInterpolator((z, x), d, bounds_error=False, fill_value=None)

def my_signed_distance_function(points):
    return interpolant(points)

# still, I use segy file
ef = sm.get_sizing_function_from_segy(
    filename,
    bbox,
    hmin=hmin,
    nz=nz,
    nx=nx,
    wl=WL,
    freq=FREQ,
    hmax=80,
    domain_pad=0,
)

sm.plot_sizing_function(ef)

# call the mesh generator
points, cells = sm.generate_mesh(
    domain=my_signed_distance_function,
    bbox=bbox,
    edge_length=ef,
    mesh_improvement=False, 
)

meshio.write_points_cells(
    "irregular_mesh.vtk",
    points[:, [1, 0]],  # FOR VISUALIZATION IN PARAVIEW
    [("triangle", cells)],
    file_format="vtk",
)

SEGY File:
vp_nz1000dz10_nx1668dx15.zip

[handechao]

sorry for pasting wrong picture of 'without spatial filter'. Here is the right one.

[krober10nd]

@handechao you should be able to install the latest branch misc and use my code above. Anyway, you arrived at the same conclusion.

I have to tune the spatial filter sigma to get the right representation of the irregular surface. I've done it before, I just need to find my code. Will take a couple days.

[krober10nd]

If you run with a small sigma, and run SeismicMesh.geometry.delete_boundary_elements() to remove some of the bad ones might get what you want. To add the pad requires more changes.

[krober10nd]

Set the spatial_filter to some multiple of the minimum edge length. so

min_edge_length = 1500 / (max_freq * cells_per_wavelength)
filter = min_edge_length / 20 

seems to work for a variety of cases. Clean up the remaining degenerate elements with SeismicMesh.geometry.delete_boundary_elements() as I said earlier.

[handechao]

Yes, I have tested a 3d model last month. The result was not well (see the figure). Following is my code. For some reason, I can't provide the model. I will design a simple 3d model for testing and upload it.

import meshio
from mpi4py import MPI
import matplotlib.pyplot as plt
import numpy as np
import skfmm
import warnings
import segyio
from scipy.interpolate import RegularGridInterpolator
from scipy.ndimage import filters

from SeismicMesh import (
    get_sizing_function_from_segy,
    generate_mesh,
    sliver_removal,
    Cube,
)
from SeismicMesh.geometry import delete_boundary_entities

def read_velocity_model(
    filename,
    nz=None,
    nx=None,
    ny=None,
    byte_order=None,
    axes_order=None,
    axes_order_sort=None,
    dtype=None,
):
    """Read a velocity model"""
    if filename.endswith(".segy"):
        return read_segy(filename)
    else:
        return read_bin(
            filename, nz, nx, ny, byte_order, axes_order, axes_order_sort, dtype
        )


def read_bin(filename, nz, nx, ny, byte_order, axes_order, axes_order_sort, dtype):
    """Read a velocity model from a binary"""
    if (nz is None) or (nx is None) or (ny is None):
        raise ValueError(
            "Please specify the number of grid points in each dimension (e.g., `nz`, `nx`, `ny`)..."
        )
    axes = [nz, nx, ny]
    ix = np.argsort(axes_order)
    axes = [axes[o] for o in ix]
    print(ix,axes)
    with open(filename, "r") as file:
        print(f"Reading binary file: {filename}")
        if byte_order == "big":
            vp = np.fromfile(file, dtype=np.dtype(dtype).newbyteorder(">"))
        elif byte_order == "little":
            vp = np.fromfile(file, dtype=np.dtype(dtype).newbyteorder("<"))
        else:
            raise ValueError("Please specify byte_order as either: little or big.")
        vp = vp.reshape(*axes, order=axes_order_sort)

        vp = np.flipud(vp.transpose((*axes_order,)))

        return vp, nz, nx, ny  # z, x and then y


def read_segy(filename):
    """Read a velocity model from a SEG-y file"""

    with segyio.open(filename, ignore_geometry=True) as f:
        nz, nx = len(f.samples), len(f.trace)
        vp = np.zeros(shape=(nz, nx))
        for index, trace in enumerate(f.trace):
            vp[:, index] = trace
        if np.amin(vp) < 1000.0:
            warnings.warn(
                "Velocity appear to be in km/s. Maybe pass `units` km-s key pair?"
            )
        return np.flipud(vp), nz, nx, 0


############################################################
############################################################
comm = MPI.COMM_WORLD
bbox=(-15300.0, 0.0, 0.0, 29880.0, 0.0, 17600.0)
hmin = 100
sigma = 2
fname = "03-Shanqiandai_3DModel_Def_02_less.bin"
cube = Cube(bbox)
nz=1531
nx=748
ny=441
dz=10
dx=40
dy=40



vp, nz, nx, ny = read_velocity_model(
    fname,
    nz=1531,
    nx=748,
    ny=441,
    byte_order="little",
    dtype="float32",
    axes_order=(0,1,2),
    )

vp = np.squeeze(vp)
print(f"THE MINIMUM VELOCITY IS {np.amin(vp)}")
print(f"THE MAXIMUM VELOCITY IS {np.amax(vp)}")

# create a signed distance field based on velocity Vp
phi = np.ones_like(vp)
phi[vp > 0] = -1
# call fast marching method
d = skfmm.distance(phi, [dz, dx, dy])
d = filters.gaussian_filter(d, sigma, mode="constant")

z = np.linspace(bbox[0], bbox[1], nz)
x = np.linspace(bbox[2], bbox[3], nx)
y = np.linspace(bbox[4], bbox[5], ny)

interpolant = RegularGridInterpolator((z, x, y), d, bounds_error=False, fill_value=None)


def my_signed_distance_function(points):
    return interpolant(points)

ef = get_sizing_function_from_segy(
    fname,
    bbox,
    hmin=hmin,
    dt=0.001,
    freq=15,
    wl=5,
    grade=0.15,
    hmax=500,
    nz=nz,
    nx=nx,
    ny=ny,
    domain_pad=0,
    pad_style="edge",
    byte_order="little",
    # axes_order_sort="F",
)

# points, cells = generate_mesh(domain=cube, edge_length=ef, max_iter=75)
points, cells = generate_mesh(domain=my_signed_distance_function, 
    bbox=bbox,
    edge_length=ef, max_iter=75)

# For 3D mesh generation, we provide an implementation to bound the minimum dihedral angle::
# We use the preserve kwarg to ensure the level-set is very accurately preserved.
points, cells = sliver_removal(
    points=points, bbox=bbox, domain=my_signed_distance_function, edge_length=ef, preserve=True
)


# Meshes can be written quickly to disk using meshio and visualized with ParaView::
if comm.rank == 0:

    # NOTE: SeismicMesh outputs assumes the domain is (z,x,y) so for visualization
    # in ParaView, we swap the axes so it appears as in the (x,y,z) plane.
    meshio.write_points_cells(
        "lms3d.vtk",
        points[:, [1, 2, 0]] / 1000.0,
        [("tetra", cells)],
    )

[krober10nd]

You can send large files to my email keithrbt0 at gmail . com

[handechao]

The copyright owner is not us. Sorry for that.

[handechao]

Hi, Here is a small 3D model for testing.
irr_X193x25m_Y186x25m_Z105x25m.zip
It warned FAILURE: Termination...maximum number of iterations reached. Try increasing max_iter when generating the mesh. I set max_iter=100 in sliver_removal.

import meshio
from mpi4py import MPI
import matplotlib.pyplot as plt
import numpy as np
import skfmm
import warnings
import segyio
from scipy.interpolate import RegularGridInterpolator
from scipy.ndimage import filters

from SeismicMesh import (
    get_sizing_function_from_segy,
    generate_mesh,
    sliver_removal,
    Cube,
)
from SeismicMesh.geometry import delete_boundary_entities

def read_velocity_model(
    filename,
    nz=None,
    nx=None,
    ny=None,
    byte_order=None,
    axes_order=None,
    axes_order_sort=None,
    dtype=None,
):
    """Read a velocity model"""
    if filename.endswith(".segy"):
        return read_segy(filename)
    else:
        return read_bin(
            filename, nz, nx, ny, byte_order, axes_order, axes_order_sort, dtype
        )

def read_bin(filename, nz, nx, ny, byte_order, axes_order, axes_order_sort, dtype):
    """Read a velocity model from a binary"""
    if (nz is None) or (nx is None) or (ny is None):
        raise ValueError(
            "Please specify the number of grid points in each dimension (e.g., `nz`, `nx`, `ny`)..."
        )
    axes = [nz, nx, ny]
    ix = np.argsort(axes_order)
    axes = [axes[o] for o in ix]
    print(ix,axes)
    with open(filename, "r") as file:
        print(f"Reading binary file: {filename}")
        if byte_order == "big":
            vp = np.fromfile(file, dtype=np.dtype(dtype).newbyteorder(">"))
        elif byte_order == "little":
            vp = np.fromfile(file, dtype=np.dtype(dtype).newbyteorder("<"))
        else:
            raise ValueError("Please specify byte_order as either: little or big.")
        vp = vp.reshape(*axes, order=axes_order_sort)

        vp = np.flipud(vp.transpose((*axes_order,)))

        return vp, nz, nx, ny  # z, x and then y

def read_segy(filename):
    """Read a velocity model from a SEG-y file"""

    with segyio.open(filename, ignore_geometry=True) as f:
        nz, nx = len(f.samples), len(f.trace)
        vp = np.zeros(shape=(nz, nx))
        for index, trace in enumerate(f.trace):
            vp[:, index] = trace
        if np.amin(vp) < 1000.0:
            warnings.warn(
                "Velocity appear to be in km/s. Maybe pass `units` km-s key pair?"
            )
        return np.flipud(vp), nz, nx, 0


############################################################
# main script #
############################################################
comm = MPI.COMM_WORLD
bbox=(-2600.0, 0.0, 0.0, 4800.0, 0.0, 4650.0)
hmin = 20
sigma = 2
fname = "irr_X193x25m_Y186x25m_Z105x25m.bin"
cube = Cube(bbox)
nz=105
nx=193
ny=186
dz=25
dx=25
dy=25

vp, nz, nx, ny = read_velocity_model(
    fname,
    nz=nz,
    nx=nx,
    ny=ny,
    byte_order="little",
    dtype="float32",
    axes_order=(0,1,2),
    )

vp = np.squeeze(vp)
print(f"THE MINIMUM VELOCITY IS {np.amin(vp)}")
print(f"THE MAXIMUM VELOCITY IS {np.amax(vp)}")

# create a signed distance field based on velocity Vp
phi = np.ones_like(vp)
phi[vp > 0] = -1
# call fast marching method
d = skfmm.distance(phi, [dz, dx, dy])
d = filters.gaussian_filter(d, sigma, mode="constant")

z = np.linspace(bbox[0], bbox[1], nz)
x = np.linspace(bbox[2], bbox[3], nx)
y = np.linspace(bbox[4], bbox[5], ny)

interpolant = RegularGridInterpolator((z, x, y), d, bounds_error=False, fill_value=None)

def my_signed_distance_function(points):
    return interpolant(points)

ef = get_sizing_function_from_segy(
    fname,
    bbox,
    hmin=hmin,
    dt=0.001,
    freq=15,
    wl=5,
    grade=0.15,
    hmax=100,
    nz=nz,
    nx=nx,
    ny=ny,
    domain_pad=0,
    pad_style="edge",
    byte_order="little",
)

points, cells = generate_mesh(domain=my_signed_distance_function, 
    bbox=bbox,
    edge_length=ef, max_iter=100)

# bound the minimum dihedral angle:
points, cells = sliver_removal(
    points=points, bbox=bbox, domain=my_signed_distance_function, edge_length=ef, preserve=True,max_iter=100
)

# Meshes can be written quickly to disk using meshio and visualized with ParaView::
if comm.rank == 0:
    meshio.write_points_cells(
        "irr3d.vtk",
        points[:, [1, 2, 0]] / 1000.0,
        [("tetra", cells)],
    )

[krober10nd]

thank you @handechao I will find some time for this over the weekend.

[krober10nd]

@handechao The binary file wasn't be read in correctly. I'd encourage you to plot slices to do a sanity check if you do this next time.

The second thing was your 0-level set didn't pass through the domain (partly the result of the fact you weren't reading the file in correctly), so I increased the threshold to 0.1 m/s instead of 0 m/s. I used a Gaussian Blur again to regularize the 0-level set with a sigma of 0.5 grid units. In the future, if you're signed distance never crosses 0 inside the domain, this indicates the 0-level set is not contained within the domain.

Now the current issue is the filtering which will need some more experimentation to get right.

import meshio
import numpy as np
import skfmm
from mpi4py import MPI
from scipy.interpolate import RegularGridInterpolator
from scipy.ndimage import filters

from SeismicMesh import (generate_mesh, get_sizing_function_from_segy,
                         sliver_removal)


def read_velocity_model(
    filename,
    nz=None,
    nx=None,
    ny=None,
    byte_order=None,
    axes_order=None,
    axes_order_sort=None,
    dtype=None,
):
    """Read a velocity model"""
    return read_bin(
        filename, nz, nx, ny, byte_order, axes_order, axes_order_sort, dtype
    )


def read_bin(filename, nz, nx, ny, byte_order, axes_order, axes_order_sort, dtype):
    """Read a velocity model from a binary"""
    if (nz is None) or (nx is None) or (ny is None):
        raise ValueError(
            "Please specify the number of grid points in each dimension (e.g., `nz`, `nx`, `ny`)..."
        )
    axes = [nz, nx, ny]
    ix = np.argsort(axes_order)
    axes = [axes[o] for o in ix]
    with open(filename, "r") as file:
        print(f"Reading binary file: {filename}")
        if byte_order == "big":
            vp = np.fromfile(file, dtype=np.dtype(dtype).newbyteorder(">"))
        elif byte_order == "little":
            vp = np.fromfile(file, dtype=np.dtype(dtype).newbyteorder("<"))
        else:
            raise ValueError("Please specify byte_order as either: little or big.")
        vp = vp.reshape(*axes, order=axes_order_sort)

        vp = np.flipud(vp.transpose((*axes_order,)))

        return vp, nz, nx, ny  # z, x and then y


############################################################
# main script #
############################################################
comm = MPI.COMM_WORLD
bbox = (-2600.0, 0.0, 0.0, 4800.0, 0.0, 4650.0)
bbox_buf = (
    -2600.0 - 250.0,
    -10.0,
    0.0 + 250.0,
    4800.0 - 250.0,
    0.0 + 250.0,
    4650.0 - 250.0,
)
hmin = 20
sigma = 0.5
fname = "irr_X193x25m_Y186x25m_Z105x25m.bin"
nz = 105
nx = 193
ny = 186
dx = dy = dz = 25

vp, nz, nx, ny = read_velocity_model(
    fname,
    nz=nz,
    nx=nx,
    ny=ny,
    byte_order="little",
    dtype="float32",
    axes_order=(2, 0, 1),
)
vp = np.squeeze(vp)

# create a signed distance field based on velocity Vp
phi = np.ones_like(vp)
phi[vp > 0.1] = -1
# call fast marching method
d = skfmm.distance(phi, [dz, dx, dy])
d = filters.gaussian_filter(d, (sigma), mode="constant")

z = np.linspace(bbox[0], bbox[1], nz)
x = np.linspace(bbox[2], bbox[3], nx)
y = np.linspace(bbox[4], bbox[5], ny)

interpolant = RegularGridInterpolator((z, x, y), d, bounds_error=False, fill_value=None)


def my_signed_distance_function(points):
    return interpolant(points)


ef = get_sizing_function_from_segy(
    fname,
    bbox,
    hmin=hmin,
    freq=15,
    wl=5,
    grade=0.15,
    hmax=1000,
    nz=nz,
    nx=nx,
    ny=ny,
    byte_order="little",  # NB: I'm passing the read_binary options here too!
    dtype="float32",
    axes_order=(2, 0, 1),
    axes_order_sort=None,
)
ef.bbox = bbox_buf

points, cells = generate_mesh(
    domain=my_signed_distance_function,
    bbox=bbox_buf,
    edge_length=ef,
    max_iter=35,
    verbose=2,
)

# bound the minimum dihedral angle:
# points, cells = sliver_removal(
#    points=points,
#    bbox=bbox,
#    domain=my_signed_distance_function,
#    edge_length=ef,
#    preserve=True,
#    max_iter=100,
#    verbose=2,
# )

# Meshes can be written quickly to disk using meshio and visualized with ParaView::
if comm.rank == 0:
    meshio.write_points_cells(
        "irr3d.vtk",
        points[:, [1, 2, 0]] / 1000.0,
        [("tetra", cells)],
    )

[krober10nd]

@handechao I think I have a good solution. In order to capture the irregular free surface, but have quality boundaries on the sides, I compute the intersection with a cube. This can also be used to add in your absorbing boundary layer if you wish.

See the code below and the image:

import meshio
import numpy as np
import skfmm
from mpi4py import MPI
from scipy.interpolate import RegularGridInterpolator
from scipy.ndimage import filters

from SeismicMesh import (Cube, generate_mesh, get_sizing_function_from_segy,
                         sliver_removal)


def read_velocity_model(
    filename,
    nz=None,
    nx=None,
    ny=None,
    byte_order=None,
    axes_order=None,
    axes_order_sort=None,
    dtype=None,
):
    """Read a velocity model"""
    return read_bin(
        filename, nz, nx, ny, byte_order, axes_order, axes_order_sort, dtype
    )


def read_bin(filename, nz, nx, ny, byte_order, axes_order, axes_order_sort, dtype):
    """Read a velocity model from a binary"""
    if (nz is None) or (nx is None) or (ny is None):
        raise ValueError(
            "Please specify the number of grid points in each dimension (e.g., `nz`, `nx`, `ny`)..."
        )
    axes = [nz, nx, ny]
    ix = np.argsort(axes_order)
    axes = [axes[o] for o in ix]
    with open(filename, "r") as file:
        print(f"Reading binary file: {filename}")
        if byte_order == "big":
            vp = np.fromfile(file, dtype=np.dtype(dtype).newbyteorder(">"))
        elif byte_order == "little":
            vp = np.fromfile(file, dtype=np.dtype(dtype).newbyteorder("<"))
        else:
            raise ValueError("Please specify byte_order as either: little or big.")
        vp = vp.reshape(*axes, order=axes_order_sort)

        vp = np.flipud(vp.transpose((*axes_order,)))

        return vp, nz, nx, ny  # z, x and then y


############################################################
# main script #
############################################################
comm = MPI.COMM_WORLD
bbox = (-2600.0, 0.0, 0.0, 4800.0, 0.0, 4650.0)
hmin = 20
sigma = 0.5
fname = "irr_X193x25m_Y186x25m_Z105x25m.bin"
nz = 105
nx = 193
ny = 186
dx = dy = dz = 25

vp, nz, nx, ny = read_velocity_model(
    fname,
    nz=nz,
    nx=nx,
    ny=ny,
    byte_order="little",
    dtype="float32",
    axes_order=(2, 0, 1),
)
vp = np.squeeze(vp)

# create a signed distance field based on velocity Vp
phi = np.ones_like(vp)
phi[vp > 0.1] = -1
# call fast marching method
d = skfmm.distance(phi, [dz, dx, dy])
d = filters.gaussian_filter(d, (sigma), mode="constant")

z = np.linspace(bbox[0], bbox[1], nz)
x = np.linspace(bbox[2], bbox[3], nx)
y = np.linspace(bbox[4], bbox[5], ny)

interpolant = RegularGridInterpolator((z, x, y), d, bounds_error=False, fill_value=None)

cube = Cube((-2600.0, 200.0, 0.0, 4800.0, 0.0, 4650.0))


def my_signed_distance_function(points):
    d = [interpolant(points), cube.eval(points)]
    return np.maximum.reduce(d)


ef = get_sizing_function_from_segy(
    fname,
    bbox,
    hmin=hmin,
    freq=15,
    wl=5,
    grade=0.15,
    hmax=1000,
    nz=nz,
    nx=nx,
    ny=ny,
    byte_order="little",  # NB: I'm passing the read_binary options here too!
    dtype="float32",
    axes_order=(2, 0, 1),
    axes_order_sort=None,
)

points, cells = generate_mesh(
    domain=my_signed_distance_function,
    bbox=bbox,
    edge_length=ef,
    max_iter=100,
    verbose=2,
)

# bound the minimum dihedral angle:
points, cells = sliver_removal(
    points=points,
    bbox=bbox,
    domain=my_signed_distance_function,
    edge_length=ef,
    preserve=True,
    max_iter=100,
    verbose=2,
    min_dh_angle_bound=5.0, 
)

# Meshes can be written quickly to disk using meshio and visualized with ParaView::
if comm.rank == 0:
    meshio.write_points_cells(
        "irr3d.vtk",
        points[:, [1, 2, 0]] / 1000.0,
        [("tetra", cells)],
    )

[handechao]

Thank you, great result. I will try to understand this code and test it on other models next week.

[krober10nd]

No problem. One thing we will have to experiment with is the sliver removal algorithm here. It's not always trivial to remove all of them and you have to try alternative strategies.