Poisson API

poissonMatrix(problem::Problem; coefficient::Union{Number,ScalarField}=1.0)
∫∇oN_c_∇oN_dΩ(problem::Problem; coefficient::Union{Number,ScalarField}=1.0)

Assembles the global stiffness matrix for a Poisson-type diffusion operator

K_ab = ∫_Ω (∇N_a · ∇N_b) α(x) dΩ

coefficient can be a constant (Number) or an elementwise ScalarField (interpolated to Gauss points using the Lagrange basis, as in the original implementation).

reactionMatrix(problem::Problem; coefficient::Union{Number,ScalarField}=1.0)
∫N_c_N_dΩ(problem::Problem; coefficient::Union{Number,ScalarField}=1.0)

Assembles the global (weighted) mass / reaction matrix

M_ab = ∫_Ω N_a N_b c(x) dΩ

coefficient can be a constant (Number) or an elementwise ScalarField, interpolated to Gauss points using the Lagrange basis (same mechanism as stiffness).

advectionMatrix(problem::Problem; 
                coefficient::Union{Number,ScalarField}=1.0, 
                dir::Int = 1)

aliases:

∫N_c_∂N∂x_dΩ(problem::Problem; coefficient::Union{Number,ScalarField}=1.0)
∫N_c_∂N∂y_dΩ(problem::Problem; coefficient::Union{Number,ScalarField}=1.0)
∫N_c_∂N∂z_dΩ(problem::Problem; coefficient::Union{Number,ScalarField}=1.0)

Assembles the global convection (advection/drift) matrix

C_ab = ∫_Ω N_a (∂N_b/∂x_dir) β(x) dΩ

• dir=1 corresponds to x, dir=2 to y, dir=3 to z (when available).
• coefficient can be a constant (Number) or an elementwise ScalarField, interpolated to Gauss points as in the stiffness implementation.

gradDivMatrix(problem::Problem; coefficient::Union{Number,ScalarField}=1.0)
∫∇N_c_∇N_dΩ(problem::Problem; coefficient::Union{Number,ScalarField}=1.0)

Assembles the global matrix corresponding to the grad-div bilinear form

G_ab = ∫_Ω (∇·u_h) (∇·v_h) α(x) dΩ

This is the weak form associated with the operator ∇(∇·u) (up to sign conventions), commonly appearing in linear elasticity and in grad-div stabilization.

Notes

• Requires problem.pdim == problem.dim (vector unknown with one component per spatial dimension).
• coefficient can be a constant (Number) or an elementwise ScalarField, interpolated to Gauss points using the Lagrange basis (same mechanism as in poissonMatrix).

symmetricGradientMatrix(problem::Problem; coefficient::Union{Number,ScalarField}=1.0)

Assembles the global matrix for the symmetric-gradient (strain) bilinear form

A(u,v) = ∫_Ω 2μ ε(u) : ε(v) dΩ
ε(u)   = 1/2 (∇u + ∇uᵀ)

Notes

• Requires problem.pdim == problem.dim.
• coefficient is typically the shear modulus μ (constant or ScalarField).

curlCurlMatrix(problem::Problem; coefficient::Union{Number,ScalarField}=1.0)
∫∇xN_c_∇xN_dΩ(problem::Problem; coefficient::Union{Number,ScalarField}=1.0)

Assembles the global curl–curl matrix

C(u,v) = ∫_Ω (∇×u) · (∇×v) α(x) dΩ

Notes

• Requires problem.pdim == problem.dim (vector field).
• Uses standard Lagrange H¹ elements (NOT H(curl)-conforming).
• Intended for operator studies, stabilization terms, and educational use.
• Not suitable as a primary operator for Maxwell-type problems.

gradMatrix(problem_u::Problem, problem_p::Problem)

Assembles the mixed gradient matrix G mapping a scalar field (pressure) to a vector field (velocity), corresponding to the weak form

(G p, v) = -∫_Ω p ∇·v dΩ

Notes

• Requires problemu.pdim == problemu.dim (vector field).
• Requires problem_p.pdim == 1 (scalar field).
• Requires problemu.dim == problemp.dim.
• The divergence matrix is obtained as D = -G'.

navierStokesAdvectionMatrix(problem::Problem, u;
                            time_index::Int = 1)

Assembles the (linearized) Navier–Stokes convection operator in Oseen/Picard form:

C(w,v) = ∫_Ω (u·∇)w · v dΩ

where u is a given velocity field (the coefficient), w is the unknown velocity, and v is the test function.

Matrix size:

• problem.pdim == problem.dim is required (vector field).
• Returns a SystemMatrix of size (dimN) × (dimN).

Arguments

• problem::Problem: Must represent the velocity unknown (vector field).
• u: Given advection velocity in the SAME dof ordering as VectorField uses: [u1x,u1y(,u1z), u2x,u2y(,u2z), ...]. Accepted forms:
  • VectorField (nodal, same mesh)
  • AbstractVector{<:Real} (length = dim*problem.non)

Keyword arguments

• time_index: If u stores multiple time steps internally (e.g. as a matrix), you can select which column to use. For plain vectors it is ignored.

Notes

• This is NOT the scalar advectionMatrix. This operator uses the given velocity u inside the Gauss integration, assembled in ONE pass (no splitting to 3 scalar fields).
• The resulting operator is block-diagonal by velocity components (no component mixing), consistent with (u·∇) acting component-wise in Cartesian coordinates.

tensorLaplaceMatrix(problem::Problem; 
                    coefficient::Union{Number,ScalarField}=1.0)

Assembles the tensor Laplace operator

∫_Ω ∇(sym σ) : ∇(sym τ) dΩ

for a nodal TensorField with pdim = 9.

Notes

• Operator is restricted to the symmetric tensor subspace.
• coefficient may be constant or an elementwise ScalarField.

traceLaplaceMatrix(problem::Problem; 
                   coefficient::Union{Number,ScalarField}=1.0)

Assembles the trace Laplace operator

∫_Ω ∇tr(σ) · ∇tr(τ) dΩ

for a nodal TensorField with pdim = 9.

beltramiMichellMatrix(problem::Problem; 
coeff_laplace::Union{Number,ScalarField}=1.0, 
coeff_trace::Union{Number,ScalarField}=1.0)

Assemble the Beltrami-Michell operator as tensorLaplaceMatrix + traceLaplaceMatrix.

Arguments

• problem::Problem: Tensor-field problem definition.
• coeff_laplace: Coefficient for tensorLaplaceMatrix.
• coeff_trace: Coefficient for traceLaplaceMatrix.

Returns

• SystemMatrix: Sum of the two assembled operators.

Example

Kbm = beltramiMichellMatrix(problem; coeff_laplace=1.0, coeff_trace=1.0)

tensorDivDivMatrix(problem::Problem; 
                   coefficient::Union{Number,ScalarField}=1.0)

Assembles the tensor div–div matrix

D(σ,τ) = ∫_Ω (∇·σ_h) · (∇·τ_h) α(x) dΩ

Notes

• Requires problem.pdim == dim^2 (second-order tensor field).
• Acts as an equilibrium-enforcing operator in stress-based formulations (e.g. Beltrami–Michell).
• coefficient may be a constant (Number) or an elementwise ScalarField.

loadTensor(problem::Problem;
           source::Union{Matrix{Float64},TensorField} = zeros(3,3))

Assembles an L2 right-hand-side vector for tensor-valued problems (pdim = 9):

f_{a,α} = ∫_Ω N_a(x) * S_α(x) dΩ

where S is either

• a constant 3×3 tensor (Matrix{Float64}), or
• a nodal TensorField.

Notes

• Intended for Beltrami–Michell and other stress-based formulations.
• No traction, pressure, or surface force interpretation.

loadVector(problem, loads; F=nothing)

Assembles the right-hand-side vector associated with external mechanical loads in weak-form finite element problems.

Depending on the problem dimension and the dimension of the physical group, this function handles concentrated forces, line/surface/volume loads, tractions, pressures, and body forces in 1D, 2D, and 3D.

The interpretation of the physical group dimension is:

• 1D problems:
  • Point → concentrated force
  • Line → distributed force
• 2D problems:
  • Point → concentrated force
  • Line → surface force
  • Surface → body force
• 3D problems:
  • Point → concentrated force
  • Line → edge force
  • Surface → surface force
  • Volume → body force

Load components may be specified as constants, spatial functions, or ScalarFields (nodal type). Axisymmetric and other weighted formulations can be handled by including the appropriate geometric factor (e.g. 2πr) in the load definition or coefficient.

Returns a VectorField for vector-valued problems (pdim > 1) and a ScalarField for scalar problems (pdim == 1).

If F is the deformation gradient (a TensorField), then the function solves the follower load for large displacement problems (in 3D).

Return: loadVec

Types:

• problem: Problem
• loads: Vector{BoundaryCondition}
• loadVec: VectorField
• F: Union{Nothing,TensorField}