PDE Operators

mutable struct PDEOperator{T<:Real, APT<:AssemblyPatternType, AT<:AssemblyType} <: AbstractPDEOperator

common structures for all finite element operators that are assembled with GradientRobustMultiPhysics; better look at the AssemblyPatternType and the constructors

LinearForm(
    operator::Type{<:??},
    action::AbstractAction;
    name,
    AT,
    regions,
    factor,
    store
) -> GradientRobustMultiPhysics.PDEOperator{Float64, LinearForm, ON_CELLS}

generates a LinearForm L(v) = (f,operator(v)) from an action

• operator : operator applied to test function
• action : action that computes a result to be multiplied with test function operator

Optional arguments:

• regions : specifies in which regions the operator should assemble, default [0] means all regions
• name : name for this LinearForm that is used in print messages
• AT : specifies on which entities of the grid the LinearForm is assembled (default: ON_CELLS)
• factor : additional factor that is multiplied during assembly
• store : stores a vector of the discretised LinearForm with the latest assembly result

Details on the action: The action is an Action consisting of a kernel function with interface (result, input, ...) and additional argument information. During assembly input is ignored (only in this constructor for LinearForms). The result computed by the kernel function is multiplied (dot product) with the operator evaluation of the test function.

LinearForm(
    operator::Type{<:??},
    f::GradientRobustMultiPhysics.AbstractUserDataType;
    name,
    kwargs...
) -> GradientRobustMultiPhysics.PDEOperator{Float64, LinearForm, ON_CELLS}

generates a LinearForm L(v) = (f,operator(v)) from a DataFunction

• operator : operator applied to test function
• action : DataFunction, evaluation is multiplied with test function operator

Optional arguments:

• regions : specifies in which regions the operator should assemble, default [0] means all regions
• name : name for this LinearForm that is used in print messages
• AT : specifies on which entities of the grid the LinearForm is assembled (default: ON_CELLS)
• factor : additional factor that is multiplied during assembly
• store : stores a vector of the discretised LinearForm with the latest assembly result

LinearForm(operator_test::Type{<:??}; ...)
LinearForm(
    operator_test::Type{<:??},
    operators_current::Vector{DataType};
    ...
) -> GradientRobustMultiPhysics.PDEOperator{Float64, LinearForm, ON_CELLS}
LinearForm(
    operator_test::Type{<:??},
    operators_current::Vector{DataType},
    coeff_from::Vector{Int64};
    ...
) -> GradientRobustMultiPhysics.PDEOperator{Float64, LinearForm, ON_CELLS}
LinearForm(
    operator_test::Type{<:??},
    operators_current::Vector{DataType},
    coeff_from::Vector{Int64},
    action::AbstractAction;
    regions,
    name,
    AT,
    factor,
    store
) -> GradientRobustMultiPhysics.PDEOperator{Float64, LinearForm, ON_CELLS}

Creates a (PDE description level) LinearForm based on:

• operator_test : operator for the test function (assumes linearity for that part)
• operators_current : additional operators for other unknowns
• coefffrom : either PDE unknown ids or block ids for CurrentSolution given to assemblyoperator! that should be used for operators_current
• action : an Action with kernel of interface (result, input, kwargs) that takes input (= all but last operator evaluations) and computes result to be dot-producted with test function evaluation (if no action is specified, the full input vector is dot-producted with the test function operator evaluation)

Optional arguments:

• regions: specifies in which regions the operator should assemble, default [0] means all regions
• name : name for this LinearForm that is used in print messages
• AT : specifies on which entities of the grid the LinearForm is assembled (default: ON_CELLS)
• factor : additional factor that is multiplied during assembly
• store : stores a vector of the LinearForm with the latest assembly result (e.g. when the operators sits in a system block that has to be reassembled in an iterative scheme)

Details on the action: The action is an Action consisting of a kernel function with interface (result, input, ...) and additional argument information. During assembly input will be filled with the operator evaluations of the other unknowns (i.e. operatorscurrent). The result computed by the kernel function is multiplied (dot product) with the operator evaluation of the test function (i.e. operatortest). If no action is given, the assembly tries to multiply the operator evaluations (that would have been given as input) directly.

BilinearForm(
    operators_linear::Vector{DataType},
    operators_current::Vector{DataType},
    coeff_from::Vector{Int64},
    action::AbstractAction;
    name,
    AT,
    APT,
    apply_action_to,
    factor,
    regions,
    transposed_assembly,
    also_transposed_block,
    transpose_factor,
    store
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

generates a BilinearForm defined by the following arguments:

• operators_linear : operator for the two linear arguments (usually ansatz and test function)
• operators_current : additional operators for other unknowns
• coefffrom : either PDE unknown ids or block ids for CurrentSolution given to assemblyoperator! that should be used for operators_current
• action : tells how to further combine the operatorscurrent+operatoransatz evaluations (=input of action) to a result that is multiplied with the test function operator (if no action is specified, the full input vector is dot-producted with the test function operator evaluation)

Optional arguments:

• applyactionto : specifies which of the two linear arguments is part of the action input ([1] = ansatz, [2] = test)
• regions : specifies in which regions the operator should assemble, default [0] means all regions
• name : name for this BilinearForm that is used in print messages
• AT : specifies on which entities of the grid the BilinearForm is assembled (default: ON_CELLS)
• APT : specifies the subtype of the APT_BilinearForm AssemblyPattern used for assembly (e.g. for lumping (wip))
• factor : additional factor that is multiplied during assembly
• transposed_assembly : transposes the resulting assembled matrix (consider true here for non-symmetric operators)
• alsotransposedblock : when true toggles assembly of transposed system matrix block
• transpose_factor : factor for transposed block (default = factor)
• store : stores a matrix of the BilinearForm with the latest assembly result (e.g. when the operators sits in a system block that has to be reassembled in an iterative scheme)

Details on the action: The action is an Action consisting of a kernel function with interface (result, input, ...) and additional argument information. During assembly input will be filled with the operator evaluations of the other unknowns (i.e. operatorcurrent, if specified) and appended to that the operator evaluation of one of the two linear argument (decided by applyaction_to). The result computed by the kernel function is multiplied (dot product) with the operator evaluation of the other linear argument. If no action is given, the assembly tries to multiply the operator evaluations (that would have been given as input) directly.

BilinearForm(
    operators_linear::Vector{DataType};
    ...
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}
BilinearForm(
    operators_linear::Vector{DataType},
    action::AbstractAction;
    kwargs...
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

same as other constructor but with operators_current =

LaplaceOperator(
;
    ...
) -> GradientRobustMultiPhysics.PDEOperator{Float64, SymmetricBilinearForm, ON_CELLS}
LaplaceOperator(
    κ;
    name,
    AT,
    ∇,
    regions,
    store
) -> Union{Nothing, GradientRobustMultiPhysics.PDEOperator{Float64, SymmetricBilinearForm, ON_CELLS}}

constructor for a bilinearform that describes a(u,v) = κ (∇u,∇v) where kappa is some constant (diffusion) coefficient.

ReactionOperator(
;
    ...
) -> GradientRobustMultiPhysics.PDEOperator{Float64, SymmetricBilinearForm, ON_CELLS}
ReactionOperator(
    α;
    ...
) -> Union{Nothing, GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}, GradientRobustMultiPhysics.PDEOperator{Float64, SymmetricBilinearForm, ON_CELLS}}
ReactionOperator(
    α,
    ncomponents;
    name,
    AT,
    id,
    regions,
    store
) -> Union{Nothing, GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}, GradientRobustMultiPhysics.PDEOperator{Float64, SymmetricBilinearForm, ON_CELLS}}

constructor for a bilinearform a(u,v) = (αu,v) or (u,αv) with some coefficient α that can be a number or an AbstractDataFunction.

LagrangeMultiplier(
    operator::Type{<:??};
    name,
    AT,
    action,
    regions,
    store,
    factor
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

constructor for a bilinearform that describes a(u,v) = (A(operator(u)), id(v)) and assembles a second transposed block at the block of the transposed PDE coordinates. It is intended to use to render one unknown of the PDE the Lagrange multiplier for another unknown by putting this operator on the coressponding subdiagonal block of the PDE description.

Example: LagrangeMultiplier(Divergence) is used to render the pressure the LagrangeMultiplier for the velocity divergence constraint in the Stokes prototype.

ConvectionOperator(
    a_from::Int64,
    a_operator::Type{<:??},
    xdim::Int64,
    ncomponents::Int64;
    name,
    AT,
    a_to,
    factor,
    ansatz_operator,
    test_operator,
    regions,
    newton,
    store,
    transposed_assembly,
    bonus_quadorder
) -> Union{GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}, GradientRobustMultiPhysics.PDEOperator{Float64, NonlinearForm, ON_CELLS}}

constructs a convection term of the form c(a,u,v) = (aoperator(a)*ansatzoperator(u),test_operator(v)) as a BilinearForm (or NonlinearForm, see newton argument)

• a_from : id of registered unknown to be used in the spot a
• a_operator : operator applied to a
• xdim : expected space dimension
• ncomponents : expected numer of components of a

optional arguments:

• newton : generates a NonlinearForm instead of a BilinearForm that triggers assembly of Newton terms for c(u,u,v)
• ato : position of a argument, set ato = 2 to trigger assembly of c(u,a,v)
• ansatz_operator : operator used in the spot u (default: Gradient)
• test_operator : operator used in the spot v (default: Identity)
• factor : additional factor multiplied in assemblxy (default: 1)
• regions : specifies in which regions the operator should assemble, default [0] means all regions
• name : name for this operator that is used in print messages
• AT : specifies on which entities of the grid the operator is assembled (default: ON_CELLS)
• store : stores a matrix of the operator with the latest assembly result

ConvectionOperator(
    β::GradientRobustMultiPhysics.AbstractUserDataType,
    ncomponents::Int64;
    name,
    store,
    AT,
    ansatz_operator,
    test_operator,
    transposed_assembly,
    regions
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

constructor for a bilinearform that describes a(u,v) = ((β ⋅ ∇) u,v) with some user-specified DataFunction β. The user also has to specify the number of components (ncomponents) the convection is applied to. The operators for u and v can be changed (if this leads to something reasonable).

ConvectionRotationFormOperator(
    beta::Int64,
    beta_operator::Type{<:??},
    xdim::Int64,
    ncomponents::Int64;
    name,
    AT,
    factor,
    ansatz_operator,
    test_operator,
    regions
)

constructor for a bilinearform a(u,v) = (beta x curl(u),v) where beta is the id of some unknown vector field of the PDEDescription, u and v are also vector-fields and x is the cross product (so far this is only implemented in 2D)

HookStiffnessOperator1D(
    μ;
    name,
    regions,
    ∇,
    store
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

constructor for a bilinearform a(u,v) = (μ ∇u,∇v) where C is the 1D stiffness tensor for given μ.

HookStiffnessOperator2D(
    μ,
    λ;
    name,
    AT,
    regions,
    ϵ,
    store
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

constructor for a bilinearform a(u,v) = (C ϵ(u), ϵ(v)) where C is the 2D stiffness tensor for isotropic media in Voigt notation, i.e. ℂ ϵ(u) = 2 μ ϵ(u) + λ tr(ϵ(u)) for Lame parameters μ and λ

In Voigt notation ℂ is a 3 x 3 matrix
ℂ = [c11,c12,  0
     c12,c11,  0
       0,  0,c33]

where c33 = μ, c12 = λ and c11 = 2*c33 + c12

Note: ϵ is the symmetric part of the gradient (in Voigt notation)

HookStiffnessOperator3D(
    μ,
    λ;
    name,
    AT,
    regions,
    ϵ,
    store
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

constructor for a bilinearform a(u,v) = (C ϵ(u), ϵ(v)) where C is the 3D stiffness tensor for isotropic media in Voigt notation, i.e. ℂ ϵ(u) = 2 μ ϵ(u) + λ tr(ϵ(u)) for Lame parameters μ and λ

In Voigt notation ℂ is a 6 x 6 matrix
ℂ = [c11,c12,c12,  0,  0,  0
     c12,c11,c12,  0,  0,  0
     c12,c12,c11,  0,  0,  0
       0,  0,  0,c44,  0,  0
       0,  0,  0,  0,c44,  0
       0,  0,  0,  0,  0,c44]   

where c44 = μ, c12 = λ and c11 = 2*c44 + c12

Note: ϵ is the symmetric part of the gradient (in Voigt notation)

NonlinearForm(
    operator_test::Type{<:??},
    operators_current::Vector{DataType},
    coeff_from::Vector{Int64},
    action_kernel,
    argsizes::Vector{Int64};
    name,
    AT,
    newton,
    sparse_jacobian,
    jacobian,
    dependencies,
    bonus_quadorder,
    store,
    factor,
    regions
) -> GradientRobustMultiPhysics.PDEOperator{Float64, NonlinearForm, ON_CELLS}

generates a NonlinearForm defined by the following arguments:

• operator_test : operator for the test function
• operators_current : additional operators for other unknowns
• coefffrom : either PDE unknown ids or block ids for CurrentSolution given to assemblyoperator! that should be used for operators_current
• action_kernel : function of interface (result, input, ...) that computes the nonlinear quantity that should be multiplied with the testfunction operator
• argsizes : dimensions of [result, input] of kernel function

Optional arguments:

• dependencies : code String for additional dependencies of the kernel/jacobians (substring of "XTI")
• jacobian : default = "auto" triggers automatic computation of jacobians by ForwardDiff, otherwise user can specify a function of interface (jacobian, input, ...) with matching dimensions and dependencies
• sparse_jacobian : use sparsity detection and sparse matrixes for local jacobians ?
• regions : specifies in which regions the operator should assemble, default [0] means all regions
• name : name for this NonlinearForm that is used in print messages
• AT : specifies on which entities of the grid the NonlinearForm is assembled (default: ON_CELLS)
• factor : additional factor that is multiplied during assembly
• store : stores a matrix of the discretised NonlinearForm with the latest assembly result
• bonus_quadorder : increases the quadrature order in assembly accordingly (additional to usual quadorder based on used FESpaces)

Some details: Given some operator G(u), the Newton iteration reads DG u_next = DG u - G(u) which is added to the rest of the (linear) operators in the PDEDescription. The local jacobians (= jacobians of the operator kernel) to build DG needed for this are computed by automatic differentation (ForwardDiff). The user can also specify a jacobian kernel function by hand (which may improve assembly times).

For default dependencies both the kernel functions for the operator and its jacobian have to satisfy the interface

function name(result,input,...)

where input is a vector of the operators of the solution and result is either what then is multiplied with operator2 of the testfunction (or the jacobian).

FVConvectionDiffusionOperator(
    beta_from::Int64;
    μ
) -> FVConvectionDiffusionOperator{Float64}

finite-volume convection diffusion operator (for cell-wise P0 rho)

• div(μ ∇ρ + β ρ)

For μ = 0, the upwind divergence div_upw(β ρ) is generated For μ > 0, TODO

DiagonalOperator(; ...) -> DiagonalOperator{Float64}
DiagonalOperator(
    value::Real;
    name,
    onlynz,
    regions
) -> DiagonalOperator{<:Real}

puts value on the diagonal entries of the cell dofs within given regions

if onlyz == true only values that are zero are changed

can only be applied in PDE LHS

CopyOperator(copy_from, factor) -> CopyOperator

copies entries from TargetVector to rhs block

can only be applied in PDE RHS