diff --git a/CMakeLists.txt b/CMakeLists.txt index 1e425bc63c..ea7af319e7 100644 --- a/CMakeLists.txt +++ b/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 ) diff --git a/cmake/thirdparty.cmake b/cmake/thirdparty.cmake index 34e5b13110..378b9add4f 100644 --- a/cmake/thirdparty.cmake +++ b/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() diff --git a/cmake/thirdparty.json b/cmake/thirdparty.json index e65aafcf41..58ab6d253c 100644 --- a/cmake/thirdparty.json +++ b/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" + } } ] -} +} \ No newline at end of file diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index e7b7355d26..3aff6b553a 100644 --- a/example/CMakeLists.txt +++ b/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 ) diff --git a/example/cmesh/t8_cmesh_mesh_deformation.cxx b/example/cmesh/t8_cmesh_mesh_deformation.cxx index 9dd882833c..c4b9342d9f 100644 --- a/example/cmesh/t8_cmesh_mesh_deformation.cxx +++ b/example/cmesh/t8_cmesh_mesh_deformation.cxx @@ -28,10 +28,10 @@ #include #include #include -#if T8CODE_ENABLE_OCC +#if T8_ENABLE_OCC && T8_ENABLE_EIGEN #include #include -#endif /* T8CODE_ENABLE_OCC */ +#endif /* T8_ENABLE_OCC and T8_ENABLE_EIGEN*/ #include #include @@ -39,11 +39,32 @@ #include #include #include +#include +#include + +#if T8_ENABLE_OCC && T8_ENABLE_EIGEN +namespace fs = std::filesystem; + +static std::vector +findBrepFiles (const char *folder) +{ + std::vector 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 (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 (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 (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; } diff --git a/example/cmesh/t8_rbf.cxx b/example/cmesh/t8_rbf.cxx new file mode 100644 index 0000000000..e97834fa37 --- /dev/null +++ b/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 +#if T8_ENABLE_EIGEN +#include +#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 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; +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5225ccc1a1..c6ff75645d 100644 --- a/src/CMakeLists.txt +++ b/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/ diff --git a/src/config.cmake.in b/src/config.cmake.in index c8fa76969b..8b97473f87 100644 --- a/src/config.cmake.in +++ b/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@ ) diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.cxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.cxx index 2356b73a60..752dbeb680 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.cxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.cxx @@ -25,6 +25,8 @@ */ #include +#include +#include #include #include #include @@ -33,15 +35,95 @@ #include #include -std::unordered_map -t8_cmesh_mesh_deformation::calculate_displacement_surface_vertices (const t8_cad_handle *cad) +static double +calculate_local_support_radius (t8_cmesh_t cmesh, const std::vector> &tree_list, + int current_global_vertex_id, const t8_3D_vec &displacements) +{ + double total_distance = 0.0; + + double displacement_magnitude = std::sqrt (displacements[0] * displacements[0] + displacements[1] * displacements[1] + + displacements[2] * displacements[2]); + /** Get the first tree where this vertex exists. */ + const auto &first_tree = tree_list.front (); + const t8_locidx_t first_tree_id = first_tree.first; + const int first_local_index = first_tree.second; + + const double *first_tree_coords = (const double *) t8_cmesh_get_attribute ( + cmesh, t8_get_package_id (), T8_CMESH_VERTICES_ATTRIBUTE_KEY, first_tree_id); + + /** Check if the coordinates are available. */ + if (first_tree_coords == nullptr) { + t8_errorf ("Error: Coordinates attribute missing for tree %d\n.", first_tree_id); + SC_ABORTF ("Vertex coordinates are missing."); + } + + const double current_vertex_x = first_tree_coords[3 * first_local_index + 0]; + const double current_vertex_y = first_tree_coords[3 * first_local_index + 1]; + const double current_vertex_z = first_tree_coords[3 * first_local_index + 2]; + + /** Set up a set to calculate the distance to every neighbors vertex just once. */ + std::set> unique_neighbors; + + /** Iterate over the tree list in which the current vertex is present. */ + for (const auto &[tree_id, local_vertex_index_of_the_current_vertex] : tree_list) { + + /** Get the vertex coordinates array of the current tree. */ + const double *tree_vertex_coords + = (const double *) t8_cmesh_get_attribute (cmesh, t8_get_package_id (), T8_CMESH_VERTICES_ATTRIBUTE_KEY, tree_id); + + /** Check if the coordinates are available. */ + if (tree_vertex_coords != nullptr) { + + int num_vertices = t8_eclass_num_vertices[t8_cmesh_get_tree_class (cmesh, tree_id)]; + + for (int neighbor_vertex = 0; neighbor_vertex < num_vertices; ++neighbor_vertex) { + /** If the vertex is the current vertex itself, we do not calculate the distance. */ + if (neighbor_vertex != local_vertex_index_of_the_current_vertex) { + + /** Get the coordinates of the neighbor vertex. */ + double x = tree_vertex_coords[3 * neighbor_vertex + 0]; + double y = tree_vertex_coords[3 * neighbor_vertex + 1]; + double z = tree_vertex_coords[3 * neighbor_vertex + 2]; + + /** Save the coordinates in the set, in which double coordinates will be filtered out. */ + unique_neighbors.insert (std::make_tuple (x, y, z)); + } + } + } + } + + /** Check if the set is empty. */ + if (unique_neighbors.empty ()) { + t8_errorf ("Error: No neighbor vertices found for global vertex %d to calculate local support radius.\n", + current_global_vertex_id); + SC_ABORTF ("Calculation of local support radius failed due to missing neighbor vertices."); + } + /** Calculate the distance to every neighbor node. */ + for (const auto &[x, y, z] : unique_neighbors) { + double distance_x = current_vertex_x - x; + double distance_y = current_vertex_y - y; + double distance_z = current_vertex_z - z; + + total_distance += std::sqrt (distance_x * distance_x + distance_y * distance_y + distance_z * distance_z); + } + + double grading_factor = 2.0; + + return (total_distance / static_cast (unique_neighbors.size ())) + * (1.0 + grading_factor * displacement_magnitude); +} + +std::unordered_map +t8_cmesh_mesh_deformation::calculate_displacement_surface_vertices (const t8_cad_handle *cad, + const t8_rbf_function_type rbf_type, + const double scale_factor_support_radius) { T8_ASSERT (t8_cmesh_is_committed (associated_cmesh)); const int mesh_dimension = t8_cmesh_get_dimension (associated_cmesh); /* Map from global vertex id -> displacement vector. */ - std::unordered_map displacements; + std::unordered_map boundary_node_data; for (const auto &global_vertex : *(associated_cmesh->vertex_connectivity)) { @@ -98,7 +180,7 @@ t8_cmesh_mesh_deformation::calculate_displacement_surface_vertices (const t8_cad /* Check if the (u,v)-parameters are available. */ if (uv_attribute == nullptr) { t8_errorf ("Error: (u,v)-parameters are missing for tree %d\n.", first_tree_id); - SC_ABORT ("(u,v)-parameters are missing.\n"); + SC_ABORTF ("(u,v)-parameters are missing.\n"); } /* Get the (u,v)-parameter of the vertex. */ const double *uv_parameter = &uv_attribute[2 * local_corner_index]; @@ -115,7 +197,7 @@ t8_cmesh_mesh_deformation::calculate_displacement_surface_vertices (const t8_cad gp_Pnt new_coords; - /* Find the new coordinates of the vertex in the cad file, based on the geometry its lying on. */ + /* Find the new coordinates of the vertex in the CAD file, based on the geometry it's lying on. */ switch (first_tree_entity_dim) { case 0: { new_coords = cad->get_cad_point (first_tree_entity_tag); @@ -140,43 +222,113 @@ t8_cmesh_mesh_deformation::calculate_displacement_surface_vertices (const t8_cad const double old_y = old_coords[3 * local_corner_index + 1]; const double old_z = old_coords[3 * local_corner_index + 2]; + t8_rbf_boundary_node node; + + node.position = { old_x, old_y, old_z }; + /* Calculate the displacement of the vertex which should be then done in the deformation. */ - displacements[global_vertex_id] = { new_coords.X () - old_x, new_coords.Y () - old_y, new_coords.Z () - old_z }; + node.displacement = { new_coords.X () - old_x, new_coords.Y () - old_y, new_coords.Z () - old_z }; + + node.weight.fill (0.0); + + if (rbf_type == T8_RBF_CP_C2) { + node.local_support_radius + = scale_factor_support_radius + * calculate_local_support_radius (associated_cmesh, tree_list, global_vertex_id, node.displacement); + } + boundary_node_data[global_vertex_id] = node; } } - return displacements; + return boundary_node_data; } void -t8_cmesh_mesh_deformation::apply_vertex_displacements (const std::unordered_map &displacements, - std::shared_ptr cad) +t8_cmesh_mesh_deformation::apply_vertex_displacements ( + std::unordered_map &boundary_node_data, std::shared_ptr cad, + const t8_rbf_function_type rbf_type) { T8_ASSERT (t8_cmesh_is_committed (associated_cmesh)); - /* Iterate over all vertices in the displacement map. */ - for (const auto &[global_vertex, displacement] : displacements) { + t8_rbf rbf_handler (rbf_type); - /* Get the list of trees where this vertex exists. */ - const auto &tree_list = associated_cmesh->vertex_connectivity->get_tree_list_of_vertex (global_vertex); + rbf_handler.set_boundary_nodes (std::move (boundary_node_data)); - /* Update the vertex coordinates in each tree. */ - for (const auto &[tree_id, local_vertex_index] : tree_list) { + /** Calculate the weights. */ + rbf_handler.solve (); + + /** Iterate over all vertices in the displacement map. */ + for (const auto &global_vertex : *(associated_cmesh->vertex_connectivity)) { + /** Get the global vertex ID of the current vertex. */ + t8_gloidx_t global_vertex_id = global_vertex.first; + + /** Get the list of trees where this vertex exists. */ + const auto &tree_list = global_vertex.second; + + /** Get the first tree where this vertex exists as a reference. */ + const auto &first_tree = tree_list.front (); - /* Get the vertex coordinates of the current tree. */ + /* Get the data of the first tree. */ + const int *first_tree_geom_attribute = static_cast (t8_cmesh_get_attribute ( + associated_cmesh, t8_get_package_id (), T8_CMESH_NODE_GEOMETRY_ATTRIBUTE_KEY, first_tree.first)); + + /* Check if the geometry attribute is available for this tree. */ + if (first_tree_geom_attribute == nullptr) { + t8_errorf ("Error: Geometry attribute missing for tree %d\n.", first_tree.first); + SC_ABORTF ("Geometry attribute is missing."); + } + + const int first_tree_entity_dim = first_tree_geom_attribute[2 * tree_list[0].second]; + + const int mesh_dimension = t8_cmesh_get_dimension (associated_cmesh); + + /** Contains the displacement whether its a boundary node or an inner node. */ + t8_3D_vec final_displacement; + + /* Check if this vertex is a boundary node. + If so, we can use the already known displacement extracted from the new CAD geometry input given.*/ + if (first_tree_entity_dim < mesh_dimension && first_tree_entity_dim >= 0) { + + final_displacement = rbf_handler.get_boundary_displacement (global_vertex_id); + } + /** If it is an inner node, we do not know the displacement right away and need to interpolate to get the new coordinates. */ + else { + + const double *vertex_coords = static_cast (t8_cmesh_get_attribute ( + associated_cmesh, t8_get_package_id (), T8_CMESH_VERTICES_ATTRIBUTE_KEY, first_tree.first)); + /** Check if the coordinates are available. */ + if (vertex_coords == nullptr) { + t8_errorf ("Error: Coordinates attribute missing for tree %d\n.", first_tree.first); + SC_ABORTF ("Vertex coordinates are missing."); + } + + /** Save the initial coordinates of the inner node. */ + t8_rbf_node inner_node; + inner_node.global_id = global_vertex_id; + inner_node.position = { vertex_coords[3 * first_tree.second], vertex_coords[3 * first_tree.second + 1], + vertex_coords[3 * first_tree.second + 2] }; + + /** Calculate the displacement of the inner node. */ + rbf_handler.interpolate (inner_node); + final_displacement = inner_node.displacement; + } + + /** Update the vertex coordinates in each tree. */ + for (const auto &[tree_id, local_vertex_index] : tree_list) { + /** Get the vertex coordinates of the current tree. */ double *tree_vertex_coords = (double *) t8_cmesh_get_attribute (associated_cmesh, t8_get_package_id (), T8_CMESH_VERTICES_ATTRIBUTE_KEY, tree_id); - /* Check if the coordinates are available. */ + /** Check if the coordinates are available. */ if (tree_vertex_coords != nullptr) { - /* Update the coordinates of the vertex. */ + /** Update the coordinates of the vertex. */ for (int coord_index = 0; coord_index < 3; ++coord_index) { - tree_vertex_coords[3 * local_vertex_index + coord_index] += displacement[coord_index]; + tree_vertex_coords[3 * local_vertex_index + coord_index] += final_displacement[coord_index]; } } } } - /* Update the cad geometry. */ + /** Update the cad geometry. */ t8_geometry_handler *geometry_handler = associated_cmesh->geometry_handler; T8_ASSERT (geometry_handler != nullptr); diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.hxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.hxx index 229f2c4fee..9f0b7196c2 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.hxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.hxx @@ -35,6 +35,8 @@ #include #include +#include + /** Struct for mesh deformation. */ struct t8_cmesh_mesh_deformation { @@ -49,10 +51,13 @@ struct t8_cmesh_mesh_deformation * Computes the displacements of the surface vertices. * * \param [in] cad A pointer to the CAD-based geometry object. - * \return Map from global vertex ID to 3D displacement vector + * \param [in] rbf_type The type of radial basis function to be used. + * \param [in] scale_factor_support_radius The scale factor for the support radius. + * \return Map from global vertex ID to RBF boundary node which contains the displacement and can than be used to calculate the weight of the boundary node. */ - std::unordered_map - calculate_displacement_surface_vertices (const t8_cad_handle *cad); + std::unordered_map + calculate_displacement_surface_vertices (const t8_cad_handle *cad, const t8_rbf_function_type rbf_type, + const double scale_factor_support_radius); /** * Apply vertex displacements to a committed cmesh. @@ -60,12 +65,13 @@ struct t8_cmesh_mesh_deformation * Iterates over the provided map of global vertex IDs to 3D displacement vectors, * updating the coordinates in each tree where the vertex appears. * - * \param [in] displacements Map from global vertex ID to 3D displacement vector [dx, dy, dz]. + * \param [in] boundary_node_data Map from global vertex ID to RBF boundary node. * \param [in] cad The shared pointer to the CAD geometry to update. + * \param [in] rbf_type The RBF type to be used. */ void - apply_vertex_displacements (const std::unordered_map &displacements, - std::shared_ptr cad); + apply_vertex_displacements (std::unordered_map &boundary_node_data, + std::shared_ptr cad, const t8_rbf_function_type rbf_type); private: /** A pointer to the cmesh for attribute retrieval */ diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx new file mode 100644 index 0000000000..6d21e2c5bb --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx @@ -0,0 +1,132 @@ +/* + 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 the Radial Basis Functions for the mesh deformation. + */ + +#include +#include +#include +#include +#include +#include + +/** + * Solves the linear system A * alpha = displacements to find the weight. + * This step is mandatory to later be able to interpolate the inner nodes. + */ +void +t8_rbf::solve () +{ + const size_t num_boundary_nodes = boundary_nodes.size (); + /** Check if there are any boundary nodes. */ + if (num_boundary_nodes == 0) { + t8_errorf ("ERROR: Boundary nodes are not added correctly. The current RBF instance has no boundary nodes\n."); + } + + /** Fill the displacement vector with the values of the boundary nodes' displacements. */ + Eigen::MatrixXd displacements (num_boundary_nodes, 3); + for (size_t i = 0; i < num_boundary_nodes; ++i) { + displacements.row (i) = Eigen::Vector3d (boundary_nodes[i].displacement[0], boundary_nodes[i].displacement[1], + boundary_nodes[i].displacement[2]); + } + /** Weight vector. */ + Eigen::MatrixXd alpha; + + if (rbf_function->is_compactly_supported ()) { + alpha = solve_compactly_supported_rbf (displacements, num_boundary_nodes); + } + else { + alpha = solve_globally_supported_rbf (displacements, num_boundary_nodes); + } + /** Copy the calculated weights back to the boundary nodes. */ + for (size_t i = 0; i < num_boundary_nodes; ++i) { + boundary_nodes[i].weight[0] = alpha (i, 0); + boundary_nodes[i].weight[1] = alpha (i, 1); + boundary_nodes[i].weight[2] = alpha (i, 2); + } +} + +Eigen::MatrixXd +t8_rbf::solve_compactly_supported_rbf (const Eigen::MatrixXd &displacements, const size_t num_boundary_nodes) const +{ + /** The RBF used is compactly supported so we can use a sparse matrix. The conjugate gradient method can be used for solving the linear equation system. */ + size_t estimation_of_entries = 50; + typedef Eigen::Triplet triplet; + std::vector coefficients; + coefficients.reserve (num_boundary_nodes * estimation_of_entries); + /** Fill the matrix A with the values of the radial basis function psi evaluated at the pairwise euclidean distances between all boundary nodes. */ + for (size_t row = 0; row < num_boundary_nodes; ++row) { + + for (size_t col = 0; col < num_boundary_nodes; ++col) { + /** Calculate the distance between the current pair of boundary nodes. */ + const double distance = t8_dist (boundary_nodes[row].position, boundary_nodes[col].position); + /** Evaluate the radial basis function for the current distance. */ + const double psi = rbf_function->evaluate (distance, boundary_nodes[col].local_support_radius); + /** Check if the basis function value (the influence factor) is not equal to zero with a numerical tolerance of 1e-12. */ + if (std::abs (psi) > 1e-12) { + coefficients.emplace_back (row, col, psi); + } + } + } + /** Fill the sparse matrix A with the triplets. */ + Eigen::SparseMatrix A (num_boundary_nodes, num_boundary_nodes); + A.setFromTriplets (coefficients.begin (), coefficients.end ()); + /** Solve the linear system using the conjugate gradient method. */ + Eigen::BiCGSTAB> solver; + solver.setTolerance (1e-10); + solver.compute (A); + + if (solver.info () != Eigen::Success) { + t8_errorf ("ERROR: Decomposition of the matrix A failed. The linear system cannot be solved.\n"); + } + + return solver.solve (displacements); +} + +Eigen::MatrixXd +t8_rbf::solve_globally_supported_rbf (const Eigen::MatrixXd &displacements, const size_t num_boundary_nodes) const +{ + /** The RBF used is globally supported so we need to use a dense matrix. The linear equation system can be solved with the built-in solver of Eigen. */ + + /** Create the matrix A for the linear system. */ + Eigen::MatrixXd A (num_boundary_nodes, num_boundary_nodes); + /** Fill the matrix A with the values of the radial basis function psi evaluated at the pairwise euclidean distances between all boundary nodes. */ + for (size_t row = 0; row < num_boundary_nodes; ++row) { + /** Because of the symmetric property of the distance between nodes, we only need to compute the upper triangular part of the matrix + * and can mirror the values to the lower triangular part. */ + for (size_t col = 0; col < num_boundary_nodes; ++col) { + /** Calculate the distance between the current pair of boundary nodes. */ + const double distance = t8_dist (boundary_nodes[row].position, boundary_nodes[col].position); + /** Evaluate the radial basis function for the current distance and write the value to the matrix A. */ + A (row, col) = rbf_function->evaluate (distance, 0.0); + } + } + Eigen::PartialPivLU solver (A); + + if (solver.info () != Eigen::Success) { + t8_errorf ("ERROR: The global RBF system could not be solved with LU-Decomposition.\n"); + } + + return solver.solve (displacements); +} diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx new file mode 100644 index 0000000000..9a7c270d79 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx @@ -0,0 +1,288 @@ +/* + 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.hxx + * Implementation of CAD-based mesh deformation. + */ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +/** + * The available RBF function types. + */ +typedef enum { + T8_RBF_CP_C2 = 0, /**< Compactly supported radial basis function (CPC2). */ + T8_RBF_TPS, /**< Globally supported radial basis function (TPS). */ + T8_RBF_FUNCTION_COUNT /**< The number of available RBF function types. */ +} t8_rbf_function_type; + +/** A struct representing a node in the mesh for the RBF interpoltaion, containing its global ID, position and displacement. */ +struct t8_rbf_node +{ + /**The global ID of the node. */ + t8_gloidx_t global_id; + /** The old position of the node. */ + t8_3D_vec position; + /** The displacement. + * For boundary nodes this is known and calculated from the different coordinates of the CAD geometries. + * For inner nodes this is the result of the interpolation. */ + t8_3D_vec displacement; +}; + +/** + * This struct will be a single boundary node from the mesh deformation. + * The new position of this node can be directly calculated from the new incoming CAD geometry given. + */ +struct t8_rbf_boundary_node: public t8_rbf_node +{ + /** The calculated RBF coefficient alpha. This is needed for the interpolation to distribute the influence of the boundary nodes. + * This weight represents the fixed influence assigned to each boundary node. */ + t8_3D_vec weight; + /** The local support radius for the boundary node. */ + double local_support_radius; +}; + +/** A base struct representing a radial basis function. */ +struct t8_rbf_function +{ + /** Destructor. */ + virtual ~t8_rbf_function () {}; + + /** + * Evaluate the radial basis function. + * \param[in] distance The euclidean distance between two points. + * \param[in] radius The support radius. + * \return The function value. + */ + virtual double + evaluate (double distance, double radius) const + = 0; + + /** + * Check if the radial basis function is compactly supported. + * \return true if the function is compactly supported, false otherwise. + */ + virtual bool + is_compactly_supported () const + = 0; +}; + +/** + * The CP C2 radial basis function which is compactly supported and will be used in the + * mesh deformation to move the inner nodes from the known movement of the boundary nodes. + */ +struct t8_rbf_cpc2: public t8_rbf_function +{ + /** Constructor. */ + t8_rbf_cpc2 () + { + } + /** Destructor. */ + ~t8_rbf_cpc2 () {}; + /** + * Solve the radial basis function. + * Formula: psi(x) = (1 - x)^4 * (4x + 1) for (1-x) > 0. + * where x = distance / support radius. The distance is the euclidean distance between two points. + * \param [in] distance The euclidean distance between two points. + * \param [in] radius The support radius. + * \return The function value psi. It returns 0.0 if the node is out of the chosen radius and so has no impact. + */ + double + evaluate (double distance, double radius) const override + { + double r = distance / radius; + if (r < 1.0) { + double d = 1.0 - r; + return (d * d * d * d) * (4.0 * r + 1.0); + } + return 0.0; + } + + /** + * The CPC2 is a compactly supported RBF. + * \return true. + */ + bool + is_compactly_supported () const override + { + return true; + } +}; + +/** + * The TPS radial basis function with global support. + */ +struct t8_rbf_tps: public t8_rbf_function +{ + /** Constructor. The TPS RBF does not have a support radius, because it is globally supported. + */ + t8_rbf_tps () + { + } + /** Destructor. */ + ~t8_rbf_tps () {}; + /** + * Solve the radial basis function. + * Formula: psi(x) = x^2 * log(x). + * where x is the euclidean distance between two points. + * \param[in] distance The euclidean distance between two points. + * \param [in] radius The support radius (not used for TPS because it is globally supported). + * \return The function value psi. + */ + double + evaluate (double distance, [[maybe_unused]] double radius) const override + { + if (distance < 1e-12) { + return 0.0; + } + return distance * distance * std::log (distance); + } + + /** + * The TPS is a globally supported RBF. + * \return false. + */ + bool + is_compactly_supported () const override + { + return false; + } +}; + +/** + * Struct for mesh deformation using Radial Basis Functions (RBF). + * It handles the interpolation of the boundary displacements to the inner nodes. + */ +struct t8_rbf +{ + public: + /** Constructor. + * \param [in] rbf_type The RBF type to be used. + */ + t8_rbf (t8_rbf_function_type rbf_type) + { + + if (rbf_type == T8_RBF_CP_C2) { + rbf_function = std::make_unique (); + } + else if (rbf_type == T8_RBF_TPS) { + rbf_function = std::make_unique (); + } + else { + SC_ABORTF ("ERROR: RBF attribute missing or not correct. Unsupported RBF type.\n"); + } + }; + + /** Destructor. */ + ~t8_rbf () {}; + + /** + * Transfers the boundary node data to the internal data structure of the RBF. + * \param [in] boundary_node_data Map of global vertex IDs to boundary nodes. + */ + void + set_boundary_nodes (std::unordered_map&& boundary_node_data) + { + boundary_nodes.clear (); + boundary_nodes.reserve (boundary_node_data.size ()); + for (auto& [global_vertex_id, boundary_node] : boundary_node_data) { + + boundary_node.global_id = global_vertex_id; + + boundary_nodes.push_back (std::move (boundary_node)); + } + } + + /** + * Search for the displacement of a specific boundary node with its associated global ID. + * \param [in] global_id The global ID of the boundary node. + * \return The displacement vector of the boundary node. + */ + t8_3D_vec + get_boundary_displacement (const t8_gloidx_t global_id) const + { + for (const auto& node : boundary_nodes) { + if (node.global_id == global_id) { + return node.displacement; + } + } + /* If the boundary node can not be found in the list of boundary nodes.*/ + t8_errorf ("ERROR: The boundary node %ld is missing in the boundary node list.\n", global_id); + SC_ABORTF ("A boundary node is not recognized as one."); + } + + /** + * Solves the linear system A * alpha = displacements to find the weight. + * This step is mandatory to later be able to interpolate the inner nodes. + */ + void + solve (); + + /** + * Interpolates the inner node. + * \param [in, out] inner_node The inner node which will be interpolated. + */ + void + interpolate (t8_rbf_node& inner_node) const + { + /** Reset the inner_node_displacement to zero. */ + inner_node.displacement.fill (0.0); + /** Iterate over all boundary nodes. */ + for (const auto& boundary_node : boundary_nodes) { + + double distance = t8_dist (inner_node.position, boundary_node.position); + + const double psi = rbf_function->evaluate (distance, boundary_node.local_support_radius); + /** Check if the basis function value (the influence factor) is not equal to zero with a numerical tolerance of 1e-12. */ + if (std::abs (psi) > 1e-12) { + + for (int coordinate = 0; coordinate < 3; ++coordinate) { + /** Update the displacement for each coordinate with the weighted influence of the boundary node. */ + inner_node.displacement[coordinate] += boundary_node.weight[coordinate] * psi; + } + } + } + } + + private: + /** The specialised solve function for either compactly supported or globally supported RBFs. + * \param [in] displacements The matrix of the boundary node displacements. Each row corresponds to a + * boundary node and the three columns correspond to the x, y, and z components of the displacement. + * \param [in] num_boundary_nodes The number of boundary nodes. + * \return The matrix of the calculated weights alpha. + */ + Eigen::MatrixXd + solve_compactly_supported_rbf (const Eigen::MatrixXd& displacements, const size_t num_boundary_nodes) const; + Eigen::MatrixXd + solve_globally_supported_rbf (const Eigen::MatrixXd& displacements, const size_t num_boundary_nodes) const; + /** List of all registered support nodes. */ + std::vector boundary_nodes; + /** Pointer to the radial basis function. */ + std::unique_ptr rbf_function; +}; diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx index f3dc67aba4..d4b178a4ea 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx @@ -158,7 +158,7 @@ struct t8_geometry_cad: public t8_geometry_with_vertices /** * Getter function for the CAD handle. - * + * * \return The CAD handle of the geometry. */ std::shared_ptr @@ -167,7 +167,7 @@ struct t8_geometry_cad: public t8_geometry_with_vertices return cad_handle; } - /** Update the CAD handle with a new one. + /** Update the CAD handle with a new one. * \param[in] new_cad_handle The new CAD handle to be used. */ void