Poisson-type operators

advectionMatrix(problem::Problem; coefficient::Union{Number,ScalarField}=1.0, dir::Int=1)
∫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.

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.

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

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'.

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.

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.

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

solveField(K, f; support=BoundaryCondition[], iterative=false,
           reltol=sqrt(eps()), maxiter=..., preconditioner=Identity(),
           ordering=true)

Solves a linear static field problem with Dirichlet boundary conditions.

The algebraic system

K * u = f

is solved for the unknown field u, where prescribed values are enforced solver-side via boundary conditions.

The solution is obtained by partitioning the degrees of freedom into free and constrained sets. On constrained DOFs the values prescribed by support override the solution, while on free DOFs the reduced system is solved:

K_ff * u_f = f_f − K_fc * u_c

The function supports both direct and iterative linear solvers and works uniformly for ScalarField and VectorField unknowns.

Arguments

• K::SystemMatrix Global system matrix.
• f::Union{ScalarField,VectorField} Right-hand-side vector.
• support::Vector{BoundaryCondition} (keyword, default = BoundaryCondition[]) Dirichlet-type boundary conditions.
• iterative::Bool (keyword, default = false) If true, use an iterative solver (conjugate gradient).
• reltol::Real (keyword, default = sqrt(eps())) Relative tolerance for the iterative solver.
• maxiter::Int (keyword, default = K.model.non * K.model.dim) Maximum number of iterations for the iterative solver.
• preconditioner (keyword, default = Identity()) Preconditioner for the iterative solver.
• ordering::Bool (keyword, default = true) If false, use an explicit LU factorization without reordering.

Returns

• u::Union{ScalarField,VectorField} Solution field with prescribed values enforced on constrained DOFs.

Notes

• Boundary conditions are applied inside the solver; the input field f is not modified.
• The algorithm itself is agnostic to the physical meaning of the field (scalar, vector, tensor), as long as K and f are dimensionally consistent.
• When iterative = true, the system is solved using conjugate gradient on the reduced matrix K_ff.

sourceVector(problem, sources)

Assembles the right-hand-side vector corresponding to a volumetric source term in Poisson-type and other scalar or vector-valued PDEs expressed in weak form.

Mathematically, this function assembles

f_a = ∫_Ω N_a f dΩ

where the source term f may be given as a constant, a spatial function, or a ScalarField.

This function is an alias of loadVector and shares the same implementation. The difference is purely interpretational: sourceVector emphasizes the role of the right-hand side as a PDE source term rather than a mechanical load.

Axisymmetric or spherically symmetric problems can be handled by including the appropriate geometric weighting (e.g. 2πr, 4πr²) as a coefficient field.

In LowLevelFEM, right-hand-side vectors are assembled independently of the governing equation. The same numerical machinery can therefore represent mechanical loads, heat sources, or generic PDE source terms.

Returns a ScalarField or VectorField, depending on the problem field dimension.

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

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.

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.