Matrix solutions and spatial tensors

Mathematical range and ordering

P0, P0g, P1 and H1 accept Math::SpatialMatrix<Scalar> ranges for Real and Complex entries. Rows and columns are explicit, independently between one and RODIN_MAXIMAL_SPACE_DIMENSION. Symmetry is not imposed. A matrix space is the tensor product of its scalar space with the matrix units: $ V_h^{m\times n}=V_h\otimes\mathbb{K}^{m\times n} $ . Matrix specializations live beside the scalar and vector specializations in each space family. Reference bases, physical mappings and DOF distribution belong to that family. The componentwise construction here describes P0, P0g, P1 and H1; future $H(\mathrm{div})$ and $H(\mathrm{curl})$ families must supply their own conforming bases and Piola mappings. Local and global coefficient indices are interleaved in row-major order: $ a(mn)+in+j $ . No Voigt scaling or symmetry reduction is applied. This coefficient numbering is independent of matrix memory storage; Matrix storage and finite-element component ordering describes both conventions and their consequences for assembly.

H1 stressSpace(std::integral_constant<size_t, 3>{}, mesh, 3, 3);
GridFunction stress(stressSpace);
TrialFunction sigma(stressSpace);
TestFunction tau(stressSpace);
BilinearForm mass(sigma, tau);
mass = Integral(sigma, tau);

Differential and constitutive operators

Inner products use $ A:B=\sum_{ij}A_{ij}\overline{B_{ij}} $ . Grad and Jacobian return SpatialTensor<Scalar,3> with $ G_{ijk}=\partial_k A_{ij} $ . Derivative selects its last axis. Div uses row-wise divergence $ (\operatorname{div}A)_i=\sum_j\partial_jA_{ij} $ and requires the column count to match the ambient spatial dimension. A matrix coefficient multiplying a matrix solution acts on its rows; a fourth-order tensor coefficient couples arbitrary matrix entries: $ (\mathcal C A)_{ij}=\sum_{kl}\mathcal C_{ijkl}A_{kl} $ . TensorFunction wraps constant or callable tensor coefficients. Tensor rank is a compile-time parameter; extents are runtime values with stack storage. Tensor-valued solution spaces are not introduced by this matrix-range API. A linear form must have a scalar integrand. Bare Integral(v) is therefore rejected for vector- and matrix-valued test functions; Integral(F,v) contracts a matching coefficient with the test function to form a scalar.

Matrix Jacobians have rank three. The vector-only P1 Jacobian handlers and the weighted, rank-one and advection H1 handlers are restricted to their supported ranges. Matrix expressions use the existing generic quadrature rule for these forms; the componentwise unweighted H1 Jacobian handler supports matrix ranges directly. An explicit quadrature order remains appropriate for rational pyramid bases and curved mappings.

Nonlocal P1 forms

The existing P1 potential quadrature handler accepts matrix trial and test functions. A kernel returning a scalar multiplies every component. A kernel filling a SpatialMatrix acts on rows, while a kernel filling a SpatialTensor<Scalar,4> couples all entries:

\[ (K A)_{ij}(x)=\int_\Omega \sum_{kl}\mathcal K_{ijkl}(x,y)A_{kl}(y)\,dy. \]

The callback signatures are kernel(x,y), kernel(matrix,x,y) and kernel(tensor,x,y), respectively. Output extents must match the field dimensions. The existing collapsed quadrature on coincident triangles and centroid rules on other panel pairs are retained; nonsingular constant kernels do not establish accuracy for arbitrary singular kernels. This handler is specific to P1; nonlocal handlers for the other space families are not introduced here. The nonlocal regressions cover local Eigen and OpenMP assembly; distributed nonlocal assembly is not validated by the Jacobian and reaction-diffusion MPI tests.

Math::SpatialTensor<Real, 4> elasticity(3, 3, 3, 3);
elasticity.setZero();
// Fill the constitutive entries, then assemble without a Voigt conversion.
auto constitutive = Integral(TensorFunction(elasticity) * sigma, tau);

Callable MatrixFunction coefficients accept explicit dimensions and setOrder for polynomial quadrature selection. Integrator setOrder specifies the integration degree. Rational pyramid bases and curved transformations require quadrature chosen for the mapped integrand; polynomial reference-degree defaults alone do not guarantee exactness. Grid-function derivatives on interfaces require a traceOf domain, as for scalar fields. Shape-function derivatives on a boundary are derivatives of the boundary element, matching the scalar finite element semantics.

Backend and format support

The same range and DOF maps are consumed by Eigen dense/sparse sequential and OpenMP assembly and by local/distributed PETSc assembly. MPI spaces expand scalar ownership and ghost maps component-wise. They cache collective global sizes at construction so local evaluations do not enter collectives.

HDF5 persistence preserves raw coefficients, matrix rows/columns and complex imaginary parts. PETSc persistence stores the owned portion per shard. XDMF writes full matrix-valued arrays, including nonsymmetric entries. MFEM uses the scalar family's node conversion for each matrix entry; matrix dimensions must be supplied by the destination space because MFEM records only component count. Its existing unsupported collections remain unsupported (in particular the global constant P0g collection). MEDIT uses separate scalar vertex fields for every entry. Vertex formats interpolate higher-order fields from vertices and cannot preserve their arbitrary interior coefficients; use HDF5 for exact persistence.

Validation coverage

Native tests exercise all eight reference geometries, every matrix shape between one and three axes, nodal duality through cubic H1 and additional fourth/fifth-order representatives. Assembly and differential tests cover every non-point cell geometry, P0/P0g/P1/H1, rectangular and square ranges, complex entries, fourth-order material coupling, physical finite differences, curved and embedded triangles, boundary traces, cache lifetime, and Dirichlet data. Manufactured reaction-diffusion tests exercise complete solve paths. MPI/PETSc tests exercise ownership, ghosts, assembly and solves at one through four ranks. IO tests check real/complex HDF5, MEDIT vertex fields, MFEM node permutations and XDMF matrix metadata. This is finite regression coverage, rather than an assertion that arbitrary order, geometry and coefficient combinations have all been enumerated.