Internals

GeometricFactors(master, mesh; T=Float64)

One-time precomputation of all mesh geometry at quadrature points, shared by every discretization (the DG residual kernels consume it directly as DGContext; the HDG/CG caches compose it): face/element connectivity resolved to plain index arrays (no runtime findfirst), face and element geometry evaluated at quadrature points, and explicit inverse mass matrices. All fields are plain arrays of eltype T (Int32 for indices), so the whole cache moves to a GPU with Adapt.adapt(CuArray, gf).

Faces 1:ni are interior, ni+1:nf are boundary (the ordering of mesh.f).

Field shapes (npl volume nodes, npf face nodes, ng/ngf quadrature points):

• facecon (npf, 2, nf): volume-node indices of face nodes; side 1 = left element, side 2 = right (unused for boundary faces).
• f_el (nf, 2): (left element, right element); column 2 is -ib (negative boundary tag) for boundary faces, as in mesh.f.
• nlg (ngf, Dim, nf), dws (ngf, nf), pfg (ngf, Dim, nf): outward unit normal (w.r.t. left element), weighted measure, and physical coordinates at face quadrature points.
• vol ::VolumeTables: volume geometry in the compact affine/curved split layout (dense per-element tables only for curved elements).
• Minv (npl, npl, nt): inverse element mass matrices (dense for every element — same cost class as the solution itself).
• shapf (npf, ngf), shap (npl, ng): shape-function values (shared).

The pre-split dense tables remain available as materializing properties: gf.shapd (npl, ng, Dim, nt), gf.wjac (ng, nt), gf.pg (ng, Dim, nt) allocate and fill the full array on access — convenient for diagnostics (sum(ctx.wjac)) and one-time CPU consumers, wasteful inside loops.

Deprecated property aliases (one release, 2D): shapx/shapy materialize shapd[:, :, 1/2, :]; sh1d, np1d, ng1d read shapf, npf, ngf.

RefTables(master)

Reference-element tabulations shared by all geometry evaluations: shape values/derivatives at volume and face quadrature points, quadrature weights, weighted derivative tables and the reference mass matrix. Dimension-generic: the face tables come from the master element's recursive face element.

SideGeometry(master, mesh; T=Float64)

Face geometry in element-local orientation, one entry per (local side, element) — what the HDG trace assembly needs (its trace basis follows the element's own canonical traversal perm[:, s, 1], not the global face's stored direction):

• nl (ngf, Dim, Dim+1, nt): outward unit normal of side s of element e,
• sw (ngf, Dim+1, nt): quadrature-weighted side measure,
• pfs (ngf, Dim, Dim+1, nt): physical coordinates of the side quadrature points.

Evaluated isoparametrically from the side's high-order nodes, so curved faces are exact to the geometric order.

VolumeTables

Volume geometry at quadrature points in the memory-compact split layout: dense per-element tables are stored only for curved elements; straight (affine) elements — the bulk of any real mesh — store just their constant Jacobian data and share one set of reference tables. This is what keeps 3D meshes affordable: the dense shapd table alone is O(npl·ng·Dim) ≈ 165 KB per element at p = 3 on tets, while the affine representation is O(Dim²).

• curved_ix (nt,): 0 for affine elements, else the element's column in the dense tables below.
• cshapd (npl, ng, Dim, ntc), cwjac (ng, ntc), cpg (ng, Dim, ntc): dense quadrature-/Jacobian-weighted derivative tables, weighted Jacobians, and quadrature-point coordinates of the ntc curved elements.
• aC, aJ (Dim, Dim, nt), av0 (Dim, nt), adetJ (nt,): per-element adjugate, Jacobian, map origin (first vertex), and det J of the affine map x(ξ) = av0 + aJ ξ (zero for curved elements).
• rshapdg (npl, ng, Dim), rgpts (ng, Dim), rgwgh (ng,): shared reference tables (quadrature-weighted reference derivatives, quadrature points and weights).

Kernels branch per element on curved_ix (see the accessors quad_coords, quad_weight and the flux-rotation pattern in the DG kernels); the trade of a few extra FLOPs for O(npl·ng) less memory traffic per affine element is a win on both backends.

adjugate(J) -> SMatrix

Adjugate (transposed cofactor matrix) of a Dim × Dim Jacobian: adjugate(J) == det(J) * inv(J) without the division — the exact-arithmetic form the affine geometry path contracts reference derivatives with.

affine_jacobian(verts, ::Val{Dim}) -> (J, detJ, C)

Constant Jacobian data of the affine map from the reference simplex to the element with vertex coordinates verts (Dim+1, Dim): J[d, k] = ∂x_d/∂ξ_k, its determinant, and its adjugate C = detJ * inv(J). The one definition shared by the curved-element evaluator and the compact affine storage, so both produce bit-identical values.

element_geometry!(shapd, wjac, pg, rt, coords; verts=nothing) -> M

Fill, for one element with high-order node coordinates coords (npl, Dim):

• shapd (npl, ng, Dim): quadrature- and Jacobian-weighted physical derivative tables (∫ ∂φ/∂x_d ⋅ f becomes shapd[:, :, d] * f(quad)),
• wjac (ng,): gwgh .* detJ,
• pg (ng, Dim): physical coordinates of the volume quadrature points,

and return the element mass matrix M (npl, npl).

For a straight (affine) element pass verts, the (Dim+1, Dim) vertex coordinate matrix; the Jacobian is then constant and M is the scaled reference mass. Otherwise the isoparametric map is evaluated per quadrature point (curved element).

face_geometry!(nlg, dws, pfg, rt, coords; tangent=nothing)

Fill, for one face with high-order node coordinates coords (npf, Dim) (ordered left-element-outward, i.e. dgnodes[perml, :, el]):

• nlg (ngf, Dim): outward unit normal w.r.t. the left element,
• dws (ngf,): quadrature-weighted face measure,
• pfg (ngf, Dim): physical coordinates of the face quadrature points.

For a straight face pass tangent, an NTuple{Dim-1, SVector{Dim}} of vertex-to-vertex edge vectors; the metric is then constant. Otherwise the metric is evaluated from the high-order nodes (curved face).

face_normal(τ) -> SVector{2}
face_normal(τ₁, τ₂) -> SVector{3}

Outward (unnormalized) face normal from the face tangent vector(s): the quarter-turn rotation of the single tangent in 2D, the cross product of the two tangents in 3D. Its norm is the face measure Jacobian. This is the one place where 2D and 3D geometry genuinely differ.

inscribed_diameter(p, t, it, ::Val{Dim}) -> h

Inscribed-circle (2D) / inscribed-sphere (3D) diameter 2r = 2·Dim·|K|/|∂K| of the vertex simplex of element it — the h that controls explicit CFL limits. p/t are the mesh vertex coordinates and element connectivity.

min_inscribed_diameter(mesh) -> h_min

Smallest inscribed_diameter over all elements of the mesh — the mesh-quality scale of CFL estimates (compute_dt, StepsizeCallback).

quad_coords(vol, g, e, ::Val{Dim}) -> SVector{Dim}

Physical coordinates of volume quadrature point g of element e: av0 + aJ ξ_g for affine elements, a dense-table read for curved ones. Device-inlineable.

quad_weight(vol, g, e) -> T

Quadrature-weighted Jacobian gwgh[g] * detJ at volume quadrature point g of element e. Device-inlineable.

Solution of a CGProblem: u (npl, 1, nt) (DG-numbered for plotting), the discrete energy, and the Krylov iteration count (0 for a direct solve).

Solution of a DGProblem: u (npl, nc, nt) at time t. The callbacks field carries whatever was passed as callback to solve (nothing by default), so callback histories — e.g. an AnalysisCallback's time/data — ride along on the solution.

Solution of an HDGProblem: u (npl, 1, nt), flux q (npl, Dim, nt), trace uhat (nps, nf) and GMRES iteration count (0 for a direct solve).

l2error(sol, exact; component=1)

L2 error of solution component component against exact(x, y).

_dg_physics(prob::DGProblem) -> DGPhysics

Assemble the physics bundle the DG kernels consume: the equation, the per-boundary condition tuple (validated against the mesh's boundary count), the numerical flux, source, and stabilization — problem-level nothings resolved to the equation's defaults.

_ordered_bcs(bc, mesh)

Boundary conditions may be given positionally (a vector/tuple indexed by boundary tag) or, when the mesh generator attached boundary names, as a NamedTuple keyed by those names, e.g. for mkmesh_square: (bottom=Dirichlet(), right=Neumann(), top=Dirichlet(), left=Neumann()).

ndims(mesh) -> Int

Spatial dimension of the mesh (the Dim type parameter): 2 for triangles, 3 for tetrahedra.

mkelcon(t2f, t2o, porder, face_plocal) -> elcon (nps, 4, nt)

3D (tetrahedral) trace connectivity: face f owns the nps = (porder+1)(porder+2)/2 trace nodes (f-1)nps+1 : f*nps in the face's canonical node order (face_plocal, the triangle face element's nodes), and each element couples through orient_perm at its orientation code.

mkelcon(t2f, t2o, porder)

Element-to-global trace-node connectivity elcon (porder+1, 3, nt) from the element-to-face map t2f and orientation codes t2o: face f owns trace nodes (f-1)*(porder+1)+1 : f*(porder+1) in the face's canonical order, and each element traverses them through orient_perm[:, t2o[it, s]].

orient_perm(nps, Val(2))          -> (nps, 2)  permutation table
orient_perm(face_plocal, Val(3))  -> (nps, 6)  permutation table

Trace-node permutations of one face under each orientation code: an element seeing a face with orientation o couples its k-th local trace node (its own canonical face traversal) to the face's global trace node op[k, o]. In 2D this is the identity and the reversal; in 3D it is built constructively from the triangle face element's node barycentrics: op[k, o] is the canonical face node m whose coordinates permuted by the o-th triangle symmetry match node k (face_plocal[m, σ_o] == face_plocal[k, :]).

simpvol(p::Matrix{T}, t::Matrix{Int}) where T<:Real

Compute signed volumes of the simplices of a mesh: areas of triangles (t (nt, 3)) or volumes of tetrahedra (t (nt, 4)), positive for counterclockwise / right-handed vertex ordering.

Parameters:

• p: N×Dim matrix of vertex coordinates
• t: M×(Dim+1) matrix of simplex vertex indices

Returns:

• Vector of signed volumes for each simplex

straight_dgnodes(p, t, plocal) -> dgnodes (npl, 2, nt)

Map the master-element nodes plocal (barycentric, (npl, 3)) affinely into every triangle of (p, t). This is the straight-element part of createnodes; use it when the high-order nodes will be transformed analytically afterwards (e.g. mkmesh_trefftz's conformal maps).

transform_points!(pnew, p, db, dt, H)

Helper function to transform points according to the duct mapping.

unique_with_inverse(arr)

Find unique rows in a matrix and return inverse mapping. Similar to np.unique(arr, return_inverse=True, axis=0) in NumPy.

build_face_perm(::Val{Dim}, plocal, face_plocal) -> perm (npf, Dim+1, norient)

Constructive face-node permutation table: for each local face j and orientation o, perm[:, j, o] lists the volume-node indices whose barycentric coordinates (restricted to the face's vertices) match the face element's nodes face_plocal under the o-th symmetry of the face simplex. No hand-tabulated index magic — works for every porder, provided the volume node set restricted to each face equals the face element's node set (true for the symmetrized localpnts/localpnts3d distributions; asserted).

gaussjacobi(n, α, β=0) -> (x, w)

Gauss–Jacobi points and weights on [-1, 1] with weight (1-x)^α (1+x)^β: n points integrate f(x) (1-x)^α (1+x)^β exactly for polynomial f of degree ≤ 2n-1. Same Vandermonde-solve construction as gaussquad1d (which is the α = β = 0 case); by orthogonality only the constant polynomial has a nonzero weighted moment, 2^(α+β+1) B(α+1, β+1).

apply_blockjacobi(B, v)

Applies a block Jacobi preconditioner to a vector.

Arguments

• B::AbstractArray: Block Jacobi preconditioner with dimensions (ncf, ncf, nf)
• v::AbstractArray: Vector to be preconditioned (can be 1D flattened or 2D array)

Returns

• w::Array: Preconditioned vector in the same format as input v

boundary_face_quad(mesh, master, i)

Quadrature data of boundary face i in global face-node order: the physical quadrature points Xq (ngf × Dim), the weighted measure wds, and the face mass matrix Tm. Used to impose boundary data through its L2 projection P∂g onto Pk(F) — nodal interpolation of the data would introduce an O(h^{k+1}) perturbation that destroys the superconvergence of the method.

compute_blockjacobi(A)

Computes block Jacobi preconditioner for HDG method.

Arguments

• A: Global matrix in dense format with dimensions (ncf, ncf, nbf, nf)

Returns

• B: Block Jacobi preconditioner with dimensions (ncf, ncf, nf)

Face matrix ⟨w μ_a, μ_b⟩_F in the nps face nodes for a quadrature weight w.

Element-local face-scalar values (nfc) of a face field Θ (one DOF per face node).

Element-local face-scalar components (nfc × Dim) of the interleaved trace vector Λ.

Assembles the global matrix and vector for a specific face.

Parameters:

AE : Array Element matrices FE : Array Element vectors f : Array Face to element connectivity t2f : Array Element to face connectivity G : Array (npf, nfe, ne) Element-local -> face-canonical trace-node gather (see hdg_densesystem) ncf : Int Number of components per face nbf : Int Number of neighboring faces nfe : Int Number of faces per element i : Int Face index

Returns:

Ai : Array Global matrix for face i Fi : Array Global vector for face i

hdg_applydbc!(ae, fe, master, mesh, dbc)

Applies the strong Dirichlet boundary condition dbc to the boundary-face rows of the element matrices/vectors ae, fe in place. The identity rows are set in the element's own face-node order, so the assembled system pins each global trace node to the boundary value at its own coordinate for any face orientation (2D or 3D).

hdg_densesystem(AE, FE, f, t2f, elcon)
hdg_densesystem(AE, FE, f, t2f, t2o, npf)

Assembles the global system in dense face-block format.

Arguments

• AE: Element matrices
• FE: Element vectors
• f: Face to element connectivity
• t2f: Element to face connectivity
• elcon (npf, nfe, ne): orientation-resolved element-to-global trace-node connectivity (mkelcon); the second form reconstructs the 2D (two-orientation) elcon from t2o and the trace size npf for backwards compatibility.

Returns

• A: Global matrix in dense format
• F: Global vector in dense format

hdg_elem_face(dg, master, s)

Face metric terms of local face s of the element with nodes dg: the volume node indices ps on the face (canonical traversal), the outward unit normal n (ngf × Dim), the face measure Jacobian ds, and the weighted measure wds at the face quadrature points. The normal comes from face_normal dispatch (tangent rotation in 2D, cross product in 3D).

hdg_elem_volume(dg, master)

Volume metric terms of one (possibly curved) element with nodes dg (npl × Dim): shape values shap (npl × ng), the quadrature- and Jacobian-weighted physical derivative tables shapd (npl × ng × Dim, shapd[m, g, d] = w_g jac_g ∂φ_m/∂x_d), the mass matrix M, the convection-type matrices C (npl × npl × Dim, C[m, n, d] = (φ_m, ∂_d φ_n)), and the weighted Jacobian wjac.

hdg_elemmats(master, mesh, source, dbc, param)

Computes the HDG element matrices/vectors ae (3nps, 3nps, nt), fe (3nps, nt) for all elements (threaded) and applies the strong Dirichlet boundary conditions dbc to the boundary-face rows.

hdg_face_ops(dg, master)

The lifting operators E (npl × nfc × Dim) of the gradient equation, E[m, ℓ, d] = ⟨μ_ℓ, φ_m n_d⟩_∂K, mapping face-scalar trace values to volume test functions.

hdg_localrecovery(master, mesh, uhath, source, param)

Recovers the element-local solution uh (npl, 1, nt) and flux qh (npl, 2, nt) from the global trace vector uhath by solving the local problems (threaded).

hdg_matvec(A, F, f2f)

Performs matrix-vector multiplication for HDG method using face-to-face connectivity.

Arguments

• A: Global matrix in dense format (ncf, ncf, nbf, nf)
• F: Vector to be multiplied (flattened)
• f2f: Face-to-face connectivity

Returns

• v: Result of matrix-vector multiplication (flattened)

hdg_ns_elemmat(dg, master, ν, τ, um, λe, fe, ffun, uolde, dtinv)

Element matrices for one Newton step of the HDG incompressible Navier-Stokes discretization, linearized about the velocity um (npl × Dim) and the trace components λe (nfc × Dim). The body force is either the DG field fe (npl × Dim, integrated exactly) or the function ffun (evaluated at the quadrature points — interpolating it at the nodes instead would spoil the superconvergence of the method); uolde is the velocity at the previous time level for backward Euler (or nothing), and dtinv = 1/Δt (0 for steady state).

The velocity gradient is eliminated analytically and the local (u, p) system is statically condensed. Returns (Ke, Ge, re, crow, Z, area): the condensed trace matrix (Dim·nfc × Dim·nfc, trace ordered [λ₁; …; λDim]), the mean-pressure coupling column, the condensed right-hand side, the element- compatibility row ⟨λ·n, 1⟩∂K, the local solution operator Z = A⁻¹[B bρ r] for the recovery of (u, p), and the element measure.

Recover the velocity gradient Lij = M⁻¹(Ej λi - Cjᵀ u_i) of one element; column layout (i - 1) Dim + j (L11, L12, …, L1Dim, L21, …).

Recover the scalar gradient qd = M⁻¹(Ed θ̂ - C_dᵀ θ) of one element.

Jacobi (diagonal) preconditioner: y = x ./ d.

Global mask of CG nodes on the (Dirichlet) boundary, from the face list.

Specific total enthalpy h = (ρE + p)/ρ.

finish!(cb, state) -> nothing

Lifecycle hook called once by the solve loop after the last step (regular end or callback-requested stop). No-op by default; built-ins use it for a final analysis row / snapshot / timing summary.

initialize!(cb, state) -> nothing

Lifecycle hook called once by the solve loop before the first step, with state.t == t0. The default is a no-op, so plain closures need nothing; built-in callbacks use it to anchor time-based schedules, record the t0 conservation reference, print headers, etc. Extend it for custom callback types (TwoDG.Callbacks.initialize!(cb::MyCallback, state) = ...).