API Reference

Domain

Date structure for the computatational domain.

• nnodes: Int64, number of nodes (each quadratical quad element has 9 nodes)

• nodes: Float64[nnodes, ndims], coordinate array of all nodes

• neles: number of elements

• elements: a list of neles element arrays, each element is a struct

• ndims: Int64, dimension of the problem space

• state: a matrix of size nnodes×ndims. Current displacement of all nodal freedoms, state[1:nnodes] are for the first direction.

• Dstate: nnodes×ndims. Previous displacement of all nodal freedoms, Dstate[1:nnodes] are for the first direction.

• LM: Int64[neles][ndims], LM(e,d) is the global equation number (active freedom number) of element e's d th freedom,

     ∘ -1 means fixed (time-independent) Dirichlet
  
     ∘ -2 means time-dependent Dirichlet
  
     ∘ >0 means the global equation number

• DOF: a matrix of size neles×ndims, DOF(e,d) is the global freedom (including both active and inactive DOFs) number of element e's d th freedom.

• ID: a matrix of size nnodes×ndims. ID(n,d) is the equation number (active freedom number) of node n's $d$-th freedom,

     ∘ -1 means fixed (time-independent) Dirichlet
  
     ∘ -2 means time-dependent Dirichlet
  
     ∘ >0 means the global equation number

• neqs: Int64, number of equations, a.k.a., active freedoms

• eq_to_dof: an integer vector of length neqs, map from to equation number (active freedom number) to the freedom number (Int64[1:nnodes] are for the first direction)

• dof_to_eq: a bolean array of size nnodes×ndims, map from freedom number(Int64[1:nnodes] are for the first direction) to booleans (active freedoms(equation number) are true)

• EBC: Int64[nnodes, ndims], EBC[n,d] is the displacement boundary condition of node n's dth freedom, -1 means fixed(time-independent) Dirichlet boundary nodes -2 means time-dependent Dirichlet boundary nodes

• g: Float64[nnodes, ndims], values for fixed(time-independent) Dirichlet boundary conditions of node n's dth freedom,

• FBC: Int64[nnodes, ndims], FBC[n,d]` is the force load boundary condition of node n's dth freedom, -1 means constant(time-independent) force load boundary nodes -2 means time-dependent force load boundary nodes

• fext: Float64[neqs], constant (time-independent) nodal forces on these freedoms

• Edge_Traction_Data: n × 3 integer matrix for natural boundary conditions.

1. Edge_Traction_Data[i,1] is the element id,

2. Edge_Traction_Data[i,2] is the local edge id in the element, where the force is exterted (should be on the boundary, but not required)

3. Edge_Traction_Data[i,3] is the force id, which should be consistent with the last component of the Edge_func in the Globdat

• time: Float64, current time
• npoints: Int64, number of points (each quadratical quad element has 4 points, npoints==nnodes, when porder==1)
• node_to_point: Int64[nnodes]:map from node number to point point, -1 means the node is not a geometry point

Domain(nodes::Array{Float64}, elements::Array, ndims::Int64 = 2,
EBC::Union{Missing, Array{Int64}} = missing, g::Union{Missing, Array{Float64}} = missing, FBC::Union{Missing, Array{Int64}} = missing, 
f::Union{Missing, Array{Float64}} = missing, edge_traction_data::Array{Int64,2}=zeros(Int64,0,3))

Creating a finite element domain.

• nodes: coordinate array of all nodes, a nnodes × 2 matrix

• elements: element array. Each element is a material struct, e.g., PlaneStrain.

• ndims: dimension of the problem space. For 2D problems, ndims = 2.

• EBC: nnodes × ndims integer matrix for essential boundary conditions EBC[n,d]is the displacement boundary condition of noden`'s $d$-th freedom,

  ∘ -1: fixed (time-independent) Dirichlet boundary nodes

  ∘ -2: time-dependent Dirichlet boundary nodes

• g: nnodes × ndims double matrix, values for fixed (time-independent) Dirichlet boundary conditions of node n's $d$-th freedom,

• FBC: nnodes × ndims integer matrix for nodal force boundary conditions.

FBC[n,d] is the force load boundary condition of node n's dth freedom,

∘ -1 means constant(time-independent) force load boundary nodes

∘ -2 means time-dependent force load boundary nodes

• f: nnodes × ndims double matrix, values for constant (time-independent) force load boundary conditions of node n's $d$-th freedom,

• Edge_Traction_Data: n × 3 integer matrix for natural boundary conditions.

Edge_Traction_Data[i,1] is the element id, Edge_Traction_Data[i,2] is the local edge id in the element, where the force is exterted (should be on the boundary, but not required) Edge_Traction_Data[i,3] is the force id, which should be consistent with the last component of the Edge_func in the Globdat

For time-dependent boundary conditions (EBC or FBC entries are -2), the corresponding f or g entries are not used.

GlobalData

Store data for finite element updates. Assume the problem has n freedoms,

• state: a vector of length $n$. Displacement array at the current time, only for active freedoms. The ordering is based on the equation number.
• Dstate: a vector of length $n$. Displacement array at the previous time.
• velo: a vector of length $n$. Velocity array at the current time.
• acce: a vector of length $n$. Acceleration array at the current time.
• time: float, current time.
• M: a matrix of size $n\times n$ spares mass matrix
• Mlumped: a vector of length $n$ lumped mass array
• MID: Float64[n, nd1] off-diagonal part of the mass matrix, between the active freedoms and the time-dependent Dirichlet freedoms, assume there are nd time-dependent Dirichlet freedoms
• EBC_func: displacement $d$, velocity $v$, and acceleration $a$ from time-dependent Dirichlet boundary conditions

The length of each output is the same as number of "-2" in EBC array. The ordering is direction major, i.e., $u_1, u_3, \ldots, v_1, v_3, \ldots$

• FBC_func: time-dependent load boundary condition.

Here $f$ is a vector. Its length is the same as number of "-2" in FBC array. The ordering is direction major, i.e., $u_1, u_3, \ldots, v_1, v_3, \ldots$

• Body_func: time-dependent/independent body force function.

Here $f$ is a vector or a matrix (second dimension is 2) depending on the dimension of state variables. The output is a $N\times n_{\text{dim}}$ matrix, where $N$ is the length of $x_{\text{array}}$ or $y_{\text{array}}$, and $n_{\text{dim}}$ is the dimension of the problem (1 or 2).

• Edge_func: time-dependent/independent traction load.

Here $f$ is a vector. Its length is the same as the length of $x_{\text{array}}$ or $y_{\text{array}}$.

GlobalData(state::Union{Array{Float64,1},Missing},
        Dstate::Union{Array{Float64,1},Missing},
        velo::Union{Array{Float64,1},Missing},
        acce::Union{Array{Float64,1},Missing}, 
        neqs::Int64,
        EBC_func::Union{Function, Nothing}=nothing, FBC_func::Union{Function, Nothing}=nothing,
        Body_func::Union{Function,Nothing}=nothing, Edge_func::Union{Function,Nothing}=nothing)

The size of state, Dstate, velo, acce must be neqs, i.e., the active DOFs. If they are missing, they are treated as zeros.

state, Dstate, velo, acce can be missing, in which case they are interpreted as zeros.

getEqns(domain::Domain, iele::Int64)

Gets the equation numbers (active freedom numbers) of the element. This excludes both the time-dependent and time-independent Dirichlet boundary conditions.

getDofs(domain::Domain, iele::Int64)

Get the global freedom numbers of the element

• domain: Domain
• iele: Int64, element number

Return: Int64[], the global freedom numbers of the element (ordering in local element ordering)

getCoords(domain::Domain, el_nodes::Array{Int64})

Get the coordinates of several nodes (possibly in one element)

• domain: Domain
• el_nodes: Int64[n], node array

Return: Float64[n, ndims], the coordinates of these nodes

getNGauss(domain::Domain)

Gets the total number of Gauss quadrature points.

getGaussPoints(elem::Continuum)

Returns the Gauss quadrature nodes of the element in the undeformed domain

getGaussPoints(domain::Domain)

Returns all Gauss points as a $n_g\times 2$ matrix, where $n_g$ is the total number of Gauss points.

getState(domain::Domain, el_dofs::Array{Int64})

Get the displacements of several nodes (possibly in one element)

• domain: Domain
• el_nodes: Int64[n], node array

Return: Float64[n, ndims], the displacements of these nodes

getStrain(domain::Domain)

Computes the strain from the domain data. The output is $n_g\times 3$ matrix, where $n_g$ is the total number of Gauss points. Each row is the strain tensor $\begin{bmatrix} \epsilon_{xx} & \epsilon_{yy} & \gamma_{xy} \end{bmatrix}$

getDStrain(domain::Domain)

Computes the strain at last time step from the domain data. The output is $n_g\times 3$ matrix, where $n_g$ is the total number of Gauss points. Each row is the strain tensor $\begin{bmatrix} \epsilon_{xx} & \epsilon_{yy} & \gamma_{xy} \end{bmatrix}$

getStress(domain::Domain, Δt::Float64 = 0.0; save_trace::Bool = false)

Returns the stress based on domain.state and domain.Dstate. If save_trace is true, the stress is also saved to domain.stress, which is useful for visualization.

The output is $n_g\times 3$ matrix, where $n_g$ is the total number of Gauss points. Each row is the strain tensor $\begin{bmatrix} \sigma_{xx} & \sigma_{yy} & \sigma_{xy} \end{bmatrix}$

getElems(domain::Domain)

Returns the element connectivity matrix $n_e \times 4$. This function implicitly assumes that all elements are quadrilateral.

getStressHistory(domain::Domain)

Returns the stress history.

getStrainHistory(domain::Domain)

Returns the strain history.

getStateHistory(domain::Domain)

Returns the state history.

getEdgeForce(elem::Continuum, iedge::Float64, fvalue::Array{Float64,2})

Returns the force imposed by boundary tractions.

fvalue is a $n_{edge_gauss}\times 2$ matrix, which is ordered the same as the Gaussian points in undeformed parent edge element.

    The element nodes are ordered as 
    #   4 ---- 3             #   4 --7-- 3
    #                        #   8   9   6 
    #   1 ---- 2             #   1 --5-- 2
    for porder=1     or          porder=2
    iedge 1, 2, 3, 4 are (1,2), (2,3), (3,4), (4,1)
                    are (1,2,5), (2,3,6), (3,4,7), (4,1,8)

Returns the nodal force due to the traction on the iedge-th edge of the element $\int_{s} \mathbf{f}(\mathbf{x})\cdot \delta \mathbf{u}(\mathbf{x}) d s = \int_{e} \mathbf{f}(\xi)\cdot \delta \mathbf{u}(\xi) |\frac{\partial \mathbf{x}}{\partial \xi}| d \xi$

getEdgeForce(domain::Domain, globdat::GlobalData, time::Float64)

Computes the edge force vector $F_\mathrm{edge}$ defined in domain.edge_traction_data

• globdat: GlobalData
• domain: Domain, finite element domain, for data structure
• time: Float64, current time step size

getBodyForce(elem::Continuum, fvalue::Array{Float64,2})

Returns the body force.

• fvalue is a $n_{gauss}\times 2$ matrix, which is ordered the same as Gaussian points in the undeformed parent element.

Returns the nodal force due to the body force $\int_{e} \mathbf{f}(\mathbf{x})\cdot \delta \mathbf{u}(\mathbf{x}) d \mathbf{x} = \int_{e} \mathbf{f}(\mathbf{\xi})\cdot \delta \mathbf{u}(\mathbf{\xi}) |\frac{\partial \mathbf{x}}{\partial \mathbf{\xi}}| d \mathbf{\xi}$

getBodyForce(elem::Continuum, fvalue::Array{Float64, 1})

Returns

on a specific element $A$ fvalue has the same length as number of Gauss points.

getBodyForce(domain::Domain, globdat::GlobalData, time::Float64)

Computes the body force vector $F_\mathrm{body}$ of length neqs

• globdat: GlobalData
• domain: Domain, finite element domain, for data structure
• Δt: Float64, current time step size

getMassMatrix(elem::Continuum)

Returns the mass matrix and lumped mass matrix of the element elem.

getNodes(elem::Continuum)

Alias for elem.elnodes

commitHistory(elem::Continuum)

Updates the historic parameters in the material properties.

commitHistory(domain::Domain)

Update current step strain and stress in the history map of the domain. This is essential for visualization and time dependent constitutive relations.

FiniteStrainContinuum(coords::Array{Float64}, elnodes::Array{Int64}, props::Dict{String, Any}, ngp::Int64=2)

FiniteStrainContinuum

Constructs a finite strain element.

• eledim: spatial dimension of the element (default = 2).
• mat: constitutive law, a length #elem vector of materials such as PlaneStress
• elnodes: the node indices in this finite element, an integer array
• props: property dictionary
• coords: coordinates of the vertices of the element
• dhdx, weights, hs: data for integral
• stress: stress at each quadrature points

Example

#   Local degrees of freedom 
#   4 ---- 3
#
#   1 ---- 2

nx = 10
ny = 5
h = 0.1
element = FiniteStrainContinuum[]
prop = Dict("name"=> "PlaneStrain", "rho"=> 0.0876584, "E"=>0.07180760098, "nu"=>0.4)
for j = 1:ny
    for i = 1:nx 
        n = (nx+1)*(j-1) + (i-1)+1
        elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        ngp = 3 # 3 x 3 Gauss points per element 
        coords = [(i-1)*h (j-1)*h
                    i*h (j-1)*h
                    i*h j*h
                    (i-1)*h j*h]
        push!(element, FiniteStrainContinuum(coords,elnodes, prop, ngp))
    end
end

SmallStrainContinuum(coords::Array{Float64}, elnodes::Array{Int64}, props::Dict{String, Any}, ngp::Int64=2)

SmallStrainContinuum

Constructs a small strain element.

• eledim: spatial dimension of the element (default = 2).
• mat: constitutive law, a length #elem vector of materials such as PlaneStress
• elnodes: the node indices in this finite element, an integer array
• props: property dictionary
• coords: coordinates of the vertices of the element
• dhdx: list of ngp shape functions for first order derivatives $\nabla \phi(x)$ (ndof×2) on the Gaussian points
• weights: weight vector of length n_gauss_points, for numerical quadrature
• hs: list of ngp shape functions for function values $\phi(x)$ (length ndof vectors) on the Gaussian points
• stress: stress at each quadrature points; this field is reserved for visualization.

Example

#   Local degrees of freedom 
#   4 ---- 3
#
#   1 ---- 2

nx = 10
ny = 5
h = 0.1
element = SmallStrainContinuum[]
prop = Dict("name"=> "PlaneStrain", "rho"=> 0.0876584, "E"=>0.07180760098, "nu"=>0.4)
for j = 1:ny
    for i = 1:nx 
        n = (nx+1)*(j-1) + (i-1)+1
        elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        ngp = 3 # 3 x 3 Gauss points per element 
        coords = [(i-1)*h (j-1)*h
                    i*h (j-1)*h
                    i*h j*h
                    (i-1)*h j*h]
        push!(element, SmallStrainContinuum(coords,elnodes, prop, ngp))
    end
end

getInternalForce(elem::SmallStrainContinuum, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)

Returns the internal force term. state and Dstate are restriction of full state variables to this element.

getStiffAndForce(elem::SmallStrainContinuum, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)

Returns the internal force term and the stiffness matrix. state and Dstate are restriction of full state variables to this element.

getStiffAndForce(elem::SmallStrainContinuum, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)

Returns the internal force term and the stiffness matrix. state and Dstate are restriction of full state variables to this element.

getStrain(elem::SmallStrainContinuum, state::Array{Float64})

Returns the strain of this element. state is restricted to this variable.

PlaneStress

Creates a plane stress element

• H: Linear elasticity matrix, $3\times3$
• E: Young's modulus
• ν: Poisson's ratio
• ρ: density
• σ0: stress at the last time step
• σ0_: (for internal use), stress to be updated in commitHistory
• ε0: strain at the last time step
• ε0_: (for internal use), strain to be updated in commitHistory

Example

prop = Dict("name"=> "PlaneStress", "rho"=> 0.0876584, "E"=>0.07180760098, "nu"=>0.4)
mat = PlaneStress(prop)

PlaneStress(prop::Dict{String, Any})

prop should contain at least the following three fields: E, nu, rho

PlaneStrain

Creates a plane strain element

• H: Linear elasticity matrix, $3\times3$
• E: Young's modulus
• ν: Poisson's ratio
• ρ: density
• σ0: stress at the last time step
• σ0_: (for internal use), stress to be updated in commitHistory
• ε0: strain at the last time step
• ε0_: (for internal use), strain to be updated in commitHistory

Example

prop = Dict("name"=> "PlaneStrain", "rho"=> 0.0876584, "E"=>0.07180760098, "nu"=>0.4)
mat = PlaneStrain(prop)

PlaneStrain(prop::Dict{String, Any})

prop should contain at least the following three fields: E, nu, rho

Pascon, João Paulo. "Large deformation analysis of plane-stress hyperelastic problems via triangular membrane finite elements." International Journal of Advanced Structural Engineering (2019): 1-20.

assembleInternalForce(domain::Domain, Δt::Float64 = 0.0)

Computes the internal force vector $F_\mathrm{int}$ of length neqs

• globdat: GlobalData
• domain: Domain, finite element domain, for data structure
• Δt: Float64, current time step size

Only the information in domain is used for computing internal force. Therefore, the boundary conditions in domain must be set appropriately.

assembleInternalForce(domain::Domain, nn::Function, E_all::PyObject, DE_all::PyObject, w∂E∂u_all::PyObject, σ0_all::PyObject)

Computes local internal force fint and then assemble to Fint, which generates inverse problem automatically.

• domain: finite element domain
• nn: constitutive relation for expressing stress = f(strain), assuming stress and strain are defined on Gauss points ((neles*nGauss) × nstrains).
• E_all: strain data of size (neles*nGauss) × nstrains at the current time step.
• DE_all: strain data of size (neles*nGauss) × nstrains at the last time step.
• w∂E∂u_all: sensitivity matrix of size (neles*nGauss) x ndofs_per_element x nstrains; neles*nGauss is the number of Gaussian quadrature points, ndofs_per_element is the number of freedoms per element, and nstrain is the number of strain components. The sensitivity matrix already considers the quadrature weights.

where $e$ is the element index, $g$ is the Gaussian index.

• σ0_all: stress data of size neles*nGauss×nstrains at the last time step.

Return:

• \[F_{\mathrm{int}}\]
  : internal force vector of length neqns
• \[\sigma_{\mathrm{all}}\]
  : predicted stress at current step, a matrix of size (neles*nGauss) × nstrains

assembleStiffAndForce(domain::Domain, Δt::Float64 = 0.0)

Computes the internal force and stiffness matrix.

• domain: Domain, finite element domain, for data structure
• Δt: Float64, current time step size

Returns a length neqs vector $F_{\mathrm{int}}$ and neqs×neqs sparse stiffness matrix.

assembleMassMatrix!(globaldat::GlobalData, domain::Domain)

Computes the constant sparse mass matrix $M_{\mathrm{mass}}$, the lumped mass matrix $M_{\mathrm{lump}}$ due to time-dependent Dirichlet boundary conditions, and store them in globaldat.

Here M is a neqns×neqns matrix, and $M_{ID}$ is a neqns×nd matrix. $M_{\mathrm{lump}}$ assumes that the local mass matrix is a diagonal matrix.

• globdat: GlobalData
• domain: Domain, finite element domain, for data structure

getExternalForce(self::Domain, globaldat::GlobalData, fext::Union{Missing,Array{Float64}}=missing)

Computes external force vector at globaldat.time, This includes all the body force, external load, and internal force caused by acceleration.

setConstantDirichletBoundary!(self::Domain, EBC::Array{Int64}, g::Array{Float64})

Bookkeepings for time-independent Dirichlet boundary conditions. Only called once in the constructor of domain. It updates the fixed (time-independent Dirichlet boundary) state entries and builds both LM and DOF arrays.

• self: Domain

• EBC: Int64[nnodes, ndims], EBC[n,d] is the displacement boundary condition of node n's dth freedom,

  ∘ -1 means fixed(time-independent) Dirichlet boundary nodes

  ∘ -2 means time-dependent Dirichlet boundary nodes

• g: Float64[nnodes, ndims], values for fixed (time-independent) Dirichlet boundary conditions of node n's dth freedom,

Bookkeepings for time-independent Nodal force boundary conditions. Only called once in the constructor of domain. It updates the fixed (time-independent Nodal forces) state entries and builds both LM and DOF arrays.

• self: Domain

• FBC: Int64[nnodes, ndims], FBC[n,d] is the displacement boundary condition of node n's dth freedom,

  ∘ -1 means fixed (time-independent) Nodal force freedoms

  ∘ -2 means time-dependent Nodal force freedoms

• f: Float64[nnodes, ndims], values for fixed (time-independent) Neumann boundary conditions of node n's dth freedom,

#The name is misleading

updateStates!(domain::Domain, globaldat::GlobalData)
update time-dependent Dirichlet boundary condition to globaldat.time

Update state and Dstate in domain. This includes

• Copy state variable values for active DOFs from globaldat

• Set time-dependent essential boundary conditions using globaldat.EBC

The time-independent boundary conditions are inherented from last time step.

updateTimeDependentEssentialBoundaryCondition!(domain::Domain, globaldat::GlobalData)

If there exists time-dependent Dirichlet boundary conditions, updateTimeDependentEssentialBoundaryCondition! must be called to update the boundaries in domain. This function is called by updateStates!

This function updates state data in domain.

ExplicitSolverStep(globdat::GlobalData, domain::Domain, Δt::Float64)

Central Difference explicit solver for M a + fint(u) = fext(u). a, v, u are acceleration, velocity and displacement.

GeneralizedAlphaSolverStep(globdat::GlobalData, domain::Domain, Δt::Float64, 
ρ::Float64 = 0.0, ε::Float64 = 1e-8, ε0::Float64 = 1e-8, maxiterstep::Int64=100, 
η::Float64 = 1.0, failsafe::Bool = false, verbose::Bool = false)

Implicit solver for $Ma + f_{int}(u) = fext$ Here $a$, $v$, $u$ are acceleration, velocity and displacement respectively.

• ρ: controls the damping effect of the α-scheme, ρ∈[0,1], ρ=1 corresponds to the maximum damping
• ε: Float64, absolute error for Newton convergence
• ε0: Float64, relative error for Newton convergence
• max_iter: Int64, maximum iteration number for Newton convergence
• η: Float64, Newton step size at the first iteration
• failsafe: Bool, if failsafe is true, when the Newton fails to converge, revert back, and return false

The nonlinear $\alpha$

'a_{n+1}' is solved by

As for \alpha_m and $\alpha_f$

use the current states a, v, u, time in globdat, and update these stetes to next time step update domain history, when failsafe is true, and Newton's solver fails, nothing will be changed.

You need to call SolverInitial! before the first time step, if f^{ext}0 != 0. SolverInitial! updates a0 in the globdat.acce a0 = M^{-1}(- f^{int}(u0) + f^{ext}_0)

We assume globdat.acce[:] = a_0 and so far initialized to 0 We also assume the external force is conservative (it does not depend on the current deformation)

ImplicitStaticSolver(globdat::GlobalData, domain::Domain; 
    N::Int64 = 10, ε::Float64 = 1.e-6, maxiterstep::Int64=100)

Solves $K(u) = F$ using the incremental load method. Specifically, at step $i$, we solve

• globdat, GlobalData

• domain, Domain

• N, an integer, load stepping steps

• ε, Float64, absolute error for Newton convergence

• maxiterstep, Int64, maximum iteration number for Newton convergence

ImplicitStaticSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject}, 
    nn::Function, θ::Union{Array{Float64, 1}, PyObject},
    Fext::Union{Array{Float64, 1}, PyObject, Missing}=missing)

Solves the static problem

using Newton's method. Users provide nn, which is a function that outputs stress and stress sensitivity given the strain tensor.

• d0: a vector of length 2domain.nnodes, the fixed Dirichlet DOF should be populated with boundary values.

• Fext: external force, a vector of length domain.neqs

SolverInitial!(Δt::Float64, globdat::GlobalData, domain::Domain)

You need to call SolverInitial! before the first time step, if $f^{ext}_0 \neq 0$

SolverInitial(Δt::Float64, globdat::GlobalData, domain::Domain)

Similar to SolverInitial!, but returns the (displacement, velocity, acceleartion) tuple.

meshread(gmshfile::String)

Reads a gmsh file gmshfile and return (nodes, elements) tuple.

psread(gmshfile::String)

Reads the physical group information from gmshfile. Returns a dictionary

PhysicalGroupName => Index

Index is a $n\times 2$ coordinate matrix. We do not output indices because there may be postprocessing procedures after creating the mesh.

Add a curve loop (a closed wire) formed by the curves curveTags.

Forms a physical group of dimension dim and tag tag with name name

visualize_boundary(domain::Domain, direction::String="x")

Visualizes the boundary conditions. The boundary configuration is shown in the direction direction.

visualize_displacement(domain::Domain)
visualize_displacement(u::Array{Float64, 2}, domain::Domain)

Animation of displacements using points.

visualize_displacement(domain::Domain)
visualize_displacement(nodes::Array{Float64,2}, elems::Array{Int64, 2})

Visualizes the mesh.

visualize_scalar_on_scoped_body(d::Array{Float64}, domain::Domain, kwargs...)
visualize_scalar_on_scoped_body(s::Array{Float64, 1}, d::Array{Float64,1}, domain::Domain;
    scale_factor::Float64 = -1.0, kwargs...)

Plot the scalar on scoped body. For example, s can be the von Mises stress tensor.

visualize_scalar_on_undeformed_body(s::Array{Float64, 1}, domain::Domain; kwargs...)
visualize_scalar_on_undeformed_body(s::Array{Float64, 2}, domain::Domain; frames::Int64 = 20, kwargs...)

Plots or animates scalar values s on the domain domain

visualize_total_deformation_on_scoped_body(d_all::Array{Float64,2}, domain::Domain;
scale_factor::Float64 = 1.0, frames::Int64 = 20)

Visualizes the total deformation

visualize_von_mises_stress(domain::Domain, t_step::Int64; kwargs...)

Plot of von Mises stress tensors at time step t_step.

visualize_von_mises_stress(domain::Domain; frames::Int64 = 20, kwargs...)

Animation of von Mises stress tensors.

visualize_von_mises_stress_on_scoped_body(d_all::Array{Float64,2}, domain::Domain;
scale_factor::Real = -1.0, frames::Int64 = 20)

Similar to visualize_von_mises_stress, but the domain can be a deformed body.

visualize_x_deformation_on_scoped_body(d_all::Array{Float64,2}, domain::Domain;
scale_factor::Real = 1.0, frames::Int64 = 20, kwargs...)

Visualizes the $x$ directional displacement.

visualize_y_deformation_on_scoped_body(d_all::Array{Float64,2}, domain::Domain;
scale_factor::Real = 1.0, frames::Int64 = 20, kwargs...)

Visualizes the $y$ directional displacement.

init_nnfem(domain::Domain)

Prepares domain for use in custom operators.

s_eval_strain_on_gauss_points(state::Union{Array{Float64,1}, PyObject})

Computes the strain on Gauss points in the small strain case. state is the full displacement vector.

s_compute_stiffness_matrix(k::Union{Array{Float64,3}, PyObject})

Computes the small strain stiffness matrix. $k$ is a $n\times 3\times 3$ matrix, where $n$ is the total number of Gauss points. Returns a SparseTensor.

s_compute_internal_force_term(stress::Union{Array{Float64,2}, PyObject})

Computes the internal force $\int_\Omega \sigma : \delta \epsilon dx$ Only active DOFs are considered.

f_eval_strain_on_gauss_points(state::Union{Array{Float64,1}, PyObject})

Computes the strain on Gauss points in the finite strain case. state is the full displacement vector.

f_compute_internal_force_term(stress::Union{Array{Float64,2}, PyObject}, 
    state::Union{Array{Float64,1}, PyObject},
    domain::Domain)

Computes the internal force for finite strain continuum

Only active DOFs are considered.

ExplicitSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject}, 
    v0::Union{Array{Float64, 1}, PyObject}, 
    a0::Union{Array{Float64, 1}, PyObject}, 
    Δt::Float64, NT::Int64, 
    H::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
    Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; strain::String = "small")

Differentiable Explicit Solver.

• d0, v0, a0: initial full displacement, velocity, and acceleration.

• Δt: time step

• Hs: linear elasticity matrix at each Gauss point

• Fext: external force, $\mathrm{NT}\times n$, where $n$ is the active dof. The external force includes all body forces, external load forces (also called edge forces in NNFEM) and boundary acceleration-induced forces.

• ubd, abd: boundary displacementt and acceleration, $\mathrm{NT}\times m$, where $m$ is time-dependent boundary DOF. Time-independent boundary conditions are extracted from domain.

• strain_type (default = "small"): small strain or finite strain

ExplicitSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject}, 
    v0::Union{Array{Float64, 1}, PyObject}, 
    a0::Union{Array{Float64, 1}, PyObject}, 
    Δt::Float64, NT::Int64, 
    nn::Function,
    Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; strain_type::String = "small"))

Similar to ExplicitSolver; however, the constituve relation from $\epsilon$ to $\sigma$ must be provided by the function nn.

ExplicitSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject}, 
    v0::Union{Array{Float64, 1}, PyObject}, 
    a0::Union{Array{Float64, 1}, PyObject}, 
    σ0::Union{Array{Float64, 1}, PyObject}, 
    ε0::Union{Array{Float64, 1}, PyObject}, 
    Δt::Float64, NT::Int64, 
    nn::Function,
    Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; strain_type::String = "small")

Similar to ExplicitSolver; however, the constitutive relation has the form

Here the strain and stress are $n \times 3$ tensors. $n$ is the total number of Gaussian points and can be obtained via getNGauss(domain).

ExplicitSolverTime(Δt::Float64, NT::Int64)

Returns the times for explicit solver. Boundary conditions and external forces should be given at these times.

GeneralizedAlphaSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject}, 
    v0::Union{Array{Float64, 1}, PyObject}, 
    a0::Union{Array{Float64, 1}, PyObject}, 
    Δt::Float64, NT::Int64, 
    Hs::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
    Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; ρ::Float64 = 0.0)

Differentiable Generalized $\alpha$ scheme. This is an extension of αscheme provided in ADCME. This function does not support damping and variable time step (for efficiency).

• d0, v0, a0: initial full displacement, velocity, and acceleration.

• Δt: time step

• Hs: linear elasticity matrix at each Gauss point

• Fext: external force, $\mathrm{NT}\times n$, where $n$ is the active dof. The external force includes all body forces, external load forces (also called edge forces in NNFEM) and boundary acceleration-induced forces.

• ubd, abd: boundary displacementt and acceleration, $\mathrm{NT}\times m$, where $m$ is boundary DOF. Time-independent boundary conditions are extracted from domain.

GeneralizedAlphaSolver does not support finite-strain continuum yet.

GeneralizedAlphaSolver(globdat::GlobalData, domain::Domain,
    d0::Union{Array{Float64, 1}, PyObject}, 
    v0::Union{Array{Float64, 1}, PyObject}, 
    a0::Union{Array{Float64, 1}, PyObject}, 
    Δt::Float64, NT::Int64, 
    Cs::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
    Hs::Union{Array{Float64, 3}, Array{Float64, 2}, PyObject},
    Fext::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    ubd::Union{Array{Float64, 2}, PyObject, Missing}=missing,
    abd::Union{Array{Float64, 2}, PyObject, Missing}=missing; ρ::Float64 = 0.0)

Solve linear dynamical structural problem using a generalized alpha solver. Cs and Hs are the corresponding linear viscosity and linear elasticity matrix.

GeneralizedAlphaSolverTime(Δt::Float64, NT::Int64;ρ::Float64 = 0.0)

Returns the times for the generalized $\alpha$ solver. Boundary conditions and external forces should be given at these times.

compute_boundary_info(domain::Domain, globdat::GlobalData, ts::Array{Float64})

Computes the boundary information ubd and abd

compute_external_force(domain::Domain, globdat::GlobalData, ts::Union{Array{Float64}, Missing} = missing)

Computes the external force (body force, edge force and force due to boundary acceleration).

compute_stress_rivlin_saunders(strain::Union{PyObject, Array{Float64,2}},c1::Union{PyObject, Float64},c2::Union{PyObject, Float64})

Computes the stress using the plane stress incompressible Rivlin Saunders model.

s_compute_stiffness_matrix1(k::Union{PyObject, Array{Float64,3}}, domain::Domain)

Computes the stiffness matrix

\int_\Omega (K\nabla u) \cdot \nabla \delta u dx