Multivariate Functional Approximation (MFA)

• Multivariate Functional Approximation is a framework for representing complex, multivariate data with continuous surrogate functions that enable evaluation beyond the sampled points.
• It encompasses a range of methods — including spline-based, neural, and polynomial approaches — each tailored for specific applications like scientific visualization, machine learning, or statistical analysis.
• MFA techniques offer practical benefits in compression and derivative evaluation, while balancing trade-offs between local simplicity and global fidelity in large-scale data modeling.

Multivariate functional approximation (MFA) denotes a family of methods that replace sampled multivariate data by a functional surrogate that can be evaluated away from the original observations. In the literature, the term has two closely related but not identical uses. In scientific computing and visualization, MFA is a specific continuous, meshless, spline-based representation of discrete data, typically built from tensor-product B-splines or rational B-splines with control points and knots, and designed to support compact storage together with high-order value and derivative queries anywhere in the domain (Sun et al., 2022). In approximation theory, machine learning, and statistics, the same term also refers more broadly to the problem of approximating a multivariate function by polynomial, Fourier, tensor, rational, neural, or stochastic-process constructions, often under low-rank, sparsity, or smoothness assumptions (Wu et al., 4 Mar 2026). Across these strands, the shared objective is to represent a function or field in a form that is more compact, more analyzable, or more adaptive than the original samples.

1. Scope and terminological range

The literature does not use MFA as the name of a single algorithm. Rather, it names a class of approximation problems and representations. One strand treats MFA as a scientific-data model for large volumetric fields, emphasizing compactness, analytical differentiation, and direct evaluation of values and gradients in structured and unstructured domains (Sun et al., 2023). A second strand uses MFA to describe the approximation of multivariate functions on domains such as $[0,1]^n$ or $[-1,1]^d$ by fixed or learned basis systems, low-rank tensor factorizations, or rational surrogates (Potts et al., 2019). A third strand embeds multivariate functional approximations inside statistical models for multivariate functional responses, random effects, and longitudinal tensors (Volkmann et al., 2021).

Usage in the literature	Core representation	Representative papers
Scientific-data MFA	Tensor-product B-splines or rational B-splines with control points and knots	(Sun et al., 2022, Sun et al., 2023, Ma et al., 2024)
Approximation-theoretic MFA	Polynomial, Fourier, Gaussian, rational, or neural expansions	(Potts et al., 2019, Rodriguez et al., 2020, Wu et al., 4 Mar 2026, Yadav et al., 12 Jan 2026)
Statistical MFA	Truncated KL/MFPCA bases, tensor decompositions, spline factor models	(Volkmann et al., 2021, Alam et al., 2024, Kowal et al., 2014)

This suggests that MFA is best understood as an umbrella concept whose concrete meaning depends on the modeling objective. In scientific visualization, the emphasis is on a continuous surrogate of discrete fields. In machine learning, the emphasis is on adaptive basis design and approximation power. In functional data analysis, the emphasis is on parsimonious representations that preserve covariance structure, supervision, or time dependence.

2. Representation families and mathematical forms

A canonical fixed-basis formulation approximates a multivariate function by a tensor-product expansion,

$$f(\mathbf{x}) \approx \sum_{k_1=1}^{K_1}\cdots\sum_{k_n=1}^{K_n} w_{\mathbf{k}} \prod_{m=1}^n \phi^{(m)}_{k_m}(x_m), \qquad \mathbf{x}=(x_1,\dots,x_n)\in[0,1]^n.$$

In this form, the approximation quality depends on the chosen basis family and the decay of the coefficients. The neural-approximation literature treats this as the baseline from which polynomial, Fourier, wavelet, and Gaussian constructions depart, and argues that predefined bases can be a poor fit for complex, nonperiodic, locally detailed, or heterogeneous signals (Wu et al., 4 Mar 2026).

A distinct but related representation arises in scientific-data MFA. There the approximation is a tensor-product spline. In $d$ dimensions, one representative form is

$$F(u_1,\cdots,u_d)= \sum_{j_1}\cdots \sum_{j_d} N^{(1)}_{j_1,p}(u_1)\cdots N^{(d)}_{j_d,p}(u_d)\, P_{j_1,\cdots,j_d},$$

where the $N^{(\ell)}_{j_\ell,p}$ are B-spline basis functions and the $P_{j_1,\cdots,j_d}$ are control points. This produces a piecewise polynomial, piecewise smooth, meshless continuous function whose derivatives can also be expressed in spline form (Ma et al., 2024).

A third family replaces fixed analytic bases by learned ones. NeuApprox introduces a neural basis function

$$f_\theta(x)=\boldsymbol W_L\Big(\cdots \sigma(\boldsymbol W_0 x+\boldsymbol b_0)\cdots\Big)+\boldsymbol b_L,$$

and organizes the multivariate approximation as a sum of block terms,

$$\psi(x_1,\dots,x_n)\approx \sum_{j=1}^{T}\mathcal{C}^j\times_1 f_{\theta_1^j}(x_1)\times_2 f_{\theta_2^j}(x_2)\cdots\times_n f_{\theta_n^j}(x_n).$$

The paper interprets each block term as a low-rank separable component, notes that the construction becomes CP-style when $R_i^j=1$ for all modes and Tucker-like when $[-1,1]^d$0, and proves that the model is dense in $[-1,1]^d$1 (Wu et al., 4 Mar 2026).

Rational approximation provides yet another form. The p-AAA framework constructs a multivariate barycentric rational approximant, for example in two variables,

$[-1,1]^d$2

The method is data-driven in the sense that it requires only function evaluations, not a full state-space model, and adapts the support points greedily from the sample set (Rodriguez et al., 2020).

Polynomial-ANOVA methods occupy a further corner of the design space. They approximate $[-1,1]^d$3 in a complete orthonormal polynomial basis and use the classical ANOVA decomposition to organize coefficients by variable-interaction support. The crucial structural claim is that many high-dimensional functions have a low-dimensional structure, meaning that only a small number of low-order interactions contribute significantly; this permits grouped least-squares reconstruction from scattered data and active-set detection through estimated global sensitivity indices (Potts et al., 2019).

3. Spline-based MFA in scientific computing

In the visualization and HPC literature, MFA is a compact, adaptive, global functional approximation of discrete scientific data. The model is built from tensor products of B-spline functions or, in the direct volume-rendering work, from tensor products of nonuniform rational B-spline functions. The basis is defined by control points and knots, and the encoded field can be queried directly at arbitrary spatial locations for values, gradients, and higher derivatives up to the polynomial degree (Sun et al., 2022).

The encoding workflow is adaptive. The data are parameterized, an initial knot distribution is established with the minimum required number of control points, and additional control points and knots are added iteratively until the approximation error in each knot span is within a specified maximum relative error, denoted $[-1,1]^d$4 or $[-1,1]^d$5 in the cited papers (Sun et al., 2023). The resulting model is global rather than local, and the knot layout reflects data complexity rather than sample spacing. This separates model structure from the distribution of input samples and distinguishes MFA from local sliding-window filters (Sun et al., 2022).

Several properties recur across the spline-based literature. MFA is described as a compact closed-form function for compression, analysis, and visualization; as a continuous, meshless surrogate model; and as a representation with high-order continuity across the domain (Mahadevan et al., 2022). Because the representation is piecewise polynomial, its derivatives are available analytically, which is central for shading, feature detection, curvature-based transfer functions, PDE residual evaluation, and topology extraction (Ma et al., 2024).

Parallelization is a major concern because single-process encoding becomes prohibitive on large datasets. A domain-decomposition approach based on overlapping additive Schwarz and an outer Restricted Additive Schwarz iteration partitions the domain into subdomains, solves constrained local MFA problems in parallel, and enforces continuity across interfaces during the solve rather than by post-processing blending (Mahadevan et al., 2022). The method uses tensor-product B-spline/NURBS representations, distinguishes clamped from floating knots, and states that floating knots allow recovery of full high-order continuity $[-1,1]^d$6 across subdomain boundaries (Mahadevan et al., 2022). The reported implementation uses Python with DIY, pyDIY, MPI-3, and MPI4Py, and the paper reports strong and weak scaling studies on the Theta Cray XC40 (Mahadevan et al., 2022).

This spline-based strand of MFA therefore combines three roles that are often separated in other literatures: approximation, compression, and differentiable surrogate modeling. A plausible implication is that its significance lies not only in reconstruction accuracy but also in replacing the original sampled field as the primary object of downstream computation.

4. Rendering, topology, and physics-informed analysis

The most direct operational use of scientific-data MFA is direct evaluation during visualization. MFA-DVR is described as the first direct volume-rendering pipeline utilizing the MFA model for both structured and unstructured volumetric datasets (Sun et al., 2022). The pipeline encodes the raw dataset offline into an MFA model, then performs ray casting by sampling points along each ray, normalizing sample positions to $[-1,1]^d$7, querying the model for scalar values and, when shading is enabled, querying the gradient as well. Sample colors are composited front-to-back, and early termination is triggered once accumulated opacity exceeds $[-1,1]^d$8 (Sun et al., 2022). Because both value and gradient come directly from the model, the renderer avoids separate gradient tables and finite-difference recomputation. The paper reports that MFA-DVR generates more faithful volume rendering than local filters, achieves the best image-quality scores among the tested methods, obtains the lowest MSE and the highest PSNR and SSIM for Marschner-Lobb gradient accuracy, is slower than trilinear interpolation, and is faster than tricubic and Catmull-Rom in the reported implementation (Sun et al., 2022).

Scalability is the central issue addressed by MFA-DVV, described as the first scalable interactive volume-visualization pipeline for large MFA models (Sun et al., 2023). Its six-stage workflow consists of volume partitioning, MFA encoding, data fetching, MFA-VV rendering, image compositing, and image merging. The system binary-partitions the volume into $[-1,1]^d$9 blocks, renders partitions in parallel, composites with the associative over operator using binary-swap compositing, and optionally uses compressed MFA models obtained by reducing the number of control points (Sun et al., 2023). One experiment reports that a compressed MFA model can be as small as $$f(\mathbf{x}) \approx \sum_{k_1=1}^{K_1}\cdots\sum_{k_n=1}^{K_n} w_{\mathbf{k}} \prod_{m=1}^n \phi^{(m)}_{k_m}(x_m), \qquad \mathbf{x}=(x_1,\dots,x_n)\in[0,1]^n.$$0 of the original dataset while still maintaining strong rendering quality (Sun et al., 2023).

Topology extraction extends the same continuous-surrogate idea beyond rendering. CPE-MFA is described as the first critical point extraction framework designed for MFA models of large-scale, high-dimensional data (Ma et al., 2024). It extracts critical points directly from the MFA model without discretization or resampling, using three stages: span filtration, Newton-based search within candidate spans, and duplicate removal via spatial hashing. Span filtration exploits the strong convex hull property of B-splines, because if the convex hull of derivative control points in a span does not contain zero, then the derivative cannot vanish there (Ma et al., 2024). On the 2D Schwefel test, the method correctly identifies all 900 critical points, partitioned into 225 minima, 225 maxima, and 450 saddles; the same paper also states that the method handles only isolated critical points, does not currently handle boundary critical points robustly, and spends more than 95% of runtime in Newton iterations (Ma et al., 2024).

A further extension extracts contours, Jacobi sets, and ridge-valley graphs directly from MFA (Ma et al., 11 Aug 2025). Contours are traced by integrating a tangent field with RK4; Jacobi set extraction is reduced to zero-contour extraction of $$f(\mathbf{x}) \approx \sum_{k_1=1}^{K_1}\cdots\sum_{k_n=1}^{K_n} w_{\mathbf{k}} \prod_{m=1}^n \phi^{(m)}_{k_m}(x_m), \qquad \mathbf{x}=(x_1,\dots,x_n)\in[0,1]^n.$$1; and ridge-valley extraction is treated as a special Jacobi set for $$f(\mathbf{x}) \approx \sum_{k_1=1}^{K_1}\cdots\sum_{k_n=1}^{K_n} w_{\mathbf{k}} \prod_{m=1}^n \phi^{(m)}_{k_m}(x_m), \qquad \mathbf{x}=(x_1,\dots,x_n)\in[0,1]^n.$$2 and $$f(\mathbf{x}) \approx \sum_{k_1=1}^{K_1}\cdots\sum_{k_n=1}^{K_n} w_{\mathbf{k}} \prod_{m=1}^n \phi^{(m)}_{k_m}(x_m), \qquad \mathbf{x}=(x_1,\dots,x_n)\in[0,1]^n.$$3, requiring derivatives up to third order (Ma et al., 11 Aug 2025). The paper states that the framework is generalizable to any continuous implicit model that supports function values and high-order derivatives. It also notes several limitations: connection ambiguity, gaps in tracing, dependence on step size, and dependence on initial sampling (Ma et al., 11 Aug 2025).

PI-MFA adds physical constraints directly to reconstruction (Jung et al., 9 Jun 2026). Using tensor-product B-splines and exact analytical derivatives, it optimizes control points against a composite loss balancing data fidelity with residuals of governing PDEs, initial conditions, boundary conditions, and, in one Navier–Stokes example, a pressure gauge constraint (Jung et al., 9 Jun 2026). The paper contrasts standard MFA, regularized MFA, and PI-MFA, and reports that PI-MFA reduces PDE residuals, improves global balance-law consistency, produces more physically faithful reconstructions, and, for physically inconsistent data, can achieve lower approximation errors than data-only spline fits while offering computational advantages over the tested physics-informed neural networks (Jung et al., 9 Jun 2026). In the reported benchmarks, fit or training times are 14 s versus 45 s for 1D convection-diffusion, 301 s versus 502 s for 2D Burgers, and 602 s versus 2520 s for 2D Navier–Stokes when PI-MFA is compared to the tested PINN baseline (Jung et al., 9 Jun 2026).

5. Statistical and low-rank models for multivariate functional data

In functional data analysis, MFA often takes the form of a low-dimensional representation embedded inside a regression or latent-variable model. The multiFAMM framework models a multivariate functional response

$$f(\mathbf{x}) \approx \sum_{k_1=1}^{K_1}\cdots\sum_{k_n=1}^{K_n} w_{\mathbf{k}} \prod_{m=1}^n \phi^{(m)}_{k_m}(x_m), \qquad \mathbf{x}=(x_1,\dots,x_n)\in[0,1]^n.$$4

and approximates each multivariate functional random process by a truncated multivariate Karhunen–Loève expansion,

$$f(\mathbf{x}) \approx \sum_{k_1=1}^{K_1}\cdots\sum_{k_n=1}^{K_n} w_{\mathbf{k}} \prod_{m=1}^n \phi^{(m)}_{k_m}(x_m), \qquad \mathbf{x}=(x_1,\dots,x_n)\in[0,1]^n.$$5

Here the low-rank approximation is learned by multivariate functional principal component analysis, which captures both auto-correlation within each dimension and cross-correlation across dimensions (Volkmann et al., 2021). The paper argues that the multivariate approach is more parsimonious than fitting independent univariate models while maintaining or improving model fit, and it emphasizes corrected confidence bands and improved random-effects estimation. It also reports implementation in the R package multifamm (Volkmann et al., 2021).

SupFTSVD addresses high-dimensional longitudinal multivariate data, especially $$f(\mathbf{x}) \approx \sum_{k_1=1}^{K_1}\cdots\sum_{k_n=1}^{K_n} w_{\mathbf{k}} \prod_{m=1}^n \phi^{(m)}_{k_m}(x_m), \qquad \mathbf{x}=(x_1,\dots,x_n)\in[0,1]^n.$$6 settings such as microbiome studies (Alam et al., 2024). Its supervised functional tensor singular value decomposition models the data as

$$f(\mathbf{x}) \approx \sum_{k_1=1}^{K_1}\cdots\sum_{k_n=1}^{K_n} w_{\mathbf{k}} \prod_{m=1}^n \phi^{(m)}_{k_m}(x_m), \qquad \mathbf{x}=(x_1,\dots,x_n)\in[0,1]^n.$$7

This is a rank-$$f(\mathbf{x}) \approx \sum_{k_1=1}^{K_1}\cdots\sum_{k_n=1}^{K_n} w_{\mathbf{k}} \prod_{m=1}^n \phi^{(m)}_{k_m}(x_m), \qquad \mathbf{x}=(x_1,\dots,x_n)\in[0,1]^n.$$8 CP-like low-rank decomposition with continuous time, subject-level covariates, feature loadings, and time-varying singular functions (Alam et al., 2024). The paper emphasizes that, unlike classical multivariate FPCA focused mainly on time reduction, SupFTSVD provides simultaneous low-dimensional representations of subjects, features, and time, as well as subject-specific trajectories (Alam et al., 2024). In simulations it reports higher $$f(\mathbf{x}) \approx \sum_{k_1=1}^{K_1}\cdots\sum_{k_n=1}^{K_n} w_{\mathbf{k}} \prod_{m=1}^n \phi^{(m)}_{k_m}(x_m), \qquad \mathbf{x}=(x_1,\dots,x_n)\in[0,1]^n.$$9 and lower MSPE than unsupervised FTSVD when auxiliary variables are informative, and in FARMM and ECAM applications it reports biologically meaningful components associated with covariates such as diet, antibiotics, breastfeeding, and delivery mode (Alam et al., 2024).

A Bayesian variant appears in the multivariate functional dynamic linear model. There, functional observations are represented as

$d$0

with cubic B-spline factor loading curves, time-invariant bases, and optional sharing of the same basis across multiple outcomes (Kowal et al., 2014). The basis is made smooth and identifiable through roughness penalties, orthogonality constraints, and an ordering of smoothing parameters. Temporal dependence is modeled in the latent factors through a hierarchical dynamic linear model, and posterior inference is performed by Gibbs sampling (Kowal et al., 2014). The paper emphasizes that the shared-basis construction improves information sharing across multivariate observations and permits exact inference up to MCMC error (Kowal et al., 2014).

These statistical uses of MFA differ from spline-based scientific-data MFA in purpose, but they share the same structural principle: a complex multivariate functional object is projected onto a lower-dimensional family of basis functions or separable components, and the approximation is judged not only by reconstruction but also by interpretability, covariance fidelity, or predictive utility.

6. Neural, rational, operator-theoretic, and probabilistic generalizations

NeuApprox reframes MFA as a learned block-separable approximation problem (Wu et al., 4 Mar 2026). Its central claim is that replacing hand-crafted bases by neural basis functions yields both stronger approximation ability and more flexible data adaptation. The paper proves a universal approximation theorem by combining the univariate universal approximation theorem with Stone–Weierstrass, and reports experiments on multispectral images, videos, light fields, traffic data, and point clouds (Wu et al., 4 Mar 2026). It also studies pretraining, LoRA-based fine-tuning, and training from scratch, and states that pretraining alone is weaker than fine-tuning while training from scratch can be slower (Wu et al., 4 Mar 2026).

The polynomial-ANOVA framework studies high-dimensional functions with low superposition dimension (Potts et al., 2019). Its core contribution is to connect orthonormal polynomial coefficients with ANOVA terms through the support of the multi-index, so that subsets of variables and groups of polynomial frequencies become directly identifiable (Potts et al., 2019). The method reconstructs from scattered data by grouped least squares, uses global sensitivity indices to estimate which variable subsets matter, and then refines the active set of interactions before recomputing the approximation (Potts et al., 2019). The paper’s central message is that approximation becomes feasible in high dimensions when the function is dominated by low-order interactions (Potts et al., 2019).

The p-AAA algorithm provides a data-driven rational approximation for parametric dynamical systems and more general multivariate sampled functions (Rodriguez et al., 2020). It generalizes the one-variable AAA barycentric form to several variables, requires only function evaluations, and selects interpolation points greedily according to the current error (Rodriguez et al., 2020). The paper also extends the construction to matrix-valued functions through tangential compression and states a connection to tangential interpolation theory (Rodriguez et al., 2020).

The multivariate max-min neural network operator framework occupies a different approximation-theoretic niche (Yadav et al., 12 Jan 2026). It approximates $d$1 by nonlinear operators built from sigmoidal activations and max-min aggregation, proves continuity, monotonicity, max-min pseudo-linearity, subadditivity, and a Lipschitz-type perturbation estimate, and establishes both pointwise and uniform convergence (Yadav et al., 12 Jan 2026). Quantitative rates are derived through a modulus of continuity and a multivariate generalized absolute moment, and the paper argues that the resulting operators provide efficient and stable approximation tools (Yadav et al., 12 Jan 2026).

A probabilistic form of multivariate functional approximation appears in Stein’s method of exchangeable pairs (Döbler et al., 2020). There the target is a random element $d$2, and the approximation is by a suitable Gaussian process via an exchangeable-pairs coupling satisfying an approximate linear regression property (Döbler et al., 2020). The paper develops an abstract process-level Gaussian approximation theorem and applies it to joint subgraph counts in Erdős–Rényi graphs and to vectors of weighted degenerate $d$3-processes, including success runs of different lengths (Döbler et al., 2020). This use of MFA is distinct from deterministic surrogate construction, but it extends the same underlying idea of replacing a complex multivariate functional object by a more tractable functional representation.

7. Recurring trade-offs, misconceptions, and open directions

A persistent theme is the trade-off between local simplicity and global fidelity. In volume rendering, local methods such as nearest neighbor, trilinear interpolation, tricubic interpolation, and Catmull-Rom are computationally straightforward, but the cited work argues that they can create artifacts, inaccurate gradients, and incorrect reconstructions near boundaries or in sparse regions. MFA-DVR attributes its gains to a global approximation optimized over the entire domain rather than to separate local fits (Sun et al., 2022). At the same time, the same papers state that MFA query time is the primary performance bottleneck for interactive visualization, which motivates distributed rendering and compressed models (Sun et al., 2023).

Another recurring trade-off is between compactness and fidelity. In spline-based MFA and MFA-DVV, the number of control points largely determines model size, rendering time, and compression level (Sun et al., 2023). In NeuApprox, performance depends on the depth and width of the basis networks, the number of block terms $d$4, and the coefficient tensor sizes; too few parameters underfit, while too many layers or too large coefficient tensors can overfit or increase computation (Wu et al., 4 Mar 2026). In SupFTSVD, tensor rank and the smoothness of time-varying singular functions govern the balance between approximation quality and parsimony (Alam et al., 2024).

Continuity enforcement is another nontrivial issue. Parallel MFA encoding can introduce discontinuities if subdomains are fitted independently. The domain-decomposition literature argues that post-processing blending only smooths symptoms, whereas constrained minimization and Schwarz-type coupling enforce continuity in the representation itself (Mahadevan et al., 2022). A plausible implication is that continuity is not merely a cosmetic property for visualization, but an operational requirement for derivative queries and topology.

The literature also distinguishes data fidelity from structural fidelity. PI-MFA explicitly states that purely data-driven reconstructions need not preserve governing physics, especially when the input data are physically inconsistent due to low-fidelity discretizations or unmodeled discrepancies (Jung et al., 9 Jun 2026). Its results suggest that a better pointwise fit is not always the most useful reconstruction if PDE residuals, conservation laws, or hidden variables such as pressure are the downstream objective.

Topological analysis introduces a parallel distinction between direct continuous analysis and sampled discrete analysis. CPE-MFA and the later topology-extraction framework argue that discretization or resampling may create or miss critical points, loops, and other descriptors, whereas direct analysis on the continuous MFA avoids those sampling artifacts (Ma et al., 2024). Yet the same papers are explicit about current limitations: isolated rather than non-isolated critical sets, imperfect boundary handling, connection ambiguity, gaps in tracing, and dependence on step size or initial sampling (Ma et al., 11 Aug 2025).

A common misconception is that MFA names a single representation. The cited literature does not support that reading. It includes spline-based global surrogates for scientific volumes, supervised tensor decompositions for longitudinal functional data, shared spline factor models for Bayesian dynamic inference, polynomial-ANOVA expansions for scattered high-dimensional data, adaptive barycentric rational surrogates for parametric systems, neural basis expansions with universal approximation guarantees, max-min neural operators with convergence theorems, and Gaussian-process approximations for random multivariate functions (Kowal et al., 2014). This suggests that the unity of the field lies less in a single basis than in a common strategy: replace a complex multivariate functional object by a lower-complexity functional representation that preserves the properties needed for the task at hand.

Several open directions are stated explicitly in the cited papers. The topology papers mention Morse or Morse–Smale complexes, 3D domains, and better guarantees for initial-point selection and connectivity (Ma et al., 11 Aug 2025). The critical-point paper mentions Bezier clipping and interval Newton methods as future alternatives to Newton’s method (Ma et al., 2024). The parallel MFA paper notes future interest in T-splines for non-matching knot locations (Mahadevan et al., 2022). The VTK integration paper notes potential use as a VTK plugin (Sun et al., 2022). The PI-MFA paper identifies structured Gauss–Newton-type methods as an interesting future direction (Jung et al., 9 Jun 2026). Taken together, these directions indicate that MFA is evolving simultaneously as an approximation theory, a differentiable surrogate-modeling paradigm, and a computational substrate for visualization, topology, and physics-aware analysis.