Feature: Feature mesh deformation

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 )

example/cmesh/t8_cmesh_mesh_deformation_example.cxx

@@ -0,0 +1,118 @@
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <algorithm>
+#include <t8_cmesh/t8_cmesh.h>
+#include <t8_cmesh/t8_cmesh.hxx>
+#include <t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_connectivity.hxx>
+#include <t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.hxx>
+#include <t8_schemes/t8_default/t8_default.hxx>
+#include <t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.h>
+#include <t8_cad/t8_cad_handle.hxx>
+#include <t8_vtk/t8_vtk_writer.h>
+#include <vector>
+#include <array>
+#include <mpi.h>
+#include <unordered_set>
+#include <sc_options.h>
+
+int
+main (int argc, char **argv)
+{
+  char usage[BUFSIZ];
+  /* Brief help message. */
+  int sreturnA = snprintf (usage, BUFSIZ, "Usage:\t%s <OPTIONS>\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<t8_cad_handle> (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;
+}

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

src/t8_cad/t8_cad_handle.hxx

@@ -28,6 +28,9 @@
 #ifndef T8_CAD_HANDLE_HXX
 #define T8_CAD_HANDLE_HXX
 
+#include <t8.h>
+#include <optional>
+#include <span>
 #include <gp_Pnt.hxx>
 #include <TopExp.hxx>
 #include <Geom_Surface.hxx>
@@ -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.

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 <unordered_map>
+#include <t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_connectivity.hxx>
+#include <t8_cad/t8_cad_handle.hxx>
+#include <t8_cmesh/t8_cmesh.h>
+#include <t8_cmesh/t8_cmesh.hxx>
+#include <t8_cmesh/t8_cmesh_mesh_deformation/t8_cmesh_mesh_deformation.hxx>
+#include <t8_geometry/t8_geometry_handler.hxx>
+#include <t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx>
+/**
+ * Calculate the displacement of vertices for CAD-based mesh deformation.
+ */
+
+std::unordered_map<t8_gloidx_t, t8_3D_vec>
+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<t8_gloidx_t, t8_3D_vec> 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<const int *> (
+      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<const int *> (
+        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<t8_gloidx_t, t8_3D_vec> &displacements,
+                            std::shared_ptr<t8_cad_handle> 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<t8_geometry_cad *> (geom->second.get ());
+      cad_geom->update_cad_manager (cad);
+      break;
+    }
+  }
+}