Cast scalar types for compatibility

cpp/basix/maps.h

@@ -12,6 +12,27 @@
 namespace basix::maps
 {
 
+namespace impl
+{
+/// @private These structs are used to get the float/value type from a
+/// template argument, including support for complex types.
+template <typename T, typename = void>
+struct scalar_value_type
+{
+  /// @internal
+  typedef T value_type;
+};
+/// @private
+template <typename T>
+struct scalar_value_type<T, std::void_t<typename T::value_type>>
+{
+  typedef typename T::value_type value_type;
+};
+/// @private Convenience typedef
+template <typename T>
+using scalar_value_type_t = typename scalar_value_type<T>::value_type;
+} // namespace impl
+
 /// Map type
 enum class type
 {
@@ -40,14 +61,15 @@ void covariant_piola(O&& r, const P& U, const Q& /*J*/, double /*detJ*/,
                      const R& K)
 {
   using T = typename std::decay_t<O>::value_type;
+  using Z = typename impl::scalar_value_type_t<T>;
   for (std::size_t p = 0; p < U.extent(0); ++p)
   {
     // r_p = K^T U_p, where p indicates the p-th row
     for (std::size_t i = 0; i < r.extent(1); ++i)
     {
       T acc = 0;
       for (std::size_t k = 0; k < K.extent(0); ++k)
-        acc += K(k, i) * U(p, k);
+        acc += static_cast<Z>(K(k, i)) * U(p, k);
       r(p, i) = acc;
     }
   }
@@ -59,19 +81,20 @@ void contravariant_piola(O&& r, const P& U, const Q& J, double detJ,
                          const R& /*K*/)
 {
   using T = typename std::decay_t<O>::value_type;
+  using Z = typename impl::scalar_value_type_t<T>;
   for (std::size_t p = 0; p < U.extent(0); ++p)
   {
     for (std::size_t i = 0; i < r.extent(1); ++i)
     {
       T acc = 0;
       for (std::size_t k = 0; k < J.extent(1); ++k)
-        acc += J(i, k) * U(p, k);
+        acc += static_cast<Z>(J(i, k)) * U(p, k);
       r(p, i) = acc;
     }
   }
 
   std::transform(r.data_handle(), r.data_handle() + r.size(), r.data_handle(),
-                 [detJ](auto ri) { return ri / detJ; });
+                 [detJ](auto ri) { return ri / static_cast<Z>(detJ); });
 }
 
 /// Double covariant Piola map
@@ -81,6 +104,7 @@ void double_covariant_piola(O&& r, const P& U, const Q& J, double /*detJ*/,
 {
   namespace stdex = std::experimental;
   using T = typename std::decay_t<O>::value_type;
+  using Z = typename impl::scalar_value_type_t<T>;
   for (std::size_t p = 0; p < U.extent(0); ++p)
   {
     stdex::mdspan<const T, stdex::dextents<std::size_t, 2>> _U(
@@ -95,7 +119,7 @@ void double_covariant_piola(O&& r, const P& U, const Q& J, double /*detJ*/,
         T acc = 0;
         for (std::size_t k = 0; k < K.extent(0); ++k)
           for (std::size_t l = 0; l < _U.extent(1); ++l)
-            acc += K(k, i) * _U(k, l) * K(l, j);
+            acc += static_cast<Z>(K(k, i)) * _U(k, l) * static_cast<Z>(K(l, j));
         _r(i, j) = acc;
       }
     }
@@ -109,7 +133,7 @@ void double_contravariant_piola(O&& r, const P& U, const Q& J, double detJ,
 {
   namespace stdex = std::experimental;
   using T = typename std::decay_t<O>::value_type;
-
+  using Z = typename impl::scalar_value_type_t<T>;
   for (std::size_t p = 0; p < U.extent(0); ++p)
   {
     stdex::mdspan<const T, stdex::dextents<std::size_t, 2>> _U(
@@ -125,14 +149,14 @@ void double_contravariant_piola(O&& r, const P& U, const Q& J, double detJ,
         T acc = 0;
         for (std::size_t k = 0; k < J.extent(1); ++k)
           for (std::size_t l = 0; l < _U.extent(1); ++l)
-            acc += J(i, k) * _U(k, l) * J(j, l);
+            acc += static_cast<Z>(J(i, k)) * _U(k, l) * static_cast<Z>(J(j, l));
         _r(i, j) = acc;
       }
     }
   }
 
   std::transform(r.data_handle(), r.data_handle() + r.size(), r.data_handle(),
-                 [detJ](auto ri) { return ri / (detJ * detJ); });
+                 [detJ](auto ri) { return ri / static_cast<Z>(detJ * detJ); });
 }
 
 } // namespace basix::maps

cpp/basix/precompute.h

@@ -7,12 +7,34 @@
 #include "mdspan.hpp"
 #include <span>
 #include <tuple>
+#include <type_traits>
 #include <vector>
 
 /// Matrix and permutation precomputation
 namespace basix::precompute
 {
 
+namespace impl
+{
+/// @private These structs are used to get the float/value type from a
+/// template argument, including support for complex types.
+template <typename T, typename = void>
+struct scalar_value_type
+{
+  /// @internal
+  typedef T value_type;
+};
+/// @private
+template <typename T>
+struct scalar_value_type<T, std::void_t<typename T::value_type>>
+{
+  typedef typename T::value_type value_type;
+};
+/// @private Convenience typedef
+template <typename T>
+using scalar_value_type_t = typename scalar_value_type<T>::value_type;
+} // namespace impl
+
 /// Prepare a permutation
 ///
 /// This computes a representation of the permutation that allows the
@@ -229,6 +251,8 @@ void apply_matrix(const std::span<const std::size_t>& v_size_t,
                   const std::span<E>& data, std::size_t offset = 0,
                   std::size_t block_size = 1)
 {
+  using U = typename impl::scalar_value_type_t<E>;
+
   const std::size_t dim = v_size_t.size();
   apply_permutation(v_size_t, data, offset, block_size);
   for (std::size_t b = 0; b < block_size; ++b)
@@ -238,16 +262,18 @@ void apply_matrix(const std::span<const std::size_t>& v_size_t,
       for (std::size_t j = i + 1; j < dim; ++j)
       {
         data[block_size * (offset + i) + b]
-            += M(i, j) * data[block_size * (offset + j) + b];
+            += static_cast<U>(M(i, j)) * data[block_size * (offset + j) + b];
       }
     }
     for (std::size_t i = 1; i <= dim; ++i)
     {
-      data[block_size * (offset + dim - i) + b] *= M(dim - i, dim - i);
+      data[block_size * (offset + dim - i) + b]
+          *= static_cast<U>(M(dim - i, dim - i));
       for (std::size_t j = 0; j < dim - i; ++j)
       {
         data[block_size * (offset + dim - i) + b]
-            += M(dim - i, j) * data[block_size * (offset + j) + b];
+            += static_cast<U>(M(dim - i, j))
+               * data[block_size * (offset + j) + b];
       }
     }
   }
@@ -267,6 +293,8 @@ void apply_matrix_to_transpose(
     const std::span<E>& data, std::size_t offset = 0,
     std::size_t block_size = 1)
 {
+  using U = typename impl::scalar_value_type_t<E>;
+
   const std::size_t dim = v_size_t.size();
   const std::size_t data_size
       = (data.size() + (dim < block_size ? block_size - dim : 0)) / block_size;
@@ -278,16 +306,17 @@ void apply_matrix_to_transpose(
       for (std::size_t j = i + 1; j < dim; ++j)
       {
         data[data_size * b + offset + i]
-            += M(i, j) * data[data_size * b + offset + j];
+            += static_cast<U>(M(i, j)) * data[data_size * b + offset + j];
       }
     }
     for (std::size_t i = 1; i <= dim; ++i)
     {
-      data[data_size * b + offset + dim - i] *= M(dim - i, dim - i);
+      data[data_size * b + offset + dim - i]
+          *= static_cast<U>(M(dim - i, dim - i));
       for (std::size_t j = 0; j < dim - i; ++j)
       {
         data[data_size * b + offset + dim - i]
-            += M(dim - i, j) * data[data_size * b + offset + j];
+            += static_cast<U>(M(dim - i, j)) * data[data_size * b + offset + j];
       }
     }
   }