assemble(gf₁₁::AbstractArray, gf₁₂::AbstractMatrix, gf₂₁::AbstractMatrix, gf₂₂::AbstractMatrix,
         pf::RateStateQuasiDynamicProperty, pa::ViscosityProperty, u0::ArrayPartition, tspan::NTuple{2};
         se::StateEvolutionLaw=DieterichStateLaw(), kwargs...)

Construct an ODEProblem for viscoelastic rate-and-state friction model with quasi-dynamic evolution.

Arguments

• gf₁₁: Greens function array for fault-fault interaction
• gf₁₂: Greens function array for fault-mantle interaction
• gf₂₁: Greens function array for mantle-fault interaction
• gf₂₂: Greens function array for mantle-mantle interaction
• pf: rate-and-state quasi-dynamic property for fault
• pa: viscosity property for mantle
• u0: initial state partition, must be an ArrayPartition with 5 components: velocity, state variable, strain, stress, and fault slip
• tspan: time span for the simulation, a tuple of two values (start, stop)
• se: state evolution law, defaults to DieterichStateLaw()

Returns

• An ODEProblem object that can be solved using OrdinaryDiffEq.jl

assemble(gf::AbstractArray, p::RateStateQuasiDynamicProperty, u0::ArrayPartition, tspan::NTuple{2};
         se::StateEvolutionLaw=DieterichStateLaw(), kwargs...)

Construct an ODEProblem for rate-and-state friction model with quasi-dynamic evolution.

Arguments

• gf: Greens function array for the fault
• p: rate-and-state quasi-dynamic property
• u0: initial state partition, must be an ArrayPartition with 3 components: velocity, state variable, and fault slip
• tspan: time span for the simulation, a tuple of two values (start, stop)
• se: state evolution law, defaults to DieterichStateLaw()

Returns

• An ODEProblem object that can be solved using OrdinaryDiffEq.jl.

assemble(gf::AbstractArray, p::RateStateQuasiDynamicProperty, dila::DilatancyProperty, u0::ArrayPartition, tspan::NTuple{2};
         se::StateEvolutionLaw=DieterichStateLaw(), kwargs...)

Construct an ODEProblem for rate-and-state friction model with dilatancy and quasi-dynamic evolution.

Arguments

• gf: Greens function array for the fault
• p: rate-and-state quasi-dynamic property
• dila: dilatancy property
• u0: initial state partition, must be an ArrayPartition with 4 components: velocity, state variable, pressure, and fault slip
• tspan: time span for the simulation, a tuple of two values (start, stop)
• se: state evolution law, defaults to DieterichStateLaw()

Returns

• An ODEProblem object that can be solved using OrdinaryDiffEq.jl

gen_gmsh_mesh(::Val{:BEMHex8Mesh},
    llx::T, lly::T, llz::T, dx::T, dy::T, dz::T, nx::I, ny::I, nz::I;
    rfx::T=one(T), rfy::T=one(T), rfzh::AbstractVector=ones(nz),
    rfxType::AbstractString="Bump", rfyType::AbstractString="Bump",
    output::AbstractString="temp.msh"
) where {T, I}

Gernate a box using 8-node hexahedron elements by vertically extruding transfinite curve on xy plane, allowing total flexibility on the mesh size in z direction, and refinement in xy plane.

Arguments

• llx, lly, llz: coordinates of low-left corner on the top surface
• dx, dy, dz: x-, y-, z-extension
• nx, ny: number of cells along x-, y-axis
• rfx, rfy: refinement coefficients along x-, y-axis using Bump algorithm, please refer gmsh.model.geo.mesh.setTransfiniteCurve
• rfzh: accumulated height of cells along z-axis which will be normalized automatically, please refer heights in gmsh.model.geo.extrude

Returns

• A GMSH mesh file of the BEMHex8Mesh` for the mantle deformation modeling.

gen_mesh(::Val{:BEMHex8Mesh}, file::AbstractString;
    rotation::Real=0.0, transpose::Bool=false)

Load a BEM Hex8 mesh from a GMSH file.

Arguments

• file::AbstractString: path to the GMSH mesh file
• rotation::Real=0.0: rotation angle in degrees to apply to the mesh
• transpose::Bool=false: if true, transpose the mesh coordinates

Returns

• A BEMHex8Mesh object representing the mantle mesh for modeling.

gen_mesh(::Val{:RectOkada},
    x::T, ξ::T, Δx::T, Δξ::T, dip::T) where T

Generate a rectangular mesh for Okada's fault model in 2D.

Arguments

• x: length of the fault along strike
• ξ: length of the fault along downdip
• Δx: cell size along strike
• Δξ: cell size along downdip
• dip: dipping angle of the fault in degrees

Returns

• A RectOkadaMesh object representing the transfinite mesh suitable for fault modeling.

gen_pvd(mf::RectOkadaMesh, mafile, solh5, tstr, ufstrs, uastrs, steps, output;
    tscale=365*86400, mafiletype=Val(:BEMHex8Mesh))

This function generates a Paraview collection file from the Okada mesh (fault plane) and the volume mesh, and solution data stored in an HDF5 file.

Arguments

• mf::RectOkadaMesh: the rectangular Okada mesh object
• mafile: the mesh file for the BEM model
• solh5: the HDF5 file containing the solution data
• tstr: the name of the time data in the HDF5 file
• ufstrs: a vector of strings representing the names of the solution components (fault) in the HDF5 file
• uastrs: a vector of strings representing the names of the solution components (mantle) in the HDF5 file
• steps: a vector of integers representing the time steps to be included in the output
• output: the output directory for the Paraview collection file
• tscale: a scaling factor for the time data (default is 365*86400 seconds, which is one year)
• mafiletype: the type of the mesh file, default is :BEMHex8Mesh, which indicates a BEM mesh with Hex8 elements

Returns

• A Paraview collection file containing the mesh and solution data at specified time steps.

gen_vtk_grid(mesh::RectOkadaMesh)

This function creates a grid suitable for visualization in Paraview or similar tools.

Arguments

• mesh::RectOkadaMesh: the rectangular Okada mesh object

Returns

• A tuple containing the x, y, and z coordinates of the grid points.

gen_vtk_grid(t::Val{:BEMHex8Mesh}, mesh)

This function uses GMSH to read the Hex8 mesh and convert it into a VTK grid format for visualization in Paraview or similar tools.

Arguments

• t::Val{:BEMHex8Mesh}: a value type indicating the mesh type
• mesh: the mesh object to be converted

Returns

• A tuple containing the points and cells of the mesh in VTK format.

gmsh_quadrature(etype::Integer, qtype::String)::QuadratureType

This function retrieves the integration points and weights for a given element type and quadrature type from GMSH.

Arguments

• etype::Integer: the element type ID in GMSH
• qtype::String: the quadrature type, e.g., "Gauss1"

Returns

• A tuple containing the local coordinates and weights for the quadrature points.

stress_greens_function(ma::BEMHex8Mesh, mf::RectOkadaMesh, λ::T, μ::T;
    ftype::FaultType=StrikeSlip(),
) where {T}

This function computes the stress Green's function for a source in Hex8 mesh to a receiver on Okada mesh.

Arguments

• ma::BEMHex8Mesh: the BEM Hex8 mesh object
• mf::RectOkadaMesh: the rectangular Okada mesh object
• λ::T: Lamé's first parameter
• μ::T: shear modulus
• ftype::FaultType=StrikeSlip(): type of fault, either StrikeSlip or DipSlip

Returns

• A 2D array of stress Green's functions, where each column corresponds to a source

mantle patch and each row corresponds to a receiver fault patch.

stress_greens_function(
    mesh::BEMHex8Mesh,
    λ::T, μ::T;
    qtype::Union{String, Tuple{AbstractVecOrMat, AbstractVector}}="Gauss1",
    checkeigvals::Bool=true,
) where {T}

This function computes the stress Green's function for a BEM Hex8 mesh.

Arguments

• mesh::BEMHex8Mesh: the BEM Hex8 mesh object
• λ::T: Lamé's first parameter
• μ::T: shear modulus
• qtype::Union{String, Tuple{AbstractVecOrMat, AbstractVector}} quadrature type for integration, can be a string or a tuple of local coordinates and weights
• checkeigvals::Bool=true: if true, check the eigenvalues of the resulting stress matrix and print the maximum real part

Returns

• A 2D array of stress Green's functions, where each column corresponds to a source

mantle patch and each row corresponds to another mantle patch.

stress_greens_function(
    mf::RectOkadaMesh, ma::BEMHex8Mesh,
    λ::T, μ::T;
    ftype::FaultType=StrikeSlip(),
    qtype::Union{String, Tuple{AbstractVecOrMat, AbstractVector}}="Gauss1",
    nrept::Integer=2, buffer_ratio::Real=0,
) where {T}

This function computes the stress Green's function for a source in Okada mesh to a receiver on Hex8 mesh.

Arguments

• mf::RectOkadaMesh: the rectangular Okada mesh object
• ma::BEMHex8Mesh: the BEM Hex8 mesh object
• λ::T: Lamé's first parameter
• μ::T: shear modulus
• ftype::FaultType=StrikeSlip(): type of fault, either StrikeSlip or DipSlip
• qtype::Union{String, Tuple{AbstractVecOrMat, AbstractVector}}: quadrature type for integration, can be a string or a tuple of local coordinates and weights
• nrept::Integer=2: number of repetitions for the dislocation
• buffer_ratio::Real=0: ratio of buffer zone around the mesh

Returns

• A 2D array of stress Green's functions, where each column corresponds to a source fault patch

and each row corresponds to a receiver mantle patch.

stress_greens_function(mesh::RectOkadaMesh, λ::T, μ::T;
    ftype::FaultType=StrikeSlip(), fourier::Bool=true,
    nrept::Integer=2, buffer_ratio::Real=0, fftw_flags::UInt32=FFTW.PATIENT
) where {T <: Real}

This function computes the stress Green's function for a rectangular Okada mesh fault.

Arguments

• mesh::RectOkadaMesh: the rectangular Okada mesh object
• λ::T: Lamé's first parameter
• μ::T: shear modulus
• ftype::FaultType=StrikeSlip(): type of fault, either StrikeSlip or DipSlip
• fourier::Bool=true: if true, return the Fourier transform of the Green's function
• nrept::Integer=2: number of repetitions for the dislocation
• buffer_ratio::Real=0: ratio of buffer zone around the mesh
• fftw_flags::UInt32=FFTW.PATIENT: flags for FFTW plan

Returns

• If fourier is true, returns the Fourier transform of the stress Green's function, otherwise, returns the stress Green's function as a 3D array. The array considers the translational symmetry of the transfinite mesh.

wsolve(prob::ODEProblem, alg::OrdinaryDiffEqAlgorithm,
    file, nstep, getu, ustrs, tstr; kwargs...)

Write the solution to HDF5 file while solving the ODE. The interface is exactly the same as solve an ODEProblem except a few more about the saving procedure. Notice, it will set save_everystep=false so to avoid memory blow up. The return code will be written as an attribute in tstr data group.

Extra Arguments

• file::AbstractString: name of file to be saved
• nstep::Integer: number of steps after which a saving operation will be performed
• getu::Function: function handler to extract desired solution for saving, which should have the signature getu(u, t, integrator), where u is the current solution, t is the current time, and integrator is the current integrator object. The output should be a tuple of arrays or vectors to be saved.
• ustr::AbstractVector: list of names to be assigned for each components, whose length must equal the length of getu output
• tstr::AbstractString: name of time data

KWARGS

• stride::Integer=1: downsampling rate for saving outputs
• append::Bool=false: if true then append solution after the end of file
• force::Bool=false: force to overwrite the existing solution file

Returns

• sol::ODESolution: the solution object of OrdinaryDiffEq.jl