From 0eead4d1bd4f5150d9900391ec3c99c01b6d3e90 Mon Sep 17 00:00:00 2001 From: Lena Radmer Date: Thu, 26 Mar 2026 14:34:20 +0100 Subject: [PATCH 01/15] include mesh deformation example (commit 78fd6b2) --- example/cmesh/t8_cmesh_mesh_deformation.cxx | 156 ++++++++++++++++++++ 1 file changed, 156 insertions(+) create mode 100644 example/cmesh/t8_cmesh_mesh_deformation.cxx diff --git a/example/cmesh/t8_cmesh_mesh_deformation.cxx b/example/cmesh/t8_cmesh_mesh_deformation.cxx new file mode 100644 index 0000000000..0e0d4aad6b --- /dev/null +++ b/example/cmesh/t8_cmesh_mesh_deformation.cxx @@ -0,0 +1,156 @@ +/* + 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_cmesh_mesh_deformation.cxx + * This file implements an example for CAD-based mesh deformation. + */ + +#include +#include +#include +#include +#if T8CODE_ENABLE_OCC +#include +#include +#endif +#include +#include + +#include +#include +#include +#include + +int +main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) +{ +#if T8CODE_ENABLE_OCC + + char usage[BUFSIZ]; + /* Brief help message. */ + int sreturnA = snprintf (usage, BUFSIZ, "Usage:\t%s \n\t%s -h\t for a brief overview of all options.", + basename (argv[0]), basename (argv[0])); + + char help[BUFSIZ]; + /* Long help message. */ + int sreturnB = snprintf ( + help, BUFSIZ, + "Deform a mesh based on a msh file with the new CAD geometry.\n" + "Required arguments are the input mesh file, the deformation geometry file, and the mesh dimension.\n\n%s\n", + usage); + + if (sreturnA > BUFSIZ || sreturnB > BUFSIZ) { + /* The usage string or help message was truncated */ + /* Note: gcc >= 7.1 prints a warning if we + * do not check the return value of snprintf. */ + t8_debugf ("WARNING: Truncated usage string and help message to '%s' and '%s'\n", usage, help); + } + + /* + * 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_ESSENTIAL); + + /* Initialize t8code with log level SC_LP_ESSENTIAL. See sc.h for more info on the log levels. */ + t8_init (SC_LP_ESSENTIAL); + + int helpme = 0; + const char *msh_file = NULL; + const char *brep_file = NULL; + int dim, level; + + /* 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)"); + 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"); + + int parsed = sc_options_parse (t8_get_package_id (), SC_LP_ERROR, opt, argc, argv); + + if (helpme) { + t8_global_productionf ("%s\n", help); + sc_options_print_usage (t8_get_package_id (), SC_LP_ERROR, opt, NULL); + } + else if (msh_file == NULL || brep_file == NULL || dim == 0) { + t8_global_errorf ("ERROR: Missing required arguments: -m, -b, and -d are mandatory.\n\n"); + 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); + sc_options_print_usage (t8_get_package_id (), SC_LP_ERROR, opt, NULL); + } + else if (parsed >= 0) { + + /* We will use MPI_COMM_WORLD as a communicator. */ + sc_MPI_Comm comm = sc_MPI_COMM_WORLD; + + /* Create cmesh from msh. */ + t8_cmesh_t cmesh = t8_cmesh_from_msh_file (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 ()); + + /* Write output. */ + t8_forest_vtk_write_file (forest, "input_forest", 1, 1, 1, 1, 0, 0, NULL); + + /* Apply displacements. */ + deformation.apply_vertex_displacements (displacements, cad); + + /* Write output. */ + t8_forest_vtk_write_file (forest, "deformed_forest", 1, 1, 1, 1, 0, 0, NULL); + + /* Cleanup. */ + t8_forest_unref (&forest); + + t8_global_productionf ("Mesh deformation completed."); + } + + sc_options_destroy (opt); + + sc_finalize (); + mpiret = sc_MPI_Finalize (); + SC_CHECK_MPI (mpiret); + +#else + t8_global_errorf ("ERROR: This example requires OpenCASCADE support to be enabled in t8code.\n\n"); +#endif + + return 0; +} From 8759747b19e591f55fb31b5a4f2b66c8ac7e8647 Mon Sep 17 00:00:00 2001 From: Lena Radmer Date: Thu, 26 Mar 2026 14:45:55 +0100 Subject: [PATCH 02/15] import other mesh deformation logic as well --- example/CMakeLists.txt | 1 + src/CMakeLists.txt | 5 + src/t8_cad/t8_cad_handle.hxx | 1 + .../t8_cmesh_mesh_deformation.cxx | 189 ++++++++++++++++++ .../t8_cmesh_mesh_deformation.hxx | 75 +++++++ .../t8_cmesh_vertex_connectivity.hxx | 21 ++ src/t8_geometry/t8_geometry_handler.hxx | 20 ++ .../t8_geometry_cad.hxx | 13 +- 8 files changed, 323 insertions(+), 2 deletions(-) create mode 100644 src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.cxx create mode 100644 src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.hxx diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index 272486619b..a897b4d314 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -71,6 +71,7 @@ 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 ) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a9a421db91..8f65b2a588 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -97,6 +97,7 @@ 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 @@ -107,6 +108,10 @@ if( T8CODE_ENABLE_OCC ) t8_cad/t8_cad_handle.hxx DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/t8_cad ) + install( FILES + t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.hxx + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/t8_cmesh/t8_cmesh_mesh_deformation/ + ) endif() if( T8CODE_BUILD_PEDANTIC ) diff --git a/src/t8_cad/t8_cad_handle.hxx b/src/t8_cad/t8_cad_handle.hxx index ec7c3c3809..f119132029 100644 --- a/src/t8_cad/t8_cad_handle.hxx +++ b/src/t8_cad/t8_cad_handle.hxx @@ -55,6 +55,7 @@ class t8_cad_handle { * \param [in] fileprefix Prefix of a .brep file from which to extract cad geometry. */ t8_cad_handle (const std::string_view fileprefix); + /** * Constructor of the cad shape. * The shape is initialized directly from an existing TopoDS_Shape. 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 new file mode 100644 index 0000000000..21e593904a --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.cxx @@ -0,0 +1,189 @@ +/* + 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_cmesh_mesh_deformation.cxx + * This file implements the routines for CAD-based mesh deformation. + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +std::unordered_map +t8_cmesh_mesh_deformation::calculate_displacement_surface_vertices (const t8_cad_handle *cad) +{ + 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; + + for (const auto &global_vertex : *(associated_cmesh->vertex_connectivity)) { + + /* Get the list of all trees associated with the vertex. */ + const auto &tree_list = global_vertex.second; + const t8_gloidx_t global_vertex_id = global_vertex.first; + + /* Get the first tree and the local corner index of the vertex in the tree. */ + const auto &first_tree = tree_list.front (); + const t8_locidx_t first_tree_id = first_tree.first; + const int local_corner_index = first_tree.second; + + /* Get the first tree as a reference. */ + 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_id)); + + /* 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_id); + SC_ABORTF ("Geometry attribute is missing."); + } + + const int first_tree_entity_dim = first_tree_geom_attribute[2 * tree_list[0].second]; + const int first_tree_entity_tag = first_tree_geom_attribute[2 * tree_list[0].second + 1]; + +#if T8_ENABLE_DEBUG + /* Iterate over all trees and compare to the reference tree. */ + for (const auto &[tree_id, local_corner_index] : tree_list) { + + const int *geom_attribute = static_cast ( + t8_cmesh_get_attribute (associated_cmesh, t8_get_package_id (), T8_CMESH_NODE_GEOMETRY_ATTRIBUTE_KEY, tree_id)); + + const int entity_dim = geom_attribute[2 * local_corner_index]; + const int entity_tag = geom_attribute[2 * local_corner_index + 1]; + + /* Check if the attribute of the vertex is the same in all trees. */ + if (!(entity_dim == first_tree_entity_dim && entity_tag == first_tree_entity_tag)) { + t8_errorf ( + "Error: Inconsistent entity info for global vertex %li: tree %d: dim=%d tag=%d, expected dim=%d tag=%d\n", + global_vertex_id, tree_id, entity_dim, entity_tag, first_tree_entity_dim, first_tree_entity_tag); + SC_ABORTF ("Inconsistency in vertex info.\n"); + } + } +#endif /*T8_ENABLE_DEBUG */ + + /* Check if this vertex is a boundary node. */ + if (first_tree_entity_dim < mesh_dimension && first_tree_entity_dim >= 0) { + + /* Get the pointer to the array of (u,v)-parameters for the CAD geometry. */ + const double *uv_attribute = (const double *) t8_cmesh_get_attribute ( + associated_cmesh, t8_get_package_id (), T8_CMESH_NODE_PARAMETERS_ATTRIBUTE_KEY, first_tree_id); + + /* 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."); + } + /* Get the (u,v)-parameter of the vertex. */ + const double *uv_parameter = &uv_attribute[2 * local_corner_index]; + + /* Get the pointer to the coordinate array as it was before the deformation. */ + const double *old_coords = (const double *) t8_cmesh_get_attribute ( + associated_cmesh, t8_get_package_id (), T8_CMESH_VERTICES_ATTRIBUTE_KEY, first_tree_id); + + /* Check if the coordinates are available. */ + if (old_coords == nullptr) { + t8_errorf ("Error: Coordinates attribute missing for tree %d\n.", first_tree_id); + SC_ABORTF ("Vertex coordinates are missing."); + } + + gp_Pnt new_coords; + + /* Find the new coordinates of the vertex in the cad file, based on the geometry its lying on. */ + switch (first_tree_entity_dim) { + case 0: { + new_coords = cad->get_cad_point (first_tree_entity_tag); + break; + } + case 1: { + Handle_Geom_Curve curve = cad->get_cad_curve (first_tree_entity_tag); + curve->D0 (uv_parameter[0], new_coords); + break; + } + case 2: { + Handle_Geom_Surface surface = cad->get_cad_surface (first_tree_entity_tag); + surface->D0 (uv_parameter[0], uv_parameter[1], new_coords); + break; + } + default: + SC_ABORT_NOT_REACHED (); + } + + /* Get the old coordinates before the deformation. */ + const double old_x = old_coords[3 * local_corner_index + 0]; + const double old_y = old_coords[3 * local_corner_index + 1]; + const double old_z = old_coords[3 * local_corner_index + 2]; + + /* 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 }; + } + } + return displacements; +} + +void +t8_cmesh_mesh_deformation::apply_vertex_displacements (const std::unordered_map &displacements, + std::shared_ptr cad) +{ + T8_ASSERT (t8_cmesh_is_committed (associated_cmesh)); + + /* Iterate over all vertices in the displacement map. */ + for (const auto &[global_vertex, displacement] : displacements) { + + /* Get the list of trees where this vertex exists. */ + const auto &tree_list = associated_cmesh->vertex_connectivity->get_tree_list_of_vertex (global_vertex); + + /*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. */ + if (tree_vertex_coords != nullptr) { + /* 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]; + } + } + } + } + + /* Update the cad geometry. */ + t8_geometry_handler *geometry_handler = associated_cmesh->geometry_handler; + T8_ASSERT (geometry_handler != nullptr); + + for (auto geom = geometry_handler->begin (); geom != geometry_handler->end (); ++geom) { + if (geom->second->t8_geom_get_type () == T8_GEOMETRY_TYPE_CAD) { + t8_geometry_cad *cad_geom = static_cast (geom->second.get ()); + cad_geom->update_cad_handle (cad); + break; + } + } +} 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 new file mode 100644 index 0000000000..229f2c4fee --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.hxx @@ -0,0 +1,75 @@ +/* + 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_cmesh_mesh_deformation.hxx + * Implementation of CAD-based mesh deformation. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +/** Struct for mesh deformation. */ +struct t8_cmesh_mesh_deformation +{ + public: + /** Constructor */ + t8_cmesh_mesh_deformation (t8_cmesh_t cmesh): associated_cmesh (cmesh), updated_geometry (nullptr) {}; + + /** Destructor */ + ~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 + */ + std::unordered_map + calculate_displacement_surface_vertices (const t8_cad_handle *cad); + + /** + * Apply vertex displacements to a committed cmesh. + * + * 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] cad The shared pointer to the CAD geometry to update. + */ + void + apply_vertex_displacements (const std::unordered_map &displacements, + std::shared_ptr cad); + + private: + /** A pointer to the cmesh for attribute retrieval */ + t8_cmesh_t associated_cmesh; + /** A shared pointer to the updated geometry which comes from a new cad file */ + std::shared_ptr updated_geometry; +}; diff --git a/src/t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_connectivity.hxx b/src/t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_connectivity.hxx index af53065c58..8ed2ff6c50 100644 --- a/src/t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_connectivity.hxx +++ b/src/t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_connectivity.hxx @@ -222,6 +222,27 @@ struct t8_cmesh_vertex_connectivity return get_tree_list_of_vertex (global_vertex_id).size (); } + /** Typedef for the iterator type. */ + using const_iterator = t8_cmesh_vertex_conn_vertex_to_tree::const_iterator; + + /** Iterator begin. + * \return const iterator pointing to the first element. + */ + inline const_iterator + begin () const + { + return vertex_to_tree.begin (); + } + + /** Iterator end. + * \return const iterator pointing behind the last element. + */ + inline const_iterator + end () const + { + return vertex_to_tree.end (); + } + private: /** The internal state. Indicating whether this structure is new and unfilled, the ttv was filled or the vertex conn was built completely. */ state current_state; diff --git a/src/t8_geometry/t8_geometry_handler.hxx b/src/t8_geometry/t8_geometry_handler.hxx index b5af9c63fe..23b83631f4 100644 --- a/src/t8_geometry/t8_geometry_handler.hxx +++ b/src/t8_geometry/t8_geometry_handler.hxx @@ -273,6 +273,26 @@ struct t8_geometry_handler } } + /** + * Return an iterator to the beginning of the geometry handler. + * \return An iterator to the first geometry. + */ + auto + begin () + { + return registered_geometries.begin (); + } + + /** + * Return an iterator to the end of the geometry handler. + * \return An iterator to the element following the last geometry. + */ + auto + end () + { + return registered_geometries.end (); + } + private: /** * Add a geometry to the geometry handler. 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 8e4d15d96a..844ce8b6fa 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx @@ -142,9 +142,9 @@ struct t8_geometry_cad: public t8_geometry_with_vertices } /** - * Getter function for the CAD manager. + * Getter function for the CAD handle. * - * \return The CAD manager of the geometry. + * \return The CAD handle of the geometry. */ std::shared_ptr get_cad_handle () const @@ -152,6 +152,15 @@ struct t8_geometry_cad: public t8_geometry_with_vertices return cad_handle; } + /** Update the CAD handle with a new one. + * \param[in] new_cad_handle The new CAD handle to be used. + */ + void + update_cad_handle (std::shared_ptr new_cad_handle) + { + cad_handle = new_cad_handle; + } + private: /** * Maps points in the reference space \f$ [0,1]^2 \f$ to \f$ \mathbb{R}^3 \f$. Only for triangle trees. From 08848677b29c6348fadc5e5475367335479b992b Mon Sep 17 00:00:00 2001 From: Lena Radmer Date: Mon, 30 Mar 2026 08:53:54 +0200 Subject: [PATCH 03/15] first steps for RBF implementation --- .../t8_cmesh_mesh_deformation.cxx | 24 ++-- .../t8_cmesh_mesh_deformation.hxx | 10 +- .../t8_cmesh_mesh_deformation_rbf.cxx | 28 +++++ .../t8_cmesh_mesh_deformation_rbf.hxx | 115 ++++++++++++++++++ 4 files changed, 165 insertions(+), 12 deletions(-) create mode 100644 src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.cxx create mode 100644 src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.hxx 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 21e593904a..20a66a7f48 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 @@ -33,7 +33,7 @@ #include #include -std::unordered_map +std::unordered_map t8_cmesh_mesh_deformation::calculate_displacement_surface_vertices (const t8_cad_handle *cad) { T8_ASSERT (t8_cmesh_is_committed (associated_cmesh)); @@ -41,7 +41,7 @@ t8_cmesh_mesh_deformation::calculate_displacement_surface_vertices (const t8_cad 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)) { @@ -139,21 +139,29 @@ 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 = { 0.0, 0.0, 0.0 }; + + 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 ( + const std::unordered_map &boundary_node_data, std::shared_ptr cad) { T8_ASSERT (t8_cmesh_is_committed (associated_cmesh)); /* Iterate over all vertices in the displacement map. */ - for (const auto &[global_vertex, displacement] : displacements) { + for (const auto &[global_vertex, rbf_boundary_node] : boundary_node_data) { /* Get the list of trees where this vertex exists. */ const auto &tree_list = associated_cmesh->vertex_connectivity->get_tree_list_of_vertex (global_vertex); @@ -169,7 +177,7 @@ t8_cmesh_mesh_deformation::apply_vertex_displacements (const std::unordered_map< if (tree_vertex_coords != nullptr) { /* 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] += rbf_boundary_node.displacement[coord_index]; } } } 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..2c684bc08f 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,9 +51,9 @@ 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 + * \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 + std::unordered_map calculate_displacement_surface_vertices (const t8_cad_handle *cad); /** @@ -60,11 +62,11 @@ 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. */ void - apply_vertex_displacements (const std::unordered_map &displacements, + apply_vertex_displacements (const std::unordered_map &boundary_node_data, std::shared_ptr cad); private: diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.cxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.cxx new file mode 100644 index 0000000000..e8d0c59436 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.cxx @@ -0,0 +1,28 @@ +/* + 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_cmesh_mesh_deformation_rbf.cxx + * This file implements the Radial Basis Functions for the mesh deformation. + */ + +#include +#include diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.hxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.hxx new file mode 100644 index 0000000000..8b96d04ae0 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.hxx @@ -0,0 +1,115 @@ +/* + 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_cmesh_mesh_deformation_rbf.hxx + * Implementation of CAD-based mesh deformation. + */ +#pragma once + +#include +#include +#include + +/** + * The available RBF function types. + */ +typedef enum { T8_RBF_CP_C2 = 0, T8_RBF_TPS, T8_RBF_FUNCTION_COUNT } t8_rbf_function_type; + +/** + * 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 +{ + /* The old position of the boundary node. */ + t8_3D_vec position; + /* The displacement known and calculated from the different coordinates of the CAD geometries. */ + t8_3D_vec displacement; + /* The calculated RBF coefficient alpha. */ + t8_3D_vec weight; +}; + +/** + * Struct for mesh deformation using Radial Basis Functions (RBF). + * It handles the interpolation of the boundary displacements to the inner nodes. + */ +struct t8_cmesh_mesh_deformation_rbf +{ + public: + /** Constructor. + * \param[in] support_radius The radius of the compact local support for the Wendland CP C2 function. + * If the TPS RBF is used, we will not use the support_radius due to its global support. + * \param[in] rbf_type The RBF type to be used. + */ + t8_cmesh_mesh_deformation_rbf (double support_radius, t8_rbf_function_type rbf_type) + : radius (support_radius), selected_type (rbf_type) {}; + + /** Destructor. */ + ~t8_cmesh_mesh_deformation_rbf () {}; + + /** + * Add a new support node to the RBF system based on the boundary node (they are the same if we do not use a greedy algorithm to get less support nodes). + * These nodes are then used in the linear system A * alpha = d. + * The coefficient alpha is the needed weight of a support node which ensures that the RBF interpolation + * exactly matches the CAD displacement because it corrects the spatial overlap of nearby basis functions. + * \param[in] position The coordinates of the node before the displement. + * \param[in] displacement The known displacement from the input geometries. + */ + void + add_node (const double position[3], const double displacement[3]); + + /** + * 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] inner_node + * \param[out] inner_node_displacement + */ + void + interpolate (const double inner_node[3], const double inner_node_displacement[3]); + + private: + /** Wendland function. psi(x) = (1 - x)_+^4 * (4x + 1) for (1-x) > 0. + * \param[in] distance The euclidean distance between two points. + * return The computed weight (when using solve()) or influence (when using interpolate()). + */ + const double + wendland_cp_c2 (const double distance); + + /** Thin Plate Spline: psi(x) = x^2 * log(x). + * \param[in] distance The euclidean distance between two points. + * return The computed weight (when using solve()) or influence (when using interpolate()). + */ + const double + thin_plate_spline (const double distance); + /** The chosen radius for the local support of the RBF CP C^2. */ + double radius; + /** List of all registered support nodes. */ + std::vector boundary_nodes; + /** Selected Radial Basis Function. */ + t8_rbf_function_type selected_type; +}; From 17f7cbdf25cd7c3fb1814ed44bfe0758da929bfe Mon Sep 17 00:00:00 2001 From: Lena Radmer Date: Wed, 8 Apr 2026 12:50:06 +0200 Subject: [PATCH 04/15] add basis RBF function class and the CPC^2 and TPS functions. add also the rbf_node struct and the interpolate function. add some minor changes to the mesh deformation to include the rbf_node struct to save the node data. --- .../t8_cmesh_mesh_deformation.cxx | 4 +- .../t8_cmesh_mesh_deformation.hxx | 2 +- .../t8_cmesh_mesh_deformation_rbf.hxx | 115 --------- ...sh_mesh_deformation_rbf.cxx => t8_rbf.cxx} | 2 +- .../t8_cmesh_mesh_deformation/t8_rbf.hxx | 236 ++++++++++++++++++ 5 files changed, 240 insertions(+), 119 deletions(-) delete mode 100644 src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.hxx rename src/t8_cmesh/t8_cmesh_mesh_deformation/{t8_cmesh_mesh_deformation_rbf.cxx => t8_rbf.cxx} (92%) create mode 100644 src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx 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 20a66a7f48..91151af1ef 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 @@ -146,7 +146,7 @@ t8_cmesh_mesh_deformation::calculate_displacement_surface_vertices (const t8_cad /* Calculate the displacement of the vertex which should be then done in the deformation. */ node.displacement = { new_coords.X () - old_x, new_coords.Y () - old_y, new_coords.Z () - old_z }; - node.weight = { 0.0, 0.0, 0.0 }; + node.weight.fill (0.0); boundary_node_data[global_vertex_id] = node; } @@ -166,7 +166,7 @@ t8_cmesh_mesh_deformation::apply_vertex_displacements ( /* Get the list of trees where this vertex exists. */ const auto &tree_list = associated_cmesh->vertex_connectivity->get_tree_list_of_vertex (global_vertex); - /*Update the vertex coordinates in each tree. */ + /* 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. */ 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 2c684bc08f..8e181db688 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,7 +35,7 @@ #include #include -#include +#include /** Struct for mesh deformation. */ struct t8_cmesh_mesh_deformation diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.hxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.hxx deleted file mode 100644 index 8b96d04ae0..0000000000 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.hxx +++ /dev/null @@ -1,115 +0,0 @@ -/* - 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_cmesh_mesh_deformation_rbf.hxx - * Implementation of CAD-based mesh deformation. - */ -#pragma once - -#include -#include -#include - -/** - * The available RBF function types. - */ -typedef enum { T8_RBF_CP_C2 = 0, T8_RBF_TPS, T8_RBF_FUNCTION_COUNT } t8_rbf_function_type; - -/** - * 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 -{ - /* The old position of the boundary node. */ - t8_3D_vec position; - /* The displacement known and calculated from the different coordinates of the CAD geometries. */ - t8_3D_vec displacement; - /* The calculated RBF coefficient alpha. */ - t8_3D_vec weight; -}; - -/** - * Struct for mesh deformation using Radial Basis Functions (RBF). - * It handles the interpolation of the boundary displacements to the inner nodes. - */ -struct t8_cmesh_mesh_deformation_rbf -{ - public: - /** Constructor. - * \param[in] support_radius The radius of the compact local support for the Wendland CP C2 function. - * If the TPS RBF is used, we will not use the support_radius due to its global support. - * \param[in] rbf_type The RBF type to be used. - */ - t8_cmesh_mesh_deformation_rbf (double support_radius, t8_rbf_function_type rbf_type) - : radius (support_radius), selected_type (rbf_type) {}; - - /** Destructor. */ - ~t8_cmesh_mesh_deformation_rbf () {}; - - /** - * Add a new support node to the RBF system based on the boundary node (they are the same if we do not use a greedy algorithm to get less support nodes). - * These nodes are then used in the linear system A * alpha = d. - * The coefficient alpha is the needed weight of a support node which ensures that the RBF interpolation - * exactly matches the CAD displacement because it corrects the spatial overlap of nearby basis functions. - * \param[in] position The coordinates of the node before the displement. - * \param[in] displacement The known displacement from the input geometries. - */ - void - add_node (const double position[3], const double displacement[3]); - - /** - * 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] inner_node - * \param[out] inner_node_displacement - */ - void - interpolate (const double inner_node[3], const double inner_node_displacement[3]); - - private: - /** Wendland function. psi(x) = (1 - x)_+^4 * (4x + 1) for (1-x) > 0. - * \param[in] distance The euclidean distance between two points. - * return The computed weight (when using solve()) or influence (when using interpolate()). - */ - const double - wendland_cp_c2 (const double distance); - - /** Thin Plate Spline: psi(x) = x^2 * log(x). - * \param[in] distance The euclidean distance between two points. - * return The computed weight (when using solve()) or influence (when using interpolate()). - */ - const double - thin_plate_spline (const double distance); - /** The chosen radius for the local support of the RBF CP C^2. */ - double radius; - /** List of all registered support nodes. */ - std::vector boundary_nodes; - /** Selected Radial Basis Function. */ - t8_rbf_function_type selected_type; -}; diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.cxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx similarity index 92% rename from src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.cxx rename to src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx index e8d0c59436..7c5b8ebcbb 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation_rbf.cxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx @@ -24,5 +24,5 @@ * This file implements the Radial Basis Functions for the mesh deformation. */ -#include +#include #include 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..f88151f171 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx @@ -0,0 +1,236 @@ +/* + 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 + +/** + * The available RBF function types. + */ +typedef enum { T8_RBF_CP_C2 = 0, T8_RBF_TPS, T8_RBF_FUNCTION_COUNT } t8_rbf_function_type; + +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. */ + t8_3D_vec weight; +}; + +struct t8_rbf_function +{ + /** Destructor. */ + virtual ~t8_rbf_function () {}; + virtual double + evaluate (double distance) const + = 0; + 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. + * \param[in] r The support radius. Beyond the radius the function is zero and so there is no impact from this node movement. + */ + t8_rbf_cpc2 (double r): radius (r) + { + } + /** 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. + * \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) 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; + } + + private: + /** The chosen support radius. */ + double radius; +}; + +/** + * 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. + * \return The function value psi. + */ + double + evaluate (double distance) 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] support_radius The radius of the compact local support for the Wendland CP C2 function. + * If the TPS RBF is used, we will not use the support_radius due to its global support. + * \param[in] rbf_type The RBF type to be used. + */ + t8_rbf (double support_radius, t8_rbf_function_type rbf_type) + { + + if (rbf_type == T8_RBF_CP_C2) { + rbf_function = std::make_unique (support_radius); + } + else if (rbf_type == T8_RBF_TPS) { + rbf_function = std::make_unique (); + } + else { + t8_errorf ("ERROR: RBF attribute missing or not correct\n."); + SC_ABORTF ("Unsupported RBF type."); + } + }; + + /** Destructor. */ + ~t8_rbf () {}; + + 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_nodes.push_back (std::move (boundary_node)); + } + } + + /** + * 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); + /** Check if the basis function value (the infleunce 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) { + inner_node.displacement[coordinate] += boundary_node.weight[coordinate] * psi; + } + } + } + } + + private: + /** List of all registered support nodes. */ + std::vector boundary_nodes; + /** */ + std::unique_ptr rbf_function; +}; From cf17193d85bf63d21418f45df942536d47b50ba5 Mon Sep 17 00:00:00 2001 From: Lena Radmer Date: Thu, 23 Apr 2026 17:24:36 +0200 Subject: [PATCH 05/15] add more functionalities to the RBF class and implement the solver with eigen next --- CMakeLists.txt | 14 ++++- src/CMakeLists.txt | 12 +++++ src/config.cmake.in | 5 ++ .../t8_cmesh_mesh_deformation/t8_rbf.cxx | 4 ++ .../t8_cmesh_mesh_deformation/t8_rbf.hxx | 52 ++++++++++++++++++- 5 files changed, 85 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 66945282af..7fcfab5ef8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -41,6 +41,14 @@ FetchContent_Declare( GIT_PROGRESS TRUE ) +FetchContent_Declare( + eigen + GIT_REPOSITORY "https://gitlab.com/libeigen/eigen.git" + GIT_TAG 5.0.1 + GIT_PROGRESS TRUE + GIT_SHALLOW TRUE + ) + set(gtest_force_shared_crt ON CACHE BOOL "" FORCE) mark_as_advanced( FORCE gtest_force_shared_crt) @@ -66,6 +74,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 ) @@ -270,7 +279,7 @@ else() mark_as_advanced( FORCE ${_new_p4est_vars} ) endif() -# Workaround: Suppress warnings for googletests, so it does not compromise the usage of -WError for t8code. +# Workaround: Suppress warnings for fetchcontent libraries, so it does not compromise the usage of -WError for t8code. # ----------- # 1. Save original flags set(ORIGINAL_CXX_FLAGS "${CMAKE_CXX_FLAGS}") @@ -281,6 +290,9 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -w" CACHE STRING "" FORCE) set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -w" CACHE STRING "" FORCE) FetchContent_MakeAvailable( googletest ) +if ( T8CODE_ENABLE_EIGEN ) + FetchContent_MakeAvailable( eigen ) +endif() # 3. Restore original flags set(CMAKE_CXX_FLAGS "${ORIGINAL_CXX_FLAGS}" CACHE STRING "" FORCE) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8f65b2a588..040c37f822 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -114,6 +114,18 @@ if( T8CODE_ENABLE_OCC ) ) 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_BUILD_PEDANTIC ) target_compile_options( T8 PUBLIC -pedantic ) set (T8_CXXFLAGS "${T8_CXXFLAGS} -Wpedantic") diff --git a/src/config.cmake.in b/src/config.cmake.in index c8fa76969b..ab1a0021da 100644 --- a/src/config.cmake.in +++ b/src/config.cmake.in @@ -25,6 +25,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_rbf.cxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx index 7c5b8ebcbb..13b02de590 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx @@ -26,3 +26,7 @@ #include #include + +#include + +/** Hier solve rein. */ diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx index f88151f171..10356aaeb2 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx @@ -93,7 +93,9 @@ struct t8_rbf_cpc2: public t8_rbf_function double evaluate (double distance) const override { - double r = distance / radius; + double r + = distance + / radius; //guten radius finden und zeit darein verschwenden ;) inspo von molekulardynamik, schranken erarbeiten und shcon schreiben, if (r < 1.0) { double d = 1.0 - r; return (d * d * d * d) * (4.0 * r + 1.0); @@ -185,12 +187,18 @@ struct t8_rbf /** Destructor. */ ~t8_rbf () {}; + /** + * Transfers the boundary node data to the internal data structure of the RBF. + */ 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)); } } @@ -202,6 +210,48 @@ struct t8_rbf void 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.\n."); + SC_ABORTF ("The current RBF instance has no boundary nodes."); + } + + if (rbf_function->is_compactly_supported ()) { + /** If the RBF used is compactly supported use a sparse matrix. Then the conjugate gradient method can be used for solving the linear equation system. */ + } + else { + /** If the RBF used is globally supported use a dense matrix. */ + } + + /** Create the matrix A for the linear system. */ + std::vector 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 = row; 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); + /** Write the value to the matrix A. */ + A[row * num_boundary_nodes + col] = psi; + /** Mirror the value to the lower triangular part of the matrix if it's not on the diagonal. */ + if (col != row) { + A[col * num_boundary_nodes + row] = psi; + } + } + } + + /** Create the right-hand side vector for the linear system. It contains the displacements of the boundary nodes. */ + std::vector displacements (num_boundary_nodes * 3); + for (size_t i = 0; i < num_boundary_nodes; ++i) { + displacements[3 * i + 0] = boundary_nodes[i].displacement[0]; + displacements[3 * i + 1] = boundary_nodes[i].displacement[1]; + displacements[3 * i + 2] = boundary_nodes[i].displacement[2]; + } } /** From 25422836ca10ba89c134df3aac1c186e22a2d826 Mon Sep 17 00:00:00 2001 From: Lena-Richter Date: Mon, 4 May 2026 07:44:03 +0200 Subject: [PATCH 06/15] include eigen library and add a little test to try if it works --- example/CMakeLists.txt | 6 ++ example/cmesh/t8_rbf.cxx | 85 +++++++++++++++++ .../t8_cmesh_mesh_deformation/t8_rbf.cxx | 94 ++++++++++++++++++- .../t8_cmesh_mesh_deformation/t8_rbf.hxx | 46 +-------- 4 files changed, 184 insertions(+), 47 deletions(-) create mode 100644 example/cmesh/t8_rbf.cxx diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index a897b4d314..048a0d818c 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -93,6 +93,12 @@ 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 ) +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_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/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx index 13b02de590..607f3a8f23 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx @@ -24,9 +24,99 @@ * This file implements the Radial Basis Functions for the mesh deformation. */ +#include #include #include - +#include #include +#include + +typedef Eigen::SparseMatrix SpMat; // declares a column-major sparse matrix type of double + +/** + * 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 () +{ + double estimation_of_entries + = 50; // brauchen noch gute Schätzung für die Anzahl der nicht-Null Element in der Matrix A für compactly supported RBFs + 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.\n."); + SC_ABORTF ("The current RBF instance has no boundary nodes."); + } + + /** 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]); + } + /** Allocate the memory for the weight Vector. */ + Eigen::MatrixXd alpha (num_boundary_nodes, 3); + + if (rbf_function->is_compactly_supported ()) { + /** If the RBF used is compactly supported use a sparse matrix. Then the conjugate gradient method can be used for solving the linear equation system. */ + 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) { + /** 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 = row; 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); + coefficients.emplace_back (row, col, psi); + if (row != col) { + coefficients.emplace_back (col, row, psi); + } + } + } + Eigen::SparseMatrix A (num_boundary_nodes, num_boundary_nodes); + A.setFromTriplets (coefficients.begin (), coefficients.end ()); + Eigen::ConjugateGradient solver; + solver.compute (A); + if (solver.info () != Eigen::Success) { + t8_errorf ("ERROR: Decomposition of the matrix A failed.\n."); + SC_ABORTF ("The linear system cannot be solved."); + } + } + else { + /** If the RBF used is globally supported use a dense matrix. */ + } -/** Hier solve rein. */ + /** Create the matrix A for the linear system. */ + std::vector 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 = row; 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); + /** Write the value to the matrix A. */ + A[row * num_boundary_nodes + col] = psi; + /** Mirror the value to the lower triangular part of the matrix if it's not on the diagonal. */ + if (col != row) { + A[col * num_boundary_nodes + row] = psi; + } + } + } +#if 0 + /** Create the right-hand side vector for the linear system. It contains the displacements of the boundary nodes. */ + std::vector displacements (num_boundary_nodes * 3); + for (size_t i = 0; i < num_boundary_nodes; ++i) { + displacements[3 * i + 0] = boundary_nodes[i].displacement[0]; + displacements[3 * i + 1] = boundary_nodes[i].displacement[1]; + displacements[3 * i + 2] = boundary_nodes[i].displacement[2]; + } +#endif +} diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx index 10356aaeb2..ec09a06f3b 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx @@ -208,51 +208,7 @@ struct t8_rbf * This step is mandatory to later be able to interpolate the inner nodes. */ void - 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.\n."); - SC_ABORTF ("The current RBF instance has no boundary nodes."); - } - - if (rbf_function->is_compactly_supported ()) { - /** If the RBF used is compactly supported use a sparse matrix. Then the conjugate gradient method can be used for solving the linear equation system. */ - } - else { - /** If the RBF used is globally supported use a dense matrix. */ - } - - /** Create the matrix A for the linear system. */ - std::vector 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 = row; 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); - /** Write the value to the matrix A. */ - A[row * num_boundary_nodes + col] = psi; - /** Mirror the value to the lower triangular part of the matrix if it's not on the diagonal. */ - if (col != row) { - A[col * num_boundary_nodes + row] = psi; - } - } - } - - /** Create the right-hand side vector for the linear system. It contains the displacements of the boundary nodes. */ - std::vector displacements (num_boundary_nodes * 3); - for (size_t i = 0; i < num_boundary_nodes; ++i) { - displacements[3 * i + 0] = boundary_nodes[i].displacement[0]; - displacements[3 * i + 1] = boundary_nodes[i].displacement[1]; - displacements[3 * i + 2] = boundary_nodes[i].displacement[2]; - } - } + solve (); /** * Interpolates the inner node. From 19f17993974d3859f28863b07232e34f0ece99e4 Mon Sep 17 00:00:00 2001 From: Lena Richter Date: Sat, 9 May 2026 15:38:09 +0200 Subject: [PATCH 07/15] add solve function --- .../t8_cmesh_mesh_deformation/t8_rbf.cxx | 103 +++++++++++------- .../t8_cmesh_mesh_deformation/t8_rbf.hxx | 17 ++- 2 files changed, 77 insertions(+), 43 deletions(-) diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx index 607f3a8f23..639e0a92d9 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx @@ -31,8 +31,6 @@ #include #include -typedef Eigen::SparseMatrix SpMat; // declares a column-major sparse matrix type of double - /** * Solves the linear system A * alpha = displacements to find the weight. * This step is mandatory to later be able to interpolate the inner nodes. @@ -40,13 +38,10 @@ typedef Eigen::SparseMatrix SpMat; // declares a column-major sparse ma void t8_rbf::solve () { - double estimation_of_entries - = 50; // brauchen noch gute Schätzung für die Anzahl der nicht-Null Element in der Matrix A für compactly supported RBFs 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.\n."); - SC_ABORTF ("The current RBF instance has no boundary nodes."); + 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. */ @@ -55,44 +50,72 @@ t8_rbf::solve () displacements.row (i) = Eigen::Vector3d (boundary_nodes[i].displacement[0], boundary_nodes[i].displacement[1], boundary_nodes[i].displacement[2]); } - /** Allocate the memory for the weight Vector. */ - Eigen::MatrixXd alpha (num_boundary_nodes, 3); + /** Weight vector. */ + Eigen::MatrixXd alpha; if (rbf_function->is_compactly_supported ()) { - /** If the RBF used is compactly supported use a sparse matrix. Then the conjugate gradient method can be used for solving the linear equation system. */ - 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) { - /** Because of the symmetric property of the distance between nodes, we only need to compute the upper triangular part of the matrix + 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. */ + double estimation_of_entries + = 50; // brauchen noch gute Schätzung für die Anzahl der nicht-Null Element in der Matrix A für compactly supported RBFs + 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) { + /** 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 = row; 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); + for (size_t col = row; 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); + /** 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); if (row != col) { coefficients.emplace_back (col, row, psi); } } } - Eigen::SparseMatrix A (num_boundary_nodes, num_boundary_nodes); - A.setFromTriplets (coefficients.begin (), coefficients.end ()); - Eigen::ConjugateGradient solver; - solver.compute (A); - if (solver.info () != Eigen::Success) { - t8_errorf ("ERROR: Decomposition of the matrix A failed.\n."); - SC_ABORTF ("The linear system cannot be solved."); - } } - else { - /** If the RBF used is globally supported use a dense matrix. */ + /** 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::ConjugateGradient, Eigen::Lower | Eigen::Upper> solver; + 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. */ - std::vector A (num_boundary_nodes * num_boundary_nodes); + 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 @@ -103,20 +126,18 @@ t8_rbf::solve () /** Evaluate the radial basis function for the current distance. */ const double psi = rbf_function->evaluate (distance); /** Write the value to the matrix A. */ - A[row * num_boundary_nodes + col] = psi; + A (row, col) = psi; /** Mirror the value to the lower triangular part of the matrix if it's not on the diagonal. */ if (col != row) { - A[col * num_boundary_nodes + row] = psi; + A (col, row) = psi; } } } -#if 0 - /** Create the right-hand side vector for the linear system. It contains the displacements of the boundary nodes. */ - std::vector displacements (num_boundary_nodes * 3); - for (size_t i = 0; i < num_boundary_nodes; ++i) { - displacements[3 * i + 0] = boundary_nodes[i].displacement[0]; - displacements[3 * i + 1] = boundary_nodes[i].displacement[1]; - displacements[3 * i + 2] = boundary_nodes[i].displacement[2]; + Eigen::LDLT solver (A); + + if (solver.info () != Eigen::Success) { + t8_errorf ("ERROR: The global RBF system could not be solved with LDLT.\n"); } -#endif + + 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 index ec09a06f3b..7755d22eb0 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx @@ -28,7 +28,10 @@ #include #include #include +#include #include +#include +#include /** * The available RBF function types. @@ -224,7 +227,7 @@ struct t8_rbf double distance = t8_dist (inner_node.position, boundary_node.position); const double psi = rbf_function->evaluate (distance); - /** Check if the basis function value (the infleunce factor) is not equal to zero with a numerical tolerance of 1e-12. */ + /** 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) { @@ -235,8 +238,18 @@ struct t8_rbf } 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; }; From 64df1b1e53264fbda1445df1067bf6050b92670e Mon Sep 17 00:00:00 2001 From: Lena Richter Date: Fri, 19 Jun 2026 14:51:48 +0200 Subject: [PATCH 08/15] add eigen as fetch content --- cmake/thirdparty.cmake | 5 +++-- cmake/thirdparty.json | 13 +++++++++++++ 2 files changed, 16 insertions(+), 2 deletions(-) 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 b7de6a4148..91e7c43084 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 From b0fcde7addf7cd6e9e0105aa43f64c72556a5de1 Mon Sep 17 00:00:00 2001 From: Lena Richter Date: Fri, 19 Jun 2026 15:13:53 +0200 Subject: [PATCH 09/15] first deformation with inner nodes --- example/CMakeLists.txt | 7 +- example/cmesh/t8_cmesh_mesh_deformation.cxx | 17 +- .../t8_cmesh_mesh_deformation.cxx | 166 ++++++++++++++++-- .../t8_cmesh_mesh_deformation.hxx | 7 +- .../t8_cmesh_mesh_deformation/t8_rbf.cxx | 17 +- .../t8_cmesh_mesh_deformation/t8_rbf.hxx | 54 +++--- 6 files changed, 212 insertions(+), 56 deletions(-) diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index 476978a94b..4ef52b0b89 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,9 +88,11 @@ if(T8CODE_ENABLE_VTK) endif() if (T8CODE_ENABLE_EIGEN) - add_t8_example( NAME t8_rbf_test SOURCES cmesh/t8_rbf.cxx ) + 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 ) - target_link_libraries( t8_rbf_test PRIVATE Eigen3::Eigen ) endif () add_t8_example( NAME t8_gmsh_to_vtk SOURCES IO/forest/gmsh/t8_gmsh_to_vtk.cxx ) diff --git a/example/cmesh/t8_cmesh_mesh_deformation.cxx b/example/cmesh/t8_cmesh_mesh_deformation.cxx index 11a8375ec3..156f52c594 100644 --- a/example/cmesh/t8_cmesh_mesh_deformation.cxx +++ b/example/cmesh/t8_cmesh_mesh_deformation.cxx @@ -28,7 +28,7 @@ #include #include #include -#if T8CODE_ENABLE_OCC +#if T8_ENABLE_OCC #include #include #endif /* T8CODE_ENABLE_OCC */ @@ -43,7 +43,7 @@ int main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) { -#if T8CODE_ENABLE_OCC +#if T8_ENABLE_OCC char usage[BUFSIZ]; /* Brief help message. */ @@ -85,6 +85,8 @@ 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]); @@ -94,6 +96,9 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) "File prefix of the deformation geometry file (without .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); @@ -124,14 +129,18 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) /* Initialize the deformation object for the given mesh. */ t8_cmesh_mesh_deformation deformation (cmesh); + /** Save the input RBF type. */ + t8_rbf_function_type rbf_type = static_cast (rbf_type_int); + /* Calculate displacements. */ - auto displacements = deformation.calculate_displacement_surface_vertices (cad.get ()); + auto displacements + = deformation.calculate_displacement_surface_vertices (cad.get (), rbf_type, scale_factor_support_radius); /* Write output. */ t8_forest_vtk_write_file (forest, "input_forest", 1, 1, 1, 1, 0, 0, NULL); /* Apply displacements. */ - deformation.apply_vertex_displacements (displacements, cad); + deformation.apply_vertex_displacements (displacements, cad, rbf_type); /* Write output. */ t8_forest_vtk_write_file (forest, "deformed_forest", 1, 1, 1, 1, 0, 0, NULL); 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 6e2b085d65..cfcb2e7605 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,8 +35,83 @@ #include #include +static double +calculate_local_support_radius (t8_cmesh_t cmesh, const std::vector> &tree_list, + int current_global_vertex_id) +{ + double total_distance = 0.0; + + /** 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); + } + /** To get the average, we divide through the amount of neighbor nodes. */ + return (total_distance / static_cast (unique_neighbors.size ())); +} + std::unordered_map -t8_cmesh_mesh_deformation::calculate_displacement_surface_vertices (const t8_cad_handle *cad) +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)); @@ -115,7 +192,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 its lying on. */ switch (first_tree_entity_dim) { case 0: { new_coords = cad->get_cad_point (first_tree_entity_tag); @@ -149,6 +226,10 @@ t8_cmesh_mesh_deformation::calculate_displacement_surface_vertices (const t8_cad 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); + } boundary_node_data[global_vertex_id] = node; } } @@ -157,34 +238,91 @@ t8_cmesh_mesh_deformation::calculate_displacement_surface_vertices (const t8_cad void t8_cmesh_mesh_deformation::apply_vertex_displacements ( - const std::unordered_map &boundary_node_data, std::shared_ptr cad) + 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, rbf_boundary_node] : boundary_node_data) { + 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 (); - /* Get the vertex coordinates of the current tree. */ + /** 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 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 a 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 a inner node we don not know the displacement right away and need to interpolate to get the new coordinates. */ + else { + + double *vertex_coords = (double *) 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] += rbf_boundary_node.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 8e181db688..f14e56cd70 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 @@ -54,7 +54,8 @@ struct t8_cmesh_mesh_deformation * \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); + 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. @@ -66,8 +67,8 @@ struct t8_cmesh_mesh_deformation * \param [in] cad The shared pointer to the CAD geometry to update. */ void - apply_vertex_displacements (const std::unordered_map &boundary_node_data, - 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 index 639e0a92d9..0844afddbe 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx @@ -71,26 +71,22 @@ 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. */ - double estimation_of_entries + size_t estimation_of_entries = 50; // brauchen noch gute Schätzung für die Anzahl der nicht-Null Element in der Matrix A für compactly supported RBFs 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) { - /** 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 = row; col < num_boundary_nodes; ++col) { + + 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); + 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); - if (row != col) { - coefficients.emplace_back (col, row, psi); - } } } } @@ -98,7 +94,8 @@ t8_rbf::solve_compactly_supported_rbf (const Eigen::MatrixXd &displacements, con 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::ConjugateGradient, Eigen::Lower | Eigen::Upper> solver; + Eigen::BiCGSTAB> solver; + solver.setTolerance (1e-10); solver.compute (A); if (solver.info () != Eigen::Success) { @@ -124,7 +121,7 @@ t8_rbf::solve_globally_supported_rbf (const Eigen::MatrixXd &displacements, cons /** 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); + const double psi = rbf_function->evaluate (distance, 0.0); /** Write the value to the matrix A. */ A (row, col) = psi; /** Mirror the value to the lower triangular part of the matrix if it's not on the diagonal. */ diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx index 7755d22eb0..4e63c3d87a 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx @@ -48,6 +48,8 @@ struct t8_rbf_node * 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; + /** The local support radius for the node. */ + double local_support_radius; }; /** @@ -56,7 +58,8 @@ struct t8_rbf_node */ 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. */ + /** 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; }; @@ -65,7 +68,7 @@ struct t8_rbf_function /** Destructor. */ virtual ~t8_rbf_function () {}; virtual double - evaluate (double distance) const + evaluate (double distance, double radius) const = 0; virtual bool is_compactly_supported () const @@ -78,10 +81,8 @@ struct t8_rbf_function */ struct t8_rbf_cpc2: public t8_rbf_function { - /** Constructor. - * \param[in] r The support radius. Beyond the radius the function is zero and so there is no impact from this node movement. - */ - t8_rbf_cpc2 (double r): radius (r) + /** Constructor. */ + t8_rbf_cpc2 () { } /** Destructor. */ @@ -94,11 +95,9 @@ struct t8_rbf_cpc2: public t8_rbf_function * \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) const override + evaluate (double distance, double radius) const override { - double r - = distance - / radius; //guten radius finden und zeit darein verschwenden ;) inspo von molekulardynamik, schranken erarbeiten und shcon schreiben, + double r = distance / radius; if (r < 1.0) { double d = 1.0 - r; return (d * d * d * d) * (4.0 * r + 1.0); @@ -115,10 +114,6 @@ struct t8_rbf_cpc2: public t8_rbf_function { return true; } - - private: - /** The chosen support radius. */ - double radius; }; /** @@ -141,7 +136,7 @@ struct t8_rbf_tps: public t8_rbf_function * \return The function value psi. */ double - evaluate (double distance) const override + evaluate (double distance, [[maybe_unused]] double radius) const override { if (distance < 1e-12) { return 0.0; @@ -168,22 +163,19 @@ struct t8_rbf { public: /** Constructor. - * \param[in] support_radius The radius of the compact local support for the Wendland CP C2 function. - * If the TPS RBF is used, we will not use the support_radius due to its global support. * \param[in] rbf_type The RBF type to be used. */ - t8_rbf (double support_radius, t8_rbf_function_type rbf_type) + t8_rbf (t8_rbf_function_type rbf_type) { if (rbf_type == T8_RBF_CP_C2) { - rbf_function = std::make_unique (support_radius); + rbf_function = std::make_unique (); } else if (rbf_type == T8_RBF_TPS) { rbf_function = std::make_unique (); } else { - t8_errorf ("ERROR: RBF attribute missing or not correct\n."); - SC_ABORTF ("Unsupported RBF type."); + SC_ABORTF ("ERROR: RBF attribute missing or not correct. Unsupported RBF type.\n"); } }; @@ -206,6 +198,22 @@ struct t8_rbf } } + /** + * Search for the displacement of a specific boundary node with its associated global ID. + */ + 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. @@ -226,11 +234,13 @@ struct t8_rbf for (const auto& boundary_node : boundary_nodes) { double distance = t8_dist (inner_node.position, boundary_node.position); - const double psi = rbf_function->evaluate (distance); + + 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; } } From f2d1315ed7c77b0744b027ef5dbd120de6db1c99 Mon Sep 17 00:00:00 2001 From: Lena Richter Date: Sat, 27 Jun 2026 22:50:45 +0200 Subject: [PATCH 10/15] updating the TPS calculation and fix some minor issuses --- example/cmesh/t8_cmesh_mesh_deformation.cxx | 30 ++++++++++++++++--- .../t8_cmesh_mesh_deformation.cxx | 26 +++++++++------- .../t8_cmesh_mesh_deformation/t8_rbf.cxx | 20 ++++--------- .../t8_cmesh_mesh_deformation/t8_rbf.hxx | 4 +-- 4 files changed, 50 insertions(+), 30 deletions(-) diff --git a/example/cmesh/t8_cmesh_mesh_deformation.cxx b/example/cmesh/t8_cmesh_mesh_deformation.cxx index 156f52c594..67db620249 100644 --- a/example/cmesh/t8_cmesh_mesh_deformation.cxx +++ b/example/cmesh/t8_cmesh_mesh_deformation.cxx @@ -98,7 +98,7 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) 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"); + "Scale factor for the support radius. Default: 5"); int parsed = sc_options_parse (t8_get_package_id (), SC_LP_ERROR, opt, argc, argv); @@ -111,7 +111,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 { @@ -132,6 +132,28 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) /** Save the input RBF type. */ t8_rbf_function_type rbf_type = static_cast (rbf_type_int); + /* Write output. */ + t8_forest_vtk_write_file (forest, "deformed_forest_step_0", 1, 1, 1, 1, 0, 0, NULL); + + int num_steps = 50; + for (int num = 1; num <= num_steps; ++num) { + + char brep_buf[256]; + snprintf (brep_buf, sizeof (brep_buf), "%s%d", brep_file, num); + + std::string current_brep (brep_buf); + + auto cad_deformed = std::make_shared (current_brep); + + 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 (num); + t8_forest_vtk_write_file (forest, output_name.c_str (), 1, 1, 1, 1, 0, 0, NULL); + } +#if 0 /* Calculate displacements. */ auto displacements = deformation.calculate_displacement_surface_vertices (cad.get (), rbf_type, scale_factor_support_radius); @@ -144,11 +166,11 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) /* Write output. */ t8_forest_vtk_write_file (forest, "deformed_forest", 1, 1, 1, 1, 0, 0, NULL); - +#endif /* Cleanup. */ t8_forest_unref (&forest); - t8_global_productionf ("Mesh deformation completed."); + t8_global_productionf ("Mesh deformation completed.\n"); } sc_options_destroy (opt); 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 cfcb2e7605..d0df4a7dbe 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 @@ -37,10 +37,12 @@ static double calculate_local_support_radius (t8_cmesh_t cmesh, const std::vector> &tree_list, - int current_global_vertex_id) + 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; @@ -104,8 +106,11 @@ calculate_local_support_radius (t8_cmesh_t cmesh, const std::vector (unique_neighbors.size ())); + + double grading_factor = 2.0; + + return (total_distance / static_cast (unique_neighbors.size ())) + * (1.0 + grading_factor * displacement_magnitude); } std::unordered_map @@ -175,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]; @@ -227,8 +232,9 @@ t8_cmesh_mesh_deformation::calculate_displacement_surface_vertices (const t8_cad 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.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; } @@ -275,7 +281,7 @@ t8_cmesh_mesh_deformation::apply_vertex_displacements ( const int mesh_dimension = t8_cmesh_get_dimension (associated_cmesh); - /** Contains the displacement whether its a boundary node or a inner node. */ + /** 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. @@ -284,11 +290,11 @@ t8_cmesh_mesh_deformation::apply_vertex_displacements ( final_displacement = rbf_handler.get_boundary_displacement (global_vertex_id); } - /** If it is a inner node we don not know the displacement right away and need to interpolate to get the new coordinates. */ + /** If it is an inner node, we do not know the displacement right away and need to interpolate to get the new coordinates. */ else { - double *vertex_coords = (double *) t8_cmesh_get_attribute (associated_cmesh, t8_get_package_id (), - T8_CMESH_VERTICES_ATTRIBUTE_KEY, first_tree.first); + 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); diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx index 0844afddbe..cca4c1abf2 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx @@ -71,8 +71,7 @@ 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; // brauchen noch gute Schätzung für die Anzahl der nicht-Null Element in der Matrix A für compactly supported RBFs + size_t estimation_of_entries = 50; typedef Eigen::Triplet triplet; std::vector coefficients; coefficients.reserve (num_boundary_nodes * estimation_of_entries); @@ -103,7 +102,6 @@ t8_rbf::solve_compactly_supported_rbf (const Eigen::MatrixXd &displacements, con } return solver.solve (displacements); - ; } Eigen::MatrixXd @@ -117,23 +115,17 @@ t8_rbf::solve_globally_supported_rbf (const Eigen::MatrixXd &displacements, cons 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 = row; col < num_boundary_nodes; ++col) { + 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, 0.0); - /** Write the value to the matrix A. */ - A (row, col) = psi; - /** Mirror the value to the lower triangular part of the matrix if it's not on the diagonal. */ - if (col != row) { - A (col, row) = psi; - } + /** 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::LDLT solver (A); + Eigen::PartialPivLU solver (A); if (solver.info () != Eigen::Success) { - t8_errorf ("ERROR: The global RBF system could not be solved with LDLT.\n"); + 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 index 4e63c3d87a..36c262649e 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx @@ -48,8 +48,6 @@ struct t8_rbf_node * 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; - /** The local support radius for the node. */ - double local_support_radius; }; /** @@ -61,6 +59,8 @@ 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; }; struct t8_rbf_function From 984da6c75eb41ce6e132371636eb366612eafa50 Mon Sep 17 00:00:00 2001 From: Lena Richter Date: Mon, 5 Oct 2026 15:01:15 +0200 Subject: [PATCH 11/15] load mesh deformation .brep files in example dynamical --- example/cmesh/t8_cmesh_mesh_deformation.cxx | 52 +++++++++++++++++---- 1 file changed, 44 insertions(+), 8 deletions(-) diff --git a/example/cmesh/t8_cmesh_mesh_deformation.cxx b/example/cmesh/t8_cmesh_mesh_deformation.cxx index 67db620249..0c62886a9d 100644 --- a/example/cmesh/t8_cmesh_mesh_deformation.cxx +++ b/example/cmesh/t8_cmesh_mesh_deformation.cxx @@ -31,7 +31,7 @@ #if T8_ENABLE_OCC #include #include -#endif /* T8CODE_ENABLE_OCC */ +#endif /* T8_ENABLE_OCC */ #include #include @@ -39,6 +39,27 @@ #include #include #include +#include +#include + +#if T8_ENABLE_OCC +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 */ int main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) @@ -93,12 +114,12 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) 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 tothe 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: 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); @@ -123,9 +144,6 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) t8_cmesh_t cmesh = t8_cmesh_from_msh_file (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); @@ -135,6 +153,23 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) /* Write output. */ t8_forest_vtk_write_file (forest, "deformed_forest_step_0", 1, 1, 1, 1, 0, 0, NULL); + auto brep_files = findBrepFiles (brep_file); + std::sort (brep_files.begin (), brep_files.end ()); + + int ifile = 0; + for (const auto &file : brep_files) { + auto file_without_ext = file.parent_path () / file.stem (); + auto cad_deformaed = std::make_shared (file_without_ext.c_str ()); + + auto displacements = deformation.calculate_displacement_surface_vertices (cad_deformaed.get (), rbf_type, + scale_factor_support_radius); + + deformation.apply_vertex_displacements (displacements, cad_deformaed, 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); + } +#if 0 int num_steps = 50; for (int num = 1; num <= num_steps; ++num) { @@ -153,6 +188,7 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) std::string output_name = "deformed_forest_step_" + std::to_string (num); t8_forest_vtk_write_file (forest, output_name.c_str (), 1, 1, 1, 1, 0, 0, NULL); } +#endif #if 0 /* Calculate displacements. */ auto displacements @@ -179,9 +215,9 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) mpiret = sc_MPI_Finalize (); SC_CHECK_MPI (mpiret); -#else /* T8CODE_ENABLE_OCC */ +#else /* T8_ENABLE_OCC */ t8_global_errorf ("ERROR: This example requires OpenCASCADE support to be enabled in t8code.\n"); -#endif /* T8CODE_ENABLE_OCC */ +#endif /* T8_ENABLE_OCC */ return 0; } From d52a3483db639a667cb4b2f00988a7a35462d3fe Mon Sep 17 00:00:00 2001 From: Lena Richter Date: Mon, 5 Oct 2026 15:39:31 +0200 Subject: [PATCH 12/15] clean up and fix some typos --- example/cmesh/t8_cmesh_mesh_deformation.cxx | 40 ++------------------- 1 file changed, 3 insertions(+), 37 deletions(-) diff --git a/example/cmesh/t8_cmesh_mesh_deformation.cxx b/example/cmesh/t8_cmesh_mesh_deformation.cxx index 0c62886a9d..3854e4b059 100644 --- a/example/cmesh/t8_cmesh_mesh_deformation.cxx +++ b/example/cmesh/t8_cmesh_mesh_deformation.cxx @@ -114,7 +114,7 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) 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, - "Path tothe folder containing the deformation geometry files (.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"); @@ -159,50 +159,16 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) int ifile = 0; for (const auto &file : brep_files) { auto file_without_ext = file.parent_path () / file.stem (); - auto cad_deformaed = std::make_shared (file_without_ext.c_str ()); - - auto displacements = deformation.calculate_displacement_surface_vertices (cad_deformaed.get (), rbf_type, - scale_factor_support_radius); - - deformation.apply_vertex_displacements (displacements, cad_deformaed, 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); - } -#if 0 - int num_steps = 50; - for (int num = 1; num <= num_steps; ++num) { - - char brep_buf[256]; - snprintf (brep_buf, sizeof (brep_buf), "%s%d", brep_file, num); - - std::string current_brep (brep_buf); - - auto cad_deformed = std::make_shared (current_brep); + 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 (num); + 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); } -#endif -#if 0 - /* Calculate displacements. */ - auto displacements - = deformation.calculate_displacement_surface_vertices (cad.get (), rbf_type, scale_factor_support_radius); - - /* Write output. */ - t8_forest_vtk_write_file (forest, "input_forest", 1, 1, 1, 1, 0, 0, NULL); - - /* Apply displacements. */ - deformation.apply_vertex_displacements (displacements, cad, rbf_type); - - /* Write output. */ - t8_forest_vtk_write_file (forest, "deformed_forest", 1, 1, 1, 1, 0, 0, NULL); -#endif /* Cleanup. */ t8_forest_unref (&forest); From d05a2e2ce20192d1f2db95c7e4a1e97e37bdb1ad Mon Sep 17 00:00:00 2001 From: Lena Richter Date: Tue, 6 Oct 2026 11:49:02 +0200 Subject: [PATCH 13/15] fix doxygen errors and missing parameter documentation --- .../t8_cmesh_mesh_deformation.hxx | 3 ++ .../t8_cmesh_mesh_deformation/t8_rbf.cxx | 2 +- .../t8_cmesh_mesh_deformation/t8_rbf.hxx | 37 +++++++++++++++---- 3 files changed, 34 insertions(+), 8 deletions(-) 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 f14e56cd70..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 @@ -51,6 +51,8 @@ struct t8_cmesh_mesh_deformation * Computes the displacements of the surface vertices. * * \param [in] cad A pointer to the CAD-based geometry object. + * \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 @@ -65,6 +67,7 @@ struct t8_cmesh_mesh_deformation * * \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 (std::unordered_map &boundary_node_data, diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx index cca4c1abf2..6d21e2c5bb 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.cxx @@ -20,7 +20,7 @@ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ -/** \file t8_cmesh_mesh_deformation_rbf.cxx +/** \file t8_rbf.cxx * This file implements the Radial Basis Functions for the mesh deformation. */ diff --git a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx index 36c262649e..9a7c270d79 100644 --- a/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_rbf.hxx @@ -36,8 +36,13 @@ /** * The available RBF function types. */ -typedef enum { T8_RBF_CP_C2 = 0, T8_RBF_TPS, T8_RBF_FUNCTION_COUNT } t8_rbf_function_type; +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. */ @@ -63,13 +68,26 @@ struct t8_rbf_boundary_node: public t8_rbf_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; @@ -91,7 +109,8 @@ struct t8_rbf_cpc2: public t8_rbf_function * 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] 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 @@ -133,6 +152,7 @@ struct t8_rbf_tps: public t8_rbf_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 @@ -163,7 +183,7 @@ struct t8_rbf { public: /** Constructor. - * \param[in] rbf_type The RBF type to be used. + * \param [in] rbf_type The RBF type to be used. */ t8_rbf (t8_rbf_function_type rbf_type) { @@ -184,6 +204,7 @@ struct 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) @@ -199,7 +220,9 @@ struct t8_rbf } /** - * Search for the displacement of a specific boundary node with its associated global ID. + * 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 @@ -223,7 +246,7 @@ struct t8_rbf /** * Interpolates the inner node. - * \param[in, out] inner_node The inner node which will be interpolated. + * \param [in, out] inner_node The inner node which will be interpolated. */ void interpolate (t8_rbf_node& inner_node) const @@ -249,9 +272,9 @@ struct t8_rbf 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 + * \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. + * \param [in] num_boundary_nodes The number of boundary nodes. * \return The matrix of the calculated weights alpha. */ Eigen::MatrixXd From 257f2e9d84b630bf9ddb561fee250bd92804b3a5 Mon Sep 17 00:00:00 2001 From: Lena Richter Date: Tue, 6 Oct 2026 13:37:19 +0200 Subject: [PATCH 14/15] build mesh deformation only with OCC and Eigen enabled --- example/cmesh/t8_cmesh_mesh_deformation.cxx | 16 ++++++++-------- src/CMakeLists.txt | 15 ++++++++++----- src/config.cmake.in | 1 + 3 files changed, 19 insertions(+), 13 deletions(-) diff --git a/example/cmesh/t8_cmesh_mesh_deformation.cxx b/example/cmesh/t8_cmesh_mesh_deformation.cxx index a11767af88..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 T8_ENABLE_OCC +#if T8_ENABLE_OCC && T8_ENABLE_EIGEN #include #include -#endif /* T8_ENABLE_OCC */ +#endif /* T8_ENABLE_OCC and T8_ENABLE_EIGEN*/ #include #include @@ -42,7 +42,7 @@ #include #include -#if T8_ENABLE_OCC +#if T8_ENABLE_OCC && T8_ENABLE_EIGEN namespace fs = std::filesystem; static std::vector @@ -59,12 +59,12 @@ findBrepFiles (const char *folder) return files; } -#endif /* T8_ENABLE_OCC */ +#endif /* T8_ENABLE_OCC && T8_ENABLE_EIGEN */ int main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) { -#if T8_ENABLE_OCC +#if T8_ENABLE_OCC && T8_ENABLE_EIGEN char usage[BUFSIZ]; /* Brief help message. */ @@ -183,9 +183,9 @@ main ([[maybe_unused]] int argc, [[maybe_unused]] char **argv) mpiret = sc_MPI_Finalize (); SC_CHECK_MPI (mpiret); -#else /* T8_ENABLE_OCC */ - t8_global_errorf ("ERROR: This example requires OpenCASCADE support to be enabled in t8code.\n"); -#endif /* T8_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/src/CMakeLists.txt b/src/CMakeLists.txt index 085147bb9a..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,10 +107,6 @@ if( T8CODE_ENABLE_OCC ) t8_cad/t8_cad_handle.hxx DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/t8_cad ) - install( FILES - t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.hxx - DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/t8_cmesh/t8_cmesh_mesh_deformation/ - ) endif() if( T8CODE_ENABLE_EIGEN ) @@ -126,6 +121,16 @@ if( T8CODE_ENABLE_EIGEN ) ) 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/ + ) +endif() + if( T8CODE_BUILD_PEDANTIC ) target_compile_options( T8 PUBLIC -pedantic ) set (T8_CXXFLAGS "${T8_CXXFLAGS} -Wpedantic") diff --git a/src/config.cmake.in b/src/config.cmake.in index ab1a0021da..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@ ) From 6458bf3a27f116163a16f7665cd123556c06eac1 Mon Sep 17 00:00:00 2001 From: Lena-Richter <132659113+Lena-Richter@users.noreply.github.com> Date: Wed, 7 Oct 2026 23:31:56 +0200 Subject: [PATCH 15/15] Update src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.cxx Co-authored-by: Ole Albers <122293607+ole-alb@users.noreply.github.com> --- .../t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 d0df4a7dbe..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 @@ -197,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);