support Cabana backend

src/CabanaController.hpp

@@ -94,8 +94,6 @@ class CabanaController {
 private:
   // all the type defenitions that are needed us to get the type of the slice
   // returned by the underlying AoSoA
-  using soa_t = Cabana::SoA<DataTypes, vecLen>;
-
   template <std::size_t index>
   using member_data_t =
       typename Cabana::MemberTypeAtIndex<index, DataTypes>::type;
@@ -112,6 +110,15 @@ class CabanaController {
   template <class T, int stride>
   using wrapper_slice_t = CabanaSliceWrapper<member_slice_t<T, stride>, T>;
 
+  template <size_t... I>
+  static auto construct(std::integer_sequence<size_t, I...>,
+                        std::vector<int> &obj) {
+    return std::tuple<
+        Cabana::AoSoA<Cabana::MemberTypes<Ts>, MemorySpace, vecLen>...>{
+        Cabana::AoSoA<Cabana::MemberTypes<Ts>, MemorySpace, vecLen>(
+            "sliceAoSoA", obj[I])...};
+  }
+
   template <typename T1, typename... Tx> void construct_sizes() {
     // There are no dynamic ranks w/ cabana controller.
     // So we only need to know extent.
@@ -134,26 +141,35 @@ class CabanaController {
       }
       break;
     }
-    extent_sizes[theta][0] = num_tuples;
+    extent_sizes[theta][0] = num_tuples[this->theta];
     this->theta += 1;
     if constexpr (sizeof...(Tx) != 0)
       construct_sizes<Tx...>();
   }
 
   // member vaiables
-  Cabana::AoSoA<DataTypes, MemorySpace, vecLen> aosoa;
-  const int num_tuples;
+  std::tuple<Cabana::AoSoA<Cabana::MemberTypes<Ts>, MemorySpace, vecLen>...>
+      aosoa;
+  // Cabana::AoSoA<DataTypes, MemorySpace, vecLen> aosoa;
+  int num_tuples[sizeof...(Ts)];
   unsigned short theta = 0;
   int extent_sizes[sizeof...(Ts)][MAX_RANK];
 
 public:
-  CabanaController() : num_tuples(0) {
+  CabanaController() {
     static_assert(sizeof...(Ts) != 0);
+    for (int i = 0; i < sizeof...(Ts); ++i) {
+      num_tuples[i] = 0;
+    }
     construct_sizes<Ts...>();
   }
-
-  CabanaController(int n) : aosoa("sliceAoSoA", n), num_tuples(n) {
+  CabanaController(const std::initializer_list<int> items) {
     static_assert(sizeof...(Ts) != 0);
+    std::vector<int> obj(items);
+    for (std::size_t i = 0; i < obj.size(); ++i) {
+      num_tuples[i] = obj[i];
+    }
+    aosoa = construct(std::make_integer_sequence<size_t, sizeof...(Ts)>{}, obj);
     construct_sizes<Ts...>();
   }
 
@@ -162,18 +178,19 @@ class CabanaController {
     assert(j >= 0);
     assert(j <= MAX_RANK);
     if (j == 0)
-      return this->tuples();
+      return this->tuples(i);
     else
       return extent_sizes[i][j];
   }
 
-  int tuples() const { return num_tuples; }
+  int tuples(int i) const { return num_tuples[i]; }
 
   template <std::size_t index> auto makeSlice() {
     // Creates wrapper object w/ meta data about the slice.
     using type = std::tuple_element_t<index, TypeTuple>;
-    const int stride = sizeof(soa_t) / sizeof(member_value_t<index>);
-    auto slice = Cabana::slice<index>(aosoa);
+    const int stride = sizeof(Cabana::SoA<Cabana::MemberTypes<type>, vecLen>) /
+                       sizeof(member_value_t<index>);
+    auto slice = Cabana::slice<0>(std::get<index>(aosoa));
     int sizes[MAX_RANK];
     for (int i = 0; i < MAX_RANK; i++)
       sizes[i] = this->size(index, i);

src/MeshField.hpp

@@ -41,8 +41,8 @@ createCoordinateField(MeshField::MeshInfo mesh_info, Omega_h::Reals coords) {
   auto coordField =
       MeshField::CreateCoordinateField<ExecutionSpace, Controller>(mesh_info);
   auto setCoordField = KOKKOS_LAMBDA(const int &i) {
-    coordField(0, 0, i, MeshField::Vertex) = coords[i * meshDim];
-    coordField(0, 1, i, MeshField::Vertex) = coords[i * meshDim + 1];
+    coordField(i, 0, 0, MeshField::Vertex) = coords[i * meshDim];
+    coordField(i, 0, 1, MeshField::Vertex) = coords[i * meshDim + 1];
   };
   MeshField::parallel_for(ExecutionSpace(), {0}, {mesh_info.numVtx},
                           setCoordField, "setCoordField");
@@ -180,6 +180,7 @@ class OmegahMeshField {
             getMeshInfo(mesh_in), mesh_in.coords())) {}
 
   template <typename DataType, size_t order, size_t dim>
+  // Ordering of field indexing changed to 'entity, node, component'
   auto CreateLagrangeField() {
     return MeshField::CreateLagrangeField<ExecutionSpace, Controller, DataType,
                                           order, dim>(meshInfo);

src/MeshField_Element.hpp

@@ -59,9 +59,20 @@ struct FieldElement {
                const ElementDofHolderAccessor elm2dofIn)
       : numMeshEnts(numMeshEntsIn), field(fieldIn), shapeFn(shapeFnIn),
         elm2dof(elm2dofIn) {}
-
-  using ValArray = Kokkos::Array<typename FieldAccessor::BaseType,
-                                 ShapeType::numComponentsPerDof>;
+  /* general template for baseType which simply sets type
+   */
+  template <typename T> struct baseType {
+    using type = T;
+  };
+  /* template specialization to recursively strip type to get base type
+   * Example: int[5][6] => int[6] => int
+   */
+  template <typename T, size_t N> struct baseType<T[N]> {
+    using type = typename baseType<T>::type;
+  };
+  using ValArray =
+      Kokkos::Array<typename baseType<typename FieldAccessor::BaseType>::type,
+                    ShapeType::numComponentsPerDof>;
   static const size_t NumComponents = ShapeType::numComponentsPerDof;
 
   /**
@@ -91,7 +102,7 @@ struct FieldElement {
         for (int ci = 0; ci < shapeFn.numComponentsPerDof; ++ci) {
           auto map = elm2dof(ni, ci, ent, topo);
           const auto fval =
-              field(map.node, map.component, map.entity, map.topo);
+              field(map.entity, map.node, map.component, map.topo);
           c[ci] += fval * shapeValues[ni];
         }
       }

src/MeshField_ShapeField.hpp

@@ -1,10 +1,10 @@
 #ifndef MESHFIELD_SHAPEFIELD_HPP
 #define MESHFIELD_SHAPEFIELD_HPP
 
-#include "KokkosController.hpp"
 #ifdef MESHFIELDS_ENABLE_CABANA
 #include "CabanaController.hpp"
 #endif
+#include "KokkosController.hpp"
 #include "MeshField_Field.hpp"
 #include "MeshField_Shape.hpp"
 #include <type_traits> //decltype
@@ -173,9 +173,29 @@
     if (meshInfo.numVtx <= 0) {
       fail("mesh has no vertices\n");
     }
-    using Ctrlr = Controller<MemorySpace, ExecutionSpace, DataType ***>;
+#ifdef MESHFIELDS_ENABLE_CABANA
+    using Ctrlr = std::conditional_t<
+        std::is_same_v<
+            Controller<ExecutionSpace, MemorySpace, DataType>,
+            MeshField::CabanaController<ExecutionSpace, MemorySpace, DataType>>,
+        Controller<ExecutionSpace, MemorySpace, DataType[1][1]>,
+        Controller<MemorySpace, ExecutionSpace, DataType ***>>;
     // 1 dof with 1 component per vtx
-    Ctrlr kk_ctrl({/*field 0*/ 1, 1, meshInfo.numVtx});
+    auto createController = [](const int numComp, auto numVtx) {
+      if constexpr (std::is_same_v<
+                        Controller<ExecutionSpace, MemorySpace, DataType>,
+                        MeshField::CabanaController<ExecutionSpace, MemorySpace,
+                                                    DataType>>) {
+        return Ctrlr({numVtx});
+      } else {
+        return Ctrlr({/*field 0*/ numVtx, 1, numComp});
+      }
+    };
+    Ctrlr kk_ctrl = createController(1, meshInfo.numVtx);
+#else
+    using Ctrlr = Controller<MemorySpace, ExecutionSpace, DataType ***>;
+    Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, 1});
+#endif
     auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
     using LA = LinearAccessor<decltype(vtxField)>;
     using LinearLagrangeShapeField = ShapeField<Ctrlr, LinearTriangleShape, LA>;
@@ -188,11 +208,32 @@
     if (meshInfo.numEdge <= 0) {
       fail("mesh has no edges\n");
     }
+#ifdef MESHFIELDS_ENABLE_CABANA
+    using Ctrlr = std::conditional_t<
+        std::is_same_v<
+            Controller<ExecutionSpace, MemorySpace, DataType>,
+            MeshField::CabanaController<ExecutionSpace, MemorySpace, DataType>>,
+        Controller<ExecutionSpace, MemorySpace, DataType[1][1], DataType[1][1]>,
+        Controller<MemorySpace, ExecutionSpace, DataType ***, DataType ***>>;
+    // 1 dof with 1 comp per vtx/edge
+    auto createController = [](const int numComp, auto numVtx, auto numEdge) {
+      if constexpr (std::is_same_v<
+                        Controller<ExecutionSpace, MemorySpace, DataType>,
+                        MeshField::CabanaController<ExecutionSpace, MemorySpace,
+                                                    DataType>>) {
+        return Ctrlr({numVtx, numEdge});
+      } else {
+        return Ctrlr({/*field 0*/ numVtx, 1, numComp,
+                      /*field 1*/ numEdge, 1, numComp});
+      }
+    };
+    Ctrlr kk_ctrl = createController(1, meshInfo.numVtx, meshInfo.numEdge);
+#else
     using Ctrlr =
         Controller<MemorySpace, ExecutionSpace, DataType ***, DataType ***>;
-    // 1 dof with 1 comp per vtx/edge
-    Ctrlr kk_ctrl({/*field 0*/ 1, 1, meshInfo.numVtx,
-                   /*field 1*/ 1, 1, meshInfo.numEdge});
+    Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, 1,
+                   /*field 1*/ meshInfo.numEdge, 1, 1});
+#endif
     auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
     auto edgeField = MeshField::makeField<Ctrlr, 1>(kk_ctrl);
     using QA = QuadraticAccessor<decltype(vtxField), decltype(edgeField)>;
@@ -224,6 +265,7 @@
  * @param meshInfo defines on-process mesh metadata
  * @return a linear ShapeField
  */
+
 template <typename ExecutionSpace, template <typename...> typename Controller =
                                        MeshField::KokkosController>
 auto CreateCoordinateField(const MeshInfo &meshInfo) {
@@ -232,9 +274,30 @@
   }
   using DataType = Real;
   using MemorySpace = typename ExecutionSpace::memory_space;
-  using Ctrlr = Controller<MemorySpace, ExecutionSpace, DataType ***>;
   const int numComp = meshInfo.dim;
-  Ctrlr kk_ctrl({/*field 0*/ 1, numComp, meshInfo.numVtx});
+// FIXME Oversized cabana datatypes when numComp = 1|2
+#ifdef MESHFIELDS_ENABLE_CABANA
+  using Ctrlr = std::conditional_t<
+      std::is_same_v<
+          Controller<ExecutionSpace, MemorySpace, DataType>,
+          MeshField::CabanaController<ExecutionSpace, MemorySpace, DataType>>,
+      Controller<ExecutionSpace, MemorySpace, DataType[1][3]>,
+      Controller<MemorySpace, ExecutionSpace, DataType ***>>;
+  auto createController = [](const int numComp, auto numVtx) {
+    if constexpr (std::is_same_v<
+                      Controller<ExecutionSpace, MemorySpace, DataType>,
+                      MeshField::CabanaController<ExecutionSpace, MemorySpace,
+                                                  DataType>>) {
+      return Ctrlr({numVtx});
+    } else {
+      return Ctrlr({/*field 0*/ numVtx, 1, numComp});
+    }
+  };
+  Ctrlr kk_ctrl = createController(numComp, meshInfo.numVtx);
+#else
+  using Ctrlr = Controller<MemorySpace, ExecutionSpace, DataType ***>;
+  Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, numComp});
+#endif
   auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
   using LA = LinearAccessor<decltype(vtxField)>;
   using LinearLagrangeShapeField = ShapeField<Ctrlr, LinearTriangleShape, LA>;