This notebooks shows the most basic mesh with all necessary commands to put material properties and boundaries ready to run with Salvus.
# initialize notebook
%matplotlib inline
from __future__ import print_function
import matplotlib.pyplot as plt
import numpy as np
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = (10, 8)
Set Basic input parameters for a rectangular domain with homogeneous material properties.
vp = 1700 # p-wave velocity in m/s
vs = 1000 # s-wave velocity in m/s
rho = 2700 # density in kg / m^3
max_x = 10000 # domain size in x direction
max_y = 5000 # domain size in y direction
fmax = 1. # maximum frequency in Hz
elements_per_wavelength = 2. # resolution criterion
Compute element size and number of elements needed according to the resolution criterion.
hmax = vs / fmax / elements_per_wavelength
nelem_x = int(max_x / hmax)
nelem_y = int(max_y / hmax)
Build a structured Grid and plot it.
from pymesher.structured_grid_2D import StructuredGrid2D
sg = StructuredGrid2D.rectangle(nelem_x=nelem_x, nelem_y=nelem_y, max_x=max_x, max_y=max_y)
sg.plot(show=True)
Convert the structured grid to an unstructured mesh (trivial in this special case, but uses the same general routines).
m = sg.get_unstructured_mesh()
Attach the material properties, constant here, but using the general array storage on element nodes.
m.attach_field('VP', vp * np.ones(m.npoint))
m.attach_field('VS', vs * np.ones(m.npoint))
m.attach_field('RHO', rho * np.ones(m.npoint))
# this is currently necessary to tell salvus that this is elastic material everywhere. Set to 1 to make it fluid.
m.attach_field('fluid', np.zeros(m.nelem))
Find outer Boundaries of the domain assuming the domain is rectangular.
m.find_side_sets(mode='cartesian')
Write to an exodus file that can be read by paraview and salvus.
m.write_exodus('basic2D.e')
Open with paraview. You propably need to enable the advanced mode clicking on the gear symbol to see the boundaries.
!~/ParaView/bin/paraview basic2D.e
Very much the same as in 2D, so less comments.
vp = 1700 # p-wave velocity in m/s
vs = 1000 # s-wave velocity in m/s
rho = 2700 # density in kg / m^3
max_x = 10000 # domain size in x direction
max_y = 5000 # domain size in y direction
max_z = 3000 # domain size in y direction
fmax = 1. # maximum frequency in Hz
elements_per_wavelength = 2. # resolution criterion
hmax = vs / fmax / elements_per_wavelength
nelem_x = int(max_x / hmax)
nelem_y = int(max_y / hmax)
nelem_z = int(max_z / hmax)
Build a structured Grid. Don't plot it, because in 3D matplotlib is to slow.
from pymesher.structured_grid_3D import StructuredGrid3D
sg = StructuredGrid3D.cube(nelem_x=nelem_x, nelem_y=nelem_y, nelem_z=nelem_z,
max_x=max_x, max_y=max_y, max_z=max_z)
m = sg.get_unstructured_mesh()
m.attach_field('VP', vp * np.ones(m.npoint))
m.attach_field('VS', vs * np.ones(m.npoint))
m.attach_field('RHO', rho * np.ones(m.npoint))
m.find_side_sets(mode='cartesian')
m.write_exodus('basic3D.e')
!~/ParaView/bin/paraview basic3D.e