HexRef class

.gitignore

@@ -4,3 +4,4 @@ Trilinos
 CMakeCache.txt
 CMakeCache
 CMakeFiles
+*~

CMakeLists.txt

@@ -12,8 +12,7 @@ set(CMAKE_CXX_STANDARD 17)
 # Build type
 if (NOT CMAKE_BUILD_TYPE)
   set(CMAKE_BUILD_TYPE RelWithDebInfo)
- # set(CMAKE_BUILD_TYPE Release)
-endif(NOT CMAKE_BUILD_TYPE)
+endif()
 
 # Set module path
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR}/cmake)
@@ -22,6 +21,28 @@ set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR}/cmake)
 ################################################################################
 # Dependencies
 ################################################################################
+# MATAR submodule
+# (Note: the way the MATAR submodule is handled may need to change in the
+# future. Currently, it is deployed as a header-only library, which is why it
+# seemed sufficient to simply include the MATAR header. However, the
+# requirements of linking Kokkos and other parts of Trilinos may change this.
+# Please update this as needed.)
+add_subdirectory(matar)
+include_directories(matar/src)
+
+# Optional Kokkos dependency for MATAR
+option(MATAR_WITH_KOKKOS "MATAR is to be compiled against Kokkos" OFF)
+if (MATAR_WITH_KOKKOS)
+  find_package(Kokkos REQUIRED)
+endif()
+
+# Optional Doxygen and Sphinx dependencies
+option(WITH_DOCS "Build the documentation locally" OFF)
+if (WITH_DOCS)
+  find_package(Doxygen REQUIRED)
+  find_package(Sphinx REQUIRED)
+endif()
+
 # Optional BLAS/LAPACK dependency (disabled by default)
 option(WITH_BLAS_LAPACK "Build ELEMENTS with BLAS/LAPACK" OFF)
 if (WITH_BLAS_LAPACK)
@@ -37,29 +58,16 @@ if (WITH_VTK)
   find_package(VTK REQUIRED NO_MODULE)
 endif()
 
-# Optional Doxygen and Sphinx dependencies
-option(WITH_DOCS "Build the documentation locally" OFF)
-if (WITH_DOCS)
-  find_package(Doxygen REQUIRED)
-  find_package(Sphinx REQUIRED)
-endif()
-
-# MATAR submodule
-# (Note: the way the MATAR submodule is handled may need to change in the
-# future. Currently, it is deployed as a header-only library, which is why it
-# seemed sufficient to simply include the MATAR header. However, the
-# requirements of linking Kokkos and other parts of Trilinos may change this.
-# Please update this as needed.)
-add_subdirectory(matar)
-include_directories(matar/src)
-
 
 ################################################################################
 # Build
 ################################################################################
 # Includes
 include_directories(common)
 include_directories(elements)
+if (MATAR_WITH_KOKKOS)
+  add_subdirectory(elements-kokkos)
+endif()
 include_directories(geometry)
 include_directories(slam)
 include_directories(swage)

build_scripts/build_matar_kokkos_openmp.sh

@@ -0,0 +1,47 @@
+#!/bin/sh
+
+# Build configuration (modify this!)
+ELEMENTS_DIR=/path/to/elements
+MATAR_DIR=/path/to/matar/within/elements
+PROJECT_BUILD_DIR=/path/to/project/build/directory
+export OMP_NUM_THREADS=sensible_number_of_openmp_threads_for_your_machine
+export OMP_PROC_BIND=spread
+
+
+# Build and install Kokkos with OpenMP
+KOKKOS_SOURCE_DIR=${MATAR_DIR}/src/Kokkos/kokkos
+KOKKOS_BUILD_DIR=${PROJECT_BUILD_DIR}/kokkos
+KOKKOS_INSTALL_DIR=${KOKKOS_BUILD_DIR}
+
+OPTIONS=(
+  -DCMAKE_BUILD_TYPE=Release
+  -DCMAKE_INSTALL_PREFIX=${KOKKOS_INSTALL_DIR}
+  -DCMAKE_CXX_STANDARD=17
+  -DKokkos_ENABLE_SERIAL=ON
+  -DKokkos_ENABLE_OPENMP=ON
+  -DKokkos_ENABLE_TESTS=OFF
+  -DBUILD_TESTING=OFF
+)
+
+cmake "${OPTIONS[@]}" -S ${KOKKOS_SOURCE_DIR} -B ${KOKKOS_BUILD_DIR}
+cmake --build ${KOKKOS_BUILD_DIR}
+cmake --install ${KOKKOS_BUILD_DIR}
+
+
+# Build and install ELEMENTS with MATAR/Kokkos
+ELEMENTS_SOURCE_DIR=${ELEMENTS_DIR}
+ELEMENTS_BUILD_DIR=${PROJECT_BUILD_DIR}/elements
+ELEMENTS_INSTALL_DIR=${ELEMENTS_BUILD_DIR}
+
+OPTIONS=(
+  -DCMAKE_BUILD_TYPE=Release
+  -DCMAKE_INSTALL_PREFIX=${ELEMENTS_INSTALL_DIR}
+  -DMATAR_WITH_KOKKOS=ON
+  -DKokkos_DIR=${KOKKOS_INSTALL_DIR}/lib/cmake/Kokkos
+  -DKOKKOS=ON
+  -DOPENMP=ON
+)
+
+cmake "${OPTIONS[@]}" -S ${ELEMENTS_SOURCE_DIR} -B ${ELEMENTS_BUILD_DIR}
+cmake --build ${ELEMENTS_BUILD_DIR} --verbose
+cmake --install ${ELEMENTS_BUILD_DIR}

common/error.h

@@ -19,3 +19,8 @@ class VtkGridError : public std::runtime_error {
  public:
   VtkGridError(const std::string info) : std::runtime_error(info) {}
 };
+
+class NotImplementedError : public std::runtime_error {
+ public:
+  NotImplementedError(const std::string info) : std::runtime_error(info) {}
+};

elements-kokkos/CMakeLists.txt

@@ -0,0 +1,3 @@
+set(ELEMENTS_KOKKOS_INCLUDES ${CMAKE_CURRENT_SOURCE_DIR}/include)
+
+add_subdirectory(src)

elements-kokkos/include/BernsteinProductBasis.h

@@ -0,0 +1,36 @@
+#pragma once
+
+
+#include "TensorProductBasis.h"
+
+
+template <typename NumType>
+struct BernsteinProductBasis : TensorProductBasis<NumType> { 
+  using TensorProductBasis<NumType>::Nd, TensorProductBasis<NumType>::Np, 
+      TensorProductBasis<NumType>::ijk, TensorProductBasis<NumType>::rad; 
+
+  SizeType N;
+  SizeType Ne;
+
+  // Work array for intermediate coefficients
+  NumType *C;
+
+  BernsteinProductBasis(const SizeType);
+  ~BernsteinProductBasis();
+
+  // Basis functions and basis function gradients
+  NumType eval_basis(const SizeType, const NumType *);
+  void eval_grad_basis(const SizeType, const NumType *, NumType *);
+
+  // Function approximation over element
+  NumType eval_approx(const NumType *, const NumType *); 
+  void eval_grad_approx(const NumType *, const NumType *, NumType *);
+
+  // Jacobian of spatial mapping
+  void eval_jac(const NumType *, const NumType *, const NumType *,
+                const NumType *, NumType *);
+  NumType eval_det_jac(const NumType *, const NumType *, const NumType *, 
+      const NumType *);
+  void eval_inv_jac(const NumType *, const NumType *, const NumType *,
+       const NumType *, NumType *);
+};

elements-kokkos/include/HexRef.h

@@ -0,0 +1,106 @@
+#pragma once
+
+
+// MATAR includes
+#include "MatarTypeDefs.h"
+
+
+enum BasisType {
+  LAGRANGE,
+  LEGENDRE,
+  BERNSTEIN
+};
+
+
+// Forward declarations
+template <typename NumType>
+struct TensorProductBasis;
+
+
+struct HexRef {
+  HexRef(int order_basis, BasisType type_basis);
+  ~HexRef();
+
+
+  // Cells (subelements used in volume integration)
+  int num_ref_cells_1D_;
+
+  int num_ref_cells_in_elem_;
+  int num_ref_cells_in_elem() const;
+
+  int cell_rid(int i, int j, int k) const;
+
+  MatarRealCArray ref_cell_positions_;
+  Real ref_cell_positions(int cell_rid, int dim) const;
+
+  MatarRealCArray ref_cell_g_weights_;
+  Real ref_cell_g_weights(int cell_rid) const;
+
+
+  // Patches (subelement faces used in surface integration)
+  int num_patches_in_elem_;
+  int num_sides_in_elem_;
+
+  const int patch_rlid_cell_neighbor_[6] = {
+    1, // side in neighbor cell mated with -xi side in this cell
+    0, // side in neighbor cell mated with +xi side in this cell
+    3, // side in neighbor cell mated with -eta side in this cell
+    2, // side in neighbor cell mated with +eta side in this cell
+    5, // side in neighbor cell mated with -zeta side in this cell
+    4  // side in neighbor cell mated with +zeta side in this cell
+  };
+  int patch_rlid_in_cell_neighbor(int patch_rlid) const;
+
+  MatarIntCArray ref_patches_in_cell_list_;
+  int ref_patches_in_cell(int cell_lid, int patch_rlid) const;
+
+  MatarRealCArray ref_patch_positions_;
+  Real ref_patch_positions(int patch_rid, int dim) const;
+
+  MatarRealCArray ref_patch_g_weights_;
+  Real ref_patch_g_weights(int patch_rid) const;
+
+  const Real cell_side_unit_normals_[18] = {
+    -1,  0,  0, // -xi 
+     1,  0,  0, // +xi 
+     0, -1,  0, // -eta 
+     0,  1,  0, // +eta 
+     0,  0, -1, // -zeta 
+     0,  0,  1  // +zeta 
+  };
+  Real cell_side_unit_normals(int side_rlid, int dim) const;
+
+
+  // Nodes (quadrature points) and vertices (degrees of freedom)
+  int num_ref_nodes_1D_;
+  int num_ref_nodes_in_elem_;
+  int num_ref_nodes() const;
+
+  int num_ref_verts_1D_;
+  int num_ref_verts_in_elem_;
+
+  Real *ref_verts_1D_;
+
+
+  // Evaluations of basis functions and their gradients at quadrature points
+  // in cells and on patches
+  static const int num_dim_ = 3;
+  int num_dim() const;
+
+  int num_basis_;
+  int num_basis() const;
+
+  MatarRealCArray ref_cell_basis_;
+  Real ref_cell_basis(int cell_rid, int basis_id) const;
+
+  MatarRealCArray ref_cell_gradient_;
+  Real ref_cell_gradient(int cell_rid, int basis_id, int dim) const;
+
+  MatarRealCArray ref_patch_basis_;
+  Real ref_patch_basis(int patch_rid, int basis_id) const;
+
+  MatarRealCArray ref_patch_gradient_;
+  Real ref_patch_gradient(int patch_rid, int basis_id, int dim) const;
+
+  TensorProductBasis<Real> *basis;
+};

elements-kokkos/include/LagrangeProductBasis.h

@@ -0,0 +1,38 @@
+#pragma once
+
+
+#include "TensorProductBasis.h"
+
+
+template <typename NumType>
+struct LagrangeProductBasis : TensorProductBasis<NumType> {
+  using TensorProductBasis<NumType>::Nd, TensorProductBasis<NumType>::Np, 
+      TensorProductBasis<NumType>::ijk, TensorProductBasis<NumType>::rad; 
+
+  SizeType N;
+  SizeType Ne;
+  const NumType *Z;
+  NumType *w;
+
+  // Work array for intermediate coefficients
+  NumType *C;
+
+  LagrangeProductBasis(const SizeType, const NumType *);
+  ~LagrangeProductBasis();
+
+  // Basis functions and basis function gradients
+  NumType eval_basis(const SizeType, const NumType *);
+  void eval_grad_basis(const SizeType, const NumType *, NumType *);
+
+  // Function approximation over element
+  NumType eval_approx(const NumType *, const NumType *);
+  void eval_grad_approx(const NumType *, const NumType *, NumType *);
+
+  // Jacobian of spatial mapping
+  void eval_jac(const NumType *, const NumType *, const NumType *, 
+      const NumType *, NumType *);
+  NumType eval_det_jac(const NumType *, const NumType *, const NumType *, 
+      const NumType *);
+  void eval_inv_jac(const NumType *, const NumType *, const NumType *, 
+      const NumType *, NumType *);
+};

elements-kokkos/include/LegendreProductBasis.h

@@ -0,0 +1,28 @@
+#pragma once
+
+
+#include "TensorProductBasis.h"
+
+
+template <typename NumType>
+struct LegendreProductBasis : TensorProductBasis<NumType> {
+  using TensorProductBasis<NumType>::Nd, TensorProductBasis<NumType>::Np, 
+      TensorProductBasis<NumType>::ijk, TensorProductBasis<NumType>::rad; 
+
+  SizeType N;
+  SizeType Ne;
+
+  // Work array for intermediate coefficients
+  NumType *C;
+
+  LegendreProductBasis(SizeType);
+  ~LegendreProductBasis();
+
+  // Basis functions and basis function gradients
+  NumType eval_basis(const SizeType, const NumType *);
+  void eval_grad_basis(const SizeType, const NumType *, NumType *);
+
+  // Function approximation over element
+  NumType eval_approx(const NumType *, const NumType *);
+  void eval_grad_approx(const NumType *, const NumType *, NumType *);
+};

elements-kokkos/include/MatarTypeDefs.h

@@ -0,0 +1,12 @@
+#include "common.h"
+#include "matar.h"
+
+#ifdef MATAR_WITH_KOKKOS
+typedef CArrayKokkos<Real> MatarRealCArray;
+typedef CArrayKokkos<SizeType> MatarUIntCArray;
+typedef CArrayKokkos<int> MatarIntCArray;
+#else
+typedef CArray<Real> MatarRealCArray;
+typedef CArray<SizeType> MatarUIntCArray;
+typedef CArray<int> MatarIntCArray;
+#endif // MATAR_WITH_KOKKOS