diff --git a/CMakeLists.txt b/CMakeLists.txt
index d4f572e..9f1c6fb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -51,7 +51,7 @@ target_include_directories(OpenABF
         $<BUILD_INTERFACE:${OPENABF_LINK_DIR}>
         $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
 )
-target_compile_features(OpenABF INTERFACE cxx_std_17)
+target_compile_features(OpenABF INTERFACE cxx_std_20)
 set_target_properties(OpenABF
     PROPERTIES
         PUBLIC_HEADER "${public_hdrs}"
diff --git a/README.md b/README.md
index d5fed8c..d01162c 100644
--- a/README.md
+++ b/README.md
@@ -5,7 +5,7 @@ The templated interface is designed for simple out-of-the-box use, and
 integration with existing geometric processing pipelines is quick and easy.
 
 ## Dependencies
-- C++17 compiler
+- C++20 compiler
 - [Eigen 3.3+](http://eigen.tuxfamily.org/)
 - CMake 3.15+ (optional)
 
@@ -161,7 +161,7 @@ path. As OpenABF depends upon the Eigen library, you will also need to add the
 Eigen headers to your include path:
 
 ```shell
-g++ -I /path/to/eigen/ -std=c++17 -DNDEBUG -O3 main.cpp -o main
+g++ -I /path/to/eigen/ -std=c++20 -DNDEBUG -O3 main.cpp -o main
 ```
 
 > [!IMPORTANT]
@@ -175,8 +175,8 @@ automatically conformant with the C++ standard in all cases. This may lead to
 the following issues when compiling against OpenABF.
 
 > [!NOTE]
-> As this project only supports C++17 and up, you should always compile 
-with at least `/std:c++17`.
+> As this project only supports C++20 and up, you should always compile 
+with at least `/std:c++20`.
 
 #### Undeclared identifier errors
 
diff --git a/conductor/code_styleguides/cpp.md b/conductor/code_styleguides/cpp.md
index f4218c3..0fcc005 100644
--- a/conductor/code_styleguides/cpp.md
+++ b/conductor/code_styleguides/cpp.md
@@ -42,7 +42,7 @@ Setters on instance configuration follow the instance-method rule: `set_<thing>(
 Pre-2026 code in several places uses `camelCase` for instance methods and `detail::` free functions (e.g., `setLevelRatio`, `detail::lscm::buildSystem`, `DecimationMesh::tryCollapse`). These are tracked for rename in issue #94. New code must follow the table above.
 
 ## C++ Standards
-- **Standard**: C++17 required; do not use C++20 features in public headers without a guard
+- **Standard**: C++20 required
 - **Headers**: All public headers must be self-contained (include what they use)
 - **`#pragma once`**: Preferred over include guards in project headers
 - **RAII**: Prefer RAII over manual resource management
diff --git a/conductor/setup_state.json b/conductor/setup_state.json
index 7dcef86..87bd669 100644
--- a/conductor/setup_state.json
+++ b/conductor/setup_state.json
@@ -12,7 +12,7 @@
     "key_goals": "Zero-to-low dependency header-only integration; multiple algorithms (ABF, ABF++, LSCM); multi-chart flattening and packing; high numerical accuracy; production-efficient implementations.",
     "voice_tone": "Concise and direct",
     "design_principles": "Simplicity over features, Performance first, Developer experience focused",
-    "languages": "C++17, Python 3.x (scripts)",
+    "languages": "C++20, Python 3.x (scripts)",
     "frontend": "None",
     "backend": "None",
     "database": "None",
diff --git a/conductor/tech-stack.md b/conductor/tech-stack.md
index 49deac9..87e86d1 100644
--- a/conductor/tech-stack.md
+++ b/conductor/tech-stack.md
@@ -1,7 +1,7 @@
 # Tech Stack
 
 ## Primary Language
-- **C++17** — all library code; required standard for consumers
+- **C++20** — all library code; required standard for consumers
 - **Python 3.x** — utility/build scripts (e.g. `thirdparty/amalgamate/amalgamate.py` for single-header generation); Python bindings (pybind11) are a future consideration
 
 ## Build System
diff --git a/docs/installation.md b/docs/installation.md
index eceee85..51c1e53 100644
--- a/docs/installation.md
+++ b/docs/installation.md
@@ -5,7 +5,7 @@ file(s) into your include directory. This page outlines the various installation
 options.
 
 ## Dependencies
-- C++17 compiler
+- C++20 compiler
 - [Eigen 3.3+](http://eigen.tuxfamily.org/)
 - CMake 3.15+ (optional)
 
@@ -74,7 +74,7 @@ path. As %OpenABF depends upon the Eigen library, you will also need to add the
 Eigen headers to your include path:
 
 ```{.sh}
-g++ -I /path/to/eigen/ -std=c++17 -DNDEBUG -O3 main.cpp -o main
+g++ -I /path/to/eigen/ -std=c++20 -DNDEBUG -O3 main.cpp -o main
 ```
 
 **Note:** For best performance, compile your application with the `-DNDEBUG -03`
@@ -86,8 +86,8 @@ For many legacy reasons, the Microsoft Visual C++ compiler (MSVC) is not
 automatically conformant with the C++ standard in all cases. This may lead to
 the following issues when compiling against OpenABF.
 
-**Note:** As this project only supports C++17 and up, you should always compile
-with at least `/std:c++17`.
+**Note:** As this project only supports C++20 and up, you should always compile
+with at least `/std:c++20`.
 
 ### Undeclared identifier errors
 
diff --git a/include/OpenABF/ABF.hpp b/include/OpenABF/ABF.hpp
index 98f5e27..01a23c6 100644
--- a/include/OpenABF/ABF.hpp
+++ b/include/OpenABF/ABF.hpp
@@ -2,6 +2,7 @@
 
 #include <cassert>
 #include <cmath>
+#include <concepts>
 #include <limits>
 #include <vector>
 
@@ -84,7 +85,7 @@ void InitializeAnglesAndWeights(MeshPtr& m)
 }
 
 /** @brief Compute ∇CTri w.r.t LambdaTri == CTri */
-template <typename T, class FacePtr, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class FacePtr>
 auto TriGrad(const FacePtr& f) -> T
 {
     T g = -PI<T>;
@@ -95,7 +96,7 @@ auto TriGrad(const FacePtr& f) -> T
 }
 
 /** @brief Compute ∇CPlan w.r.t LambdaPlan == CPlan */
-template <typename T, class VertPtr, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class VertPtr>
 auto PlanGrad(const VertPtr& v) -> T
 {
     T g = -2 * PI<T>;
@@ -106,7 +107,7 @@ auto PlanGrad(const VertPtr& v) -> T
 }
 
 /** @brief Compute ∇CLen w.r.t LambdaLen == CLen */
-template <typename T, class VertPtr, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class VertPtr>
 auto LenGrad(const VertPtr& vertex) -> T
 {
     T p1{1};
@@ -119,8 +120,7 @@ auto LenGrad(const VertPtr& vertex) -> T
 }
 
 /** @brief Compute ∇CLen w.r.t edge->alpha */
-template <typename T, class VertPtr, class EdgePtr,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class VertPtr, class EdgePtr>
 auto LenGrad(const VertPtr& vertex, const EdgePtr& edge) -> T
 {
     T p1{1};
@@ -146,7 +146,7 @@ auto LenGrad(const VertPtr& vertex, const EdgePtr& edge) -> T
 }
 
 /** @brief Compute ∇F w.r.t an edge's alpha */
-template <typename T, class EdgePtr, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class EdgePtr>
 auto AlphaGrad(const EdgePtr& edge) -> T
 {
     // δE/δα
@@ -172,7 +172,7 @@ auto AlphaGrad(const EdgePtr& edge) -> T
 }
 
 /** @brief Compute ∇F w.r.t all parameters */
-template <typename T, class MeshPtr, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class MeshPtr>
 auto Gradient(const MeshPtr& mesh) -> T
 {
     T g{0};
@@ -225,9 +225,8 @@ auto Gradient(const MeshPtr& mesh) -> T
  * concept](https://eigen.tuxfamily.org/dox-devel/group__TopicSparseSystems.html)
  * and templated on Eigen::SparseMatrix<T>
  */
-template <typename T, class MeshType = detail::ABF::Mesh<T>,
-          class Solver = Eigen::SparseLU<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int>>,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class MeshType = detail::ABF::Mesh<T>,
+          class Solver = Eigen::SparseLU<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int>>>
 class ABF
 {
 public:
diff --git a/include/OpenABF/ABFPlusPlus.hpp b/include/OpenABF/ABFPlusPlus.hpp
index 6999d9a..59050a0 100644
--- a/include/OpenABF/ABFPlusPlus.hpp
+++ b/include/OpenABF/ABFPlusPlus.hpp
@@ -2,6 +2,7 @@
 
 #include <cassert>
 #include <cmath>
+#include <concepts>
 #include <limits>
 #include <vector>
 
@@ -38,9 +39,8 @@ namespace OpenABF
  * concept](https://eigen.tuxfamily.org/dox-devel/group__TopicSparseSystems.html)
  * and templated on Eigen::SparseMatrix<T>
  */
-template <typename T, class MeshType = detail::ABF::Mesh<T>,
-          class Solver = Eigen::SparseLU<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int>>,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class MeshType = detail::ABF::Mesh<T>,
+          class Solver = Eigen::SparseLU<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int>>>
 class ABFPlusPlus
 {
 public:
diff --git a/include/OpenABF/AngleBasedLSCM.hpp b/include/OpenABF/AngleBasedLSCM.hpp
index 1ce626a..47c071c 100644
--- a/include/OpenABF/AngleBasedLSCM.hpp
+++ b/include/OpenABF/AngleBasedLSCM.hpp
@@ -2,6 +2,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <concepts>
 #include <cstddef>
 #include <iterator>
 #include <optional>
@@ -32,9 +33,8 @@ template <template <class...> class U, class... Vs>
 constexpr bool is_instance_of_v<U<Vs...>, U> = std::true_type{};
 
 /** Solve least squares using A'Ab  */
-template <
-    class SparseMatrix, class DenseMatrix, class Solver,
-    std::enable_if_t<!is_instance_of_v<Solver, Eigen::LeastSquaresConjugateGradient>, bool> = false>
+template <class SparseMatrix, class DenseMatrix, class Solver>
+    requires(!is_instance_of_v<Solver, Eigen::LeastSquaresConjugateGradient>)
 auto SolveLeastSquares(SparseMatrix A, SparseMatrix b) -> DenseMatrix
 {
     // Setup AtA and solver
@@ -56,9 +56,8 @@ auto SolveLeastSquares(SparseMatrix A, SparseMatrix b) -> DenseMatrix
 }
 
 /** Solve least squares with LeastSquaresConjugateGradient */
-template <
-    class SparseMatrix, class DenseMatrix, class Solver,
-    std::enable_if_t<is_instance_of_v<Solver, Eigen::LeastSquaresConjugateGradient>, bool> = true>
+template <class SparseMatrix, class DenseMatrix, class Solver>
+    requires(is_instance_of_v<Solver, Eigen::LeastSquaresConjugateGradient>)
 auto SolveLeastSquares(SparseMatrix A, SparseMatrix b) -> DenseMatrix
 {
     // Solve
@@ -118,9 +117,8 @@ auto SolveLeastSquares(SparseMatrix A, SparseMatrix b) -> DenseMatrix
  * at the mesh sizes we target. Use IC only if you have profiled it favorably
  * against Diagonal for your specific mesh class.
  */
-template <typename T, class MeshType = HalfEdgeMesh<T>,
-          class Solver = Eigen::SparseLU<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int>>,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class MeshType = HalfEdgeMesh<T>,
+          class Solver = Eigen::SparseLU<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int>>>
 class AngleBasedLSCM
 {
 public:
diff --git a/include/OpenABF/HalfEdgeMesh.hpp b/include/OpenABF/HalfEdgeMesh.hpp
index 5e58baa..9738490 100644
--- a/include/OpenABF/HalfEdgeMesh.hpp
+++ b/include/OpenABF/HalfEdgeMesh.hpp
@@ -430,14 +430,16 @@ class HalfEdgeMesh
 
         /** Dereference operator */
         template <bool Const_ = Const>
-        std::enable_if_t<Const_, reference> operator*() const
+            requires(Const_)
+        reference operator*() const
         {
             return current_;
         }
 
         /** Dereference operator */
         template <bool Const_ = Const>
-        std::enable_if_t<not Const_, reference> operator*()
+            requires(not Const_)
+        reference operator*()
         {
             return current_;
         }
@@ -507,13 +509,15 @@ class HalfEdgeMesh
 
         /** Dereference (const overload) */
         template <bool C = Const>
-        auto operator*() const -> std::enable_if_t<C, reference>
+            requires(C)
+        auto operator*() const -> reference
         {
             return current_;
         }
         /** Dereference (non-const overload) */
         template <bool C = Const>
-        auto operator*() -> std::enable_if_t<!C, reference>
+            requires(!C)
+        auto operator*() -> reference
         {
             return current_;
         }
@@ -598,13 +602,15 @@ class HalfEdgeMesh
 
         /** Dereference (const overload) */
         template <bool C = Const>
-        auto operator*() const -> std::enable_if_t<C, reference>
+            requires(C)
+        auto operator*() const -> reference
         {
             return *edgeIt_;
         }
         /** Dereference (non-const overload) */
         template <bool C = Const>
-        auto operator*() -> std::enable_if_t<!C, reference>
+            requires(!C)
+        auto operator*() -> reference
         {
             return *edgeIt_;
         }
diff --git a/include/OpenABF/HierarchicalLSCM.hpp b/include/OpenABF/HierarchicalLSCM.hpp
index 55b2ad5..5c48c18 100644
--- a/include/OpenABF/HierarchicalLSCM.hpp
+++ b/include/OpenABF/HierarchicalLSCM.hpp
@@ -3,6 +3,7 @@
 #include <algorithm>
 #include <array>
 #include <cmath>
+#include <concepts>
 #include <limits>
 #include <numeric>
 #include <optional>
@@ -969,10 +970,9 @@ auto SolveLSCMLevel(const typename HalfEdgeMesh<T>::Pointer& levelMesh,
  *       The result is numerically equivalent but may differ in convergence
  *       behavior from a plain AngleBasedLSCM call.
  */
-template <typename T, class MeshType = HalfEdgeMesh<T>,
+template <std::floating_point T, class MeshType = HalfEdgeMesh<T>,
           class Solver =
-              Eigen::ConjugateGradient<Eigen::SparseMatrix<T>, Eigen::Lower | Eigen::Upper>,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+              Eigen::ConjugateGradient<Eigen::SparseMatrix<T>, Eigen::Lower | Eigen::Upper>>
 class HierarchicalLSCM
 {
 public:
diff --git a/include/OpenABF/Math.hpp b/include/OpenABF/Math.hpp
index 0a75d05..6c8692a 100644
--- a/include/OpenABF/Math.hpp
+++ b/include/OpenABF/Math.hpp
@@ -2,6 +2,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <concepts>
 #include <limits>
 #include <numeric>
 
@@ -80,16 +81,14 @@ auto interior_angle(const Vector1& a, const Vector2& b)
 }
 
 /** @brief Convert degrees to radians */
-template <typename T = float, typename T2,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T = float, typename T2>
 constexpr auto to_radians(T2 deg) -> T
 {
     return deg * PI<T> / T(180);
 }
 
 /** @brief Convert radians to degrees */
-template <typename T = float, typename T2,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T = float, typename T2>
 constexpr auto to_degrees(T2 rad) -> T
 {
     return rad * T(180) / PI<T>;
diff --git a/include/OpenABF/Vec.hpp b/include/OpenABF/Vec.hpp
index c40ac8d..1643829 100644
--- a/include/OpenABF/Vec.hpp
+++ b/include/OpenABF/Vec.hpp
@@ -16,7 +16,8 @@ namespace OpenABF
  * @tparam T Element type
  * @tparam Dims Number of elements
  */
-template <typename T, std::size_t Dims, std::enable_if_t<std::is_arithmetic<T>::value, bool> = true>
+template <typename T, std::size_t Dims>
+    requires(std::is_arithmetic_v<T>)
 class Vec
 {
     /** Underlying element storage */
@@ -203,7 +204,8 @@ class Vec
     }
 
     /** @brief Multiplication assignment operator */
-    template <typename T2, std::enable_if_t<std::is_arithmetic<T2>::value, bool> = true>
+    template <typename T2>
+        requires(std::is_arithmetic_v<T2>)
     Vec& operator*=(const T2& b)
     {
         for (auto& v : val_) {
@@ -213,7 +215,8 @@ class Vec
     }
 
     /** @brief Multiplication operator */
-    template <typename T2, std::enable_if_t<std::is_arithmetic<T2>::value, bool> = true>
+    template <typename T2>
+        requires(std::is_arithmetic_v<T2>)
     friend Vec operator*(Vec lhs, const T2& rhs)
     {
         lhs *= rhs;
@@ -221,7 +224,8 @@ class Vec
     }
 
     /** @brief Division assignment operator */
-    template <typename T2, std::enable_if_t<std::is_arithmetic<T2>::value, bool> = true>
+    template <typename T2>
+        requires(std::is_arithmetic_v<T2>)
     Vec& operator/=(const T2& b)
     {
         for (auto& v : val_) {
@@ -231,7 +235,8 @@ class Vec
     }
 
     /** @brief Division operator */
-    template <typename T2, std::enable_if_t<std::is_arithmetic<T2>::value, bool> = true>
+    template <typename T2>
+        requires(std::is_arithmetic_v<T2>)
     friend Vec operator/(Vec lhs, const T2& rhs)
     {
         lhs /= rhs;
@@ -247,7 +252,8 @@ class Vec
 
     /** @brief Compute the vector cross product */
     template <class Vector, std::size_t D = Dims>
-    std::enable_if_t<D == 3, Vec> cross(const Vector& v)
+        requires(D == 3)
+    Vec cross(const Vector& v)
     {
         return OpenABF::cross(*this, v);
     }
diff --git a/single_include/OpenABF/OpenABF.hpp b/single_include/OpenABF/OpenABF.hpp
index b21d488..81213a5 100644
--- a/single_include/OpenABF/OpenABF.hpp
+++ b/single_include/OpenABF/OpenABF.hpp
@@ -54,6 +54,7 @@ class MeshException : public std::runtime_error
 
 #include <algorithm>
 #include <cmath>
+#include <concepts>
 #include <limits>
 #include <numeric>
 
@@ -132,16 +133,14 @@ auto interior_angle(const Vector1& a, const Vector2& b)
 }
 
 /** @brief Convert degrees to radians */
-template <typename T = float, typename T2,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T = float, typename T2>
 constexpr auto to_radians(T2 deg) -> T
 {
     return deg * PI<T> / T(180);
 }
 
 /** @brief Convert radians to degrees */
-template <typename T = float, typename T2,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T = float, typename T2>
 constexpr auto to_degrees(T2 rad) -> T
 {
     return rad * T(180) / PI<T>;
@@ -168,7 +167,8 @@ namespace OpenABF
  * @tparam T Element type
  * @tparam Dims Number of elements
  */
-template <typename T, std::size_t Dims, std::enable_if_t<std::is_arithmetic<T>::value, bool> = true>
+template <typename T, std::size_t Dims>
+    requires(std::is_arithmetic_v<T>)
 class Vec
 {
     /** Underlying element storage */
@@ -355,7 +355,8 @@ class Vec
     }
 
     /** @brief Multiplication assignment operator */
-    template <typename T2, std::enable_if_t<std::is_arithmetic<T2>::value, bool> = true>
+    template <typename T2>
+        requires(std::is_arithmetic_v<T2>)
     Vec& operator*=(const T2& b)
     {
         for (auto& v : val_) {
@@ -365,7 +366,8 @@ class Vec
     }
 
     /** @brief Multiplication operator */
-    template <typename T2, std::enable_if_t<std::is_arithmetic<T2>::value, bool> = true>
+    template <typename T2>
+        requires(std::is_arithmetic_v<T2>)
     friend Vec operator*(Vec lhs, const T2& rhs)
     {
         lhs *= rhs;
@@ -373,7 +375,8 @@ class Vec
     }
 
     /** @brief Division assignment operator */
-    template <typename T2, std::enable_if_t<std::is_arithmetic<T2>::value, bool> = true>
+    template <typename T2>
+        requires(std::is_arithmetic_v<T2>)
     Vec& operator/=(const T2& b)
     {
         for (auto& v : val_) {
@@ -383,7 +386,8 @@ class Vec
     }
 
     /** @brief Division operator */
-    template <typename T2, std::enable_if_t<std::is_arithmetic<T2>::value, bool> = true>
+    template <typename T2>
+        requires(std::is_arithmetic_v<T2>)
     friend Vec operator/(Vec lhs, const T2& rhs)
     {
         lhs /= rhs;
@@ -399,7 +403,8 @@ class Vec
 
     /** @brief Compute the vector cross product */
     template <class Vector, std::size_t D = Dims>
-    std::enable_if_t<D == 3, Vec> cross(const Vector& v)
+        requires(D == 3)
+    Vec cross(const Vector& v)
     {
         return OpenABF::cross(*this, v);
     }
@@ -874,14 +879,16 @@ class HalfEdgeMesh
 
         /** Dereference operator */
         template <bool Const_ = Const>
-        std::enable_if_t<Const_, reference> operator*() const
+            requires(Const_)
+        reference operator*() const
         {
             return current_;
         }
 
         /** Dereference operator */
         template <bool Const_ = Const>
-        std::enable_if_t<not Const_, reference> operator*()
+            requires(not Const_)
+        reference operator*()
         {
             return current_;
         }
@@ -951,13 +958,15 @@ class HalfEdgeMesh
 
         /** Dereference (const overload) */
         template <bool C = Const>
-        auto operator*() const -> std::enable_if_t<C, reference>
+            requires(C)
+        auto operator*() const -> reference
         {
             return current_;
         }
         /** Dereference (non-const overload) */
         template <bool C = Const>
-        auto operator*() -> std::enable_if_t<!C, reference>
+            requires(!C)
+        auto operator*() -> reference
         {
             return current_;
         }
@@ -1042,13 +1051,15 @@ class HalfEdgeMesh
 
         /** Dereference (const overload) */
         template <bool C = Const>
-        auto operator*() const -> std::enable_if_t<C, reference>
+            requires(C)
+        auto operator*() const -> reference
         {
             return *edgeIt_;
         }
         /** Dereference (non-const overload) */
         template <bool C = Const>
-        auto operator*() -> std::enable_if_t<!C, reference>
+            requires(!C)
+        auto operator*() -> reference
         {
             return *edgeIt_;
         }
@@ -2307,6 +2318,7 @@ class HalfEdgeMesh
 
 #include <cassert>
 #include <cmath>
+#include <concepts>
 #include <limits>
 #include <vector>
 
@@ -2392,7 +2404,7 @@ void InitializeAnglesAndWeights(MeshPtr& m)
 }
 
 /** @brief Compute ∇CTri w.r.t LambdaTri == CTri */
-template <typename T, class FacePtr, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class FacePtr>
 auto TriGrad(const FacePtr& f) -> T
 {
     T g = -PI<T>;
@@ -2403,7 +2415,7 @@ auto TriGrad(const FacePtr& f) -> T
 }
 
 /** @brief Compute ∇CPlan w.r.t LambdaPlan == CPlan */
-template <typename T, class VertPtr, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class VertPtr>
 auto PlanGrad(const VertPtr& v) -> T
 {
     T g = -2 * PI<T>;
@@ -2414,7 +2426,7 @@ auto PlanGrad(const VertPtr& v) -> T
 }
 
 /** @brief Compute ∇CLen w.r.t LambdaLen == CLen */
-template <typename T, class VertPtr, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class VertPtr>
 auto LenGrad(const VertPtr& vertex) -> T
 {
     T p1{1};
@@ -2427,8 +2439,7 @@ auto LenGrad(const VertPtr& vertex) -> T
 }
 
 /** @brief Compute ∇CLen w.r.t edge->alpha */
-template <typename T, class VertPtr, class EdgePtr,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class VertPtr, class EdgePtr>
 auto LenGrad(const VertPtr& vertex, const EdgePtr& edge) -> T
 {
     T p1{1};
@@ -2454,7 +2465,7 @@ auto LenGrad(const VertPtr& vertex, const EdgePtr& edge) -> T
 }
 
 /** @brief Compute ∇F w.r.t an edge's alpha */
-template <typename T, class EdgePtr, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class EdgePtr>
 auto AlphaGrad(const EdgePtr& edge) -> T
 {
     // δE/δα
@@ -2480,7 +2491,7 @@ auto AlphaGrad(const EdgePtr& edge) -> T
 }
 
 /** @brief Compute ∇F w.r.t all parameters */
-template <typename T, class MeshPtr, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class MeshPtr>
 auto Gradient(const MeshPtr& mesh) -> T
 {
     T g{0};
@@ -2533,9 +2544,8 @@ auto Gradient(const MeshPtr& mesh) -> T
  * concept](https://eigen.tuxfamily.org/dox-devel/group__TopicSparseSystems.html)
  * and templated on Eigen::SparseMatrix<T>
  */
-template <typename T, class MeshType = detail::ABF::Mesh<T>,
-          class Solver = Eigen::SparseLU<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int>>,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class MeshType = detail::ABF::Mesh<T>,
+          class Solver = Eigen::SparseLU<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int>>>
 class ABF
 {
 public:
@@ -2769,6 +2779,7 @@ class ABF
 
 #include <cassert>
 #include <cmath>
+#include <concepts>
 #include <limits>
 #include <vector>
 
@@ -2809,9 +2820,8 @@ namespace OpenABF
  * concept](https://eigen.tuxfamily.org/dox-devel/group__TopicSparseSystems.html)
  * and templated on Eigen::SparseMatrix<T>
  */
-template <typename T, class MeshType = detail::ABF::Mesh<T>,
-          class Solver = Eigen::SparseLU<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int>>,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class MeshType = detail::ABF::Mesh<T>,
+          class Solver = Eigen::SparseLU<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int>>>
 class ABFPlusPlus
 {
 public:
@@ -3069,6 +3079,7 @@ class ABFPlusPlus
 
 #include <algorithm>
 #include <cmath>
+#include <concepts>
 #include <cstddef>
 #include <iterator>
 #include <optional>
@@ -3391,9 +3402,8 @@ template <template <class...> class U, class... Vs>
 constexpr bool is_instance_of_v<U<Vs...>, U> = std::true_type{};
 
 /** Solve least squares using A'Ab  */
-template <
-    class SparseMatrix, class DenseMatrix, class Solver,
-    std::enable_if_t<!is_instance_of_v<Solver, Eigen::LeastSquaresConjugateGradient>, bool> = false>
+template <class SparseMatrix, class DenseMatrix, class Solver>
+    requires(!is_instance_of_v<Solver, Eigen::LeastSquaresConjugateGradient>)
 auto SolveLeastSquares(SparseMatrix A, SparseMatrix b) -> DenseMatrix
 {
     // Setup AtA and solver
@@ -3415,9 +3425,8 @@ auto SolveLeastSquares(SparseMatrix A, SparseMatrix b) -> DenseMatrix
 }
 
 /** Solve least squares with LeastSquaresConjugateGradient */
-template <
-    class SparseMatrix, class DenseMatrix, class Solver,
-    std::enable_if_t<is_instance_of_v<Solver, Eigen::LeastSquaresConjugateGradient>, bool> = true>
+template <class SparseMatrix, class DenseMatrix, class Solver>
+    requires(is_instance_of_v<Solver, Eigen::LeastSquaresConjugateGradient>)
 auto SolveLeastSquares(SparseMatrix A, SparseMatrix b) -> DenseMatrix
 {
     // Solve
@@ -3477,9 +3486,8 @@ auto SolveLeastSquares(SparseMatrix A, SparseMatrix b) -> DenseMatrix
  * at the mesh sizes we target. Use IC only if you have profiled it favorably
  * against Diagonal for your specific mesh class.
  */
-template <typename T, class MeshType = HalfEdgeMesh<T>,
-          class Solver = Eigen::SparseLU<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int>>,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+template <std::floating_point T, class MeshType = HalfEdgeMesh<T>,
+          class Solver = Eigen::SparseLU<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int>>>
 class AngleBasedLSCM
 {
 public:
@@ -3635,6 +3643,7 @@ class AngleBasedLSCM
 #include <algorithm>
 #include <array>
 #include <cmath>
+#include <concepts>
 #include <limits>
 #include <numeric>
 #include <optional>
@@ -4605,10 +4614,9 @@ auto SolveLSCMLevel(const typename HalfEdgeMesh<T>::Pointer& levelMesh,
  *       The result is numerically equivalent but may differ in convergence
  *       behavior from a plain AngleBasedLSCM call.
  */
-template <typename T, class MeshType = HalfEdgeMesh<T>,
+template <std::floating_point T, class MeshType = HalfEdgeMesh<T>,
           class Solver =
-              Eigen::ConjugateGradient<Eigen::SparseMatrix<T>, Eigen::Lower | Eigen::Upper>,
-          std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
+              Eigen::ConjugateGradient<Eigen::SparseMatrix<T>, Eigen::Lower | Eigen::Upper>>
 class HierarchicalLSCM
 {
 public: