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.start_assemble
(AssemblerSparsityPattern
orAssemblerSymmetricSparsityPattern
)KeReAssembler
ReAssembler
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.ReAssembler
— TypeReAssembler(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:
scaling
: Calculate a scaling measure locally at the residual level, see e.g.,ElementResidualScaling
.
reset_scaling!
is called when creating ReAssembler
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.KeReAssembler
— TypeKeReAssembler(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 aConstraintHandler
is given, local applications of constraints will be applied, usingFerrite.apply_assemble
.apply_zero
: Required if a constraint handler is given, and forwarded toFerrite.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
Functions to overload
FerriteAssembly.element_routine!
— Functionelement_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 NaN
s unless a
is passed to work!
.
The following variables can be obtained from buffer
.
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
.
FerriteAssembly.element_residual!
— Functionelement_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.
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.)
FerriteAssembly.facet_routine!
— Functionfacet_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.
FerriteAssembly.facet_residual!
— Functionfacet_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.
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.ElementResidualScaling
— TypeElementResidualScaling(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.