Program Listing for File t8_vtk_writer.hxx
↰ Return to documentation for file (src/t8_vtk/t8_vtk_writer.hxx)
/*
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) 2024 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.
*/
#ifndef T8_VTK_WRITER_HXX
#define T8_VTK_WRITER_HXX
#include <t8_cmesh/t8_cmesh.h>
#include <t8_forest/t8_forest.h>
#include <t8_forest/t8_forest_general.h>
#include <t8_forest/t8_forest_ghost.h>
#include <t8.h>
#include "t8_forest/t8_forest_types.h"
#include "t8_vtk/t8_vtk_writer_helper.hxx"
#include "t8_vtk/t8_vtk_write_ASCII.hxx"
#include <string>
#include <t8_vtk.h>
#include <t8_types/t8_vec.hxx>
#if T8_ENABLE_VTK
#include <vtkUnstructuredGrid.h>
#include <vtkXMLPUnstructuredGridWriter.h>
#include <vtkDoubleArray.h>
#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkMergePoints.h>
#include <vtkTypeInt64Array.h>
#if T8_ENABLE_MPI
#include <vtkMPI.h>
#include <vtkMPICommunicator.h>
#include <vtkMPIController.h>
#endif /* T8_ENABLE_MPI */
#endif /* T8_ENABLE_VTK */
template <typename grid_t>
class vtk_writer {
public:
vtk_writer (const bool write_treeid, const bool write_mpirank, const bool write_level, const bool write_element_id,
const bool write_ghosts, const bool curved_flag, std::string fileprefix, const int num_data,
t8_vtk_data_field_t *data, sc_MPI_Comm comm, const bool merge_points = true)
: write_treeid (write_treeid), write_mpirank (write_mpirank), write_level (write_level),
write_element_id (write_element_id), write_ghosts (write_ghosts), curved_flag (curved_flag),
fileprefix (fileprefix), num_data (num_data), data (data), comm (comm), merge_points (merge_points)
{
}
vtk_writer (std::string fileprefix, sc_MPI_Comm comm): fileprefix (fileprefix), comm (comm)
{
}
#if T8_ENABLE_VTK
void
grid_to_vtkUnstructuredGrid (const grid_t grid, vtkSmartPointer<vtkUnstructuredGrid> unstructuredGrid)
{
this->t8_grid_to_vtkUnstructuredGrid (grid, unstructuredGrid);
}
#endif /* T8_ENABLE_VTK */
bool
write_with_API (const grid_t grid)
{
return write_vtk (grid);
}
bool
write_ASCII (const grid_t grid);
inline void
set_write_treeid (const bool write_treeid)
{
this->write_treeid = write_treeid;
}
inline void
set_write_mpirank (const bool write_mpirank)
{
this->write_mpirank = write_mpirank;
}
inline void
set_write_level (const bool write_level)
{
this->write_level = write_level;
}
inline void
set_write_element_id (const bool write_element_id)
{
this->write_element_id = write_element_id;
}
inline void
set_write_ghosts (const bool write_ghosts)
{
this->write_ghosts = write_ghosts;
}
inline void
set_curved_flag (const bool curved_flag)
{
this->curved_flag = curved_flag;
}
inline void
set_fileprefix (std::string fileprefix)
{
this->fileprefix = fileprefix;
}
private:
#if T8_ENABLE_VTK
void
t8_grid_element_to_vtk_cell (const grid_t grid, const t8_element_t *element, const t8_locidx_t itree,
const t8_gloidx_t offset, const bool is_ghost, const int elem_id, long int *point_id,
int *cellTypes, vtkSmartPointer<vtkMergePoints> points,
vtkSmartPointer<vtkCellArray> cellArray,
vtkSmartPointer<t8_vtk_gloidx_array_type_t> vtk_treeid,
vtkSmartPointer<t8_vtk_gloidx_array_type_t> vtk_mpirank,
vtkSmartPointer<t8_vtk_gloidx_array_type_t> vtk_level,
vtkSmartPointer<t8_vtk_gloidx_array_type_t> vtk_element_id, bool merge_points = true)
{
/* Get the shape of the current element and the respective shape of the vtk_cell. */
const t8_element_shape_t element_shape = grid_element_shape (grid, itree, element);
const int num_node = t8_get_number_of_vtk_nodes (element_shape, curved_flag);
/* Create an array with the correct number of indices per cell */
vtkIdType *vecCellIds = new vtkIdType[num_node];
/* Compute the coordinates of the element/tree. */
double *coordinates = T8_ALLOC (double, 3 * num_node);
grid_element_to_coords (grid, itree, element, curved_flag, coordinates, num_node, element_shape);
vtkIdType ptId = -1;
for (int ivertex = 0; ivertex < num_node; ivertex++, (*point_id)++) {
const size_t offset_3d = 3 * ivertex;
/* Insert the point in the points array. */
double vtkCoords[3] = { coordinates[offset_3d], coordinates[offset_3d + 1], coordinates[offset_3d + 2] };
if (merge_points) {
points->InsertUniquePoint (vtkCoords, ptId);
}
else {
ptId = points->InsertNextPoint (vtkCoords);
}
/* Add the returned point id to the cell ids*/
vecCellIds[ivertex] = ptId;
}
T8_FREE (coordinates);
/* Fill the cell array. */
cellArray->InsertNextCell (num_node, vecCellIds);
delete[] vecCellIds;
/* Write additional information if desired. */
if (curved_flag == 0) {
cellTypes[elem_id - offset] = t8_eclass_vtk_type[element_shape];
}
else {
cellTypes[elem_id - offset] = t8_curved_eclass_vtk_type[element_shape];
}
if (write_treeid == 1) {
const t8_gloidx_t gtree_id = tree_local_to_global_id (grid, itree);
if (is_ghost) {
vtk_treeid->InsertNextValue (-1);
}
else {
vtk_treeid->InsertNextValue (gtree_id);
}
}
if (write_mpirank == 1) {
int mpirank;
int mpiret = sc_MPI_Comm_rank (this->comm, &mpirank);
SC_CHECK_MPI (mpiret);
vtk_mpirank->InsertNextValue (mpirank);
}
if (write_level == 1) {
vtk_level->InsertNextValue (grid_element_level (grid, itree, element));
}
if (write_element_id == 1) {
vtk_element_id->InsertNextValue (elem_id);
}
}
void
t8_grid_tree_to_vtk_cells (const grid_t grid, vtkSmartPointer<vtkUnstructuredGrid> unstructuredGrid,
vtkSmartPointer<t8_vtk_gloidx_array_type_t> vtk_treeid,
vtkSmartPointer<t8_vtk_gloidx_array_type_t> vtk_mpirank,
vtkSmartPointer<t8_vtk_gloidx_array_type_t> vtk_level,
vtkSmartPointer<t8_vtk_gloidx_array_type_t> vtk_element_id,
vtkSmartPointer<vtkCellArray> cellArray, vtkSmartPointer<vtkMergePoints> points,
int *cellTypes, const t8_locidx_t num_local_trees, t8_gloidx_t *elem_id,
long int *point_id, const t8_gloidx_t offset, const bool ghosts, const t8_locidx_t itree,
const bool merge_points = true);
void
t8_grid_to_vtkUnstructuredGrid (const grid_t grid, vtkSmartPointer<vtkUnstructuredGrid> unstructuredGrid,
bool merge_points = true)
{
T8_ASSERT (grid != NULL);
vtkSmartPointer<vtkCellArray> cellArray = vtkSmartPointer<vtkCellArray>::New ();
vtkSmartPointer<vtkPoints> points_store = vtkSmartPointer<vtkPoints>::New ();
vtkSmartPointer<vtkMergePoints> points = vtkSmartPointer<vtkMergePoints>::New ();
vtkSmartPointer<t8_vtk_gloidx_array_type_t> vtk_treeid = vtkSmartPointer<t8_vtk_gloidx_array_type_t>::New ();
vtkSmartPointer<t8_vtk_gloidx_array_type_t> vtk_mpirank = vtkSmartPointer<t8_vtk_gloidx_array_type_t>::New ();
vtkSmartPointer<t8_vtk_gloidx_array_type_t> vtk_level = vtkSmartPointer<t8_vtk_gloidx_array_type_t>::New ();
vtkSmartPointer<t8_vtk_gloidx_array_type_t> vtk_element_id = vtkSmartPointer<t8_vtk_gloidx_array_type_t>::New ();
/*
* We need the dataArray for writing double valued user defined data in the vtu files.
* We want to write num_data many timesteps/arrays.
* We need num_data many vtkDoubleArrays, so we need to allocate storage.
* Later we call the constructor with: dataArrays[idata]=vtkDoubleArray::New().
*/
vtkDoubleArray **dataArrays;
dataArrays = T8_ALLOC (vtkDoubleArray *, num_data);
long int point_id = 0; /* The id of the point in the points Object. */
const t8_gloidx_t offset
= grid_first_local_id (grid); /* offset to take the elements of the previous processes into account*/
t8_gloidx_t elem_id = offset;
/* Check if we have to write ghosts on this process. */
const bool do_ghosts = grid_do_ghosts (grid, write_ghosts);
/* Compute the number of cells on this process. */
const t8_locidx_t num_cells = num_cells_to_write (grid, do_ghosts);
int *cellTypes = T8_ALLOC (int, num_cells);
T8_ASSERT (cellTypes != NULL);
/* Get the local bounds of the forest or cmesh */
double bounds[6] = { 0 };
if (num_cells > 0 && merge_points) {
// This function expects a non-empty partition, so we only call it
// if we have cells. If the local bounds cannot be computed,
// merge points will be disabled.
merge_points = grid_get_local_bounds (grid, bounds);
if (!merge_points) {
t8_errorf (
"Computation of bounding box failed. VTK merge points will be disabled. VTK output is still valid. \n");
}
}
/* Allocate VTK Memory for the arrays */
const int grid_dim = grid_get_dim (grid);
const int maximum_cell_size = curved_flag ? t8_curved_dim_max_nodes[grid_dim] : t8_dim_max_nodes[grid_dim];
cellArray->AllocateEstimate (num_cells, maximum_cell_size);
points_store->Allocate (num_cells * maximum_cell_size);
vtk_treeid->Allocate (num_cells);
vtk_mpirank->Allocate (num_cells);
vtk_level->Allocate (num_cells);
vtk_element_id->Allocate (num_cells);
/* Initialize vtkMergePoints for insertion. It needs the bounds to create the internal acceleration structure */
points->InitPointInsertion (points_store, bounds);
/* Iterate over all trees and translate them. */
const t8_locidx_t num_local_trees = grid_local_num_trees (grid);
for (t8_locidx_t itree = 0; itree < num_local_trees; itree++) {
t8_grid_tree_to_vtk_cells (grid, unstructuredGrid, vtk_treeid, vtk_mpirank, vtk_level, vtk_element_id, cellArray,
points, cellTypes, num_local_trees, &elem_id, &point_id, offset, false, itree);
}
if (do_ghosts) {
/* Iterate over all ghost trees and translate them. */
const t8_locidx_t num_ghost_trees = grid_local_num_ghost_trees (grid);
for (t8_locidx_t itree_ghost = 0; itree_ghost < num_ghost_trees; itree_ghost++) {
t8_grid_tree_to_vtk_cells (grid, unstructuredGrid, vtk_treeid, vtk_mpirank, vtk_level, vtk_element_id,
cellArray, points, cellTypes, num_local_trees, &elem_id, &point_id, offset, true,
itree_ghost);
}
}
/* Construct the unstructuredGrid. */
unstructuredGrid->SetPoints (points->GetPoints ());
unstructuredGrid->SetCells (cellTypes, cellArray);
if (this->write_treeid) {
vtk_treeid->SetName ("treeid");
unstructuredGrid->GetCellData ()->AddArray (vtk_treeid);
}
if (this->write_mpirank) {
vtk_mpirank->SetName ("mpirank");
unstructuredGrid->GetCellData ()->AddArray (vtk_mpirank);
}
if (this->write_level) {
vtk_level->SetName ("level");
unstructuredGrid->GetCellData ()->AddArray (vtk_level);
}
if (this->write_element_id) {
vtk_element_id->SetName ("element_id");
unstructuredGrid->GetCellData ()->AddArray (vtk_element_id);
}
/* Write the user defined data fields. For that we iterate over the idata, set the name, the array and then give
* this data to the unstructured Grid Object.We differentiate between scalar and vector data.
*/
for (int idata = 0; idata < num_data; idata++) {
dataArrays[idata] = vtkDoubleArray::New ();
const int num_components = this->data[idata].type == T8_VTK_SCALAR ? 1 : 3;
dataArrays[idata]->SetName (this->data[idata].description); /* Set the name of the array. */
dataArrays[idata]->SetNumberOfTuples (num_cells); /* We want number of tuples=number of cells. */
dataArrays[idata]->SetNumberOfComponents (num_components); /* Each tuples has 3 values. */
dataArrays[idata]->SetVoidArray (this->data[idata].data, num_cells * num_components, 1);
unstructuredGrid->GetCellData ()->AddArray (dataArrays[idata]);
}
/* We have to free the allocated memory for the cellTypes Array and the other arrays we allocated memory for. */
for (int idata = 0; idata < num_data; idata++) {
dataArrays[idata]->Delete ();
}
/* Release unused memory in the arrays */
unstructuredGrid->Squeeze ();
T8_FREE (cellTypes);
T8_FREE (dataArrays);
return;
}
#endif /* T8_ENABLE_VTK */
bool
write_vtk ([[maybe_unused]] const grid_t grid)
{
#if T8_ENABLE_VTK
T8_ASSERT (!fileprefix.empty ());
/*
* Write file: First we construct the unstructured Grid that will store the points and elements. It requires
* information about the points(coordinates, stored in the points object) and the cells(cellTypes and which points
* belong to this cell).
*/
vtkSmartPointer<vtkUnstructuredGrid> unstructuredGrid = vtkSmartPointer<vtkUnstructuredGrid>::New ();
std::string mpifilename = fileprefix + std::string (".pvtu");
vtkSmartPointer<vtkXMLPUnstructuredGridWriter> pwriterObj = vtkSmartPointer<vtkXMLPUnstructuredGridWriter>::New ();
t8_grid_to_vtkUnstructuredGrid (grid, unstructuredGrid);
/*
* Get/Set whether the appended data section is base64 encoded.
* If encoded, reading and writing will be slower, but the file
* will be fully valid XML and text-only.
* If not encoded, the XML specification will be violated,
* but reading and writing will be fast. The default is to do the encoding.
* Documentation: https://vtk.org/doc/release/5.0/html/a02260.html#z3560_2
*/
pwriterObj->EncodeAppendedDataOff ();
/* We set the filename of the pvtu file. The filenames of the vtu files
* are given based on the name of the pvtu file and the process number.
*/
pwriterObj->SetFileName (mpifilename.c_str ());
/*
* Since we want to write multiple files, the processes
* have to communicate. Therefore, we define the communicator
* vtk_comm and set it as the communicator.
* We have to set a controller for the pwriterObj,
* therefore we define the controller vtk_mpi_ctrl.
*/
#if T8_ENABLE_MPI
vtkSmartPointer<vtkMPICommunicator> vtk_comm = vtkSmartPointer<vtkMPICommunicator>::New ();
vtkMPICommunicatorOpaqueComm vtk_opaque_comm (&comm);
vtk_comm->InitializeExternal (&vtk_opaque_comm);
vtkSmartPointer<vtkMPIController> vtk_mpi_ctrl = vtkSmartPointer<vtkMPIController>::New ();
vtk_mpi_ctrl->SetCommunicator (vtk_comm);
pwriterObj->SetController (vtk_mpi_ctrl);
#endif /* T8_ENABLE_MPI */
/*
* We set the number of pieces as the number of mpi processes,
* since we want to write a file for each process. We also
* need to define a Start and EndPiece for the current
* process. Then we can set the inputData for the writer:
* We want to write the unstructured Grid, update the writer
* and then write.
*
* Note: We could write more than one file per process here, if desired.
*/
int mpisize;
int mpirank;
int mpiret = sc_MPI_Comm_size (comm, &mpisize);
SC_CHECK_MPI (mpiret);
mpiret = sc_MPI_Comm_rank (comm, &mpirank);
SC_CHECK_MPI (mpiret);
pwriterObj->SetNumberOfPieces (mpisize);
pwriterObj->SetStartPiece (mpirank);
pwriterObj->SetEndPiece (mpirank);
/* We set the input data and write the vtu files. */
pwriterObj->SetInputData (unstructuredGrid);
pwriterObj->Update ();
if (pwriterObj->Write ()) {
return true;
}
else {
t8_errorf ("Error when writing vtk file.\n");
}
/* Return whether writing was successful */
return false;
#else
t8_global_errorf ("Warning: t8code is not linked against vtk library. Vtk output will not be generated.\n");
t8_global_productionf ("Consider calling 't8_forest_write_vtk' or 't8_forest_vtk_write_file' instead.\n");
return false;
#endif
}
bool write_treeid = false;
bool write_mpirank = false;
bool write_level = false;
bool write_element_id = false;
bool write_ghosts = false;
bool curved_flag = false;
std::string fileprefix;
int num_data = 0;
t8_vtk_data_field_t *data = NULL;
sc_MPI_Comm comm;
bool merge_points = true;
};
#endif /* T8_VTK_WRITER_HXX */