Path: blob/main/tutorials/general/t8_step5_element_data.cxx
903 views
/*
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) 2015 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.
*/
/* See also: https://github.com/DLR-AMR/t8code/wiki/Step-5---Store-element-data
*
* This is step5 of the t8code tutorials.
* In the following we will store data in the individual elements of our forest.
* To do this, we will again create a uniform forest, which will get adapted as in step4,
* with the difference that we partition, balance and create ghost elements all in the same step.
* After adapting the forest we will learn how to build a data array and gather data for
* the local elements. Furthermore, we exchange the data values of the ghost elements and
* output the volume data to vtu.
*
* How you can experiment here:
* - Look at the paraview output files of the adapted forest.
* You can apply a clip filter to look into the cube. Also you can apply (in addition)
* the threshold filter to display only elements with certain properties.
* But at first you may just want to enter the tooltip selection mode 'Hover Cells On'
* to display cell information when hover over them.
* - Change the adaptation criterion as you wish to adapt elements or families as desired.
* - Store even more data per element, for instance the coordinates of its midpoint.
* You can again apply the threshold filter to your new data. Don't forget to write the
* data into the output file.
* */
#include <t8.h> /* General t8code header, always include this. */
#include <t8_cmesh/t8_cmesh.h> /* cmesh definition and basic interface. */
#include <t8_cmesh/t8_cmesh_examples.h> /* A collection of exemplary cmeshes */
#include <t8_forest/t8_forest_general.h> /* forest definition and basic interface. */
#include <t8_forest/t8_forest_io.h> /* save forest */
#include <t8_forest/t8_forest_geometrical.h> /* geometrical information */
#include <t8_schemes/t8_default/t8_default.hxx> /* default refinement scheme. */
#include <tutorials/general/t8_step3.h>
T8_EXTERN_C_BEGIN ();
/* 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;
};
static t8_forest_t
t8_step5_build_forest (sc_MPI_Comm comm, int level)
{
t8_cmesh_t cmesh = t8_cmesh_new_hypercube_hybrid (comm, 0, 0);
const t8_scheme *scheme = t8_scheme_new_default ();
struct t8_step3_adapt_data adapt_data = {
{ 0.5, 0.5, 1 }, /* Midpoints of the sphere. */
0.2, /* Refine if inside this radius. */
0.4 /* Coarsen if outside this radius. */
};
/* Start with a uniform forest. */
t8_forest_t forest = t8_forest_new_uniform (cmesh, scheme, level, 0, comm);
t8_forest_t forest_apbg;
/* Adapt, partition, balance and create ghost elements all in the same step. */
t8_forest_init (&forest_apbg);
t8_forest_set_user_data (forest_apbg, &adapt_data);
t8_forest_set_adapt (forest_apbg, forest, t8_step3_adapt_callback, 0);
t8_forest_set_partition (forest_apbg, NULL, 0);
t8_forest_set_balance (forest_apbg, NULL, 0);
t8_forest_set_ghost (forest_apbg, 1, T8_GHOST_FACES);
t8_forest_commit (forest_apbg);
return forest_apbg;
}
static struct t8_step5_data_per_element *
t8_step5_create_element_data (t8_forest_t forest)
{
t8_locidx_t num_local_elements;
t8_locidx_t num_ghost_elements;
struct t8_step5_data_per_element *element_data;
/* Check that forest is a committed, that is valid and usable, forest. */
T8_ASSERT (t8_forest_is_committed (forest));
/* Get the number of local elements of forest. */
num_local_elements = t8_forest_get_local_num_leaf_elements (forest);
/* Get the number of ghost elements of forest. */
num_ghost_elements = t8_forest_get_num_ghosts (forest);
/* Now we need to build an array of our data that is as long as the number
* of elements plus the number of ghosts. You can use any allocator such as
* new, malloc or the t8code provide allocation macro T8_ALLOC.
* Note that in the latter case you need
* to use T8_FREE in order to free the memory.
*/
element_data = T8_ALLOC (struct t8_step5_data_per_element, num_local_elements + num_ghost_elements);
/* Note: We will later need to associate this data with an sc_array in order to exchange the values for
* the ghost elements, which we can do with sc_array_new_data (see t8_step5_exchange_ghost_data).
* We could also have directly allocated the data here in an sc_array with
* sc_array_new_count (sizeof (struct data_per_element), num_local_elements + num_ghost_elements);
*/
/* Let us now fill the data with something.
* For this, we iterate through all trees and for each tree through all its elements, calling
* t8_forest_get_leaf_element_in_tree to get a pointer to the current element.
* This is the recommended and most performant way.
* An alternative is to iterate over the number of local elements and use
* t8_forest_get_leaf_element. However, this function needs to perform a binary search
* for the element and the tree it is in, while t8_forest_get_leaf_element_in_tree has a
* constant look up time. You should only use t8_forest_get_leaf_element if you do not know
* in which tree an element is.
*/
{
t8_locidx_t itree, num_local_trees;
t8_locidx_t current_index;
t8_locidx_t ielement, num_elements_in_tree;
const t8_element_t *element;
/* Get the scheme of the forest for later usage. */
const t8_scheme *scheme = t8_forest_get_scheme (forest);
/* Get the number of trees that have elements of this process. */
num_local_trees = t8_forest_get_num_local_trees (forest);
for (itree = 0, current_index = 0; itree < num_local_trees; ++itree) {
/* This loop iterates through all local trees in the forest. */
/* Each tree may have a different tree class (quad/tri/hex/tet etc.) and therefore
* also a different way to interpret its elements. In order to be able to handle elements
* of a tree, we need the scheme we obtained earlier and in order to use it we get the eclass of the tree. */
const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, itree);
/* Get the number of elements of this tree. */
num_elements_in_tree = t8_forest_get_tree_num_leaf_elements (forest, itree);
for (ielement = 0; ielement < num_elements_in_tree; ++ielement, ++current_index) {
/* This loop iterates through all the local elements of the forest in the current tree. */
/* We can now write to the position current_index into our array in order to store
* data for this element. */
/* Since in this example we want to compute the data based on the element in question,
* we need to get a pointer to this element. */
element = t8_forest_get_leaf_element_in_tree (forest, itree, ielement);
/* We want to store the elements level and its volume as data. We compute these
* via the scheme and the forest_element interface. */
element_data[current_index].level = scheme->element_get_level (tree_class, element);
element_data[current_index].volume = t8_forest_element_volume (forest, itree, element);
}
}
}
return element_data;
}
/* Each process has computed the data entries for its local elements.
* In order to get the values for the ghost elements, we use t8_forest_ghost_exchange_data.
* Calling this function will fill all the ghost entries of our element data array with the
* value on the process that owns the corresponding element. */
static void
t8_step5_exchange_ghost_data (t8_forest_t forest, struct t8_step5_data_per_element *data)
{
sc_array *sc_array_wrapper;
t8_locidx_t num_elements = t8_forest_get_local_num_leaf_elements (forest);
t8_locidx_t num_ghosts = t8_forest_get_num_ghosts (forest);
/* t8_forest_ghost_exchange_data expects an sc_array (of length num_local_elements + num_ghosts).
* We wrap our data array to an sc_array. */
sc_array_wrapper = sc_array_new_data (data, sizeof (struct t8_step5_data_per_element), num_elements + num_ghosts);
/* Carry out the data exchange. The entries with indices > num_local_elements will get overwritten.
*/
t8_forest_ghost_exchange_data (forest, sc_array_wrapper);
/* Destroy the wrapper array. This will not free the data memory since we used sc_array_new_data. */
sc_array_destroy (sc_array_wrapper);
}
/* Write the forest 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
*/
static void
t8_step5_output_data_to_vtu (t8_forest_t forest, struct t8_step5_data_per_element *data, const char *prefix)
{
t8_locidx_t num_elements = t8_forest_get_local_num_leaf_elements (forest);
t8_locidx_t ielem;
/* 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 (ielem = 0; ielem < num_elements; ++ielem) {
element_volumes[ielem] = data[ielem].volume;
}
{
/* To write user defined data, we need the extended output function t8_forest_vtk_write_file
* from t8_forest_vtk.h. Despite writing user data, it also offers more control over which
* properties of the forest to write. */
int write_treeid = 1;
int write_mpirank = 1;
int write_level = 1;
int write_element_id = 1;
int write_ghosts = 0;
t8_forest_write_vtk_ext (forest, prefix, write_treeid, write_mpirank, write_level, write_element_id, write_ghosts,
0, 0, num_data, &vtk_data);
}
T8_FREE (element_volumes);
}
int
t8_step5_main (int argc, char **argv)
{
int mpiret;
sc_MPI_Comm comm;
t8_forest_t forest;
/* The prefix for our output files. */
const char *prefix_forest = "t8_step5_forest";
const char *prefix_forest_with_data = "t8_step5_forest_with_volume_data";
/* The uniform refinement level of the forest. */
const int level = 3;
/* The array that will hold our per element data. */
t8_step5_data_per_element *data;
/* Initialize MPI. This has to happen before we initialize sc or t8code. */
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);
/* Print a message on the root process. */
t8_global_productionf (" [step5] \n");
t8_global_productionf (" [step5] Hello, this is the step5 example of t8code.\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");
/* We will use MPI_COMM_WORLD as a communicator. */
comm = sc_MPI_COMM_WORLD;
/*
* Setup.
* Build cmesh and uniform forest.
* Adapt forest similar to step 3 & 4.
*/
t8_global_productionf (" [step5] \n");
t8_global_productionf (" [step5] Creating an adapted forest as in step3.\n");
t8_global_productionf (" [step5] \n");
forest = t8_step5_build_forest (comm, level);
t8_forest_write_vtk (forest, prefix_forest);
t8_global_productionf (" [step5] Wrote forest to vtu files: %s*\n", prefix_forest);
/*
* Build data array and gather data for the local elements.
*/
data = t8_step5_create_element_data (forest);
t8_global_productionf (" [step5] Computed level and volume data for local elements.\n");
if (t8_forest_get_local_num_leaf_elements (forest) > 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", data[0].level, data[0].volume);
}
/*
* Exchange the data values of the ghost elements
*/
t8_step5_exchange_ghost_data (forest, data);
t8_global_productionf (" [step5] Exchanged ghost data.\n");
if (t8_forest_get_num_ghosts (forest) > 0) {
/* output the data of the first ghost element (if it exists) */
t8_locidx_t first_ghost_index = t8_forest_get_local_num_leaf_elements (forest);
t8_global_productionf (" [step5] Ghost 0 has level %i and volume %e.\n", data[first_ghost_index].level,
data[first_ghost_index].volume);
}
/*
* Output the volume data to vtu.
*/
t8_step5_output_data_to_vtu (forest, data, prefix_forest_with_data);
t8_global_productionf (" [step5] Wrote forest and volume data to %s*.\n", prefix_forest_with_data);
/*
* clean-up
*/
/* Free the data array. */
T8_FREE (data);
/* Destroy the forest. */
t8_forest_unref (&forest);
t8_global_productionf (" [step5] Destroyed forest.\n");
sc_finalize ();
mpiret = sc_MPI_Finalize ();
SC_CHECK_MPI (mpiret);
return 0;
}
T8_EXTERN_C_END ();