Implementing Basis Module

[zasexton]

Problem

In the solver, basis functions are hard-coded through a Fortran-style set of per-element dispatch tables, cached directly inside mesh and function-space data structures, and then copied, recomputed, or patched by special-case logic whenever we need mixed spaces or NURBS behavior. That makes basis evaluation tightly coupled to mesh storage, quadrature setup, and physics assembly. Adding new basis families or extending derivative support tends to spread changes across unrelated parts of the solver. A modular Basis folder therefore is necessary because it gives us a single mesh-agnostic interface for basis families, plus centralized construction, caching, and evaluation machinery, so different basis families (Lagrange, Bernstein, B-spline/NURBS, and vector bases) can be implemented once, tested once, and reused consistently across elements, geometry mappings, and assembly kernels instead of proliferating element-specific and physics-specific code paths.

• Basis evaluation is implemented as a monolithic Fortran-port dispatch layer. Code/Source/solver/nn.cpp pulls in large element tables, and Code/Source/solver/nn_elem_gnn.h hard-codes closed-form N/Nx formulas per element type instead of exposing a reusable basis abstraction.
• Basis state is stored directly on mesh and function-space containers rather than encapsulated in its own module. Code/Source/solver/ComMod.h and Code/Source/solver/ComMod.h embed N, Nx, Nxx, bounds, and quadrature data into fsType and mshType, while Code/Source/solver/nn.cpp and Code/Source/solver/fs.cpp allocate and populate those arrays procedurally.
• Mixed spaces and extensibility are handled by special cases. Code/Source/solver/fs.cpp and Code/Source/solver/fs.cpp build Taylor-Hood spaces by copying or regenerating alternate shape tables, while NURBS support is still scattered as repeated update hooks or TODOs in places like Code/Source/solver/fluid.cpp, Code/Source/solver/eq_assem.cpp, and Code/Source/solver/nn.cpp.

Solution

To address the current solver’s basis-function, we should move all reference-element basis logic into a dedicated, mesh-agnostic Basis module with a common interface for evaluating values, gradients, and higher derivatives, and create those bases through a factory keyed by element family, order, continuity, and field type. Instead of storing handwritten N, Nx, and Nxx tables directly on mshType and fsType and rebuilding them through scattered element-specific branches, the mesh and element layers should hold lightweight basis and quadrature objects, while geometry mapping handles only the transformation from reference to physical space. Mixed discretizations such as Taylor-Hood should be represented as explicit mixed/composite elements rather than ad hoc control flow, and any future implemented NURBS/IGA support should be implemented as another basis family rather than a repeated special case. With caching and batched evaluation centralized in the basis layer, the same implementations can be reused consistently by physics assembly kernels, which reduces duplication, makes new basis families easier to add, and localizes both testing and performance optimization to a single part of the codebase.

This issue will outline Basis module design and API/ABI infrastructure.

[zasexton]

Currently, the solver evaluates shape functions through a set of hand-coded dispatch maps in nn_elem_gnn.h, nn_elem_gip.h, and nn_elem_gnnxx.h — per-element-type lambdas that compute basis function values, reference gradients, second derivatives, and Gauss quadrature points. Each supported element type (Tet4, Hex8, Tri6, etc.) has its own lambda with the shape function formulas written out explicitly. These maps are called through the nn::get_gnn and nn::get_gip free functions, which populate fsType arrays (N, Nx, Nxx, w, xi) at setup time. The physics assembly files then read from these arrays during element loops, with nn::gnn() handling the Jacobian transform from reference to physical gradients.

We want to replace the dispatch maps with calls through an FE Basis module's BasisFunction / BasisFactory / BasisCache API (this is probably the cleanest and most impactful route though it will require some extra data copying at setup time). The nn::get_gnn and nn::get_gip implementations will be rewritten to create a BasisFunction via a factory, evaluate it at the requested reference coordinates, and copy the results into the same fsType arrays that the current solver already expects. The Jacobian transform (nn::gnn), the fsType struct, and all physics files will not require any changes — they will continue to read fs[k].N(a,g) and fs[k].Nx(d,a,g) exactly as before.

The two prerequisites are an ElementType enum converter (the current solver uses a different enums for the same element shapes) and a quadrature point-count-to-polynomial-order mapping; the current code selects quadrature by number of Gauss points, but a more robust strategy for Basis would select by integration order. Specifying by point count is ambiguous. 4 points on a tetrahedron implies a specific rule, but 4 points on a triangle could be multiple different rules with different polynomial exactness. It also couples the caller to implementation details. For example, if someone wants to switch from a Gauss-Legendre rule to a more efficient rule with fewer points at the same accuracy, every call site that hardcodes nG has to change. The current code works around this by having the nn_elem_gip.h lambdas implicitly encode the order-to-point-count relationship per element type, but that's just a manual version of what an eventual QuadratureFactory would do automatically.

Both of these prerequisites should be small translations. Once these are in place the solver automatically gains access to every basis family the module supports — higher-order Lagrange, Serendipity, etc — without touching any physics code or other infrastructural code.

[alisonmarsden]

@mrp089 it would be great to get your input on this - it's the first of several issues @zasexton and @ktbolt are jointly working on related to overall code refactoring plans. I'm also wondering if we should consider how interfacing with MFEM or deal.II might be worked into our plan, either now or in the future. (FEBio is migrating to MFEM now and LifeX uses deal.II)

[alisonmarsden]

@zasexton can we generate a pull request on this so others can review?

[zasexton]

Yes, still waiting on one more review for PR draft #537 for the exception infrastructure that we will use within the basis module

[alisonmarsden]

@zasexton are we ready to go with this one now?

[ktbolt]

@zasexton The current nn*.h code is about 3000 lines of code; the svMultiPhysics/Code/Source/solver/FE/Basis code is quite a bit larger with the 20,000 lines of code. Let's try to reduce the size of that code by removing functionality that is not relevant to our current requirements.