diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt
index 5e50e9834..aa5eb47f5 100644
--- a/Code/Source/solver/CMakeLists.txt
+++ b/Code/Source/solver/CMakeLists.txt
@@ -177,15 +177,15 @@ if(ENABLE_COVERAGE)
 endif()
 
 add_subdirectory(Core)
+add_subdirectory(common/math)
 add_subdirectory(FE/Common)
-add_subdirectory(FE/Math)
 add_subdirectory(FE/Quadrature)
 add_subdirectory(FE/Basis)
 
 set(SVMP_ARCHITECTURE_TARGETS
   svmp_core
   svmp_fe_common
-  svmp_fe_math
+  svmp_math
   svmp_fe_basis
   svmp_fe_quadrature
 )
diff --git a/Code/Source/solver/FE/Basis/BasisFunction.h b/Code/Source/solver/FE/Basis/BasisFunction.h
index b8e52778d..6fc6198b9 100644
--- a/Code/Source/solver/FE/Basis/BasisFunction.h
+++ b/Code/Source/solver/FE/Basis/BasisFunction.h
@@ -6,8 +6,8 @@
 
 #include "BasisExceptions.h"
 #include "BasisTraits.h"
-#include "Math/Matrix.h"
-#include "Math/Vector.h"
+#include "math/Matrix.h"
+#include "math/Vector.h"
 #include "Types.h"
 
 #include <cstddef>
diff --git a/Code/Source/solver/FE/Basis/CMakeLists.txt b/Code/Source/solver/FE/Basis/CMakeLists.txt
index 3d0ee337f..6678ff91d 100644
--- a/Code/Source/solver/FE/Basis/CMakeLists.txt
+++ b/Code/Source/solver/FE/Basis/CMakeLists.txt
@@ -19,6 +19,6 @@ target_include_directories(svmp_fe_basis PUBLIC
 
 target_link_libraries(svmp_fe_basis PUBLIC
   svmp_fe_common
-  svmp_fe_math
+  svmp_math
   svmp_fe_quadrature
 )
diff --git a/Code/Source/solver/FE/Basis/NodeOrderingConventions.h b/Code/Source/solver/FE/Basis/NodeOrderingConventions.h
index 201e3e3c4..0f857727e 100644
--- a/Code/Source/solver/FE/Basis/NodeOrderingConventions.h
+++ b/Code/Source/solver/FE/Basis/NodeOrderingConventions.h
@@ -5,7 +5,7 @@
 #define SVMP_FE_BASIS_NODEORDERINGCONVENTIONS_H
 
 #include "BasisTraits.h"
-#include "Math/Vector.h"
+#include "math/Vector.h"
 #include "Types.h"
 
 #include <array>
diff --git a/Code/Source/solver/FE/Basis/SerendipityBasis.cpp b/Code/Source/solver/FE/Basis/SerendipityBasis.cpp
index bf9027136..35ab72bd5 100644
--- a/Code/Source/solver/FE/Basis/SerendipityBasis.cpp
+++ b/Code/Source/solver/FE/Basis/SerendipityBasis.cpp
@@ -3,7 +3,8 @@
 
 #include "SerendipityBasis.h"
 #include "NodeOrderingConventions.h"
-#include "Math/DenseLinearAlgebra.h"
+#include "Core/Exception.h"
+#include "math/DenseLinearAlgebra.h"
 
 #include <algorithm>
 #include <array>
@@ -195,21 +196,21 @@ std::vector<double> build_inverse_vandermonde(
     // too.
     const double norm_v = matrix_norm_inf(vandermonde, n);
 
-    // invert_dense_matrix raises a generic FEException if the Vandermonde is
+    // invert_dense_matrix raises a generic CoreException if the Vandermonde is
     // exactly singular (a rank-deficient pivot). For a serendipity family that
     // means the node set is not unisolvent at this order -- a construction failure
     // in basis terms -- so translate it to BasisConstructionException, presenting
     // the singular and the ill-conditioned cases (below) as one catchable type in
     // one vocabulary. The matrix was just built n-by-n from n modes, so a
-    // size-mismatch FEException cannot originate here; rank deficiency is the only
-    // FEException this call can raise. (Defensive: the supported node sets are
+    // size-mismatch CoreException cannot originate here; rank deficiency is the only
+    // CoreException this call can raise. (Defensive: the supported node sets are
     // provably unisolvent, so this branch is not reachable for the shipped orders.)
     std::vector<double> inverse;
     try {
         inverse = math::invert_dense_matrix(
             std::move(vandermonde), n,
             "SerendipityBasis interpolation matrix for " + label);
-    } catch (const FEException&) {
+    } catch (const CoreException&) {
         svmp::raise<BasisConstructionException>("SerendipityBasis: " + label +
                 " interpolation matrix is singular; the serendipity node set is not "
                 "unisolvent at the requested order");
diff --git a/Code/Source/solver/FE/FE.h b/Code/Source/solver/FE/FE.h
index c51c4a787..852ef25eb 100644
--- a/Code/Source/solver/FE/FE.h
+++ b/Code/Source/solver/FE/FE.h
@@ -11,7 +11,7 @@
  * This header intentionally contains no declarations. It gives Doxygen a
  * header-based home for the top-level FE group; submodule groups attach to it
  * from their own headers, including FE_Basis (Basis/BasisFunction.h),
- * FE_Common (Common/Types.h), FE_Math (Math/Vector.h), and FE_Quadrature
+ * FE_Common (Common/Types.h), shared Math (math/Vector.h), and FE_Quadrature
  * (Quadrature/QuadratureRule.h).
  */
 
diff --git a/Code/Source/solver/FE/Math/CMakeLists.txt b/Code/Source/solver/FE/Math/CMakeLists.txt
deleted file mode 100644
index c3a1c0d3b..000000000
--- a/Code/Source/solver/FE/Math/CMakeLists.txt
+++ /dev/null
@@ -1,23 +0,0 @@
-file(GLOB SVMP_FE_MATH_FILES CONFIGURE_DEPENDS
-  "${CMAKE_CURRENT_SOURCE_DIR}/*.cpp"
-  "${CMAKE_CURRENT_SOURCE_DIR}/*.h"
-  "${CMAKE_CURRENT_SOURCE_DIR}/*.hpp"
-  "${CMAKE_CURRENT_SOURCE_DIR}/*.inl"
-)
-
-add_library(svmp_fe_math OBJECT
-  ${SVMP_FE_MATH_FILES}
-)
-
-target_compile_features(svmp_fe_math PUBLIC cxx_std_20)
-
-target_include_directories(svmp_fe_math PUBLIC
-  ${CMAKE_CURRENT_SOURCE_DIR}
-  ${CMAKE_CURRENT_SOURCE_DIR}/../..
-  ${SV_SOURCE_DIR}/ThirdParty/eigen/include
-  ${SV_SOURCE_DIR}/ThirdParty/eigen/include/eigen3
-)
-
-target_link_libraries(svmp_fe_math PUBLIC
-  svmp_fe_common
-)
diff --git a/Code/Source/solver/FE/Quadrature/CMakeLists.txt b/Code/Source/solver/FE/Quadrature/CMakeLists.txt
index 92f1c9b2d..e8b434446 100644
--- a/Code/Source/solver/FE/Quadrature/CMakeLists.txt
+++ b/Code/Source/solver/FE/Quadrature/CMakeLists.txt
@@ -18,5 +18,5 @@ target_include_directories(svmp_fe_quadrature PUBLIC
 
 target_link_libraries(svmp_fe_quadrature PUBLIC
   svmp_fe_common
-  svmp_fe_math
+  svmp_math
 )
diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.h b/Code/Source/solver/FE/Quadrature/QuadratureRule.h
index c1a8e0b7b..a066279f4 100644
--- a/Code/Source/solver/FE/Quadrature/QuadratureRule.h
+++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.h
@@ -40,7 +40,7 @@
  */
 
 #include "FE/Common/Types.h"
-#include "FE/Math/Vector.h"
+#include "math/Vector.h"
 
 #include <cstddef>
 #include <span>
diff --git a/Code/Source/solver/common/math/CMakeLists.txt b/Code/Source/solver/common/math/CMakeLists.txt
new file mode 100644
index 000000000..bc8121f72
--- /dev/null
+++ b/Code/Source/solver/common/math/CMakeLists.txt
@@ -0,0 +1,22 @@
+file(GLOB_RECURSE SVMP_MATH_FILES CONFIGURE_DEPENDS
+  "${CMAKE_CURRENT_SOURCE_DIR}/include/*.h"
+  "${CMAKE_CURRENT_SOURCE_DIR}/include/*.hpp"
+  "${CMAKE_CURRENT_SOURCE_DIR}/include/*.inl"
+  "${CMAKE_CURRENT_SOURCE_DIR}/src/*.cpp"
+)
+
+add_library(svmp_math OBJECT
+  ${SVMP_MATH_FILES}
+)
+
+target_compile_features(svmp_math PUBLIC cxx_std_20)
+
+target_include_directories(svmp_math PUBLIC
+  ${CMAKE_CURRENT_SOURCE_DIR}/include
+  ${SV_SOURCE_DIR}/ThirdParty/eigen/include
+  ${SV_SOURCE_DIR}/ThirdParty/eigen/include/eigen3
+)
+
+target_link_libraries(svmp_math PUBLIC
+  svmp_core
+)
diff --git a/Code/Source/solver/FE/Math/DenseLinearAlgebra.h b/Code/Source/solver/common/math/include/math/DenseLinearAlgebra.h
similarity index 89%
rename from Code/Source/solver/FE/Math/DenseLinearAlgebra.h
rename to Code/Source/solver/common/math/include/math/DenseLinearAlgebra.h
index 5c7d6cc6e..b54dab5e5 100644
--- a/Code/Source/solver/FE/Math/DenseLinearAlgebra.h
+++ b/Code/Source/solver/common/math/include/math/DenseLinearAlgebra.h
@@ -1,10 +1,8 @@
 // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others.
 // SPDX-License-Identifier: BSD-3-Clause
 
-#ifndef SVMP_FE_MATH_DENSELINEARALGEBRA_H
-#define SVMP_FE_MATH_DENSELINEARALGEBRA_H
-
-#include "Types.h"
+#ifndef SVMP_COMMON_MATH_DENSELINEARALGEBRA_H
+#define SVMP_COMMON_MATH_DENSELINEARALGEBRA_H
 
 #include <cstddef>
 #include <limits>
@@ -14,7 +12,7 @@
 #include <string_view>
 #include <vector>
 
-namespace svmp::FE::math {
+namespace svmp::math {
 
 // Dense solve, inverse, rank, and pseudo-inverse support for FE construction
 // utilities. Matrices are row-major: matrix[row * cols + col]. The
@@ -23,7 +21,7 @@ namespace svmp::FE::math {
 
 /**
  * @brief Largest absolute entry of a dense matrix.
- * @ingroup FE_Math
+ * @ingroup SVMP_Math
  * @param matrix Row-major matrix entries.
  * @return Maximum of |entry| over all entries, or 0 for an empty matrix.
  */
@@ -31,7 +29,7 @@ namespace svmp::FE::math {
 
 /**
  * @brief Scale-aware pivot tolerance for dense factorization.
- * @ingroup FE_Math
+ * @ingroup SVMP_Math
  *
  * @details Proportional to machine epsilon scaled by the matrix size and
  * magnitude; pivots below it are treated as rank-deficient.
@@ -49,7 +47,7 @@ namespace svmp::FE::math {
 
 /**
  * @brief Scale-aware singular-value tolerance for rank decisions.
- * @ingroup FE_Math
+ * @ingroup SVMP_Math
  *
  * @details Singular values at or below the returned tolerance are treated as
  * zero when computing rank or a pseudo-inverse.
@@ -65,7 +63,7 @@ namespace svmp::FE::math {
                                                          double largest_singular_value,
                                                          double multiplier = double(64)) noexcept;
 
-/** @brief Result of a rank-revealing pseudo-inverse. @ingroup FE_Math */
+/** @brief Result of a rank-revealing pseudo-inverse. @ingroup SVMP_Math */
 struct DensePseudoInverseResult {
     std::vector<double> inverse;                 ///< Row-major pseudo-inverse.
     std::size_t rank{0};                         ///< Numerical rank at the chosen tolerance.
@@ -74,7 +72,7 @@ struct DensePseudoInverseResult {
     double smallest_retained_singular_value{0};  ///< Smallest singular value kept.
 };
 
-/** @brief SVD-based conditioning and rank diagnostics for a dense matrix. @ingroup FE_Math */
+/** @brief SVD-based conditioning and rank diagnostics for a dense matrix. @ingroup SVMP_Math */
 struct DenseMatrixDiagnostics {
     std::size_t rank{0};                         ///< Numerical rank at @ref tolerance.
     double tolerance{0};                         ///< Singular-value tolerance used.
@@ -83,7 +81,7 @@ struct DenseMatrixDiagnostics {
     double condition_estimate{std::numeric_limits<double>::infinity()};  ///< Condition estimate; infinite when rank-deficient.
 };
 
-/** @brief A dense inverse together with its diagnostics. @ingroup FE_Math */
+/** @brief A dense inverse together with its diagnostics. @ingroup SVMP_Math */
 struct DenseInverseResult {
     std::vector<double> inverse;        ///< Row-major inverse.
     DenseMatrixDiagnostics diagnostics; ///< Conditioning/rank diagnostics of the input.
@@ -92,20 +90,20 @@ struct DenseInverseResult {
 
 /**
  * @brief Condition estimate above which the inverse switches to an SVD fallback.
- * @ingroup FE_Math
+ * @ingroup SVMP_Math
  * @return The fallback condition-number threshold.
  */
 [[nodiscard]] double dense_matrix_condition_fallback_threshold() noexcept;
 /**
  * @brief Condition estimate above which validation rejects a dense inverse.
- * @ingroup FE_Math
+ * @ingroup SVMP_Math
  * @return The error condition-number threshold.
  */
 [[nodiscard]] double dense_matrix_condition_error_threshold() noexcept;
 
 /**
  * @brief LU factorization of a dense square matrix with a cached pivot summary.
- * @ingroup FE_Math
+ * @ingroup SVMP_Math
  *
  * @details Produced by factor_dense_matrix(); move-only because it owns the Eigen
  * factorization. @ref error_message_label prefixes the messages of exceptions
@@ -159,13 +157,13 @@ struct DenseLUSolver {
 
 /**
  * @brief SVD-based rank and conditioning diagnostics for a dense matrix.
- * @ingroup FE_Math
+ * @ingroup SVMP_Math
  * @param matrix Row-major matrix of size rows * cols.
  * @param rows Row count.
  * @param cols Column count.
  * @param error_message_label Prefix for the message of any exception thrown.
  * @return Rank, tolerance, singular-value, and condition diagnostics.
- * @throws FEException If the matrix size is inconsistent or the matrix is empty.
+ * @throws CoreException If the matrix size is inconsistent or the matrix is empty.
  */
 [[nodiscard]] DenseMatrixDiagnostics dense_matrix_diagnostics(
     std::span<const double> matrix,
@@ -175,12 +173,12 @@ struct DenseLUSolver {
 
 /**
  * @brief LU-factor a dense square matrix.
- * @ingroup FE_Math
+ * @ingroup SVMP_Math
  * @param matrix Row-major n * n matrix (consumed).
  * @param n Matrix dimension.
  * @param error_message_label Prefix for the message of any exception thrown.
  * @return The factorization.
- * @throws FEException If the size is inconsistent or the matrix is rank-deficient.
+ * @throws CoreException If the size is inconsistent or the matrix is rank-deficient.
  */
 [[nodiscard]] DenseLUSolver factor_dense_matrix(std::vector<double> matrix,
                                                 std::size_t n,
@@ -188,12 +186,12 @@ struct DenseLUSolver {
 
 /**
  * @brief Invert a dense square matrix.
- * @ingroup FE_Math
+ * @ingroup SVMP_Math
  * @param matrix Row-major n * n matrix (consumed).
  * @param n Matrix dimension.
  * @param error_message_label Prefix for the message of any exception thrown.
  * @return Row-major inverse of size n * n.
- * @throws FEException If the size is inconsistent or the matrix is singular.
+ * @throws CoreException If the size is inconsistent or the matrix is singular.
  */
 [[nodiscard]] std::vector<double> invert_dense_matrix(std::vector<double> matrix,
                                                     std::size_t n,
@@ -202,12 +200,12 @@ struct DenseLUSolver {
 /**
  * @brief Invert a dense square matrix with diagnostics, using an SVD fallback for
  * high-condition matrices.
- * @ingroup FE_Math
+ * @ingroup SVMP_Math
  * @param matrix Row-major n * n matrix (consumed).
  * @param n Matrix dimension.
  * @param error_message_label Prefix for the message of any exception thrown.
  * @return Inverse plus diagnostics and whether the SVD fallback was used.
- * @throws FEException If the size is inconsistent or the matrix is rank-deficient.
+ * @throws CoreException If the size is inconsistent or the matrix is rank-deficient.
  */
 [[nodiscard]] DenseInverseResult invert_dense_matrix_with_diagnostics(
     std::vector<double> matrix,
@@ -216,12 +214,12 @@ struct DenseLUSolver {
 
 /**
  * @brief Validate that a dense inverse has full rank and acceptable conditioning.
- * @ingroup FE_Math
+ * @ingroup SVMP_Math
  * @param result Result from invert_dense_matrix_with_diagnostics().
  * @param expected_rank Required (full) rank.
  * @param error_message_label Prefix for the message of any exception thrown.
  * @param max_condition Largest acceptable condition estimate.
- * @throws FEException If the rank is below expected_rank or the condition exceeds max_condition.
+ * @throws CoreException If the rank is below expected_rank or the condition exceeds max_condition.
  */
 void validate_dense_inverse_diagnostics(
     const DenseInverseResult& result,
@@ -231,12 +229,12 @@ void validate_dense_inverse_diagnostics(
 
 /**
  * @brief Numerical rank of a dense matrix from its singular values.
- * @ingroup FE_Math
+ * @ingroup SVMP_Math
  * @param matrix Row-major matrix of size rows * cols (consumed).
  * @param rows Row count.
  * @param cols Column count.
  * @return Number of singular values above the scale-aware tolerance.
- * @throws FEException If the matrix size is inconsistent.
+ * @throws CoreException If the matrix size is inconsistent.
  */
 [[nodiscard]] std::size_t dense_matrix_rank(std::vector<double> matrix,
                                             std::size_t rows,
@@ -244,13 +242,13 @@ void validate_dense_inverse_diagnostics(
 
 /**
  * @brief Moore-Penrose pseudo-inverse via a rank-revealing SVD.
- * @ingroup FE_Math
+ * @ingroup SVMP_Math
  * @param matrix Row-major matrix of size rows * cols.
  * @param rows Row count.
  * @param cols Column count.
  * @param error_message_label Prefix for the message of any exception thrown.
  * @return Row-major pseudo-inverse (cols * rows) plus rank/tolerance diagnostics.
- * @throws FEException If the matrix size is inconsistent or the matrix is empty.
+ * @throws CoreException If the matrix size is inconsistent or the matrix is empty.
  */
 [[nodiscard]] DensePseudoInverseResult rank_revealing_pseudo_inverse(
     std::span<const double> matrix,
@@ -258,6 +256,6 @@ void validate_dense_inverse_diagnostics(
     std::size_t cols,
     std::string_view error_message_label = "dense matrix");
 
-} // namespace svmp::FE::math
+} // namespace svmp::math
 
-#endif // SVMP_FE_MATH_DENSELINEARALGEBRA_H
+#endif // SVMP_COMMON_MATH_DENSELINEARALGEBRA_H
diff --git a/Code/Source/solver/FE/Math/Matrix.h b/Code/Source/solver/common/math/include/math/Matrix.h
similarity index 82%
rename from Code/Source/solver/FE/Math/Matrix.h
rename to Code/Source/solver/common/math/include/math/Matrix.h
index e61fb4f42..641f4300b 100644
--- a/Code/Source/solver/FE/Math/Matrix.h
+++ b/Code/Source/solver/common/math/include/math/Matrix.h
@@ -1,8 +1,8 @@
 // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others.
 // SPDX-License-Identifier: BSD-3-Clause
 
-#ifndef SVMP_FE_MATH_MATRIX_H
-#define SVMP_FE_MATH_MATRIX_H
+#ifndef SVMP_COMMON_MATH_MATRIX_H
+#define SVMP_COMMON_MATH_MATRIX_H
 
 /**
  * @file Matrix.h
@@ -23,16 +23,16 @@
 #include <cstddef>
 
 /**
- * @defgroup FE_MatrixMath Matrix
- * @ingroup FE_Math
+ * @defgroup SVMP_Math_Matrix Matrix
+ * @ingroup SVMP_Math
  * @brief Fixed-size matrix type aliases.
  */
 
-namespace svmp::FE::math {
+namespace svmp::math {
 
 /**
  * @brief Fixed-size matrix for element-level computations
- * @ingroup FE_MatrixMath
+ * @ingroup SVMP_Math_Matrix
  * @tparam T Scalar type (float, double)
  * @tparam M Number of rows
  * @tparam N Number of columns
@@ -40,6 +40,6 @@ namespace svmp::FE::math {
 template<typename T, std::size_t M, std::size_t N>
 using Matrix = Eigen::Matrix<T, static_cast<int>(M), static_cast<int>(N)>;
 
-} // namespace svmp::FE::math
+} // namespace svmp::math
 
-#endif // SVMP_FE_MATH_MATRIX_H
+#endif // SVMP_COMMON_MATH_MATRIX_H
diff --git a/Code/Source/solver/FE/Math/Vector.h b/Code/Source/solver/common/math/include/math/Vector.h
similarity index 81%
rename from Code/Source/solver/FE/Math/Vector.h
rename to Code/Source/solver/common/math/include/math/Vector.h
index 6f09e6ac9..d3ab5bf0c 100644
--- a/Code/Source/solver/FE/Math/Vector.h
+++ b/Code/Source/solver/common/math/include/math/Vector.h
@@ -1,8 +1,8 @@
 // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others.
 // SPDX-License-Identifier: BSD-3-Clause
 
-#ifndef SVMP_FE_MATH_VECTOR_H
-#define SVMP_FE_MATH_VECTOR_H
+#ifndef SVMP_COMMON_MATH_VECTOR_H
+#define SVMP_COMMON_MATH_VECTOR_H
 
 /**
  * @file Vector.h
@@ -15,7 +15,7 @@
  * zeroed value is required.
  *
  * This is a small, fixed-size (compile-time length) vector for element-level FE
- * kernels in namespace svmp::FE::math. It is distinct from, and not a replacement
+ * kernels in namespace svmp::math. It is distinct from, and not a replacement
  * for, the legacy dynamically sized container in solver/Vector.h: the two differ
  * in namespace, size model (compile-time vs runtime), and memory management, and
  * coexist deliberately.
@@ -26,30 +26,29 @@
 #include <cstddef>
 
 /**
- * @defgroup FE_Math Math
- * @ingroup FE
+ * @defgroup SVMP_Math Math
  * @brief Linear algebra vocabulary types and dense utilities for finite-element computations.
  *
  * @details The Math module defines the fixed-size vector and matrix types
  * used in element-level kernels (as aliases of Eigen types) and dense linear
  * algebra utilities used by basis construction and local transforms.
  *
- * @defgroup FE_VectorMath Vector
- * @ingroup FE_Math
+ * @defgroup SVMP_Math_Vector Vector
+ * @ingroup SVMP_Math
  * @brief Fixed-size vector type aliases.
  */
 
-namespace svmp::FE::math {
+namespace svmp::math {
 
 /**
  * @brief Fixed-size column vector for element-level computations
- * @ingroup FE_VectorMath
+ * @ingroup SVMP_Math_Vector
  * @tparam T Scalar type (float, double)
  * @tparam N Vector dimension
  */
 template<typename T, std::size_t N>
 using Vector = Eigen::Matrix<T, static_cast<int>(N), 1>;
 
-} // namespace svmp::FE::math
+} // namespace svmp::math
 
-#endif // SVMP_FE_MATH_VECTOR_H
+#endif // SVMP_COMMON_MATH_VECTOR_H
diff --git a/Code/Source/solver/FE/Math/DenseLinearAlgebra.cpp b/Code/Source/solver/common/math/src/DenseLinearAlgebra.cpp
similarity index 95%
rename from Code/Source/solver/FE/Math/DenseLinearAlgebra.cpp
rename to Code/Source/solver/common/math/src/DenseLinearAlgebra.cpp
index d98a9279b..f7ed794ca 100644
--- a/Code/Source/solver/FE/Math/DenseLinearAlgebra.cpp
+++ b/Code/Source/solver/common/math/src/DenseLinearAlgebra.cpp
@@ -1,9 +1,9 @@
 // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others.
 // SPDX-License-Identifier: BSD-3-Clause
 
-#include "DenseLinearAlgebra.h"
+#include "math/DenseLinearAlgebra.h"
 
-#include "FEException.h"
+#include "Core/Exception.h"
 
 #include <Eigen/Dense>
 
@@ -13,7 +13,7 @@
 #include <string>
 #include <utility>
 
-namespace svmp::FE::math {
+namespace svmp::math {
 
 namespace {
 
@@ -90,11 +90,11 @@ void DenseLUSolver::solve_in_place(std::span<double> rhs) const {
 
 void DenseLUSolver::solve_in_place(std::span<double> rhs,
                                    std::size_t rhs_count) const {
-    ::svmp::check<FEException>(
+    check<CoreException>(
         rhs_count > 0, error_message_label + ": dense solve requires at least one right-hand side");
-    ::svmp::check<FEException>(
+    check<CoreException>(
         rhs.size() == n * rhs_count, error_message_label + ": dense multi-RHS solve size mismatch");
-    ::svmp::check<FEException>(
+    check<CoreException>(
         impl && impl->lu.rows() == static_cast<Eigen::Index>(n), error_message_label + ": dense solver is not factorized");
     if (n == 0) {
         return;
@@ -119,9 +119,9 @@ DenseMatrixDiagnostics dense_matrix_diagnostics(
     std::size_t rows,
     std::size_t cols,
     std::string_view error_message_label) {
-    ::svmp::check<FEException>(
+    check<CoreException>(
         matrix.size() == rows * cols, std::string(error_message_label) + ": diagnostic size mismatch");
-    ::svmp::check<FEException>(
+    check<CoreException>(
         rows > 0 && cols > 0, std::string(error_message_label) + ": diagnostics require a nonempty matrix");
 
     const DenseMatrix dense = map_row_major(matrix, rows, cols);
@@ -157,7 +157,7 @@ DenseMatrixDiagnostics dense_matrix_diagnostics(
 DenseLUSolver factor_dense_matrix(std::vector<double> matrix,
                                   std::size_t n,
                                   std::string_view error_message_label) {
-    ::svmp::check<FEException>(
+    check<CoreException>(
         matrix.size() == n * n, std::string(error_message_label) + ": dense factorization size mismatch");
 
     DenseLUSolver solver;
@@ -177,7 +177,7 @@ DenseLUSolver factor_dense_matrix(std::vector<double> matrix,
     const auto diagonal = solver.impl->lu.matrixLU().diagonal();
     for (Eigen::Index col = 0; col < diagonal.size(); ++col) {
         const double pivot_magnitude = std::abs(diagonal[col]);
-        ::svmp::check<FEException>(
+        check<CoreException>(
             pivot_magnitude > solver.pivot_tolerance, solver.error_message_label + ": rank-deficient matrix (rank " +
                 std::to_string(col) + " of " + std::to_string(n) +
                 ", pivot below scale-aware tolerance " +
@@ -200,7 +200,7 @@ DenseInverseResult invert_dense_matrix_with_diagnostics(
     std::vector<double> matrix,
     std::size_t n,
     std::string_view error_message_label) {
-    ::svmp::check<FEException>(
+    check<CoreException>(
         matrix.size() == n * n, std::string(error_message_label) + ": dense inverse size mismatch");
     std::vector<double> matrix_for_lu = matrix;
     const DenseLUSolver solver =
@@ -225,7 +225,7 @@ DenseInverseResult invert_dense_matrix_with_diagnostics(
             // singular value (rank < full_rank). A sub-tolerance singular value
             // therefore cannot reach here in current code; the guard protects
             // against future refactors that derive the fallback condition differently.
-            ::svmp::check<FEException>(
+            check<CoreException>(
                 singular_values[i] > result.diagnostics.tolerance, std::string(error_message_label) + ": high-condition SVD fallback encountered a dropped singular value");
             sigma_inverse(i, i) = double(1) / singular_values[i];
         }
@@ -245,7 +245,7 @@ void validate_dense_inverse_diagnostics(
     std::size_t expected_rank,
     std::string_view error_message_label,
     double max_condition) {
-    ::svmp::check<FEException>(
+    check<CoreException>(
         result.diagnostics.rank == expected_rank, std::string(error_message_label) + ": rank-deficient matrix (rank " +
             std::to_string(result.diagnostics.rank) + " of " +
             std::to_string(expected_rank) + ")");
@@ -254,7 +254,7 @@ void validate_dense_inverse_diagnostics(
         return;
     }
 
-    ::svmp::check<FEException>(
+    check<CoreException>(
         result.diagnostics.condition_estimate <= max_condition, std::string(error_message_label) + ": condition estimate " +
             std::to_string(result.diagnostics.condition_estimate) +
             " exceeds supported threshold " + std::to_string(max_condition));
@@ -273,7 +273,7 @@ std::vector<double> invert_dense_matrix(std::vector<double> matrix,
 std::size_t dense_matrix_rank(std::vector<double> matrix,
                               std::size_t rows,
                               std::size_t cols) {
-    ::svmp::check<FEException>(
+    check<CoreException>(
         matrix.size() == rows * cols, "dense_matrix_rank: size mismatch");
 
     const DenseMatrix dense =
@@ -300,9 +300,9 @@ DensePseudoInverseResult rank_revealing_pseudo_inverse(
     std::size_t rows,
     std::size_t cols,
     std::string_view error_message_label) {
-    ::svmp::check<FEException>(
+    check<CoreException>(
         matrix.size() == rows * cols, std::string(error_message_label) + ": pseudo-inverse size mismatch");
-    ::svmp::check<FEException>(
+    check<CoreException>(
         rows > 0 && cols > 0, std::string(error_message_label) + ": pseudo-inverse requires a nonempty matrix");
 
     const DenseMatrix dense = map_row_major(matrix, rows, cols);
@@ -334,4 +334,4 @@ DensePseudoInverseResult rank_revealing_pseudo_inverse(
     return result;
 }
 
-} // namespace svmp::FE::math
+} // namespace svmp::math
diff --git a/Code/Source/solver/nn.cpp b/Code/Source/solver/nn.cpp
index 5e2893fda..b42ce607a 100644
--- a/Code/Source/solver/nn.cpp
+++ b/Code/Source/solver/nn.cpp
@@ -213,7 +213,7 @@ std::size_t basis_index_for_solver_node(consts::ElementType eType, const int sol
 /// xi array, zero-filling the trailing components that are inactive for
 /// lower-dimensional elements. Throws BasisConfigurationException when xi has
 /// fewer rows than the basis reference dimension.
-fe::math::Vector<double, 3> make_basis_point(const febasis::BasisFunction& basis,
+svmp::math::Vector<double, 3> make_basis_point(const febasis::BasisFunction& basis,
                                                const int g,
                                                const Array<double>& xi)
 {
@@ -225,7 +225,7 @@ fe::math::Vector<double, 3> make_basis_point(const febasis::BasisFunction& basis
 
   // Inactive trailing components must be zero for lower-dimensional elements;
   // Eigen-backed vectors are not zero-initialized by default construction.
-  fe::math::Vector<double, 3> point = fe::math::Vector<double, 3>::Zero();
+  svmp::math::Vector<double, 3> point = svmp::math::Vector<double, 3>::Zero();
   for (int d = 0; d < basis.dimension(); ++d) {
     point[static_cast<std::size_t>(d)] = xi(d, g);
   }
diff --git a/Code/Source/solver/post.cpp b/Code/Source/solver/post.cpp
index ffc9b9c8d..b6765816f 100644
--- a/Code/Source/solver/post.cpp
+++ b/Code/Source/solver/post.cpp
@@ -6,7 +6,7 @@
 #include "ActiveStress.h"
 #include "Core/Exception.h"
 #include "FE/Common/FEException.h"
-#include "FE/Math/DenseLinearAlgebra.h"
+#include "math/DenseLinearAlgebra.h"
 #include "all_fun.h"
 #include "porous/darcy.h"
 #include "fluid/fluid.h"
@@ -1209,7 +1209,7 @@ void post(Simulation* simulation, const mshType& lM, Array<double>& res, const S
       for (double& value : mass) {
         value /= element_volume;
       }
-      svmp::FE::math::factor_dense_matrix(mass, eNoN, "Darcy flux mass matrix").solve_in_place(flux_rhs, nsd);
+      svmp::math::factor_dense_matrix(mass, eNoN, "Darcy flux mass matrix").solve_in_place(flux_rhs, nsd);
 
       for (int a = 0; a < eNoN; ++a) {
         const int Ac = lM.IEN(a,e);
diff --git a/tests/unitTests/FE/Basis/test_BasisErrorPaths.cpp b/tests/unitTests/FE/Basis/test_BasisErrorPaths.cpp
index e81a9c377..536e9907a 100644
--- a/tests/unitTests/FE/Basis/test_BasisErrorPaths.cpp
+++ b/tests/unitTests/FE/Basis/test_BasisErrorPaths.cpp
@@ -48,7 +48,7 @@ class MinimalScalarBasis : public BasisFunction {
     int order() const noexcept override { return 1; }
     std::size_t size() const noexcept override { return 2u; }
 
-    void evaluate_values_to(const math::Vector<double, 3>&,
+    void evaluate_values_to(const svmp::math::Vector<double, 3>&,
                             std::span<double> values_out) const override
     {
         std::fill(values_out.begin(), values_out.end(), double(0));
@@ -70,7 +70,7 @@ class ExactQuadraticBasis : public BasisFunction {
     int order() const noexcept override { return 2; }
     std::size_t size() const noexcept override { return 2u; }
 
-    void evaluate_values_to(const math::Vector<double, 3>& xi,
+    void evaluate_values_to(const svmp::math::Vector<double, 3>& xi,
                             std::span<double> values_out) const override
     {
         const double x = xi[0];
@@ -84,7 +84,7 @@ class ExactQuadraticBasis : public BasisFunction {
                         x * y + x * z - y * z;
     }
 
-    void evaluate_gradients_to(const math::Vector<double, 3>& xi,
+    void evaluate_gradients_to(const svmp::math::Vector<double, 3>& xi,
                                std::span<Gradient> gradients_out) const override
     {
         const double x = xi[0];
@@ -122,14 +122,14 @@ class SpanPrimitiveBasis : public BasisFunction {
     int order() const noexcept override { return 1; }
     std::size_t size() const noexcept override { return 2u; }
 
-    void evaluate_values_to(const math::Vector<double, 3>& xi,
+    void evaluate_values_to(const svmp::math::Vector<double, 3>& xi,
                             std::span<double> values_out) const override
     {
         values_out[0] = double(1) + xi[0];
         values_out[1] = double(2) + xi[1];
     }
 
-    void evaluate_gradients_to(const math::Vector<double, 3>&,
+    void evaluate_gradients_to(const svmp::math::Vector<double, 3>&,
                                std::span<Gradient> gradients_out) const override
     {
         gradients_out[0] = Gradient::Zero();
@@ -138,7 +138,7 @@ class SpanPrimitiveBasis : public BasisFunction {
         gradients_out[1][1] = double(1);
     }
 
-    void evaluate_hessians_to(const math::Vector<double, 3>& xi,
+    void evaluate_hessians_to(const svmp::math::Vector<double, 3>& xi,
                               std::span<Hessian> hessians_out) const override
     {
         for (std::size_t d = 0; d < size(); ++d) {
@@ -350,7 +350,7 @@ TEST(BasisErrorPaths, NodeOrderingInvalidNodeThrows) {
 
 TEST(BasisErrorPaths, BasisFunctionDefaultsThrowForMissingDerivatives) {
     MinimalScalarBasis basis;
-    const math::Vector<double, 3> xi{double(0), double(0), double(0)};
+    const svmp::math::Vector<double, 3> xi{double(0), double(0), double(0)};
     std::vector<Gradient> gradients;
     std::vector<Hessian> hessians;
 
@@ -360,7 +360,7 @@ TEST(BasisErrorPaths, BasisFunctionDefaultsThrowForMissingDerivatives) {
 
 TEST(BasisErrorPaths, NumericalDerivativeHelpersMatchAnalyticDerivatives) {
     ExactQuadraticBasis basis;
-    const math::Vector<double, 3> xi{double(0.2), double(-0.35), double(0.4)};
+    const svmp::math::Vector<double, 3> xi{double(0.2), double(-0.35), double(0.4)};
 
     // On a quadratic, centered differences are exact except for the round-off
     // floor ~ eps_machine/step. The tolerances below are a few times
@@ -401,7 +401,7 @@ TEST(BasisErrorPaths, NumericalDerivativeHelpersMatchAnalyticDerivatives) {
 
 TEST(BasisErrorPaths, BasisFunctionVectorOverloadsForwardToSpanPrimitives) {
     SpanPrimitiveBasis basis;
-    const math::Vector<double, 3> point{double(0.25), double(0.5), double(-0.25)};
+    const svmp::math::Vector<double, 3> point{double(0.25), double(0.5), double(-0.25)};
 
     // Reference results taken directly from the span primitives the basis defines.
     std::vector<double> span_values(basis.size());
diff --git a/tests/unitTests/FE/Basis/test_BasisHessians.cpp b/tests/unitTests/FE/Basis/test_BasisHessians.cpp
index fd66a9e74..8eb39a5ed 100644
--- a/tests/unitTests/FE/Basis/test_BasisHessians.cpp
+++ b/tests/unitTests/FE/Basis/test_BasisHessians.cpp
@@ -26,7 +26,7 @@ namespace {
 constexpr double kHessianInvariantTol = double(1e-12);
 
 void numerical_gradient_helper(const BasisFunction& basis,
-                               const math::Vector<double, 3>& xi,
+                               const svmp::math::Vector<double, 3>& xi,
                                std::vector<Gradient>& gradients,
                                double eps = double(1e-6))
 {
@@ -36,8 +36,8 @@ void numerical_gradient_helper(const BasisFunction& basis,
 
     for (int d = 0; d < basis.dimension(); ++d) {
         const std::size_t sd = static_cast<std::size_t>(d);
-        math::Vector<double, 3> xi_p = xi;
-        math::Vector<double, 3> xi_m = xi;
+        svmp::math::Vector<double, 3> xi_p = xi;
+        svmp::math::Vector<double, 3> xi_m = xi;
         xi_p[sd] += eps;
         xi_m[sd] -= eps;
 
@@ -53,7 +53,7 @@ void numerical_gradient_helper(const BasisFunction& basis,
 }
 
 void numerical_hessian_helper(const BasisFunction& basis,
-                              const math::Vector<double, 3>& xi,
+                              const svmp::math::Vector<double, 3>& xi,
                               std::vector<Hessian>& hessians,
                               double eps = double(1e-5))
 {
@@ -62,8 +62,8 @@ void numerical_hessian_helper(const BasisFunction& basis,
 
     for (int i = 0; i < dim; ++i) {
         for (int j = 0; j < dim; ++j) {
-            math::Vector<double, 3> xi_p = xi;
-            math::Vector<double, 3> xi_m = xi;
+            svmp::math::Vector<double, 3> xi_p = xi;
+            svmp::math::Vector<double, 3> xi_m = xi;
             const std::size_t sj = static_cast<std::size_t>(j);
             xi_p[sj] += eps;
             xi_m[sj] -= eps;
@@ -81,7 +81,7 @@ void numerical_hessian_helper(const BasisFunction& basis,
     }
 }
 
-std::vector<math::Vector<double, 3>> sample_points_for(BasisTopology topology) {
+std::vector<svmp::math::Vector<double, 3>> sample_points_for(BasisTopology topology) {
     switch (topology) {
         case BasisTopology::Line:
             return {{double(-0.35), double(0), double(0)}, {double(0.2), double(0), double(0)}};
@@ -101,7 +101,7 @@ std::vector<math::Vector<double, 3>> sample_points_for(BasisTopology topology) {
 }
 
 void expect_gradients_match_numerical(const BasisFunction& basis,
-                                      const std::vector<math::Vector<double, 3>>& points,
+                                      const std::vector<svmp::math::Vector<double, 3>>& points,
                                       double tol,
                                       double eps = double(1e-6))
 {
@@ -125,7 +125,7 @@ void expect_gradients_match_numerical(const BasisFunction& basis,
 }
 
 void expect_hessians_match_numerical(const BasisFunction& basis,
-                                     const std::vector<math::Vector<double, 3>>& points,
+                                     const std::vector<svmp::math::Vector<double, 3>>& points,
                                      double tol,
                                      double eps = double(1e-5))
 {
@@ -152,7 +152,7 @@ void expect_hessians_match_numerical(const BasisFunction& basis,
 }
 
 void expect_partition_hessian_sum_zero(const BasisFunction& basis,
-                                       const math::Vector<double, 3>& xi,
+                                       const svmp::math::Vector<double, 3>& xi,
                                        double tol)
 {
     std::vector<Hessian> hessians;
@@ -179,7 +179,7 @@ void expect_partition_hessian_sum_zero(const BasisFunction& basis,
 }
 
 void expect_hessians_symmetric(const BasisFunction& basis,
-                               const math::Vector<double, 3>& xi,
+                               const svmp::math::Vector<double, 3>& xi,
                                double tol)
 {
     std::vector<Hessian> hessians;
@@ -197,7 +197,7 @@ void expect_hessians_symmetric(const BasisFunction& basis,
 }
 
 void expect_inactive_z_derivatives_zero(const BasisFunction& basis,
-                                        const std::vector<math::Vector<double, 3>>& points,
+                                        const std::vector<svmp::math::Vector<double, 3>>& points,
                                         double tol)
 {
     ASSERT_EQ(basis.dimension(), 2);
@@ -228,7 +228,7 @@ void expect_inactive_z_derivatives_zero(const BasisFunction& basis,
     }
 }
 
-std::vector<math::Vector<double, 3>> serendipity_sample_points(BasisTopology topology) {
+std::vector<svmp::math::Vector<double, 3>> serendipity_sample_points(BasisTopology topology) {
     if (topology == BasisTopology::Quadrilateral) {
         return {{double(0.17), double(-0.31), double(0)}, {double(-0.45), double(0.25), double(0)}};
     }
@@ -265,7 +265,7 @@ TEST(BasisHessians, LagrangeHessiansSumToZeroAndAreSymmetric) {
     const struct Case {
         BasisTopology topology;
         int order;
-        math::Vector<double, 3> xi;
+        svmp::math::Vector<double, 3> xi;
     } cases[] = {
         {BasisTopology::Line, 3, {double(0.15), double(0), double(0)}},
         {BasisTopology::Triangle, 3, {double(0.2), double(0.25), double(0)}},
@@ -286,7 +286,7 @@ TEST(BasisHessians, SerendipityHessiansSumToZeroAndAreSymmetric) {
     const struct Case {
         ElementType type;
         int order;
-        math::Vector<double, 3> xi;
+        svmp::math::Vector<double, 3> xi;
     } cases[] = {
         {ElementType::Quad8, 2, {double(0.17), double(-0.31), double(0)}},
         {ElementType::Hex20, 2, {double(0.2), double(-0.1), double(0.3)}},
@@ -309,7 +309,7 @@ TEST(BasisHessians, SolverMappedVolumeSelectionsSatisfyInvariants) {
         ElementType type;
         BasisType basis_type;
         int order;
-        math::Vector<double, 3> xi;
+        svmp::math::Vector<double, 3> xi;
     } cases[] = {
         {ElementType::Line2, BasisType::Lagrange, 1, {double(0.15), double(0), double(0)}},
         {ElementType::Line3, BasisType::Lagrange, 2, {double(-0.25), double(0), double(0)}},
diff --git a/tests/unitTests/FE/Basis/test_HigherOrderWedge.cpp b/tests/unitTests/FE/Basis/test_HigherOrderWedge.cpp
index dd4511b18..490009718 100644
--- a/tests/unitTests/FE/Basis/test_HigherOrderWedge.cpp
+++ b/tests/unitTests/FE/Basis/test_HigherOrderWedge.cpp
@@ -16,8 +16,8 @@ using namespace svmp::FE::basis;
 
 namespace {
 
-void expect_nodes_close(const std::vector<math::Vector<double, 3>>& lhs,
-                        const std::vector<math::Vector<double, 3>>& rhs,
+void expect_nodes_close(const std::vector<svmp::math::Vector<double, 3>>& lhs,
+                        const std::vector<svmp::math::Vector<double, 3>>& rhs,
                         double tol)
 {
     ASSERT_EQ(lhs.size(), rhs.size());
@@ -46,7 +46,7 @@ void expect_kronecker_at_nodes(const LagrangeBasis& basis, double tol)
 }
 
 void expect_partition_gradient_hessian_sums(const LagrangeBasis& basis,
-                                            const std::vector<math::Vector<double, 3>>& points,
+                                            const std::vector<svmp::math::Vector<double, 3>>& points,
                                             double value_tol,
                                             double derivative_tol)
 {
@@ -83,7 +83,7 @@ void expect_partition_gradient_hessian_sums(const LagrangeBasis& basis,
 }
 
 void expect_all_entries_finite(const LagrangeBasis& basis,
-                               const math::Vector<double, 3>& xi)
+                               const svmp::math::Vector<double, 3>& xi)
 {
     std::vector<double> values;
     std::vector<Gradient> gradients;
diff --git a/tests/unitTests/FE/Basis/test_LagrangeBasis.cpp b/tests/unitTests/FE/Basis/test_LagrangeBasis.cpp
index a24bcb96a..ff8bb35e1 100644
--- a/tests/unitTests/FE/Basis/test_LagrangeBasis.cpp
+++ b/tests/unitTests/FE/Basis/test_LagrangeBasis.cpp
@@ -21,7 +21,7 @@ using namespace svmp::FE::basis;
 
 namespace {
 
-using Point = math::Vector<double, 3>;
+using Point = svmp::math::Vector<double, 3>;
 
 // Tolerance convention (set from measured residuals, not guessed). Identities that
 // share a generator -- node coordinates, the lattice forward-image, span vs vector
diff --git a/tests/unitTests/FE/Basis/test_SerendipityBasis.cpp b/tests/unitTests/FE/Basis/test_SerendipityBasis.cpp
index b764c8a8a..1aff8fff7 100644
--- a/tests/unitTests/FE/Basis/test_SerendipityBasis.cpp
+++ b/tests/unitTests/FE/Basis/test_SerendipityBasis.cpp
@@ -8,7 +8,7 @@
 #include "FE/Basis/LagrangeBasis.h"
 #include "FE/Basis/NodeOrderingConventions.h"
 #include "FE/Basis/SerendipityBasis.h"
-#include "FE/Math/DenseLinearAlgebra.h"
+#include "math/DenseLinearAlgebra.h"
 
 #include <algorithm>
 #include <array>
@@ -21,7 +21,7 @@ using namespace svmp::FE::basis;
 namespace {
 
 void expect_partition_of_unity(const SerendipityBasis& basis,
-                               const math::Vector<double, 3>& xi,
+                               const svmp::math::Vector<double, 3>& xi,
                                double tolerance = double(1e-10))
 {
     std::vector<double> values;
@@ -47,7 +47,7 @@ void expect_partition_of_unity(const SerendipityBasis& basis,
 }
 
 void expect_nodal_delta(const SerendipityBasis& basis,
-                        const std::vector<math::Vector<double, 3>>& nodes,
+                        const std::vector<svmp::math::Vector<double, 3>>& nodes,
                         double tolerance)
 {
     ASSERT_EQ(nodes.size(), basis.size());
@@ -62,10 +62,10 @@ void expect_nodal_delta(const SerendipityBasis& basis,
     }
 }
 
-std::vector<math::Vector<double, 3>> reference_nodes(ElementType type,
+std::vector<svmp::math::Vector<double, 3>> reference_nodes(ElementType type,
                                                    std::size_t count)
 {
-    std::vector<math::Vector<double, 3>> nodes;
+    std::vector<svmp::math::Vector<double, 3>> nodes;
     nodes.reserve(count);
     for (std::size_t i = 0; i < count; ++i) {
         nodes.push_back(ReferenceNodeLayout::node_coord_at(type, i));
@@ -75,7 +75,7 @@ std::vector<math::Vector<double, 3>> reference_nodes(ElementType type,
 
 template<typename Function>
 double interpolate_nodal_function(const SerendipityBasis& basis,
-                                const math::Vector<double, 3>& xi,
+                                const svmp::math::Vector<double, 3>& xi,
                                 Function&& nodal_function)
 {
     std::vector<double> values;
@@ -129,14 +129,14 @@ double integer_power_for_test(double base, int exponent) {
     return result;
 }
 
-double monomial_value_for_test(const math::Vector<double, 3>& p,
+double monomial_value_for_test(const svmp::math::Vector<double, 3>& p,
                              const std::array<int, 2>& exponent) {
     return integer_power_for_test(p[0], exponent[0]) *
            integer_power_for_test(p[1], exponent[1]);
 }
 
 std::vector<double> quadrilateral_vandermonde_for_test(
-    const std::vector<math::Vector<double, 3>>& nodes,
+    const std::vector<svmp::math::Vector<double, 3>>& nodes,
     const std::vector<std::array<int, 2>>& exponents)
 {
     const std::size_t n = nodes.size();
@@ -150,7 +150,7 @@ std::vector<double> quadrilateral_vandermonde_for_test(
     return vandermonde;
 }
 
-void expect_no_duplicate_nodes(const std::vector<math::Vector<double, 3>>& nodes,
+void expect_no_duplicate_nodes(const std::vector<svmp::math::Vector<double, 3>>& nodes,
                                double tolerance)
 {
     for (std::size_t a = 0; a < nodes.size(); ++a) {
@@ -163,8 +163,8 @@ void expect_no_duplicate_nodes(const std::vector<math::Vector<double, 3>>& nodes
     }
 }
 
-void expect_nodes_near(const std::vector<math::Vector<double, 3>>& actual,
-                       const std::vector<math::Vector<double, 3>>& expected,
+void expect_nodes_near(const std::vector<svmp::math::Vector<double, 3>>& actual,
+                       const std::vector<svmp::math::Vector<double, 3>>& expected,
                        double tolerance)
 {
     ASSERT_EQ(actual.size(), expected.size());
@@ -178,7 +178,7 @@ void expect_nodes_near(const std::vector<math::Vector<double, 3>>& actual,
 
 // Every monomial here has superlinear degree at most three, so it lies in the
 // order-three quadrilateral serendipity space.
-double cubic_serendipity_function(const math::Vector<double, 3>& p) {
+double cubic_serendipity_function(const svmp::math::Vector<double, 3>& p) {
     const double x = p[0];
     const double y = p[1];
     return double(1) + double(2) * x - y + double(3) * x * y +
@@ -186,7 +186,7 @@ double cubic_serendipity_function(const math::Vector<double, 3>& p) {
            double(0.5) * x * x * x * y - double(0.25) * x * y * y * y;
 }
 
-double bilinear_function(const math::Vector<double, 3>& p) {
+double bilinear_function(const svmp::math::Vector<double, 3>& p) {
     return double(2) - double(3) * p[0] + double(4) * p[1] + double(0.5) * p[0] * p[1];
 }
 
@@ -197,7 +197,7 @@ double bilinear_function(const math::Vector<double, 3>& p) {
 // each serendipity space, independently of the production tables in
 // SerendipityBasis.cpp
 
-double monomial_value_3d_for_test(const math::Vector<double, 3>& p,
+double monomial_value_3d_for_test(const svmp::math::Vector<double, 3>& p,
                                 const std::array<int, 3>& exponent) {
     return integer_power_for_test(p[0], exponent[0]) *
            integer_power_for_test(p[1], exponent[1]) *
@@ -205,7 +205,7 @@ double monomial_value_3d_for_test(const math::Vector<double, 3>& p,
 }
 
 std::vector<double> vandermonde_3d_for_test(
-    const std::vector<math::Vector<double, 3>>& nodes,
+    const std::vector<svmp::math::Vector<double, 3>>& nodes,
     const std::vector<std::array<int, 3>>& exponents) {
     const std::size_t n = nodes.size();
     std::vector<double> vandermonde(n * n, double(0));
@@ -300,8 +300,8 @@ std::vector<std::array<int, 3>> wedge15_serendipity_exponents_for_test() {
 
 // Independent Hex20 reference layout: the eight cube corners followed by the
 // twelve edge midpoints, in the corner/edge order the reference layout uses.
-std::vector<math::Vector<double, 3>> hex20_reference_nodes_for_test() {
-    std::vector<math::Vector<double, 3>> corners;
+std::vector<svmp::math::Vector<double, 3>> hex20_reference_nodes_for_test() {
+    std::vector<svmp::math::Vector<double, 3>> corners;
     corners.push_back({double(-1), double(-1), double(-1)});
     corners.push_back({double(1), double(-1), double(-1)});
     corners.push_back({double(1), double(1), double(-1)});
@@ -317,9 +317,9 @@ std::vector<math::Vector<double, 3>> hex20_reference_nodes_for_test() {
         {0, 4}, {1, 5}, {2, 6}, {3, 7},
     };
 
-    std::vector<math::Vector<double, 3>> nodes = corners;
+    std::vector<svmp::math::Vector<double, 3>> nodes = corners;
     for (const auto& edge : edges) {
-        const math::Vector<double, 3> midpoint =
+        const svmp::math::Vector<double, 3> midpoint =
             (corners[static_cast<std::size_t>(edge[0])] +
              corners[static_cast<std::size_t>(edge[1])]) * double(0.5);
         nodes.push_back(midpoint);
@@ -329,8 +329,8 @@ std::vector<math::Vector<double, 3>> hex20_reference_nodes_for_test() {
 
 // Independent Wedge15 reference layout: the six prism corners followed by the
 // nine edge midpoints, in reference-layout order.
-std::vector<math::Vector<double, 3>> wedge15_reference_nodes_for_test() {
-    std::vector<math::Vector<double, 3>> corners;
+std::vector<svmp::math::Vector<double, 3>> wedge15_reference_nodes_for_test() {
+    std::vector<svmp::math::Vector<double, 3>> corners;
     corners.push_back({double(0), double(0), double(-1)});
     corners.push_back({double(1), double(0), double(-1)});
     corners.push_back({double(0), double(1), double(-1)});
@@ -344,9 +344,9 @@ std::vector<math::Vector<double, 3>> wedge15_reference_nodes_for_test() {
         {0, 3}, {1, 4}, {2, 5},
     };
 
-    std::vector<math::Vector<double, 3>> nodes = corners;
+    std::vector<svmp::math::Vector<double, 3>> nodes = corners;
     for (const auto& edge : edges) {
-        const math::Vector<double, 3> midpoint =
+        const svmp::math::Vector<double, 3> midpoint =
             (corners[static_cast<std::size_t>(edge[0])] +
              corners[static_cast<std::size_t>(edge[1])]) * double(0.5);
         nodes.push_back(midpoint);
@@ -357,8 +357,8 @@ std::vector<math::Vector<double, 3>> wedge15_reference_nodes_for_test() {
 // Independent Quad8 reference layout: the four quad corners followed by the four
 // edge midpoints, in the corner/edge order the reference layout uses (the VTK
 // quad boundary traversal). Mirrors the Hex20/Wedge15 anchors above.
-std::vector<math::Vector<double, 3>> quad8_reference_nodes_for_test() {
-    std::vector<math::Vector<double, 3>> corners;
+std::vector<svmp::math::Vector<double, 3>> quad8_reference_nodes_for_test() {
+    std::vector<svmp::math::Vector<double, 3>> corners;
     corners.push_back({double(-1), double(-1), double(0)});
     corners.push_back({double(1), double(-1), double(0)});
     corners.push_back({double(1), double(1), double(0)});
@@ -366,9 +366,9 @@ std::vector<math::Vector<double, 3>> quad8_reference_nodes_for_test() {
 
     const int edges[4][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}};
 
-    std::vector<math::Vector<double, 3>> nodes = corners;
+    std::vector<svmp::math::Vector<double, 3>> nodes = corners;
     for (const auto& edge : edges) {
-        const math::Vector<double, 3> midpoint =
+        const svmp::math::Vector<double, 3> midpoint =
             (corners[static_cast<std::size_t>(edge[0])] +
              corners[static_cast<std::size_t>(edge[1])]) * double(0.5);
         nodes.push_back(midpoint);
@@ -394,7 +394,7 @@ double legendre_value_for_test(double x, int degree) {
     return p_k;
 }
 
-double legendre_mode_for_test(const math::Vector<double, 3>& p,
+double legendre_mode_for_test(const svmp::math::Vector<double, 3>& p,
                               const std::array<int, 3>& mode) {
     return legendre_value_for_test(p[0], mode[0]) *
            legendre_value_for_test(p[1], mode[1]) *
@@ -416,7 +416,7 @@ double matrix_norm_inf_for_test(const std::vector<double>& matrix, std::size_t n
 // Infinity-norm condition number of the Legendre generalized Vandermonde the
 // production basis inverts, rebuilt from the basis nodes and the re-derived
 // serendipity modes (an independent check that Source B is fixed).
-double legendre_vandermonde_condition(const std::vector<math::Vector<double, 3>>& nodes,
+double legendre_vandermonde_condition(const std::vector<svmp::math::Vector<double, 3>>& nodes,
                                       const std::vector<std::array<int, 3>>& modes) {
     const std::size_t n = nodes.size();
     std::vector<double> v(n * n, double(0));
@@ -426,7 +426,7 @@ double legendre_vandermonde_condition(const std::vector<math::Vector<double, 3>>
         }
     }
     const double norm_v = matrix_norm_inf_for_test(v, n);
-    const auto inverse = math::invert_dense_matrix(v, n, "test Legendre Vandermonde");
+    const auto inverse = svmp::math::invert_dense_matrix(v, n, "test Legendre Vandermonde");
     return norm_v * matrix_norm_inf_for_test(inverse, n);
 }
 
@@ -452,7 +452,7 @@ double serendipity_lebesgue_constant(const SerendipityBasis& basis, int samples)
         for (int j = 0; j <= samples; ++j) {
             const int kmax = (dim >= 3) ? samples : 0;
             for (int k = 0; k <= kmax; ++k) {
-                const math::Vector<double, 3> xi{
+                const svmp::math::Vector<double, 3> xi{
                     axis(i), axis(j), dim >= 3 ? axis(k) : double(0)};
                 basis.evaluate_values(xi, values);
                 double sum = double(0);
@@ -690,7 +690,7 @@ TEST(SerendipityBasis, QuadrilateralNodesFollowDocumentedConstructionThroughOrde
 TEST(SerendipityBasis, QuadrilateralOrderOneReproducesBilinearFunctions) {
     SerendipityBasis basis(BasisTopology::Quadrilateral, 1);
 
-    const std::vector<math::Vector<double, 3>> points = {
+    const std::vector<svmp::math::Vector<double, 3>> points = {
         {double(0.25), double(-0.4), double(0)},
         {double(-0.7), double(0.6), double(0)},
     };
@@ -704,7 +704,7 @@ TEST(SerendipityBasis, QuadrilateralOrderOneReproducesBilinearFunctions) {
 TEST(SerendipityBasis, QuadrilateralOrderThreeReproducesSerendipityCubics) {
     SerendipityBasis basis(BasisTopology::Quadrilateral, 3);
 
-    const std::vector<math::Vector<double, 3>> points = {
+    const std::vector<svmp::math::Vector<double, 3>> points = {
         {double(0.25), double(-0.4), double(0)},
         {double(-0.7), double(0.6), double(0)},
     };
@@ -716,7 +716,7 @@ TEST(SerendipityBasis, QuadrilateralOrderThreeReproducesSerendipityCubics) {
 }
 
 TEST(SerendipityBasis, QuadrilateralOrdersReproduceEverySerendipityMonomial) {
-    const std::vector<math::Vector<double, 3>> points = {
+    const std::vector<svmp::math::Vector<double, 3>> points = {
         {double(0.25), double(-0.4), double(0)},
         {double(-0.7), double(0.6), double(0)},
         {double(0.11), double(0.23), double(0)},
@@ -737,7 +737,7 @@ TEST(SerendipityBasis, QuadrilateralOrdersReproduceEverySerendipityMonomial) {
                     interpolate_nodal_function(
                         basis,
                         xi,
-                        [&exponent](const math::Vector<double, 3>& node) {
+                        [&exponent](const svmp::math::Vector<double, 3>& node) {
                             return monomial_value_for_test(node, exponent);
                         });
                 const double expected = monomial_value_for_test(xi, exponent);
@@ -760,7 +760,7 @@ TEST(SerendipityBasis, QuadrilateralVandermondeHasFullRankThroughOrderTen) {
 
         ASSERT_EQ(exponents.size(), n) << "order=" << order;
         ASSERT_EQ(vandermonde.size(), n * n) << "order=" << order;
-        EXPECT_EQ(math::dense_matrix_rank(vandermonde, n, n), n)
+        EXPECT_EQ(svmp::math::dense_matrix_rank(vandermonde, n, n), n)
             << "order=" << order;
     }
 }
@@ -779,7 +779,7 @@ TEST(SerendipityBasis, TrilinearHexMatchesLagrangeHex8) {
                        reference_nodes(ElementType::Hex8, serendipity.size()),
                        double(1e-12));
 
-    const std::vector<math::Vector<double, 3>> points = {
+    const std::vector<svmp::math::Vector<double, 3>> points = {
         {double(0.2), double(-0.1), double(0.3)},
         {double(-0.35), double(0.25), double(-0.15)},
     };
@@ -818,7 +818,7 @@ TEST(SerendipityBasis, Hex20VandermondeHasFullRank) {
     ASSERT_EQ(exponents.size(), n);
     const auto vandermonde = vandermonde_3d_for_test(basis.nodes(), exponents);
     ASSERT_EQ(vandermonde.size(), n * n);
-    EXPECT_EQ(math::dense_matrix_rank(vandermonde, n, n), n);
+    EXPECT_EQ(svmp::math::dense_matrix_rank(vandermonde, n, n), n);
 }
 
 TEST(SerendipityBasis, Wedge15VandermondeHasFullRank) {
@@ -828,7 +828,7 @@ TEST(SerendipityBasis, Wedge15VandermondeHasFullRank) {
     ASSERT_EQ(exponents.size(), n);
     const auto vandermonde = vandermonde_3d_for_test(basis.nodes(), exponents);
     ASSERT_EQ(vandermonde.size(), n * n);
-    EXPECT_EQ(math::dense_matrix_rank(vandermonde, n, n), n);
+    EXPECT_EQ(svmp::math::dense_matrix_rank(vandermonde, n, n), n);
 }
 
 // V * C == I guard: independently invert the Vandermonde and confirm the basis
@@ -841,9 +841,9 @@ TEST(SerendipityBasis, Hex20MatchesIndependentlyInvertedVandermonde) {
     ASSERT_EQ(exponents.size(), n);
     auto vandermonde = vandermonde_3d_for_test(basis.nodes(), exponents);
     const auto coefficients =
-        math::invert_dense_matrix(std::move(vandermonde), n, "Hex20 test Vandermonde");
+        svmp::math::invert_dense_matrix(std::move(vandermonde), n, "Hex20 test Vandermonde");
 
-    const std::vector<math::Vector<double, 3>> points = {
+    const std::vector<svmp::math::Vector<double, 3>> points = {
         {double(0.2), double(-0.1), double(0.3)},
         {double(-0.35), double(0.25), double(-0.15)},
         {double(0.11), double(0.23), double(-0.42)},
@@ -870,9 +870,9 @@ TEST(SerendipityBasis, Wedge15MatchesIndependentlyInvertedVandermonde) {
     ASSERT_EQ(exponents.size(), n);
     auto vandermonde = vandermonde_3d_for_test(basis.nodes(), exponents);
     const auto coefficients =
-        math::invert_dense_matrix(std::move(vandermonde), n, "Wedge15 test Vandermonde");
+        svmp::math::invert_dense_matrix(std::move(vandermonde), n, "Wedge15 test Vandermonde");
 
-    const std::vector<math::Vector<double, 3>> points = {
+    const std::vector<svmp::math::Vector<double, 3>> points = {
         {double(0.2), double(0.3), double(0.1)},
         {double(0.25), double(0.25), double(-0.4)},
         {double(0.1), double(0.6), double(0.5)},
@@ -899,7 +899,7 @@ TEST(SerendipityBasis, Hex20ReproducesEverySerendipityMonomial) {
     const auto exponents = hex20_serendipity_exponents_for_test();
     ASSERT_EQ(exponents.size(), basis.size());
 
-    const std::vector<math::Vector<double, 3>> points = {
+    const std::vector<svmp::math::Vector<double, 3>> points = {
         {double(0.2), double(-0.1), double(0.3)},
         {double(-0.35), double(0.25), double(-0.15)},
         {double(0.11), double(0.23), double(-0.42)},
@@ -908,7 +908,7 @@ TEST(SerendipityBasis, Hex20ReproducesEverySerendipityMonomial) {
         for (const auto& xi : points) {
             const double interpolated = interpolate_nodal_function(
                 basis, xi,
-                [&exponent](const math::Vector<double, 3>& node) {
+                [&exponent](const svmp::math::Vector<double, 3>& node) {
                     return monomial_value_3d_for_test(node, exponent);
                 });
             EXPECT_NEAR(interpolated, monomial_value_3d_for_test(xi, exponent),
@@ -924,7 +924,7 @@ TEST(SerendipityBasis, Wedge15ReproducesEverySerendipityMonomial) {
     const auto exponents = wedge15_serendipity_exponents_for_test();
     ASSERT_EQ(exponents.size(), basis.size());
 
-    const std::vector<math::Vector<double, 3>> points = {
+    const std::vector<svmp::math::Vector<double, 3>> points = {
         {double(0.2), double(0.3), double(0.1)},
         {double(0.25), double(0.25), double(-0.4)},
         {double(0.1), double(0.6), double(0.5)},
@@ -933,7 +933,7 @@ TEST(SerendipityBasis, Wedge15ReproducesEverySerendipityMonomial) {
         for (const auto& xi : points) {
             const double interpolated = interpolate_nodal_function(
                 basis, xi,
-                [&exponent](const math::Vector<double, 3>& node) {
+                [&exponent](const svmp::math::Vector<double, 3>& node) {
                     return monomial_value_3d_for_test(node, exponent);
                 });
             EXPECT_NEAR(interpolated, monomial_value_3d_for_test(xi, exponent),
@@ -1044,12 +1044,12 @@ TEST(SerendipityBasis, HexahedralSerendipityVandermondeHasFullRankThroughOrderSi
         ASSERT_EQ(exponents.size(), n) << "order=" << order;
         const auto vandermonde = vandermonde_3d_for_test(basis.nodes(), exponents);
         ASSERT_EQ(vandermonde.size(), n * n) << "order=" << order;
-        EXPECT_EQ(math::dense_matrix_rank(vandermonde, n, n), n) << "order=" << order;
+        EXPECT_EQ(svmp::math::dense_matrix_rank(vandermonde, n, n), n) << "order=" << order;
     }
 }
 
 TEST(SerendipityBasis, HexahedralTopologyIsNodalAndPartitionsUnity) {
-    const std::vector<math::Vector<double, 3>> points = {
+    const std::vector<svmp::math::Vector<double, 3>> points = {
         {double(0.2), double(-0.1), double(0.3)},
         {double(-0.35), double(0.25), double(-0.15)},
     };
@@ -1071,7 +1071,7 @@ TEST(SerendipityBasis, HexahedralTopologyIsNodalAndPartitionsUnity) {
 // monomial in its span at interior points, pinning the production monomial space
 // against the re-derived verification.
 TEST(SerendipityBasis, HexahedralTopologyReproducesEverySerendipityMonomial) {
-    const std::vector<math::Vector<double, 3>> points = {
+    const std::vector<svmp::math::Vector<double, 3>> points = {
         {double(0.2), double(-0.1), double(0.3)},
         {double(-0.35), double(0.25), double(-0.15)},
         {double(0.11), double(0.23), double(-0.42)},
@@ -1084,7 +1084,7 @@ TEST(SerendipityBasis, HexahedralTopologyReproducesEverySerendipityMonomial) {
             for (const auto& xi : points) {
                 const double interpolated = interpolate_nodal_function(
                     basis, xi,
-                    [&exponent](const math::Vector<double, 3>& node) {
+                    [&exponent](const svmp::math::Vector<double, 3>& node) {
                         return monomial_value_3d_for_test(node, exponent);
                     });
                 EXPECT_NEAR(interpolated, monomial_value_3d_for_test(xi, exponent),
@@ -1102,7 +1102,7 @@ TEST(SerendipityBasis, NamedHexLayoutsMatchTopologyConstruction) {
         {ElementType::Hex8, 1},
         {ElementType::Hex20, 2},
     };
-    const std::vector<math::Vector<double, 3>> points = {
+    const std::vector<svmp::math::Vector<double, 3>> points = {
         {double(0.2), double(-0.1), double(0.3)},
         {double(-0.35), double(0.25), double(-0.15)},
         {double(0.11), double(0.23), double(-0.42)},
diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp
index 6d2edf3c7..5055e7c66 100644
--- a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp
+++ b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp
@@ -60,7 +60,7 @@ void expect_invalid_argument_with_message(
 
 TEST(QuadPointContract, UsesFixedSizeFEVectorRepresentation)
 {
-    static_assert(std::is_same_v<QuadPoint, math::Vector<double, 3>>);
+    static_assert(std::is_same_v<QuadPoint, svmp::math::Vector<double, 3>>);
     static_assert(QuadPoint::RowsAtCompileTime == 3);
     static_assert(QuadPoint::ColsAtCompileTime == 1);
 
diff --git a/tests/unitTests/FE/Math/test_DenseLinearAlgebra.cpp b/tests/unitTests/common/math/test_DenseLinearAlgebra.cpp
similarity index 97%
rename from tests/unitTests/FE/Math/test_DenseLinearAlgebra.cpp
rename to tests/unitTests/common/math/test_DenseLinearAlgebra.cpp
index b029a55b4..672fff401 100644
--- a/tests/unitTests/FE/Math/test_DenseLinearAlgebra.cpp
+++ b/tests/unitTests/common/math/test_DenseLinearAlgebra.cpp
@@ -5,15 +5,15 @@
 
 #include <gtest/gtest.h>
 
-#include "FE/Common/FEException.h"
-#include "FE/Math/DenseLinearAlgebra.h"
+#include "Core/Exception.h"
+#include "math/DenseLinearAlgebra.h"
 
 #include <cmath>
 #include <span>
 #include <vector>
 
-using namespace svmp::FE;
-using namespace svmp::FE::math;
+using namespace svmp::math;
+using svmp::CoreException;
 
 namespace {
 
@@ -195,19 +195,19 @@ TEST(DenseLinearAlgebra, SolveInPlaceValidatesInputs) {
 
     std::vector<double> rhs{double(1), double(2)};
     EXPECT_THROW(solver.solve_in_place(std::span<double>(rhs.data(), rhs.size()), 0u),
-                 FEException);
+                 CoreException);
 
     std::vector<double> wrong_size{double(1), double(2), double(3)};
     EXPECT_THROW(
         solver.solve_in_place(std::span<double>(wrong_size.data(), wrong_size.size()), 1u),
-        FEException);
+        CoreException);
 
     DenseLUSolver unfactored;
     unfactored.n = 2u;
     unfactored.error_message_label = "unfactored";
     EXPECT_FALSE(unfactored.empty());
     EXPECT_THROW(unfactored.solve_in_place(std::span<double>(rhs.data(), rhs.size()), 1u),
-                 FEException);
+                 CoreException);
 }
 
 TEST(DenseLinearAlgebra, DiagnosticValidationRejectsRankMismatch) {
@@ -215,7 +215,7 @@ TEST(DenseLinearAlgebra, DiagnosticValidationRejectsRankMismatch) {
     result.diagnostics.rank = 1u;
 
     EXPECT_THROW(validate_dense_inverse_diagnostics(result, 2u, "rank mismatch"),
-                 FEException);
+                 CoreException);
 }
 
 TEST(DenseLinearAlgebra, RankHandlesNonSquareMatrices) {
@@ -266,7 +266,7 @@ TEST(DenseLinearAlgebra, DiagnosticValidationRejectsUnsupportedCondition) {
               dense_matrix_condition_error_threshold());
     EXPECT_THROW(validate_dense_inverse_diagnostics(
                      result, 2u, "excessive-condition diagonal"),
-                 FEException);
+                 CoreException);
 }
 
 TEST(DenseLinearAlgebra, ThrowsForScaleAwareSingularPivot) {
@@ -276,7 +276,7 @@ TEST(DenseLinearAlgebra, ThrowsForScaleAwareSingularPivot) {
     };
 
     EXPECT_THROW((void)invert_dense_matrix(singular, 2u, "singular 2x2"),
-                 FEException);
+                 CoreException);
 }
 
 TEST(DenseLinearAlgebra, FactorizationThrowsForRankDeficientMatrix) {
@@ -286,7 +286,7 @@ TEST(DenseLinearAlgebra, FactorizationThrowsForRankDeficientMatrix) {
     };
 
     EXPECT_THROW((void)factor_dense_matrix(singular, 2u, "rank-one 2x2"),
-                 FEException);
+                 CoreException);
 }
 
 TEST(DenseLinearAlgebra, RankUsesScaleAwareTolerance) {
diff --git a/tests/unitTests/linear_solver_tests/test_precond_rcs.cpp b/tests/unitTests/linear_solver_tests/test_precond_rcs.cpp
index 74c9ce2c2..28bd07699 100644
--- a/tests/unitTests/linear_solver_tests/test_precond_rcs.cpp
+++ b/tests/unitTests/linear_solver_tests/test_precond_rcs.cpp
@@ -40,7 +40,7 @@
  */
 
 #include "precond.h"
-#include "FE/Math/DenseLinearAlgebra.h"
+#include "math/DenseLinearAlgebra.h"
 
 #include <gtest/gtest.h>
 
@@ -152,7 +152,7 @@ TEST_P(RcsPreconditioner, PreservesConstrainedSystem) {
     }
   }
 
-  svmp::FE::math::factor_dense_matrix(matrix, n, "RCS regression").solve_in_place(solution);
+  svmp::math::factor_dense_matrix(matrix, n, "RCS regression").solve_in_place(solution);
   for (int row = 0; row < n; ++row) {
     solution[row] *= right(row % dof, row / dof);
     EXPECT_NEAR(solution[row], exact[row], 1e-11);