Stokes, Brinkman and Oseen
Incompressible-flow formulations, residual estimators and adaptive policies.
User-defined formulations compose these public numerical owners with
LocalEquations and Equation; their editable providers are listed in the
variational guide.
Physical field records provide optional interpretation of the resulting
coefficients. Named FieldDefinition views require no physical solver object.
The volume operator is the vector Laplacian. The normal-oriented hybrid
multiplier is \(\lambda_K=(-\nu\nabla u+pI+u\otimes\beta/2)n_K\);
the traction convenience arguments prescribe its negative. With zero
convection the half-advection term vanishes. Velocity/pressure coefficient
arrays require their executed bases to recover these physical fields.
The API below documents its kernels, pressure gauge and residual signs.
For the analytical hypotheses, consult
Araya et al. (2017) and
Araya et al. (2025).
Successful matrix assembly does not establish those hypotheses.
Physical velocity-pressure volume forms on explicitly supplied triangular meshes.
FlowOperators
dataclass
Velocity first, pressure second; literal volume operator and physical moments.
minimum_resistance(drag: Any, supplied: float | None) -> float
Determine and validate an explicit global lower eigenvalue bound for 2017 USFEM.
advection_contract(advection: Any, divergence: Any, bound: float | None, *, stabilized: bool) -> tuple[Any, Any, float | None]
Validate exact divergence data and an optional certified velocity bound.
A callable needs its analytical divergence, including explicit zero for a solenoidal field. Stabilized Oseen also needs a supplied global upper bound on its Euclidean magnitude; quadrature sampling cannot certify a supremum.
triangle_flow_operators(fine: TriangleMesh, *, degree: int, formulation: str, viscosity: float, drag: Any, beta: Any, source: Any, order: int, gamma_min: float | None = None, pointwise: bool = False, residual_weight: Any = None, beta_divergence: Any = 0.0, beta_bound: float | None = None) -> FlowOperators
Integrate grad-grad flow and the declared residual test/trial operators.
Coordinates are interleaved continuous Pk velocity then continuous pressure, P(k-1) for Taylor--Hood or Pk for USFEM/Oseen. Pressure equations have negative divergence. Skew advection includes -div(beta)/2; Oseen residuals use opposite test/trial convection and positive grad-div. Cartesian material cut integration preserves physical resistance data; inverse bounds use independent Gaussian polynomial integration. No trace, boundary, kernel elimination or solve is chosen.
Three-dimensional grad-grad flow forms and residual stabilization on tetrahedra.
Flow3DOperators
dataclass
Assembled velocity-pressure algebra and physical moments on one local mesh.
Velocity coefficients are interleaved Cartesian components at continuous Pk nodes, followed by scalar continuous pressure coefficients. Taylor--Hood uses P(k-1) pressure; stabilized forms use Pk. Material moments support exact-zero translation tests; sampled positive resistances are never classified by a small absolute threshold.
resistance_values_3d(field: Any, points: FloatArray) -> FloatArray
Evaluate finite symmetric nonnegative scalar or 3-by-3 resistance tensors.
flow_contract_3d(viscosity: Any, drag: Any, advection: Any, divergence: Any, bound: Any, formulation: str, stabilization: str, gamma_min: Any, degree: int | None) -> tuple[float, Any, Any, float | None, float | None, int]
Validate the common physical and polynomial contracts for 3D flow forms.
tetra_flow_operators(mesh: TetraMesh, *, viscosity: float = 1.0, drag: Any = 0.0, advection: Any = (0.0, 0.0, 0.0), advection_divergence: Any = None, advection_bound: float | None = None, source: Any = (0.0, 0.0, 0.0), degree: int | None = None, formulation: str = 'taylor-hood', stabilization: str = 'tensor-2025', gamma_min: float | None = None, order: int = 5) -> Flow3DOperators
Assemble -nuDelta(u)+beta.grad(u)+Gammau+grad(p)=f, div(u)=0.
The pressure test equation is multiplied by -1. Advection uses its skew form minus div(beta)u.v/2; the natural traction is (nugrad(u)-pI)n-beta.nu/2. USFEM subtracts the full trial/test momentum residual product and its forcing term. Oseen uses the adjoint convection sign in the test residual, plus positive grad-div stabilization.
The polynomial inverse bound m=min(1/3,C) is computed per tetrahedron from Ch^2||Delta(v)||^2 <= ||grad(v)||^2. The three USFEM coefficient choices preserve the declared 2D conventions. For variable materials, Gaussian quadrature is used; discontinuities must be resolved by the local mesh. No tetrahedral Cartesian cut integration is implied.
Finite symmetric nonnegative Brinkman resistance evaluation in two dimensions.
resistance_values(field: Any, points: FloatArray) -> FloatArray
Evaluate a symmetric nonnegative scalar or tensor Brinkman resistance.
Two-level residual estimators for Stokes, Brinkman and Oseen velocity-pressure fields.
FlowEstimator
dataclass
Quadrature-evaluated first/second-level residual contributions.
face_squared stores eta_1,Ftilde², with the original macroface length
in the denominator. Interior faces count twice in eta_1, as in the sum
over macrocell boundaries. local_squared stores unscaled eta_2,K².
fine_squared partitions each eta_2,K² among its fine cells, splitting
an interior fine-face contribution equally between its two neighbors.
Oseen's reported eta_2 includes 2*(-2ell); the Stokes/Brinkman estimator
does not. These are estimators up to the constants/higher-order terms in
their respective theorems, not certified numerical upper bounds.
eta1: float
property
Return the jump/boundary mismatch indicator with macrocell multiplicity.
eta2: float
property
Return the literature-specific scaled second-level residual indicator.
total: float
property
Return eta_1+eta_2, not the root of their squared sum.
mixed_error(velocity: Any, velocity_gradient: Any, pressure: Any, *, order: int = 8) -> float
Integrate ||(u-uh,p-ph)||_(V×Q), V²=diam(Omega)^(-2)L2²+H1².
estimate_flow_error(solution: VectorSolution, *, viscosity: float = 1.0, drag: Any = 0.0, advection: Any = (0.0, 0.0), source: Any = (0.0, 0.0), dirichlet: Any = (0.0, 0.0), full_dirichlet: bool, variant: Literal['oseen-2021', 'stokes-brinkman-2021'] = 'oseen-2021', quadrature_order: int = 8) -> FlowEstimator
Evaluate the two-level Oseen or Stokes--Brinkman estimator.
Use equations (33)--(34) of Araya et al. (2021) or equations (4.1)--(4.7) of Araya, Rebolledo and Valentin (2021).
Supply the same physical coefficients and Dirichlet data used in the solve. The formula assumes full Dirichlet boundaries, constant positive viscosity, and uniform skeletal polynomial degree ell. Piecewise trace partitions and nonmatching local submesh refinements are supported; edge integration is cut at every fine and skeletal breakpoint. The supplied full_dirichlet declaration is explicit because VectorSolution does not retain boundary data. Volume integration is Gaussian; a discontinuous material must therefore be aligned with local cells for this estimator. The Oseen theorem additionally assumes gamma-div(beta)/2>0. A finite estimator for a case outside that hypothesis is not a reliability proof.
Fixed-macro-mesh face and local refinement for the two-level flow estimators.
AdaptiveFlowResult
dataclass
Solved fields and indicators at each adaptive state, with an explicit stop reason.
mark_flow_faces(estimator: FlowEstimator, theta: float = 0.5, *, variant: Literal['oseen-2021', 'stokes-brinkman-2021'] = 'oseen-2021') -> tuple[tuple[np.ndarray, ...], np.ndarray]
Apply face marking with local-error dominance.
Use Algorithm 1 of Araya et al. (2021) or Algorithm 2 of Araya, Rebolledo and Valentin (2021).
Both algorithms use eta_F=eta_1,F+sum eta_2,K, with eta_1,F=sqrt(sum_segments eta_1,Ftilde²). Mark faces above theta times the maximum and, within each marked face, every segment tied for its largest indicator. Refine incident local meshes when eta_1,F < sum eta_2,K. Values of eta_2,K are unscaled, exactly as in those algorithms.
adapt_flow(mesh: TriangleMesh, *, skeleton: SkeletonSpace | None = None, local_refinement: int | tuple[int, ...] = 1, iterations: int = 4, theta: float = 0.5, tolerance: float = 0.0, max_local_refinement: int = 64, estimator_variant: Literal['oseen-2021', 'stokes-brinkman-2021'] = 'oseen-2021', estimator_order: int = 8, local_refiner: Literal['uniform', 'longest-edge'] = 'uniform', max_local_cells: int = 65536, local_error_marking: Literal['uniform', 'maximum'] = 'uniform', solve_step: Callable[..., VectorSolution] | None = None, **solve_options: Any) -> AdaptiveFlowResult
Solve, estimate and refine traces/local meshes with a fixed macro topology.
Defaults are local P3/P3 Oseen with discontinuous P1 traces. Supply physical coefficients and
forcing through solve_flow keyword arguments. Full Dirichlet data are required by these
published estimators. The initial local meshes must resolve the supplied skeleton. Each
subsequent local grid doubles on local-error dominance and is additionally refined to contain
every new skeletal breakpoint, as in the matching submesh construction in Section 5.3 of
Araya et al. (2021). This rational uniform-grid
closure is explicit, not a minimal unstructured refinement. The maximum local resolution bounds
memory consumption and returns a distinct stop reason; it never changes tolerances or solver
choice. local_refiner='longest-edge' instead closes new boundary breakpoints with local
longest-edge propagation; local-error dominance red-refines that local mesh once.
max_local_cells limits this nonuniform route. Its refinements result is empty because a
uniform subdivision count does not describe these meshes; actual meshes are stored in each
solution. With longest-edge closure, local_error_marking='maximum' selects cells whose eta_2
contribution is at least theta times the local maximum, within every local mesh triggered by the
published macroface marking. Interior fine-face terms are split equally between their two
neighboring cells. This within-local policy is explicit; it is not prescribed by
Araya, Rebolledo and Valentin (2021).
solve_step(mesh, skeleton=current_space, **problem_data) may assemble
the user's own mathematical equations at each state. It receives the actual
transferred boundary data and local partitions; its returned physical field
must satisfy the estimator's existing hypotheses. The default retains the
built-in solve. Callback exceptions propagate without changing refinement.
Macroelement marking for the Stokes--Brinkman estimator.
The marking rule follows Algorithm 1 of Araya, Rebolledo and Valentin (2021).
AdaptiveFlowMacroResult
dataclass
Solved macro meshes, original indicators and conforming refinement ancestry.
mark_flow_cells(estimator: FlowEstimator, theta: float = 0.5) -> np.ndarray
Mark eta_K >= theta*max(eta_K), where eta_K=eta_2,K+sum_F eta_1,F.
The first-level contribution on each macroface is the square root of its summed segment contributions. Interior macrofaces contribute to both cells. Local indicators are unscaled; this follows Algorithm 1 of Araya, Rebolledo and Valentin (2021), rather than bulk marking. A zero estimator gives an empty marking set.
adapt_flow_macros(mesh: TriangleMesh, *, trace_degree: int = 0, local_refinement: int | tuple[int, ...] = 1, iterations: int = 4, theta: float = 0.5, tolerance: float = 0.0, maximum_cells: int = 10000, estimator_order: int = 8, macro_refiner: Literal['red-green', 'longest-edge'] = 'red-green', solve_step: Callable[..., VectorSolution] | None = None, **solve_options: Any) -> AdaptiveFlowMacroResult
Apply macro marking with an explicitly selected conforming refinement.
The marking rule follows Algorithm 1 of Araya, Rebolledo and Valentin (2021).
iterations counts refinement steps after the initial solve. Marked triangles use red--green
closure by default. macro_refiner='longest-edge' selects longest-edge propagation,
preserving the same marking rule and controlling the minimum angle under repeated local
refinement. Each child inherits its parent's local subdivision count, so no additional local
refinement is induced. Local physical cells shrink only with their parent macrocell. Traces are
unsplit discontinuous P_ell on every current macroface. This fully specified closure is an
original realization of the published strategy; the paper does not provide its historical
connectivity.
The estimator assumes full Dirichlet Stokes/Brinkman data. Defaults are local P3/P3 USFEM and P0
traces. solve_options carries the same physical data to both solver and estimator; advection
and prescribed tractions are excluded. A cell cap stops before solving an oversized mesh. No
convergence rate or unit reliability constant is asserted by this adaptive loop.
solve_step(mesh, skeleton=current_space, **problem_data) may assemble
the user's own mathematical equations at each state. It receives the actual
transferred boundary data and local partitions; its returned physical field
must satisfy the estimator's existing hypotheses. The default retains the
built-in solve. Callback exceptions propagate without changing refinement.
Broken velocity, pressure and grad-grad pseudostress on tetrahedral local meshes.
The scalar skeleton supplies face geometry and scalar modes. Its hybrid
trace interleaves three canonical negative-traction coefficients per mode.
Velocity and pressure are continuous only inside each macrocell. The raw
pseudostress is not claimed to be globally H(div)-conforming.
evaluate(cell: int, bary: FloatArray) -> tuple[FloatArray, FloatArray]
Evaluate one-sided velocity and pressure on all fine tetrahedra of a macrocell.
gradient(cell: int, bary: FloatArray) -> FloatArray
Evaluate grad(u)[component, derivative] without averaging macro interfaces.
pseudostress(cell: int, bary: FloatArray) -> FloatArray
Return -nugrad(u)+pI+u tensor beta/2, whose normal is the MHM multiplier.
l2_error(exact: Any, order: int = 7) -> float
Return the complete broken velocity L2 error.
pressure_l2_error(exact: Any, order: int = 7) -> float
Return pressure L2 error with the explicitly imposed physical gauge.
h1_seminorm_error(exact_gradient: Any, order: int = 7) -> float
Return the complete broken velocity-gradient L2 error.
divergence_l2(order: int = 7) -> float
Integrate the raw pointwise velocity divergence; no H(div) recovery is implied.
References
-
Rodolfo Araya, Ramiro Rebolledo, and Frédéric Valentin (2021). On a multiscale a posteriori error estimator for the Stokes and Brinkman equations, IMA Journal of Numerical Analysis 41(1), 344–380. DOI: 10.1093/imanum/drz053. An earlier version is HAL: hal-01945934v1.
-
Rodolfo Araya, Cristian Cárcamo, Abner H. Poza, and Frédéric Valentin (2021). An adaptive multiscale hybrid-mixed method for the Oseen equations, Advances in Computational Mathematics 47, article 15. DOI: 10.1007/s10444-020-09833-8.