Feature: mesh deformation with RBFs (TPS and CPC2)

CMakeLists.txt

@@ -54,6 +54,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_OCC "Enable t8code's features which rely on OpenCASCADE" OFF )
+option( T8CODE_ENABLE_EIGEN "Enable t8code's features which rely on eigen" OFF )
 
 option( T8CODE_USE_SYSTEM_SC "Use system-installed sc library" OFF )
 option( T8CODE_USE_SYSTEM_P4EST "Use system-installed p4est library" OFF )

cmake/thirdparty.cmake

@@ -40,13 +40,14 @@ foreach(INDEX RANGE ${DEPS_RANGE})
 
     # If the DEP_CMAKE_OPTION field is non-empty, check the CMake option.
     if(NOT DEP_CMAKE_OPTION STREQUAL "")
+        string(JSON DEP_CMAKE_OPTION_VALUE GET "${DEPS_JSON}" "thirdparty" ${INDEX} "cmake_option_install_value")
         # If the named option variable does not exist in the CMake cache, abort.
         if(NOT DEFINED ${DEP_CMAKE_OPTION})
             message(FATAL_ERROR "Loading thirdparty library ${DEP_NAME} at index ${INDEX} references unknown CMake option '${DEP_CMAKE_OPTION}'. Aborting.")
         else()
             # If the named option is defined but set to ON, skip this thirdparty library.
-            if(${${DEP_CMAKE_OPTION}})
-                message(STATUS "Skipping FetchContent-step for thirdparty library ${DEP_NAME} because CMake option '${DEP_CMAKE_OPTION}' is '${${DEP_CMAKE_OPTION}}'")
+            if(NOT ${DEP_CMAKE_OPTION} STREQUAL ${DEP_CMAKE_OPTION_VALUE})
+                message(STATUS "Skipping FetchContent-step for thirdparty library ${DEP_NAME} because CMake option '${DEP_CMAKE_OPTION}' is set to '${${DEP_CMAKE_OPTION}}'")
                 continue()
             endif()
         endif()

cmake/thirdparty.json

@@ -13,6 +13,7 @@
         {
             "name": "SC",
             "depends_on_cmake_option": "T8CODE_USE_SYSTEM_SC",
+            "cmake_option_install_value": "OFF",
             "source": {
                 "type": "git",
                 "url": "https://github.com/cburstedde/libsc.git",
@@ -23,12 +24,24 @@
         {
             "name": "P4EST",
             "depends_on_cmake_option": "T8CODE_USE_SYSTEM_P4EST",
+            "cmake_option_install_value": "OFF",
             "source": {
                 "type": "git",
                 "url": "https://github.com/cburstedde/p4est.git",
                 "ref": "2296a990d8b6b54731a63be0ba5bc17b08cd1f3d",
                 "shallow": "TRUE"
             }
+        },
+        {
+            "name": "EIGEN",
+            "depends_on_cmake_option": "T8CODE_ENABLE_EIGEN",
+            "cmake_option_install_value": "ON",
+            "source": {
+                "type": "git",
+                "url": "https://gitlab.com/libeigen/eigen.git",
+                "ref": "5.0.1",
+                "shallow": "TRUE"
+            }
         }
     ]
-}
+}

example/CMakeLists.txt

@@ -71,7 +71,6 @@ add_t8_example( NAME t8_cmesh_set_join_by_vertices      SOURCES cmesh/t8_cmesh_s
 add_t8_example( NAME t8_cmesh_geometry_examples         SOURCES cmesh/t8_cmesh_geometry_examples.cxx )
 add_t8_example( NAME t8_cmesh_create_partitioned        SOURCES cmesh/t8_cmesh_create_partitioned.cxx )
 add_t8_example( NAME t8_cmesh_hypercube_pad             SOURCES cmesh/t8_cmesh_hypercube_pad.cxx )
-add_t8_example( NAME t8_cmesh_mesh_deformation          SOURCES cmesh/t8_cmesh_mesh_deformation.cxx )
 
 add_t8_example( NAME t8_test_ghost                      SOURCES forest/t8_test_ghost.cxx )
 add_t8_example( NAME t8_test_face_iterate               SOURCES forest/t8_test_face_iterate.cxx )
@@ -89,6 +88,14 @@ if(T8CODE_ENABLE_VTK)
     add_t8_example( NAME t8_cmesh_read_from_vtk             SOURCES IO/cmesh/vtk/t8_cmesh_read_from_vtk.cxx )
 endif()
 
+if (T8CODE_ENABLE_EIGEN)
+    add_t8_example( NAME t8_rbf_test                        SOURCES cmesh/t8_rbf.cxx )
+    target_link_libraries( t8_rbf_test                      PRIVATE Eigen3::Eigen )
+    add_t8_example( NAME t8_cmesh_mesh_deformation          SOURCES cmesh/t8_cmesh_mesh_deformation.cxx )
+    target_link_libraries( t8_cmesh_mesh_deformation        PRIVATE Eigen3::Eigen )
+
+endif ()
+
 add_t8_example( NAME t8_gmsh_to_vtk                     SOURCES IO/forest/gmsh/t8_gmsh_to_vtk.cxx )
 
 add_t8_example( NAME t8_example_spheres                 SOURCES remove/t8_example_spheres.cxx )

example/cmesh/t8_cmesh_mesh_deformation.cxx

@@ -28,22 +28,43 @@
 #include <t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_connectivity.hxx>
 #include <t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.h>
 #include <t8_schemes/t8_default/t8_default.hxx>
-#if T8CODE_ENABLE_OCC
+#if T8_ENABLE_OCC && T8_ENABLE_EIGEN
 #include <t8_cad/t8_cad_handle.hxx>
 #include <t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.hxx>
-#endif /* T8CODE_ENABLE_OCC */
+#endif /* T8_ENABLE_OCC and T8_ENABLE_EIGEN*/
 #include <t8_vtk/t8_vtk_writer.h>
 #include <sc_options.h>
 
 #include <iostream>
 #include <string>
 #include <vector>
 #include <array>
+#include <filesystem>
+#include <algorithm>
+
+#if T8_ENABLE_OCC && T8_ENABLE_EIGEN
+namespace fs = std::filesystem;
+
+static std::vector<fs::path>
+findBrepFiles (const char *folder)
+{
+  std::vector<fs::path> files;
+
+  for (const auto &entry : fs::directory_iterator (folder)) {
+    if (entry.is_regular_file () && entry.path ().extension () == ".brep") {
+      files.push_back (entry.path ());
+    }
+  }
+
+  return files;
+}
+
+#endif /* T8_ENABLE_OCC && T8_ENABLE_EIGEN */
 
 int
 main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv)
 {
-#if T8CODE_ENABLE_OCC
+#if T8_ENABLE_OCC && T8_ENABLE_EIGEN
 
   char usage[BUFSIZ];
   /* Brief help message. */
@@ -76,7 +97,7 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv)
   SC_CHECK_MPI (mpiret);
 
   /* Initialize the sc library, has to happen before we initialize t8code. */
-  sc_init (sc_MPI_COMM_WORLD, 1, 1, NULL, SC_LP_PRODUCTION);
+  sc_init (sc_MPI_COMM_WORLD, 1, 1, NULL, SC_LP_ESSENTIAL);
 
   /* Initialize t8code with log level SC_LP_PRODUCTION. See sc.h for more info on the log levels. */
   t8_init (SC_LP_PRODUCTION);
@@ -85,15 +106,20 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv)
   const char *msh_file = NULL;
   const char *brep_file = NULL;
   int dim, level;
+  int rbf_type_int = 0;
+  double scale_factor_support_radius = 1.5;
 
   /* Initialize command line argument parser. */
   sc_options_t *opt = sc_options_new (argv[0]);
   sc_options_add_switch (opt, 'h', "help", &helpme, "Display a short help message.");
   sc_options_add_string (opt, 'm', "mshfile", &msh_file, NULL, "File prefix of the input mesh file (without .msh)");
   sc_options_add_string (opt, 'b', "brepfile", &brep_file, NULL,
-                         "File prefix of the deformation geometry file (without .brep)");
+                         "Path to the folder containing the deformation geometry files (.brep)");
   sc_options_add_int (opt, 'd', "dimension", &dim, 0, "Dimension of the mesh (1, 2 or 3)");
   sc_options_add_int (opt, 'l', "level", &level, 2, "Uniform refinement level for the input mesh. Default: 2");
+  sc_options_add_int (opt, 't', "rbftype", &rbf_type_int, 0, "RBF type (0 for CP_C2, 1 for TPS). Default: 0");
+  sc_options_add_double (opt, 's', "scalefactor", &scale_factor_support_radius, 1.5,
+                         "Scale factor for the support radius. Default: 1.5");
 
   int parsed = sc_options_parse (t8_get_package_id (), SC_LP_ERROR, opt, argc, argv);
 
@@ -106,7 +132,7 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv)
     sc_options_print_usage (t8_get_package_id (), SC_LP_ERROR, opt, NULL);
   }
   else if (dim < 1 || dim > 3) {
-    t8_global_errorf ("ERROR: Invalid mesh dimension: dim=%d. Dimension must be 1, 2 or 3.\n\n", dim);
+    t8_global_errorf ("ERROR: Invalid mesh dimension: dim=%d. Dimension must be 1, 2 or 3.\n", dim);
     sc_options_print_usage (t8_get_package_id (), SC_LP_ERROR, opt, NULL);
   }
   else {
@@ -120,28 +146,35 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv)
     t8_cmesh_from_msh_file (&cmesh, msh_file, 0, comm, dim, 0, 1);
     t8_forest_t forest = t8_forest_new_uniform (cmesh, t8_scheme_new_default (), level, 0, comm);
 
-    /* Load CAD geometry from .brep file. */
-    auto cad = std::make_shared<t8_cad_handle> (brep_file);
-
     /* Initialize the deformation object for the given mesh. */
     t8_cmesh_mesh_deformation deformation (cmesh);
 
-    /* Calculate displacements. */
-    auto displacements = deformation.calculate_displacement_surface_vertices (cad.get ());
+    /** Save the input RBF type. */
+    t8_rbf_function_type rbf_type = static_cast<t8_rbf_function_type> (rbf_type_int);
 
     /* Write output. */
-    t8_forest_vtk_write_file (forest, "input_forest", 1, 1, 1, 1, 0, 0, NULL);
+    t8_forest_vtk_write_file (forest, "deformed_forest_step_0", 1, 1, 1, 1, 0, 0, NULL);
 
-    /* Apply displacements. */
-    deformation.apply_vertex_displacements (displacements, cad);
+    auto brep_files = findBrepFiles (brep_file);
+    std::sort (brep_files.begin (), brep_files.end ());
 
-    /* Write output. */
-    t8_forest_vtk_write_file (forest, "deformed_forest", 1, 1, 1, 1, 0, 0, NULL);
+    int ifile = 0;
+    for (const auto &file : brep_files) {
+      auto file_without_ext = file.parent_path () / file.stem ();
+      auto cad_deformed = std::make_shared<t8_cad_handle> (file_without_ext.c_str ());
+
+      auto displacements = deformation.calculate_displacement_surface_vertices (cad_deformed.get (), rbf_type,
+                                                                                scale_factor_support_radius);
+
+      deformation.apply_vertex_displacements (displacements, cad_deformed, rbf_type);
 
+      std::string output_name = "deformed_forest_step_" + std::to_string (ifile++);
+      t8_forest_vtk_write_file (forest, output_name.c_str (), 1, 1, 1, 1, 0, 0, NULL);
+    }
     /* Cleanup. */
     t8_forest_unref (&forest);
 
-    t8_global_productionf ("Mesh deformation completed.");
+    t8_global_productionf ("Mesh deformation completed.\n");
   }
 
   sc_options_destroy (opt);
@@ -150,9 +183,9 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv)
   mpiret = sc_MPI_Finalize ();
   SC_CHECK_MPI (mpiret);
 
-#else  /* T8CODE_ENABLE_OCC */
-  t8_global_errorf ("ERROR: This example requires OpenCASCADE support to be enabled in t8code.\n");
-#endif /* T8CODE_ENABLE_OCC */
+#else  /* T8_ENABLE_OCC and T8_ENABLE_EIGEN*/
+  t8_global_errorf ("ERROR: This example requires OpenCASCADE and Eigen support to be enabled in t8code.\n");
+#endif /* T8_ENABLE_OCC and T8_ENABLE_EIGEN*/
 
   return 0;
 }

example/cmesh/t8_rbf.cxx

@@ -0,0 +1,85 @@
+/*
+  This file is part of t8code.
+  t8code is a C library to manage a collection (a forest) of multiple
+  connected adaptive space-trees of general element classes in parallel.
+
+  Copyright (C) 2026 the developers
+
+  t8code is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2 of the License, or
+  (at your option) any later version.
+
+  t8code is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with t8code; if not, write to the Free Software Foundation, Inc.,
+  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+*/
+
+/** \file t8_rbf.cxx
+ *  This file implements an example for CAD-based mesh deformation.
+ */
+#include <t8.h>
+#if T8_ENABLE_EIGEN
+#include <Eigen/Dense>
+#endif
+
+int
+main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv)
+{
+#if T8_ENABLE_EIGEN
+  t8_global_productionf ("--- STARTING EIGEN BASIC TEST --- \n");
+  /*
+   * Initialization.
+   */
+
+  /* Initialize MPI. This has to happen before we initialize sc or t8code. */
+  int mpiret = sc_MPI_Init (&argc, &argv);
+
+  /* Error check the MPI return value. */
+  SC_CHECK_MPI (mpiret);
+
+  /* Initialize the sc library, has to happen before we initialize t8code. */
+  sc_init (sc_MPI_COMM_WORLD, 1, 1, NULL, SC_LP_PRODUCTION);
+
+  /* Initialize t8code with log level SC_LP_ESSENTIAL. See sc.h for more info on the log levels. */
+  t8_init (SC_LP_PRODUCTION);
+
+  /** Little Test. */
+  Eigen::Matrix3d A;
+  A << 4.0, 1.2, 0.5, 1.2, 5.0, 2.1, 0.5, 2.1, 6.0;
+
+  /** vector b. */
+  Eigen::Vector3d b;
+  b << 1.0, 2.0, 3.0;
+
+  /* We will test a dense matrix example. Given that the matrix is symmetric and positive semidefinite we will test out the LDLT decomposition. */
+  Eigen::LDLT<Eigen::Matrix3d> solver (A);
+
+  if (solver.info () == Eigen::Success) {
+    Eigen::Vector3d x = solver.solve (b);
+
+    t8_global_productionf ("The solution is: x = [%f, %f, %f]\n", x (0), x (1), x (2));
+
+    /* Test if A * x = b for a quick check. */
+    Eigen::Vector3d check = A * x;
+    t8_global_productionf ("(A*x): [%f, %f, %f] (should be the same as b)\n", check (0), check (1), check (2));
+  }
+  else {
+    t8_global_productionf ("ERROR: the system could not be solved.\n");
+  }
+
+  sc_finalize ();
+  mpiret = sc_MPI_Finalize ();
+  SC_CHECK_MPI (mpiret);
+
+#else
+  t8_global_productionf ("ERROR: This example requires Eigen support to be enabled in t8code.\n");
+#endif
+
+  return 0;
+}

src/CMakeLists.txt

@@ -97,7 +97,6 @@ if( T8CODE_ENABLE_OCC )
     target_sources(T8 PRIVATE
       t8_geometry/t8_geometry_implementations/t8_geometry_cad.cxx
       t8_cad/t8_cad_handle.cxx
-      t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.cxx
     )
     install( FILES
       t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx
@@ -108,6 +107,24 @@ if( T8CODE_ENABLE_OCC )
       t8_cad/t8_cad_handle.hxx
       DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/t8_cad
     )
+endif()
+
+if( T8CODE_ENABLE_EIGEN )
+    target_compile_definitions( T8 PUBLIC T8_ENABLE_EIGEN=1 )
+    target_link_libraries( T8 PRIVATE Eigen3::Eigen )
+    target_sources(T8 PRIVATE
+      t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx
+    )
+    install( FILES
+      t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx
+      DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/t8_cmesh/t8_cmesh_mesh_deformation/
+    )
+endif()
+
+if( T8CODE_ENABLE_OCC AND T8CODE_ENABLE_EIGEN )
+    target_sources(T8 PRIVATE
+      t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.cxx
+    )
     install( FILES
       t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.hxx
       DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/t8_cmesh/t8_cmesh_mesh_deformation/

src/config.cmake.in

@@ -13,6 +13,7 @@ set( T8CODE_BUILD_MESH_HANDLE @T8CODE_BUILD_MESH_HANDLE@ )
 
 set( T8CODE_ENABLE_MPI @T8CODE_ENABLE_MPI@ )
 set( T8CODE_ENABLE_VTK @T8CODE_ENABLE_VTK@ )
+set( T8CODE_ENABLE_EIGEN @T8CODE_ENABLE_EIGEN@ )
 
 set( T8CODE_USE_SYSTEM_SC @T8CODE_USE_SYSTEM_SC@ )
 set( T8CODE_USE_SYSTEM_P4EST @T8CODE_USE_SYSTEM_P4EST@ )
@@ -25,6 +26,11 @@ if(T8CODE_ENABLE_VTK)
     set (T8CODE_VTK_VERSION_USED "@T8CODE_VTK_VERSION_USED@")
 endif()
 
+# Ensure that external libraries using for example find_package ( t8code REQUIRED) link automatically against eigen
+if(T8CODE_ENABLE_EIGEN)
+    find_dependency(Eigen3)
+endif()
+
 include( "${CMAKE_CURRENT_LIST_DIR}/@PROJECT_NAME@-targets.cmake" )
 
 check_required_components( @PROJECT_NAME@ )