Merge pull request #53 from SCOREC/cws/integration

add support for linear and quadratic triangle integration

Author: cwsmith

CMakeLists.txt

@@ -7,6 +7,7 @@ include(CMakePackageConfigHelpers)
 
 option(MeshFields_USE_Cabana "Build with the Cabana storage backend" OFF)
 
+find_package(Kokkos REQUIRED)
 find_package(Omega_h REQUIRED)
 #Clear the omegah compilation flags that it passes to cuda.  Using the
 # kokkos target, and nvcc_wrapper, provide sufficient flags.
@@ -17,11 +18,14 @@ if(MeshFields_USE_Cabana)
 endif()
 
 set(MESHFIELD_HEADERS
+  src/MeshField_Defines.hpp
   src/MeshField_Utility.hpp
   src/MeshField_Macros.hpp
   src/KokkosController.hpp
   src/MeshField_Element.hpp
+  src/MeshField_Integrate.hpp
   src/MeshField_Shape.hpp
+  src/MeshField_ShapeField.hpp
   src/MeshField_Fail.hpp
   src/MeshField_For.hpp
   src/MeshField_Reduce.hpp
@@ -42,12 +46,13 @@ set(MESHFIELD_SOURCES
 
 add_library(meshfields ${MESHFIELD_SOURCES})
 
-target_compile_features(meshfields INTERFACE cxx_std_17)
+target_compile_features(meshfields INTERFACE cxx_std_20)
 target_include_directories(meshfields
   PUBLIC
   "$<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>" # for MeshField_Config.hpp
   "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/src>"
   "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>")
+target_link_libraries(meshfields INTERFACE Kokkos::kokkos)
 target_link_libraries(meshfields INTERFACE Omega_h::omega_h)
 if(Kokkos_ENABLE_CUDA)
   target_compile_options(meshfields INTERFACE "--expt-relaxed-constexpr")
@@ -159,8 +164,10 @@ endfunction()
 meshfields_add_exe(KokkosTests test/testKokkos.cpp)
 meshfields_add_exe(SerializationTests test/testSerialize.cpp)
 meshfields_add_exe(ElementTests test/testElement.cpp)
+meshfields_add_exe(ElementJacobian1d test/testElementJacobian1d.cpp)
+meshfields_add_exe(ElementJacobian2d test/testElementJacobian2d.cpp)
+meshfields_add_exe(CountIntegrator test/testCountIntegrator.cpp)
 meshfields_add_exe(OmegahElementTests test/testOmegahElement.cpp)
-meshfields_add_exe(OmegahCoordFieldTest test/testOmegahCoordField.cpp)
 meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
 
 if(MeshFields_USE_Cabana)
@@ -173,8 +180,10 @@ endif()
 test_func(KokkosTests ./KokkosTests)
 test_func(SerializationTests ./SerializationTests)
 test_func(ElementTests ./ElementTests)
+test_func(ElementJacobian1d ./ElementJacobian1d)
+test_func(ElementJacobian2d ./ElementJacobian2d)
+test_func(CountIntegrator ./CountIntegrator)
 test_func(OmegahElementTests ./OmegahElementTests)
-test_func(OmegahCoordFieldTest ./OmegahCoordFieldTest)
 if(MeshFields_USE_EXCEPTIONS)
   # exception caught - no error
   test_func(ExceptionTest ./ExceptionTest)

src/MeshField.hpp

@@ -36,7 +36,8 @@ template <typename ExecutionSpace, template <typename...> typename Controller =
                                        MeshField::KokkosController>
 decltype(MeshField::CreateCoordinateField<ExecutionSpace, Controller>(
     MeshField::MeshInfo()))
-createCoordinateField(MeshField::MeshInfo mesh_info, Omega_h::Reals coords) {
+createCoordinateField(const MeshField::MeshInfo &mesh_info,
+                      Omega_h::Reals coords) {
   const auto meshDim = mesh_info.dim;
   auto coordField =
       MeshField::CreateCoordinateField<ExecutionSpace, Controller>(mesh_info);
@@ -49,6 +50,11 @@ createCoordinateField(MeshField::MeshInfo mesh_info, Omega_h::Reals coords) {
   return coordField;
 }
 
+} // anonymous namespace
+
+namespace MeshField {
+
+namespace Omegah {
 struct LinearTriangleToVertexField {
   Omega_h::LOs triVerts;
   LinearTriangleToVertexField(Omega_h::Mesh &mesh)
@@ -159,9 +165,7 @@ template <int ShapeOrder> auto getTriangleElement(Omega_h::Mesh &mesh) {
   }
 }
 
-} // namespace
-
-namespace MeshField {
+} // namespace Omegah
 
 template <typename ExecutionSpace, template <typename...> typename Controller =
                                        MeshField::KokkosController>
@@ -210,9 +214,9 @@ class OmegahMeshField {
   }
 
   // evaluate a field at the specified local coordinate for each triangle
-  template <typename AnalyticFunction, typename ViewType, typename ShapeField>
+  template <typename ViewType, typename ShapeField>
   auto triangleLocalPointEval(ViewType localCoords, size_t NumPtsPerElem,
-                              AnalyticFunction func, ShapeField field) {
+                              ShapeField field) {
     const auto MeshDim = 2;
     if (mesh.dim() != MeshDim) {
       MeshField::fail("input mesh must be 2d\n");
@@ -222,7 +226,7 @@ class OmegahMeshField {
       MeshField::fail("input field order must be 1 or 2\n");
     }
 
-    const auto [shp, map] = getTriangleElement<ShapeOrder>(mesh);
+    const auto [shp, map] = Omegah::getTriangleElement<ShapeOrder>(mesh);
 
     MeshField::FieldElement f(meshInfo.numTri, field, shp, map);
     auto offsets = createOffsets(meshInfo.numTri, NumPtsPerElem);