diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index c4e98a14929b4daaad9b5fcf1b65a1b400f9318f..f4c66232816b33d8ee7e8ac897436507181d5644 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -282,7 +282,7 @@ class [[nodiscard]] TinyMatrix
   constexpr TinyMatrix& operator=(TinyMatrix&& A) noexcept = default;
 
   template <typename... Args>
-  PUGS_INLINE constexpr TinyMatrix(const T& t, Args&&... args) noexcept
+  PUGS_INLINE explicit constexpr TinyMatrix(const T& t, Args&&... args) noexcept
   {
     static_assert(sizeof...(args) == M * N - 1, "wrong number of parameters");
     this->_unpackVariadicInput(t, std::forward<Args>(args)...);
diff --git a/src/algebra/TinyVector.hpp b/src/algebra/TinyVector.hpp
index 029b259f568866208919f5864980fdd3079dc975..3ca3d2a38550b5a31bb16843d25a35cebdb10fd0 100644
--- a/src/algebra/TinyVector.hpp
+++ b/src/algebra/TinyVector.hpp
@@ -203,7 +203,7 @@ class [[nodiscard]] TinyVector
   constexpr TinyVector& operator=(TinyVector&& v) noexcept = default;
 
   template <typename... Args>
-  PUGS_INLINE constexpr TinyVector(const T& t, Args&&... args) noexcept
+  PUGS_INLINE explicit constexpr TinyVector(const T& t, Args&&... args) noexcept
   {
     static_assert(sizeof...(args) == N - 1, "wrong number of parameters");
     this->_unpackVariadicInput(t, std::forward<Args>(args)...);
diff --git a/src/analysis/CubeGaussQuadrature.cpp b/src/analysis/CubeGaussQuadrature.cpp
index 1d35bd875ee46c0faf5d492328c34d9ec7e6466c..f1aaf9fd23dfdd55cd5ac3eced298a2b83101204 100644
--- a/src/analysis/CubeGaussQuadrature.cpp
+++ b/src/analysis/CubeGaussQuadrature.cpp
@@ -25,7 +25,7 @@ CubeGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
       switch (id) {
       case 1: {
         Assert(position_list.size() == 0);
-        point_list[k]  = {0, 0, 0};
+        point_list[k]  = R3{0, 0, 0};
         weight_list[k] = w;
         ++k;
         break;
@@ -34,12 +34,12 @@ CubeGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
         Assert(position_list.size() == 1);
         const double a = position_list[0];
 
-        point_list[k + 0] = {+a, 0, 0};
-        point_list[k + 1] = {-a, 0, 0};
-        point_list[k + 2] = {0, +a, 0};
-        point_list[k + 3] = {0, -a, 0};
-        point_list[k + 4] = {0, 0, +a};
-        point_list[k + 5] = {0, 0, -a};
+        point_list[k + 0] = R3{+a, 0, 0};
+        point_list[k + 1] = R3{-a, 0, 0};
+        point_list[k + 2] = R3{0, +a, 0};
+        point_list[k + 3] = R3{0, -a, 0};
+        point_list[k + 4] = R3{0, 0, +a};
+        point_list[k + 5] = R3{0, 0, -a};
 
         for (size_t l = 0; l < 6; ++l) {
           weight_list[k + l] = w;
@@ -52,14 +52,14 @@ CubeGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
         Assert(position_list.size() == 1);
         const double a = position_list[0];
 
-        point_list[k + 0] = {+a, +a, +a};
-        point_list[k + 1] = {+a, +a, -a};
-        point_list[k + 2] = {+a, -a, +a};
-        point_list[k + 3] = {+a, -a, -a};
-        point_list[k + 4] = {-a, +a, +a};
-        point_list[k + 5] = {-a, +a, -a};
-        point_list[k + 6] = {-a, -a, +a};
-        point_list[k + 7] = {-a, -a, -a};
+        point_list[k + 0] = R3{+a, +a, +a};
+        point_list[k + 1] = R3{+a, +a, -a};
+        point_list[k + 2] = R3{+a, -a, +a};
+        point_list[k + 3] = R3{+a, -a, -a};
+        point_list[k + 4] = R3{-a, +a, +a};
+        point_list[k + 5] = R3{-a, +a, -a};
+        point_list[k + 6] = R3{-a, -a, +a};
+        point_list[k + 7] = R3{-a, -a, -a};
 
         for (size_t l = 0; l < 8; ++l) {
           weight_list[k + l] = w;
@@ -72,18 +72,18 @@ CubeGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
         Assert(position_list.size() == 1);
         const double a = position_list[0];
 
-        point_list[k + 0]  = {+a, +a, 0};
-        point_list[k + 1]  = {+a, -a, 0};
-        point_list[k + 2]  = {-a, +a, 0};
-        point_list[k + 3]  = {-a, -a, 0};
-        point_list[k + 4]  = {+a, 0, +a};
-        point_list[k + 5]  = {+a, 0, -a};
-        point_list[k + 6]  = {-a, 0, +a};
-        point_list[k + 7]  = {-a, 0, -a};
-        point_list[k + 8]  = {0, +a, +a};
-        point_list[k + 9]  = {0, +a, -a};
-        point_list[k + 10] = {0, -a, +a};
-        point_list[k + 11] = {0, -a, -a};
+        point_list[k + 0]  = R3{+a, +a, 0};
+        point_list[k + 1]  = R3{+a, -a, 0};
+        point_list[k + 2]  = R3{-a, +a, 0};
+        point_list[k + 3]  = R3{-a, -a, 0};
+        point_list[k + 4]  = R3{+a, 0, +a};
+        point_list[k + 5]  = R3{+a, 0, -a};
+        point_list[k + 6]  = R3{-a, 0, +a};
+        point_list[k + 7]  = R3{-a, 0, -a};
+        point_list[k + 8]  = R3{0, +a, +a};
+        point_list[k + 9]  = R3{0, +a, -a};
+        point_list[k + 10] = R3{0, -a, +a};
+        point_list[k + 11] = R3{0, -a, -a};
 
         for (size_t l = 0; l < 12; ++l) {
           weight_list[k + l] = w;
@@ -97,30 +97,30 @@ CubeGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
         const double a = position_list[0];
         const double b = position_list[1];
 
-        point_list[k + 0]  = {+a, +b, 0};
-        point_list[k + 1]  = {+a, -b, 0};
-        point_list[k + 2]  = {-a, +b, 0};
-        point_list[k + 3]  = {-a, -b, 0};
-        point_list[k + 4]  = {+b, +a, 0};
-        point_list[k + 5]  = {-b, +a, 0};
-        point_list[k + 6]  = {+b, -a, 0};
-        point_list[k + 7]  = {-b, -a, 0};
-        point_list[k + 8]  = {+a, 0, +b};
-        point_list[k + 9]  = {+a, 0, -b};
-        point_list[k + 10] = {-a, 0, +b};
-        point_list[k + 11] = {-a, 0, -b};
-        point_list[k + 12] = {+b, 0, +a};
-        point_list[k + 13] = {-b, 0, +a};
-        point_list[k + 14] = {+b, 0, -a};
-        point_list[k + 15] = {-b, 0, -a};
-        point_list[k + 16] = {0, +a, +b};
-        point_list[k + 17] = {0, +a, -b};
-        point_list[k + 18] = {0, -a, +b};
-        point_list[k + 19] = {0, -a, -b};
-        point_list[k + 20] = {0, +b, +a};
-        point_list[k + 21] = {0, -b, +a};
-        point_list[k + 22] = {0, +b, -a};
-        point_list[k + 23] = {0, -b, -a};
+        point_list[k + 0]  = R3{+a, +b, 0};
+        point_list[k + 1]  = R3{+a, -b, 0};
+        point_list[k + 2]  = R3{-a, +b, 0};
+        point_list[k + 3]  = R3{-a, -b, 0};
+        point_list[k + 4]  = R3{+b, +a, 0};
+        point_list[k + 5]  = R3{-b, +a, 0};
+        point_list[k + 6]  = R3{+b, -a, 0};
+        point_list[k + 7]  = R3{-b, -a, 0};
+        point_list[k + 8]  = R3{+a, 0, +b};
+        point_list[k + 9]  = R3{+a, 0, -b};
+        point_list[k + 10] = R3{-a, 0, +b};
+        point_list[k + 11] = R3{-a, 0, -b};
+        point_list[k + 12] = R3{+b, 0, +a};
+        point_list[k + 13] = R3{-b, 0, +a};
+        point_list[k + 14] = R3{+b, 0, -a};
+        point_list[k + 15] = R3{-b, 0, -a};
+        point_list[k + 16] = R3{0, +a, +b};
+        point_list[k + 17] = R3{0, +a, -b};
+        point_list[k + 18] = R3{0, -a, +b};
+        point_list[k + 19] = R3{0, -a, -b};
+        point_list[k + 20] = R3{0, +b, +a};
+        point_list[k + 21] = R3{0, -b, +a};
+        point_list[k + 22] = R3{0, +b, -a};
+        point_list[k + 23] = R3{0, -b, -a};
 
         for (size_t l = 0; l < 24; ++l) {
           weight_list[k + l] = w;
@@ -134,30 +134,30 @@ CubeGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
         const double a = position_list[0];
         const double b = position_list[1];
 
-        point_list[k + 0]  = {+a, +a, +b};
-        point_list[k + 1]  = {+a, +a, -b};
-        point_list[k + 2]  = {+a, -a, +b};
-        point_list[k + 3]  = {+a, -a, -b};
-        point_list[k + 4]  = {-a, +a, +b};
-        point_list[k + 5]  = {-a, +a, -b};
-        point_list[k + 6]  = {-a, -a, +b};
-        point_list[k + 7]  = {-a, -a, -b};
-        point_list[k + 8]  = {+a, +b, +a};
-        point_list[k + 9]  = {+a, -b, +a};
-        point_list[k + 10] = {+a, +b, -a};
-        point_list[k + 11] = {+a, -b, -a};
-        point_list[k + 12] = {-a, +b, +a};
-        point_list[k + 13] = {-a, -b, +a};
-        point_list[k + 14] = {-a, +b, -a};
-        point_list[k + 15] = {-a, -b, -a};
-        point_list[k + 16] = {+b, +a, +a};
-        point_list[k + 17] = {-b, +a, +a};
-        point_list[k + 18] = {+b, +a, -a};
-        point_list[k + 19] = {-b, +a, -a};
-        point_list[k + 20] = {+b, -a, +a};
-        point_list[k + 21] = {-b, -a, +a};
-        point_list[k + 22] = {+b, -a, -a};
-        point_list[k + 23] = {-b, -a, -a};
+        point_list[k + 0]  = R3{+a, +a, +b};
+        point_list[k + 1]  = R3{+a, +a, -b};
+        point_list[k + 2]  = R3{+a, -a, +b};
+        point_list[k + 3]  = R3{+a, -a, -b};
+        point_list[k + 4]  = R3{-a, +a, +b};
+        point_list[k + 5]  = R3{-a, +a, -b};
+        point_list[k + 6]  = R3{-a, -a, +b};
+        point_list[k + 7]  = R3{-a, -a, -b};
+        point_list[k + 8]  = R3{+a, +b, +a};
+        point_list[k + 9]  = R3{+a, -b, +a};
+        point_list[k + 10] = R3{+a, +b, -a};
+        point_list[k + 11] = R3{+a, -b, -a};
+        point_list[k + 12] = R3{-a, +b, +a};
+        point_list[k + 13] = R3{-a, -b, +a};
+        point_list[k + 14] = R3{-a, +b, -a};
+        point_list[k + 15] = R3{-a, -b, -a};
+        point_list[k + 16] = R3{+b, +a, +a};
+        point_list[k + 17] = R3{-b, +a, +a};
+        point_list[k + 18] = R3{+b, +a, -a};
+        point_list[k + 19] = R3{-b, +a, -a};
+        point_list[k + 20] = R3{+b, -a, +a};
+        point_list[k + 21] = R3{-b, -a, +a};
+        point_list[k + 22] = R3{+b, -a, -a};
+        point_list[k + 23] = R3{-b, -a, -a};
 
         for (size_t l = 0; l < 24; ++l) {
           weight_list[k + l] = w;
@@ -172,54 +172,54 @@ CubeGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
         const double b = position_list[1];
         const double c = position_list[2];
 
-        point_list[k + 0]  = {+a, +b, +c};
-        point_list[k + 1]  = {+a, +b, -c};
-        point_list[k + 2]  = {+a, -b, +c};
-        point_list[k + 3]  = {+a, -b, -c};
-        point_list[k + 4]  = {-a, +b, +c};
-        point_list[k + 5]  = {-a, +b, -c};
-        point_list[k + 6]  = {-a, -b, +c};
-        point_list[k + 7]  = {-a, -b, -c};
-        point_list[k + 8]  = {+a, +c, +b};
-        point_list[k + 9]  = {+a, -c, +b};
-        point_list[k + 10] = {+a, +c, -b};
-        point_list[k + 11] = {+a, -c, -b};
-        point_list[k + 12] = {-a, +c, +b};
-        point_list[k + 13] = {-a, -c, +b};
-        point_list[k + 14] = {-a, +c, -b};
-        point_list[k + 15] = {-a, -c, -b};
-        point_list[k + 16] = {+b, +a, +c};
-        point_list[k + 17] = {+b, +a, -c};
-        point_list[k + 18] = {-b, +a, +c};
-        point_list[k + 19] = {-b, +a, -c};
-        point_list[k + 20] = {+b, -a, +c};
-        point_list[k + 21] = {+b, -a, -c};
-        point_list[k + 22] = {-b, -a, +c};
-        point_list[k + 23] = {-b, -a, -c};
-        point_list[k + 24] = {+b, +c, +a};
-        point_list[k + 25] = {+b, -c, +a};
-        point_list[k + 26] = {-b, +c, +a};
-        point_list[k + 27] = {-b, -c, +a};
-        point_list[k + 28] = {+b, +c, -a};
-        point_list[k + 29] = {+b, -c, -a};
-        point_list[k + 30] = {-b, +c, -a};
-        point_list[k + 31] = {-b, -c, -a};
-        point_list[k + 32] = {+c, +a, +b};
-        point_list[k + 33] = {-c, +a, +b};
-        point_list[k + 34] = {+c, +a, -b};
-        point_list[k + 35] = {-c, +a, -b};
-        point_list[k + 36] = {+c, -a, +b};
-        point_list[k + 37] = {-c, -a, +b};
-        point_list[k + 38] = {+c, -a, -b};
-        point_list[k + 39] = {-c, -a, -b};
-        point_list[k + 40] = {+c, +b, +a};
-        point_list[k + 41] = {-c, +b, +a};
-        point_list[k + 42] = {+c, -b, +a};
-        point_list[k + 43] = {-c, -b, +a};
-        point_list[k + 44] = {+c, +b, -a};
-        point_list[k + 45] = {-c, +b, -a};
-        point_list[k + 46] = {+c, -b, -a};
-        point_list[k + 47] = {-c, -b, -a};
+        point_list[k + 0]  = R3{+a, +b, +c};
+        point_list[k + 1]  = R3{+a, +b, -c};
+        point_list[k + 2]  = R3{+a, -b, +c};
+        point_list[k + 3]  = R3{+a, -b, -c};
+        point_list[k + 4]  = R3{-a, +b, +c};
+        point_list[k + 5]  = R3{-a, +b, -c};
+        point_list[k + 6]  = R3{-a, -b, +c};
+        point_list[k + 7]  = R3{-a, -b, -c};
+        point_list[k + 8]  = R3{+a, +c, +b};
+        point_list[k + 9]  = R3{+a, -c, +b};
+        point_list[k + 10] = R3{+a, +c, -b};
+        point_list[k + 11] = R3{+a, -c, -b};
+        point_list[k + 12] = R3{-a, +c, +b};
+        point_list[k + 13] = R3{-a, -c, +b};
+        point_list[k + 14] = R3{-a, +c, -b};
+        point_list[k + 15] = R3{-a, -c, -b};
+        point_list[k + 16] = R3{+b, +a, +c};
+        point_list[k + 17] = R3{+b, +a, -c};
+        point_list[k + 18] = R3{-b, +a, +c};
+        point_list[k + 19] = R3{-b, +a, -c};
+        point_list[k + 20] = R3{+b, -a, +c};
+        point_list[k + 21] = R3{+b, -a, -c};
+        point_list[k + 22] = R3{-b, -a, +c};
+        point_list[k + 23] = R3{-b, -a, -c};
+        point_list[k + 24] = R3{+b, +c, +a};
+        point_list[k + 25] = R3{+b, -c, +a};
+        point_list[k + 26] = R3{-b, +c, +a};
+        point_list[k + 27] = R3{-b, -c, +a};
+        point_list[k + 28] = R3{+b, +c, -a};
+        point_list[k + 29] = R3{+b, -c, -a};
+        point_list[k + 30] = R3{-b, +c, -a};
+        point_list[k + 31] = R3{-b, -c, -a};
+        point_list[k + 32] = R3{+c, +a, +b};
+        point_list[k + 33] = R3{-c, +a, +b};
+        point_list[k + 34] = R3{+c, +a, -b};
+        point_list[k + 35] = R3{-c, +a, -b};
+        point_list[k + 36] = R3{+c, -a, +b};
+        point_list[k + 37] = R3{-c, -a, +b};
+        point_list[k + 38] = R3{+c, -a, -b};
+        point_list[k + 39] = R3{-c, -a, -b};
+        point_list[k + 40] = R3{+c, +b, +a};
+        point_list[k + 41] = R3{-c, +b, +a};
+        point_list[k + 42] = R3{+c, -b, +a};
+        point_list[k + 43] = R3{-c, -b, +a};
+        point_list[k + 44] = R3{+c, +b, -a};
+        point_list[k + 45] = R3{-c, +b, -a};
+        point_list[k + 46] = R3{+c, -b, -a};
+        point_list[k + 47] = R3{-c, -b, -a};
 
         for (size_t l = 0; l < 48; ++l) {
           weight_list[k + l] = w;
diff --git a/src/analysis/PrismGaussQuadrature.cpp b/src/analysis/PrismGaussQuadrature.cpp
index c2cde582ba1c4e2b6b5077e976a43dea5deaca77..10c8fab3720ab794299522452e5528d25f38e394 100644
--- a/src/analysis/PrismGaussQuadrature.cpp
+++ b/src/analysis/PrismGaussQuadrature.cpp
@@ -19,9 +19,9 @@ PrismGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
 
     Assert(point_list.size() == weight_list.size());
 
-    const R2 A = {0, 0};
-    const R2 B = {1, 0};
-    const R2 C = {0, 1};
+    const R2 A{0, 0};
+    const R2 B{1, 0};
+    const R2 C{0, 1};
 
     size_t k = 0;
     for (size_t i = 0; i < descriptor_list.size(); ++i) {
@@ -31,7 +31,7 @@ PrismGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
       case 1: {
         Assert(value_list.size() == 0);
 
-        point_list[k]  = {1. / 3, 1. / 3, 0.};
+        point_list[k]  = R3{1. / 3, 1. / 3, 0.};
         weight_list[k] = w;
 
         k += 1;
@@ -40,8 +40,8 @@ PrismGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
       case 2: {
         Assert(value_list.size() == 1);
         const double c    = value_list[0];
-        point_list[k + 0] = {1. / 3, 1. / 3, +c};
-        point_list[k + 1] = {1. / 3, 1. / 3, -c};
+        point_list[k + 0] = R3{1. / 3, 1. / 3, +c};
+        point_list[k + 1] = R3{1. / 3, 1. / 3, -c};
 
         for (size_t l = 0; l < 2; ++l) {
           weight_list[k + l] = w;
diff --git a/src/analysis/PyramidGaussQuadrature.cpp b/src/analysis/PyramidGaussQuadrature.cpp
index 869408a914c58c2402e2b9b8c7954549466a97f7..4b4d17359486f0f44bff24dfed013c790f50a36e 100644
--- a/src/analysis/PyramidGaussQuadrature.cpp
+++ b/src/analysis/PyramidGaussQuadrature.cpp
@@ -27,7 +27,7 @@ PyramidGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
         Assert(value_list.size() == 1);
         const double z = value_list[0];
 
-        point_list[k]  = {0, 0, z};
+        point_list[k]  = R3{0, 0, z};
         weight_list[k] = w;
         ++k;
         break;
@@ -37,10 +37,10 @@ PyramidGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
         const double a = value_list[0];
         const double z = value_list[1];
 
-        point_list[k + 0] = {+a, 0, z};
-        point_list[k + 1] = {-a, 0, z};
-        point_list[k + 2] = {0, +a, z};
-        point_list[k + 3] = {0, -a, z};
+        point_list[k + 0] = R3{+a, 0, z};
+        point_list[k + 1] = R3{-a, 0, z};
+        point_list[k + 2] = R3{0, +a, z};
+        point_list[k + 3] = R3{0, -a, z};
 
         for (size_t l = 0; l < 4; ++l) {
           weight_list[k + l] = w;
@@ -54,10 +54,10 @@ PyramidGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
         const double a = value_list[0];
         const double z = value_list[1];
 
-        point_list[k + 0] = {+a, +a, z};
-        point_list[k + 1] = {+a, -a, z};
-        point_list[k + 2] = {-a, +a, z};
-        point_list[k + 3] = {-a, -a, z};
+        point_list[k + 0] = R3{+a, +a, z};
+        point_list[k + 1] = R3{+a, -a, z};
+        point_list[k + 2] = R3{-a, +a, z};
+        point_list[k + 3] = R3{-a, -a, z};
 
         for (size_t l = 0; l < 4; ++l) {
           weight_list[k + l] = w;
@@ -72,14 +72,14 @@ PyramidGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
         const double b = value_list[1];
         const double z = value_list[2];
 
-        point_list[k + 0] = {+a, +b, z};
-        point_list[k + 1] = {+a, -b, z};
-        point_list[k + 2] = {-a, +b, z};
-        point_list[k + 3] = {-a, -b, z};
-        point_list[k + 4] = {+b, +a, z};
-        point_list[k + 5] = {-b, +a, z};
-        point_list[k + 6] = {+b, -a, z};
-        point_list[k + 7] = {-b, -a, z};
+        point_list[k + 0] = R3{+a, +b, z};
+        point_list[k + 1] = R3{+a, -b, z};
+        point_list[k + 2] = R3{-a, +b, z};
+        point_list[k + 3] = R3{-a, -b, z};
+        point_list[k + 4] = R3{+b, +a, z};
+        point_list[k + 5] = R3{-b, +a, z};
+        point_list[k + 6] = R3{+b, -a, z};
+        point_list[k + 7] = R3{-b, -a, z};
 
         for (size_t l = 0; l < 8; ++l) {
           weight_list[k + l] = w;
diff --git a/src/analysis/SquareGaussQuadrature.cpp b/src/analysis/SquareGaussQuadrature.cpp
index f0df6af9f2cbbc23318e861403dbaadaa844d749..20e27c8c02d9ebce1c2305834aeda3b6430a0d27 100644
--- a/src/analysis/SquareGaussQuadrature.cpp
+++ b/src/analysis/SquareGaussQuadrature.cpp
@@ -24,7 +24,7 @@ SquareGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
       case 1: {
         Assert(position_list.size() == 0);
 
-        point_list[k]  = {0, 0};
+        point_list[k]  = R2{0, 0};
         weight_list[k] = w;
 
         k += 1;
@@ -34,10 +34,10 @@ SquareGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
         Assert(position_list.size() == 1);
         const double& a = position_list[0];
 
-        point_list[k + 0] = {+a, 0};
-        point_list[k + 1] = {-a, 0};
-        point_list[k + 2] = {0, +a};
-        point_list[k + 3] = {0, -a};
+        point_list[k + 0] = R2{+a, 0};
+        point_list[k + 1] = R2{-a, 0};
+        point_list[k + 2] = R2{0, +a};
+        point_list[k + 3] = R2{0, -a};
 
         for (size_t l = 0; l < 4; ++l) {
           weight_list[k + l] = w;
@@ -50,10 +50,10 @@ SquareGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
         Assert(position_list.size() == 1);
         const double& a = position_list[0];
 
-        point_list[k + 0] = {+a, +a};
-        point_list[k + 1] = {+a, -a};
-        point_list[k + 2] = {-a, +a};
-        point_list[k + 3] = {-a, -a};
+        point_list[k + 0] = R2{+a, +a};
+        point_list[k + 1] = R2{+a, -a};
+        point_list[k + 2] = R2{-a, +a};
+        point_list[k + 3] = R2{-a, -a};
 
         for (size_t l = 0; l < 4; ++l) {
           weight_list[k + l] = w;
@@ -67,14 +67,14 @@ SquareGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
         const double& a = position_list[0];
         const double& b = position_list[1];
 
-        point_list[k + 0] = {+a, +b};
-        point_list[k + 1] = {+a, -b};
-        point_list[k + 2] = {-a, +b};
-        point_list[k + 3] = {-a, -b};
-        point_list[k + 4] = {+b, +a};
-        point_list[k + 5] = {+b, -a};
-        point_list[k + 6] = {-b, +a};
-        point_list[k + 7] = {-b, -a};
+        point_list[k + 0] = R2{+a, +b};
+        point_list[k + 1] = R2{+a, -b};
+        point_list[k + 2] = R2{-a, +b};
+        point_list[k + 3] = R2{-a, -b};
+        point_list[k + 4] = R2{+b, +a};
+        point_list[k + 5] = R2{+b, -a};
+        point_list[k + 6] = R2{-b, +a};
+        point_list[k + 7] = R2{-b, -a};
 
         for (size_t l = 0; l < 8; ++l) {
           weight_list[k + l] = w;
diff --git a/src/analysis/TensorialGaussLegendreQuadrature.cpp b/src/analysis/TensorialGaussLegendreQuadrature.cpp
index 9a4b524c0f70cbc3a68088bd5b14ad17d6c521ea..634ac9f9559ee87a6926dccf28d6c966aede0acc 100644
--- a/src/analysis/TensorialGaussLegendreQuadrature.cpp
+++ b/src/analysis/TensorialGaussLegendreQuadrature.cpp
@@ -5,23 +5,25 @@ template <>
 void
 TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degree)
 {
+  using R1 = TinyVector<1>;
+
   const size_t nb_points = degree / 2 + 1;
 
-  SmallArray<TinyVector<1>> point_list(nb_points);
+  SmallArray<R1> point_list(nb_points);
   SmallArray<double> weight_list(nb_points);
 
   switch (degree) {
   case 0:
   case 1: {
-    point_list[0] = 0;
+    point_list[0] = R1{0};
 
     weight_list[0] = 2;
     break;
   }
   case 2:
   case 3: {
-    point_list[0] = -0.5773502691896257645091488;
-    point_list[1] = +0.5773502691896257645091488;
+    point_list[0] = R1{-0.5773502691896257645091488};
+    point_list[1] = R1{+0.5773502691896257645091488};
 
     weight_list[0] = 1;
     weight_list[1] = 1;
@@ -29,9 +31,9 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degr
   }
   case 4:
   case 5: {
-    point_list[0] = -0.7745966692414833770358531;
-    point_list[1] = 0;
-    point_list[2] = +0.7745966692414833770358531;
+    point_list[0] = R1{-0.7745966692414833770358531};
+    point_list[1] = R1{0};
+    point_list[2] = R1{+0.7745966692414833770358531};
 
     weight_list[0] = 0.5555555555555555555555556;
     weight_list[1] = 0.8888888888888888888888889;
@@ -40,10 +42,10 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degr
   }
   case 6:
   case 7: {
-    point_list[0] = -0.8611363115940525752239465;
-    point_list[1] = -0.3399810435848562648026658;
-    point_list[2] = +0.3399810435848562648026658;
-    point_list[3] = +0.8611363115940525752239465;
+    point_list[0] = R1{-0.8611363115940525752239465};
+    point_list[1] = R1{-0.3399810435848562648026658};
+    point_list[2] = R1{+0.3399810435848562648026658};
+    point_list[3] = R1{+0.8611363115940525752239465};
 
     weight_list[0] = 0.3478548451374538573730639;
     weight_list[1] = 0.6521451548625461426269361;
@@ -53,11 +55,11 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degr
   }
   case 8:
   case 9: {
-    point_list[0] = -0.9061798459386639927976269;
-    point_list[1] = -0.5384693101056830910363144;
-    point_list[2] = 0;
-    point_list[3] = +0.5384693101056830910363144;
-    point_list[4] = +0.9061798459386639927976269;
+    point_list[0] = R1{-0.9061798459386639927976269};
+    point_list[1] = R1{-0.5384693101056830910363144};
+    point_list[2] = R1{0};
+    point_list[3] = R1{+0.5384693101056830910363144};
+    point_list[4] = R1{+0.9061798459386639927976269};
 
     weight_list[0] = 0.2369268850561890875142640;
     weight_list[1] = 0.4786286704993664680412915;
@@ -68,12 +70,12 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degr
   }
   case 10:
   case 11: {
-    point_list[0] = -0.9324695142031520278123016;
-    point_list[1] = -0.6612093864662645136613996;
-    point_list[2] = -0.2386191860831969086305017;
-    point_list[3] = +0.2386191860831969086305017;
-    point_list[4] = +0.6612093864662645136613996;
-    point_list[5] = +0.9324695142031520278123016;
+    point_list[0] = R1{-0.9324695142031520278123016};
+    point_list[1] = R1{-0.6612093864662645136613996};
+    point_list[2] = R1{-0.2386191860831969086305017};
+    point_list[3] = R1{+0.2386191860831969086305017};
+    point_list[4] = R1{+0.6612093864662645136613996};
+    point_list[5] = R1{+0.9324695142031520278123016};
 
     weight_list[0] = 0.1713244923791703450402961;
     weight_list[1] = 0.3607615730481386075698335;
@@ -85,13 +87,13 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degr
   }
   case 12:
   case 13: {
-    point_list[0] = -0.9491079123427585245261897;
-    point_list[1] = -0.7415311855993944398638648;
-    point_list[2] = -0.4058451513773971669066064;
-    point_list[3] = 0;
-    point_list[4] = +0.4058451513773971669066064;
-    point_list[5] = +0.7415311855993944398638648;
-    point_list[6] = +0.9491079123427585245261897;
+    point_list[0] = R1{-0.9491079123427585245261897};
+    point_list[1] = R1{-0.7415311855993944398638648};
+    point_list[2] = R1{-0.4058451513773971669066064};
+    point_list[3] = R1{0};
+    point_list[4] = R1{+0.4058451513773971669066064};
+    point_list[5] = R1{+0.7415311855993944398638648};
+    point_list[6] = R1{+0.9491079123427585245261897};
 
     weight_list[0] = 0.1294849661688696932706114;
     weight_list[1] = 0.2797053914892766679014678;
@@ -104,14 +106,14 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degr
   }
   case 14:
   case 15: {
-    point_list[0] = -0.9602898564975362316835609;
-    point_list[1] = -0.7966664774136267395915539;
-    point_list[2] = -0.5255324099163289858177390;
-    point_list[3] = -0.1834346424956498049394761;
-    point_list[4] = +0.1834346424956498049394761;
-    point_list[5] = +0.5255324099163289858177390;
-    point_list[6] = +0.7966664774136267395915539;
-    point_list[7] = +0.9602898564975362316835609;
+    point_list[0] = R1{-0.9602898564975362316835609};
+    point_list[1] = R1{-0.7966664774136267395915539};
+    point_list[2] = R1{-0.5255324099163289858177390};
+    point_list[3] = R1{-0.1834346424956498049394761};
+    point_list[4] = R1{+0.1834346424956498049394761};
+    point_list[5] = R1{+0.5255324099163289858177390};
+    point_list[6] = R1{+0.7966664774136267395915539};
+    point_list[7] = R1{+0.9602898564975362316835609};
 
     weight_list[0] = 0.1012285362903762591525314;
     weight_list[1] = 0.2223810344533744705443560;
@@ -125,15 +127,15 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degr
   }
   case 16:
   case 17: {
-    point_list[0] = -0.9681602395076260898355762;
-    point_list[1] = -0.8360311073266357942994298;
-    point_list[2] = -0.6133714327005903973087020;
-    point_list[3] = -0.3242534234038089290385380;
-    point_list[4] = 0;
-    point_list[5] = +0.3242534234038089290385380;
-    point_list[6] = +0.6133714327005903973087020;
-    point_list[7] = +0.8360311073266357942994298;
-    point_list[8] = +0.9681602395076260898355762;
+    point_list[0] = R1{-0.9681602395076260898355762};
+    point_list[1] = R1{-0.8360311073266357942994298};
+    point_list[2] = R1{-0.6133714327005903973087020};
+    point_list[3] = R1{-0.3242534234038089290385380};
+    point_list[4] = R1{0};
+    point_list[5] = R1{+0.3242534234038089290385380};
+    point_list[6] = R1{+0.6133714327005903973087020};
+    point_list[7] = R1{+0.8360311073266357942994298};
+    point_list[8] = R1{+0.9681602395076260898355762};
 
     weight_list[0] = 0.0812743883615744119718922;
     weight_list[1] = 0.1806481606948574040584720;
@@ -148,16 +150,16 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degr
   }
   case 18:
   case 19: {
-    point_list[0] = -0.9739065285171717200779640;
-    point_list[1] = -0.8650633666889845107320967;
-    point_list[2] = -0.6794095682990244062343274;
-    point_list[3] = -0.4333953941292471907992659;
-    point_list[4] = -0.1488743389816312108848260;
-    point_list[5] = +0.1488743389816312108848260;
-    point_list[6] = +0.4333953941292471907992659;
-    point_list[7] = +0.6794095682990244062343274;
-    point_list[8] = +0.8650633666889845107320967;
-    point_list[9] = +0.9739065285171717200779640;
+    point_list[0] = R1{-0.9739065285171717200779640};
+    point_list[1] = R1{-0.8650633666889845107320967};
+    point_list[2] = R1{-0.6794095682990244062343274};
+    point_list[3] = R1{-0.4333953941292471907992659};
+    point_list[4] = R1{-0.1488743389816312108848260};
+    point_list[5] = R1{+0.1488743389816312108848260};
+    point_list[6] = R1{+0.4333953941292471907992659};
+    point_list[7] = R1{+0.6794095682990244062343274};
+    point_list[8] = R1{+0.8650633666889845107320967};
+    point_list[9] = R1{+0.9739065285171717200779640};
 
     weight_list[0] = 0.0666713443086881375935688;
     weight_list[1] = 0.1494513491505805931457763;
@@ -173,17 +175,17 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degr
   }
   case 20:
   case 21: {
-    point_list[0]  = -0.9782286581460569928039380;
-    point_list[1]  = -0.8870625997680952990751578;
-    point_list[2]  = -0.7301520055740493240934163;
-    point_list[3]  = -0.5190961292068118159257257;
-    point_list[4]  = -0.2695431559523449723315320;
-    point_list[5]  = 0;
-    point_list[6]  = +0.2695431559523449723315320;
-    point_list[7]  = +0.5190961292068118159257257;
-    point_list[8]  = +0.7301520055740493240934163;
-    point_list[9]  = +0.8870625997680952990751578;
-    point_list[10] = +0.9782286581460569928039380;
+    point_list[0]  = R1{-0.9782286581460569928039380};
+    point_list[1]  = R1{-0.8870625997680952990751578};
+    point_list[2]  = R1{-0.7301520055740493240934163};
+    point_list[3]  = R1{-0.5190961292068118159257257};
+    point_list[4]  = R1{-0.2695431559523449723315320};
+    point_list[5]  = R1{0};
+    point_list[6]  = R1{+0.2695431559523449723315320};
+    point_list[7]  = R1{+0.5190961292068118159257257};
+    point_list[8]  = R1{+0.7301520055740493240934163};
+    point_list[9]  = R1{+0.8870625997680952990751578};
+    point_list[10] = R1{+0.9782286581460569928039380};
 
     weight_list[0]  = 0.0556685671161736664827537;
     weight_list[1]  = 0.1255803694649046246346940;
@@ -200,18 +202,18 @@ TensorialGaussLegendreQuadrature<1>::_buildPointAndWeightLists(const size_t degr
   }
   case 22:
   case 23: {
-    point_list[0]  = -0.9815606342467192506905491;
-    point_list[1]  = -0.9041172563704748566784659;
-    point_list[2]  = -0.7699026741943046870368938;
-    point_list[3]  = -0.5873179542866174472967024;
-    point_list[4]  = -0.3678314989981801937526915;
-    point_list[5]  = -0.1252334085114689154724414;
-    point_list[6]  = +0.1252334085114689154724414;
-    point_list[7]  = +0.3678314989981801937526915;
-    point_list[8]  = +0.5873179542866174472967024;
-    point_list[9]  = +0.7699026741943046870368938;
-    point_list[10] = +0.9041172563704748566784659;
-    point_list[11] = +0.9815606342467192506905491;
+    point_list[0]  = R1{-0.9815606342467192506905491};
+    point_list[1]  = R1{-0.9041172563704748566784659};
+    point_list[2]  = R1{-0.7699026741943046870368938};
+    point_list[3]  = R1{-0.5873179542866174472967024};
+    point_list[4]  = R1{-0.3678314989981801937526915};
+    point_list[5]  = R1{-0.1252334085114689154724414};
+    point_list[6]  = R1{+0.1252334085114689154724414};
+    point_list[7]  = R1{+0.3678314989981801937526915};
+    point_list[8]  = R1{+0.5873179542866174472967024};
+    point_list[9]  = R1{+0.7699026741943046870368938};
+    point_list[10] = R1{+0.9041172563704748566784659};
+    point_list[11] = R1{+0.9815606342467192506905491};
 
     weight_list[0]  = 0.0471753363865118271946160;
     weight_list[1]  = 0.1069393259953184309602547;
diff --git a/src/analysis/TensorialGaussLobattoQuadrature.cpp b/src/analysis/TensorialGaussLobattoQuadrature.cpp
index 3ab97460b55dad3a2a48a24db56f5c575542361e..d848943751a6c5e3487b5fcc24394c66985a271a 100644
--- a/src/analysis/TensorialGaussLobattoQuadrature.cpp
+++ b/src/analysis/TensorialGaussLobattoQuadrature.cpp
@@ -5,6 +5,8 @@ template <>
 void
 TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists(const size_t degree)
 {
+  using R1 = TinyVector<1>;
+
   const size_t nb_points = degree / 2 + 2;
 
   SmallArray<TinyVector<1>> point_list(nb_points);
@@ -13,8 +15,8 @@ TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists(const size_t degre
   switch (degree) {
   case 0:
   case 1: {
-    point_list[0] = -1;
-    point_list[1] = +1;
+    point_list[0] = R1{-1};
+    point_list[1] = R1{+1};
 
     weight_list[0] = 1;
     weight_list[1] = 1;
@@ -22,9 +24,9 @@ TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists(const size_t degre
   }
   case 2:
   case 3: {
-    point_list[0] = -1;
-    point_list[1] = 0;
-    point_list[2] = +1;
+    point_list[0] = R1{-1};
+    point_list[1] = R1{0};
+    point_list[2] = R1{+1};
 
     weight_list[0] = 0.3333333333333333333333333;
     weight_list[1] = 1.3333333333333333333333333;
@@ -33,10 +35,10 @@ TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists(const size_t degre
   }
   case 4:
   case 5: {
-    point_list[0] = -1;
-    point_list[1] = -0.4472135954999579392818347;
-    point_list[2] = +0.4472135954999579392818347;
-    point_list[3] = +1;
+    point_list[0] = R1{-1};
+    point_list[1] = R1{-0.4472135954999579392818347};
+    point_list[2] = R1{+0.4472135954999579392818347};
+    point_list[3] = R1{+1};
 
     weight_list[0] = 0.1666666666666666666666667;
     weight_list[1] = 0.8333333333333333333333333;
@@ -46,11 +48,11 @@ TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists(const size_t degre
   }
   case 6:
   case 7: {
-    point_list[0] = -1;
-    point_list[1] = -0.6546536707079771437982925;
-    point_list[2] = 0;
-    point_list[3] = +0.6546536707079771437982925;
-    point_list[4] = +1;
+    point_list[0] = R1{-1};
+    point_list[1] = R1{-0.6546536707079771437982925};
+    point_list[2] = R1{0};
+    point_list[3] = R1{+0.6546536707079771437982925};
+    point_list[4] = R1{+1};
 
     weight_list[0] = 0.1;
     weight_list[1] = 0.5444444444444444444444444;
@@ -61,12 +63,12 @@ TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists(const size_t degre
   }
   case 8:
   case 9: {
-    point_list[0] = -1;
-    point_list[1] = -0.7650553239294646928510030;
-    point_list[2] = -0.2852315164806450963141510;
-    point_list[3] = +0.2852315164806450963141510;
-    point_list[4] = +0.7650553239294646928510030;
-    point_list[5] = +1;
+    point_list[0] = R1{-1};
+    point_list[1] = R1{-0.7650553239294646928510030};
+    point_list[2] = R1{-0.2852315164806450963141510};
+    point_list[3] = R1{+0.2852315164806450963141510};
+    point_list[4] = R1{+0.7650553239294646928510030};
+    point_list[5] = R1{+1};
 
     weight_list[0] = 0.0666666666666666666666667;
     weight_list[1] = 0.3784749562978469803166128;
@@ -78,13 +80,13 @@ TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists(const size_t degre
   }
   case 10:
   case 11: {
-    point_list[0] = -1;
-    point_list[1] = -0.8302238962785669298720322;
-    point_list[2] = -0.4688487934707142138037719;
-    point_list[3] = 0;
-    point_list[4] = +0.4688487934707142138037719;
-    point_list[5] = +0.8302238962785669298720322;
-    point_list[6] = +1;
+    point_list[0] = R1{-1};
+    point_list[1] = R1{-0.8302238962785669298720322};
+    point_list[2] = R1{-0.4688487934707142138037719};
+    point_list[3] = R1{0};
+    point_list[4] = R1{+0.4688487934707142138037719};
+    point_list[5] = R1{+0.8302238962785669298720322};
+    point_list[6] = R1{+1};
 
     weight_list[0] = 0.0476190476190476190476190;
     weight_list[1] = 0.2768260473615659480107004;
@@ -97,14 +99,14 @@ TensorialGaussLobattoQuadrature<1>::_buildPointAndWeightLists(const size_t degre
   }
   case 12:
   case 13: {
-    point_list[0] = -1;
-    point_list[1] = -0.8717401485096066153374457;
-    point_list[2] = -0.5917001814331423021445107;
-    point_list[3] = -0.2092992179024788687686573;
-    point_list[4] = +0.2092992179024788687686573;
-    point_list[5] = +0.5917001814331423021445107;
-    point_list[6] = +0.8717401485096066153374457;
-    point_list[7] = +1;
+    point_list[0] = R1{-1};
+    point_list[1] = R1{-0.8717401485096066153374457};
+    point_list[2] = R1{-0.5917001814331423021445107};
+    point_list[3] = R1{-0.2092992179024788687686573};
+    point_list[4] = R1{+0.2092992179024788687686573};
+    point_list[5] = R1{+0.5917001814331423021445107};
+    point_list[6] = R1{+0.8717401485096066153374457};
+    point_list[7] = R1{+1};
 
     weight_list[0] = 0.0357142857142857142857143;
     weight_list[1] = 0.2107042271435060393829921;
diff --git a/src/analysis/TetrahedronGaussQuadrature.cpp b/src/analysis/TetrahedronGaussQuadrature.cpp
index 644e2e29387821a274fa2ae8dd2861c60f2c4abc..e5ee82d3fb149ab7c50647d98092754f146dcab1 100644
--- a/src/analysis/TetrahedronGaussQuadrature.cpp
+++ b/src/analysis/TetrahedronGaussQuadrature.cpp
@@ -16,10 +16,10 @@ TetrahedronGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
   auto fill_quadrature_points = [](auto descriptor_list, auto& point_list, auto& weight_list) {
     Assert(point_list.size() == weight_list.size());
 
-    const R3 A = {0, 0, 0};
-    const R3 B = {1, 0, 0};
-    const R3 C = {0, 1, 0};
-    const R3 D = {0, 0, 1};
+    const R3 A{0, 0, 0};
+    const R3 B{1, 0, 0};
+    const R3 C{0, 1, 0};
+    const R3 D{0, 0, 1};
 
     size_t k = 0;
     for (size_t i = 0; i < descriptor_list.size(); ++i) {
diff --git a/src/analysis/TriangleGaussQuadrature.cpp b/src/analysis/TriangleGaussQuadrature.cpp
index df6a9a1ced3aa15d098aaa7696eb67af41c6090b..19013da54f79268329416b5499fd80748612b4cc 100644
--- a/src/analysis/TriangleGaussQuadrature.cpp
+++ b/src/analysis/TriangleGaussQuadrature.cpp
@@ -16,9 +16,9 @@ TriangleGaussQuadrature::_buildPointAndWeightLists(const size_t degree)
   auto fill_quadrature_points = [](auto descriptor_list, auto& point_list, auto& weight_list) {
     Assert(point_list.size() == weight_list.size());
 
-    const R2 A = {0, 0};
-    const R2 B = {1, 0};
-    const R2 C = {0, 1};
+    const R2 A{0, 0};
+    const R2 B{1, 0};
+    const R2 C{0, 1};
 
     size_t k = 0;
     for (size_t i = 0; i < descriptor_list.size(); ++i) {
diff --git a/src/geometry/CubeTransformation.hpp b/src/geometry/CubeTransformation.hpp
index 7bc12c324e4a55750b41b1584633cc789f8a7b2d..48582eef10148522adb2e4cb60f5c6b0f29df9dc 100644
--- a/src/geometry/CubeTransformation.hpp
+++ b/src/geometry/CubeTransformation.hpp
@@ -21,8 +21,7 @@ class CubeTransformation
   TinyVector<Dimension>
   operator()(const TinyVector<Dimension>& x) const
   {
-    const TinyVector<NumberOfPoints - 1> X =
-      {x[0], x[1], x[2], x[0] * x[1], x[1] * x[2], x[0] * x[2], x[0] * x[1] * x[2]};
+    const TinyVector<NumberOfPoints - 1> X{x[0], x[1], x[2], x[0] * x[1], x[1] * x[2], x[0] * x[2], x[0] * x[1] * x[2]};
     return m_coefficients * X + m_shift;
   }
 
@@ -35,17 +34,17 @@ class CubeTransformation
     const double& y = X[1];
     const double& z = X[2];
 
-    const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 3) * y + T(0, 5) * z + T(0, 6) * y * z,
-                                                T(0, 1) + T(0, 3) * x + T(0, 4) * z + T(0, 6) * x * z,
-                                                T(0, 2) + T(0, 4) * y + T(0, 5) * x + T(0, 6) * x * y,
-                                                //
-                                                T(1, 0) + T(1, 3) * y + T(1, 5) * z + T(1, 6) * y * z,
-                                                T(1, 1) + T(1, 3) * x + T(1, 4) * z + T(1, 6) * x * z,
-                                                T(1, 2) + T(1, 4) * y + T(1, 5) * x + T(1, 6) * x * y,
-                                                //
-                                                T(2, 0) + T(2, 3) * y + T(2, 5) * z + T(2, 6) * y * z,
-                                                T(2, 1) + T(2, 3) * x + T(2, 4) * z + T(2, 6) * x * z,
-                                                T(2, 2) + T(2, 4) * y + T(2, 5) * x + T(2, 6) * x * y};
+    const TinyMatrix<Dimension, Dimension> J{T(0, 0) + T(0, 3) * y + T(0, 5) * z + T(0, 6) * y * z,
+                                             T(0, 1) + T(0, 3) * x + T(0, 4) * z + T(0, 6) * x * z,
+                                             T(0, 2) + T(0, 4) * y + T(0, 5) * x + T(0, 6) * x * y,
+                                             //
+                                             T(1, 0) + T(1, 3) * y + T(1, 5) * z + T(1, 6) * y * z,
+                                             T(1, 1) + T(1, 3) * x + T(1, 4) * z + T(1, 6) * x * z,
+                                             T(1, 2) + T(1, 4) * y + T(1, 5) * x + T(1, 6) * x * y,
+                                             //
+                                             T(2, 0) + T(2, 3) * y + T(2, 5) * z + T(2, 6) * y * z,
+                                             T(2, 1) + T(2, 3) * x + T(2, 4) * z + T(2, 6) * x * z,
+                                             T(2, 2) + T(2, 4) * y + T(2, 5) * x + T(2, 6) * x * y};
 
     return det(J);
   }
diff --git a/src/geometry/PrismTransformation.hpp b/src/geometry/PrismTransformation.hpp
index 5e6e509ed0c92023b81ea13a3e686392dec52fcd..fc83e1911e196e665a8d0954c7df3c81bb2abdf4 100644
--- a/src/geometry/PrismTransformation.hpp
+++ b/src/geometry/PrismTransformation.hpp
@@ -18,7 +18,7 @@ class PrismTransformation
   TinyVector<Dimension>
   operator()(const TinyVector<Dimension>& x) const
   {
-    const TinyVector<NumberOfPoints - 1> X = {x[0], x[1], x[2], x[0] * x[2], x[1] * x[2]};
+    const TinyVector<NumberOfPoints - 1> X{x[0], x[1], x[2], x[0] * x[2], x[1] * x[2]};
     return m_coefficients * X + m_shift;
   }
 
@@ -30,17 +30,17 @@ class PrismTransformation
     const double& y = X[1];
     const double& z = X[2];
 
-    const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 3) * z,   //
-                                                T(0, 1) + T(0, 4) * z,   //
-                                                T(0, 2) + T(0, 3) * x + T(0, 4) * y,
-                                                //
-                                                T(1, 0) + T(1, 3) * z,   //
-                                                T(1, 1) + T(1, 4) * z,   //
-                                                T(1, 2) + T(1, 3) * x + T(1, 4) * y,
-                                                //
-                                                T(2, 0) + T(2, 3) * z,   //
-                                                T(2, 1) + T(2, 4) * z,   //
-                                                T(2, 2) + T(2, 3) * x + T(2, 4) * y};
+    const TinyMatrix<Dimension, Dimension> J{T(0, 0) + T(0, 3) * z,   //
+                                             T(0, 1) + T(0, 4) * z,   //
+                                             T(0, 2) + T(0, 3) * x + T(0, 4) * y,
+                                             //
+                                             T(1, 0) + T(1, 3) * z,   //
+                                             T(1, 1) + T(1, 4) * z,   //
+                                             T(1, 2) + T(1, 3) * x + T(1, 4) * y,
+                                             //
+                                             T(2, 0) + T(2, 3) * z,   //
+                                             T(2, 1) + T(2, 4) * z,   //
+                                             T(2, 2) + T(2, 3) * x + T(2, 4) * y};
     return det(J);
   }
 
diff --git a/src/geometry/PyramidTransformation.hpp b/src/geometry/PyramidTransformation.hpp
index 2fa3eaa8d8b010dda5ac3d591b1183db707ea3f8..5b2a79a1b46a95b18b773b2aacb6448081453441 100644
--- a/src/geometry/PyramidTransformation.hpp
+++ b/src/geometry/PyramidTransformation.hpp
@@ -18,7 +18,7 @@ class PyramidTransformation
   TinyVector<Dimension>
   operator()(const TinyVector<Dimension>& x) const
   {
-    const TinyVector<NumberOfPoints - 1> X = {x[0], x[1], x[2], x[0] * x[1]};
+    const TinyVector<NumberOfPoints - 1> X{x[0], x[1], x[2], x[0] * x[1]};
     return m_coefficients * X + m_shift;
   }
 
@@ -29,17 +29,17 @@ class PyramidTransformation
     const double& x = X[0];
     const double& y = X[1];
 
-    const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 3) * y,   //
-                                                T(0, 1) + T(0, 3) * x,   //
-                                                T(0, 2),
-                                                //
-                                                T(1, 0) + T(1, 3) * y,   //
-                                                T(1, 1) + T(1, 3) * x,   //
-                                                T(1, 2),
-                                                //
-                                                T(2, 0) + T(2, 3) * y,   //
-                                                T(2, 1) + T(2, 3) * x,   //
-                                                T(2, 2)};
+    const TinyMatrix<Dimension, Dimension> J{T(0, 0) + T(0, 3) * y,   //
+                                             T(0, 1) + T(0, 3) * x,   //
+                                             T(0, 2),
+                                             //
+                                             T(1, 0) + T(1, 3) * y,   //
+                                             T(1, 1) + T(1, 3) * x,   //
+                                             T(1, 2),
+                                             //
+                                             T(2, 0) + T(2, 3) * y,   //
+                                             T(2, 1) + T(2, 3) * x,   //
+                                             T(2, 2)};
     return det(J);
   }
 
diff --git a/src/geometry/SquareTransformation.hpp b/src/geometry/SquareTransformation.hpp
index d79ae12c96c696f6be0bcf6245ede0b16c5f2a83..d6474f33e4ecbaa4bcfba525dadc9fd38e95962d 100644
--- a/src/geometry/SquareTransformation.hpp
+++ b/src/geometry/SquareTransformation.hpp
@@ -25,7 +25,7 @@ class SquareTransformation<2>
   TinyVector<Dimension>
   operator()(const TinyVector<2>& x) const
   {
-    const TinyVector<NumberOfPoints - 1> X = {x[0], x[1], x[0] * x[1]};
+    const TinyVector<NumberOfPoints - 1> X{x[0], x[1], x[0] * x[1]};
     return m_coefficients * X + m_shift;
   }
 
@@ -36,11 +36,11 @@ class SquareTransformation<2>
     const double& x = X[0];
     const double& y = X[1];
 
-    const TinyMatrix<Dimension, Dimension> J = {T(0, 0) + T(0, 2) * y,   //
-                                                T(0, 1) + T(0, 2) * x,
-                                                //
-                                                T(1, 0) + T(1, 2) * y,   //
-                                                T(1, 1) + T(1, 2) * x};
+    const TinyMatrix<Dimension, Dimension> J{T(0, 0) + T(0, 2) * y,   //
+                                             T(0, 1) + T(0, 2) * x,
+                                             //
+                                             T(1, 0) + T(1, 2) * y,   //
+                                             T(1, 1) + T(1, 2) * x};
     return det(J);
   }
 
@@ -80,7 +80,7 @@ class SquareTransformation
   TinyVector<Dimension>
   operator()(const TinyVector<2>& x) const
   {
-    const TinyVector<NumberOfPoints - 1> X = {x[0], x[1], x[0] * x[1]};
+    const TinyVector<NumberOfPoints - 1> X{x[0], x[1], x[0] * x[1]};
     return m_coefficients * X + m_shift;
   }
 
@@ -91,13 +91,13 @@ class SquareTransformation
     const double& x = X[0];
     const double& y = X[1];
 
-    const TinyVector<Dimension> dxT = {T(0, 0) + T(0, 2) * y,   //
-                                       T(1, 0) + T(1, 2) * y,   //
-                                       T(2, 0) + T(2, 2) * y};
+    const TinyVector<Dimension> dxT{T(0, 0) + T(0, 2) * y,   //
+                                    T(1, 0) + T(1, 2) * y,   //
+                                    T(2, 0) + T(2, 2) * y};
 
-    const TinyVector<Dimension> dyT = {T(0, 1) + T(0, 2) * x,   //
-                                       T(1, 1) + T(1, 2) * x,   //
-                                       T(2, 1) + T(2, 2) * x};
+    const TinyVector<Dimension> dyT{T(0, 1) + T(0, 2) * x,   //
+                                    T(1, 1) + T(1, 2) * x,   //
+                                    T(2, 1) + T(2, 2) * x};
 
     return l2Norm(crossProduct(dxT, dyT));
   }
diff --git a/src/language/node_processor/AffectationProcessor.hpp b/src/language/node_processor/AffectationProcessor.hpp
index f24e8473be7717f4ddd2e8f20d3d2ab666f19001..002c0e4ed7aa427a20734689089de590cfd96bc0 100644
--- a/src/language/node_processor/AffectationProcessor.hpp
+++ b/src/language/node_processor/AffectationProcessor.hpp
@@ -165,7 +165,7 @@ class AffectationExecutor final : public IAffectationExecutor
                 [&](auto&& v) {
                   using Vi_T = std::decay_t<decltype(v)>;
                   if constexpr (std::is_convertible_v<Vi_T, double>) {
-                    m_lhs = v;
+                    m_lhs = TinyVector<1>(v);
                   } else {
                     // LCOV_EXCL_START
                     throw UnexpectedError("unexpected rhs type in affectation");
@@ -178,7 +178,7 @@ class AffectationExecutor final : public IAffectationExecutor
                 [&](auto&& v) {
                   using Vi_T = std::decay_t<decltype(v)>;
                   if constexpr (std::is_convertible_v<Vi_T, double>) {
-                    m_lhs = v;
+                    m_lhs = TinyMatrix<1>(v);
                   } else {
                     // LCOV_EXCL_START
                     throw UnexpectedError("unexpected rhs type in affectation");
diff --git a/src/language/node_processor/FunctionProcessor.hpp b/src/language/node_processor/FunctionProcessor.hpp
index cead55c20fd2e6b34edb4df510840d683fc8b5b8..c6c81d2938e9f774f0315f6a9af5c2e03c6a19e5 100644
--- a/src/language/node_processor/FunctionProcessor.hpp
+++ b/src/language/node_processor/FunctionProcessor.hpp
@@ -58,6 +58,10 @@ class FunctionExpressionProcessor final : public INodeProcessor
       return ReturnType{ZeroType::zero};
     } else if constexpr (std::is_convertible_v<ExpressionValueType, ReturnType>) {
       return static_cast<ReturnType>(std::get<ExpressionValueType>(m_function_expression.execute(exec_policy)));
+    } else if constexpr (std::is_arithmetic_v<ExpressionValueType> and
+                         (is_tiny_vector_v<ReturnType> or is_tiny_matrix_v<ReturnType>)) {
+      static_assert(ReturnType::Dimension == 1, "invalid conversion");
+      return ReturnType(std::get<ExpressionValueType>(m_function_expression.execute(exec_policy)));
     } else {
       throw UnexpectedError("invalid conversion");
     }
diff --git a/src/mesh/CartesianMeshBuilder.cpp b/src/mesh/CartesianMeshBuilder.cpp
index 6b7ef2b26f8aa8fbe1d78ab24f606eb8c335340d..c0cf86a73886b59e56f417fbf1731c78e00d15d5 100644
--- a/src/mesh/CartesianMeshBuilder.cpp
+++ b/src/mesh/CartesianMeshBuilder.cpp
@@ -23,7 +23,7 @@ CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<1>& a,
                                           const TinyVector<1, uint64_t>& cell_size,
                                           const IConnectivity& connectivity) const
 {
-  const double h = (b[0] - a[0]) / cell_size[0];
+  const TinyVector<1> h{(b[0] - a[0]) / cell_size[0]};
   NodeValue<TinyVector<1>> xr(connectivity);
   parallel_for(
     connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId r) { xr[r] = a + r * h; });
diff --git a/src/mesh/LogicalConnectivityBuilder.cpp b/src/mesh/LogicalConnectivityBuilder.cpp
index bfeddbeca5471469ce5dee98382464a7cedc9746..db14ff1df24cede3afcc7d77e96950b185274546 100644
--- a/src/mesh/LogicalConnectivityBuilder.cpp
+++ b/src/mesh/LogicalConnectivityBuilder.cpp
@@ -42,31 +42,29 @@ LogicalConnectivityBuilder::_buildBoundaryNodeList(
 {
   const TinyVector<2, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1};
 
-  const auto node_number = [&](const TinyVector<2, uint64_t> node_logic_id) {
-    return node_logic_id[0] * node_size[1] + node_logic_id[1];
-  };
+  const auto node_number = [&](const uint64_t i, const uint64_t j) { return i * node_size[1] + j; };
 
   {   // xminymin
     Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({0, 0});
+    boundary_nodes[0] = node_number(0, 0);
     descriptor.addRefItemList(RefNodeList{RefId{10, "XMINYMIN"}, boundary_nodes});
   }
 
   {   // xmaxymin
     Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({cell_size[0], 0});
+    boundary_nodes[0] = node_number(cell_size[0], 0);
     descriptor.addRefItemList(RefNodeList{RefId{11, "XMAXYMIN"}, boundary_nodes});
   }
 
   {   // xmaxymax
     Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({cell_size[0], cell_size[1]});
+    boundary_nodes[0] = node_number(cell_size[0], cell_size[1]);
     descriptor.addRefItemList(RefNodeList{RefId{12, "XMAXYMAX"}, boundary_nodes});
   }
 
   {   // xminymax
     Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({0, cell_size[1]});
+    boundary_nodes[0] = node_number(0, cell_size[1]);
     descriptor.addRefItemList(RefNodeList{RefId{13, "XMINYMAX"}, boundary_nodes});
   }
 }
@@ -78,55 +76,55 @@ LogicalConnectivityBuilder::_buildBoundaryNodeList(const TinyVector<3, uint64_t>
 {
   const TinyVector<3, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1, cell_size[2] + 1};
 
-  const auto node_number = [&](const TinyVector<3, uint64_t>& node_logic_id) {
-    return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
+  const auto node_number = [&](const uint64_t i, const uint64_t j, const uint64_t k) {
+    return (i * node_size[1] + j) * node_size[2] + k;
   };
 
   {   // xminyminzmin
     Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({0, 0, 0});
+    boundary_nodes[0] = node_number(0, 0, 0);
     descriptor.addRefItemList(RefNodeList{RefId{10, "XMINYMINZMIN"}, boundary_nodes});
   }
 
   {   // xmaxyminzmin
     Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({cell_size[0], 0, 0});
+    boundary_nodes[0] = node_number(cell_size[0], 0, 0);
     descriptor.addRefItemList(RefNodeList{RefId{11, "XMAXYMINZMIN"}, boundary_nodes});
   }
 
   {   // xmaxymaxzmin
     Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({cell_size[0], cell_size[1], 0});
+    boundary_nodes[0] = node_number(cell_size[0], cell_size[1], 0);
     descriptor.addRefItemList(RefNodeList{RefId{12, "XMAXYMAXZMIN"}, boundary_nodes});
   }
 
   {   // xminymaxzmin
     Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({0, cell_size[1], 0});
+    boundary_nodes[0] = node_number(0, cell_size[1], 0);
     descriptor.addRefItemList(RefNodeList{RefId{13, "XMINYMAXZMIN"}, boundary_nodes});
   }
 
   {   // xminyminzmax
     Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({0, 0, cell_size[2]});
+    boundary_nodes[0] = node_number(0, 0, cell_size[2]);
     descriptor.addRefItemList(RefNodeList{RefId{14, "XMINYMINZMAX"}, boundary_nodes});
   }
 
   {   // xmaxyminzmax
     Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({cell_size[0], 0, cell_size[2]});
+    boundary_nodes[0] = node_number(cell_size[0], 0, cell_size[2]);
     descriptor.addRefItemList(RefNodeList{RefId{15, "XMAXYMINZMAX"}, boundary_nodes});
   }
 
   {   // xmaxymaxzmax
     Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({cell_size[0], cell_size[1], cell_size[2]});
+    boundary_nodes[0] = node_number(cell_size[0], cell_size[1], cell_size[2]);
     descriptor.addRefItemList(RefNodeList{RefId{16, "XMAXYMAXZMAX"}, boundary_nodes});
   }
 
   {   // xminymaxzmax
     Array<NodeId> boundary_nodes(1);
-    boundary_nodes[0] = node_number({0, cell_size[1], cell_size[2]});
+    boundary_nodes[0] = node_number(0, cell_size[1], cell_size[2]);
     descriptor.addRefItemList(RefNodeList{RefId{17, "XMINYMAXZMAX"}, boundary_nodes});
   }
 }
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 4a3ab6939b02dffef76dbdda5e99273c6653bce9..596ab38bd3713686e3b69b96e7b854b170bca5e0 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -365,8 +365,8 @@ class MeshData : public IMeshData
       NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
       parallel_for(
         m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-          Cjr(j, 0) = -1;
-          Cjr(j, 1) = 1;
+          Cjr(j, 0) = Rd{-1};
+          Cjr(j, 1) = Rd{1};
         });
       m_Cjr = Cjr;
     } else if constexpr (Dimension == 2) {
diff --git a/tests/test_Array.cpp b/tests/test_Array.cpp
index d6382d60699c670033a4ae64cd5970f6dcb86c21..fc5a76735fcb78c81ca7635751f39e85f8d24040 100644
--- a/tests/test_Array.cpp
+++ b/tests/test_Array.cpp
@@ -286,16 +286,16 @@ TEST_CASE("Array", "[utils]")
     {
       using N2 = TinyVector<2, int>;
       Array<N2> b(10);
-      b[0] = {13, 2};
-      b[1] = {1, 3};
-      b[2] = {8, -2};
-      b[3] = {-3, 2};
-      b[4] = {23, 4};
-      b[5] = {-1, -3};
-      b[6] = {13, 17};
-      b[7] = {0, 9};
-      b[8] = {12, 13};
-      b[9] = {9, -17};
+      b[0] = N2{13, 2};
+      b[1] = N2{1, 3};
+      b[2] = N2{8, -2};
+      b[3] = N2{-3, 2};
+      b[4] = N2{23, 4};
+      b[5] = N2{-1, -3};
+      b[6] = N2{13, 17};
+      b[7] = N2{0, 9};
+      b[8] = N2{12, 13};
+      b[9] = N2{9, -17};
 
       REQUIRE((sum(b) == N2{75, 28}));
     }
@@ -304,16 +304,16 @@ TEST_CASE("Array", "[utils]")
     {
       using N22 = TinyMatrix<2, 2, int>;
       Array<N22> b(10);
-      b[0] = {13, 2, 0, 1};
-      b[1] = {1, 3, 6, 3};
-      b[2] = {8, -2, -1, 21};
-      b[3] = {-3, 2, 5, 12};
-      b[4] = {23, 4, 7, 1};
-      b[5] = {-1, -3, 33, 11};
-      b[6] = {13, 17, 12, 13};
-      b[7] = {0, 9, 1, 14};
-      b[8] = {12, 13, -3, -71};
-      b[9] = {9, -17, 0, 16};
+      b[0] = N22{13, 2, 0, 1};
+      b[1] = N22{1, 3, 6, 3};
+      b[2] = N22{8, -2, -1, 21};
+      b[3] = N22{-3, 2, 5, 12};
+      b[4] = N22{23, 4, 7, 1};
+      b[5] = N22{-1, -3, 33, 11};
+      b[6] = N22{13, 17, 12, 13};
+      b[7] = N22{0, 9, 1, 14};
+      b[8] = N22{12, 13, -3, -71};
+      b[9] = N22{9, -17, 0, 16};
 
       REQUIRE((sum(b) == N22{75, 28, 60, 21}));
     }
diff --git a/tests/test_BuiltinFunctionEmbedder.cpp b/tests/test_BuiltinFunctionEmbedder.cpp
index ba88ac0b0a0e3d6e7c3a16977e03437a1e14a05c..51b2462b44a696e9b5214678b17e96d888c6dd27 100644
--- a/tests/test_BuiltinFunctionEmbedder.cpp
+++ b/tests/test_BuiltinFunctionEmbedder.cpp
@@ -89,7 +89,7 @@ TEST_CASE("BuiltinFunctionEmbedder", "[language]")
 
     BuiltinFunctionEmbedder<TinyVector<2>(TinyMatrix<2>, TinyVector<2>)> embedded_c{c};
 
-    TinyMatrix<2> a_arg = {2.3, 1, -2, 3};
+    TinyMatrix<2> a_arg{2.3, 1, -2, 3};
     TinyVector<2> x_arg{3, 2};
 
     std::vector<DataVariant> args;
diff --git a/tests/test_BuiltinFunctionEmbedderUtils.cpp b/tests/test_BuiltinFunctionEmbedderUtils.cpp
index 965d709b97187fbab363b4cacad7cde4f8ad7b7c..da2fbdd1dc8466377621b077af3ef790ec15e24e 100644
--- a/tests/test_BuiltinFunctionEmbedderUtils.cpp
+++ b/tests/test_BuiltinFunctionEmbedderUtils.cpp
@@ -226,7 +226,7 @@ foo(3,1);
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyVector<2>(double, const TinyVector<1>&)>>(
                                             [](double a, const TinyVector<1>& x) -> TinyVector<2> {
-                                              return {a, x[0]};
+                                              return TinyVector<2>{a, x[0]};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
@@ -253,7 +253,7 @@ foo(3.1,f(1));
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyVector<2>(double, const TinyVector<1>&)>>(
                                             [](double a, const TinyVector<1>& x) -> TinyVector<2> {
-                                              return {a, x[0]};
+                                              return TinyVector<2>{a, x[0]};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
@@ -280,7 +280,7 @@ foo(f(2));
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyVector<2>(double, const TinyVector<1>&)>>(
                                             [](double a, const TinyVector<1>& x) -> TinyVector<2> {
-                                              return {a, x[0]};
+                                              return TinyVector<2>{a, x[0]};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
@@ -334,7 +334,7 @@ foo(3,0);
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyVector<3>(double, const TinyVector<2>&)>>(
                                             [](double a, const TinyVector<2>& x) -> TinyVector<3> {
-                                              return {a, x[0], x[1]};
+                                              return TinyVector<3>{a, x[0], x[1]};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
@@ -360,7 +360,7 @@ foo(3.1,(1,2.3));
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyVector<2>(double, const TinyVector<2>&)>>(
                                             [](double a, const TinyVector<2>& x) -> TinyVector<2> {
-                                              return {a * x[1], x[0]};
+                                              return TinyVector<2>{a * x[1], x[0]};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
@@ -390,7 +390,7 @@ foo(3,x);
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyVector<2>(double, const TinyVector<3>&)>>(
                                             [](double a, const TinyVector<3>& x) -> TinyVector<2> {
-                                              return {a * x[0], x[1] + x[2]};
+                                              return TinyVector<2>{a * x[0], x[1] + x[2]};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
@@ -442,7 +442,7 @@ foo(3.1,(1,2.3,4));
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyVector<2>(double, const TinyVector<3>&)>>(
                                             [](double a, const TinyVector<3>& x) -> TinyVector<2> {
-                                              return {a * x[1], x[0] * x[2]};
+                                              return TinyVector<2>{a * x[1], x[0] * x[2]};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
@@ -496,7 +496,7 @@ foo(3,1);
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyMatrix<2>(double, const TinyMatrix<1>&)>>(
                                             [](double a, const TinyMatrix<1>& x) -> TinyMatrix<2> {
-                                              return {a, x(0, 0), a * x(0, 0), x(0, 0)};
+                                              return TinyMatrix<2>{a, x(0, 0), a * x(0, 0), x(0, 0)};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
@@ -523,7 +523,7 @@ foo(3.1,f(1));
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyMatrix<2>(double, const TinyMatrix<1>&)>>(
                                             [](double a, const TinyMatrix<1>& x) -> TinyMatrix<2> {
-                                              return {a, x(0, 0), 2 + x(0, 0), 1};
+                                              return TinyMatrix<2>{a, x(0, 0), 2 + x(0, 0), 1};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
@@ -550,7 +550,7 @@ foo(f(2));
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyMatrix<2>(double, const TinyMatrix<1>&)>>(
                                             [](double a, const TinyMatrix<1>& x) -> TinyMatrix<2> {
-                                              return {a, x(0, 0), a + x(0, 0), 1};
+                                              return TinyMatrix<2>{a, x(0, 0), a + x(0, 0), 1};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
@@ -604,9 +604,9 @@ foo(3,0);
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyMatrix<3>(double, const TinyMatrix<2>&)>>(
                                             [](double a, const TinyMatrix<2>& x) -> TinyMatrix<3> {
-                                              return {a, x(0, 0), x(1, 0),   //
-                                                      a, x(1, 0), x(1, 1),   //
-                                                      a, x(0, 0), x(1, 1)};
+                                              return TinyMatrix<3>{a, x(0, 0), x(1, 0),   //
+                                                                   a, x(1, 0), x(1, 1),   //
+                                                                   a, x(0, 0), x(1, 1)};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
@@ -632,7 +632,7 @@ foo(3.1,(1,2.3,0,3));
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyMatrix<2>(double, const TinyMatrix<2>&)>>(
                                             [](double a, const TinyMatrix<2>& x) -> TinyMatrix<2> {
-                                              return {a * x(0, 0), x(1, 1), a, x(1, 0)};
+                                              return TinyMatrix<2>{a * x(0, 0), x(1, 1), a, x(1, 0)};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
@@ -662,8 +662,8 @@ foo(3,x);
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyMatrix<2>(double, const TinyMatrix<3>&)>>(
                                             [](double a, const TinyMatrix<3>& x) -> TinyMatrix<2> {
-                                              return {a * x(0, 0), x(1, 0) + x(0, 1),   //
-                                                      a, x(1, 1)};
+                                              return TinyMatrix<2>{a * x(0, 0), x(1, 0) + x(0, 1),   //
+                                                                   a, x(1, 1)};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
@@ -715,8 +715,8 @@ foo(3.1,(1, 2.3, 4, 0.3, 2.5, 4.6, 2.7, 8.1, -9));
                                           std::make_shared<
                                             BuiltinFunctionEmbedder<TinyMatrix<2>(double, const TinyMatrix<3>&)>>(
                                             [](double a, const TinyMatrix<3>& x) -> TinyMatrix<2> {
-                                              return {a * x(1, 1), x(1, 0),   //
-                                                      x(0, 1), a * x(0, 0)};
+                                              return TinyMatrix<2>{a * x(1, 1), x(1, 0),   //
+                                                                   x(0, 1), a * x(0, 0)};
                                             }))});
 
       auto function_embedder = getBuiltinFunctionEmbedder(*root_node->children[0]);
diff --git a/tests/test_BuiltinFunctionRegister.hpp b/tests/test_BuiltinFunctionRegister.hpp
index 0ac5226fcac218e5ad195e5f2ecbdef6aa058cb3..49bceabb565ec13b04350f083009c85702886fc2 100644
--- a/tests/test_BuiltinFunctionRegister.hpp
+++ b/tests/test_BuiltinFunctionRegister.hpp
@@ -53,7 +53,7 @@ class test_BuiltinFunctionRegister
 
     m_name_builtin_function_map.insert(
       std::make_pair("RtoR1:R", std::make_shared<BuiltinFunctionEmbedder<TinyVector<1>(double)>>(
-                                  [](double r) -> TinyVector<1> { return {r}; })));
+                                  [](double r) -> TinyVector<1> { return TinyVector<1>{r}; })));
 
     m_name_builtin_function_map.insert(
       std::make_pair("R1toR:R^1", std::make_shared<BuiltinFunctionEmbedder<double(TinyVector<1>)>>(
@@ -74,7 +74,7 @@ class test_BuiltinFunctionRegister
 
     m_name_builtin_function_map.insert(
       std::make_pair("RtoR11:R", std::make_shared<BuiltinFunctionEmbedder<TinyMatrix<1>(double)>>(
-                                   [](double r) -> TinyMatrix<1> { return {r}; })));
+                                   [](double r) -> TinyMatrix<1> { return TinyMatrix<1>{r}; })));
 
     m_name_builtin_function_map.insert(
       std::make_pair("R11toR:R^1x1", std::make_shared<BuiltinFunctionEmbedder<double(TinyMatrix<1>)>>(
diff --git a/tests/test_CubeTransformation.cpp b/tests/test_CubeTransformation.cpp
index 7d060ccd6faf409fb08813607209da1d9da616b1..346dd3c4e9b1f3412bb53503c6880ea07a1508e6 100644
--- a/tests/test_CubeTransformation.cpp
+++ b/tests/test_CubeTransformation.cpp
@@ -15,25 +15,25 @@ TEST_CASE("CubeTransformation", "[geometry]")
 
   SECTION("invertible transformation")
   {
-    const R3 a_hat = {-1, -1, -1};
-    const R3 b_hat = {+1, -1, -1};
-    const R3 c_hat = {+1, +1, -1};
-    const R3 d_hat = {-1, +1, -1};
-    const R3 e_hat = {-1, -1, +1};
-    const R3 f_hat = {+1, -1, +1};
-    const R3 g_hat = {+1, +1, +1};
-    const R3 h_hat = {-1, +1, +1};
+    const R3 a_hat{-1, -1, -1};
+    const R3 b_hat{+1, -1, -1};
+    const R3 c_hat{+1, +1, -1};
+    const R3 d_hat{-1, +1, -1};
+    const R3 e_hat{-1, -1, +1};
+    const R3 f_hat{+1, -1, +1};
+    const R3 g_hat{+1, +1, +1};
+    const R3 h_hat{-1, +1, +1};
 
     const R3 m_hat = zero;
 
-    const R3 a = {0, 0, 0};
-    const R3 b = {3, 1, 3};
-    const R3 c = {2, 5, 2};
-    const R3 d = {0, 3, 1};
-    const R3 e = {1, 2, 5};
-    const R3 f = {3, 1, 7};
-    const R3 g = {2, 5, 5};
-    const R3 h = {0, 3, 6};
+    const R3 a{0, 0, 0};
+    const R3 b{3, 1, 3};
+    const R3 c{2, 5, 2};
+    const R3 d{0, 3, 1};
+    const R3 e{1, 2, 5};
+    const R3 f{3, 1, 7};
+    const R3 g{2, 5, 5};
+    const R3 h{0, 3, 6};
 
     const R3 m = 0.125 * (a + b + c + d + e + f + g + h);
 
@@ -153,10 +153,10 @@ TEST_CASE("CubeTransformation", "[geometry]")
   {
     using R3 = TinyVector<3>;
 
-    const R3 a = {1, 2, 1};
-    const R3 b = {3, 1, 3};
-    const R3 c = {2, 5, 2};
-    const R3 d = {2, 3, 4};
+    const R3 a{1, 2, 1};
+    const R3 b{3, 1, 3};
+    const R3 c{2, 5, 2};
+    const R3 d{2, 3, 4};
 
     const CubeTransformation t(a, b, c, c, d, d, d, d);
 
@@ -209,12 +209,12 @@ TEST_CASE("CubeTransformation", "[geometry]")
 
   SECTION("degenerate to prism")
   {
-    const R3 a = {1, 2, 0};
-    const R3 b = {3, 1, 3};
-    const R3 c = {2, 5, 2};
-    const R3 d = {0, 3, 1};
-    const R3 e = {1, 2, 5};
-    const R3 f = {3, 1, 7};
+    const R3 a{1, 2, 0};
+    const R3 b{3, 1, 3};
+    const R3 c{2, 5, 2};
+    const R3 d{0, 3, 1};
+    const R3 e{1, 2, 5};
+    const R3 f{3, 1, 7};
 
     const CubeTransformation t(a, b, c, c, d, e, f, f);
 
@@ -269,11 +269,11 @@ TEST_CASE("CubeTransformation", "[geometry]")
 
   SECTION("degenerate to pyramid")
   {
-    const R3 a = {1, 2, 0};
-    const R3 b = {3, 1, 3};
-    const R3 c = {2, 5, 2};
-    const R3 d = {0, 3, 1};
-    const R3 e = {1, 2, 5};
+    const R3 a{1, 2, 0};
+    const R3 b{3, 1, 3};
+    const R3 c{2, 5, 2};
+    const R3 d{0, 3, 1};
+    const R3 e{1, 2, 5};
 
     const CubeTransformation t(a, b, c, d, e, e, e, e);
 
diff --git a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
index 251424f6a7340244b6ac60879cc85653b3115351..a7f0cc0d74eb79f0897e2ea2a2a4b7212dbc306e 100644
--- a/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
+++ b/tests/test_EmbeddedIDiscreteFunctionMathFunctions.cpp
@@ -175,11 +175,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
 
         constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
           if constexpr (Dimension == 1) {
-            return {x[0], 1 + x[0] * x[0]};
+            return TinyVector<2>{x[0], 1 + x[0] * x[0]};
           } else if constexpr (Dimension == 2) {
-            return {x[0], x[1]};
+            return TinyVector<2>{x[0], x[1]};
           } else if constexpr (Dimension == 3) {
-            return {x[0], x[1] + x[2]};
+            return TinyVector<2>{x[0], x[1] + x[2]};
           }
         };
 
@@ -188,7 +188,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           parallel_for(
             uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<2> x = to_2d(xj[cell_id]);
-              uj[cell_id]           = {2 * x[0] + 1, 1 - x[1]};
+              uj[cell_id]           = TinyVector<2>{2 * x[0] + 1, 1 - x[1]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
@@ -199,7 +199,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           parallel_for(
             vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<2> x = to_2d(xj[cell_id]);
-              vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1]};
+              vj[cell_id]           = TinyVector<2>{x[0] * x[1] + 1, 2 * x[1]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, vj);
@@ -210,11 +210,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
 
         constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
           if constexpr (Dimension == 1) {
-            return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+            return TinyVector<3>{x[0], 1 + x[0] * x[0], 2 - x[0]};
           } else if constexpr (Dimension == 2) {
-            return {x[0], x[1], x[0] + x[1]};
+            return TinyVector<3>{x[0], x[1], x[0] + x[1]};
           } else if constexpr (Dimension == 3) {
-            return {x[0], x[1], x[2]};
+            return TinyVector<3>{x[0], x[1], x[2]};
           }
         };
 
@@ -223,7 +223,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           parallel_for(
             uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<3> x = to_3d(xj[cell_id]);
-              uj[cell_id]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+              uj[cell_id]           = TinyVector<3>{2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
@@ -234,7 +234,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           parallel_for(
             vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<3> x = to_3d(xj[cell_id]);
-              vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
+              vj[cell_id]           = TinyVector<3>{x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, vj);
@@ -246,7 +246,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
         std::shared_ptr p_R1x1_u = [=] {
           CellValue<TinyMatrix<1>> uj{mesh->connectivity()};
           parallel_for(
-            uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = {2 * xj[cell_id][0] + 1}; });
+            uj.numberOfItems(),
+            PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = TinyMatrix<1>{2 * xj[cell_id][0] + 1}; });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
         }();
@@ -257,8 +258,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
             uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<2> x = to_2d(xj[cell_id]);
 
-              uj[cell_id] = {2 * x[0] + 1, 1 - x[1],   //
-                             2 * x[1], -x[0]};
+              uj[cell_id] = TinyMatrix<2>{2 * x[0] + 1, 1 - x[1],   //
+                                          2 * x[1], -x[0]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
@@ -270,9 +271,9 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
             uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<3> x = to_3d(xj[cell_id]);
 
-              uj[cell_id] = {2 * x[0] + 1,    1 - x[1],        3,             //
-                             2 * x[1],        -x[0],           x[0] - x[1],   //
-                             3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
+              uj[cell_id] = TinyMatrix<3>{2 * x[0] + 1,    1 - x[1],        3,             //
+                                          2 * x[1],        -x[0],           x[0] - x[1],   //
+                                          3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
@@ -684,11 +685,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
 
         constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
           if constexpr (Dimension == 1) {
-            return {x[0], 1 + x[0] * x[0]};
+            return TinyVector<2>{x[0], 1 + x[0] * x[0]};
           } else if constexpr (Dimension == 2) {
-            return {x[0], x[1]};
+            return TinyVector<2>{x[0], x[1]};
           } else if constexpr (Dimension == 3) {
-            return {x[0], x[1] + x[2]};
+            return TinyVector<2>{x[0], x[1] + x[2]};
           }
         };
 
@@ -697,7 +698,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           parallel_for(
             uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<2> x = to_2d(xj[cell_id]);
-              uj[cell_id]           = {2 * x[0] + 1, 1 - x[1]};
+              uj[cell_id]           = TinyVector<2>{2 * x[0] + 1, 1 - x[1]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
@@ -708,7 +709,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           parallel_for(
             vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<2> x = to_2d(xj[cell_id]);
-              vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1]};
+              vj[cell_id]           = TinyVector<2>{x[0] * x[1] + 1, 2 * x[1]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, vj);
@@ -719,11 +720,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
 
         constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
           if constexpr (Dimension == 1) {
-            return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+            return TinyVector<3>{x[0], 1 + x[0] * x[0], 2 - x[0]};
           } else if constexpr (Dimension == 2) {
-            return {x[0], x[1], x[0] + x[1]};
+            return TinyVector<3>{x[0], x[1], x[0] + x[1]};
           } else if constexpr (Dimension == 3) {
-            return {x[0], x[1], x[2]};
+            return TinyVector<3>{x[0], x[1], x[2]};
           }
         };
 
@@ -732,7 +733,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           parallel_for(
             uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<3> x = to_3d(xj[cell_id]);
-              uj[cell_id]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+              uj[cell_id]           = TinyVector<3>{2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
@@ -743,7 +744,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           parallel_for(
             vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<3> x = to_3d(xj[cell_id]);
-              vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
+              vj[cell_id]           = TinyVector<3>{x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, vj);
@@ -755,7 +756,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
         std::shared_ptr p_R1x1_u = [=] {
           CellValue<TinyMatrix<1>> uj{mesh->connectivity()};
           parallel_for(
-            uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = {2 * xj[cell_id][0] + 1}; });
+            uj.numberOfItems(),
+            PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = TinyMatrix<1>{2 * xj[cell_id][0] + 1}; });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
         }();
@@ -766,8 +768,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
             uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<2> x = to_2d(xj[cell_id]);
 
-              uj[cell_id] = {2 * x[0] + 1, 1 - x[1],   //
-                             2 * x[1], -x[0]};
+              uj[cell_id] = TinyMatrix<2>{2 * x[0] + 1, 1 - x[1],   //
+                                          2 * x[1], -x[0]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
@@ -779,9 +781,9 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
             uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<3> x = to_3d(xj[cell_id]);
 
-              uj[cell_id] = {2 * x[0] + 1,    1 - x[1],        3,             //
-                             2 * x[1],        -x[0],           x[0] - x[1],   //
-                             3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
+              uj[cell_id] = TinyMatrix<3>{2 * x[0] + 1,    1 - x[1],        3,             //
+                                          2 * x[1],        -x[0],           x[0] - x[1],   //
+                                          3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
@@ -1193,11 +1195,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
 
         constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
           if constexpr (Dimension == 1) {
-            return {x[0], 1 + x[0] * x[0]};
+            return TinyVector<2>{x[0], 1 + x[0] * x[0]};
           } else if constexpr (Dimension == 2) {
-            return {x[0], x[1]};
+            return TinyVector<2>{x[0], x[1]};
           } else if constexpr (Dimension == 3) {
-            return {x[0], x[1] + x[2]};
+            return TinyVector<2>{x[0], x[1] + x[2]};
           }
         };
 
@@ -1206,7 +1208,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           parallel_for(
             uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<2> x = to_2d(xj[cell_id]);
-              uj[cell_id]           = {2 * x[0] + 1, 1 - x[1]};
+              uj[cell_id]           = TinyVector<2>{2 * x[0] + 1, 1 - x[1]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
@@ -1217,7 +1219,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           parallel_for(
             vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<2> x = to_2d(xj[cell_id]);
-              vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1]};
+              vj[cell_id]           = TinyVector<2>{x[0] * x[1] + 1, 2 * x[1]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, vj);
@@ -1228,11 +1230,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
 
         constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
           if constexpr (Dimension == 1) {
-            return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+            return TinyVector<3>{x[0], 1 + x[0] * x[0], 2 - x[0]};
           } else if constexpr (Dimension == 2) {
-            return {x[0], x[1], x[0] + x[1]};
+            return TinyVector<3>{x[0], x[1], x[0] + x[1]};
           } else if constexpr (Dimension == 3) {
-            return {x[0], x[1], x[2]};
+            return TinyVector<3>{x[0], x[1], x[2]};
           }
         };
 
@@ -1241,7 +1243,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           parallel_for(
             uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<3> x = to_3d(xj[cell_id]);
-              uj[cell_id]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+              uj[cell_id]           = TinyVector<3>{2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
@@ -1252,7 +1254,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
           parallel_for(
             vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<3> x = to_3d(xj[cell_id]);
-              vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
+              vj[cell_id]           = TinyVector<3>{x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, vj);
@@ -1264,7 +1266,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
         std::shared_ptr p_R1x1_u = [=] {
           CellValue<TinyMatrix<1>> uj{mesh->connectivity()};
           parallel_for(
-            uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = {2 * xj[cell_id][0] + 1}; });
+            uj.numberOfItems(),
+            PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = TinyMatrix<1>{2 * xj[cell_id][0] + 1}; });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
         }();
@@ -1275,8 +1278,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
             uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<2> x = to_2d(xj[cell_id]);
 
-              uj[cell_id] = {2 * x[0] + 1, 1 - x[1],   //
-                             2 * x[1], -x[0]};
+              uj[cell_id] = TinyMatrix<2>{2 * x[0] + 1, 1 - x[1],   //
+                                          2 * x[1], -x[0]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
@@ -1288,9 +1291,9 @@ TEST_CASE("EmbeddedIDiscreteFunctionMathFunctions", "[scheme]")
             uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
               const TinyVector<3> x = to_3d(xj[cell_id]);
 
-              uj[cell_id] = {2 * x[0] + 1,    1 - x[1],        3,             //
-                             2 * x[1],        -x[0],           x[0] - x[1],   //
-                             3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
+              uj[cell_id] = TinyMatrix<3>{2 * x[0] + 1,    1 - x[1],        3,             //
+                                          2 * x[1],        -x[0],           x[0] - x[1],   //
+                                          3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
             });
 
           return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
diff --git a/tests/test_EmbeddedIDiscreteFunctionOperators.cpp b/tests/test_EmbeddedIDiscreteFunctionOperators.cpp
index 5d91447d90147537c3428addfea8c2710defed84..b955ba24b785a83c3be6f73360fa9ea0026cc56e 100644
--- a/tests/test_EmbeddedIDiscreteFunctionOperators.cpp
+++ b/tests/test_EmbeddedIDiscreteFunctionOperators.cpp
@@ -254,11 +254,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
 
           constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
             if constexpr (Dimension == 1) {
-              return {x[0], 1 + x[0] * x[0]};
+              return TinyVector<2>{x[0], 1 + x[0] * x[0]};
             } else if constexpr (Dimension == 2) {
-              return {x[0], x[1]};
+              return TinyVector<2>{x[0], x[1]};
             } else if constexpr (Dimension == 3) {
-              return {x[0], x[1] + x[2]};
+              return TinyVector<2>{x[0], x[1] + x[2]};
             }
           };
 
@@ -267,7 +267,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
-                uj[cell_id]           = {2 * x[0] + 1, 1 - x[1]};
+                uj[cell_id]           = TinyVector<2>{2 * x[0] + 1, 1 - x[1]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
@@ -278,7 +278,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
-                vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1]};
+                vj[cell_id]           = TinyVector<2>{x[0] * x[1] + 1, 2 * x[1]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, vj);
@@ -289,11 +289,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
 
           constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
             if constexpr (Dimension == 1) {
-              return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+              return TinyVector<3>{x[0], 1 + x[0] * x[0], 2 - x[0]};
             } else if constexpr (Dimension == 2) {
-              return {x[0], x[1], x[0] + x[1]};
+              return TinyVector<3>{x[0], x[1], x[0] + x[1]};
             } else if constexpr (Dimension == 3) {
-              return {x[0], x[1], x[2]};
+              return TinyVector<3>{x[0], x[1], x[2]};
             }
           };
 
@@ -302,7 +302,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
-                uj[cell_id]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+                uj[cell_id]           = TinyVector<3>{2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
@@ -313,7 +313,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
-                vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
+                vj[cell_id]           = TinyVector<3>{x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, vj);
@@ -325,7 +325,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
           std::shared_ptr p_R1x1_u = [=] {
             CellValue<TinyMatrix<1>> uj{mesh->connectivity()};
             parallel_for(
-              uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = {2 * xj[cell_id][0] + 1}; });
+              uj.numberOfItems(),
+              PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = TinyMatrix<1>{2 * xj[cell_id][0] + 1}; });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
           }();
@@ -336,7 +337,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
           std::shared_ptr p_R1x1_v = [=] {
             CellValue<TinyMatrix<1>> vj{mesh->connectivity()};
             parallel_for(
-              vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { vj[cell_id] = {0.3 - xj[cell_id][0]}; });
+              vj.numberOfItems(),
+              PUGS_LAMBDA(const CellId cell_id) { vj[cell_id] = TinyMatrix<1>{0.3 - xj[cell_id][0]}; });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, vj);
           }();
@@ -347,8 +349,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
 
-                uj[cell_id] = {2 * x[0] + 1, 1 - x[1],   //
-                               2 * x[1], -x[0]};
+                uj[cell_id] = TinyMatrix<2>{2 * x[0] + 1, 1 - x[1],   //
+                                            2 * x[1], -x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
@@ -363,8 +365,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
 
-                vj[cell_id] = {x[0] + 0.3, 1 - x[1] - x[0],   //
-                               2 * x[1] + x[0], x[1] - x[0]};
+                vj[cell_id] = TinyMatrix<2>{x[0] + 0.3, 1 - x[1] - x[0],   //
+                                            2 * x[1] + x[0], x[1] - x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, vj);
@@ -376,9 +378,9 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
 
-                uj[cell_id] = {2 * x[0] + 1,    1 - x[1],        3,             //
-                               2 * x[1],        -x[0],           x[0] - x[1],   //
-                               3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
+                uj[cell_id] = TinyMatrix<3>{2 * x[0] + 1,    1 - x[1],        3,             //
+                                            2 * x[1],        -x[0],           x[0] - x[1],   //
+                                            3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
@@ -393,9 +395,9 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
 
-                vj[cell_id] = {0.2 * x[0] + 1,  2 + x[1],          3 - x[2],      //
-                               2.3 * x[2],      x[1] - x[0],       x[2] - x[1],   //
-                               2 * x[2] + x[0], x[1] + 0.2 * x[2], x[2] - 2 * x[0]};
+                vj[cell_id] = TinyMatrix<3>{0.2 * x[0] + 1,  2 + x[1],          3 - x[2],      //
+                                            2.3 * x[2],      x[1] - x[0],       x[2] - x[1],   //
+                                            2 * x[2] + x[0], x[1] + 0.2 * x[2], x[2] - 2 * x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, vj);
@@ -951,11 +953,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
 
           constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
             if constexpr (Dimension == 1) {
-              return {x[0], 1 + x[0] * x[0]};
+              return TinyVector<2>{x[0], 1 + x[0] * x[0]};
             } else if constexpr (Dimension == 2) {
-              return {x[0], x[1]};
+              return TinyVector<2>{x[0], x[1]};
             } else if constexpr (Dimension == 3) {
-              return {x[0], x[1] + x[2]};
+              return TinyVector<2>{x[0], x[1] + x[2]};
             }
           };
 
@@ -964,7 +966,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
-                uj[cell_id]           = {2 * x[0] + 1, 1 - x[1]};
+                uj[cell_id]           = TinyVector<2>{2 * x[0] + 1, 1 - x[1]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
@@ -975,7 +977,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
-                vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1]};
+                vj[cell_id]           = TinyVector<2>{x[0] * x[1] + 1, 2 * x[1]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, vj);
@@ -986,11 +988,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
 
           constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
             if constexpr (Dimension == 1) {
-              return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+              return TinyVector<3>{x[0], 1 + x[0] * x[0], 2 - x[0]};
             } else if constexpr (Dimension == 2) {
-              return {x[0], x[1], x[0] + x[1]};
+              return TinyVector<3>{x[0], x[1], x[0] + x[1]};
             } else if constexpr (Dimension == 3) {
-              return {x[0], x[1], x[2]};
+              return TinyVector<3>{x[0], x[1], x[2]};
             }
           };
 
@@ -999,7 +1001,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
-                uj[cell_id]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+                uj[cell_id]           = TinyVector<3>{2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
@@ -1010,7 +1012,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
-                vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
+                vj[cell_id]           = TinyVector<3>{x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, vj);
@@ -1022,7 +1024,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
           std::shared_ptr p_R1x1_u = [=] {
             CellValue<TinyMatrix<1>> uj{mesh->connectivity()};
             parallel_for(
-              uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = {2 * xj[cell_id][0] + 1}; });
+              uj.numberOfItems(),
+              PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = TinyMatrix<1>{2 * xj[cell_id][0] + 1}; });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
           }();
@@ -1033,7 +1036,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
           std::shared_ptr p_R1x1_v = [=] {
             CellValue<TinyMatrix<1>> vj{mesh->connectivity()};
             parallel_for(
-              vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { vj[cell_id] = {0.3 - xj[cell_id][0]}; });
+              vj.numberOfItems(),
+              PUGS_LAMBDA(const CellId cell_id) { vj[cell_id] = TinyMatrix<1>{0.3 - xj[cell_id][0]}; });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, vj);
           }();
@@ -1044,8 +1048,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
 
-                uj[cell_id] = {2 * x[0] + 1, 1 - x[1],   //
-                               2 * x[1], -x[0]};
+                uj[cell_id] = TinyMatrix<2>{2 * x[0] + 1, 1 - x[1],   //
+                                            2 * x[1], -x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
@@ -1060,8 +1064,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
 
-                vj[cell_id] = {x[0] + 0.3, 1 - x[1] - x[0],   //
-                               2 * x[1] + x[0], x[1] - x[0]};
+                vj[cell_id] = TinyMatrix<2>{x[0] + 0.3, 1 - x[1] - x[0],   //
+                                            2 * x[1] + x[0], x[1] - x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, vj);
@@ -1073,9 +1077,9 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
 
-                uj[cell_id] = {2 * x[0] + 1,    1 - x[1],        3,             //
-                               2 * x[1],        -x[0],           x[0] - x[1],   //
-                               3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
+                uj[cell_id] = TinyMatrix<3>{2 * x[0] + 1,    1 - x[1],        3,             //
+                                            2 * x[1],        -x[0],           x[0] - x[1],   //
+                                            3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
@@ -1090,9 +1094,9 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
 
-                vj[cell_id] = {0.2 * x[0] + 1,  2 + x[1],          3 - x[2],      //
-                               2.3 * x[2],      x[1] - x[0],       x[2] - x[1],   //
-                               2 * x[2] + x[0], x[1] + 0.2 * x[2], x[2] - 2 * x[0]};
+                vj[cell_id] = TinyMatrix<3>{0.2 * x[0] + 1,  2 + x[1],          3 - x[2],      //
+                                            2.3 * x[2],      x[1] - x[0],       x[2] - x[1],   //
+                                            2 * x[2] + x[0], x[1] + 0.2 * x[2], x[2] - 2 * x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, vj);
@@ -1648,11 +1652,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
 
           constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
             if constexpr (Dimension == 1) {
-              return {x[0], 1 + x[0] * x[0]};
+              return TinyVector<2>{x[0], 1 + x[0] * x[0]};
             } else if constexpr (Dimension == 2) {
-              return {x[0], x[1]};
+              return TinyVector<2>{x[0], x[1]};
             } else if constexpr (Dimension == 3) {
-              return {x[0], x[1] + x[2]};
+              return TinyVector<2>{x[0], x[1] + x[2]};
             }
           };
 
@@ -1661,7 +1665,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
-                uj[cell_id]           = {2 * x[0] + 1, 1 - x[1]};
+                uj[cell_id]           = TinyVector<2>{2 * x[0] + 1, 1 - x[1]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
@@ -1672,7 +1676,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
-                vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1]};
+                vj[cell_id]           = TinyVector<2>{x[0] * x[1] + 1, 2 * x[1]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, vj);
@@ -1683,11 +1687,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
 
           constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
             if constexpr (Dimension == 1) {
-              return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+              return TinyVector<3>{x[0], 1 + x[0] * x[0], 2 - x[0]};
             } else if constexpr (Dimension == 2) {
-              return {x[0], x[1], x[0] + x[1]};
+              return TinyVector<3>{x[0], x[1], x[0] + x[1]};
             } else if constexpr (Dimension == 3) {
-              return {x[0], x[1], x[2]};
+              return TinyVector<3>{x[0], x[1], x[2]};
             }
           };
 
@@ -1696,7 +1700,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
-                uj[cell_id]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+                uj[cell_id]           = TinyVector<3>{2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
@@ -1707,7 +1711,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
-                vj[cell_id]           = {x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
+                vj[cell_id]           = TinyVector<3>{x[0] * x[1] + 1, 2 * x[1], x[2] * x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, vj);
@@ -1719,7 +1723,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
           std::shared_ptr p_R1x1_u = [=] {
             CellValue<TinyMatrix<1>> uj{mesh->connectivity()};
             parallel_for(
-              uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = {2 * xj[cell_id][0] + 1}; });
+              uj.numberOfItems(),
+              PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = TinyMatrix<1>{2 * xj[cell_id][0] + 1}; });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
           }();
@@ -1730,7 +1735,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
           std::shared_ptr p_R1x1_v = [=] {
             CellValue<TinyMatrix<1>> vj{mesh->connectivity()};
             parallel_for(
-              vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { vj[cell_id] = {0.3 - xj[cell_id][0]}; });
+              vj.numberOfItems(),
+              PUGS_LAMBDA(const CellId cell_id) { vj[cell_id] = TinyMatrix<1>{0.3 - xj[cell_id][0]}; });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, vj);
           }();
@@ -1741,8 +1747,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
 
-                uj[cell_id] = {2 * x[0] + 1, 1 - x[1],   //
-                               2 * x[1], -x[0]};
+                uj[cell_id] = TinyMatrix<2>{2 * x[0] + 1, 1 - x[1],   //
+                                            2 * x[1], -x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
@@ -1757,8 +1763,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
 
-                vj[cell_id] = {x[0] + 0.3, 1 - x[1] - x[0],   //
-                               2 * x[1] + x[0], x[1] - x[0]};
+                vj[cell_id] = TinyMatrix<2>{x[0] + 0.3, 1 - x[1] - x[0],   //
+                                            2 * x[1] + x[0], x[1] - x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, vj);
@@ -1770,9 +1776,9 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
 
-                uj[cell_id] = {2 * x[0] + 1,    1 - x[1],        3,             //
-                               2 * x[1],        -x[0],           x[0] - x[1],   //
-                               3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
+                uj[cell_id] = TinyMatrix<3>{2 * x[0] + 1,    1 - x[1],        3,             //
+                                            2 * x[1],        -x[0],           x[0] - x[1],   //
+                                            3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
@@ -1787,9 +1793,9 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               vj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
 
-                vj[cell_id] = {0.2 * x[0] + 1,  2 + x[1],          3 - x[2],      //
-                               2.3 * x[2],      x[1] - x[0],       x[2] - x[1],   //
-                               2 * x[2] + x[0], x[1] + 0.2 * x[2], x[2] - 2 * x[0]};
+                vj[cell_id] = TinyMatrix<3>{0.2 * x[0] + 1,  2 + x[1],          3 - x[2],      //
+                                            2.3 * x[2],      x[1] - x[0],       x[2] - x[1],   //
+                                            2 * x[2] + x[0], x[1] + 0.2 * x[2], x[2] - 2 * x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, vj);
@@ -2322,11 +2328,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
 
           constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
             if constexpr (Dimension == 1) {
-              return {x[0], 1 + x[0] * x[0]};
+              return TinyVector<2>{x[0], 1 + x[0] * x[0]};
             } else if constexpr (Dimension == 2) {
-              return {x[0], x[1]};
+              return TinyVector<2>{x[0], x[1]};
             } else if constexpr (Dimension == 3) {
-              return {x[0], x[1] + x[2]};
+              return TinyVector<2>{x[0], x[1] + x[2]};
             }
           };
 
@@ -2335,7 +2341,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
-                uj[cell_id]           = {2 * x[0] + 1, 1 - x[1]};
+                uj[cell_id]           = TinyVector<2>{2 * x[0] + 1, 1 - x[1]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
@@ -2343,11 +2349,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
 
           constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
             if constexpr (Dimension == 1) {
-              return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+              return TinyVector<3>{x[0], 1 + x[0] * x[0], 2 - x[0]};
             } else if constexpr (Dimension == 2) {
-              return {x[0], x[1], x[0] + x[1]};
+              return TinyVector<3>{x[0], x[1], x[0] + x[1]};
             } else if constexpr (Dimension == 3) {
-              return {x[0], x[1], x[2]};
+              return TinyVector<3>{x[0], x[1], x[2]};
             }
           };
 
@@ -2356,7 +2362,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
-                uj[cell_id]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+                uj[cell_id]           = TinyVector<3>{2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
@@ -2365,7 +2371,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
           std::shared_ptr p_R1x1_u = [=] {
             CellValue<TinyMatrix<1>> uj{mesh->connectivity()};
             parallel_for(
-              uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = {2 * xj[cell_id][0] + 1}; });
+              uj.numberOfItems(),
+              PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = TinyMatrix<1>{2 * xj[cell_id][0] + 1}; });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
           }();
@@ -2376,8 +2383,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
 
-                uj[cell_id] = {2 * x[0] + 1, 1 - x[1],   //
-                               2 * x[1], -x[0]};
+                uj[cell_id] = TinyMatrix<2>{2 * x[0] + 1, 1 - x[1],   //
+                                            2 * x[1], -x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
@@ -2389,9 +2396,9 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
 
-                uj[cell_id] = {2 * x[0] + 1,    1 - x[1],        3,             //
-                               2 * x[1],        -x[0],           x[0] - x[1],   //
-                               3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
+                uj[cell_id] = TinyMatrix<3>{2 * x[0] + 1,    1 - x[1],        3,             //
+                                            2 * x[1],        -x[0],           x[0] - x[1],   //
+                                            3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
@@ -2466,11 +2473,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
 
           constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
             if constexpr (Dimension == 1) {
-              return {x[0], 1 + x[0] * x[0]};
+              return TinyVector<2>{x[0], 1 + x[0] * x[0]};
             } else if constexpr (Dimension == 2) {
-              return {x[0], x[1]};
+              return TinyVector<2>{x[0], x[1]};
             } else if constexpr (Dimension == 3) {
-              return {x[0], x[1] + x[2]};
+              return TinyVector<2>{x[0], x[1] + x[2]};
             }
           };
 
@@ -2479,7 +2486,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
-                uj[cell_id]           = {2 * x[0] + 1, 1 - x[1]};
+                uj[cell_id]           = TinyVector<2>{2 * x[0] + 1, 1 - x[1]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
@@ -2487,11 +2494,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
 
           constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
             if constexpr (Dimension == 1) {
-              return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+              return TinyVector<3>{x[0], 1 + x[0] * x[0], 2 - x[0]};
             } else if constexpr (Dimension == 2) {
-              return {x[0], x[1], x[0] + x[1]};
+              return TinyVector<3>{x[0], x[1], x[0] + x[1]};
             } else if constexpr (Dimension == 3) {
-              return {x[0], x[1], x[2]};
+              return TinyVector<3>{x[0], x[1], x[2]};
             }
           };
 
@@ -2500,7 +2507,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
-                uj[cell_id]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+                uj[cell_id]           = TinyVector<3>{2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
@@ -2509,7 +2516,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
           std::shared_ptr p_R1x1_u = [=] {
             CellValue<TinyMatrix<1>> uj{mesh->connectivity()};
             parallel_for(
-              uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = {2 * xj[cell_id][0] + 1}; });
+              uj.numberOfItems(),
+              PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = TinyMatrix<1>{2 * xj[cell_id][0] + 1}; });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
           }();
@@ -2520,8 +2528,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
 
-                uj[cell_id] = {2 * x[0] + 1, 1 - x[1],   //
-                               2 * x[1], -x[0]};
+                uj[cell_id] = TinyMatrix<2>{2 * x[0] + 1, 1 - x[1],   //
+                                            2 * x[1], -x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
@@ -2533,9 +2541,9 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
 
-                uj[cell_id] = {2 * x[0] + 1,    1 - x[1],        3,             //
-                               2 * x[1],        -x[0],           x[0] - x[1],   //
-                               3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
+                uj[cell_id] = TinyMatrix<3>{2 * x[0] + 1,    1 - x[1],        3,             //
+                                            2 * x[1],        -x[0],           x[0] - x[1],   //
+                                            3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
@@ -2610,11 +2618,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
 
           constexpr auto to_2d = [&](const TinyVector<Dimension>& x) -> TinyVector<2> {
             if constexpr (Dimension == 1) {
-              return {x[0], 1 + x[0] * x[0]};
+              return TinyVector<2>{x[0], 1 + x[0] * x[0]};
             } else if constexpr (Dimension == 2) {
-              return {x[0], x[1]};
+              return TinyVector<2>{x[0], x[1]};
             } else if constexpr (Dimension == 3) {
-              return {x[0], x[1] + x[2]};
+              return TinyVector<2>{x[0], x[1] + x[2]};
             }
           };
 
@@ -2623,7 +2631,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
-                uj[cell_id]           = {2 * x[0] + 1, 1 - x[1]};
+                uj[cell_id]           = TinyVector<2>{2 * x[0] + 1, 1 - x[1]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(mesh, uj);
@@ -2631,11 +2639,11 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
 
           constexpr auto to_3d = [&](const TinyVector<Dimension>& x) -> TinyVector<3> {
             if constexpr (Dimension == 1) {
-              return {x[0], 1 + x[0] * x[0], 2 - x[0]};
+              return TinyVector<3>{x[0], 1 + x[0] * x[0], 2 - x[0]};
             } else if constexpr (Dimension == 2) {
-              return {x[0], x[1], x[0] + x[1]};
+              return TinyVector<3>{x[0], x[1], x[0] + x[1]};
             } else if constexpr (Dimension == 3) {
-              return {x[0], x[1], x[2]};
+              return TinyVector<3>{x[0], x[1], x[2]};
             }
           };
 
@@ -2644,7 +2652,7 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
             parallel_for(
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
-                uj[cell_id]           = {2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
+                uj[cell_id]           = TinyVector<3>{2 * x[0] + 1, 1 - x[1] * x[2], x[0] + x[2]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(mesh, uj);
@@ -2653,7 +2661,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
           std::shared_ptr p_R1x1_u = [=] {
             CellValue<TinyMatrix<1>> uj{mesh->connectivity()};
             parallel_for(
-              uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = {2 * xj[cell_id][0] + 1}; });
+              uj.numberOfItems(),
+              PUGS_LAMBDA(const CellId cell_id) { uj[cell_id] = TinyMatrix<1>{2 * xj[cell_id][0] + 1}; });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(mesh, uj);
           }();
@@ -2664,8 +2673,8 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<2> x = to_2d(xj[cell_id]);
 
-                uj[cell_id] = {2 * x[0] + 1, 1 - x[1],   //
-                               2 * x[1], -x[0]};
+                uj[cell_id] = TinyMatrix<2>{2 * x[0] + 1, 1 - x[1],   //
+                                            2 * x[1], -x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(mesh, uj);
@@ -2677,9 +2686,9 @@ TEST_CASE("EmbeddedIDiscreteFunctionOperators", "[scheme]")
               uj.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
                 const TinyVector<3> x = to_3d(xj[cell_id]);
 
-                uj[cell_id] = {2 * x[0] + 1,    1 - x[1],        3,             //
-                               2 * x[1],        -x[0],           x[0] - x[1],   //
-                               3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
+                uj[cell_id] = TinyMatrix<3>{2 * x[0] + 1,    1 - x[1],        3,             //
+                                            2 * x[1],        -x[0],           x[0] - x[1],   //
+                                            3 * x[2] - x[1], x[1] - 2 * x[2], x[2] - x[0]};
               });
 
             return std::make_shared<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(mesh, uj);
diff --git a/tests/test_LineTransformation.cpp b/tests/test_LineTransformation.cpp
index 908d2d693973d13d48288cda7f7d021ee187c7a7..82559ac50cb17e77ff6916c5d066685649c69986 100644
--- a/tests/test_LineTransformation.cpp
+++ b/tests/test_LineTransformation.cpp
@@ -13,14 +13,14 @@ TEST_CASE("LineTransformation", "[geometry]")
   {
     using R1 = TinyVector<1>;
 
-    const R1 a = 0.3;
-    const R1 b = 2.7;
+    const R1 a{0.3};
+    const R1 b{2.7};
 
     const LineTransformation<1> t(a, b);
 
-    REQUIRE(t(-1)[0] == Catch::Approx(0.3));
-    REQUIRE(t(0)[0] == Catch::Approx(1.5));
-    REQUIRE(t(1)[0] == Catch::Approx(2.7));
+    REQUIRE(t(R1{-1})[0] == Catch::Approx(0.3));
+    REQUIRE(t(R1{0})[0] == Catch::Approx(1.5));
+    REQUIRE(t(R1{1})[0] == Catch::Approx(2.7));
 
     REQUIRE(t.jacobianDeterminant() == Catch::Approx(1.2));
 
@@ -46,19 +46,20 @@ TEST_CASE("LineTransformation", "[geometry]")
 
   SECTION("2D")
   {
+    using R1 = TinyVector<1>;
     using R2 = TinyVector<2>;
 
-    const R2 a = {0.3, 1.2};
-    const R2 b = {2.7, 0.7};
+    const R2 a{0.3, 1.2};
+    const R2 b{2.7, 0.7};
 
     const LineTransformation<2> t(a, b);
 
-    REQUIRE(t(-1)[0] == Catch::Approx(0.3));
-    REQUIRE(t(-1)[1] == Catch::Approx(1.2));
-    REQUIRE(t(0)[0] == Catch::Approx(1.5));
-    REQUIRE(t(0)[1] == Catch::Approx(0.95));
-    REQUIRE(t(1)[0] == Catch::Approx(2.7));
-    REQUIRE(t(1)[1] == Catch::Approx(0.7));
+    REQUIRE(t(R1{-1})[0] == Catch::Approx(0.3));
+    REQUIRE(t(R1{-1})[1] == Catch::Approx(1.2));
+    REQUIRE(t(R1{0})[0] == Catch::Approx(1.5));
+    REQUIRE(t(R1{0})[1] == Catch::Approx(0.95));
+    REQUIRE(t(R1{1})[0] == Catch::Approx(2.7));
+    REQUIRE(t(R1{1})[1] == Catch::Approx(0.7));
 
     REQUIRE(t.velocityNorm() == Catch::Approx(l2Norm(0.5 * (b - a))));
 
@@ -80,22 +81,23 @@ TEST_CASE("LineTransformation", "[geometry]")
 
   SECTION("3D")
   {
+    using R1 = TinyVector<1>;
     using R3 = TinyVector<3>;
 
-    const R3 a = {0.3, 1.2, -1};
-    const R3 b = {2.7, 0.7, 0.3};
+    const R3 a{0.3, 1.2, -1};
+    const R3 b{2.7, 0.7, 0.3};
 
     const LineTransformation<3> t(a, b);
 
-    REQUIRE(t(-1)[0] == Catch::Approx(0.3));
-    REQUIRE(t(-1)[1] == Catch::Approx(1.2));
-    REQUIRE(t(-1)[2] == Catch::Approx(-1.));
-    REQUIRE(t(0)[0] == Catch::Approx(1.5));
-    REQUIRE(t(0)[1] == Catch::Approx(0.95));
-    REQUIRE(t(0)[2] == Catch::Approx(-0.35));
-    REQUIRE(t(1)[0] == Catch::Approx(2.7));
-    REQUIRE(t(1)[1] == Catch::Approx(0.7));
-    REQUIRE(t(1)[2] == Catch::Approx(0.3));
+    REQUIRE(t(R1{-1})[0] == Catch::Approx(0.3));
+    REQUIRE(t(R1{-1})[1] == Catch::Approx(1.2));
+    REQUIRE(t(R1{-1})[2] == Catch::Approx(-1.));
+    REQUIRE(t(R1{0})[0] == Catch::Approx(1.5));
+    REQUIRE(t(R1{0})[1] == Catch::Approx(0.95));
+    REQUIRE(t(R1{0})[2] == Catch::Approx(-0.35));
+    REQUIRE(t(R1{1})[0] == Catch::Approx(2.7));
+    REQUIRE(t(R1{1})[1] == Catch::Approx(0.7));
+    REQUIRE(t(R1{1})[2] == Catch::Approx(0.3));
 
     REQUIRE(t.velocityNorm() == Catch::Approx(l2Norm(0.5 * (b - a))));
 
diff --git a/tests/test_MathModule.cpp b/tests/test_MathModule.cpp
index 4d71927f13d3229d450c941250bb4202d071b109..6846b7c4e24a17dde8747afd9e4fcb1982678985 100644
--- a/tests/test_MathModule.cpp
+++ b/tests/test_MathModule.cpp
@@ -386,8 +386,8 @@ TEST_CASE("MathModule", "[language]")
   {
     SECTION("dot:R^1*R^1")
     {
-      TinyVector<1> arg0 = 3;
-      TinyVector<1> arg1 = 2;
+      TinyVector<1> arg0{3};
+      TinyVector<1> arg1{2};
 
       DataVariant arg0_variant = arg0;
       DataVariant arg1_variant = arg1;
diff --git a/tests/test_PrismTransformation.cpp b/tests/test_PrismTransformation.cpp
index f933fe308fbc64915faf04449aeca364f0904053..8ffa0516c653820357eaf3b5f1f7cea3c8b757ef 100644
--- a/tests/test_PrismTransformation.cpp
+++ b/tests/test_PrismTransformation.cpp
@@ -11,21 +11,21 @@ TEST_CASE("PrismTransformation", "[geometry]")
 {
   using R3 = TinyVector<3>;
 
-  const R3 a_hat = {0, 0, -1};
-  const R3 b_hat = {1, 0, -1};
-  const R3 c_hat = {0, 1, -1};
-  const R3 d_hat = {0, 0, +1};
-  const R3 e_hat = {1, 0, +1};
-  const R3 f_hat = {0, 1, +1};
-
-  const R3 m_hat = {1. / 3, 1. / 3, 0};
-
-  const R3 a = {1, 2, 0};
-  const R3 b = {3, 1, 3};
-  const R3 c = {2, 5, 2};
-  const R3 d = {0, 3, 1};
-  const R3 e = {1, 2, 5};
-  const R3 f = {3, 1, 7};
+  const R3 a_hat{0, 0, -1};
+  const R3 b_hat{1, 0, -1};
+  const R3 c_hat{0, 1, -1};
+  const R3 d_hat{0, 0, +1};
+  const R3 e_hat{1, 0, +1};
+  const R3 f_hat{0, 1, +1};
+
+  const R3 m_hat{1. / 3, 1. / 3, 0};
+
+  const R3 a{1, 2, 0};
+  const R3 b{3, 1, 3};
+  const R3 c{2, 5, 2};
+  const R3 d{0, 3, 1};
+  const R3 e{1, 2, 5};
+  const R3 f{3, 1, 7};
 
   const PrismTransformation t(a, b, c, d, e, f);
 
diff --git a/tests/test_PugsFunctionAdapter.cpp b/tests/test_PugsFunctionAdapter.cpp
index cb0dbc7c5e1b76dfeeb6390ddca70524e22d7dc1..d8708c49a9601d38602251a5179cad939382f291 100644
--- a/tests/test_PugsFunctionAdapter.cpp
+++ b/tests/test_PugsFunctionAdapter.cpp
@@ -221,7 +221,7 @@ let R33toR33zero: R^3x3 -> R^3x3, x -> 0;
         FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
         TinyVector<1> result = tests_adapter::TestBinary<TinyVector<1>(bool)>::one_arg(function_symbol_id, b);
 
-        REQUIRE(result == not b);
+        REQUIRE(result == TinyVector<1>{not b});
       }
 
       {
@@ -230,7 +230,7 @@ let R33toR33zero: R^3x3 -> R^3x3, x -> 0;
         FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
         TinyVector<1> result = tests_adapter::TestBinary<TinyVector<1>(bool)>::one_arg(function_symbol_id, b);
 
-        REQUIRE(result == not b);
+        REQUIRE(result == TinyVector<1>{not b});
       }
     }
 
@@ -245,7 +245,7 @@ let R33toR33zero: R^3x3 -> R^3x3, x -> 0;
         FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
         TinyMatrix<1> result = tests_adapter::TestBinary<TinyMatrix<1>(bool)>::one_arg(function_symbol_id, b);
 
-        REQUIRE(result == not b);
+        REQUIRE(result == TinyMatrix<1>{not b});
       }
 
       {
@@ -254,7 +254,7 @@ let R33toR33zero: R^3x3 -> R^3x3, x -> 0;
         FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
         TinyMatrix<1> result = tests_adapter::TestBinary<TinyMatrix<1>(bool)>::one_arg(function_symbol_id, b);
 
-        REQUIRE(result == not b);
+        REQUIRE(result == TinyMatrix<1>{not b});
       }
     }
 
@@ -268,7 +268,7 @@ let R33toR33zero: R^3x3 -> R^3x3, x -> 0;
       FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
       TinyVector<1> result = tests_adapter::TestBinary<TinyVector<1>(uint64_t)>::one_arg(function_symbol_id, n);
 
-      REQUIRE(result == n * n);
+      REQUIRE(result == TinyVector<1>{n * n});
     }
 
     {
@@ -281,7 +281,7 @@ let R33toR33zero: R^3x3 -> R^3x3, x -> 0;
       FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
       TinyMatrix<1> result = tests_adapter::TestBinary<TinyMatrix<1>(uint64_t)>::one_arg(function_symbol_id, n);
 
-      REQUIRE(result == n * n);
+      REQUIRE(result == TinyMatrix<1>{n * n});
     }
 
     {
@@ -294,7 +294,7 @@ let R33toR33zero: R^3x3 -> R^3x3, x -> 0;
       FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
       TinyVector<1> result = tests_adapter::TestBinary<TinyVector<1>(int64_t)>::one_arg(function_symbol_id, z);
 
-      REQUIRE(result == -z);
+      REQUIRE(result == TinyVector<1>{-z});
     }
 
     {
@@ -307,7 +307,7 @@ let R33toR33zero: R^3x3 -> R^3x3, x -> 0;
       FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
       TinyMatrix<1> result = tests_adapter::TestBinary<TinyMatrix<1>(int64_t)>::one_arg(function_symbol_id, z);
 
-      REQUIRE(result == -z);
+      REQUIRE(result == TinyMatrix<1>{-z});
     }
 
     {
@@ -320,7 +320,7 @@ let R33toR33zero: R^3x3 -> R^3x3, x -> 0;
       FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
       TinyVector<1> result = tests_adapter::TestBinary<TinyVector<1>(double)>::one_arg(function_symbol_id, x);
 
-      REQUIRE(result == x * x);
+      REQUIRE(result == TinyVector<1>{x * x});
     }
 
     {
@@ -333,7 +333,7 @@ let R33toR33zero: R^3x3 -> R^3x3, x -> 0;
       FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
       TinyMatrix<1> result = tests_adapter::TestBinary<TinyMatrix<1>(double)>::one_arg(function_symbol_id, x);
 
-      REQUIRE(result == x * x);
+      REQUIRE(result == TinyMatrix<1>{x * x});
     }
 
     {
diff --git a/tests/test_PyramidTransformation.cpp b/tests/test_PyramidTransformation.cpp
index 23b7369e625b20307ddab31d222b1b2ae3d607b4..8e1a41889036674e8afcd2e6c12be4cec4157825 100644
--- a/tests/test_PyramidTransformation.cpp
+++ b/tests/test_PyramidTransformation.cpp
@@ -13,19 +13,19 @@ TEST_CASE("PyramidTransformation", "[geometry]")
 {
   using R3 = TinyVector<3>;
 
-  const R3 a_hat = {-1, -1, +0};
-  const R3 b_hat = {+1, -1, +0};
-  const R3 c_hat = {+1, +1, +0};
-  const R3 d_hat = {-1, +1, +0};
-  const R3 e_hat = {+0, +0, +1};
-
-  const R3 m_hat = {0, 0, 1. / 5};
-
-  const R3 a = {1, 2, 0};
-  const R3 b = {3, 1, 3};
-  const R3 c = {2, 5, 2};
-  const R3 d = {0, 3, 1};
-  const R3 e = {1, 2, 5};
+  const R3 a_hat{-1, -1, +0};
+  const R3 b_hat{+1, -1, +0};
+  const R3 c_hat{+1, +1, +0};
+  const R3 d_hat{-1, +1, +0};
+  const R3 e_hat{+0, +0, +1};
+
+  const R3 m_hat{0, 0, 1. / 5};
+
+  const R3 a{1, 2, 0};
+  const R3 b{3, 1, 3};
+  const R3 c{2, 5, 2};
+  const R3 d{0, 3, 1};
+  const R3 e{1, 2, 5};
 
   const PyramidTransformation t(a, b, c, d, e);
 
diff --git a/tests/test_SmallArray.cpp b/tests/test_SmallArray.cpp
index 251134047754e6c2b8f5012c4d1713d973fa7823..2565ac0d7eb82beb828076840fda74d9b3ffbfa5 100644
--- a/tests/test_SmallArray.cpp
+++ b/tests/test_SmallArray.cpp
@@ -317,16 +317,16 @@ TEST_CASE("SmallArray", "[utils]")
     {
       using N2 = TinyVector<2, int>;
       SmallArray<N2> b(10);
-      b[0] = {13, 2};
-      b[1] = {1, 3};
-      b[2] = {8, -2};
-      b[3] = {-3, 2};
-      b[4] = {23, 4};
-      b[5] = {-1, -3};
-      b[6] = {13, 17};
-      b[7] = {0, 9};
-      b[8] = {12, 13};
-      b[9] = {9, -17};
+      b[0] = N2{13, 2};
+      b[1] = N2{1, 3};
+      b[2] = N2{8, -2};
+      b[3] = N2{-3, 2};
+      b[4] = N2{23, 4};
+      b[5] = N2{-1, -3};
+      b[6] = N2{13, 17};
+      b[7] = N2{0, 9};
+      b[8] = N2{12, 13};
+      b[9] = N2{9, -17};
 
       REQUIRE((sum(b) == N2{75, 28}));
     }
@@ -335,16 +335,16 @@ TEST_CASE("SmallArray", "[utils]")
     {
       using N22 = TinyMatrix<2, 2, int>;
       SmallArray<N22> b(10);
-      b[0] = {13, 2, 0, 1};
-      b[1] = {1, 3, 6, 3};
-      b[2] = {8, -2, -1, 21};
-      b[3] = {-3, 2, 5, 12};
-      b[4] = {23, 4, 7, 1};
-      b[5] = {-1, -3, 33, 11};
-      b[6] = {13, 17, 12, 13};
-      b[7] = {0, 9, 1, 14};
-      b[8] = {12, 13, -3, -71};
-      b[9] = {9, -17, 0, 16};
+      b[0] = N22{13, 2, 0, 1};
+      b[1] = N22{1, 3, 6, 3};
+      b[2] = N22{8, -2, -1, 21};
+      b[3] = N22{-3, 2, 5, 12};
+      b[4] = N22{23, 4, 7, 1};
+      b[5] = N22{-1, -3, 33, 11};
+      b[6] = N22{13, 17, 12, 13};
+      b[7] = N22{0, 9, 1, 14};
+      b[8] = N22{12, 13, -3, -71};
+      b[9] = N22{9, -17, 0, 16};
 
       REQUIRE((sum(b) == N22{75, 28, 60, 21}));
     }
diff --git a/tests/test_SquareTransformation.cpp b/tests/test_SquareTransformation.cpp
index a1ead02fc75d25f5ff0194b4157e22497d6005dc..5b3723e766886d8d6883900463a803ebbb94c216 100644
--- a/tests/test_SquareTransformation.cpp
+++ b/tests/test_SquareTransformation.cpp
@@ -15,17 +15,17 @@ TEST_CASE("SquareTransformation", "[geometry]")
   {
     using R2 = TinyVector<2>;
 
-    const R2 a_hat = {-1, -1};
-    const R2 b_hat = {+1, -1};
-    const R2 c_hat = {+1, +1};
-    const R2 d_hat = {-1, +1};
+    const R2 a_hat{-1, -1};
+    const R2 b_hat{+1, -1};
+    const R2 c_hat{+1, +1};
+    const R2 d_hat{-1, +1};
 
     const R2 m_hat = zero;
 
-    const R2 a = {0, 0};
-    const R2 b = {8, -2};
-    const R2 c = {12, 7};
-    const R2 d = {3, 7};
+    const R2 a{0, 0};
+    const R2 b{8, -2};
+    const R2 c{12, 7};
+    const R2 d{3, 7};
 
     const R2 m = 0.25 * (a + b + c + d);
 
@@ -113,9 +113,9 @@ TEST_CASE("SquareTransformation", "[geometry]")
     {
       using R2 = TinyVector<2>;
 
-      const R2 a = {1, 2};
-      const R2 b = {3, 1};
-      const R2 c = {2, 5};
+      const R2 a{1, 2};
+      const R2 b{3, 1};
+      const R2 c{2, 5};
 
       const SquareTransformation<2> t(a, b, c, c);
 
@@ -170,19 +170,19 @@ TEST_CASE("SquareTransformation", "[geometry]")
   {
     using R2 = TinyVector<2>;
 
-    const R2 a_hat = {-1, -1};
-    const R2 b_hat = {+1, -1};
-    const R2 c_hat = {+1, +1};
-    const R2 d_hat = {-1, +1};
+    const R2 a_hat{-1, -1};
+    const R2 b_hat{+1, -1};
+    const R2 c_hat{+1, +1};
+    const R2 d_hat{-1, +1};
 
     const R2 m_hat = zero;
 
     using R3 = TinyVector<3>;
 
-    const R3 a = {0, 0, -1};
-    const R3 b = {8, -2, 3};
-    const R3 c = {12, 7, 2};
-    const R3 d = {3, 7, 1};
+    const R3 a{0, 0, -1};
+    const R3 b{8, -2, 3};
+    const R3 c{12, 7, 2};
+    const R3 d{3, 7, 1};
 
     const R3 m = 0.25 * (a + b + c + d);
 
diff --git a/tests/test_TensorialGaussLegendreQuadrature.cpp b/tests/test_TensorialGaussLegendreQuadrature.cpp
index b2ac6c80c57009ea16ac8911e0609e59f60cb178..aa8708ec610e0c43fa9ad2ded1f50c78e5461b91 100644
--- a/tests/test_TensorialGaussLegendreQuadrature.cpp
+++ b/tests/test_TensorialGaussLegendreQuadrature.cpp
@@ -41,7 +41,7 @@ TEST_CASE("TensorialGaussLegendreQuadrature", "[analysis]")
       double alpha = 0.5 * (b - a);
       double beta  = 0.5 * (a + b);
 
-      auto x = [&alpha, &beta](auto x_hat) { return alpha * x_hat + beta; };
+      auto x = [&alpha, &beta](auto x_hat) { return TinyVector<1>{alpha * x_hat[0] + beta}; };
 
       auto value = weight_list[0] * f(x(point_list[0]));
       for (size_t i = 1; i < weight_list.size(); ++i) {
diff --git a/tests/test_TensorialGaussLobattoQuadrature.cpp b/tests/test_TensorialGaussLobattoQuadrature.cpp
index 1c08404ecaaf98fedb4dfa841a5ed8d845cfca09..e038b692ea6dce069e037640307e9d57864fffef 100644
--- a/tests/test_TensorialGaussLobattoQuadrature.cpp
+++ b/tests/test_TensorialGaussLobattoQuadrature.cpp
@@ -41,7 +41,7 @@ TEST_CASE("TensorialGaussLobattoQuadrature", "[analysis]")
       double alpha = 0.5 * (b - a);
       double beta  = 0.5 * (a + b);
 
-      auto x = [&alpha, &beta](auto x_hat) { return alpha * x_hat + beta; };
+      auto x = [&alpha, &beta](auto x_hat) { return TinyVector<1>{alpha * x_hat[0] + beta}; };
 
       auto value = weight_list[0] * f(x(point_list[0]));
       for (size_t i = 1; i < weight_list.size(); ++i) {
diff --git a/tests/test_TetrahedronTransformation.cpp b/tests/test_TetrahedronTransformation.cpp
index 91739b99e7030e1b8490869fb770f1978cce9334..504f48b34fabcae986f2e24000366066548a5c6c 100644
--- a/tests/test_TetrahedronTransformation.cpp
+++ b/tests/test_TetrahedronTransformation.cpp
@@ -11,32 +11,32 @@ TEST_CASE("TetrahedronTransformation", "[geometry]")
 {
   using R3 = TinyVector<3>;
 
-  const R3 a = {1, 2, 1};
-  const R3 b = {3, 1, 3};
-  const R3 c = {2, 5, 2};
-  const R3 d = {2, 3, 4};
+  const R3 a{1, 2, 1};
+  const R3 b{3, 1, 3};
+  const R3 c{2, 5, 2};
+  const R3 d{2, 3, 4};
 
   const TetrahedronTransformation t(a, b, c, d);
 
-  REQUIRE(t({0, 0, 0})[0] == Catch::Approx(1));
-  REQUIRE(t({0, 0, 0})[1] == Catch::Approx(2));
-  REQUIRE(t({0, 0, 0})[2] == Catch::Approx(1));
+  REQUIRE(t(R3{0, 0, 0})[0] == Catch::Approx(1));
+  REQUIRE(t(R3{0, 0, 0})[1] == Catch::Approx(2));
+  REQUIRE(t(R3{0, 0, 0})[2] == Catch::Approx(1));
 
-  REQUIRE(t({1, 0, 0})[0] == Catch::Approx(3));
-  REQUIRE(t({1, 0, 0})[1] == Catch::Approx(1));
-  REQUIRE(t({1, 0, 0})[2] == Catch::Approx(3));
+  REQUIRE(t(R3{1, 0, 0})[0] == Catch::Approx(3));
+  REQUIRE(t(R3{1, 0, 0})[1] == Catch::Approx(1));
+  REQUIRE(t(R3{1, 0, 0})[2] == Catch::Approx(3));
 
-  REQUIRE(t({0, 1, 0})[0] == Catch::Approx(2));
-  REQUIRE(t({0, 1, 0})[1] == Catch::Approx(5));
-  REQUIRE(t({0, 1, 0})[2] == Catch::Approx(2));
+  REQUIRE(t(R3{0, 1, 0})[0] == Catch::Approx(2));
+  REQUIRE(t(R3{0, 1, 0})[1] == Catch::Approx(5));
+  REQUIRE(t(R3{0, 1, 0})[2] == Catch::Approx(2));
 
-  REQUIRE(t({0, 0, 1})[0] == Catch::Approx(2));
-  REQUIRE(t({0, 0, 1})[1] == Catch::Approx(3));
-  REQUIRE(t({0, 0, 1})[2] == Catch::Approx(4));
+  REQUIRE(t(R3{0, 0, 1})[0] == Catch::Approx(2));
+  REQUIRE(t(R3{0, 0, 1})[1] == Catch::Approx(3));
+  REQUIRE(t(R3{0, 0, 1})[2] == Catch::Approx(4));
 
-  REQUIRE(t({0.25, 0.25, 0.25})[0] == Catch::Approx(2));
-  REQUIRE(t({0.25, 0.25, 0.25})[1] == Catch::Approx(11. / 4));
-  REQUIRE(t({0.25, 0.25, 0.25})[2] == Catch::Approx(2.5));
+  REQUIRE(t(R3{0.25, 0.25, 0.25})[0] == Catch::Approx(2));
+  REQUIRE(t(R3{0.25, 0.25, 0.25})[1] == Catch::Approx(11. / 4));
+  REQUIRE(t(R3{0.25, 0.25, 0.25})[2] == Catch::Approx(2.5));
 
   REQUIRE(t.jacobianDeterminant() == Catch::Approx(14));
 
diff --git a/tests/test_TriangleTransformation.cpp b/tests/test_TriangleTransformation.cpp
index e9561e859d727ea6c4019eaa78de125fbab8ba00..69803b7081a4ddbd529ce9288b2c14973832ed48 100644
--- a/tests/test_TriangleTransformation.cpp
+++ b/tests/test_TriangleTransformation.cpp
@@ -13,23 +13,23 @@ TEST_CASE("TriangleTransformation", "[geometry]")
   {
     using R2 = TinyVector<2>;
 
-    const R2 a = {1, 2};
-    const R2 b = {3, 1};
-    const R2 c = {2, 5};
+    const R2 a{1, 2};
+    const R2 b{3, 1};
+    const R2 c{2, 5};
 
     const TriangleTransformation<2> t(a, b, c);
 
-    REQUIRE(t({0, 0})[0] == Catch::Approx(1));
-    REQUIRE(t({0, 0})[1] == Catch::Approx(2));
+    REQUIRE(t(R2{0, 0})[0] == Catch::Approx(1));
+    REQUIRE(t(R2{0, 0})[1] == Catch::Approx(2));
 
-    REQUIRE(t({1, 0})[0] == Catch::Approx(3));
-    REQUIRE(t({1, 0})[1] == Catch::Approx(1));
+    REQUIRE(t(R2{1, 0})[0] == Catch::Approx(3));
+    REQUIRE(t(R2{1, 0})[1] == Catch::Approx(1));
 
-    REQUIRE(t({0, 1})[0] == Catch::Approx(2));
-    REQUIRE(t({0, 1})[1] == Catch::Approx(5));
+    REQUIRE(t(R2{0, 1})[0] == Catch::Approx(2));
+    REQUIRE(t(R2{0, 1})[1] == Catch::Approx(5));
 
-    REQUIRE(t({1. / 3, 1. / 3})[0] == Catch::Approx(2));
-    REQUIRE(t({1. / 3, 1. / 3})[1] == Catch::Approx(8. / 3));
+    REQUIRE(t(R2{1. / 3, 1. / 3})[0] == Catch::Approx(2));
+    REQUIRE(t(R2{1. / 3, 1. / 3})[1] == Catch::Approx(8. / 3));
 
     REQUIRE(t.jacobianDeterminant() == Catch::Approx(7));
 
@@ -56,29 +56,30 @@ TEST_CASE("TriangleTransformation", "[geometry]")
   {
     SECTION("Data")
     {
+      using R2 = TinyVector<2>;
       using R3 = TinyVector<3>;
 
-      const R3 a = {1, 2, 2};
-      const R3 b = {4, 1, 3};
-      const R3 c = {2, 5, 1};
+      const R3 a{1, 2, 2};
+      const R3 b{4, 1, 3};
+      const R3 c{2, 5, 1};
 
       const TriangleTransformation<3> t(a, b, c);
 
-      REQUIRE(t({0, 0})[0] == Catch::Approx(1));
-      REQUIRE(t({0, 0})[1] == Catch::Approx(2));
-      REQUIRE(t({0, 0})[2] == Catch::Approx(2));
+      REQUIRE(t(R2{0, 0})[0] == Catch::Approx(1));
+      REQUIRE(t(R2{0, 0})[1] == Catch::Approx(2));
+      REQUIRE(t(R2{0, 0})[2] == Catch::Approx(2));
 
-      REQUIRE(t({1, 0})[0] == Catch::Approx(4));
-      REQUIRE(t({1, 0})[1] == Catch::Approx(1));
-      REQUIRE(t({1, 0})[2] == Catch::Approx(3));
+      REQUIRE(t(R2{1, 0})[0] == Catch::Approx(4));
+      REQUIRE(t(R2{1, 0})[1] == Catch::Approx(1));
+      REQUIRE(t(R2{1, 0})[2] == Catch::Approx(3));
 
-      REQUIRE(t({0, 1})[0] == Catch::Approx(2));
-      REQUIRE(t({0, 1})[1] == Catch::Approx(5));
-      REQUIRE(t({0, 1})[2] == Catch::Approx(1));
+      REQUIRE(t(R2{0, 1})[0] == Catch::Approx(2));
+      REQUIRE(t(R2{0, 1})[1] == Catch::Approx(5));
+      REQUIRE(t(R2{0, 1})[2] == Catch::Approx(1));
 
-      REQUIRE(t({1. / 3, 1. / 3})[0] == Catch::Approx(7. / 3));
-      REQUIRE(t({1. / 3, 1. / 3})[1] == Catch::Approx(8. / 3));
-      REQUIRE(t({1. / 3, 1. / 3})[2] == Catch::Approx(2));
+      REQUIRE(t(R2{1. / 3, 1. / 3})[0] == Catch::Approx(7. / 3));
+      REQUIRE(t(R2{1. / 3, 1. / 3})[1] == Catch::Approx(8. / 3));
+      REQUIRE(t(R2{1. / 3, 1. / 3})[2] == Catch::Approx(2));
 
       REQUIRE(t.areaVariationNorm() == Catch::Approx(2 * std::sqrt(30)));
     }
@@ -87,9 +88,9 @@ TEST_CASE("TriangleTransformation", "[geometry]")
     {
       using R3 = TinyVector<3>;
 
-      const R3 a_hat = {0, 0, 0};
-      const R3 b_hat = {1, 0, 0};
-      const R3 c_hat = {0, 1, 0};
+      const R3 a_hat{0, 0, 0};
+      const R3 b_hat{1, 0, 0};
+      const R3 c_hat{0, 1, 0};
 
       const double theta = 2.3;
       TinyMatrix<3> r_theta(1, 0, 0,                               //
@@ -115,9 +116,9 @@ TEST_CASE("TriangleTransformation", "[geometry]")
     {
       using R3 = TinyVector<3>;
 
-      const R3 a = {1, 2, 2};
-      const R3 b = {4, 1, 3};
-      const R3 c = {2, 5, 1};
+      const R3 a{1, 2, 2};
+      const R3 b{4, 1, 3};
+      const R3 c{2, 5, 1};
 
       const TriangleTransformation<3> t(a, b, c);