Merge branch 'main' into mesh-fields-tolerances

Author: jacobmerson

.github/workflows/cmake.yml

@@ -27,6 +27,13 @@ jobs:
         sudo apt-get install -yq cmake
         cmake --version
 
+    # There is a known bug using MPICH ubuntu-24.04 (ubuntu-latest as of 11/10/2025)
+    # https://bugs.launchpad.net/ubuntu/+source/mpich/+bug/2072338
+    - name: Install mpi
+      run: |
+        sudo apt-get update -yq
+        sudo apt-get install -yq openmpi-bin libopenmpi-dev
+
     ## Kokkos
     - name: Kokkos Checkout repo
       uses: actions/checkout@v4
@@ -51,7 +58,7 @@ jobs:
     - name: Omega_h Checkout repo
       uses: actions/checkout@v4
       with:
-        repository: sandialabs/omega_h
+        repository: SCOREC/omega_h
         path: omegah
 
     - name: Omega_h Configure CMake

.github/workflows/doxygen-gh-pages.yml

@@ -44,7 +44,7 @@ jobs:
       - name: Omega_h Checkout repo
         uses: actions/checkout@v4
         with:
-            repository: sandialabs/omega_h
+            repository: SCOREC/omega_h
             path: omegah
 
       - name: Omega_h Configure CMake

.github/workflows/pr_comment_trigger_self_hosted.yml

@@ -5,9 +5,6 @@ on:
   issue_comment:                                     
     types: [created]
 
-  # Allows you to run this workflow manually from the Actions tab
-  workflow_dispatch:
-
 jobs:
   build_and_test:
     permissions:
@@ -34,11 +31,11 @@ jobs:
         run: |
           set +e #avoid exiting when lua modules return non-zero on 'warning' messages
           source /etc/profile #provides module command
-          module use /opt/scorec/spack/rhel9/v0201_4/lmod/linux-rhel9-x86_64/Core/
-          module load gcc/12.3.0-iil3lno 
-          module load mpich/4.1.1-xpoyz4t
-          module load cmake/3.26.3-2duxfcd
-          module load cuda/12.1.1-zxa4msk
+          module use /opt/scorec/spack/rhel9/v0222_2/lmod/linux-rhel9-x86_64/Core/
+          module load gcc/13.2.0-4eahhas
+          module load mpich/4.2.3-62uy3hd
+          module load cuda/12.6.2-gqq65nw
+          module load cmake
           set -e
 
           echo "github.workspace ${{github.workspace}}"
@@ -93,15 +90,16 @@ jobs:
             -DCMAKE_CUDA_ARCHITECTURES="80" \
             -DKokkos_PREFIX=$kkbdir/install \
             -DBUILD_SHARED_LIBS=on
-          cmake --build $ohbdir -j 4 --target install
+          cmake --build $ohbdir -j 6 --target install
 
           #meshfields
           mfbdir=${workDir}/build-meshfields
           cmake -S meshfields_${{ github.event.issue.number }} -B $mfbdir \
             -DMeshFields_USE_Cabana=on \
             -DOmega_h_ROOT=$ohbdir/install \
             -DKokkos_ROOT=$kkbdir/install \
-            -DCabana_ROOT=$cabbdir/install
+            -DCabana_ROOT=$cabbdir/install \
+            -DMeshFields_IS_TESTING=on
           cmake --build $mfbdir -j4
           ctest --test-dir $mfbdir
 

CMakeLists.txt

@@ -9,9 +9,12 @@ 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.
-set_property(TARGET Omega_h::omega_h PROPERTY INTERFACE_COMPILE_OPTIONS "")
+find_package(MPI REQUIRED)
+
+# Verify Omega_h was built with Kokkos enabled
+if(NOT Omega_h_USE_Kokkos)
+  message(FATAL_ERROR "Omega_h must be built with Kokkos enabled (Omega_h_USE_Kokkos=ON)")
+endif()
 
 if(MeshFields_USE_Cabana)
   find_package(Cabana 0.7.0 REQUIRED)
@@ -32,6 +35,7 @@ set(MESHFIELD_HEADERS
   src/MeshField_Scan.hpp
   src/MeshField_Field.hpp
   src/MeshField.hpp
+  src/MeshField_SPR_ErrorEstimator.hpp
   "${CMAKE_CURRENT_BINARY_DIR}/MeshField_Config.hpp"
 )
 if(MeshFields_USE_Cabana)
@@ -112,8 +116,8 @@ endif()
 enable_testing()
 include(CTest)
 
-option(IS_TESTING "Build for CTest" OFF)
-message(STATUS "IS_TESTING: ${IS_TESTING}")
+option(MeshFields_IS_TESTING "Build for CTest" OFF)
+message(STATUS "MeshFields_IS_TESTING: ${MeshFields_IS_TESTING}")
 
 #check for valgrind
 find_program(VALGRIND_CMD valgrind DOC "Location of the valgrind program")
@@ -124,79 +128,125 @@ function(test_func_impl TEST_NAME)
   # need to run as a cmake script to capture assert and other 'system failures'
   # https://cmake.org/cmake/help/latest/prop_test/WILL_FAIL.html#prop_test:WILL_FAIL
   add_test(NAME ${TEST_NAME} COMMAND ${CMAKE_COMMAND} -E env ${TEST_STR})
+  if(TEST ${TEST_NAME})
+    set_property(TEST ${TEST_NAME} PROPERTY LABELS "meshfields::base")
+  endif()
 endfunction(test_func_impl)
 
-function(test_func TEST_NAME)
+#smoke tests that are always enabled
+function(smoke_test_func TEST_NAME)
   test_func_impl(${TEST_NAME} ${ARGN})
   if(TEST ${TEST_NAME})
-    set_property(TEST ${TEST_NAME} PROPERTY LABELS "base")
+    set_property(TEST ${TEST_NAME} PROPERTY LABELS "meshfields::smoke")
+  endif()
+endfunction(smoke_test_func)
+
+function(test_func TEST_NAME)
+  if(MeshFields_IS_TESTING)
+    test_func_impl(${TEST_NAME} ${ARGN})
   endif()
 endfunction(test_func)
 
+function(add_test_property TEST_NAME PROP)
+  if(TEST ${TEST_NAME})
+    set_property(TEST ${TEST_NAME} PROPERTY ${PROP} ${ARGN})
+  endif()
+endfunction(add_test_property)
+
 # Unlike test_func, will_fail_test_func assumes the command for the test will fail
 function(will_fail_test_func TEST_NAME)
-  test_func_impl(${TEST_NAME} ${ARGN})
-  set_property(TEST ${TEST_NAME} PROPERTY WILL_FAIL TRUE)
-  if(TEST ${TEST_NAME})
-    set_property(TEST ${TEST_NAME} PROPERTY LABELS "base")
+  if(MeshFields_IS_TESTING)
+    test_func_impl(${TEST_NAME} ${ARGN})
+    set_property(TEST ${TEST_NAME} PROPERTY WILL_FAIL TRUE)
   endif()
 endfunction()
 
 function(will_fail_valgrind_test_func TEST_NAME)
-  if(VALGRIND_CMD)
-    test_func_impl(${TEST_NAME} ${VALGRIND_CMD} ${ARGN})
-    set_property(TEST ${TEST_NAME} PROPERTY
-      FAIL_REGULAR_EXPRESSION "Invalid read;Invalid write"
-    )
-    set_property(TEST ${TEST_NAME} PROPERTY WILL_FAIL TRUE)
-    if(TEST ${TEST_NAME})
-      set_property(TEST ${TEST_NAME} PROPERTY LABELS "base")
+  if(MeshFields_IS_TESTING)
+    if(VALGRIND_CMD)
+      test_func_impl(${TEST_NAME} ${VALGRIND_CMD} ${ARGN})
+      set_property(TEST ${TEST_NAME} PROPERTY
+        FAIL_REGULAR_EXPRESSION "Invalid read;Invalid write"
+      )
+      set_property(TEST ${TEST_NAME} PROPERTY WILL_FAIL TRUE)
     endif()
   endif()
 endfunction()
 
-function(meshfields_add_exe EXE_NAME EXE_SRC)
+#smoke executables are always built
+function(meshfields_add_smoke_exe EXE_NAME EXE_SRC)
   add_executable(${EXE_NAME} ${EXE_SRC})
   target_link_libraries(${EXE_NAME} PRIVATE meshfields)
 endfunction()
 
-# Creating minimal reproduction of error
-meshfields_add_exe(KokkosTests test/testKokkos.cpp)
+function(meshfields_add_exe EXE_NAME EXE_SRC)
+  if(MeshFields_IS_TESTING)
+    add_executable(${EXE_NAME} ${EXE_SRC})
+    target_link_libraries(${EXE_NAME} PRIVATE meshfields)
+  endif()
+endfunction()
+
+meshfields_add_smoke_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(ElementJacobian3d test/testElementJacobian3d.cpp)
 meshfields_add_exe(CountIntegrator test/testCountIntegrator.cpp)
-meshfields_add_exe(OmegahTriTests test/testOmegahTri.cpp)
+meshfields_add_smoke_exe(OmegahTriTests test/testOmegahTri.cpp)
 meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
 meshfields_add_exe(PointMapping test/testPointMapping.cpp)
 meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp)
+meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp)
 
 if(MeshFields_USE_Cabana)
   meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp)
-  meshfields_add_exe(CabanaTests test/testCabana.cpp)
-  test_func(CabanaTests ./CabanaTests)
+  meshfields_add_smoke_exe(CabanaTests test/testCabana.cpp)
+  smoke_test_func(CabanaTests ./CabanaTests)
   test_func(ControllerPerformance ./ControllerPerformance)
 endif()
 
-test_func(KokkosTests ./KokkosTests)
+smoke_test_func(KokkosTests ./KokkosTests)
 test_func(SerializationTests ./SerializationTests)
 test_func(ElementTests ./ElementTests)
 test_func(ElementJacobian1d ./ElementJacobian1d)
 test_func(ElementJacobian2d ./ElementJacobian2d)
 test_func(ElementJacobian3d ./ElementJacobian3d)
 test_func(CountIntegrator ./CountIntegrator)
-test_func(OmegahTriTests ./OmegahTriTests)
+smoke_test_func(OmegahTriTests ./OmegahTriTests)
 test_func(PointMapping ./PointMapping)
-test_func(OmegahTetTest, ./OmegahTetTest)
+test_func(OmegahTetTest ./OmegahTetTest)
 if(MeshFields_USE_EXCEPTIONS)
   # exception caught - no error
   test_func(ExceptionTest ./ExceptionTest)
 else()
   will_fail_test_func(ExceptionTest ./ExceptionTest)
 endif()
 
+meshfields_add_exe(OmegahSprAdapt test/testSprThwaitesAdapt.cpp)
+test_func(OmegahSprAdapt ./OmegahSprAdapt
+  ${CMAKE_SOURCE_DIR}/meshes/thwaites_basal_effectiveStrain.osh/
+  thwaitesAdapted #output prefix
+  0.1 # adapt ratio
+  1.0 # min size
+  16.0 # max size
+)
+add_test_property(OmegahSprAdapt 
+  PASS_REGULAR_EXPRESSION
+  "afterAdapt nTri: 14439; solution_1 min -4084.78 max 213.322")
+
+test_func(OmegahSprAdapt_p2
+  ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 ${MPIEXEC_PREFLAGS}
+  ./OmegahSprAdapt ${CMAKE_SOURCE_DIR}/meshes/thwaites_basal_effectiveStrain.osh/
+  thwaitesAdapted_p2 #output prefix
+  0.1 # adapt ratio
+  1.0 # min size
+  16.0 # max size
+)
+add_test_property(OmegahSprAdapt_p2
+  PASS_REGULAR_EXPRESSION
+  "afterAdapt nTri: 14468; solution_1 min -4098.9 max 213.322")
+
 #Code Coverage set up -------------------------------------------------------
 
 option(meshfields_ENABLE_COVERAGE_BUILD "Do a coverage build" OFF)

docs/addingShapeFunctions.md

@@ -0,0 +1,18 @@
+In order to add shape function, you must create a struct in
+https://github.com/SCOREC/meshFields/blob/main/src/MeshField_Shape.hpp.
+
+You must define numNodes, meshEntDim, and Order as const variables as well as the functions
+getValues and getLocalGradients. You then need to define a mapping from how you store the
+nodes to the ordering of the shape function defined earlier. The Omegah mappings are in
+https://github.com/SCOREC/meshFields/blob/main/src/MeshField.hpp#L63. In the Jacobian
+tests we also an identity mapping defined as another example:
+https://github.com/SCOREC/meshFields/blob/main/test/testElementJacobian2d.cpp#L12.
+
+In order to use integration with the added shape function it must be added here:
+https://github.com/SCOREC/meshFields/blob/main/src/MeshField_Integrate.hpp#L44.
+
+Finally, when trying to compute an integral using the shape function, a class derived from the
+class Integrator (defined here:
+https://github.com/SCOREC/meshFields/blob/main/src/MeshField_Integrate.hpp#L204) must be
+declared and should implement the atPoints function. An example of this can be found here:
+https://github.com/SCOREC/meshFields/blob/main/test/testCountIntegrator.cpp#L41.

meshes/thwaites_basal_effectiveStrain.osh/0.osh

@@ -0,0 +1 @@
+1

meshes/thwaites_basal_effectiveStrain.osh/nparts

@@ -0,0 +1 @@
+10

meshes/thwaites_basal_effectiveStrain.osh/version

@@ -312,7 +312,8 @@ class OmegahMeshField {
     // HACK assumes there is a vertex field.. in the Field Mixin object
     auto field_view = field.vtxField.serialize();
     Omega_h::Write<FieldDataType> field_write(field_view);
-    mesh.add_tag(0, "field", 1, Omega_h::read(field_write));
+    mesh.add_tag(0, "field", 1, Omega_h::read(field_write), false,
+                 Omega_h::ArrayType::VectorND);
     Omega_h::vtk::write_parallel("foo.vtk", &mesh, mesh.dim());
   }
 

src/MeshField.hpp

@@ -170,14 +170,15 @@ struct FieldElement {
    */
   KOKKOS_INLINE_FUNCTION ValArray
   getValue(int ent, Kokkos::Array<Real, MeshEntDim + 1> localCoord) const {
-    assert(ent < numMeshEnts);
+    assert(ent >= 0);
+    assert(static_cast<size_t>(ent) < numMeshEnts);
     ValArray c;
     const auto shapeValues = shapeFn.getValues(localCoord);
-    for (int ci = 0; ci < NumComponents; ++ci)
+    for (size_t ci = 0; ci < NumComponents; ++ci)
       c[ci] = 0;
     for (auto topo : elm2dof.getTopology()) { // element topology
-      for (int ni = 0; ni < shapeFn.numNodes; ++ni) {
-        for (int ci = 0; ci < NumComponents; ++ci) {
+      for (size_t ni = 0; ni < shapeFn.numNodes; ++ni) {
+        for (size_t ci = 0; ci < NumComponents; ++ci) {
           auto map = elm2dof(ni, ci, ent, topo);
           const auto fval =
               field(map.entity, map.node, map.component, map.topo);
@@ -194,8 +195,8 @@ struct FieldElement {
   KOKKOS_INLINE_FUNCTION NodeArray getNodeValues(int ent) const {
     NodeArray c;
     for (auto topo : elm2dof.getTopology()) { // element topology
-      for (int ni = 0; ni < ShapeType::numNodes; ++ni) {
-        for (int d = 0; d < ShapeType::meshEntDim; ++d) {
+      for (size_t ni = 0; ni < ShapeType::numNodes; ++ni) {
+        for (size_t d = 0; d < ShapeType::meshEntDim; ++d) {
           auto map = elm2dof(ni, d, ent, topo);
           const auto fval =
               field(map.entity, map.node, map.component, map.topo);
@@ -219,11 +220,12 @@ struct FieldElement {
    * @return the result of evaluation
    */
   KOKKOS_INLINE_FUNCTION Real getJacobian1d(int ent) const {
-    assert(ent < numMeshEnts);
+    assert(ent >= 0);
+    assert(static_cast<size_t>(ent) < numMeshEnts);
     const auto nodalGradients = shapeFn.getLocalGradients();
     const auto nodeValues = getNodeValues(ent);
     auto g = nodalGradients[0] * nodeValues[0];
-    for (int i = 1; i < shapeFn.numNodes; ++i) {
+    for (size_t i = 1; i < shapeFn.numNodes; ++i) {
       g = g + nodalGradients[i] * nodeValues[i];
     }
     return g;
@@ -328,7 +330,7 @@ struct FieldElement {
           KOKKOS_LAMBDA(const int &ent, LO &lerrors) {
             Real sum = 0;
             LO isError = 0;
-            for (int i = 0; i < localCoords.extent(1); i++) {
+            for (size_t i = 0; i < localCoords.extent(1); i++) {
               if (localCoords(ent, i) < 0)
                 isError++;
               sum += localCoords(ent, i);
@@ -449,7 +451,7 @@ evaluate(FieldElement &fes, Kokkos::View<Real **> localCoords,
         KOKKOS_LAMBDA(const int &ent, LO &lerrors) {
           Real sum = 0;
           LO isError = 0;
-          for (int i = 0; i < localCoords.extent(1); i++) {
+          for (size_t i = 0; i < localCoords.extent(1); i++) {
             if (localCoords(ent, i) < 0)
               isError++;
             sum += localCoords(ent, i);
@@ -477,7 +479,7 @@ evaluate(FieldElement &fes, Kokkos::View<Real **> localCoords,
   LO numLocalCoords;
   Kokkos::deep_copy(numLocalCoords,
                     Kokkos::subview(offsets, offsets.size() - 1));
-  if (localCoords.extent(0) != numLocalCoords) {
+  if (localCoords.extent(0) != static_cast<size_t>(numLocalCoords)) {
     fail("The size of dimension 0 of the local coordinates input array (%zu) "
          "does not match the last entry of the offsets array (%d).\n",
          localCoords.extent(0), numLocalCoords);
@@ -491,10 +493,10 @@ evaluate(FieldElement &fes, Kokkos::View<Real **> localCoords,
         Kokkos::Array<Real, FieldElement::MeshEntDim + 1> lc;
         // TODO use nested parallel for?
         for (auto pt = offsets(ent); pt < offsets(ent + 1); pt++) {
-          for (int i = 0; i < localCoords.extent(1); i++) // better way?
+          for (size_t i = 0; i < localCoords.extent(1); i++) // better way?
             lc[i] = localCoords(pt, i);
           const auto val = fes.getValue(ent, lc);
-          for (int i = 0; i < numComponents; i++)
+          for (size_t i = 0; i < numComponents; i++)
             res(pt, i) = val[i];
         }
       });