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/example/cmesh/t8_cmesh_mesh_deformation.cxx b/example/cmesh/t8_cmesh_mesh_deformation.cxx new file mode 100644 index 0000000000..4d8e236f9f --- /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 /* T8CODE_ENABLE_OCC */ +#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 /* T8CODE_ENABLE_OCC */ + t8_global_errorf ("ERROR: This example requires OpenCASCADE support to be enabled in t8code.\n"); +#endif /* T8CODE_ENABLE_OCC */ + + return 0; +} 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..8f1bfd0c06 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.cxx @@ -0,0 +1,190 @@ +/* + 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]; + + /* Check if all trees sharing this vertex have consistent geometry attributes. */ +#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.