Edit on GitHub

API

Data Structures

AdFem.PoreData — Type

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
source
AdFem.Mesh — Type

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
  • elem_type: type of the element (P1, P2 or BDM1)

Internally, the mesh mmesh is represented by a collection of NNFEM_Element object with some other attributes

int nelem; // total number of elements
int nnode; // total number of nodes
int ngauss; // total number of Gauss points
int ndof; // total number of dofs
int order; // quadrature integration order
int degree; // Degree of Polynomials, 1 - P1 element, 2 - P2 element 
int elem_ndof; // 3 for P1, 6 for P2
MatrixXd GaussPts; // coordinates of Gauss quadrature points
std::vector<NNFEM_Element*> elements; // array of elements

The NNFEM_Element has data

VectorXd h;   // basis functions, elem_ndof × ng  
VectorXd hx;  // x-directional basis functions, elem_ndof × ng  
VectorXd hy;  // y-directional basis functions, elem_ndof × ng  
MatrixXd hs;  // shape functions for linear element, 3 × ng
VectorXd w;   // weight vectors, ng  
double area;  // area of the triangle
MatrixXd coord; // coordinates array, 3 × 2
int nnode; // total number of nodes 
int ngauss; // total number of Gauss points
int dof[6]; // global indices for both nodes and edges, note that edge indices are offset by `nv`
int node[3]; // global indices of local vertices
int edge[3]; // global indices of local edges
int ndof; // DOF, 3 for P1 element, 6 for P2 element
source

Matrix Assembling Functions

AdFem.compute_fem_stiffness_matrix — Function
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}\]
source
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\]
source
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)
source
AdFem.compute_interaction_matrix — Function
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.

source
compute_interaction_matrix(mesh::Mesh)
source
AdFem.compute_fvm_tpfa_matrix — Function
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)\]

Warning

No flow boundary condition is assumed.

source
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.

source
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.

Info

K can also be missing, in which case K is treated as a all-one vector.

source
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.

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

A differentiable kernel for compute_fvm_tpfa_matrix.

source
AdFem.compute_fem_mass_matrix — Function
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)$.

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

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

source
AdFem.compute_fem_mass_matrix1 — Function
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.

source
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.

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

A differentiable kernel.

source
compute_fem_mass_matrix1(rho::Union{PyObject, Array{Float64, 1}}, mesh::Mesh)
source
AdFem.compute_fem_stiffness_matrix1 — Function
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

source
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\]
source
AdFem.compute_fvm_advection_matrix — Function
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}\]
source
AdFem.compute_fem_laplace_matrix1 — Function
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.

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

K is duplicated on each Gauss point.

source
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.

source
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.

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

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

source
compute_fem_laplace_matrix1(kappa::PyObject, mesh::Mesh)
source
compute_fem_laplace_matrix1(kappa::Array{Float64,1}, mesh::Mesh)
source
AdFem.compute_fem_laplace_matrix — Function
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.

source
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$

source
compute_fem_laplace_matrix(kappa::Union{PyObject, Array{Float64, 1}}, mesh::Mesh)
source
AdFem.compute_fem_advection_matrix1 — Function
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)$

source
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)
source
AdFem.compute_fem_bdm_mass_matrix — Function
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.

source
compute_fem_bdm_mass_matrix(mmesh::Mesh)

Same as compute_fem_bdm_mass_matrix

source
compute_fem_bdm_mass_matrix(alpha::Array{Float64,1},beta::Array{Float64,1}, mmesh::Mesh)
source
AdFem.compute_fem_bdm_mass_matrix1 — Function
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.

source
compute_fem_bdm_mass_matrix1(mmesh::Mesh)

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

source
AdFem.compute_fem_bdm_div_matrix — Function
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.

source
AdFem.compute_fem_bdm_div_matrix1 — Function
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.

source
AdFem.compute_fem_bdm_skew_matrix — Function
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.

source

Vector Assembling Functions

AdFem.compute_fem_source_term — Function
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.

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

A differentiable kernel.

source
compute_fem_source_term(f1::Union{PyObject,Array{Float64,2}}, f2::Union{PyObject,Array{Float64,2}}, mesh::Mesh)
source
AdFem.compute_fvm_source_term — Function
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.

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

A differentiable kernel.

source
compute_fvm_source_term(f::Array{Float64, 1}, mmesh::Mesh)
source
AdFem.compute_fvm_mechanics_term — Function
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\]
source
compute_fvm_mechanics_term(u::PyObject, m::Int64, n::Int64, h::Float64)
source
AdFem.compute_fem_normal_traction_term — Function
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.

source
AdFem.compute_fem_traction_term — Function
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.

source
compute_fem_traction_term(t::Array{Float64, 2},
bdedge::Array{Int64,2}, mesh::Mesh)
source
compute_fem_traction_term(t1::Array{Float64, 1}, t2::Array{Float64, 1},
bdedge::Array{Int64,2}, mesh::Mesh)
source
AdFem.compute_von_mises_stress_term — Function
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.

source
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.

source
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)
source
AdFem.compute_fem_source_term1 — Function
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.

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

A differentiable kernel.

source
compute_fem_source_term1(f::PyObject, mesh::Mesh)
source
compute_fem_source_term1(f::Array{Float64,1}, mesh::Mesh)
source
AdFem.compute_fem_flux_term1 — Function
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}\]
source
AdFem.compute_strain_energy_term — Function
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.

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

A differentiable kernel.

source
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.

source
compute_strain_energy_term(Sigma::PyObject, mmesh::Mesh)
source
AdFem.compute_strain_energy_term1 — Function
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.

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

A differentiable operator.

source
AdFem.compute_fem_viscoelasticity_strain_energy_term — Function
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$

source
AdFem.compute_fvm_advection_term — Function
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$.

source
AdFem.compute_interaction_term — Function
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.

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

A differentiable kernel.

source
compute_interaction_term(p::Union{PyObject,Array{Float64, 1}}, mesh::Mesh)
source
AdFem.compute_fem_laplace_term1 — Function
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\]
source
compute_fem_laplace_term1(u::Array{Float64, 1},nu::Array{Float64, 1}, mesh::Mesh)
source
compute_fem_laplace_term1(u::Union{PyObject, Array{Float64, 1}},
                            nu::Union{PyObject, Array{Float64, 1}},
                            mesh::Mesh)
source
AdFem.compute_fem_traction_term1 — Function
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.

source
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.

source

Evaluation Functions

AdFem.eval_f_on_gauss_pts — Function
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.

source
eval_f_on_gauss_pts(f::Function, mesh::Mesh; tensor_input::Bool = false)
source
AdFem.eval_f_on_dof_pts — Function
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.

source
AdFem.eval_f_on_boundary_node — Function
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
source
eval_f_on_boundary_node(f::Function, bdnode::Array{Int64}, mesh::Mesh)
source
AdFem.eval_f_on_boundary_edge — Function
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
source
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

source
AdFem.eval_strain_on_gauss_pts — Function
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.

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

A differentiable kernel.

source
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.

source
eval_strain_on_gauss_pts(u::PyObject, mmesh::Mesh)
source
AdFem.eval_f_on_fvm_pts — Function
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.

source
eval_f_on_fvm_pts(f::Function, mesh::Mesh; tensor_input::Bool = false)
source
AdFem.eval_f_on_fem_pts — Function
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.

source
eval_f_on_fem_pts(f::Function, mesh::Mesh; tensor_input::Bool = false)
source
AdFem.eval_grad_on_gauss_pts1 — Function
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$.

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

A differentiable kernel.

source
eval_grad_on_gauss_pts1(u::Union{Array{Float64,1}, PyObject}, mesh::Mesh)
source
AdFem.eval_grad_on_gauss_pts — Function
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$.

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

A differentiable kernel.

source

Boundary Conditions

AdFem.fem_impose_Dirichlet_boundary_condition — Function
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}\]
source
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.

source
AdFem.fem_impose_Dirichlet_boundary_condition1 — 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.

source
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.

source
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.

source
AdFem.impose_Dirichlet_boundary_conditions — 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}\]
source
AdFem.impose_bdm_traction_boundary_condition1 — Function
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)
source
AdFem.impose_bdm_traction_boundary_condition — Function
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.

source

Visualization

In visualize_scalar_on_XXX_points, the first argument is the data matrix. When the data matrix is 1D, one snapshot is plotted. When the data matrix is 2D, it is understood as multiple snapshots at different time steps (each row is a snapshot). When the data matrix is 3D, it is understood as time step × height × width.

AdFem.visualize_pressure — Function
visualize_pressure(U::Array{Float64, 2}, m::Int64, n::Int64, h::Float64)

Visualizes pressure. U is the solution vector.

source
AdFem.visualize_displacement — Function
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)}$.

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

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

source
AdFem.visualize_stress — Function
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$).

source
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.

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

Visualizes the Von Mises stress.

source
visualize_von_mises_stress(K::Array{Float64}, u::Array{Float64, 1}, mmesh::Mesh, args...; kwargs...)
source
AdFem.visualize_scalar_on_gauss_points — Function
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

source
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.

source
AdFem.visualize_scalar_on_fem_points — Function
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

source
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.

source
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.
source
AdFem.visualize_scalar_on_fvm_points — Function
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}$.

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

Modeling Tools

AdFem.layer_model — Function
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.

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

A differential kernel for layer_model.

source
AdFem.compute_vel — Function
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\]
source
AdFem.compute_plane_strain_matrix — Function
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}\]
source
compute_plane_strain_matrix(E::Union{PyObject, Array{Float64, 1}}, nu::Union{PyObject, Array{Float64, 1}})
source
AdFem.compute_plane_stress_matrix — Function
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}\]
source
compute_plane_stress_matrix(E::Union{PyObject, Array{Float64, 1}}, nu::Union{PyObject, Array{Float64, 1}})
source
AdFem.compute_space_varying_tangent_elasticity_matrix — Function
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}\]
source
AdFem.mantle_viscosity — Function
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.

Info

See Towards adjoint-based inversion of time-dependent mantle convection with nonlinear viscosity for details.

source
AdFem.antiplane_viscosity — Function
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})\]
source
AdFem.update_stress_viscosity — Function
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.

source
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})
source

Mesh

AdFem.get_edge_dof — Function
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.

source
AdFem.get_boundary_edge_orientation — Function
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.

source
AdFem.bcnode — Function
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)
source
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
source
bcnode(f::Function, mesh::Mesh; by_dof::Bool = true)

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

Note

For BDM1 element and by_dof = true, because the degrees of freedoms are associated with edges, f has the signature

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

bcnode only returns DOFs on edges such that f(x1, y1, x2, y2)=true.

source
AdFem.bcedge — Function
bcedge(desc::String, m::Int64, n::Int64, h::Float64)

Returns the edge indices for description. See bcnode

source
bcedge(mesh::Mesh)

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

source
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)$.

source
AdFem.interior_node — Function
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 thosee on the boundary.

source
AdFem.femidx — Function
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
source
AdFem.fvmidx — Function
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
source
AdFem.subdomain — Function
subdomain(f::Function, m::Int64, n::Int64, h::Float64)

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

source
AdFem.gauss_nodes — Function
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$

source
gauss_nodes(mesh::Mesh)
source
AdFem.fem_nodes — Function
fem_nodes(m::Int64, n::Int64, h::Float64)

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

source
fem_nodes(mesh::Mesh)
source
AdFem.fvm_nodes — Function
fvm_nodes(m::Int64, n::Int64, h::Float64)

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

source
fvm_nodes(mesh::Mesh)
source

Misc

AdFem.trim_coupled — Function
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.

source
AdFem.coupled_impose_pressure — Function
coupled_impose_pressure(A::SparseMatrixCSC{Float64,Int64}, pnode::Array{Int64}, 
m::Int64, n::Int64, h::Float64)

Returns a trimmed matrix.

source
AdFem.cholesky_outproduct — Function
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$

source
AdFem.fem_to_fvm — Function
fem_to_fvm(u::Union{PyObject, Array{Float64}}, m::Int64, n::Int64, h::Float64)

Interpolates the nodal values of u to cell values.

source
AdFem.fem_to_gauss_points — Function
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$.

source
fem_to_gauss_points(u::PyObject, mesh::Mesh)
source
fem_to_gauss_points(u::Array{Float64,1}, mesh::Mesh)
source
AdFem.dof_to_gauss_points — Function
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.

source