GitHub

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 — Type
BilinearOperator(
    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 to u_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.

source
ExtendableFEM.BilinearOperator — Type
function 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 to oa_test.

Keyword Arguments

  • bonus_quadorder: additional quadrature order added to quadorder. Default: 0

  • entities: assemble operator on these grid entities (default = ONCELLS). Default: ONCELLS

  • entry_tolerance: threshold to add entry to sparse matrix. Default: 0

  • factor: factor that should be multiplied during assembly. Default: 1

  • lump: diagonal lumping (=0 no lumping, =1 only keep diagonal entry, =2 accumulate full row to diagonal). Default: 0

  • name: name for operator used in printouts. Default: ''BilinearOperator''

  • parallel: assemble operator in parallel using colors/partitions information (assembles into full matrix directly). Default: false

  • parallel_groups: assemble operator in parallel using CellAssemblyGroups (assembles separated matrices that are added together sequantially). Default: false

  • params: array of parameters that should be made available in qpinfo argument of kernel function. Default: nothing

  • quadorder: 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: false

  • time_dependent: operator is time-dependent ?. Default: false

  • transposed_copy: assemble a transposed copy of that operator into the transposed matrix block(s), 0 = no, 1 = symmetric, -1 = skew-symmetric. Default: 0

  • use_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)
source
ExtendableFEM.BilinearOperator — Method
BilinearOperator(
    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: 0

  • entities: assemble operator on these grid entities (default = ONCELLS). Default: ONCELLS

  • entry_tolerance: threshold to add entry to sparse matrix. Default: 0

  • factor: factor that should be multiplied during assembly. Default: 1

  • lump: diagonal lumping (=0 no lumping, =1 only keep diagonal entry, =2 accumulate full row to diagonal). Default: 0

  • name: name for operator used in printouts. Default: ''BilinearOperator''

  • parallel: assemble operator in parallel using colors/partitions information (assembles into full matrix directly). Default: false

  • parallel_groups: assemble operator in parallel using CellAssemblyGroups (assembles separated matrices that are added together sequantially). Default: false

  • params: array of parameters that should be made available in qpinfo argument of kernel function. Default: nothing

  • quadorder: 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: false

  • time_dependent: operator is time-dependent ?. Default: false

  • transposed_copy: assemble a transposed copy of that operator into the transposed matrix block(s), 0 = no, 1 = symmetric, -1 = skew-symmetric. Default: 0

  • use_sparsity_pattern: read sparsity pattern of jacobian of kernel to find out which components couple. Default: ''auto''

  • verbosity: verbosity level. Default: 0

source
ExtendableFEM.assemble! — Method
assemble!(
    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 to nothing.

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.).
source

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 — Type
function 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 to oa_test.
  • kwargs...: Additional keyword arguments for operator assembly (see below).

Keyword arguments

  • bonus_quadorder: additional quadrature order added to quadorder. Default: 0

  • callback!: function with interface (A, b, sol) that is called in each assembly step. Default: nothing

  • entities: assemble operator on these grid entities (default = ONFACES). Default: ONFACES

  • entry_tolerance: threshold to add entry to sparse matrix. Default: 0

  • factor: factor that should be multiplied during assembly. Default: 1

  • lump: lump the operator (= only assemble the diagonal). Default: false

  • name: name for operator used in printouts. Default: ''BilinearOperatorDG''

  • parallel: assemble operator in parallel using colors/partitions information (assembles into full matrix directly). Default: false

  • parallel_groups: assemble operator in parallel using CellAssemblyGroups (assembles separated matrices that are added together sequantially). Default: false

  • params: array of parameters that should be made available in qpinfo argument of kernel function. Default: nothing

  • quadorder: 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: false

  • time_dependent: operator is time-dependent ?. Default: false

  • transposed_copy: assemble a transposed copy of that operator into the transposed matrix block(s), 0 = no, 1 = symmetric, -1 = skew-symmetric. Default: 0

  • use_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)
source
ExtendableFEM.BilinearOperatorDG — Method
function 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: 0

  • entities: assemble operator on these grid entities (default = ONCELLS). Default: ONCELLS

  • entry_tolerance: threshold to add entry to sparse matrix. Default: 0

  • factor: factor that should be multiplied during assembly. Default: 1

  • lump: diagonal lumping (=0 no lumping, =1 only keep diagonal entry, =2 accumulate full row to diagonal). Default: 0

  • name: name for operator used in printouts. Default: ''BilinearOperator''

  • parallel: assemble operator in parallel using colors/partitions information (assembles into full matrix directly). Default: false

  • parallel_groups: assemble operator in parallel using CellAssemblyGroups (assembles separated matrices that are added together sequantially). Default: false

  • params: array of parameters that should be made available in qpinfo argument of kernel function. Default: nothing

  • quadorder: 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: false

  • time_dependent: operator is time-dependent ?. Default: false

  • transposed_copy: assemble a transposed copy of that operator into the transposed matrix block(s), 0 = no, 1 = symmetric, -1 = skew-symmetric. Default: 0

  • use_sparsity_pattern: read sparsity pattern of jacobian of kernel to find out which components couple. Default: ''auto''

  • verbosity: verbosity level. Default: 0

source

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))