Issue 525 FE Basis refactor implementation#561
Conversation
updating fork repo with upstream head
update main to upstream branch
… unit tests for refactored basis functions
|
From the CI test cases we got the following:
The affected meshes are |
|
In the FSI fluid assembly,
That means:
Therefore the Hessian is not just participating in the tangent assembly and may not necessarily converge on the same solution as the current reference vtu. This is also why we do not see such CI failures for test cases that use Tet elements because |
|
Overall, we see that providing the non-zero Hessian values slightly increases the number of linear iterations required but does not impact the number of nonlinear iterations for the FSI test cases considered. More importantly, this difference in how the Hessian is used in the formulated residual suggested that the problem being solved is slightly different depending on if a user is providing brick elements versus 
|
|
Confirmed. The difference in accuracy for the FSI test cases when employing the new FE Basis functions is due to the non-zero Hessian values for the reference elements. We can see that when we zero-out the values for the Hessian on the 
|
|
@zasexton Were Hessians also zeroed out in the Fortran svFSI? I think we should update reference solutions unless there is a good reason to use zero Hessians. @sujaldave, @hanzhao2020 Is there any reason (in VMS or something) to force |
|
@aabrown100-git yes, the Hessian was also effectively zero for the |
Regenerate affected FSI and FSI-ustruct HEX8 result_005.vtu references for the FE Basis path with nonzero HEX8 Hessian contributions. Update the pipe_3d PETSc and Trilinos references to match the base pipe_3d reference, preserving the existing shared-reference pattern across linear algebra variants.
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #561 +/- ##
==========================================
+ Coverage 69.44% 72.98% +3.54%
==========================================
Files 181 237 +56
Lines 34072 37670 +3598
Branches 5930 6539 +609
==========================================
+ Hits 23662 27495 +3833
+ Misses 10273 9940 -333
- Partials 137 235 +98 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
… reference tables and related unit tests
…from the single ReferenceNodeLayout generator
…r bases while preserving named-element request behavior and covering both paths with FE basis tests
…he shared tensor-Lagrange node distribution) well-conditioned at high order
ReferenceNodeLayout::serendipity_node_coords owner
michelebucelli
left a comment
michelebucelli
left a comment
There was a problem hiding this comment.
Thank you @zasexton, this is lots of work! I added some other comments.
Also good job with the documentation, it's very detailed.
| @@ -32,10 +34,27 @@ struct BasisRequest { | |||
| std::vector<std::vector<double>> axis_weights{}; | |||
| std::vector<int> tensor_extents{}; | |||
| std::string custom_id{}; | |||
| // Reference topology for arbitrary-order bases. This field is intentionally | |||
| // last so existing aggregate initializers for named elements keep their | |||
| // positional meaning. | |||
There was a problem hiding this comment.
I suggest turning these comments into Doxygen comments. It would also be nice to add Doxygen comments to the other members as well.
(Maybe the comment on topology being intentionally last can be omitted from Doxygen, since it's more of a low-level implementation detail that shouldn't affect the users of this type).
| namespace svmp { | ||
| namespace FE { | ||
| namespace basis { |
There was a problem hiding this comment.
Minor: I find the following more compact syntax a bit better (less brackets to keep track of):
| namespace svmp { | |
| namespace FE { | |
| namespace basis { | |
| namespace svmp::FE::basis { |
| /** | ||
| * @brief Invalid Basis request or configuration | ||
| */ | ||
| class BasisConfigurationException : public BasisException { | ||
| public: | ||
| BasisConfigurationException(const std::string& message, | ||
| const char* file, | ||
| int line, | ||
| const char* function) | ||
| : BasisException(message, file, line, function, StatusCode::InvalidArgument) {} | ||
| }; | ||
|
|
||
| /** | ||
| * @brief Requested element topology is incompatible with the basis family | ||
| */ | ||
| class BasisElementCompatibilityException : public BasisException { | ||
| public: | ||
| BasisElementCompatibilityException(const std::string& message, | ||
| const char* file, | ||
| int line, | ||
| const char* function) | ||
| : BasisException(message, file, line, function, StatusCode::InvalidArgument) {} | ||
| }; | ||
|
|
||
| /** | ||
| * @brief Basis evaluation request cannot be satisfied | ||
| */ | ||
| class BasisEvaluationException : public BasisException { | ||
| public: | ||
| BasisEvaluationException(const std::string& message, | ||
| const char* file, | ||
| int line, | ||
| const char* function) | ||
| : BasisException(message, file, line, function, StatusCode::InvalidArgument) {} | ||
| }; | ||
|
|
||
| /** | ||
| * @brief Public-to-canonical node ordering or coordinate lookup failure | ||
| */ | ||
| class BasisNodeOrderingException : public BasisException { | ||
| public: | ||
| BasisNodeOrderingException(const std::string& message, | ||
| const char* file, | ||
| int line, | ||
| const char* function) | ||
| : BasisException(message, file, line, function, StatusCode::InvalidArgument) {} | ||
| }; | ||
|
|
||
| /** | ||
| * @brief Internal basis construction or transform setup failure | ||
| */ | ||
| class BasisConstructionException : public BasisException { | ||
| public: | ||
| BasisConstructionException(const std::string& message, | ||
| const char* file, | ||
| int line, | ||
| const char* function) | ||
| : BasisException(message, file, line, function, StatusCode::InternalError) {} | ||
| }; |
There was a problem hiding this comment.
In the same vein as the previous comment, I am not sure I fully understand the use cases for these classes. While I imagine that the typical developer wouldn't have to raise these exceptions, I still think it would be nice to be a bit more pedagogic for future reference.
I think that adding a more detailed explanatory sentence to the documentation of these classes would help. Additionally, it would be nice to provide one clarifying example for each (e.g. for BasisElementCompatibilityException, you could write that it's thrown if requesting a tetrahedral type of element on a hexahedral grid, or something in that vein).
| void LagrangeBasis::evaluate_values(const Vec3& xi, | ||
| std::vector<double>& values) const { | ||
| values.resize(size()); | ||
| evaluate_values_to(xi, std::span<double>(values.data(), values.size())); | ||
| } |
There was a problem hiding this comment.
Related to my previous comment on BasisFunction::evaluate_values, I think this could be the default implementation in BasisFunction (perhaps even marked as final).
The same suggestion applies to the other evaluate_ methods.
| #if defined(SVMP_FE_WITH_MESH) && SVMP_FE_WITH_MESH | ||
| # include "Mesh/Core/MeshTypes.h" | ||
| /** Nonzero when FE shares scalar/index types with the Mesh library. */ | ||
| # define SVMP_FE_HAS_MESH_TYPES 1 | ||
| #else | ||
| // Build FE without Mesh types unless explicitly enabled. | ||
| /** Nonzero when FE shares scalar/index types with the Mesh library. */ | ||
| # define SVMP_FE_HAS_MESH_TYPES 0 | ||
| #endif |
There was a problem hiding this comment.
Is this Mesh library an external (optional) module of svmp? Or is it part of the planned refactor?
There was a problem hiding this comment.
In my opinion, it will be part of the planned refactor as another library that will be used by the FE library. The trait interface will allow for topology and cell identification and mapping to the corresponding FE elements. Otherwise how mesh data is processed and augmented I imagine we can detangle from much of the FE code. These traits and cmake definitions are only the precursors. I can also remove them here to make the PR more narrow
| // Linear elements | ||
| Line2 = 0, ///< 2-node line | ||
| Triangle3 = 1, ///< 3-node triangle | ||
| Quad4 = 2, ///< 4-node quadrilateral | ||
| Tetra4 = 3, ///< 4-node tetrahedron | ||
| Hex8 = 4, ///< 8-node hexahedron | ||
| Wedge6 = 5, ///< 6-node wedge/prism | ||
| Pyramid5 = 6, ///< 5-node pyramid | ||
|
|
||
| // Quadratic elements | ||
| Line3 = 10, ///< 3-node line | ||
| Triangle6 = 11, ///< 6-node triangle | ||
| Quad9 = 12, ///< 9-node quadrilateral (bi-quadratic) | ||
| Quad8 = 13, ///< 8-node quadrilateral (serendipity) | ||
| Tetra10 = 14, ///< 10-node tetrahedron | ||
| Hex27 = 15, ///< 27-node hexahedron (tri-quadratic) | ||
| Hex20 = 16, ///< 20-node hexahedron (serendipity) | ||
| Wedge15 = 17, ///< 15-node wedge | ||
| Wedge18 = 18, ///< 18-node wedge (complete quadratic) | ||
| Pyramid13 = 19, ///< 13-node pyramid | ||
| Pyramid14 = 20, ///< 14-node pyramid | ||
|
|
||
| // Special elements | ||
| Point1 = 30, ///< 1-node point element | ||
|
|
||
| Unknown = 255 ///< Unrecognized or uninitialized element type |
There was a problem hiding this comment.
Is there any meaning attached to the numbers used here? If not, I would suggest removing them (after all, the point of enum class is that the underlying numbers are unimportant as long as they are different).
| /** | ||
| * @brief DOF-specific index type | ||
| * | ||
| * Strong type alias to prevent mixing DOF indices with other indices. | ||
| * Provides type safety at compile time. | ||
| */ | ||
| struct DofIndex { | ||
| GlobalIndex value; ///< Underlying global DOF index; negative values are invalid. | ||
|
|
||
| /** | ||
| * @brief Construct a DOF index, defaulting to the invalid sentinel. | ||
| * @param v Global DOF index value. | ||
| */ | ||
| constexpr explicit DofIndex(GlobalIndex v = -1) noexcept : value(v) {} | ||
| /** | ||
| * @brief Convert to the underlying global index value. | ||
| * @return The stored global index. | ||
| */ | ||
| constexpr operator GlobalIndex() const noexcept { return value; } | ||
| /** | ||
| * @brief Check whether this index refers to a valid DOF. | ||
| * @return True when the stored value is non-negative. | ||
| */ | ||
| constexpr bool is_valid() const noexcept { return value >= 0; } | ||
| }; |
There was a problem hiding this comment.
Would it make sense for this to be replaced by StrongType<GlobalIndex, ...>?
There was a problem hiding this comment.
I agree. I'll make this change.
| /// (basis_factory::default_basis_request). The switch deliberately has no | ||
| /// default case so that compilers building with -Wswitch flag any newly added | ||
| /// solver element type that is missing a mapping here. Returns std::nullopt for |
There was a problem hiding this comment.
I would have a mild preference towards a default case throwing an exception, rather than relying on the warning being turned on (although I recognize this would be caught at runtime, so the behavior wouldn't be exactly the same).
- include the basis family in the documented basis identity (topology, order, family) - replace the order() normalization wording with a plain description and example - add cross-reference links to the module entities in the group documentation - add topology and family examples to the basis object contract - refer to shape function values in the BasisFunction class description
Add an explanatory sentence and a concrete example to each derived basis exception, drawn from its real throw sites, so the specific types guide developers on when to raise or expect them.
…ation history Name the orders and elements that use equispaced nodes directly, instead of referring to what previous element layouts kept, so the note stays accurate as the implementation evolves.
create() and create_default_for() now return std::unique_ptr<BasisFunction>, which has simpler ownership semantics and is cheaper than shared_ptr. The single caching consumer keeps a shared_ptr cache; the returned unique_ptr converts to it on insertion. All FE/Basis unit tests pass.
…erloads in the base The span *_to evaluators are now the primitives a concrete basis implements: evaluate_values_to is pure virtual, evaluate_gradients_to/evaluate_hessians_to default to reporting not-implemented, and a protected evaluate_all_to provides the single-pass combined evaluation. The std::vector overloads (evaluate_values, evaluate_gradients, evaluate_hessians, evaluate_all) are implemented once on the base class, so LagrangeBasis and SerendipityBasis no longer duplicate those wrappers. This inverts the previous arrangement where the vector form was the primitive and the base supplied an allocate-and-copy span fallback. A minimal basis now overrides the span primitive instead of the vector form. The BasisFunction test helpers are updated accordingly, and the former fallback test is rewritten to verify the base vector overloads forward to the span primitives. All FE/Basis unit tests pass.
…tation - note that the vector evaluators use an output argument so callers can reuse a buffer across evaluations rather than allocating per call - document that the BasisTraits classifiers return -1 / Unknown as sentinels that callers validate into exceptions (the constexpr noexcept helpers cannot throw) - keep the BasisRequest::topology field-ordering note out of the rendered docs
…roup Replace the @cond INTERNAL exclusion with a documented internal group whose @warning states the declarations are internal (use basis_factory and BasisFunction::nodes() instead) and may change. Core developers now get the rendered documentation while model-level callers are clearly steered away.
… Vector Document that FE/Math/Vector.h is a fixed-size, compile-time-length element-level vector in svmp::FE::math, distinct from and not a replacement for the legacy dynamic global ::Vector in solver/Vector.h.
…nd ElementType bands - explain that the Mesh library is an optional external module and FE builds standalone with fallback types when it is absent - record why GlobalIndex stays a plain alias (raw PETSc/Trilinos interop) and that DofIndex is the strong wrapper - document that the explicit ElementType values are intentional bands with a uint8_t Unknown sentinel
…e_label Add Doxygen to the dense linear-algebra functions and the DenseLUSolver type, and rename the label argument (and the DenseLUSolver member) to error_message_label so its role as an error-message prefix is self-documenting. All FE math unit tests pass.
Add a default case that raises BasisElementCompatibilityException for a solver element type with no FE mapping, so a missing mapping fails loudly instead of relying on the unhandled-enum compiler warning. The deliberate NA/PNT/NRB cases still return std::nullopt.
MathJax 2 was pinned as a default when the FE documentation was added and is now end-of-life. Move to the maintained MathJax 3 (chtml output, mathjax@3 CDN, ams extension); doxygen emits the v3 bootstrap and the AMS math used in the FE docs is supported.
…mbers Generating HTML (not just parsing) surfaced reference/doc issues the no-output checks missed: - qualify the @ref class/struct/enum links in the Basis group doc so they resolve - use plain text for the basis_factory namespace mention (no resolvable @ref) - avoid the ::Vector explicit-link request in the Math vector doc - add @return tags and document the NodeOrderingConventions entities that became visible when the @cond INTERNAL exclusion was removed Net FE Doxygen warnings drop from 47 to 19 (the rest pre-existing).
…assertion helpers - add svmp::NotImplementedException and svmp::IndexOutOfRangeException (CoreException-derived) - default the ExceptionT template parameter of not_implemented() and check_index() to those types, so they can be used without naming an exception while existing explicit callers are unaffected - document the helpers, including that check vs check_arg are the same check named for intent and that throw_if is the logical inverse of check, and when to use each not_implemented overload Note: svmp::FE::NotImplementedException stays for FE code needing FEException ancestry. All FE/exception-helper unit tests pass.
…ption infrastructure
4 participants
Current situation
Tracks #525.
The solver currently relies on legacy table-driven shape-function paths in
nn.cpp, which makes basis evaluation, Hessian support, node-ordering validation, and parity testing difficult to extend. This PR introduces a self-contained FE Basis layer while preserving the existing solver-facing storage contracts.Release Notes
nn::get_gnnandnn::get_gn_nxxvolume/face evaluations through the new FE Basis adapter.SVMP_BASIS_MODE=auto|legacy|feto compare legacy and FE Basis evaluation paths at runtime.Documentation
SVMP_BASIS_MODE.Testing
tests/unitTests/FE/Basisfor Lagrange/Serendipity basis evaluation, Hessians, cache/factory behavior, error paths, node ordering, solver adapter parity, and supported mapped element coverage.tests/unitTests/FE/Mathfor matrix/vector operations, expression helpers, math constants, and dense linear algebra.Code of Conduct & Contributing Guidelines