Multiresolution Based Grid Adaptation for Triangles and Cartesian Grids

.gitignore

@@ -4,6 +4,7 @@ build/
 
 .deps/
 .libs/
+.cache/
 
 libtool
 config.log

CMakeLists.txt

@@ -29,6 +29,8 @@ project(
     VERSION "${T8CODE_VERSION_MAJOR}.${T8CODE_VERSION_MINOR}.${T8CODE_VERSION_PATCH}"
     )
 
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
 include( GNUInstallDirs)
 include( CTest )
 include( CMakeDependentOption )
@@ -53,6 +55,7 @@ option( T8CODE_BUILD_FORTRAN_INTERFACE "Build t8code's Fortran interface" OFF )
 
 option( T8CODE_ENABLE_MPI "Enable t8code's features which rely on MPI" ON )
 option( T8CODE_ENABLE_VTK "Enable t8code's features which rely on VTK" OFF )
+option( T8CODE_ENABLE_MRA "Enable t8code's features which rely on MRA" OFF )
 option( T8CODE_ENABLE_OCC "Enable t8code's features which rely on OpenCASCADE" OFF )
 
 option( T8CODE_USE_SYSTEM_SC "Use system-installed sc library" OFF )
@@ -218,6 +221,13 @@ if( T8CODE_ENABLE_VTK )
     endif (VTK_FOUND)
 endif( T8CODE_ENABLE_VTK )
 
+if( T8CODE_ENABLE_MRA )
+  find_package(GSL REQUIRED)
+  if(GSL_FOUND)
+    message("Found GSL")
+  endif(GSL_FOUND)
+endif( T8CODE_ENABLE_MRA )
+
 if( T8CODE_ENABLE_OCC )
     find_package( OpenCASCADE REQUIRED COMPONENTS
     TKBO TKPrim TKTopAlgo

cmake/thirdparty.cmake

@@ -72,6 +72,22 @@ foreach(INDEX RANGE ${DEPS_RANGE})
         message(WARNING "Unknown thirdparty library type '${DEP_TYPE}' for ${DEP_NAME}")
     endif()
 
-    # 3. Populate the thirdparty library
-    FetchContent_MakeAvailable(${DEP_NAME})
+    # 3. Populate the thirdparty library.
+    # Build googletest statically even when t8code is a shared library: a shared
+    # libgtest.so collides at load time with a system libgtest.so of the same
+    # soname but different symbols (e.g. the DLR-SC MPI fork), causing runtime
+    # "undefined symbol" errors. A static gtest is baked into each test binary.
+    if(DEP_NAME STREQUAL "googletest")
+        set(_t8_saved_build_shared "${BUILD_SHARED_LIBS}")
+        set(_t8_saved_pic "${CMAKE_POSITION_INDEPENDENT_CODE}")
+        set(BUILD_SHARED_LIBS OFF)
+        set(CMAKE_POSITION_INDEPENDENT_CODE ON)
+        FetchContent_MakeAvailable(${DEP_NAME})
+        set(BUILD_SHARED_LIBS "${_t8_saved_build_shared}")
+        set(CMAKE_POSITION_INDEPENDENT_CODE "${_t8_saved_pic}")
+        unset(_t8_saved_build_shared)
+        unset(_t8_saved_pic)
+    else()
+        FetchContent_MakeAvailable(${DEP_NAME})
+    endif()
 endforeach()

cmake/thirdparty.json

@@ -29,6 +29,16 @@
                 "ref": "2296a990d8b6b54731a63be0ba5bc17b08cd1f3d",
                 "shallow": "TRUE"
             }
+        },
+        {
+            "name": "unordered_dense",
+            "depends_on_cmake_option": "",
+            "source": {
+                "type": "git",
+                "url": "https://github.com/martinus/unordered_dense.git",
+                "ref": "v4.5.0",
+                "shallow": "TRUE"
+            }
         }
     ]
 }

src/CMakeLists.txt

@@ -114,6 +114,109 @@ if( T8CODE_ENABLE_OCC )
     )
 endif()
 
+if( T8CODE_ENABLE_MRA )
+  target_compile_definitions( T8 PUBLIC T8_ENABLE_MRA )
+
+  # Add GSL library 
+  target_link_libraries( T8 PUBLIC GSL::gsl GSL::gslcblas )
+
+  # Add unordered_dense library https://github.com/martinus/unordered_dense
+  # Fetched via cmake/thirdparty.json; unordered_dense_SOURCE_DIR is set by FetchContent.
+  set( UNORDERED_DENSE_DIR "${unordered_dense_SOURCE_DIR}/include" )
+
+  target_include_directories(T8 PUBLIC
+    $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/src>
+    $<BUILD_INTERFACE:${UNORDERED_DENSE_DIR}>
+    $<INSTALL_INTERFACE:include>
+  )
+
+  target_sources( T8 PRIVATE
+    t8_mra/num/basis/dubiner.cxx
+    t8_mra/num/quadrature/dunavant.cxx
+  )
+
+  install( FILES
+    t8_mra/t8_mra.hxx
+    DESTINATION include/t8_mra )
+
+  install( FILES
+    t8_mra/core/multiscale.hxx
+    t8_mra/core/forest_backend.hxx
+    t8_mra/core/mst.hxx
+    t8_mra/core/shape_traits.hxx
+    t8_mra/core/face_neighbor.hxx
+    DESTINATION include/t8_mra/core )
+
+  install( FILES
+    t8_mra/core/shape/mst_policy.hxx
+    t8_mra/core/shape/cartesian.hxx
+    t8_mra/core/shape/triangle.hxx
+    DESTINATION include/t8_mra/core/shape )
+
+  install( FILES
+    t8_mra/core/adapt/grading.hxx
+    t8_mra/core/adapt/coarsen.hxx
+    t8_mra/core/adapt/refine.hxx
+    t8_mra/core/adapt/balance.hxx
+    DESTINATION include/t8_mra/core/adapt )
+
+  install( FILES
+    t8_mra/io/vtk.hxx
+    t8_mra/io/dataset.hxx
+    DESTINATION include/t8_mra/io )
+
+  install( FILES
+    t8_mra/dg/dg_base.hxx
+    t8_mra/dg/cartesian.hxx
+    t8_mra/dg/triangle.hxx
+    DESTINATION include/t8_mra/dg )
+
+  install( FILES
+    t8_mra/criteria/coarsening_criterion.hxx
+    t8_mra/criteria/refinement_criterion.hxx
+    DESTINATION include/t8_mra/criteria )
+
+  install( FILES
+    t8_mra/data/element_data.hxx
+    t8_mra/data/levelindex_map.hxx
+    t8_mra/data/levelindex_set.hxx
+    t8_mra/data/levelmultiindex.hxx
+    t8_mra/data/triangle_order.hxx
+    DESTINATION include/t8_mra/data )
+
+  install( FILES
+    t8_mra/data/shape/cartesian.hxx
+    t8_mra/data/shape/triangle.hxx
+    DESTINATION include/t8_mra/data/shape )
+
+  install( FILES
+    t8_mra/num/cell_geometry.hxx
+    t8_mra/num/dg_basis.hxx
+    t8_mra/num/geometry.hxx
+    t8_mra/num/mask_coefficients.hxx
+    t8_mra/num/mat.hxx
+    t8_mra/num/nodal_to_modal.hxx
+    DESTINATION include/t8_mra/num )
+
+  install( FILES
+    t8_mra/num/shape/cartesian.hxx
+    t8_mra/num/shape/triangle.hxx
+    DESTINATION include/t8_mra/num/shape )
+
+  install( FILES
+    t8_mra/num/basis/basis.hxx
+    t8_mra/num/basis/legendre.hxx
+    t8_mra/num/basis/dubiner.hxx
+    DESTINATION include/t8_mra/num/basis )
+
+  install( FILES
+    t8_mra/num/quadrature/quadrature.hxx
+    t8_mra/num/quadrature/gauss_legendre.hxx
+    t8_mra/num/quadrature/dunavant.hxx
+    t8_mra/num/quadrature/dunavant_table.hxx
+    DESTINATION include/t8_mra/num/quadrature )
+endif()
+
 if( T8CODE_BUILD_PEDANTIC )
   target_compile_options( T8 PUBLIC -pedantic )
   set (T8_CXXFLAGS "${T8_CXXFLAGS} -Wpedantic")

src/t8_mra/core/adapt/balance.hxx

@@ -0,0 +1,90 @@
+#pragma once
+
+#ifdef T8_ENABLE_MRA
+
+#include <cstddef>
+#include <vector>
+
+#include "sc_mpi.h"
+
+#include "t8.h"
+#include "t8_eclass/t8_eclass.h"
+#include "t8_element/t8_element.h"
+
+#include "t8_mra/core/adapt/grading.hxx"
+
+namespace t8_mra::adapt
+{
+
+/**
+ * @brief One balancing round: refine covering leaves across faces
+ *
+ * Every leaf resolves its face neighbours to their covering leaf; a covering
+ * leaf more than one level coarser is refined one level (children data =
+ * inverse two-scale with zero details, so the data is unchanged).
+ *
+ * @return Number of leaves marked in this round
+ */
+template <typename TMultiscale>
+[[nodiscard]] unsigned int
+balance_round (TMultiscale &mra)
+{
+  clear_state (mra);
+
+  int mpirank;
+  int mpisize;
+  sc_MPI_Comm_rank (mra.grid.comm, &mpirank);
+  sc_MPI_Comm_size (mra.grid.comm, &mpisize);
+  std::vector<std::vector<size_t>> outgoing (mpisize);
+
+  mra.grid.for_each_face_neigh (
+    [] (const auto &lmi) { return lmi.level () >= 2; },
+    [&] (const auto &, t8_eclass_t tree_class, t8_gloidx_t neigh_gtreeid, t8_element_t *neigh_element,
+         const auto &neigh_lmi) {
+      if (refine_covering_leaf (mra, neigh_lmi, 0, 1u, no_prior_marks {}) < 0 && mpisize > 1) {
+        const auto owner = mra.grid.find_owner (neigh_gtreeid, neigh_element, tree_class);
+
+        if (owner != mpirank)
+          outgoing[owner].push_back (neigh_lmi.index);
+      }
+    });
+
+  if (mpisize > 1)
+    exchange_refine_requests (mra, outgoing, 0, 1u, no_prior_marks {});
+
+  const auto num_marked = mra.grid.global_num_marks (num_refinement_marks (mra, 0, mra.grid.maximum_level));
+
+  if (num_marked == 0) {
+    clear_state (mra);
+    return 0;
+  }
+
+  apply_refinement (mra, 0, mra.grid.maximum_level, 0);
+  clear_state (mra);
+
+  return num_marked;
+}
+
+/**
+ * @brief Restore the 2:1 face balance of the grid
+ *
+ * Rounds iterate until no leaf has a face neighbour more than one level coarser;
+ * a jump of k levels resolves in k-1 rounds. Terminates: every round refines at
+ * least one leaf and levels are bounded by max_level.
+ */
+template <typename TMultiscale>
+void
+balance (TMultiscale &mra)
+{
+  auto rounds = 0;
+
+  while (balance_round (mra) > 0)
+    t8_debugf ("MRA balance round %d\n", rounds++);
+
+  if (rounds > 0)
+    mra.grid.repartition ();
+}
+
+}  // namespace t8_mra::adapt
+
+#endif  // T8_ENABLE_MRA