Assemblers

Assemblers are used to assemble the system matrices and residual vectors. To calculate the contribution to these matrices, one or more of the following functions should be overloaded for the specific material:

Available assemblers

The following assemblers can be used to assemble the system matrix and vector:

Ferrite assemblers

Ferrite assemblers are supported and are used to assemble both a matrix and a vector. They are created by calling Ferrite's start_assemble function. They will primarily call element_routine! or facet_routine!. However, if not defined for the given material, automatic differentiation will be used to get the matrix contribution by calling element_residual! or facet_residual! instead.

ReAssembler

The ReAssembler assembles only the residual vector, and will thus call element_residual! or facet_residual!.

FerriteAssembly.ReAssemblerType
ReAssembler(r; fillzero=true, scaling=NoScaling())

ReAssembler will call element_residual! and assemble re into the system vector r. fillzero=true implies that r is zeroed upon construction of ReAssembler. Furthermore, remaining keyword arguments can enable the following features:

reset_scaling! is called when creating ReAssembler

source

KeReAssembler

The KeReAssembler assembles both the residual vector and the system matrix, similar to Ferrite's assemblers. It will also primarily call element_routine! or facet_routine! and fall back to automatic differentiation of element_residual! and facet_residual! if the former methods are not defined.

FerriteAssembly.KeReAssemblerType
KeReAssembler(K, r; fillzero=true, kwargs...)
KeReAssembler(a::Ferrite.AbstractAssembler; kwargs...)

The default KeReAssembler works just like a Ferrite.AbstractAssembler: It will call element_routine! and assemble Ke and re into the global system matrix K and vector r. However, it comes with the additional possible features, that are controllable via the keyword arguments:

  • ch::Union{Nothing,ConstraintHandler}=nothing: If a ConstraintHandler is given, local applications of constraints will be applied, using Ferrite.apply_assemble.
  • apply_zero: Required if a constraint handler is given, and forwarded to Ferrite.apply_assemble.
  • scaling: Calculate a scaling measure locally at the residual level, see e.g., ElementResidualScaling

a=Ferrite.start_assemble(K, r; fillzero=fillzero) is passed to the second definition if a matrix, K, and vector, r, are given as input.

reset_scaling! is called when creating KeReAssembler

source

Functions to overload

FerriteAssembly.element_routine!Function
element_routine!(
    Ke::AbstractMatrix, re::AbstractVector, state,
    ae::AbstractVector, material, cellvalues, buffer)

The main function to be overloaded for the specific material. In most cases, the same implementation can be used for different cellvalues (e.g. for different interpolation orders) This function should modify the element stiffness matrix Ke, the residual re, and potentially state. The element degree of freedom values, ae, are filled by NaNs unless a is passed to work!.

The following variables can be obtained from buffer.

source
FerriteAssembly.element_routine!(
    Ke, re, state::Vector{<:MMB.AbstractMaterialState}, ae, 
    m::MMB.AbstractMaterial, cv::AbstractCellValues, buffer)

Solve the weak form

\[ \int_\Omega [\boldsymbol{\delta u}\otimes\nabla]^\mathrm{sym} : \boldsymbol{\sigma}\ \mathrm{d}\Omega = \int_\Gamma \boldsymbol{\delta u} \cdot \boldsymbol{t}\ \mathrm{d}\Gamma + \int_\Omega \boldsymbol{\delta u} \cdot \boldsymbol{b}\ \mathrm{d}\Omega\]

where $\sigma$ is calculated with the material_response function from MaterialModelsBase.jl. Note that create_cell_state is already implemented for <:AbstractMaterial.

source
FerriteAssembly.element_residual!Function
element_residual!(
    re::AbstractVector, state,
    ae::AbstractVector, material, cellvalues, buffer)

To calculate the element tangent stiffness Ke automatically by using ForwardDiff, it is possible to overload element_residual! instead of element_routine!. See element_routine! for a description of the input parameters.

Note

MethodError with ForwardDiff.Dual
When using automatic differentiation for elements with state variables (or other mutating values in e.g. cache), an error will be thrown if trying to change the type in many cases. When mutating state, call ForwardDiff.value() on the value to be assigned after it will no longer be used to calculate re. (If done on a value which is later affects re inside the element, the tangent, Ke, will be wrong.)

source
FerriteAssembly.facet_routine!Function
facet_routine!(Ke, re, ae, material, facetvalues, facetbuffer)

Calculate contributions to the stiffness matrix and residual vector from a facet domain. It can be used, for example, to implement Robin boundary conditions.

source
FerriteAssembly.facet_residual!Function
facet_residual!(re, ae, material, facetvalues, facetbuffer)

Calculate contributions to the residual vector from a facet domain. Internally, this is used for calculating Neumann boundary conditions.

source

Residual scaling

There are many options for how to scale the residual in finite element simulations. This package does not intend to implement many different options, but does give the user the option to calculate scaling contributions from each cell, which may be useful. By defining a scaling that is passed to KeReAssembler or ReAssembler, it can be updated based on the output from each cell.

One type of scaling, ElementResidualScaling, is included as described below. Its code can be used as a template for how to include custom scaling that works on the element level.

ElementResidualScaling

FerriteAssembly.ElementResidualScalingType
ElementResidualScaling(dh::AbstractDofHandler, p=Val(2), T=Float64)

Create tolerance scaling based on the sum of the p-norm of each cell's residual vector, separately for each field. I.e. pseudo-code for field :u the scaling factor::T is

for each cell
    element_routine!(Ke, re, args...) # Calculate re
    factor += sum(abs.(re[dof_range(dh, :u)]).^p)^(1/p)

Note that p=Val(2) is a special case that is supported, for different values it should be given as a Real. p=2 is equivalent to Val(2), but is less efficient.

source