import numpy as np
import xarray as xr
import salvus.namespace as sn
vp = xr.DataArray(
np.r_[np.full(11, 1.0), np.linspace(2.0, 7.5, 51)],
[("y", np.linspace(1.0, 0.0, 62))],
)
vs = xr.DataArray(
np.r_[np.full(11, 0.5), np.linspace(1.0, 3.5, 51)],
[("y", np.linspace(1.0, 0.0, 62))],
)
rho = xr.DataArray(
np.r_[np.full(11, 1.0), np.linspace(2.0, 7.5, 51)],
[("y", np.linspace(1.0, 0.0, 62))],
)
vs.plot()
[<matplotlib.lines.Line2D at 0x7af55fc5f450>]
xc = np.linspace(0, 1, 101)
i0 = sn.layered_meshing.interface.Curve.from_points(
xc, np.sin(2 * np.pi * xc) * 0.05 - 0.05, reference_elevation=1.0, axis="x"
)
d = sn.domain.dim2.BoxDomain(x0=0, x1=1, y0=0, y1=1)
mr = sn.MeshResolution(reference_frequency=10.0, elements_per_wavelength=2.0)
m = sn.layered_meshing.LayeredModel(
[i0, sn.material.from_params(rho=rho, vp=vp, vs=vs)]
)
sn.layered_meshing.mesh_from_domain(domain=d, model=m, mesh_resolution=mr)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7af55b1774d0>
m_p1 = sn.layered_meshing.partition(
d, m, model_predicate=lambda x: x.ds.VS == 0.5, mode="discontinuous"
)
m_p1.n_layers
2
sn.layered_meshing.mesh_from_domain(domain=d, model=m_p1, mesh_resolution=mr)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7af55b021590>
"filter_doubling_monotonic_top_down"
. This will identify regions in the
model where doubling can safely occur. We choose a "continuous" partitioning
mode here as there is no true discontinuity we're looking for, simply a
logical one that will allow the mesh to add a doubling layer.m_p2 = sn.layered_meshing.partition(
d,
m_p1,
model_predicate=(
sn.layered_meshing.filters.filter_doubling_monotonic_top_down
),
mode="continuous",
)
m_p2.n_layers
3
sn.layered_meshing.mesh_from_domain(domain=d, model=m_p2, mesh_resolution=mr)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7af552de0790>
sn.layered_meshing.mesh_from_domain(
domain=d,
model=sn.layered_meshing.MeshingProtocol(
m_p2,
intralayer_coarsening_policy=[
# Allow for a variable number of elements in the vertical
# direction. doubling: Place refinements at the top of the layer.
sn.layered_meshing.meshing_protocol.IntralayerVerticalRefine(
refinement_type="doubling"
),
# Enforce a constant element size.
sn.layered_meshing.meshing_protocol.IntralayerConstant(),
],
),
mesh_resolution=mr,
)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7af551b724d0>
InterlayerDoubling()
. We can change
this by explicitly setting a different policy, i.e. InterlayerConstant()
as
below. One policy is allowed per internal interface, and the last policy will
be applied to all subsequent interfaces. Let's quickly inspect the effect of
changing the policy to InterlayerConstant()
below, a more thorough look at
the different options will be explored in a future tutorial.sn.layered_meshing.mesh_from_domain(
domain=d,
model=sn.layered_meshing.MeshingProtocol(
m_p2,
intralayer_coarsening_policy=[
# Allow for a variable number of elements in the vertical
# direction. doubling: Place refinements at the top of the layer.
sn.layered_meshing.meshing_protocol.IntralayerVerticalRefine(
refinement_type="doubling"
),
# Enforce a constant element size.
sn.layered_meshing.meshing_protocol.IntralayerConstant(),
],
# Enforce a constant number of elements in the horizontal direction
# across layer boundaries.
interlayer_coarsening_policy=(
sn.layered_meshing.meshing_protocol.InterlayerConstant()
),
),
mesh_resolution=mr,
)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7af552d18a50>