diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index 272486619b..95628142b7 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_example SOURCES cmesh/t8_cmesh_mesh_deformation_example.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_example.cxx b/example/cmesh/t8_cmesh_mesh_deformation_example.cxx new file mode 100644 index 0000000000..7a4154c85c --- /dev/null +++ b/example/cmesh/t8_cmesh_mesh_deformation_example.cxx @@ -0,0 +1,118 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +int +main (int argc, char **argv) +{ + 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_PRODUCTION); + + /* Initialize t8code with log level SC_LP_PRODUCTION. See sc.h for more info on the log levels. */ + t8_init (SC_LP_PRODUCTION); + + int helpme = 0; + const char *msh_file = NULL; + const char *brep_file = NULL; + int dim = 0; + + /* 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 (.brep)"); + sc_options_add_int (opt, 'd', "dimension", &dim, 0, "Dimension of the mesh (1, 2 or 3)"); + + 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_productionf ("\n\t 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_productionf ("\n\t 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 (), 2, 0, comm); + + /* Load CAD geometry from .brep file. */ + auto cad = std::make_shared (brep_file); + + /* Calculate displacements. */ + auto displacements = calculate_displacement_surface_vertices (cmesh, cad.get ()); + + /* Write output. */ + t8_forest_vtk_write_file (forest, "input_forest", 1, 1, 1, 1, 0, 0, NULL); + + /* Apply displacements. */ + apply_vertex_displacements (cmesh, displacements, cad); + + /* Write output. */ + t8_forest_vtk_write_file (forest, "deformed_forest", 1, 1, 1, 1, 0, 0, NULL); + + /* Cleanup. */ + t8_forest_unref (&forest); + + std::cout << "Mesh deformation completed." << std::endl; + } + + MPI_Finalize (); + sc_options_destroy (opt); + return 0; +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index cd7b1e9db1..872a786bb5 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 @@ -153,6 +154,7 @@ target_sources( T8 PRIVATE t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.cxx t8_cmesh/t8_cmesh_io/t8_cmesh_save.cxx t8_cmesh/t8_cmesh_io/t8_cmesh_triangle.cxx + t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.cxx t8_data/t8_shmem.c t8_data/t8_containers.cxx t8_eclass/t8_eclass.c @@ -162,7 +164,7 @@ target_sources( T8 PRIVATE t8_forest/t8_forest_partition.cxx t8_forest/t8_forest_partition_for_coarsening.cxx t8_forest/t8_forest_pfc_helper.cxx - t8_forest/t8_forest.cxx + t8_forest/t8_forest.cxx t8_forest/t8_forest_private.cxx t8_forest/t8_forest_ghost.cxx t8_forest/t8_forest_iterate.cxx diff --git a/src/t8_cad/t8_cad_handle.hxx b/src/t8_cad/t8_cad_handle.hxx index ec7c3c3809..ef14f72654 100644 --- a/src/t8_cad/t8_cad_handle.hxx +++ b/src/t8_cad/t8_cad_handle.hxx @@ -28,6 +28,9 @@ #ifndef T8_CAD_HANDLE_HXX #define T8_CAD_HANDLE_HXX +#include +#include +#include #include #include #include @@ -55,6 +58,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..72838f8f33 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.cxx @@ -0,0 +1,198 @@ + +/* + 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) 2025 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 +/** + * Calculate the displacement of vertices for CAD-based mesh deformation. + */ + +std::unordered_map +calculate_displacement_surface_vertices (t8_cmesh_t cmesh, const t8_cad_handle *cad) +{ + T8_ASSERT (t8_cmesh_is_committed (cmesh)); + + /* Get number of global vertices on this rank */ + //const t8_gloidx_t num_global_vertices = t8_cmesh_get_num_global_vertices (cmesh); + + const int mesh_dimension = t8_cmesh_get_dimension (cmesh); + + /* Map from global vertex id -> displacement vector. */ + std::unordered_map displacements; + + for (const auto &global_vertex : *(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 (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 (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 ( + 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 ( + 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; +} +/** + * Apply vertex displacements to a cmesh and update the CAD geometry. + */ +void +apply_vertex_displacements (t8_cmesh_t cmesh, const std::unordered_map &displacements, + std::shared_ptr cad) +{ + T8_ASSERT (t8_cmesh_is_committed (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 = 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 (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 = 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..2e9e3148b2 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.hxx @@ -0,0 +1,80 @@ + +/* + 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) 2025 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. + */ + +#ifndef T8_CMESH_DEFORMATION_HXX +#define T8_CMESH_DEFORMATION_HXX + +#include +#include +#include + +#include +#include +#include +#include + +/** Struct for mesh deformation. */ +struct t8_cmesh_mesh_deformation +{ + public: + /** Constructor */ + t8_cmesh_mesh_deformation (): associated_cmesh (nullptr), updated_geometry (nullptr) {}; + + /** Destructor */ + ~t8_cmesh_mesh_deformation () {}; + + 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; +}; + +/** + * Computes the displacements of the surface vertices. + * + * \param [in] cmesh The coarse mesh structure. + * \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 (t8_cmesh_t cmesh, 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] cmesh The committed coarse mesh structure. + * \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 (t8_cmesh_t cmesh, const std::unordered_map &displacements, + std::shared_ptr cad); +#endif /* !T8_CMESH_DEFORMATION_HXX */ 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..b1f9237b7c 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,23 @@ 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. */ + inline const_iterator + begin () const + { + return vertex_to_tree.begin (); + } + + /** Iterator end. */ + 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..4d3172073f 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx @@ -152,6 +152,15 @@ struct t8_geometry_cad: public t8_geometry_with_vertices return cad_handle; } + /** Update the CAD manager 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. diff --git a/tutorials/features/t8_features_curved_meshes.cxx b/tutorials/features/t8_features_curved_meshes.cxx index 58eb3a622c..684aa3fcd6 100644 --- a/tutorials/features/t8_features_curved_meshes.cxx +++ b/tutorials/features/t8_features_curved_meshes.cxx @@ -45,7 +45,9 @@ #include /* default refinement scheme. */ #include /* Linear geometry calculation of trees */ #if T8_ENABLE_OCC -#include /* Curved geometry calculation of trees */ +#include /* Curved geometry calculation of trees */ +#include /* Mesh deformation struct */ +#include /* CAD data structure */ #endif #include /* msh file reader */ #include /* std::string */ @@ -300,7 +302,7 @@ t8_naca_plane_refinement (t8_forest_t forest, const std::string &fileprefix, int /* Moving plane loop */ while (adapt_data.t < steps) { - /* Adapt and balance the forest. + /* Adapt and balance the forest. * Note, that we have to hand the adapt data to the forest before the commit. */ t8_forest_init (&forest_new); t8_forest_set_adapt (forest_new, forest, t8_naca_plane_adapt_callback, 1);