API

PoreData is a collection of physical parameters for coupled geomechanics and flow simulation

• M: Biot modulus
• b: Biot coefficient
• ρb: Bulk density
• ρf: Fluid density
• kp: Permeability
• E: Young modulus
• ν: Poisson ratio
• μ: Fluid viscosity
• Pi: Initial pressure
• Bf: formation volume, $B_f=\frac{\rho_{f,0}}{\rho_f}$
• g: Gravity acceleration

Mesh holds data structures for an unstructured mesh.

• nodes: a $n_v \times 2$ coordinates array
• edges: a $n_{\text{edge}} \times 2$ integer array for edges
• elems: a $n_e \times 3$ connectivity matrix, 1-based.
• nnode, nedge, nelem: number of nodes, edges, and elements
• ndof: total number of degrees of freedoms
• conn: connectivity matrix, nelems × 3 or nelems × 6, depending on whether a linear element or a quadratic element is used.
• lorder: order of quadrature rule for line integrals (default = 6, 4 gauss points per line segment)
• elem_type: type of the element (P1, P2 or BDM1)

Constructors

Mesh(m::Int64, n::Int64, h::Float64; order::Int64 = -1, 
            degree::Union{FiniteElementType, Int64} = 1, lorder::Int64 = -1, 
            version::Int64 = 1)

Constructs a mesh of a rectangular domain. The rectangle is split into $m\times n$ cells, and each cell is further split into two triangles. order specifies the quadrature rule order. degree determines the degree for finite element basis functions.

Integration Order : order and lorder

TODO

CrackMesh(m::Int64, n::Int64, h::Float64, k::Int64 = 1)

Creates a crack mesh.

mmesh = CrackMesh(20, 10, 0.1, 4)
visualize_mesh(mmesh)

To access the underlying Mesh object, use mmesh.mesh.

compute_fem_stiffness_matrix(K::Array{Float64,2}, m::Int64, n::Int64, h::Float64)

Computes the term

\[\int_{A}\delta \varepsilon :\sigma\mathrm{d}x = \int_A u_AB^TKB\delta u_A\mathrm{d}x\]

where the constitutive relation is given by

\[\begin{bmatrix}\sigma_{xx}\\\sigma_{yy}\\\sigma_{xy}\end{bmatrix} = K \begin{bmatrix}\varepsilon_{xx}\\\varepsilon_{yy}\\2\varepsilon_{xy}\end{bmatrix}\]

compute_fem_stiffness_matrix(hmat::PyObject,m::Int64, n::Int64, h::Float64)

A differentiable kernel. hmat has one of the following sizes

• \[3\times 3\]

• \[4mn \times 3 \times 3\]

compute_fem_stiffness_matrix(kappa::PyObject, mesh::Mesh)
compute_fem_stiffness_matrix(kappa::Array{Float64, 3}, mesh::Mesh)
compute_fem_stiffness_matrix(kappa::Array{Float64, 2}, mesh::Mesh)

Computes the stiffness matrix. Here, the acceptable sizes of $\kappa$

• a $3\times 3$ matrix, which is the pointwise uniform stiffness matrix
• a $N_g \times 3 \times 3$ tensor, which includes a specific stiffness matrix at each Gauss node

compute_interaction_matrix(m::Int64, n::Int64, h::Float64)

Computes the interaction term

\[\int_A p \delta \varepsilon_v\mathrm{d}x = \int_A p [1,1,0]B^T\delta u_A\mathrm{d}x\]

Here $\varepsilon_v = \text{tr}\; \varepsilon = \text{div}\; \mathbf{u}$.

The output is a $mn \times 2(m+1)(n+1)$ matrix.

compute_interaction_matrix(mesh::Mesh)

compute_fvm_tpfa_matrix(m::Int64, n::Int64, h::Float64)

Computes the term with two-point flux approximation

\[\int_{A_i} \Delta p \mathrm{d}x = \sum_{j=1}^{n_{\mathrm{faces}}} (p_j-p_i)\]

compute_fvm_tpfa_matrix(K::Array{Float64}, m::Int64, n::Int64, h::Float64)

Computes the term with two-point flux approximation with distinct permeability at each cell

\[\int_{A_i} K_i \Delta p \mathrm{d}x = K_i\sum_{j=1}^{n_{\mathrm{faces}}} (p_j-p_i)\]

Note that $K$ is a length $mn$ vector, representing values per cell.

compute_fvm_tpfa_matrix(K::Array{Float64}, bc::Array{Int64,2}, pval::Array{Float64,1}, m::Int64, n::Int64, h::Float64)

Computes the term with two-point flux approximation with distinct permeability at each cell

\[\int_{A_i} K_i \Delta p \mathrm{d}x = K_i\sum_{j=1}^{n_{\mathrm{faces}}} (p_j-p_i)\]

Here $K$ is a length $mn$ vector, representing values per cell.

Additionally, Dirichlet boundary conditions are imposed on the boundary edges bc (a $N\times 2$ integer matrix), i.e., the $i$-th edge has value pval. The ghost node method is used for imposing the Dirichlet boundary condition. The other boundaries are no-blow boundaries, i.e., $\frac{\partial T}{\partial n} = 0$. The function outputs a length $mn$ vector and $mn\times mn$ matrix $M$.

\[\int_{A_i} K_i \Delta p \mathrm{d}x = f_i + M_{i,:}\mathbf{p}\]

Returns both the sparse matrix A and the right hand side rhs.

compute_fvm_tpfa_matrix(K::PyObject, bc::Array{Int64,2}, pval::Union{Array{Float64},PyObject}, m::Int64, n::Int64, h::Float64)

A differentiable kernel for compute_fvm_tpfa_matrix.

compute_fvm_tpfa_matrix(K::PyObject, m::Int64, n::Int64, h::Float64)

A differentiable kernel for compute_fvm_tpfa_matrix.

compute_fem_mass_matrix(m::Int64, n::Int64, h::Float64)

Computes the finite element mass matrix

\[\int_{\Omega} u \delta u \mathrm{d}x\]

The matrix size is $2(m+1)(n+1) \times 2(m+1)(n+1)$.

compute_fem_mass_matrix(ρ::PyObject,m::Int64, n::Int64, h::Float64)

A differentiable kernel. $\rho$ is a vector of length $4mn$ or $8mn$

compute_fvm_mass_matrix(m::Int64, n::Int64, h::Float64)

Returns the FVM mass matrix

\[\int_{A_i} p_i \mathrm{d}x = h^2 p_i \]

compute_fem_mass_matrix1(ρ::Array{Float64,1}, m::Int64, n::Int64, h::Float64)

Computes the mass matrix for a scalar value $u$

\[\int_A \rho u \delta u \mathrm{d} x\]

The output is a $(m+1)(n+1)\times (m+1)(n+1)$ sparse matrix.

compute_fem_mass_matrix1(m::Int64, n::Int64, h::Float64)

Computes the mass matrix for a scalar value $u$

\[\int_A u \delta u \mathrm{d} x\]

The output is a $(m+1)(n+1)\times (m+1)(n+1)$ sparse matrix.

compute_fem_mass_matrix1(ρ::PyObject,m::Int64, n::Int64, h::Float64)

A differentiable kernel.

compute_fem_mass_matrix1(rho::Union{PyObject, Array{Float64, 1}}, mesh::Mesh)

compute_fem_mass_matrix1(ρ::Union{PyObject, Array{Float64, 1}}, 
        mmesh::Mesh3)

compute_fem_mass_matrix1(mmesh::Mesh3)

compute_fem_stiffness_matrix1(K::Array{Float64,2}, m::Int64, n::Int64, h::Float64)

Computes the term

\[\int_{A} (K \nabla u) \cdot \nabla \delta u \mathrm{d}x = \int_A u_A B^T K B \delta u_A\mathrm{d}x\]

Returns a $(m+1)\times (n+1)$ matrix

compute_fem_stiffness_matrix1(hmat::PyObject, m::Int64, n::Int64, h::Float64)

A differentiable kernel for computing the stiffness matrix. Two possible shapes for hmat are supported:

• \[4mn \times 2\times 2\]

• \[2 \times 2\]

compute_fvm_advection_matrix(v::Union{PyObject, Array{Float64, 2}},
    bc::Array{Int64, 2},bcval::Union{PyObject, Array{Float64}},m::Int64,n::Int64,h::Float64)

Computes the advection matrix for use in the implicit scheme

\[\int_A \mathbf{v} \cdot \nabla u dx \]

Here v is a $2mn$ vector, where the first $mn$ entries corresponds to the first dimension of $\mathbf{v}$ and the remaining $mn$ entries corresponds to the second dimension.

It returns a matrix $mn\times mn$ matrix $K$ and an auxilliary term $\mathbf{f}$ due to boundary conditions.

\[\int_\Omega \mathbf{v} \cdot \nabla u dx = K \mathbf{u} + \mathbf{f}\]

compute_fem_laplace_matrix1(K::Array{Float64, 1}, m::Int64, n::Int64, h::Float64)

Computes the coefficient matrix for

\[\int_\Omega K \nabla u \cdot \nabla (\delta u) \; dx \]

Here $K\in \mathbf{R}^{2\times 2}$, $u$ is a scalar variable, and K is a $4mn \times 2 \times 2$ matrix.

Returns a $(m+1)(n+1)\times (m+1)(n+1)$ sparse matrix.

compute_fem_laplace_matrix1(K::Array{Float64, 2}, m::Int64, n::Int64, h::Float64)

K is duplicated on each Gauss point.

compute_fem_laplace_matrix1(K::Array{Float64, 1}, m::Int64, n::Int64, h::Float64)

Computes the coefficient matrix for

\[\int_\Omega K\nabla u \cdot \nabla (\delta u) \; dx \]

Here K is a vector with length $4mn$ (defined on Gauss points).

Returns a $(m+1)(n+1)\times (m+1)(n+1)$ sparse matrix.

compute_fem_laplace_matrix1(m::Int64, n::Int64, h::Float64)

Computes the coefficient matrix for

\[\int_\Omega \nabla u \cdot \nabla (\delta u) \; dx \]

Returns a $(m+1)(n+1)\times (m+1)(n+1)$ sparse matrix.

compute_fem_laplace_matrix1(K::PyObject, m::Int64, n::Int64, h::Float64)

A differentiable kernel. Only $K\in \mathbb{R}^{4mn}$ is supported.

compute_fem_laplace_matrix1(kappa::PyObject, mesh::Mesh)

compute_fem_laplace_matrix1(kappa::Array{Float64,1}, mesh::Mesh)

compute_fem_laplace_matrix1(kappa::PyObject, mesh::Mesh3)

compute_fem_laplace_matrix1(kappa::Array{Float64,1}, mesh::Mesh3)

compute_fem_laplace_matrix(m::Int64, n::Int64, h::Float64)

Computes the coefficient matrix for

\[\int_\Omega \nabla \mathbf{u} \cdot \nabla (\delta \mathbf{u}) \; dx \]

Here

\[\mathbf{u} = \begin{bmatrix} u \\ v \end{bmatrix}\]

and

\[\nabla \mathbf{u} = \begin{bmatrix}u_x & u_y \\ v_x & v_y \end{bmatrix}\]

Returns a $2(m+1)(n+1)\times 2(m+1)(n+1)$ sparse matrix.

compute_fem_laplace_matrix(K::Array{Float64, 1}, m::Int64, n::Int64, h::Float64)

Computes the coefficient matrix for

\[\int_\Omega K \nabla \mathbf{u} \cdot \nabla (\delta \mathbf{u}) \; dx \]

Here $K$ is a scalar defined on Gauss points. K is a vector of length $4mn$

compute_fem_laplace_matrix(kappa::Union{PyObject, Array{Float64, 1}}, mesh::Mesh)

compute_fem_laplace_matrix(kappa::Union{PyObject, Array{Float64, 1}}, mesh::Mesh3)

compute_fem_advection_matrix1(u0::PyObject,v0::PyObject,m::Int64,n::Int64,h::Float64)

Computes the advection term for a scalar function $u$ defined on an FEM grid. The weak form is

\[\int_\Omega (\mathbf{u}_0 \cdot \nabla u) \delta u \; d\mathbf{x} = \int_\Omega \left(u_0 \frac{\partial u}{\partial x} \delta u + v_0 \frac{\partial u}{\partial y} \delta u\right)\; d\mathbf{x}\]

Here $u_0$ and $v_0$ are both vectors of length $4mn$.

Returns a sparse matrix of size $(m+1)(n+1)\times (m+1)(n+1)$

compute_fem_advection_matrix1(u::Union{Array{Float64,1}, PyObject},v::Union{Array{Float64,1}, PyObject}, mesh::Mesh)
compute_fem_advection_matrix1(u::Array{Float64,1}, v::Array{Float64,1}, mesh::Mesh)

compute_fem_bdm_mass_matrix(alpha::Union{Array{Float64,1}, PyObject},beta::Union{Array{Float64,1}, PyObject}, mmesh::Mesh)

Computes

\[\int_\Omega A\sigma : \delta \tau dx\]

Here

\[A\sigma = \alpha \sigma + \beta \text{tr} \sigma I\]

Here $\sigma$ and $\tau$ are both fourth-order tensors. The output is a 4mmesh.nedge × 4mmesh.nedge matrix.

compute_fem_bdm_mass_matrix(mmesh::Mesh)

Same as compute_fem_bdm_mass_matrix

compute_fem_bdm_mass_matrix(alpha::Array{Float64,1},beta::Array{Float64,1}, mmesh::Mesh)

compute_fem_bdm_mass_matrix1(alpha::Array{Float64,1}, mmesh::Mesh)

Computes

\[\int_\Omega \alpha\sigma \cdot \delta \tau dx\]

Here $\alpha$ is a scalar, and $\sigma$ and $\delta \tau$ are second order tensors.

The returned value is a 2mmesh.nedge × 2mmesh.nedge matrix.

compute_fem_bdm_mass_matrix1(mmesh::Mesh)

Same as compute_fem_bdm_mass_matrix1, except that $\alpha\equiv 1$

compute_fem_bdm_div_matrix(mmesh::Mesh)

Computes the coefficient matrix for

\[\int_\Omega \text{div} \tau \delta u dx\]

Here $\tau \in \mathbb{R}^{2\times 2}$ is a fourth-order tensor (not necessarily symmetric). mmesh uses the BDM1 finite element. The output is a 2mmesh.nelem × 4mmesh.nedge matrix.

compute_fem_bdm_div_matrix1(mmesh::Mesh)

Computes the coefficient matrix for

\[\int_\Omega \text{div} \tau \delta u dx\]

Here $\tau \in \mathbb{R}^{2}$ is a second-order tensor (not necessarily symmetric). mmesh uses the BDM1 finite element. The output is a mmesh.nelem × 2mmesh.nedge matrix.

compute_fem_bdm_skew_matrix(mmesh::Mesh)

Computes $\int_\Omega \sigma : v dx$ where

\[v = \begin{bmatrix}0 & \rho \\-\rho & 0 \end{bmatrix}\]

Here $\sigma$ is a fourth-order tensor.

The returned value is a mmesh.nelem × 4mmesh.nedge matrix.

compute_fem_boundary_mass_matrix1(c::Union{Array{Float64}, PyObject}, bdedge::Array{Int64, 2}, mmesh::Mesh)

Computes the matrix

\[\int_\Gamma cu \delta u ds\]

The parameters are

• bdedge: a $N_e \times 2$ integer array, the boundary edge to integrate on
• c: given by a vector of length $4N_e$; currently, each edge has 4 quadrature points;

The output is a $N_v\times N_v$ sparse matrix.

compute_fem_source_term(f1::Array{Float64}, f2::Array{Float64},
m::Int64, n::Int64, h::Float64)

Computes the term

\[\int_\Omega \mathbf{f}\cdot\delta u \mathrm{d}x\]

Returns a $2(m+1)(n+1)$ vector.

compute_fem_source_term(f1::PyObject, f2::PyObject,
m::Int64, n::Int64, h::Float64)

A differentiable kernel.

compute_fem_source_term(f1::Union{PyObject,Array{Float64,2}}, f2::Union{PyObject,Array{Float64,2}}, mesh::Mesh)

compute_fem_source_term(f1::Union{PyObject,Array{Float64,2}}, 
    f2::Union{PyObject,Array{Float64,2}}, mesh::Mesh3)

compute_fvm_source_term(f::Array{Float64}, m::Int64, n::Int64, h::Float64)

Computes the source term

\[\int_{A_i} f\mathrm{d}x\]

Here $f$ has length $4mn$ or $mn$. In the first case, an average value of four quadrature nodal values of $f$ is used per cell.

compute_fvm_source_term(f::PyObject, m::Int64, n::Int64, h::Float64)

A differentiable kernel.

compute_fvm_source_term(f::Array{Float64, 1}, mmesh::Mesh)

compute_fvm_mechanics_term(u::Array{Float64}, m::Int64, n::Int64, h::Float64)

Computes the mechanic interaction term

\[\int_{A_i} \varepsilon_v\mathrm{d}x\]

Here

\[\varepsilon_v = \mathrm{tr} \varepsilon = \varepsilon_{xx} + \varepsilon_{yy}\]

Numerically, we have

\[\varepsilon_v = [1 \ 1 \ 0] B^T \delta u_A\]

compute_fvm_mechanics_term(u::PyObject, m::Int64, n::Int64, h::Float64)

compute_fem_normal_traction_term(t::Array{Float64,1}, bdedge::Array{Int64},
m::Int64, n::Int64, h::Float64)
compute_fem_normal_traction_term(t::Float64, bdedge::Array{Int64},
m::Int64, n::Int64, h::Float64)

Computes the normal traction term

\[\int_{\Gamma} t(\mathbf{n})\cdot\delta u \mathrm{d}\]

Here $t(\mathbf{n})\parallel\mathbf{n}$ points outward to the domain and the magnitude is given by t. bdedge is a $N\times2$ matrix and each row denotes the indices of two endpoints of the boundary edge.

See compute_fem_traction_term for graphical illustration.

compute_fem_traction_term(t::Array{Float64, 2},
bdedge::Array{Int64,2}, m::Int64, n::Int64, h::Float64)

Computes the traction term

\[\int_{\Gamma} t(\mathbf{n})\cdot\delta u \mathrm{d}\]

The number of rows of t is equal to the number of edges in bdedge. The first component of t describes the $x$ direction traction, while the second component of t describes the $y$ direction traction.

Also see compute_fem_normal_traction_term.

compute_fem_traction_term(t1::PyObject, t2::PyObject, bdedge::Array{Int64, 2}, mmesh::Mesh)

Computes the boundary integral

\[\int_\Gamma \mathbf{t}(x,y)\cdot \delta \mathbf{u} d\mathbf{x}\]

Here $\mathbf{t}(x,y)$ is the boundary traction term

\[\mathbf{t}(x,y) = \begin{bmatrix} t_1(x,y) \\ t_2(x,y) \end{bmatrix}\]

and

\[\mathbf{u} = \begin{bmatrix} u_1(x,y) \\ u_2(x,y) \end{bmatrix}\]

\[t_1\]

and $t_2$ are defined on boundary Gauss points. See eval_f_on_boundary_edge.

compute_fem_traction_term(t::Array{Float64, 2},
bdedge::Array{Int64,2}, mesh::Mesh)

compute_fem_traction_term(t1::Array{Float64, 1}, t2::Array{Float64, 1},
bdedge::Array{Int64,2}, mesh::Mesh)

compute_von_mises_stress_term(K::Array{Float64}, u::Array{Float64}, m::Int64, n::Int64, h::Float64)

Compute the von Mises stress on the Gauss quadrature nodes.

compute_von_mises_stress_term(Se::Array{Float64,2},  m::Int64, n::Int64, h::Float64)

Se is a $4mn\times3$ array that stores the stress data at each Gauss point.

compute_von_mises_stress_term(K::Array{Float64, 3}, u::Array{Float64, 1}, mesh::Mesh)
compute_von_mises_stress_term(K::Array{Float64, 2}, u::Array{Float64, 1}, mesh::Mesh)

Computes the von Mises stress term

\[\sigma_v = \sqrt{\sigma_1^2 - \sigma_1 \sigma_2 + \sigma_2^2 + 3\sigma_{12}^2}\]

compute_fem_source_term1(f::Array{Float64},
m::Int64, n::Int64, h::Float64)

Computes the term

\[\int_\Omega f \delta u dx\]

Returns a $(m+1)\times (n+1)$ vector. f is a length $4mn$ vector, given by its values on Gauss points.

compute_fem_source_term1(f::PyObject,
m::Int64, n::Int64, h::Float64)

A differentiable kernel.

compute_fem_source_term1(f::PyObject, mesh::Mesh)

compute_fem_source_term1(f::Array{Float64,1}, mesh::Mesh)

compute_fem_source_term1(f::PyObject, mesh::Mesh3)

compute_fem_source_term1(f::Array{Float64,1}, mesh::Mesh3)

compute_fem_flux_term1(t::Array{Float64},
bdedge::Array{Int64,2}, m::Int64, n::Int64, h::Float64)

Computes the traction term

\[\int_{\Gamma} q \delta u \mathrm{d}\]

compute_strain_energy_term(S::Array{Float64, 2}, m::Int64, n::Int64, h::Float64)

Computes the strain energy

\[\int_{A} \sigma : \delta \varepsilon \mathrm{d}x\]

where $\sigma$ is provided by S, a $4mn \times 3$ matrix. The values $\sigma_{11}, \sigma_{22}, \sigma_{12}$ are defined on 4 Gauss points per element.

compute_strain_energy_term(S::PyObject,m::Int64, n::Int64, h::Float64)

A differentiable kernel.

compute_strain_energy_term(Sigma::Array{Float64, 2}, mmesh::Mesh)

Computes the strain energy term

\[\int_A \sigma : \varepsilon (\delta u) dx\]

Here $\sigma$ is a fourth-order tensor. Sigma is a ngauss × 3 matrix, each row represents $[\sigma_{11}, \sigma_{22}, \sigma_{12}]$ at each Gauss point.

The output is a length 2mmesh.ndof vector.

compute_strain_energy_term(Sigma::PyObject, mmesh::Mesh)

compute_strain_energy_term1(S::PyObject, m::Int64, n::Int64, h::Float64)

Computes the strain energy

\[\int_{A} \sigma : \delta \varepsilon \mathrm{d}x\]

where $\sigma$ is provided by S, a $4mn \times 2$ matrix. The values $\sigma_{31}, \sigma_{32}$ are defined on 4 Gauss points per element.

compute_strain_energy_term1(sigma::PyObject, m::Int64, n::Int64, h::Float64)

A differentiable operator.

compute_fem_viscoelasticity_strain_energy_term(ε0, σ0, ε, A, B, m, n, h)

Given the constitutive relation

\[\sigma^{n+1} = S \sigma^n + H (\varepsilon^{n+1}-\varepsilon^n),\]

this function computes

\[\int_A {\sigma:\delta \varepsilon}\mathrm{d} x = \underbrace{\int_A { B \varepsilon^{n+1}:\delta \varepsilon}\mathrm{d} x} + \underbrace{ \int_A { A \sigma^{n+1}:\delta \varepsilon}\mathrm{d} x - \int_A { B \varepsilon^{n+1}:\delta \varepsilon}\mathrm{d} x }_f\]

and returns $f$

compute_fvm_advection_term(v::Union{PyObject, Array{Float64, 2}},
u::Union{PyObject, Array{Float64,1}},m::Int64,n::Int64,h::Float64)

Computes the advection term using upwind schemes

\[\int_A \mathbf{v} \cdot \nabla u dx \]

Here $\mathbf{v}$ is a $mn\times 2$ matrix and $u$ is a length $mn$ vector. Zero boundary conditions are assumed. $u$ is a vector of length $m\times n$.

compute_interaction_term(p::Array{Float64, 1}, m::Int64, n::Int64, h::Float64)

Computes the FVM-FEM interaction term

\[ \begin{bmatrix} \int p \frac{\partial \delta u}{\partial x} dx \\ \int p \frac{\partial \delta v}{\partial y} dy \end{bmatrix} \]

The input is a vector of length $mn$. The output is a $2(m+1)(n+1)$ vector.

compute_interaction_term(p::PyObject, m::Int64, n::Int64, h::Float64)

A differentiable kernel.

compute_interaction_term(p::Union{PyObject,Array{Float64, 1}}, mesh::Mesh)

compute_fem_laplace_term1(u::PyObject,κ::PyObject,m::Int64,n::Int64,h::Float64)
compute_fem_laplace_term1(u::PyObject,m::Int64,n::Int64,h::Float64)
compute_fem_laplace_term1(u::Array{Float64},κ::PyObject, m::Int64,n::Int64,h::Float64)
compute_fem_laplace_term1(u::PyObject,κ::Array{Float64}, m::Int64,n::Int64,h::Float64)

Computes the Laplace term for a scalar function $u$

\[\int_\Omega K\nabla u \cdot \nabla (\delta u) \mathrm{d}x\]

Here κ is a vector of length $4mn$, and u is a vector of length $(m+1)(n+1)$.

When κ is not provided, the following term is calculated:

\[\int_\Omega \nabla u \cdot \nabla (\delta u) \mathrm{d}x\]

compute_fem_laplace_term1(u::Array{Float64, 1},nu::Array{Float64, 1}, mesh::Mesh)

compute_fem_laplace_term1(u::Union{PyObject, Array{Float64, 1}},
                            nu::Union{PyObject, Array{Float64, 1}},
                            mesh::Mesh)

compute_fem_laplace_term1(u::Union{PyObject, Array{Float64, 1}},
    nu::Union{PyObject, Array{Float64, 1}},
    mesh::Mesh3)

compute_fem_traction_term1(t::Array{Float64, 2},
bdedge::Array{Int64,2}, m::Int64, n::Int64, h::Float64)

Computes the traction term

\[\int_{\Gamma} t(n) \delta u \mathrm{d}\]

The number of rows of t is equal to the number of edges in bdedge. The output is a length $(m+1)*(n+1)$ vector.

Also see compute_fem_traction_term.

compute_fem_traction_term1(t::Array{Float64, 1},
bdedge::Array{Int64,2}, mesh::Mesh)

Computes the boundary integral

\[\int_{\Gamma} t(x, y) \delta u dx\]

Returns a vector of size dof.

compute_fem_traction_term1(t::PyObject,bdedge::Array{Int64, 2}, mmesh::Mesh)

Computes the boundary integral

\[\int_{\Gamma} t(x, y) \delta u dx\]

Returns a vector of size dof.

Here $t$ is defined on boundary Gauss points, e.g., the quantity computed from eval_f_on_boundary_edge.

Example

Here is an example to evaluate the boundary integral on a domain $[0,1]^2$

\[\int_{0}^1 t(x, 1) \delta u dx\]

mmesh = Mesh(10,10,0.1)
a = constant(1.0)
f = (x,y)->x+y*a
top = bcedge((x1,y1,x2,y2)->(y1>0.99) && (y2>0.99), mmesh)
t = eval_f_on_boundary_edge(f, top, mmesh; tensor_input = true)
T = compute_fem_traction_term1(t, top, mmesh)

compute_fem_boundary_mass_term1(u::Union{Array{Float64}, PyObject}, 
    c::Union{Array{Float64}, PyObject}, bdedge::Array{Int64, 2}, mmesh::Mesh)

Computes the term

\[\int_\Gamma cu \delta u ds\]

The parameters are

• u : a vector of length $N_v$
• bdedge : a $N_e \times 2$ integer array, the boundary edge to integrate on
• c: given by a vector of length $4N_e$; currently, each edge has 4 quadrature points;

The output is a $N_v\times N_v$ sparse matrix.

eval_f_on_gauss_pts(f::Function, m::Int64, n::Int64, h::Float64; tensor_input::Bool = false)

Evaluates f at Gaussian points and return the result as $4mn$ vector out (4 Gauss points per element)

If tensor_input = true, the function f is assumed to map a tensor to a tensor output.

eval_f_on_gauss_pts(f::Function, mesh::Mesh; tensor_input::Bool = false)

eval_f_on_dof_pts(f::Function, mesh::Mesh)

Evaluates f on the DOF points.

• For P1 element, the DOF points are FEM points and therefore eval_f_on_dof_pts is equivalent to eval_on_on_fem_pts.
• For P2 element, the DOF points are FEM points plus the middle point for each edge.

Returns a vector of length mesh.ndof.

eval_f_on_boundary_node(f::Function, bdnode::Array{Int64}, m::Int64, n::Int64, h::Float64)

Returns a vector of the same length as bdnode whose entries corresponding to bdnode nodes are filled with values computed from f.

f has the following signature

f(x::Float64, y::Float64)::Float64

eval_f_on_boundary_node(f::Function, bdnode::Array{Int64}, mesh::Mesh)

eval_strain_on_gauss_pts(u::Array{Float64}, m::Int64, n::Int64, h::Float64)

Computes the strain on Gauss points. Returns a $4mn\times3$ matrix, where each row denotes $(\varepsilon_{11}, \varepsilon_{22}, 2\varepsilon_{12})$ at the corresponding Gauss point.

eval_strain_on_gauss_pts(u::PyObject, m::Int64, n::Int64, h::Float64)

A differentiable kernel.

eval_strain_on_gauss_pts(u::Array{Float64}, mmesh::Mesh)

Evaluates the strain on Gauss points. u is a vector of size 2mmesh.ndof.

The output is a ngauss × 3 vector.

eval_strain_on_gauss_pts(u::PyObject, mmesh::Mesh)

eval_strain_on_gauss_pts1(u::PyObject, m::Int64, n::Int64, h::Float64)

A differentiable kernel.

eval_f_on_fvm_pts(f::Function, m::Int64, n::Int64, h::Float64; tensor_input::Bool = false)

Returns $f(x_i, y_i)$ where $(x_i,y_i)$ are FVM nodes.

If tensor_input = true, the function f is assumed to map a tensor to a tensor output.

eval_f_on_fvm_pts(f::Function, mesh::Mesh; tensor_input::Bool = false)

eval_f_on_fem_pts(f::Function, m::Int64, n::Int64, h::Float64; tensor_input::Bool = false)

Returns $f(x_i, y_i)$ where $(x_i,y_i)$ are FEM nodes.

If tensor_input = true, the function f is assumed to map a tensor to a tensor output.

eval_f_on_fem_pts(f::Function, mesh::Mesh; tensor_input::Bool = false)

eval_grad_on_gauss_pts1(u::Array{Float64,1}, m::Int64, n::Int64, h::Float64)

Evaluates $\nabla u$ on each Gauss point. Here $u$ is a scalar function.

The input u is a vector of length $(m+1)*(n+1)$. The output is a matrix of size $4mn\times 2$.

eval_grad_on_gauss_pts1(u::PyObject, m::Int64, n::Int64, h::Float64)

A differentiable kernel.

eval_grad_on_gauss_pts1(u::Union{Array{Float64,1}, PyObject}, mesh::Mesh)

eval_grad_on_gauss_pts(u::Array{Float64,1}, m::Int64, n::Int64, h::Float64)

Evaluates $\nabla u$ on each Gauss point. Here $\mathbf{u} = (u, v)$.

\[\texttt{g[i,:,:]} = \begin{bmatrix} u_x & u_y\\ v_x & v_y \end{bmatrix}\]

The input u is a vector of length $2(m+1)*(n+1)$. The output is a matrix of size $4mn\times 2 \times 2$.

eval_grad_on_gauss_pts(u::PyObject, m::Int64, n::Int64, h::Float64)

A differentiable kernel.

eval_scalar_on_boundary_edge(u::Union{PyObject, Array{Float64, 1}},
        edge::Array{Int64, 1}, mmesh::Mesh)

Returns an array of values on the Gauss quadrature nodes for each edge.

• u: nodal values
• edge: a $N\times 2$ integer array; each row represents an edge $(x_1, x_2)$

The returned array consists of $(y_1, y_2, \ldots)$

eval_f_on_boundary_edge(f::Function, bdedge::Array{Int64,2}, m::Int64, n::Int64, h::Float64)

Returns a vector of the same length as bdedge whose entries corresponding to bdedge nodes are filled with values computed from f.

f has the following signature

f(x::Float64, y::Float64)::Float64

eval_f_on_boundary_edge(f::Function, bdedge::Array{Int64, 2}, mesh::Mesh; tensor_input::Bool = false)

Evaluates f on the boundary Gauss points. Here f has the signature

f(Float64, Float64)::Float64

or

f(PyObject, PyObject)::PyObject

eval_strain_on_boundary_edge(u::Union{PyObject, Array{Float64, 1}},
edge::Array{Int64, 2}, mmesh::Mesh)

Returns an array of strain tensors on the Gauss quadrature nodes for each edge.

The returned value has size $N_e N_g\times 3$. Here $N_e$ is the number of edges, and $N_g$ is the number of Gauss points on the edge. Each row of edge, $(l,r)$, has the following property: $l < r$

eval_normal_and_shear_stress_on_boundary_edge(sigma::Union{PyObject, Array{Float64, 2}},
    edge::Array{Int64, 2}, mmesh::Mesh)

Calculates the normal and shear stress on boundary edges.

\[\sigma_n = (\sigma n)\cdot n, \tau = (\sigma n) \cdot m\]

Here $m = \begin{bmatrix}n_2\\-n_1\end{bmatrix}$

fem_impose_Dirichlet_boundary_condition(A::SparseMatrixCSC{Float64,Int64}, 
bd::Array{Int64}, m::Int64, n::Int64, h::Float64)

Imposes the Dirichlet boundary conditions on the matrix A.

Returns 2 matrix,

\[\begin{bmatrix} A_{BB} & A_{BI} \\ A_{IB} & A_{II} \end{bmatrix} \Rightarrow \begin{bmatrix} I & 0 \\ 0 & A_{II} \end{bmatrix}, \quad \begin{bmatrix} 0 \\ A_{IB} \end{bmatrix}\]

fem_impose_Dirichlet_boundary_condition(L::SparseTensor, bdnode::Array{Int64}, m::Int64, n::Int64, h::Float64)

A differentiable kernel for imposing the Dirichlet boundary of a vector-valued function.

fem_impose_Dirichlet_boundary_condition1(A::SparseMatrixCSC{Float64,Int64}, 
    bd::Array{Int64}, m::Int64, n::Int64, h::Float64)

Imposes the Dirichlet boundary conditions on the matrix A Returns 2 matrix,

\[\begin{bmatrix} A_{BB} & A_{BI} \\ A_{IB} & A_{II} \end{bmatrix} \Rightarrow \begin{bmatrix} I & 0 \\ 0 & A_{II} \end{bmatrix}, \quad \begin{bmatrix} 0 \\ A_{IB} \end{bmatrix}\]

bd must NOT have duplicates.

fem_impose_Dirichlet_boundary_condition1(L::SparseTensor, bdnode::Array{Int64}, m::Int64, n::Int64, h::Float64)

A differentiable kernel for imposing the Dirichlet boundary of a scalar-valued function.

fem_impose_Dirichlet_boundary_condition1(L::SparseTensor, bdnode::Array{Int64}, mesh::Mesh)

A differentiable kernel for imposing the Dirichlet boundary of a scalar-valued function.

impose_Dirichlet_boundary_conditions(A::Union{SparseArrays, Array{Float64, 2}}, rhs::Array{Float64,1}, bdnode::Array{Int64, 1}, 
    bdval::Array{Float64,1})
impose_Dirichlet_boundary_conditions(A::SparseTensor, rhs::Union{Array{Float64,1}, PyObject}, bdnode::Array{Int64, 1}, 
    bdval::Union{Array{Float64,1}, PyObject})

Algebraically impose the Dirichlet boundary conditions. We want the solutions at indices bdnode to be bdval. Given the matrix and the right hand side

\[\begin{bmatrix} A_{II} & A_{IB} \\ A_{BI} & A_{BB} \end{bmatrix}, \begin{bmatrix}r_I \\ r_B \end{bmatrix}\]

The function returns

\[\begin{bmatrix} A_{II} & 0 \\ 0 & I \end{bmatrix}, \begin{bmatrix}r_I - A_{IB} u_B \\ r_B \end{bmatrix}\]

impose_Dirichlet_boundary_conditions(A::SparseTensor, bdnode::Array{Int64, 1})

A helper function to impose homogeneous Dirichlet boundary condition.

impose_bdm_traction_boundary_condition1(gN::Array{Float64, 1}, bdedge::Array{Int64, 2}, mesh::Mesh)

Imposes the BDM traction boundary condition

\[\int_{\Gamma} \sigma \mathbf{n} g_N ds\]

Here $\sigma$ is a second-order tensor. gN is defined on the Gauss points, e.g.

gN = eval_f_on_boundary_edge(func, bdedge, mesh)

impose_bdm_traction_boundary_condition(g1:Array{Float64, 1}, g2:Array{Float64, 1},
bdedge::Array{Int64, 2}, mesh::Mesh)

Imposes the BDM traction boundary condition

\[\int_{\Gamma} \sigma \mathbf{n} \cdot \mathbf{g}_N ds\]

Here $\sigma$ is a fourth-order tensor. $\mathbf{g}_N = \begin{bmatrix}g_{N1}\\ g_{N2}\end{bmatrix}$ See impose_bdm_traction_boundary_condition1.

Returns a dof vector and a val vector.

visualize_mesh(mesh::Mesh)

Visualizes the unstructured meshes.

visualize_mesh(mmesh::Mesh3; filename::Union{Missing, String} = missing)

Visualizes the unstructured mesh mmesh. When filename is provided, a screenshot is saved to filename

visualize_pressure(U::Array{Float64, 2}, m::Int64, n::Int64, h::Float64)

Visualizes pressure. U is the solution vector.

visualize_displacement(u::Array{Float64, 2}, m::Int64, n::Int64, h::Float64)

Generates scattered plot animation for displacement $u\in \mathbb{R}^{(NT+1)\times 2(m+1)(n+1)}$.

visualize_displacement(u::Array{Float64, 1}, mmesh::Mesh)

visualize_displacement(u::Array{Float64, 2}, mmesh::Mesh)

Generates scattered plot animation for displacement $u\in \mathbb{R}^{(NT+1)\times 2N_{\text{dof}}}$.

visualize_stress(K::Array{Float64, 2}, U::Array{Float64, 2}, m::Int64, n::Int64, h::Float64; name::String="")

Visualizes displacement. U is the solution vector, K is the elasticity matrix ($3\times 3$).

visualize_stress(Se::Array{Float64, 2}, m::Int64, n::Int64, h::Float64; name::String="")

Visualizes the Von Mises stress. Se is the Von Mises at the cell center.

visualize_von_mises_stress(Se::Array{Float64, 2}, m::Int64, n::Int64, h::Float64; name::String="")

Visualizes the Von Mises stress.

visualize_von_mises_stress(K::Array{Float64}, u::Array{Float64, 1}, mmesh::Mesh, args...; kwargs...)

visualize_scalar_on_gauss_points(u::Array{Float64,1}, m::Int64, n::Int64, h::Float64, args...;kwargs...)

Visualizes the scalar u using pcolormesh. Here u is a length $4mn$ vector and the values are defined on the Gauss points

visualize_scalar_on_gauss_points(u::Array{Float64,1}, mesh::Mesh, args...;kwargs...)

Visualizes scalar values on Gauss points. For unstructured meshes, the values on each element are averaged to produce a uniform value for each element.

visualize_scalar_on_fvm_points(u::Array{Float64,1}, mmesh::Mesh3; filename::Union{Missing, String} = missing)

Visualizes the unstructured mesh mmesh with elements colored with values u. When filename is provided, a screenshot is saved to filename

visualize_scalar_on_fem_points(u::Array{Float64,1}, m::Int64, n::Int64, h::Float64, args...;kwargs...)

Visualizes the scalar u using pcolormesh. Here u is a length $(m+1)(n+1)$ vector and the values are defined on the FEM points

visualize_scalar_on_fem_points(u::Array{Float64,2}, m::Int64, n::Int64, h::Float64, args...;kwargs...)

Visualizes the scalar u using pcolormesh. Here u is a matrix of size $NT \times (m+1)(n+1)$ ($NT$ is the number of time steps) and the values are defined on the FEM points.

visualize_scalar_on_fem_points(u::Array{Float64,1}, mesh::Mesh, args...;
    with_mesh::Bool = false, kwargs...)

Visualizes the nodal values u on the unstructured mesh mesh.

• with_mesh: if true, the unstructured mesh is also plotted.

visualize_scalar_on_fem_points(U::Array{Float64,2}, mesh::Mesh)

visualize_scalar_on_fem_points(u::Array{Float64,1}, mmesh::Mesh3; filename::Union{Missing, String} = missing)

Visualizes the unstructured mesh mmesh with nodes colored with values u. When filename is provided, a screenshot is saved to filename

visualize_scalar_on_fvm_points(φ::Array{Float64, 3}, m::Int64, n::Int64, h::Float64;
vmin::Union{Real, Missing} = missing, vmax::Union{Real, Missing} = missing)

Generates scattered potential animation for the potential $\phi\in \mathbb{R}^{(NT+1)\times n \times m}$.

visualize_scalar_on_fvm_points(u::Array{Float64,1}, mesh::Mesh, args...;kwargs...)

visualize_scalar_on_fvm_points(u::Array{Float64,1}, mmesh::Mesh3; filename::Union{Missing, String} = missing)

Visualizes the unstructured mesh mmesh with elements colored with values u. When filename is provided, a screenshot is saved to filename

visualize_vector_on_fem_points(u1::Array{Float64,1}, u2::Array{Float64,1}, mesh::Mesh, args...;kwargs...)

Visualizes a vector on the mesh mesh. Here u1 and u2 are the $x$ and $y$ component of the vector. They are defined on the FEM nodes (not DOFs).

layer_model(u::Array{Float64, 1}, m::Int64, n::Int64, h::Float64)

Convert the vertical profile of a quantity to a layer model. The input u is a length $n$ vector, the output is a length $4mn$ vector, representing the $4mn$ Gauss points.

layer_model(u::PyObject, m::Int64, n::Int64, h::Float64)

A differential kernel for layer_model.

compute_vel(a::Union{PyObject, Array{Float64, 1}},
v0::Union{PyObject, Float64},psi::Union{PyObject, Array{Float64, 1}},
sigma::Union{PyObject, Array{Float64, 1}},
tau::Union{PyObject, Array{Float64, 1}},eta::Union{PyObject, Float64})

Computes $x = u_3(x_1, x_2)$ from rate and state friction. The governing equation is

\[a \sinh^{-1}\left( \frac{x - u}{\Delta t} \frac{1}{2V_0} e^{\frac{\Psi}{a}} \right) \sigma - \tau + \eta \frac{x-u}{\Delta t} = 0\]

compute_plane_strain_matrix(E::Float64, ν::Float64)

Computes the stiffness matrix for 2D plane strain. The matrix is given by

\[\frac{E(1-\nu)}{(1+\nu)(1-2\nu)}\begin{bmatrix} 1 & \frac{\nu}{1-\nu} & \frac{\nu}{1-\nu}\\ \frac{\nu}{1-\nu} & 1 & \frac{\nu}{1-\nu} \\ \frac{\nu}{1-\nu} & \frac{\nu}{1-\nu} & 1 \end{bmatrix}\]

compute_plane_strain_matrix(E::Union{PyObject, Array{Float64, 1}}, nu::Union{PyObject, Array{Float64, 1}})

Returns the pointwise plane strain matrix of size $N\times 3 \times 3$. Here $N$ is the length of $E$ or $\nu$. $N$ can be any number.

compute_plane_stress_matrix(E::Float64, ν::Float64)

Computes the stiffness matrix for 2D plane stress. The matrix is given by

\[\frac{E}{(1+\nu)(1-2\nu)}\begin{bmatrix} 1-\nu & \nu & 0\\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-2\nu}{2} \end{bmatrix}\]

compute_plane_stress_matrix(E::Union{PyObject, Array{Float64, 1}}, nu::Union{PyObject, Array{Float64, 1}})

Returns the pointwise plane stress matrix of size $N\times 3 \times 3$. Here $N$ is the length of $E$ or $\nu$. $N$ can be any number.

compute_space_varying_tangent_elasticity_matrix(mu::Union{PyObject, Array{Float64,1}},m::Int64,n::Int64,h::Float64,type::Int64=1)

Computes the space varying tangent elasticity matrix given $\mu$. It returns a matrix of size $4mn\times 2\times 2$

• If type==1, the $i$-th matrix will be

\[\begin{bmatrix}\mu_i & 0 \\ 0 & \mu_i \end{bmatrix}\]

• If type==2, the $i$-th matrix will be

\[\begin{bmatrix}\mu_i & 0 \\ 0 & \mu_{i+4mn} \end{bmatrix}\]

• If type==3, the $i$-th matrix will be

\[\begin{bmatrix}\mu_i & \mu_{i+8mn} \\ \mu_{i+8mn} & \mu_{i+4mn}\end{bmatrix}\]

mantle_viscosity(u::Union{Array{Float64}, PyObject},
    T::Union{Array{Float64}, PyObject}, m::Int64, n::Int64, h::Float64;
    σ_yield::Union{Float64, PyObject} = 300e6, 
    ω::Union{Float64, PyObject}, 
    η_min::Union{Float64, PyObject} = 1e18, 
    η_max::Union{Float64, PyObject} = 1e23, 
    E::Union{Float64, PyObject} = 9.0, 
    C::Union{Float64, PyObject} = 1000., N::Union{Float64, PyObject} = 2.)

\[\eta = \eta_{\min} + \min\left( \frac{\sigma_{\text{yield}}}{2\sqrt{\epsilon_{II}}}, \omega\min(\eta_{\max}, \eta) \right)\]

with

\[\epsilon_{II} = \frac{1}{2} \epsilon(u)\qquad \eta = C e^{E(0.5-T)} (\epsilon_{II})^{(1-n)/2n}\]

Here $\epsilon_{II}$ is the second invariant of the strain rate tensor, $C > 0$ is a viscosity pre-factor, $E > 0$ is the non-dimensional activation energy, $n > 0$ is the nonlinear exponent, $η_\min$, $η_\max$ act as minimum and maximum bounds for the effective viscosity, and $σ_{\text{yield}} > 0$ is the yield stress. $w\in (0, 1]$ is the weakening factor, which is used to incorporate phenomenological aspects that cannot be represented in a purely viscous flow model, such as processes which govern mega-thrust faults along the subduction interface, or partial melting near a mid-ocean ridge.

The viscosity of the mantle is governed by the high-temperature creep of silicates, for which laboratory experiments show that the creep strength is temperature-, pressure-, compositional- and stress-dependent.

The output is a length $4mn$ vector.

antiplane_viscosity(ε::Union{PyObject, Array{Float64}}, σ::Union{PyObject, Array{Float64}}, 
μ::Union{PyObject, Float64}, η::Union{PyObject, Float64}, Δt::Float64)

Calculates the stress at time $t_{n+1}$ given the strain at $t_{n+1}$ and stress at $t_{n}$. The governing equation is

\[\dot\sigma + \frac{\mu}{\eta}\sigma = 2\mu \dot\epsilon\]

The discretization form is

\[\sigma^{n+1} = \frac{1}{\frac{1}{\Delta t}+\frac{\mu}{\eta}}(2\mu\dot\epsilon^{n+1} + \frac{\sigma^n}{\Delta t})\]

update_stress_viscosity(ε2::Array{Float64,2}, ε1::Array{Float64,2}, σ1::Array{Float64,2}, 
invη::Array{Float64,1}, μ::Array{Float64,1}, λ::Array{Float64,1}, Δt::Float64)

Updates the stress for the Maxwell model

\[\begin{bmatrix} \dot\sigma_{xx}\ \dot\sigma_{yy}\ \dot\sigma_{xy} \end{bmatrix} + \frac{\mu}{\eta}\begin{bmatrix} 2/3 & -1/3 & 0 \ -1/3 & 2/3 & 0 \ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} \sigma_{xx}\ \sigma_{yy}\ \sigma_{xy}\end{bmatrix} = \begin{bmatrix} 2\mu + \lambda & \lambda & 0 \ \lambda & 2\mu + \lambda & 0 \ 0 & 0 & \mu \end{bmatrix}\begin{bmatrix} \dot\epsilon_{xx}\ \dot\epsilon_{yy}\ \dot\gamma_{xy} \end{bmatrix}\]

See here for details.

update_stress_viscosity(ε2::Union{PyObject,Array{Float64,2}}, ε1::Union{PyObject,Array{Float64,2}}, σ1::Union{PyObject,Array{Float64,2}}, 
invη::Union{PyObject,Array{Float64,1}}, μ::Union{PyObject,Array{Float64,1}}, λ::Union{PyObject,Array{Float64,1}}, Δt::Union{PyObject,Float64})

compute_pml_term(u::Union{Array{Float64,1}, PyObject},βprime::Union{Array{Float64,1}, PyObject},
c::Union{Array{Float64,1}, PyObject},nv::Union{Array{Float64,2}, PyObject}, mmesh::Mesh)

• u: a vector of length mmesh.ndof
• βprime: a vector of length $n_{\text{gauss}}$
• c: a tensor of size $n_{\text{gauss}}$
• nv: a matrix of size $n_{\text{gauss}}\times 2$

This function returns four outputs

\[\begin{aligned}k_1&=(c^2 n\partial_n u, n\partial_n\delta u)\\k_2&=(\beta'n\cdot(c^2 n\partial_n u), \delta u)\\ k_3&=(c^2\nabla^\parallel u, n\partial_n \delta u) + (c^2 n\partial_ n u, \nabla^\parallel \delta u)\\ k_4&= (c^2 \nabla^\parallel u, \nabla^\parallel \delta u)\end{aligned}\]

compute_pml_term(u::Union{Array{Float64,1}, PyObject},βprime::Union{Array{Float64,1}, PyObject},
λ::Union{Array{Float64,1}, PyObject},μ::Union{Array{Float64,1}, PyObject}, nv::Union{Array{Float64,2}, PyObject}, mmesh::Mesh)

• u: a vector of length mmesh.ndof
• βprime: a vector of length $n_{\text{gauss}}$
• λ, μ: Lam\'e parameters
• nv: a matrix of size $n_{\text{gauss}}\times 2$

This is a 2D version of compute_pml_term.

compute_pml_term(u::Union{Array{Float64,1}, PyObject},
    βprime::Union{Array{Float64,1}, PyObject},
    λ::Union{Array{Float64,1}, PyObject},μ::Union{Array{Float64,1}, PyObject}, 
    nv::Union{Array{Float64,2}, PyObject}, mmesh::Mesh3)

solve_slip_law(v, ψ, dc, v0, a, b, f0, Δt::Float64)

Solves one step of the slip law equation

\[\dot \psi = - \frac{|V|}{d_C}\left( a \sinh^{-1}\left( \frac{|V|}{2V_0}e^{\frac{\psi}{a}} \right) - f_0 + (b-a)*\log \frac{|V|}{V_0} \right)\]

We discretize the equation with a central difference scheme

\[\frac{\psi^{n+1} - \psi^{n-1}}{2\Delta t} = - \frac{|V|}{d_c}\left( a \sinh^{-1} \left( \frac{|V|}{2V_0} e^{\frac{\psi^{n+1} + \psi^{n-1}}{2a}{}} \right) - f_0 + (b-a) \log \frac{|V|}{V_0} \right)\]

• dc, v0, a, b, and f0 are scalars

solve_slip_law( 
    A::Union{Array{Float64,1}, PyObject}, 
    B::Union{Array{Float64,1}, PyObject}, 
    C::Union{Array{Float64,1}, PyObject},
    X0::Union{Array{Float64,1}, PyObject})

A helper function for solve_slip_law the nonlinear equation

\[x - A\sinh^{-1}(Bx) - C = 0\]

A, B, and C are vectors of the same length. X0 is the initial guess.

get_edge_dof(edges::Array{Int64, 2}, mesh::Mesh)
get_edge_dof(edges::Array{Int64, 1}, mesh::Mesh)

Returns the DOFs for edges, which is a K × 2 array containing vertex indices. The DOFs are not offset by nnode, i.e., the smallest edge DOF could be 1.

When the input is a length 2 vector, it returns a single index for the corresponding edge DOF.

get_edge_normal(edge::Array{Int64,1}, m::Int64, n::Int64, h::Float64)

Returns the normal vector given edge edge.

get_edge_normal(mmesh::Mesh)

get_edge_normal(bdedge::Array{Int64, 2}, mmesh::Mesh)

Returns the outer normal for all edges in bdedge.

get_boundary_edge_orientation(bdedge::Array{Int64, 2}, mmesh::Mesh)

Returns the orientation of the edges in bdedge. For example, if for a boundary element [1,2,3], assume [1,2] is the boundary edge, then

get_boundary_edge_orientation([1 2;2 1], mmesh) = [1.0;-1.0]

The return values for non-boundary edges in bdedge is undefined.

get_ngauss(mesh::Mesh)

Return the areas of triangles as an array.

get_ngauss(mesh::Mesh)

Return the total number of Gauss points.

get_ngauss(mesh::Mesh3)

Return the total number of Gauss points.

bcnode(desc::String, m::Int64, n::Int64, h::Float64)

Returns the node indices for the description. Multiple descriptions can be concatented via |

                upper
        |------------------|
left    |                  | right
        |                  |
        |__________________|

                lower

Example

bcnode("left|upper", m, n, h)

bcnode(mesh::Mesh; by_dof::Bool = true)

Returns all boundary node indices.

If by_dof = true, bcnode returns the global indices for boundary DOFs.

• For P2 elements, the returned values are boundary node DOFs + boundary edge DOFs (offseted by mesh.nnode)
• For BDM1 elements, the returned values are boundary edge DOFs + boundary edge DOFs offseted by mesh.nedge

bcnode(f::Function, mesh::Mesh; by_dof::Bool = true)

Returns the boundary node DOFs that satisfies f(x,y) = true.

f(x1::Float64, y1::Float64, x2::Float64, y2::Float64)::Bool

bcnode(mmesh::Mesh3)

bcedge(desc::String, m::Int64, n::Int64, h::Float64)

Returns the edge indices for description. See bcnode

bcedge(mesh::Mesh)

Returns all boundary edges as a set of integer pairs (edge vertices).

bcedge(f::Function, mesh::Mesh)

Returns all edge indices that satisfies f(x1, y1, x2, y2) = true Here the edge endpoints are given by $(x_1, y_1)$ and $(x_2, y_2)$.

interior_node(desc::String, m::Int64, n::Int64, h::Float64)

In contrast to bcnode, interior_node returns the nodes that are not specified by desc, including those on the boundary.

femidx(d::Int64, m::Int64)

Returns the FEM index of the dof d. Basically, femidx is the inverse of

(i,j) → d = (j-1)*(m+1) + i

fvmidx(d::Int64, m::Int64)

Returns the FVM index of the dof d. Basically, femidx is the inverse of

(i,j) → d = (j-1)*m + i

subdomain(f::Function, m::Int64, n::Int64, h::Float64)

Returns the subdomain defined by f(x, y)==true.

gauss_nodes(m::Int64, n::Int64, h::Float64)

Returns the node matrix of Gauss points for all elements. The matrix has a size $4mn\times 2$

gauss_nodes(mesh::Mesh)

gauss_nodes(mesh::Mesh3)

gauss_weights(mmesh::Mesh)

Returns the weights for each Gauss points.

fem_nodes(m::Int64, n::Int64, h::Float64)

Returns the FEM node matrix of size $(m+1)(n+1)\times 2$

fem_nodes(mesh::Mesh)

fem_nodes(mesh::Mesh3)

fvm_nodes(m::Int64, n::Int64, h::Float64)

Returns the FVM node matrix of size $(m+1)(n+1)\times 2$

fvm_nodes(mesh::Mesh)

fvm_nodes(mesh::Mesh3)

pcl_impose_Dirichlet_boundary_conditions(indices::Array{Int64, 2}, bdnode::Array{Int64, 1}, outdof::Int64)

Computes the Jacobian matrix for impose_Dirichlet_boundary_conditions. Assume that impose_Dirichlet_boundary_conditions transforms a sparse matrix A to B, and v_A and v_B are the nonzero entries, this function computes

\[J_{ij} = \frac{\partial (v_B)_j}{\partial (v_A)_i}\]

• indices: the 1-based $n_A\times 2$ index array for A;
• bdnode: boundary DOF (the same as inputs to impose_Dirichlet_boundary_conditions)
• outdof: the number of nonzero entries in B, i.e., length of v_B

pcl_compute_fem_laplace_matrix1(mmesh::Mesh)

Computes the Jacobian matrix for pcl_compute_fem_laplace_matrix1. Assume we contruct a sparse matrix via

A = compute_fem_laplace_matrix1(κ, mmesh)
v = A.o.values 

Then the function returns

\[J = \frac{\partial v_j}{\partial κ_i}\]

trim_coupled(pd::PoreData, Q::SparseMatrixCSC{Float64,Int64}, L::SparseMatrixCSC{Float64,Int64}, 
M::SparseMatrixCSC{Float64,Int64}, 
bd::Array{Int64}, Δt::Float64, m::Int64, n::Int64, h::Float64)

Assembles matrices from mechanics and flow and assemble the coupled matrix

\[\begin{bmatrix} \hat M & -\hat L^T\\ \hat L & \hat Q \end{bmatrix}\]

Q is obtained from compute_fvm_tpfa_matrix, M is obtained from compute_fem_stiffness_matrix, and L is obtained from compute_interaction_matrix.

coupled_impose_pressure(A::SparseMatrixCSC{Float64,Int64}, pnode::Array{Int64}, 
m::Int64, n::Int64, h::Float64)

Returns a trimmed matrix.

cholesky_factorize(A::Union{Array{<:Real,2}, PyObject})

Returns the cholesky factor of A. See cholesky_outproduct for details.

cholesky_outproduct(L::Union{Array{<:Real,2}, PyObject})

Returns $A = LL'$ where L (length=6) is a vectorized form of $L$ $L = \begin{matrix} l_1 & 0 & 0\\ l_4 & l_2 & 0 \\ l_5 & l_6 & l_3 \end{matrix}$ and A (length=9) is also a vectorized form of $A$

fem_to_fvm(u::Union{PyObject, Array{Float64}}, m::Int64, n::Int64, h::Float64)

Interpolates the nodal values of u to cell values.

fem_to_gauss_points(u::PyOject, m::Int64, n::Int64, h::Float64)
fem_to_gauss_points(u::Array{Float64,1}, m::Int64, n::Int64, h::Float64)

Given a vector of length $(m+1)(n+1)$, u, returns the function values at each Gauss point.

Returns a vector of length $4mn$.

fem_to_gauss_points(u::PyObject, mesh::Mesh)

fem_to_gauss_points(u::Array{Float64,1}, mesh::Mesh)

dof_to_gauss_points(u::PyObject, mesh::Mesh)
dof_to_gauss_points(u::Array{Float64,1}, mesh::Mesh)

Similar to fem_to_gauss_points. The only difference is that the function uses all DOFs–-which means, for quadratic elements, the nodal values on the edges are also used.