import numpy as np
import xarray as xr
import salvus.namespace as sn
sn.layered_meshing.mesh_from_domain(
domain=sn.domain.dim2.BoxDomain(x0=0.0, x1=1.0, y0=0.0, y1=1.0),
model=sn.material.elastic.Velocity.from_params(rho=1.0, vp=1.0, vs=0.5),
mesh_resolution=sn.MeshResolution(reference_frequency=10.0),
)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7a4d1a04da90>
# We cant actually access any material much faster by letting salvus figure
# out which one we want by the arguments supplied
# m0 = sn.material.elastic.Velocity.from_params(rho=1.0, vp=1.0, vs=0.5)
m0 = sn.material.from_params(rho=1.0, vp=1.0, vs=0.5)
xarray
Dataset,
which can be queried by the materials .ds
variable.m0.ds
<xarray.Dataset> Size: 24B Dimensions: () Data variables: RHO float64 8B 1.0 VP float64 8B 1.0 VS float64 8B 0.5
<xarray.Dataset> Size: 24B Dimensions: () Data variables: RHO float64 8B 1.0 VP float64 8B 1.0 VS float64 8B 0.5
xarray
representation of the material
parameters we specified. You might notice that no coordinates are specified
here. This is no problem, as the coordinates will be "realized" when the
model is eventually used in the meshing process. Realization in this context
refers to ensuring that all coordinates present in a layered model are
consistent with the spatial dimension and domain they are being used within.
As no coordinates are present here, the realization stage will thus ensure
that the model values are propagated to all locations in the host domain. In
the future we will see how realization can be exploited to define models on
only a subset of a domain's coordinate axes, and additionally define them
relative to either the domain bounds, or to the bounds of a host layer.m0_lame = sn.material.from_params(lam=0.5, mu=0.25, rho=1.0)
print(sn.material.elastic.Velocity.from_material(m0_lame))
Velocity(RHO=ConstantParameter(p=1.0), VP=ConstantParameter(p=1.0), VS=ConstantParameter(p=0.5))
.to_wavelength_oracle()
that will automatically compute the critical
parameter required to determine the required mesh size. For isotropic elastic
models this is simply the model's shear wave velocity, but for less symmetric
materials the oracle may instead be a combination of the material's
parameters.print(f"{m0.to_wavelength_oracle()}\n{m0_lame.to_wavelength_oracle()}")
ConstantParameter(p=0.5) ConstantParameter(p=0.5)
m0_lame
in
terms of , , and , Salvus automatically computed it for us.
The presence of this oracle for each parameterization means that the tedious
process of manually computing the minimum wavelength is a thing of the past
-- Salvus will take care of this for you.LayeredModel
to tell Salvus that we are ready to proceed to the next stage
of mesh generation.lm_0 = sn.layered_meshing.LayeredModel(m0)
LayeredModel
object has some interesting properties. First off, we can
query the models it contains and see our material stored therein. Note that
the models are presented as a list with one element; again, the utility of
this will be apparent later on.lm_0.models
[Velocity(RHO=ConstantParameter(p=1.0), VP=ConstantParameter(p=1.0), VS=ConstantParameter(p=0.5))]
.interfaces
property. Here, we see an empty list.lm_0.interfaces
[]
.complete()
method on
the layered model and inspect the interfaces added.print("\n\n".join(str(s) for s in lm_0.complete().interfaces))
Hyperplane(da=<xarray.DataArray ()> Size: 8B array(0.) Attributes: reference_elevation: Depth(value=0.0), extender=<cyfunction extrude_like_and_pad at 0x7a4d1bd79560>, interpolation_method='linear') Hyperplane(da=<xarray.DataArray ()> Size: 8B array(0.) Attributes: reference_elevation: Height(value=0.0), extender=<cyfunction extrude_like_and_pad at 0x7a4d1bd79560>, interpolation_method='linear')
strata
.print("\n\n".join(str(s) for s in lm_0.complete().strata))
Hyperplane(da=<xarray.DataArray ()> Size: 8B array(0.) Attributes: reference_elevation: Depth(value=0.0), extender=<cyfunction extrude_like_and_pad at 0x7a4d1bd79560>, interpolation_method='linear') Velocity(RHO=ConstantParameter(p=1.0), VP=ConstantParameter(p=1.0), VS=ConstantParameter(p=0.5)) Hyperplane(da=<xarray.DataArray ()> Size: 8B array(0.) Attributes: reference_elevation: Height(value=0.0), extender=<cyfunction extrude_like_and_pad at 0x7a4d1bd79560>, interpolation_method='linear')
i0, i1 = lm_0.complete().interfaces
print(i0.da.reference_elevation)
print(i1.da.reference_elevation)
Depth(value=0.0) Height(value=0.0)
Domain
object. For now, let's
just use a simple 2-D Box domain.d_2d = sn.domain.dim2.BoxDomain(x0=0, x1=1, y0=0, y1=1)
d_2d.plot()
MeshResolution
object.
Here we can set the desired elements per wavelength, reference frequency, and
order of model representation as follows:mr = sn.MeshResolution(
reference_frequency=2.0, elements_per_wavelength=1.5, model_order=2
)
mesh_from_domain
function as follows:mesh = sn.layered_meshing.mesh_from_domain(
domain=d_2d, model=lm_0, mesh_resolution=mr
)
mesh