diff --git a/mesh_handle/README.md b/mesh_handle/README.md index 20f220981d..fc578bf32e 100644 --- a/mesh_handle/README.md +++ b/mesh_handle/README.md @@ -12,5 +12,6 @@ The folder's most important files are: - The [element.hxx](element.hxx) defines the elements (mesh or ghost elements) of the mesh handle. - The [competences.hxx](competences.hxx) defines additional competences/functionality of an element to access additional data. - The [competence_pack.hxx](competence_pack.hxx) is needed to pack element or mesh competences to pass it as template parameter to the mesh. +- The [mesh_io.hxx](mesh_io.hxx) allows to output results in vtk format. Headers in the [internal/](internal/) folder are only intended to implement details of the mesh handle, so they should not need to be included for any other purpose. \ No newline at end of file diff --git a/mesh_handle/mesh_io.hxx b/mesh_handle/mesh_io.hxx new file mode 100644 index 0000000000..6fa2c16846 --- /dev/null +++ b/mesh_handle/mesh_io.hxx @@ -0,0 +1,92 @@ +/* + 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 mesh_io.hxx + * In- and output of meshes. For example, writing a mesh to vtk format. + */ + +#pragma once +#include + +namespace t8_mesh_handle +{ +/** + * Write the mesh in a parallel vtu format. Extended version. + * See \see write_mesh_to_vtk for the standard version of this function. + * Writes one master .pvtu file and each process writes in its own .vtu file. + * If linked and not otherwise specified, the VTK API is used. + * If the VTK library is not linked, an ASCII file is written. + * This may change in accordance with \a write_ghosts, \a write_curved and + * \a do_not_use_API, because the export of ghosts is not yet available with + * the VTK API and the export of curved elements is not available with the + * inbuilt function to write ASCII files. The function will for example + * still use the VTK API to satisfy \a write_curved, even if \a do_not_use_API + * is set to true. + * This function is collective and must be called on each process. + * \param [in] mesh The mesh to write. + * \param [in] fileprefix The prefix of the files where the vtk will be stored. + * The master file is then fileprefix.pvtu and the process with rank r writes in the file fileprefix_r.vtu. + * \param [in] num_data Number of user defined double valued data fields to write. + * \param [in] data Array of t8_vtk_data_field_t of length \a num_data providing the user defined + * per element data. If scalar and vector fields are used, all scalar fields must come first in the array. + * \param [in] write_treeid If true, the global tree id of the underlying forest is written for each element. + * \param [in] write_mpirank If true, the mpirank is written for each element. + * \param [in] write_level If true, the refinement level is written for each element. + * \param [in] write_element_id If true, the global element id is written for each element. + * \param [in] write_ghosts If true, each process additionally writes its ghost elements. + * For ghost element the treeid of the underlying forest is -1. + * \param [in] write_curved If true, write the elements as curved element types from vtk. + * \param [in] do_not_use_API Do not use the VTK API, even if linked and available. + * \return True if successful, false if not (process local). + */ +template +int +write_mesh_to_vtk_ext (TMeshClass &mesh, const char *fileprefix, const int num_data, t8_vtk_data_field_t *data, + bool write_treeid = false, bool write_mpirank = true, bool write_level = true, + bool write_element_id = true, bool write_ghosts = false, bool write_curved = false, + bool do_not_use_API = false) +{ + return t8_forest_write_vtk_ext (mesh->get_forest (), fileprefix, write_treeid, write_mpirank, write_level, + write_element_id, write_ghosts, write_curved, do_not_use_API, num_data, data); +} + +/** + * Write the mesh in a parallel vtu format. Writes one master + * .pvtu file and each process writes in its own .vtu file. + * If linked, the VTK API is used. + * If the VTK library is not linked, an ASCII file is written. + * This function writes the elements, level, mpirank, element id and the tree id of the underlying forest as data. + * This function is collective and must be called on each process. + * For more options use \see write_mesh_to_vtk_ext. + * \param [in] mesh The mesh to write. + * \param [in] fileprefix The prefix of the files where the vtk will be stored. + * The master file is then fileprefix.pvtu and the process with rank r writes in the file fileprefix_r.vtu. + * \return True if successful, false if not (process local). + */ +template +int +write_mesh_to_vtk (TMeshClass &mesh, const char *fileprefix) +{ + return t8_forest_write_vtk (mesh->get_forest (), fileprefix); +} + +} // namespace t8_mesh_handle diff --git a/tutorials/CMakeLists.txt b/tutorials/CMakeLists.txt index 5d3a1d1d87..276c1759b3 100644 --- a/tutorials/CMakeLists.txt +++ b/tutorials/CMakeLists.txt @@ -33,6 +33,10 @@ function( add_t8_tutorial ) add_executable( ${ADD_T8_TUTORIAL_NAME} ${ADD_T8_TUTORIAL_SOURCES} ) target_link_libraries( ${ADD_T8_TUTORIAL_NAME} PRIVATE T8 SC::SC ) + string( FIND ${ADD_T8_TUTORIAL_NAME} "mesh_handle" is_mesh_handle_file ) + if ( (${is_mesh_handle_file} GREATER_EQUAL 0) ) + target_link_libraries(${ADD_T8_TUTORIAL_NAME} PRIVATE T8_MESH_HANDLE ) + endif () target_include_directories( ${ADD_T8_TUTORIAL_NAME} PRIVATE ${CMAKE_CURRENT_LIST_DIR}/.. ) set_target_properties(${ADD_T8_TUTORIAL_NAME} PROPERTIES @@ -69,3 +73,7 @@ copy_tutorial_file (features/t8_features_curved_meshes_generate_cmesh_hex.geo) copy_tutorial_file (features/t8_features_curved_meshes_generate_cmesh_quad.geo) copy_tutorial_file (features/t8_features_curved_meshes_generate_cmesh_tet.geo) copy_tutorial_file (features/t8_features_curved_meshes_generate_cmesh_tri.geo) + +if( T8CODE_BUILD_MESH_HANDLE ) + add_t8_tutorial( NAME t8_step5_mesh_handle SOURCES general/t8_step5_mesh_handle.cxx ) +endif() diff --git a/tutorials/general/t8_step5_mesh_handle.cxx b/tutorials/general/t8_step5_mesh_handle.cxx new file mode 100644 index 0000000000..9467532c21 --- /dev/null +++ b/tutorials/general/t8_step5_mesh_handle.cxx @@ -0,0 +1,239 @@ +/* + 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 types 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_step5_mesh_handle.cxx + * + * This is the same as t8_step5_element_data.cxx but using the mesh handle interface instead of the forest interface. + */ + +#include + +#include +#include +#include +#include +#include +#include +#include + +/* The data that we want to store for each element. + * In this example we want to store the element's level and volume. */ +struct t8_step5_data_per_element +{ + int level; + double volume; +}; + +/** User data type we will pass to the adapt callback. */ +struct user_data +{ + t8_3D_vec midpoint; /**< The midpoint of our sphere. */ + double refine_if_inside_radius; /**< If an element's center is smaller than this value, we refine the element. */ + double coarsen_if_outside_radius; /**< If an element's center is larger this value, we coarsen its family. */ +}; + +/** The adaptation callback function. This will refine elements inside of a given sphere and coarsen the elements + * outside of a given sphere. + * \tparam TMeshClass The mesh handle class. + * \param [in] mesh The mesh that should be adapted. + * \param [in] elements One element or a family of elements to consider for adaptation. + * \param [in] user_data The user data to be used during the adaptation process. + * \return 1 if the first entry in \a elements should be refined, + * -1 if the family \a elements shall be coarsened, + * 0 else. + */ +template +int +adapt_callback ([[maybe_unused]] const TMeshClass &mesh, std::span elements, + const user_data &user_data) +{ + auto element_centroid = elements[0].get_centroid (); + double dist = t8_dist (element_centroid, user_data.midpoint); + if (dist < user_data.refine_if_inside_radius) { + return 1; + } + // Check if we got a family and if yes, if we should coarsen. + if ((elements.size () > 1) && (dist > user_data.coarsen_if_outside_radius)) { + return -1; + } + return 0; +} + +/** Build a mesh with initial uniform refinement level \a level which is adapted according to \ref adapt_callback, + * partitioned and balanced afterwards and ghost elements are set. + * \tparam TMeshClass The mesh handle class. + * \param [in] comm MPI communicator to use. + * \param [in] level Initial refinement level. + * \return Unique pointer to the mesh created. + */ +template +std::unique_ptr +t8_step5_build_mesh (sc_MPI_Comm comm, int level) +{ + auto mesh_handle = t8_mesh_handle::handle_hypercube_hybrid_uniform_default (level, comm); + struct user_data adapt_data = { + { 0.5, 0.5, 1 }, /* Midpoint of the sphere. */ + 0.2, /* Refine if inside this radius. */ + 0.4 /* Coarsen if outside this radius. */ + }; + /* Adapt, partition, balance and create ghost elements. */ + mesh_handle->set_balance (); + mesh_handle->set_partition (); + mesh_handle->set_adapt (TMesh::template mesh_adapt_callback_wrapper (adapt_callback, adapt_data), + false); + mesh_handle->set_ghost (); + mesh_handle->commit (); + return mesh_handle; +} + +/** Set element data to the mesh handle. */ +template +void +t8_step5_set_element_data (TMeshClass &mesh) +{ + for (auto &elem : *mesh) { + elem.set_element_data ({ elem.get_level (), elem.get_volume () }); + } +} + +/** Exchange element data set in \ref t8_step5_set_element_data for ghost elements. */ +template +void +t8_step5_exchange_ghost_data (TMeshClass &mesh) +{ + mesh->exchange_ghost_data (); +} + +/** Write the mesh as vtu and also write the element's volumes in the file. + * t8code supports writing element based data to vtu as long as its stored + * as doubles. Each of the data fields to write has to be provided in its own + * array of length num_local_elements. + * We support two types: T8_VTK_SCALAR - One double per element + * and T8_VTK_VECTOR - 3 doubles per element + */ +template +static void +t8_step5_output_data_to_vtu (TMeshClass &mesh, const char *prefix) +{ + t8_locidx_t num_elements = mesh->get_num_local_elements (); + /* We need to allocate a new array to store the volumes on their own. + * This array has one entry per local element. */ + double *element_volumes = T8_ALLOC (double, num_elements); + /* The number of user defined data fields to write. */ + int num_data = 1; + /* For each user defined data field we need one t8_vtk_data_field_t variable */ + t8_vtk_data_field_t vtk_data; + /* Set the type of this variable. Since we have one value per element, we pick T8_VTK_SCALAR. */ + vtk_data.type = T8_VTK_SCALAR; + /* The name of the field as should be written to the file. */ + strcpy (vtk_data.description, "Element volume"); + vtk_data.data = element_volumes; + /* Copy the element's volumes from our data array to the output array. */ + for (t8_locidx_t ielem = 0; ielem < num_elements; ++ielem) { + element_volumes[ielem] = ((*mesh)[ielem]).get_element_data ().volume; + } + /* To write user defined data, we need the extended output function write_mesh_to_vtk_ext. + * Despite writing user data, it also offers more control over which properties to write. */ + t8_mesh_handle::write_mesh_to_vtk_ext (mesh, prefix, num_data, &vtk_data); + T8_FREE (element_volumes); +} + +/** Main function. */ +int +t8_step5_main (int argc, char **argv) +{ + /* The prefix for our output files. */ + const char *prefix_mesh = "t8_step5_handle"; + const char *prefix_mesh_with_data = "t8_step5_handle_with_volume_data"; + /* The initial uniform refinement level. */ + const int level = 3; + + /* 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_PRODUCTION. See sc.h for more info on the log levels. */ + t8_init (SC_LP_PRODUCTION); + /* We will use MPI_COMM_WORLD as a communicator. */ + sc_MPI_Comm comm = sc_MPI_COMM_WORLD; + + /* Print a message on the root process. */ + t8_global_productionf (" [step5] \n"); + t8_global_productionf (" [step5] Hello, this is the step5 example of t8code using the mesh handle.\n"); + t8_global_productionf ( + " [step5] In this example we will store data on our elements and exchange the data of ghost elements.\n"); + t8_global_productionf (" [step5] \n"); + + /* Setup: Build cmesh and adapt uniformly. Adapt similar to step 3 & 4. */ + t8_global_productionf (" [step5] \n"); + t8_global_productionf (" [step5] Creating an adapted mesh as in step3.\n"); + t8_global_productionf (" [step5] \n"); + { /* We put the mesh in its own scope so that it is automatically destroyed at the end of the scope. + * This is only necessary because sc_finalize checks if there are leftover references. + * This unique pointer would have been destroyed automatically at the end of the programme. */ + using mesh_class = t8_mesh_handle::mesh, t8_step5_data_per_element>; + auto mesh = t8_step5_build_mesh (comm, level); + + t8_mesh_handle::write_mesh_to_vtk (mesh, prefix_mesh); + t8_global_productionf (" [step5] Wrote mesh to vtu files: %s*\n", prefix_mesh); + + t8_step5_set_element_data (mesh); + t8_global_productionf (" [step5] Computed level and volume data for local elements.\n"); + if (mesh->get_num_local_elements () > 0) { + /* Output the stored data of the first local element (if it exists). */ + t8_global_productionf (" [step5] Element 0 has level %i and volume %e.\n", ((*mesh)[0]).get_element_data ().level, + ((*mesh)[0]).get_element_data ().volume); + } + + /* Exchange the data values of the ghost elements. */ + t8_step5_exchange_ghost_data (mesh); + t8_global_productionf (" [step5] Exchanged ghost data.\n"); + if (mesh->get_num_ghosts () > 0) { + /* Output the data of the first ghost element (if it exists). */ + t8_locidx_t first_ghost_index = mesh->get_num_local_elements (); + t8_global_productionf (" [step5] Ghost 0 has level %i and volume %e.\n", + ((*mesh)[first_ghost_index]).get_element_data ().level, + ((*mesh)[first_ghost_index]).get_element_data ().volume); + } + + /* Output the volume data to vtu. */ + t8_step5_output_data_to_vtu (mesh, prefix_mesh_with_data); + t8_global_productionf (" [step5] Wrote mesh and volume data to %s*.\n", prefix_mesh_with_data); + + /* Cleanup. */ + } // End scope of mesh + sc_finalize (); + + mpiret = sc_MPI_Finalize (); + SC_CHECK_MPI (mpiret); + + return 0; +} + +/** Entry point of the program. */ +int +main (int argc, char **argv) +{ + return t8_step5_main (argc, argv); +}