BilinearOperator
BilinearOperator
provides a flexible interface for defining bilinear forms (matrices) in finite element problems. These operators typically represent (linearizations of) PDE operators, stabilization terms, or custom user-defined couplings. The interface supports both standard and discontinuous Galerkin (DG) forms, and allows for custom assembly on cells, faces, or other grid entities.
Overview
- Standard Bilinear Operators: Assemble contributions to the system matrix, including diffusion, mass, and convection terms, as well as custom couplings. Supports both cell- and face-based assembly, including jump terms for conforming and piecewise spaces that only require the dofs on that entity.
- Discontinuous Galerkin (DG) Operators: Use
BilinearOperatorDG
for forms that require evaluation of jumps or averages involving all degrees of freedom on neighboring cells (e.g., interior penalty stabilization, gradient jumps for H1 elements). - Custom Matrices: You can also assign a user-assembled matrix as a
BilinearOperator
.
Constructors
ExtendableFEM.BilinearOperator
— TypeBilinearOperator(
A::AbstractMatrix,
u_test::Vector{<:Union{Int64, Unknown}};
...
) -> Union{ExtendableFEM.BilinearOperatorFromMatrix{_A, <:AbstractMatrix{T}} where T, ExtendableFEM.BilinearOperatorFromMatrix{_A, MT} where {T, MT<:AbstractMatrix{T}}} where _A<:Union{Integer, Unknown}
BilinearOperator(
A::AbstractMatrix,
u_test::Vector{<:Union{Int64, Unknown}},
u_ansatz;
...
) -> Union{ExtendableFEM.BilinearOperatorFromMatrix{_A, <:AbstractMatrix{T}} where T, ExtendableFEM.BilinearOperatorFromMatrix{_A, MT} where {T, MT<:AbstractMatrix{T}}} where _A<:Union{Integer, Unknown}
BilinearOperator(
A::AbstractMatrix,
u_test::Vector{<:Union{Int64, Unknown}},
u_ansatz,
u_args;
kwargs...
) -> Union{ExtendableFEM.BilinearOperatorFromMatrix{_A, <:AbstractMatrix{T}} where T, ExtendableFEM.BilinearOperatorFromMatrix{_A, MT} where {T, MT<:AbstractMatrix{T}}} where _A<:Union{Integer, Unknown}
Constructs a bilinear operator from a user-supplied matrix A
, which can be a sparse matrix or a block-structured FEMatrix. The arguments u_test
and u_ansatz
specify the test and ansatz (trial) unknowns or indices, determining where the (blocks of the) matrix are inserted in the global system. If u_ansatz
is not provided, it defaults to u_test
.
Arguments
A::AbstractMatrix
: The matrix representing the bilinear form, e.g., a sparse matrix.u_test::Vector{<:Union{Unknown, Int}}
: Identifiers or indices for the test functions.u_ansatz::Vector{<:Union{Unknown, Int}}
(optional): Identifiers or indices for the ansatz (trial) functions. Defaults tou_test
.u_args::Vector{<:Union{Unknown, Int}}
(optional): Identifiers or indices for unknowns on which the matrix depends in a nonlinear way (this tells the solver which solution blocks trigger reassembly).
Returns
A BilinearOperatorFromMatrix
object that can be used in the assembly process.
ExtendableFEM.BilinearOperator
— Typefunction BilinearOperator(
[kernel!::Function],
oa_test::Array{<:Tuple{Union{Unknown,Int}, DataType},1},
oa_ansatz::Array{<:Tuple{Union{Unknown,Int}, DataType},1} = oa_test;
kwargs...)
Constructs a bilinear form by evaluating the vector product of operator evaluations of the test and ansatz functions. Operator evaluations are specified as tuples pairing an unknown identifier (or integer) with a function operator (such as those generated by grad(u)
, id(u)
, etc).
If a kernel function is provided as the first argument, it customizes the ansatz function evaluation. The kernel must have the signature:
kernel!(result, eval_ansatz, qpinfo)
where result
is the output vector, eval_ansatz
is the ansatz function evaluation at the quadrature point, and qpinfo
provides quadrature point information.
If no kernel is provided, the standard kernel is used that produces the dot product.
Arguments
oa_test
: Array of tuples(unknown, operator)
for test functions.oa_ansatz
: Array of tuples(unknown, operator)
for ansatz (trial) functions. Defaults tooa_test
.
Keyword Arguments
bonus_quadorder
: additional quadrature order added to quadorder. Default: 0entities
: assemble operator on these grid entities (default = ONCELLS). Default: ONCELLSentry_tolerance
: threshold to add entry to sparse matrix. Default: 0factor
: factor that should be multiplied during assembly. Default: 1lump
: diagonal lumping (=0 no lumping, =1 only keep diagonal entry, =2 accumulate full row to diagonal). Default: 0name
: name for operator used in printouts. Default: ''BilinearOperator''parallel
: assemble operator in parallel using colors/partitions information (assembles into full matrix directly). Default: falseparallel_groups
: assemble operator in parallel using CellAssemblyGroups (assembles separated matrices that are added together sequantially). Default: falseparams
: array of parameters that should be made available in qpinfo argument of kernel function. Default: nothingquadorder
: quadrature order. Default: ''auto''regions
: subset of regions where operator should be assembly only. Default: Any[]store
: store matrix separately (and copy from there when reassembly is triggered). Default: falsetime_dependent
: operator is time-dependent ?. Default: falsetransposed_copy
: assemble a transposed copy of that operator into the transposed matrix block(s), 0 = no, 1 = symmetric, -1 = skew-symmetric. Default: 0use_sparsity_pattern
: read sparsity pattern of jacobian of kernel to find out which components couple. Default: ''auto''verbosity
: verbosity level. Default: 0
Example
BilinearOperator([grad(1)], [grad(1)]; entities=ON_CELLS)
ExtendableFEM.BilinearOperator
— MethodBilinearOperator(
kernel,
oa_test::Vector{<:Tuple{Union{Int64, Unknown}, DataType}},
oa_ansatz::Vector{<:Tuple{Union{Int64, Unknown}, DataType}},
oa_args::Vector{<:Tuple{Union{Int64, Unknown}, DataType}};
kwargs...
) -> BilinearOperator{Float64}
Constructs a nonlinear bilinear form by evaluating a user-supplied kernel function that depends on both the operator evaluations of the ansatz (trial) functions and additional argument functions (e.g., the current solution). The result of the kernel is then contracted with the operator evaluations of the test functions. This is typically used for linearizations of nonlinear operators.
The kernel function must have the following signature:
kernel!(result, eval_ansatz, eval_args, qpinfo)
where
result
: Output vector for the kernel evaluation.eval_ansatz
: Operator evaluation(s) of the ansatz functions at the quadrature point.eval_args
: Operator evaluation(s) of the argument functions (e.g., current solution) at the quadrature point.qpinfo
: Quadrature point information structure.
Operator evaluations are specified as tuples pairing an unknown identifier (or integer) with a function operator (such as those generated by grad(u)
, id(u)
, etc).
Arguments
kernel
: User-supplied kernel function as described above.oa_test
: Array of tuples(unknown, operator)
for test functions.oa_ansatz
: Array of tuples(unknown, operator)
for ansatz (trial) functions.oa_args
: Array of tuples(unknown, operator)
for argument functions (e.g., current solution).kwargs...
: Additional keyword arguments for operator assembly (see below).
Keyword Arguments
bonus_quadorder
: additional quadrature order added to quadorder. Default: 0entities
: assemble operator on these grid entities (default = ONCELLS). Default: ONCELLSentry_tolerance
: threshold to add entry to sparse matrix. Default: 0factor
: factor that should be multiplied during assembly. Default: 1lump
: diagonal lumping (=0 no lumping, =1 only keep diagonal entry, =2 accumulate full row to diagonal). Default: 0name
: name for operator used in printouts. Default: ''BilinearOperator''parallel
: assemble operator in parallel using colors/partitions information (assembles into full matrix directly). Default: falseparallel_groups
: assemble operator in parallel using CellAssemblyGroups (assembles separated matrices that are added together sequantially). Default: falseparams
: array of parameters that should be made available in qpinfo argument of kernel function. Default: nothingquadorder
: quadrature order. Default: ''auto''regions
: subset of regions where operator should be assembly only. Default: Any[]store
: store matrix separately (and copy from there when reassembly is triggered). Default: falsetime_dependent
: operator is time-dependent ?. Default: falsetransposed_copy
: assemble a transposed copy of that operator into the transposed matrix block(s), 0 = no, 1 = symmetric, -1 = skew-symmetric. Default: 0use_sparsity_pattern
: read sparsity pattern of jacobian of kernel to find out which components couple. Default: ''auto''verbosity
: verbosity level. Default: 0
ExtendableFEM.assemble!
— Methodassemble!(
A::FEMatrix,
O::BilinearOperator{Tv, UT};
...
) -> Any
assemble!(
A::FEMatrix,
O::BilinearOperator{Tv, UT},
sol;
assemble_matrix,
kwargs...
) -> Any
Assembles the given BilinearOperator
into the provided FEMatrix
A
. This function computes the matrix representation of the bilinear form defined by O
and inserts the resulting blocks into A
. If the operator depends on additional argument functions (e.g., for nonlinear problems), these are provided via the sol
argument.
Arguments
A::FEMatrix
: The finite element matrix to assemble into.O::BilinearOperator
: The bilinear operator to assemble.sol
: (Optional) Array of solution vectors or argument blocks required for nonlinear operators. Defaults tonothing
.
Keyword Arguments
assemble_matrix::Bool
: Whether to assemble the matrix (default:true
).kwargs...
: Additional keyword arguments passed to the assembly routines (e.g., for setting the time, verbosity, etc.).
BilinearOperatorDG
BilinearOperatorDG
is intended for bilinear forms that involve jumps or averages of discontinuous quantities on faces, requiring access to all degrees of freedom on neighboring cells. This is essential for DG methods and certain stabilization techniques.
ExtendableFEM.BilinearOperatorDG
— Typefunction BilinearOperatorDG(
[kernel!::Function],
oa_test::Array{<:Tuple{Union{Unknown,Int}, DataType},1},
oa_ansatz::Array{<:Tuple{Union{Unknown,Int}, DataType},1} = oa_test;
kwargs...)
Generates a bilinear form that evaluates the vector product of the (discontinuous) operator evaluation(s) of the test function(s) with the (discontinuous) operator evaluation(s) of the ansatz function(s). If a function is provided in the first argument, the ansatz function evaluations can be customized by the kernel function, and its result vector is then used in a dot product with the test function evaluations. In this case, the header of the kernel function must conform to the interface:
kernel!(result, eval_ansatz, qpinfo)
where qpinfo
allows access to information at the current quadrature point.
Operator evaluations are tuples that pair an unknown identifier or integer with a function operator (such as those generated by jump(grad(u))
, id(u)
, etc).
Arguments
oa_test
: Array of tuples(unknown, operator)
for test functions.oa_ansatz
: Array of tuples(unknown, operator)
for ansatz (trial) functions. Defaults tooa_test
.kwargs...
: Additional keyword arguments for operator assembly (see below).
Keyword arguments
bonus_quadorder
: additional quadrature order added to quadorder. Default: 0callback!
: function with interface (A, b, sol) that is called in each assembly step. Default: nothingentities
: assemble operator on these grid entities (default = ONFACES). Default: ONFACESentry_tolerance
: threshold to add entry to sparse matrix. Default: 0factor
: factor that should be multiplied during assembly. Default: 1lump
: lump the operator (= only assemble the diagonal). Default: falsename
: name for operator used in printouts. Default: ''BilinearOperatorDG''parallel
: assemble operator in parallel using colors/partitions information (assembles into full matrix directly). Default: falseparallel_groups
: assemble operator in parallel using CellAssemblyGroups (assembles separated matrices that are added together sequantially). Default: falseparams
: array of parameters that should be made available in qpinfo argument of kernel function. Default: nothingquadorder
: quadrature order. Default: ''auto''regions
: subset of regions where operator should be assembly only. Default: Any[]store
: store matrix separately (and copy from there when reassembly is triggered). Default: falsetime_dependent
: operator is time-dependent ?. Default: falsetransposed_copy
: assemble a transposed copy of that operator into the transposed matrix block(s), 0 = no, 1 = symmetric, -1 = skew-symmetric. Default: 0use_sparsity_pattern
: read sparsity pattern of jacobian of kernel to find out which components couple. Default: ''auto''verbosity
: verbosity level. Default: 0
Example
BilinearOperatorDG([jump(grad(1))], [jump(grad(1))]; entities=ON_FACES)
ExtendableFEM.BilinearOperatorDG
— Methodfunction BilinearOperatorDG(
kernel::Function,
oa_test::Array{<:Tuple{Union{Unknown,Int}, DataType},1},
oa_ansatz::Array{<:Tuple{Union{Unknown,Int}, DataType},1},
oa_args::Array{<:Tuple{Union{Unknown,Int}, DataType},1};
kwargs...)
Generates a nonlinear bilinear form that evaluates a kernel function that depends on the (discontinuou) operator evaluation(s) of the ansatz function(s) and the (discontinuous) operator evaluations of the current solution. The result of the kernel function is used in a vector product with the operator evaluation(s) of the test function(s). Hence, this can be used as a linearization of a nonlinear operator. The header of the kernel functions needs to be conform to the interface
kernel!(result, eval_ansatz, eval_args, qpinfo)
where qpinfo allows to access information at the current quadrature point.
Operator evaluations are tuples that pair an unknown identifier or integer with a Function operator.
Keyword arguments:
bonus_quadorder
: additional quadrature order added to quadorder. Default: 0entities
: assemble operator on these grid entities (default = ONCELLS). Default: ONCELLSentry_tolerance
: threshold to add entry to sparse matrix. Default: 0factor
: factor that should be multiplied during assembly. Default: 1lump
: diagonal lumping (=0 no lumping, =1 only keep diagonal entry, =2 accumulate full row to diagonal). Default: 0name
: name for operator used in printouts. Default: ''BilinearOperator''parallel
: assemble operator in parallel using colors/partitions information (assembles into full matrix directly). Default: falseparallel_groups
: assemble operator in parallel using CellAssemblyGroups (assembles separated matrices that are added together sequantially). Default: falseparams
: array of parameters that should be made available in qpinfo argument of kernel function. Default: nothingquadorder
: quadrature order. Default: ''auto''regions
: subset of regions where operator should be assembly only. Default: Any[]store
: store matrix separately (and copy from there when reassembly is triggered). Default: falsetime_dependent
: operator is time-dependent ?. Default: falsetransposed_copy
: assemble a transposed copy of that operator into the transposed matrix block(s), 0 = no, 1 = symmetric, -1 = skew-symmetric. Default: 0use_sparsity_pattern
: read sparsity pattern of jacobian of kernel to find out which components couple. Default: ''auto''verbosity
: verbosity level. Default: 0
Examples
Example: Stokes Operator
For the linear operator of a Stokes problem, a kernel could look like:
μ = 0.1 # viscosity parameter
function kernel!(result, input, qpinfo)
∇u, p = view(input,1:4), view(input, 5)
result[1] = μ*∇u[1] - p[1]
result[2] = μ*∇u[2]
result[3] = μ*∇u[3]
result[4] = μ*∇u[4] - p[1]
result[5] = -(∇u[1] + ∇u[4])
return nothing
end
The corresponding BilinearOperator
constructor call reads:
u = Unknown("u"; name = "velocity")
p = Unknown("p"; name = "pressure")
BilinearOperator(kernel!, [grad(u), id(p)]; use_sparsity_pattern = true)
The use_sparsity_pattern
argument ensures that the zero pressure-pressure block of the matrix is not assembled, since input[5]
does not couple with result[5]
.
Example: Interior Penalty Stabilization
A popular stabilization for convection-dominated problems is based on jumps of the gradient, which can be realized with the following kernel:
function stab_kernel!(result, input, qpinfo)
result .= input .* qpinfo.volume^2
end
and the BilinearOperatorDG
constructor call:
u = Unknown("u")
assign_operator!(PD, BilinearOperatorDG(stab_kernel!, [jump(grad(u))]; entities = ON_IFACES, factor = 0.01))