Public API

High-level problem/solve API: thin orchestration over the validated low-level solvers. Physics (equations, boundary conditions, numerical fluxes) lives in TwoDG.Equations and is passed through to the kernels unchanged — there is no translation layer.

using TwoDG

mesh = mkmesh_square(17, 17, 3, 0, 1)
eq   = EulerEquations(γ = 1.4)
prob = DGProblem(eq, mesh; bc = [FarField(uinf), SlipWall(), FarField(uinf), FarField(uinf)],
                 u0 = [ρ0, ρu0, ρv0, ρE0])
sol  = solve(prob, RK4(); dt = 1e-3, tfinal = 1.0)
err  = l2error(sol, exact)      # or scaplot(sol.prob.mesh, sol.u[:, 1, :])

CGProblem(equation, mesh; source, reaction=0.0)

Continuous-Galerkin discretization of a Poisson / convection-diffusion(- reaction) problem with homogeneous Dirichlet boundaries. source is a function (x, y) -> value. Solve with Direct (the default).

ConjugateGradient(; tol=1e-10, maxit=5000, preconditioner=true)

Jacobi-preconditioned conjugate-gradient solve of the (SPD) CG stiffness system, matrix-free on any KA backend (ArrayT=CuArray runs the iteration on the GPU). Only valid for CGProblems without convection; use GMRES otherwise.

DGProblem(equation, mesh; bc, u0, source=nothing, numerical_flux=nothing,
          stabilization=nothing, T=Float64)

Explicit (L)DG semidiscretization of equation on mesh.

• bc — one BoundaryCondition per boundary tag (a vector/tuple in tag order, or a NamedTuple keyed by the mesh's boundary names).
• u0 — initial condition: an (npl, nc, nt) array or a vector of nc constants / (x, y) -> value functions.
• source — nothing or a pointwise source (u, x, t) -> SVector.
• numerical_flux — any callable (eq, uL, uR, n, x, t) -> SVector; defaults to default_numerical_flux(equation).
• stabilization — LDG penalty policy (LDGStabilization) for diffusive equations; defaults to default_stabilization(equation).

Solve with RK4:

solve(prob, RK4(); dt, tfinal (or nstep), t0=0.0, ArrayT=Array)

ArrayT=CuArray (with CUDA.jl loaded) runs the whole time loop on the GPU.

Direct (sparse LU) solve of the condensed system.

GMRES(; restart=80, tol=1e-6, maxit=2000, preconditioner=true, batched=true)

Restarted, block-Jacobi-preconditioned GMRES (Krylov.jl) on the HDG trace system. batched=true uses the batched KA assembly/recovery path (hdg_parsolve_batched); otherwise the per-element threaded assembly (hdg_parsolve).

HDGProblem(equation, mesh; bc, source=nothing, stabilization=1.0)

Steady convection-diffusion (or Poisson) problem discretized with HDG, on triangular (2D) or tetrahedral (3D) meshes. bc is a single Dirichlet applied on the whole boundary (constant, (x, y) -> g, or (x, y, z) -> g); source is nothing or a function of the quadrature-point coordinate matrix p -> values. Solve with Direct (sparse LU) or GMRES (default).

Classic four-stage Runge-Kutta time integration (the GPU-tight internal stepper).

solve(prob::DGProblem, RK4(); dt, tfinal (or nstep), t0=0.0,
      ArrayT=Array, ngauss=nothing, callback=nothing, restart=nothing)

Run the internal RK4 time loop.

callback — any callable cb(state::SolveState) -> Union{Nothing, Bool}: a plain closure, a built-in callback (ProgressCallback, AnalysisCallback, SaveSolutionCallback, SteadyStateCallback, CheckpointCallback, StepsizeCallback), or a CallbackSet composing several. It is called after every step; state.u is the live solution array (device-resident under ArrayT=CuArray; copy or Array(state.u) it if you keep it). Return true to stop early. Callbacks are observers: with a fixed dt the computed solution is bit-identical with and without them. A StepsizeCallback may write state.dt; the loop then advances at the new step size and, when tfinal is given, clamps the last step to land on tfinal exactly. The callback rides back on the solution as sol.callbacks.

restart — path to a CheckpointCallback file; the solve resumes from its u/t/step (the problem's u0 and the t0 keyword are ignored, and nstep counts the additional steps to take).

compute_dt(prob::DGProblem; cfl=0.3) -> dt

CFL-limited explicit time step: over all elements,

dt = cfl * min 1 / ( λ (2p+1)/h + κ ((2p+1)/h)² )

with h the smallest inscribed-circle (2D) / inscribed-sphere (3D) diameter (min_inscribed_diameter), p the polynomial order, λ the maximum characteristic speed of the equation — for state-dependent equations (Euler, user equations) the maximum of wavespeed(eq, u) over the initial state — and κ its diffusivity (LDG diffusion limits the step quadratically). cfl = 0.3 is a conservative default for the internal RK4 stepper. StepsizeCallback applies the same formula to the running solution.

semidiscretize(prob::DGProblem, tspan; ArrayT=Array, ngauss=nothing) -> ODEProblem

Spatial semidiscretization of prob as a SciML ODEProblem, so any OrdinaryDiffEq.jl time integrator (adaptive, SSP, IMEX, …) can drive the DG right-hand side. Requires SciMLBase to be loaded (it usually is, via any using OrdinaryDiffEq... package; otherwise using SciMLBase).

using TwoDG, OrdinaryDiffEqTsit5
ode = semidiscretize(prob, (0.0, 1.0))
sol = solve(ode, Tsit5())

Composable run-time layer for the internal time-stepping loop, following Trixi.jl's callback catalog and Oceananigans.jl's schedule/callback split: schedules decide when to fire, callbacks decide what to do, and a SolveState is the documented view of the running solve they act on.

Any callable cb(state::SolveState) -> Union{Nothing, Bool} is a valid callback — a plain closure, one of the built-ins (ProgressCallback, AnalysisCallback, SteadyStateCallback, SaveSolutionCallback, CheckpointCallback, StepsizeCallback), or a CallbackSet composing several. Returning true stops the solve. Custom callback types may additionally extend initialize! and finish!, which the solve loop calls once before the first step and once after the last.

AbstractSchedule

A schedule decides when a callback fires; the callback decides what happens — the two concerns never mix (Oceananigans' design). A schedule is a small callable object: sched(state::SolveState) -> Bool. Stateful schedules are anchored at t0 by initialize!(sched, state) and advance their own next-fire state when they report true.

Built-in schedules: EveryStep, IterationInterval, TimeInterval, SpecifiedTimes, WallTimeInterval. All built-in callbacks accept schedule = ..., with interval = n as sugar for IterationInterval(n).

AnalysisCallback(; schedule=IterationInterval(100), interval=nothing,
                 integrals=NamedTuple(), errors=nothing, io=stdout)

In-loop solution analysis (Trixi's AnalysisCallback): at t0, at each firing of schedule, and at the end of the run, computes

• the component-wise conserved integrals ∫ u_c dx and their maximum drift against the t0 values (a conservative scheme should keep this at round-off),
• the per-component min/max of the nodal values,
• the user integrals — a NamedTuple of pointwise functionals (eq, u::SVector) -> Real (e.g. (ek = energy_kinetic, s = entropy) or any closure), each integrated with integrate,
• optionally the L² errors against errors — a function or tuple of functions exact(x::SVector{Dim}, t) -> Real, one per solution component starting from the first, integrated with the context's quadrature,

prints them as a table row to io, and appends them to the history carried by the callback: cb.time::Vector{Float64}, cb.steps::Vector{Int}, and cb.data::Dict{Symbol, Vector{Float64}} with keys :conservation_drift, Symbol(:integral_, var), Symbol(:min_, var) / Symbol(:max_, var) per component name, each user-integral name, and Symbol(:l2error_, var) per error function. The history also rides on the returned solution as sol.callbacks.

All reductions run on the device holding state.u; the callback is an observer — the computed solution is bit-identical with and without it.

CallbackSet(callbacks...)

Compose callbacks: calling the set calls every member in order (all run even if an earlier one requests a stop) and returns true — stop the solve — if any member returned true. initialize!/finish! forward to every member. Any mix of built-in callbacks and plain closures is allowed.

CheckpointCallback(; path="checkpoint.jls", schedule=WallTimeInterval(600), interval=nothing)

Restart-file writer: at each firing (every 10 wall-clock minutes by default) serializes the full solver state — u (host copy), t, step, dt — to the single file path, atomically (written to a temporary file first, then renamed, so a killed run never leaves a torn checkpoint). Resume with

solve(prob, RK4(); dt, nstep (or tfinal), restart = path, ...)

which restores u/t/step from the file and continues; a resumed run reproduces the uninterrupted one to floating-point tolerance.

EveryStep() — fire after every time step.

IterationInterval(n)

Fire every n steps (step % n == 0).

NaNCheckCallback(; schedule=EveryStep(), interval=nothing, io=stdout)

Early-abort guard for diverging runs (in the spirit of Oceananigans' NaNChecker): at each firing checks the live solution for non-finite values (NaN/Inf) with a single device reduction and, if any are found, reports the step, the time, and the offending solution component(s) to io, then returns true to stop the solve. The partial solution is still returned, so the incipient blow-up can be inspected (sol.u, plotting, save_vtk) instead of burning the rest of the wall-clock budget on a field of NaNs.

The healthy-path cost is one isfinite reduction over the field per firing; on a GPU each firing synchronizes the device, so pass interval = n to check every n steps if that matters.

ProgressCallback(; schedule=IterationInterval(100), interval=nothing, io=stdout)

Heartbeat for long runs (Trixi's AliveCallback): at each firing prints one line with the step, time, current dt, wall-clock rate since the last firing, and an ETA (from nstep or tfinal, whichever drives the run); at the end of the solve prints the total wall time. Never touches the solution; allocation-free on steps where it does not fire.

SaveSolutionCallback(; schedule=IterationInterval(100), interval=nothing,
                     path="out", prefix="solution", fields=(:u,),
                     save_initial=true, save_final=true)

Snapshot writer: at t0 (if save_initial), at each firing, and at the end of the run (if save_final), copies the requested fields to the host and serializes them to path/prefix_<step>.jls (stdlib Serialization; read back with Serialization.deserialize(file), which returns a NamedTuple with t, step, and one entry per field). Written paths are collected in cb.files.

fields entries are either :u (the full conserved array) or name => f pairs of pointwise functions with the derived_field contract, e.g. fields = (:u, :p => pressure, :M => mach).

SolveState

The view of a running solve that a callback receives after every step (and at t0 through initialize!). Fields:

• u — the live solution array (npl, nc, nt). Under ArrayT = CuArray it is device-resident; copy it (Array(state.u)) if you keep a reference.
• t::Float64, step::Int — current time and step count.
• dt::Float64 — the time step. Writable: a callback that assigns state.dt (see StepsizeCallback) controls the loop from the next step on.
• nstep::Int, tfinal::Float64 — the planned end of the run (typemax(Int) / NaN when the run is driven by the other criterion); used for progress/ETA reporting.
• prob, ctx, phys — the problem, the GeometricFactors geometry cache (for quadrature-exact diagnostics, see integrate), and the compiled physics bundle.

These fields are the documented callback API; anything else is internal.

SpecifiedTimes(ts...)
SpecifiedTimes(ts::AbstractVector)

Fire on the first step whose time reaches each of the given times (sorted internally). Times beyond the end of the run simply never fire; a step that crosses several times fires once.

SteadyStateCallback(; abstol=1e-8, reltol=1e-6,
                    schedule=IterationInterval(10), interval=nothing)

Terminate the solve when the solution stops changing: at each firing compares the finite-difference rate ‖u - u_prev‖₂ / Δt between consecutive firings against abstol + reltol * ‖u‖₂ and returns true (stop) once it falls below. The RK4 stepper does not expose its stage residuals, so the criterion is this between-firings rate — fire it every few steps, not every step, for a meaningful Δt. The comparison state is kept on the same device as u; nothing crosses the bus but scalars.

StepsizeCallback(; cfl=0.3, schedule=EveryStep(), interval=nothing)

Adaptive CFL time step (Trixi's StepsizeCallback): recomputes

dt = cfl / ( λ (2p+1)/h + κ ((2p+1)/h)² )

from the current solution and writes it to state.dt, which the RK4 loop reads before every step. λ is the maximum of the pointwise characteristic-speed bound wavespeed(eq, u) over all solution nodes (a device reduction), κ the equation's diffusivity, and h the smallest inscribed-simplex diameter of the mesh (cached at initialize!) — the same formula as compute_dt, re-evaluated as the solution develops.

The initial dt is set at initialize!, so the dt passed to solve is only a placeholder. Solve with tfinal (the loop then clamps the last step to land on tfinal exactly); with nstep, the run simply takes nstep CFL-sized steps. This is the one callback that controls the loop rather than observing it, and it does not translate to the SciML bridge — OrdinaryDiffEq owns adaptivity there.

TimeInterval(Δt)

Fire on the first step whose time reaches the next multiple t0 + k·Δt (fixed-step loop: no interpolation to the exact time — the firing t may overshoot by up to one dt). A step that crosses several multiples fires once.

WallTimeInterval(seconds)

Fire when at least seconds of wall-clock time have elapsed since the last firing (or since initialize!). For heartbeats and CheckpointCallbacks on long runs.

integrate(u, ctx) -> Vector

Component-wise integrals ∫_Ω u_c dx of a DG field u (npl, nc, nt) — the conserved totals a scheme should preserve. AnalysisCallback records them at t0 and reports the drift. Device-capable; returns a host Vector of length nc.

integrate(f, eq, u, ctx) -> Real

Integral ∫_Ω f(eq, u(x)) dx of a pointwise functional over a DG field u (npl, nc, nt), using the quadrature of the geometry cache ctx::GeometricFactors. f is any callable (eq, u::SVector) -> Real — the same contract as derived_field — e.g. pressure, energy_kinetic, entropy, or a user closure. Runs as a KernelAbstractions kernel plus a device reduction, so u may live on the GPU; only the scalar result crosses to the host.

l2norm(u, ctx; component=nothing) -> Real

‖u‖_{L²(Ω)} of one solution component (or, with component = nothing, the root-sum-square over all components) of a DG field u (npl, nc, nt), integrated with the quadrature of ctx::GeometricFactors. Device-capable.

Mesh

Solver-ready triangular mesh: vertex geometry, DG connectivity, and high-order (possibly curved) nodes. Built by the mkmesh_* generators or by discretizeing a MeshGeometry; every mesh those return has p–f2f fully populated (pcg/tcg are filled by cgmesh, which discretize runs for you). The keyword constructor's nothing defaults are construction-internal staging only — solver code must never receive a mesh with nothing fields.

Fields and conventions

• p :: (np, 2) — vertex coordinates.
• t :: (nt, 3) — triangles as vertex indices, counterclockwise.
• f :: (nf, 4) — faces: f[i, 1:2] are the endpoint vertices, f[i, 3] the element to the left, f[i, 4] the element to the right or -k when face i lies on boundary k (see boundary_names). Interior faces are listed first, then boundary faces grouped by tag.
• t2f :: (nt, 3) — element-to-face map (unsigned): local face j of element it is t2f[it, j].
• t2o :: (nt, 3) — orientation code of each local face w.r.t. the face's stored traversal (1 matching, 2 reversed in 2D; up to 6 codes on tetrahedra), indexing master.perm[:, s, o]. Replaces the former sign of t2f; meshes constructed with a signed t2f and no t2o are migrated automatically for one release.
• fcurved :: (nf,), tcurved :: (nt,) — flags marking faces/elements that touch a curved boundary (their high-order nodes are projected during createnodes).
• porder — polynomial order of the elements.
• plocal :: (npl, 3), tlocal — master-element node positions (barycentric) and their triangulation, npl = (porder+1)(porder+2)/2.
• dgnodes :: (npl, 2, nt) — coordinates of the high-order DG nodes of each element (isoparametric on curved elements).
• pcg, tcg — deduplicated continuous (CG) node coordinates and connectivity, filled by cgmesh.
• elcon :: (porder+1, 3, nt) — element-to-global trace-node connectivity used by HDG: global numbering of the face nodes of each local face, already orientation-corrected via the sign of t2f.
• f2f :: (nf, 5) — face-to-face connectivity (mkf2f), used by the block-Jacobi preconditioner.
• boundary_names :: Vector{Symbol} — names of the boundary tags, in tag order (empty when the generator did not attach names).

MeshGeometry(p, t; boundaries, curved=Symbol[], fd=nothing)

Geometry-stage description of a triangulated domain, before any choice of polynomial order or node placement.

• p (np, 2): vertex coordinates; t (nt, 3): triangle vertex indices.
• boundaries: a NamedTuple (or pairs) of name => classifier, where classifier(midpoints) -> Bool selects the boundary faces belonging to name. The order defines the boundary tags (-1, -2, … in mesh.f).
• curved: names of boundaries that are curved (their high-order face nodes are projected onto the true boundary during discretize).
• fd: nothing, or one signed-distance function per boundary (same order as boundaries); required when curved is nonempty.

See discretize.

boundary_names(mesh) -> Vector{Symbol}

Names of the mesh's boundary segments, in boundary-tag order (tag k is stored as -k in mesh.f[:, 4]). Empty when the generator did not attach names.

box_geometry(m=2, n=2, o=2; lengths=(1.0, 1.0, 1.0)) -> MeshGeometry{3}

Geometry of the box [0,L₁]×[0,L₂]×[0,L₃] on an m × n × o vertex grid (Kuhn 6-tet cells), with boundaries named :left/:right (x), :front/ :back (y), :bottom/:top (z), tags 1–6.

cgmesh(mesh, ptol=2e-13) -> Mesh

Build the continuous-Galerkin numbering of mesh by deduplicating the DG nodes shared between elements: fills pcg (unique node coordinates) and tcg (nt, npl) (element-to-CG-node connectivity). Required before cg_solve. Nodes closer than ptol (after rounding) are merged.

circle_geometry(p, t) -> MeshGeometry

Geometry of the unit circle from a raw triangulation (p, t) (deduplicated with fixmesh). The single boundary is named :boundary and is curved: high-order nodes are projected onto the unit circle during discretize.

createnodes(mesh, fd=nothing) -> Mesh

Compute the high-order DG node coordinates dgnodes (npl, 2, nt) of mesh by mapping the master-element nodes (mesh.plocal) into each triangle. When fd is a vector of signed-distance functions (one per boundary tag), nodes on curved boundary edges are projected onto the true boundary, making those elements isoparametric.

Also fills the HDG trace connectivity elcon and the face-to-face map f2f when the mesh already has face connectivity (mkt2f). Returns a new Mesh with the extra fields populated.

discretize(geo::MeshGeometry, porder; nodetype=0) -> Mesh

Discretization stage: build the solver-ready Mesh from a MeshGeometry — face/element connectivity (f, t2f, t2o, elcon, f2f), high-order nodes (dgnodes, projected onto curved boundaries in 2D), and CG numbering (pcg, tcg). nodetype = 0 places nodes uniformly, 1/2 use extended Chebyshev points (2D only).

For a MeshGeometry{3} the mesh is tetrahedral: face orientations carry the 6-code t2o and the HDG trace connectivity goes through the triangle face element. Curved boundaries project the boundary-face nodes onto the signed- distance zero set; nodes on edges shared by two different curved boundaries are projected onto the intersection curve by alternating projections and written back to every element that shares them (the edge-before-face rule).

face_orientation(stored, mine, ::Val{Dim}) -> Int

Orientation code o of an element's face traversal mine (tuple of global vertex ids, the element's outward ordering, Dim of them) relative to the face's canonical stored traversal stored: the index of the face-simplex symmetry σ with σ[j] = position of mine[j] in stored, so that master.perm[:, s, o] lists the element's volume nodes in the canonical face-node order.

face_vertices(::Val{Dim}, j) -> NTuple{Dim, Int}

Local vertices of local face j (the face opposite vertex j, where barycentric coordinate j vanishes), in the traversal that makes the face outward-oriented for the reference element: counterclockwise edges in 2D, right-hand-rule outward triangles in 3D.

fixmesh(p::Matrix{T}, t::Matrix{Int}, ptol::Real=2e-13) where T<:Real

Remove duplicated nodes and fix element orientation in a simplicial mesh (triangles or tetrahedra).

Parameters:

• p: N×Dim matrix of vertex coordinates
• t: M×(Dim+1) matrix of simplex vertex indices
• ptol: tolerance for identifying duplicate vertices (default: 2e-13)

Returns:

• Tuple of (cleaned vertex matrix, fixed simplex matrix with positive orientation)

gmsh_geometry(filepath; boundaries, curved=Symbol[], fd=nothing) -> MeshGeometry

Read a linear Gmsh .msh file into a MeshGeometry: tetrahedra give a MeshGeometry{3}, triangles a MeshGeometry{2}. Duplicate nodes are merged and elements positively oriented (fixmesh); boundaries, curved, fd as in MeshGeometry (curved boundaries are projected onto their signed-distance zero sets by discretize). Implemented in the TwoDGGmshExt package extension: available once Gmsh.jl is loaded (using TwoDG, Gmsh).

lshape_geometry(m=2; parity=0) -> MeshGeometry

Geometry of the unit L-shape assembled from three subsquares of m × m vertices each. The whole boundary carries the single name :boundary.

make_box_mesh(m=2, n=2, o=2; lengths=(1.0, 1.0, 1.0)) -> (p, t)

Structured tetrahedral mesh of the box [0,L₁]×[0,L₂]×[0,L₃] on an m × n × o vertex grid: each grid cell is split into 6 tetrahedra around its main diagonal (Kuhn/Freudenthal triangulation) — unlike the 5-tet split it is conforming without parity alternation and refines uniformly. Returns vertices p (np, 3) and positively oriented tetrahedra t (nt, 4).

make_circle_mesh(size, boundary_refinement) -> (p, t)

Raw unstructured triangulation of the unit circle with element size size, generated by shelling out to Python's distmesh (pyscripts/; requires python with numpy on the path, distmesh is vendored). Pass boundary_refinement = h_boundary for a mesh graded toward the boundary, or nothing for uniform sizing. Returns vertices p (np, 2) and triangles t (nt, 3); use mkmesh_circle for a solver-ready mesh.

make_circle_nodes(p, t, porder, nodetype) -> Mesh

Discretize a raw unit-circle triangulation (p, t) (e.g. from make_circle_mesh) at order porder, projecting the high-order boundary nodes onto the unit circle. Equivalent to discretize(circle_geometry(p, t), porder; nodetype).

makesquaremesh 2-d regular triangle mesh generator for the unit square p, t = squaremesh(m, n, parity)

p: node positions (np,2) t: triangle indices (nt,3) parity: flag determining the triangular pattern flag = 0 (diagonals sw - ne) (default) flag = 1 (diagonals nw - se)

mkf2f(f, t2f)

Create face to face connectivity.

Arguments

• f: Face to element connectivity
• t2f: Element to face connectivity

Returns

• f2f::Matrix{Int}: Face to face connectivity

mkmesh_box(m=2, n=2, o=2, porder=1; lengths=(1.0, 1.0, 1.0)) -> Mesh

Solver-ready tetrahedral mesh of the box: discretize(box_geometry(m, n, o; lengths), porder).

mkmesh_circle(siz=0.4, porder=3, nodetype=0; boundary_refinement=nothing) -> Mesh

Unstructured curved mesh of the unit circle (element size siz), generated with distmesh via Python (see make_circle_mesh) and discretized at order porder.

mkmesh_distort!(mesh, wig=0.05)

Distort a unit-square mesh (from mkmesh_square) in place with a smooth sinusoidal warp of amplitude wig, keeping the boundary fixed. Useful for testing solver accuracy on non-affine element mappings.

tfiduct map a unit square mesh to a cos^2 duct mesh = mkmeshduct(mesh,db,dt,H) mesh: mesh structure generated by mkamesh_square flag = 0 (diagonals sw - ne) (default) flag = 1 (diagonals nw - se) db: height of bottom bump dt: height of top bump H: height of channel

mkmesh_lshape(m=2, porder=1, parity=0, nodetype=1) -> Mesh

Solver-ready mesh of the unit L-shape: discretize(lshape_geometry(m; parity), porder; nodetype).

mkmesh_naca(t_naca=10, porder=2, name="naca0012", display_gmsh=false)

Generate a curved high-order mesh around a NACA 4-digit symmetric airfoil (thickness t_naca in percent of chord). Implemented in the TwoDGGmshExt package extension: it becomes available once Gmsh.jl is loaded (using TwoDG, Gmsh).

mkmesh_square(m=2, n=2, porder=1, parity=0, nodetype=0) -> Mesh

Solver-ready mesh of the unit square: discretize(square_geometry(m, n; parity), porder; nodetype).

mkmesh_trefftz(m=15, n=30, porder=3, node_spacing_type=0,
               tparam=[0.1, 0.05, 1.98]) -> Mesh

Curved O-mesh around a Trefftz (Karman–Trefftz) airfoil, built by conformally mapping an m × n structured rectangle through exp and the K–T transform. tparam = [x0, y0, exponent] are the airfoil parameters: circle-center shifts x0/y0 and K–T exponent (2 gives a Joukowski airfoil). Boundary tag 1 (:airfoil) is the airfoil surface, tag 2 (:farfield) the outer circle. See also trefftz for the potential-flow driver.

mkt2f(t) -> (f, t2f, t2o)

Face connectivity of a simplicial mesh: triangles t (nt, 3) or tetrahedra t (nt, 4) (routed by the column count):

• f (nf, Dim+2): face rows [v..., left element, right element] — Dim vertices, then the two adjacent elements, the right entry 0 for boundary faces (later replaced by -tag, see setbndnbrs). The stored vertex order is the left element's outward traversal (counterclockwise edge in 2D; right-hand-rule outward triangle in 3D).
• t2f (nt, Dim+1): face index of each local face (unsigned; local face j is opposite local vertex j).
• t2o (nt, Dim+1): orientation code of each local face — the index into master.perm[:, s, o] mapping the face's canonical (stored) node ordering to the element's own traversal: 1 when they match (always the case for the left element), 2 when reversed in 2D; 1:6 over the triangle symmetry group in 3D. This replaces the former sign of t2f (an explicit small integer is the representation that scales to 3D).

norient(Dim) -> Int

Number of relative orientations in which two elements can see a shared face: the order of the symmetry group of the Dim-1 face simplex — 1 in 1D (point face), 2 in 2D (segment: identity + reversal), 6 in 3D (triangle: 3 rotations × optional reflection).

orientation_permutations(::Val{Dim}) -> NTuple{norient, NTuple{Dim, Int}}

The symmetry group of the Dim-element's face simplex, as permutations of the face's vertex positions; orientation code o refers to the o-th entry. 2D face = segment: identity, reversal. 3D face = triangle: the 3 rotations, then the 3 reflected traversals.

p: Node positions (:,2) f: Face Array (:,4) bndexpr: Cell Array of boundary expressions. The number of elements in BNDEXPR determines the number of different boundaries

Example: (Setting boundary types for a unit square mesh - 4 types) bndexpr = [lambda p: np.all(p[:,0]<1e-3, lambda p: np.all(p[:,0]>1-1e-3), lambda p: np.all(p[:,1]<1e-3, lambda p: np.all(p[:,1]>1-1e-3)] f = setbndnbrs(p,f,bndexpr);

Example: (Setting boundary types for the unit circle - 1 type) bndexpr = [lambda p: np.all(np.sqrt((p**2).sum(1))>1.0-1e-3)] f = setbndnbrs(p,f,bndexpr);

square_geometry(m=2, n=2; parity=0) -> MeshGeometry

Geometry of the unit square on an m × n vertex grid, with boundaries named :bottom, :right, :top, :left (tags 1–4). parity selects the diagonal direction (see make_square_mesh).

uniref(p, t, nref=1)

Uniformly refine a simplicial mesh: each triangle into 4 (edge midpoints), each tetrahedron into 8 (Bey's red refinement — 4 corner tets plus the interior octahedron split around its shortest diagonal for shape regularity; Bey, Computing 55, 1995).

Arguments

• p::Matrix{T}: Nodes as a matrix of size (np, dim)
• t::Matrix{Int}: Triangulation as a matrix of size (nt, dim+1)
• nref::Int=1: Number of uniform refinements

Returns

• p::Matrix{T}: Refined nodes
• t::Matrix{Int}: Refined triangulation (positive orientation preserved)

Master

Deprecated alias for ReferenceElement; will be removed one release after the rename (see NEWS.md).

ReferenceElement(porder; pgauss=max(4porder, 1), nodetype=0)
ReferenceElement(plocal, porder; pgauss=max(4porder, 1))
ReferenceElement(mesh::Mesh, pgauss=nothing)

Reference (master) simplex of polynomial order porder in Dim dimensions: tabulated shape functions, quadrature rules, and the local node orderings the solvers need. The face of the element is itself the Dim-1-dimensional reference simplex, stored whole in the face field (the triangle's face is the 1D segment; the tetrahedron's face is the full triangle element), so face quadrature, face shape functions, and the HDG trace basis come from one recursive structure. The element is mesh-independent; the Mesh convenience constructor just reads mesh.porder/mesh.plocal so the element's nodes match the mesh's dgnodes.

pgauss is the polynomial degree integrated exactly by the quadrature rules; nodetype selects the node distribution (0 uniform; 1 extended Chebyshev in 2D, warp-and-blend in 3D; 2 extended Chebyshev of the second kind, 2D only — see localpnts and localpnts3d); plocal may be passed directly to use a custom node set.

Fields and conventions (nv = Dim + 1 vertices/faces, npf nodes per face,

norient face orientations: 1 in 1D, 2 in 2D, 6 in 3D)

• porder — polynomial order; npl nodes ((porder+1)(porder+2)/2 for the triangle, (porder+1)(porder+2)(porder+3)/6 for the tetrahedron).
• plocal :: (npl, nv) — node positions in barycentric coordinates.
• corner :: Vector{Int} — indices of the nv vertex nodes in plocal.
• perm :: (npf, nv, norient) — face-node orderings: perm[:, j, o] lists the volume nodes on local face j matching the face's canonical traversal under orientation o (o = 1 canonical, o = 2 reversed in 2D; the 6 triangle symmetries in 3D). Used with mesh.t2o when a neighboring element sees the shared face in a different orientation.
• gpts :: (ng, Dim), gwgh :: (ng,) — volume quadrature points/weights.
• shap :: (npl, Dim+1, ng) — shape functions ([:, 1, :]) and their reference-coordinate derivatives ([:, 1+d, :] for direction d) at the volume quadrature points.
• mass :: (npl, npl) — reference-element mass matrix.
• conv :: (npl, npl, Dim) — reference convection matrices ($∫ φᵢ ∂φⱼ/∂ξ_d$).
• face — the Dim-1-dimensional ReferenceElement of the faces (nothing for the 1D segment), built with the matching node distribution and quadrature degree.

Deprecated property aliases (2D)

The pre-Dim 1D face tables remain readable as properties for one release: master.sh1d == master.face.shap, master.ma1d == master.face.mass, master.gw1d == master.face.gwgh, master.gp1d == vec(master.face.gpts), master.ploc1d == master.face.plocal.

ReferenceElement(porder; dim=2, nodetype=0, pgauss=max(4porder, 1))

Reference simplex of order porder: the triangle for dim = 2, the tetrahedron for dim = 3 (whose face field is the full triangle element).

gaussquad1d(pgauss::Integer)

Calculate Gauss integration points and weights for the interval [0,1] using Gauss-Legendre quadrature (special case of Gauss-Jacobi with α = β = 0).

Arguments

• pgauss::Integer: Order of the polynomial to be integrated exactly

Returns

• x::Vector{Float64}: Coordinates of the integration points
• w::Vector{Float64}: Integration weights

Example

x, w = gaussquad1d(3)

gaussquad2d(pgauss::Integer)

Calculate Gauss integration points and weights in 2D for [0,1]×[0,1].

Arguments

• pgauss::Integer: Order of the polynomial to be integrated exactly

Returns

• x::Matrix{Float64}: Coordinates of the integration points (N×2 matrix)
• w::Vector{Float64}: Integration weights

gaussquad3d(pgauss::Integer)

Gauss integration points and weights on the unit tetrahedron {ξ ≥ 0, ξ₁+ξ₂+ξ₃ ≤ 1}, exact for total polynomial degree pgauss.

For pgauss ≤ 10 this returns the tabulated symmetric rules of Witherden & Vincent (Comput. Math. Appl. 69, 2015) — interior points, positive weights, and 2–3× fewer points than a product rule of the same degree (e.g. 46 vs 125 at degree 8). Above degree 10 it falls back to collapsed-coordinate Gauss–Jacobi product rules (degree-exact for any pgauss; Hesthaven & Warburton 2008, §10.2): the cube [-1,1]³ maps to the tetrahedron with Jacobian (1-η₂)(1-η₃)²/64, absorbed into Gauss–Jacobi weights with α = 1 and α = 2 — no quadrature point lies on the collapsed edge.

Returns

• x::Matrix{Float64}: integration points (ng, 3)
• w::Vector{Float64}: integration weights (summing to the volume 1/6)

get_local_face_nodes(mesh, master, face_number, flip_face_direction=false)

Local (element) indices of the porder+1 nodes lying on global face face_number of the element to its left (mesh.f[face_number, 3]), in counterclockwise face order — or reversed with flip_face_direction=true, which matches how the right element sees the same face (master.perm[..., 2]).

koornwinder1d vandermonde matrix for legenedre polynomials in [0,1] [f,fx]=koornwinder(x,p)

x: coordinates of the points wherethe polynomials are to be evaluated (npoints) p: maximum order of the polynomials consider. that is all polynomials of degree up to p, npoly=p+1 f: vandermonde matrix (npoints,npoly) fx: vandermonde matrix for the derivative of the koornwinder polynomials w.r.t. x (npoints,npoly)

koornwinder2d(x, p) -> (f, fx, fy)

Vandermonde matrices of the orthonormal PKD (Proriol–Koornwinder–Dubiner) basis on the master triangle [0,0]-[1,0]-[0,1] and of its two derivatives, at the points x (npoints, 2). All polynomials of complete degree up to p, npoly = (p+1)(p+2)/2, orthonormal w.r.t. the unit-triangle measure (collapsed-coordinate + Jacobi-polynomial construction; Hesthaven & Warburton 2008, §6.1). koornwinder3d is the tetrahedral analog.

koornwinder3d(x, p) -> (f, fx, fy, fz)

Vandermonde matrices of the orthonormal PKD (Proriol–Koornwinder–Dubiner) basis on the unit tetrahedron [0,0,0]-[1,0,0]-[0,1,0]-[0,0,1] and of its three derivatives, at the points x (npoints, 3). All polynomials of complete degree up to p, npoly = (p+1)(p+2)(p+3)/6, orthonormal w.r.t. the unit-tetrahedron measure. Direct extension of koornwinder2d's collapsed-coordinate + Jacobi-polynomial recipe (Hesthaven & Warburton 2008, §10.1):

φ_pqr = c · P_p(a) · ((1-b)/2)^p P_q^{(2p+1,0)}(b) · ((1-c)/2)^{p+q} P_r^{(2p+2q+2,0)}(c)

with (a, b, c) the collapsed coordinates of the bi-unit tetrahedron and c = sqrt(2(2p+1)(p+q+1)(2(p+q+r)+3)).

localpnts 2-d mesh generator for the master element. Returns (plocal, tlocal) where: plocal: node positions (npl,3) in barycentric coordinates tlocal: triangle indices (nt,3) porder: order of the complete polynomial npl = (porder+1)*(porder+2)/2

Compute node positions on the master volume element.

Arguments

• porder::Integer: Polynomial order
• nodetype::Integer=0: Flag determining node distribution:
  • nodetype = 0: Uniform distribution (default)
  • nodetype = 1: Extended Chebyshev nodes of the first kind
  • nodetype = 2: Extended Chebyshev nodes of the second kind

Returns

• plocal::Vector{Float64}: Vector of node positions on the master volume element

localpnts3d(porder, nodetype=0) -> plocal (npl, 4)

Node positions on the master tetrahedron in barycentric coordinates, npl = (porder+1)(porder+2)(porder+3)/6.

• nodetype = 0: uniform barycentric lattice (adequate conditioning for porder ≤ 4).
• nodetype = 1: warp-and-blend nodes (Warburton, J. Eng. Math. 56, 2006; Hesthaven & Warburton 2008 §10.5) — much better Vandermonde conditioning at high order.

Both distributions are invariant under the tetrahedron's symmetry group and restrict to a symmetric triangle node set on every face — the two properties perm/t2o depend on (asserted in the constructor). The warp-and-blend set is symmetrized exactly by averaging over the 24 barycentric-coordinate permutations (the 3D analog of localpnts' rotate-and-average).

shape1d calculates the nodal shapefunctions and its derivatives for the master 1d element [0,1]

Arguments:

• porder: polynomial order
• plocal: node positions (np,2) (np=porder+1)
• pts: coordinates of the points where the shape functions and derivatives are to be evaluated (npoints)

Returns:

• nsf: shape function and derivatives (np,2,npoints) nsf[:,1,:] shape functions nsf[:,2,:] shape functions derivatives w.r.t. x

shape2d calculates the nodal shapefunctions and its derivatives for the master triangle [0,0]-[1,0]-[0,1] nfs=shape2d(porder,plocal,pts)

porder: polynomial order plocal: node positions (np,2) (np=(porder+1)*(porder+2)/2) pts: coordinates of the points where the shape fucntions and derivatives are to be evaluated (npoints,2) nfs: shape function adn derivatives (np,3,npoints) nsf[:,0,:] shape functions nsf[:,1,:] shape fucntions derivatives w.r.t. x nsf[:,2,:] shape fucntions derivatives w.r.t. y

shape3d(porder, plocal, pts)

Nodal shape functions and derivatives on the master tetrahedron [0,0,0]-[1,0,0]-[0,1,0]-[0,0,1]: nfs (npl, 4, npoints) with values in [:, 1, :] and the ξ/η/ζ derivatives in [:, 2:4, :]. Same Vandermonde-solve pattern as shape2d, on the koornwinder3d basis.

uniformlocalpnts 2-d mesh generator for the master element. [plocal,tlocal]=uniformlocalpnts(porder)

plocal: node positions (npl,3) tlocal: triangle indices (nt,3) porder: order of the complete polynomial npl = (porder+1)*(porder+2)/2

DGContext

Alias for GeometricFactors — the one-time geometry/connectivity cache the KernelAbstractions residual path consumes. Construct with DGContext(master, mesh; T=Float64) and move to a GPU with Adapt.adapt(CuArray, ctx).

DGPhysics(equation; boundary_conditions,
          numerical_flux=default_numerical_flux(equation),
          source=nothing,
          stabilization=default_stabilization(equation))

Everything the DG residual kernels need to know about the physics:

• equation :: AbstractEquation — implements flux/nvariables/….
• boundary_conditions — one BoundaryCondition per boundary tag (any collection; stored as a Tuple so kernels dispatch statically per boundary).
• numerical_flux — any callable (eq, uL, uR, n, x, t) -> SVector.
• source — nothing or a callable (u, x, t) -> SVector.
• stabilization — the LDG penalty policy for diffusive equations.

All components must be isbits (GPU-movable); adapt distributes over the fields.

RinvWorkspace(ctx, nc)

Staging buffers reused across residual evaluations (e.g. RK stages): fng (ngf, nc, nf) weighted face fluxes, fdg (ng, nc, Dim, nt) volume fluxes, srcg (ng, nc, nt) source, rtmp (npl, nc, nt) pre-mass-solve residual. Allocated on the backend of ctx.

RldgWorkspace(ctx, nc)

Staging buffers for the LDG residual path: everything in RinvWorkspace plus ug (ng, nc, nt) (u at volume quadrature points), qfd (ngf, nc, Dim, nf) (weighted gradient face fluxes), and qtmp/q (npl, Dim, nc, nt) (pre-mass-solve and final LDG gradients). Allocated on the backend of ctx.

LDG gradient computation — primary name for getq!.

getq!(q, ctx, phys, u, time; ws=RldgWorkspace(ctx, nvariables(phys)))

LDG gradient q = ∇u ((npl, Dim, nc, nt), already multiplied by the inverse mass matrix). Also leaves u interpolated to volume quadrature points in ws.ug (reused by rldgexpl!).

getq_ka(ctx, phys, u, time)

Allocating convenience wrapper around getq!.

Inviscid DG residual — primary name for rinvexpl!.

rinvexpl!(r, ctx, phys::DGPhysics, u, time; ws=RinvWorkspace(ctx, nvariables(phys)))

Inviscid DG residual r = du/dt (already multiplied by the inverse mass matrix). r, u are (npl, nc, nt) arrays on the same backend as ctx. Pass a pre-built ws to avoid re-allocating staging buffers across calls.

rinvexpl_ka(ctx, phys, u, time)

Allocating convenience wrapper around rinvexpl!.

rk4_ka!([residual!,] ctx, phys, u, time, dt, nstep; ws, stages)

In-place RK4 time integrator driving a KA residual (rinvexpl! by default; pass rldgexpl! for the LDG viscous path). All stage buffers are allocated once up front; the stage updates are plain broadcasts, so the stepper runs unchanged on CPU and GPU arrays. Returns u after nstep steps.

stages — the five stage buffers (each similar(u)); pass a preallocated tuple to make repeated single-step calls (the callback-driven solve loop) allocation-free.

rldgexpl!(r, ctx, phys, u, time; ws=RldgWorkspace(ctx, nvariables(phys)))

LDG (viscous) DG residual r = du/dt (already multiplied by the inverse mass matrix) of the LDG discretization with viscous terms. Falls back to the inviscid rinvexpl! kernels when the equation has no diffusion.

rldgexpl_ka(ctx, phys, u, time)

Allocating convenience wrapper around rldgexpl!.

LDG viscous DG residual — primary name for rldgexpl!.

HDGBatch(master, mesh, source, param; T=Float64)

Per-element operator matrices for the batched HDG assembly/recovery path, built once on the CPU. All fields are plain arrays (Adapt.@adapt_structure), so adapt(CuArray, batch) moves everything to the GPU. Dimension-generic: Dim is the direction axis of C/Le/R (2 on triangles, 3 on tetrahedra).

Shapes (npl volume nodes, nps trace nodes per face — porder + 1 in 2D, (porder+1)(porder+2)/2 in 3D — ndf = (Dim+1) nps trace DOFs per element, nc1 = ndf + 1 local-solve columns — the unit-trace columns plus the source column):

• MinvM (npl, npl, nt): κ × inverse mass (the inverse of localprob's M).
• C (npl, npl, Dim, nt), D (npl, npl, nt): coupling matrices and convection + stabilization operator (stabilization with tau = κ + |c·n|, as in localprob).
• Le (npl, nc1, Dim, nt): face lifts of the trace (F_d = -Le_d m), last column zero so the source column threads through the same GEMMs.
• B0 (npl, nc1, nt): local-solve RHS seed [-Lu | F_src].
• Alam (ndf, ndf, nt), R (ndf, npl, Dim, nt), Ru (ndf, npl, nt): trace-test matrices of the numerical flux (tau = param[:taud] + |c·n|, as in elemmat_hdg).
• elcon (ndf, nt): element trace DOFs -> global trace DOFs (orientation resolved through the mesh's mkelcon connectivity — reversal in 2D, the 6 triangle symmetries in 3D), for the assembly scatter and recovery gather.

HDGCDBatch(master, mesh, κ, τ; T=Float64)

Per-element constants of the batched scalar HDG transport step (hdg_cd_step): like HDGNSBatch but scalar-valued — A0 (npl, npl, nt) and B0 (npl, nfc + 1, nt) carry the κ-viscous composites and the constant −τ⟨θ̂, ·⟩ trace stabilization; the state-dependent convection (volume u, trace λ) is added per call. Dimension-generic like the NS batch.

HDGCDCache(master, mesh, κ, τ, tbc; ArrayT=Array, T=Float64)

Preallocated workspace for hdg_cd_step_batched: the batched element data (batch::HDGCDBatch, on the backend of ArrayT), the trace sparsity pattern, the device work arrays of the local solves/recovery, the host-side quadrature geometry used to evaluate function sources, and the reused sparse-LU factorization of the trace system (F).

Built automatically on the first hdg_cd_step_batched call; pass the returned cache back in on subsequent calls with the same master, mesh, κ, τ, and boundary-condition types to skip all setup and reuse the factorization (this reuse is what makes implicit time stepping cheap).

HDGNSBatch(master, mesh, ν, τ; T=Float64)

Per-element constant operator matrices of the batched HDG Navier-Stokes Newton step, built once on the CPU (all fields plain arrays, Adapt.@adapt_structure). Dimension-generic: with npl volume nodes, nps face nodes per face, nfc = (Dim+1) nps, the local solve size nv = (Dim+1) npl (u₁, …, u_Dim, p) and ncB = Dim·nfc + 2 right-hand-side columns (trace columns + bρ + r):

• shap (npl, ng), shf (nps, nqf): shared shape values at quadrature.
• shapd (npl, ng, Dim, nt): weighted physical derivative tables.
• M (npl, npl, nt): element mass (the dtinv·M term is added per call).
• A0 (nv, nv, nt): constant part of the local Newton matrix, including the viscous composite and the mean-pressure gauge row.
• Bhat0 (nv, ncB, nt): constant part of the local RHS block [B | bρ | r] (r column zero).
• Hx (Dim·nfc, nv, nt), Hlam0 (Dim·nfc, Dim·nfc, nt): flux-continuity test blocks (fully constant, resp. constant part).
• MiE (npl, nfc, Dim, nt), MiC (npl, npl, Dim, nt): gradient-recovery maps M⁻¹E_d, M⁻¹C_dᵀ.
• wds (nqf, Dim+1, nt), fn (nqf, Dim, Dim+1, nt): face quadrature measure and unit normal per local face.
• perm (nps, Dim+1), elcon (nps, Dim+1, nt) (Int32): face-to-volume node map and global face-node connectivity.
• crow (Dim·nfc, nt), area (nt): compatibility row and element measure (consumed on the host during global assembly).

HDGNSCache(master, mesh, ν, τ; ArrayT=Array, T=Float64)

Preallocated workspace for hdg_ns_step_batched (and through it hdg_ns_solve): the batched element data (batch::HDGNSBatch, on the backend of ArrayT), the trace saddle-point sparsity pattern, the device work arrays of the Newton assembly/local solves/recovery, the host-side quadrature geometry used to evaluate function sources, and the reused UMFPACK factorization of the trace system (F).

Built automatically on the first hdg_ns_step_batched call; pass the returned cache back in on subsequent calls with the same master, mesh, ν, τ to skip all setup and refactorize in place — across Newton iterations and time steps this dominates the cost of an unsteady NS run.

HDGSystem(ae, fe, mesh; T=Float64)

The statically condensed HDG trace system in face-block format, ready for the KA/Krylov iterative solver. All fields are plain arrays, so the whole struct moves to a GPU with Adapt.adapt(CuArray, sys).

Fields (ncf trace DOFs per face, nbf = 2nfe - 1 neighbor blocks, nf faces):

• A (ncf, ncf, nbf, nf): global matrix; block k of face i couples face i to face f2f[i, k] (block 1 is the self/diagonal block).
• B (ncf, ncf, nf): inverted diagonal blocks (block-Jacobi preconditioner).
• b (ncf * nf): right-hand side.
• f2f (nf, nbf): face-to-face connectivity, 0 = no neighbor.

elemmat_hdg(dg, master, source, param)

Calculates the element and force vectors for the HDG method.

Arguments

• dg: DG nodes
• master: Master element structure
• source: Source term function or nothing
• param: Dictionary with parameters :kappa (diffusivity) and :c (convective velocity)

Returns

• ae: Element matrix
• fe: Element force vector

hdg_cd_step(master, mesh, κ, tbc; τ=1.0, u=nothing, Λ=nothing, θold=nothing,
            dtinv=0.0, source=nothing)

Solves one (linear) implicit step of the scalar HDG advection-diffusion equation ∂θ/∂t + ∇·(uθ) - κΔθ = s with the DG velocity field u (npl × Dim × nt) and the velocity trace Λ from the HDG Navier-Stokes solver.

tbc(p, tag) prescribes the boundary condition on a boundary face with tag tag (-mesh.f[i, end]) at node coordinates p: return (:d, value) for a Dirichlet condition on θ̂ or (:n, flux) for a prescribed total normal flux (e.g. (:n, 0.0) for an insulated wall).

Returns a named tuple (θ, q, Θ) with the scalar field (npl × nt), its gradient q = ∇θ (npl × Dim × nt), and the face trace.

hdg_cd_step_batched(master, mesh, κ, tbc; τ=1.0, u=nothing, Λ=nothing,
                    θold=nothing, dtinv=0.0, source=nothing,
                    ArrayT=Array, T=Float64, cache=nothing)

Batched/KA counterpart of hdg_cd_step (same arguments, same returned fields plus cache). The per-element assembly, local solves and (θ, q) recovery run on the backend of ArrayT; the trace system is a CPU sparse LU with the pattern and factorization reused across calls. The boundary-condition types returned by tbc must not change between calls with the same cache.

hdg_direct_batched(master, mesh, source, dbc, param; T=Float64)

Direct (sparse LU/Cholesky-free) counterpart of hdg_parsolve_batched: the same batched element assembly (hdg_local_solves) and local recovery, with the condensed trace system assembled to a SparseMatrixCSC and factorized instead of solved iteratively — the Direct/GMRES choice is an algorithm choice over one assembly engine, not a separate code path.

Returns (uh, qh, uhath) with the same shapes as hdg_parsolve_batched.

hdg_gmres_ka(sys::HDGSystem; restart=80, tol=1e-6, maxit=2000,
             preconditioner=true, verbose=false)

Solves the HDG trace system with restarted GMRES (Krylov.jl), left- preconditioned by the block-Jacobi blocks of sys, on whatever backend sys lives on. Returns (x, stats) with stats::Krylov.SimpleStats.

hdg_local_solves(batch::HDGBatch)

Runs the batched local solves and static condensation on the backend of batch. Returns (; ae, fe, U, Q): the element matrices/vectors ae (ndf, ndf, nt), fe (ndf, nt) (Dirichlet BCs not yet applied) and the local solution maps U (npl, ndf + 1, nt) / Q (npl, ndf + 1, Dim, nt) (unit-trace columns plus the source column) consumed by hdg_recover.

hdg_ns_postprocess(master, mesh, master1, mesh1, result)

Element-by-element postprocessing of the HDG Navier-Stokes solution (paper §3.3): finds u* ∈ [P_{k+1}(K)]² such that, on every edge F of K,

⟨(u* - û)·n, μ⟩_F = 0                          ∀μ ∈ P_k(F),
⟨t·∇(u*·n) - nᵀ{{L}}t, t·∇μ⟩_F = 0             ∀μ ∈ P_{k+1}(F)^⊥,

and, in the element interior,

(u* - u, ∇w)_K = 0                             ∀w ∈ P_k(K),
(∇×u* - w_h, w b_K)_K = 0                      ∀w ∈ P_{k-1}(K),

where {{L}} is the single-valued average of the velocity gradient across the face, wh = L21 - L12 the discrete vorticity, and bK the cubic bubble. The resulting velocity is exactly divergence-free, H(div)-conforming, and converges with order k+2 for k ≥ 1.

Two-dimensional only: the 3D variant needs the vector vorticity and a basis of the face-tangential space and is not implemented yet.

master1/mesh1 hold the degree k+1 approximation on the same triangulation and must be built with the same quadrature order as master/mesh (see hdg_postprocess for the same convention; on curved meshes mesh1 must also carry the same discrete geometry — see match_geometry!). Returns u* as a (npl1 × 2 × nt) array on the nodes of mesh1.

hdg_ns_solve(master, mesh, ν, dbc; τ=1.0, source=nothing, maxiter=20,
             tol=1e-10, verbose=true, u0=nothing, Λ0=nothing, callback=nothing)

Solves the steady incompressible Navier-Stokes equations (2D or 3D) with the HDG method by Newton iteration (the first iteration, started from rest, is a Stokes solve). See hdg_ns_step for the arguments and the returned fields.

callback, if given, is called after every Newton iteration as callback((; iter, residual, u, p)) — the iteration count, the relative trace update ‖ΔΛ‖/‖Λ‖ (Inf on the first iteration), and the current velocity/pressure fields. Return true to stop the iteration early. This is the Newton-granularity analog of the time-loop callback convention.

hdg_ns_step(master, mesh, ν, dbc; τ=1.0, source=nothing, u=nothing, Λ=nothing,
            uold=nothing, dtinv=0.0)

Performs one Newton step of the HDG discretization of the incompressible Navier-Stokes equations (2D or 3D, from the mesh), linearized about the velocity u (npl × Dim × nt) and the velocity trace Λ (vector of length Dim·nps·nf).

Arguments

• ν: kinematic viscosity
• dbc: Dirichlet boundary velocity, called as dbc(p) with p the node coordinates, returning the velocity vector [g1, …, gDim]
• τ: HDG stabilization parameter (τ ≈ ν/ℓ + |u|)
• source: body force; nothing, a function p -> [f1, …, fDim], or a nodal array (npl × Dim × nt)
• uold, dtinv: previous-time-level velocity and 1/Δt for backward Euler (use dtinv = 0 for steady state)

Returns

Named tuple (u, gradu, p, Λ, ρ) with the new velocity (npl × Dim × nt), velocity gradient (npl × Dim² × nt, column (i-1)Dim + j holding Lij = ∂ui/∂xj), pressure (npl × nt, zero global mean), trace vector, and element mean pressures.

hdg_ns_step_batched(master, mesh, ν, dbc; τ=1.0, source=nothing, u=nothing,
                    Λ=nothing, uold=nothing, dtinv=0.0,
                    ArrayT=Array, T=Float64, cache=nothing)

Batched/KA counterpart of hdg_ns_step (same arguments, same returned fields plus cache): the per-element Newton assembly, local solves and (u, p, L) recovery run on the backend of ArrayT (e.g. CuArray); the condensed trace saddle-point system is solved on the CPU with a sparse LU whose sparsity pattern and factorization are reused across calls. Works on triangles and tetrahedra alike.

Pass the returned cache back in on subsequent calls (same master, mesh, ν, τ) to skip all setup and reuse the factorization.

hdg_parsolve(master, mesh, source, dbc, param;
             ArrayT=Array, T=Float64, restart=80, tol=1e-6, maxit=2000,
             preconditioner=true, verbose=false)

Solves the convection-diffusion equation using the HDG method with restarted, block-Jacobi-preconditioned GMRES (Krylov.jl) on the statically condensed trace system. The iteration runs on the KernelAbstractions backend of ArrayT (pass ArrayT=CuArray with CUDA.jl loaded for a GPU solve); element assembly and local recovery stay on the CPU.

Arguments

• master: master structure
• mesh: mesh structure
• source: source term
• dbc: dirichlet data
• param: dictionary with parameters:
  • param[:kappa]: diffusivity coefficient
  • param[:c]: convective velocity
  • param[:taud]: stabilization parameter

Returns

• uh (npl, 1, nt): approximate scalar variable
• qh (npl, 2, nt): approximate flux
• uhath (nps, nf): approximate trace
• niter: number of GMRES iterations

hdg_parsolve_batched(master, mesh, source, dbc, param;
                     ArrayT=Array, T=Float64, kwargs...)

Fully device-resident counterpart of hdg_parsolve: batched element assembly (hdg_local_solves) and local recovery (hdg_recover) run on the backend of ArrayT alongside the hdg_gmres_ka trace solve. Only the Dirichlet-BC application and the face-block global assembly (hdg_densesystem) remain on the CPU. kwargs are forwarded to hdg_gmres_ka. Works on triangles and tetrahedra alike.

Returns (uh, qh, uhath, niter) with the same shapes as hdg_parsolve (uh (npl, 1, nt), qh (npl, Dim, nt), uhath (nps, nf)).

hdg_parsolve_ka(master, mesh, source, dbc, param; kwargs...)

Alias of hdg_parsolve, kept for backwards compatibility: hdg_parsolve itself now solves the trace system with the KA/Krylov GMRES (ArrayT=CuArray for a GPU solve). Returns (uh, qh, uhath, niter).

hdg_postprocess(master, mesh, master1, mesh1, uh, qh)

Postprocesses the HDG solution to obtain a superconvergent solution (Nguyen, Peraire & Cockburn, JCP 230, 2011). Dimension-generic: works on triangles and tetrahedra alike. On meshes with curved elements, mesh1 must carry the same discrete geometry as mesh — see match_geometry!.

Arguments

• master, mesh: master structure and mesh of order porder
• master1, mesh1: master structure and mesh of order porder + 1 (built with the same quadrature rule as master)
• uh (npl, 1, nt): approximate scalar variable
• qh (npl, Dim, nt): approximate flux q = -∇u (pass q ./ κ when the solver returns κ∇u-scaled fluxes)

Returns

• ustarh (npl1, 1, nt): postprocessed scalar variable

HDG Postprocessing

Solves, element by element, the local Neumann problem: find u* in P_{p+1}(K) such that

• (∇u*, ∇v)_K = -(qh, ∇v)_K for all v in P_{p+1}(K),
• the mean of u* over K equals the mean of uh,

which converges at order p+2 for diffusion problems.

hdg_recover(batch, loc, x) -> (uh, qh)

Recovers uh (npl, nt) and qh (npl, Dim, nt) from the global trace vector x (on the same backend as batch) using the local solution maps from hdg_local_solves — a batched matvec, no per-element re-solve.

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

Solves the convection-diffusion equation using the HDG method with a direct (sparse LU) solve of the statically condensed trace system. Builds exactly the same system as hdg_parsolve (element matrices via hdg_elemmats, strong Dirichlet rows on boundary faces), so the two solvers agree to solver — not just discretization — accuracy.

Arguments

• master: Master element structure
• mesh: Mesh structure
• source: Source term function or nothing
• dbc: Dirichlet boundary condition data
• param: Dictionary with parameters :kappa (diffusivity), :c (convective velocity) and :taud (stabilization)

Returns

• uh (npl, 1, nt): Approximate scalar variable
• qh (npl, 2, nt): Approximate flux
• uhath (nps, nf): Approximate trace

hdg_trace_system(ae, fe, elcon) -> (K, F)

Assemble the global condensed trace system from the element trace blocks ae (ndf, ndf, nt), fe (ndf, nt) via the orientation-resolved trace connectivity elcon (ndf, nt): triplet-based sparse assembly (duplicate entries sum, exactly like the matrix-free face matvec).

localprob(dg, master, m, source, param)

Solves the local convection-diffusion problems for the HDG method.

Arguments

• dg: DG nodes
• master: Master element structure
• m: Values of uhat at element edges
• source: Source term function or nothing
• param: Dictionary with parameters :kappa (diffusivity) and :c (convective velocity)

Returns

• umf: Local solution uh
• qmf: Local solution qh

match_geometry!(master, mesh, master1, mesh1) -> mesh1

Replace the node coordinates of mesh1's curved elements with the degree-p isoparametric map of mesh evaluated at mesh1's reference nodes (master1.plocal) — exact, since $P_p ⊂ P_{p+1}$. After the call both meshes describe the same discrete geometry.

This is a prerequisite for the p+2 superconvergence of hdg_postprocess (and hdg_ns_postprocess) on curved elements: when mesh and mesh1 project their boundary nodes independently, their isoparametric maps differ by O(h^{p+1}), the p- and (p+1)-fields are composed with different maps of the reference element, and u* degrades to O(h^{p+1}) (no gain over u). Straight elements are unaffected (both maps are already the same affine map).

CGMatVecOp(ae, tcg, dirichlet)

Matrix-free CG stiffness operator for Krylov.jl (implements eltype, size, 3-arg mul!): y = K x as gather → batched row-matvec → atomic scatter, with the identity on Dirichlet rows. Lives on whatever backend its arrays live on.

cg_assemble(ae, fe, tcg, dirichlet) -> (K, F)

Triplet-based global assembly (sparse(I, J, V) sums duplicates): the sparse matrix K with unit diagonal on Dirichlet rows and the load vector F. Replaces the legacy incremental K[i, j] += ... insertion, which re-shuffled the CSC structure per entry.

cg_element_system(mesh, master, source, param; T=Float64, eliminate=true)

Batched CG element matrices and load vectors: ae (npl, npl, nt), fe (npl, nt), the global Dirichlet-node mask dirichlet (nn), and whether the (eliminated) operator is symmetric (c == 0).

Geometry is evaluated isoparametrically on mesh.pcg[mesh.tcg[e, :], :] — the same coordinates elemmat_cg uses (they can differ from dgnodes by cgmesh's deduplication rounding), so parity with elemmat_cg is exact. eliminate=false skips the symmetric Dirichlet elimination (used by the parity tests).

cg_parsolve(mesh, master, source, param; T=Float64, ArrayT=Array,
            tol=1e-10, maxit=5000, restart=80, preconditioner=true,
            verbose=false) -> (uh, energy, niter)

Matrix-free iterative CG solve on any KA backend: conjugate gradients when the operator is symmetric (no convection), restarted GMRES otherwise, both Jacobi-preconditioned. ArrayT=CuArray (with CUDA.jl loaded) runs the whole iteration on the GPU. Returns the same uh (npl, nt) / energy as cg_solve, plus the Krylov iteration count.

cg_solve(mesh, master, source, param) -> (uh, energy)

Solve the steady convection-diffusion(-reaction) equation with continuous Galerkin finite elements and homogeneous Dirichlet boundaries, on triangles or tetrahedra. mesh must carry the CG numbering (see cgmesh); the source term is called with the coordinates splatted (source(x, y) in 2D, source(x, y, z) in 3D); param is a named tuple (; κ, c, s) with diffusivity κ, convective velocity c (length Dim), and reaction coefficient s.

Element matrices are built batched (cg_element_system), assembled from triplets (cg_assemble), and factorized directly: Cholesky when the operator is SPD (no convection, s ≥ 0), sparse LU otherwise. For an iterative, GPU-capable solve see cg_parsolve.

Returns the solution uh (npl, nt) in DG (element-local) numbering, ready for scaplot, and the discrete energy ½ uᵀKu - uᵀF.

The high-level equivalent is solve(CGProblem(equation, mesh; source)).

elemmat_cg(pcg, master, source, param) -> (A, F)

Elemental stiffness matrix and force vector of the continuous-Galerkin convection-diffusion(-reaction) operator — the per-element reference implementation the batched path (cg_element_system) is verified against. Dimension-generic (triangles and tetrahedra).

• pcg (npl, Dim): element node coordinates
• master: reference element
• source: forcing function, called with the coordinates splatted (source(x, y) in 2D, source(x, y, z) in 3D)
• param: named tuple with diffusivity κ, convective velocity c (length Dim), and reaction coefficient s

Returns the local element matrix A (npl, npl) and force vector F (npl).

Computes the normal flux on each edge ∇u · n by averaging the gradients from the two neighboring elements and then correcting them to ensure that the integral of qn along the element boundary is equal to the integral of f on the volume.

Parameters:

• master: master structure
• mesh: mesh structure
• guh: solution gradient, size (npl, 2, nt)
• forcing: forcing function

Returns:

• qn: equilibrated normal fluxes to the face. The direction is pointing into element mesh.f(i,3)
• qn0: non-equilibrated (averaged) normal fluxes to the face

Computes the gradient of a scalar field on a mesh.

Parameters:

• master: master structure
• mesh: mesh structure
• uh: approximate scalar variable with local numbering, size (npl, nt)

Returns:

• guh: solution gradient, size (npl, 2, nt)

l2error(mesh, uh, exact)

L2 norm of the error between a numerical solution and an exact solution, computed with high-order (4p) quadrature.

Arguments

• mesh: Mesh structure
• uh: Scalar field with local (DG) numbering — either (npl, nt) or a single-component (npl, 1, nt) array
• exact: Exact solution function of the coordinates, (x, y) -> value in 2D, (x, y, z) -> value in 3D

Returns

• ‖uh - exact‖_{L²(Ω)}

Given normal equilibrated fluxes, it computes the fluxes q in the element interior satisfying - Div q = f and evaluates the lower energy bound.

Parameters:

• master: master structure
• mesh: mesh structure
• forcing: forcing function
• qn: normal gradient to the face

Returns:

• q: reconstructed element fluxes

AbstractEquation{Dim}

Supertype of all equations in Dim space dimensions (Trixi's NDIMS pattern: the dimension is a type parameter the kernels and problem constructors dispatch on). A concrete equation is a small immutable struct holding its physical parameters as typed fields (e.g. EulerEquations(γ)), and implements the dispatch contract:

States u are SVector{nc}, positions x and normals n are SVector{Dim}, gradients q are SMatrix{Dim, nc} (row d = ∂/∂x_d). All methods must be @inline, allocation-free, and generic in the element type T so they compile inside CPU/GPU kernels at any precision.

Scalar equations are dimension-polymorphic: ConvectionEquation(v) infers Dim = length(v) from its velocity; equations without a directional parameter (e.g. PoissonEquation()) default to 2D, with PoissonEquation{3} available explicitly — a user never spells Dim twice.

BoundaryCondition

Supertype of boundary conditions (Dirichlet, Neumann, SlipWall, FarField, IncomingWave, or a user type). Implement, on your own type,

boundary_state(bc::MyBC, eq, uL, n, x, t) -> SVector    # ghost state

(the numerical flux is then evaluated against it), or override

boundary_flux(bc::MyBC, eq, numerical_flux, uL, n, x, t) -> SVector

for non-ghost-state conditions. Fields must be isbits (numbers, SVectors, or closures over them) so the condition can move to a GPU.

ConvectionDiffusionEquation(velocity, κ)
ConvectionDiffusionEquation{Dim}(velocity, κ)

Linear convection-diffusion

∂u/∂t + ∇ ⋅ (v u) = ∇ ⋅ (κ ∇u) + s,

with velocity a constant Dim-vector (dimension inferred from its length) or a callable x -> SVector{Dim} (2D by default; use the explicit {Dim} form in 3D), and diffusivity κ. Discretized with LDG viscous fluxes in DG (stabilized by an LDGStabilization policy passed to the problem) and with HDG for steady problems.

ConvectionEquation(velocity)
ConvectionEquation{Dim}(velocity)

Linear scalar convection

∂u/∂t + ∇ ⋅ (v u) = s,

with velocity either a constant Dim-vector (the dimension is inferred from its length) or a callable x::SVector{Dim} -> SVector{Dim} for a spatially varying field (2D by default; use ConvectionEquation{3}(f) in 3D).

Dirichlet(value=0.0)

Prescribed solution value on a boundary: a constant, an SVector (one entry per component), or a function (x, t) -> value.

EulerEquations(; γ=1.4)
EulerEquations{Dim}(; γ=1.4)

Compressible Euler equations for the conserved state (ρ, ρu, ρv[, ρw], ρE) (Dim + 2 components; 2D by default), with ideal-gas ratio of specific heats γ. Default numerical flux: RoeFlux.

FarField(state)

Far-field boundary carrying the free-stream state (one value per component); the numerical flux against it upwinds automatically.

Incoming-wave boundary for WaveEquation (uses the equation's k, f).

LDGStabilization(c11=1.0, c11int=0.0)

Penalty coefficients of the LDG viscous fluxes (Cockburn & Shu, SINUM 35, 1998): c11 on Dirichlet boundary faces, c11int on interior faces. A swappable policy object — pass to DGProblem(...; stabilization=...).

LaxFriedrichs()

Local Lax–Friedrichs (Rusanov) flux,

F̂ = ½ (F(uL) + F(uR)) ⋅ n − ½ λ (uR − uL),   λ = max(|λ(uL)|, |λ(uR)|),

using the equation's flux and max_abs_speed. Works for any equation that implements those two methods; exact upwinding for linear scalar convection.

Neumann(flux=0.0)

Prescribed (currently homogeneous) normal flux: the solution trace is taken from the interior, and the viscous boundary flux is flux.

PoissonEquation(κ=1.0)
PoissonEquation{Dim}(κ=1.0)

Poisson / pure-diffusion equation -∇ ⋅ (κ ∇u) = s, for HDGProblem and CGProblem. Dimension-polymorphic: 2D by default (there is no directional parameter to infer from), PoissonEquation{3}() in 3D.

RoeFlux()

Roe's approximate Riemann solver (Roe, JCP 43:357, 1981),

F̂ = ½ (F(uL) + F(uR)) ⋅ n − ½ |Â(ũ)| (uR − uL),

with Â the flux Jacobian evaluated at the Roe average ũ. Implemented for EulerEquations (density-weighted Roe average) and WaveEquation (exact characteristic decomposition of the linear system).

Impermeable slip wall (normal-velocity reflection) for the wave and Euler systems.

WaveEquation(c; k=nothing, f=nothing)

First-order acoustic wave system for u = (q₁, …, q_Dim, p),

∂q/∂t + c ∇p = 0,    ∂p/∂t + c ∇ ⋅ q = 0,

with speed c (Dim + 1 components; 2D by default, or inferred from the wave vector k). k and f(c, k, x, t) prescribe the incident field and are only needed when an IncomingWave boundary is used.

boundary_flux(bc, eq, numerical_flux, uL, n, x, t) -> SVector

Normal boundary flux: by default the numerical flux evaluated against the condition's boundary_state. Override for conditions that prescribe the flux itself.

boundary_state(bc, eq, uL, n, x, t) -> SVector

Ghost (exterior) state a boundary condition presents to the numerical flux: the boundary data for Dirichlet/FarField, the interior state for Neumann, the reflected state for SlipWall, ….

boundary_trace(bc, eq, uL, n, x, t) -> SVector

Solution trace û a boundary condition prescribes for the LDG gradient: the boundary data for Dirichlet, the interior trace for Neumann.

boundary_viscous_flux(bc, stab, eq, uL, qL, n, x, t) -> SVector

Viscous flux on a boundary face: for Dirichlet data g, -κ qL ⋅ n + c11 (uL - g); for Neumann, the prescribed flux (currently homogeneous).

canonical_to_riemann(ρ, ρu₁, ρu₂, ρE, γ) -> (v, s, J⁺, J⁻)

Convert conserved Euler variables to the 1D Riemann-invariant variables: tangential velocity v, entropy s = p/ρ^γ, and invariants J± = u ± 2c/(γ-1). See riemann_to_canonical.

default_numerical_flux(eq)

The numerical flux solve uses when the problem does not specify one (e.g. RoeFlux() for EulerEquations, LaxFriedrichs for scalar convection). Any callable (eq, uL, uR, n, x, t) -> SVector{nc} can replace it per problem.

default_stabilization(eq)

The viscous (LDG) stabilization solve uses when the problem does not specify one; nothing for purely hyperbolic equations.

Density ρ of a conserved Euler state.

derived_field(f, eq, u) -> Matrix

Evaluate a pointwise derived quantity f(eq, u::SVector) -> Real (e.g. pressure, mach) over a solution field u (npl, nc, nt); returns (npl, nt).

diffusivity(eq) -> Real

Scalar diffusivity κ of the equation's second-order term; 0 for purely hyperbolic equations (the default). CFL estimates charge it against the quadratic (2p+1)²/h² limit.

Internal energy density ρe = ρE - ρ|v|²/2 of a conserved Euler state.

Kinetic energy density ρ|v|²/2 of a conserved Euler state.

Total energy density ρE of a conserved Euler state.

Entropy function s = p/ρ^γ.

eulereval(u, str, γ) -> Matrix

String-keyed Euler derived quantities for plotting scripts — a thin lookup over the dispatched primitives (density, pressure, mach, …): "r", "u", "v", "p", "c", "M", "s", and the 1D characteristic combinations "Jp"/"Jm" (which, as in the original plotting scripts, read component 2 as the Riemann u rather than dividing by density).

flux(eq, u, x, t) -> NTuple{Dim, SVector{nc}}

Physical (volume) flux of the equation at state u::SVector{nc} and position x::SVector{Dim}: one SVector{nc} flux per space direction, such that the conservation law reads ∂u/∂t + Σ_d ∂f_d/∂x_d = source (in 2D this is the familiar (fx, fy) pair). Defined once per equation; the numerical fluxes, boundary fluxes, and volume terms all call it.

has_diffusion(eq) -> Bool

Whether the equation has second-order (viscous/diffusive) terms, i.e. whether the LDG gradient and viscous-flux path must run. Defaults to false.

Mach number |v|/c.

max_abs_speed(eq, u, n, x, t) -> Real

Maximum absolute characteristic speed of eq at state u in direction n (e.g. |v ⋅ n| for convection, |u ⋅ n| + c for Euler). Used by LaxFriedrichs-type dissipation and CFL estimates.

normal_flux(eq, u, n, x, t) -> SVector{nc}

Physical flux projected on the unit normal n: flux(eq, u, x, t) ⋅ n. Numerical fluxes use it for their central part.

nvariables(eq) -> Int

Number of conserved components of an equation (or of a physics bundle that wraps one).

Ideal-gas pressure p = (γ-1)(ρE - ρ|v|²/2) of a conserved Euler state.

riemann_to_canonical(v, s, J⁺, J⁻, γ) -> (ρ, ρu₁, ρu₂, ρE)

Convert the 1D Riemann-invariant variables — tangential velocity v, entropy s = p/ρ^γ, and invariants J± = u ± 2c/(γ-1) — to conserved Euler variables. Inverse of canonical_to_riemann; used to prescribe characteristic far-field states.

Speed of sound c = √(γp/ρ).

varnames(eq) -> NTuple{nc, Symbol}

Names of the conserved components, in storage order.

Velocity vector ρv/ρ of a conserved Euler state.

viscous_flux(eq, u, q, x, t) -> (fx, fy)

Viscous physical flux at state u and gradient q::SMatrix{2, nc} (rows are the x/y derivatives). Only needed when has_diffusion is true.

viscous_numerical_flux(stab, eq, uL, uR, qL, qR, n, x, t) -> SVector

LDG viscous interface flux: the alternating choice takes the trace û from the left and the gradient from the right element, plus the interior penalty: F̂ᵥ = -κ qR ⋅ n + c11int (uL - uR).

wavespeed(eq, u) -> Real

Direction-independent bound on the characteristic speed at state u::SVector — |v| + c for Euler, |v| for constant-velocity convection, |c| for the wave system, 0 for pure diffusion. CFL estimates (compute_dt, StepsizeCallback) take its maximum over the solution nodes. Equations whose speed depends on position (a velocity given as a function of x) cannot provide a state-only bound and throw.

meshplot(mesh; nodes=false, annotate="")

Plot a simplicial mesh. nodes=true marks the DG nodes; annotate may contain 'p' (number the vertices) and/or 't' (number the elements). Requires a Makie backend to be loaded (e.g. using CairoMakie).

meshplot_curved(mesh; nodes=false, annotate="", pplot=0, figure_size=(800, 800), title="")

Plot a curved (isoparametric) mesh, subdividing each element with order pplot for display. Requires a Makie backend to be loaded (e.g. using CairoMakie).

save_vtk(mesh, u, filename; names=nothing) -> saved file paths
save_vtk(sol, filename)

Write a solution as high-order Lagrange VTK cells for ParaView (the 3D visualization path; works for 2D meshes too). u is a (npl, nc, nt) (or (npl, nt) scalar) field; names optionally labels the components. The solution-object form names components by varnames(sol.prob.equation). Curved (isoparametric) elements are written with their true high-order geometry, which ParaView renders natively.

Implemented in the TwoDGWriteVTKExt package extension: available once WriteVTK.jl is loaded (using TwoDG, WriteVTK).

scaplot(mesh, c; limits=nothing, show_mesh=false, figure_size=(800, 800), title="", cmap=:turbo)

Plot contours of a scalar field c of shape (npl, nt) on mesh. Requires a Makie backend to be loaded (e.g. using CairoMakie).

initu(mesh, nc::Integer, value)

Initialize the vector of unknowns u (npl, nc, nt) from per-component constants or functions (see interpolate, which this wraps). Pass nvariables(equation) for nc.

interpolate(mesh, values)

Nodal interpolation of initial/reference data onto the DG nodes. values is a collection with one entry per solution component; each entry is either a number (constant field) or a function of the coordinates ((x, y) -> value in 2D, (x, y, z) -> value in 3D). Returns u (npl, nc, nt) with nc = length(values).

newton_raphson(f, fgrad, x₀; abstol=1e-8, reltol=1e-8, maxiter=100)

Newton-Raphson method for root finding.

unique_rows(A; return_index=false, return_inverse=false)

Sorted unique rows of the matrix A, like NumPy's np.unique(A, axis=0). With return_index=true also returns the row indices I of the first occurrences (B = A[I, :]); with return_inverse=true also returns the inverse map J reconstructing A (A = B[J, :]).

areacircle calcualte the area and perimeter of a unit circle [area,perim]=areacircle(sp,porder)

siz: desired element size porder: polynomial order of approximation (default=1) area1: area of the circle (π) area2: area of the circle (π) perim: perimeter of the circumference (2π)

potential_trefftz(x, y; V=1.0, alpha=0.0, tparam=[0.1, 0.05, 1.98])
    -> (psi, velx, vely, gamma)

Analytic 2D potential-flow solution around a Karman–Trefftz airfoil, evaluated at the points (x, y): stream function psi, velocity components velx/vely, and circulation gamma (lift force = V * gamma). V is the free-stream speed, alpha the angle of attack in degrees, and tparam = [x0, y0, exponent] the airfoil parameters (circle-center shifts, K–T exponent ≤ 2 with 2 giving a Joukowski airfoil).

trefftz(V∞, α, m=15, n=30, porder=3, node_spacing_type=0,
        tparam=[0.1, 0.05, 1.98])

Full Karman–Trefftz potential-flow study: build the curved airfoil mesh (mkmesh_trefftz), evaluate the analytic stream function (potential_trefftz), differentiate it for the velocity field, and integrate surface pressure for the aerodynamic coefficients.

Returns the 20-tuple (mesh, master, xs_foil, ys_foil, chord, ψ, vx, vy, Γ, CP, CF, Clift, Cdrag, CM, vx_analytical, vy_analytical, Γ_analytical, CP_analytical, CL, CL_analytical) — numerical and analytic velocities, circulation, pressure/force/moment coefficients. See examples/dg/runtrefftz.jl.

trefftz_points(tparam=[0.1, 0.05, 1.98], np=120) -> (x, y, chord)

Coordinates of np points on a Karman–Trefftz airfoil surface and its chord length. tparam = [x0, y0, exponent]: circle-center shifts and K–T exponent (2 gives a Joukowski airfoil, trailing edge at (exponent, 0)).