Program Listing for File vtk_utils.h¶
↰ Return to documentation for file (include/vtk_utils.h)
// ---------------------------------------------------------------------
//
// Copyright (C) 2024 by Luca Heltai
//
// This file is part of the reduced_lagrange_multipliers application, based on
// the deal.II library.
//
// The reduced_lagrange_multipliers application is free software; you can use
// it, redistribute it, and/or modify it under the terms of the Apache-2.0
// License WITH LLVM-exception as published by the Free Software Foundation;
// either version 3.0 of the License, or (at your option) any later version. The
// full text of the license can be found in the file LICENSE.md at the top level
// of the reduced_lagrange_multipliers distribution.
//
// ---------------------------------------------------------------------
#ifndef rdl_vtk_utils_h
#define rdl_vtk_utils_h
#include <deal.II/base/config.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/point.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/vector.h>
#include <iostream>
#include <map>
#include <unordered_map>
#ifdef DEAL_II_WITH_VTK
# include <deal.II/grid/tria.h>
# include <deal.II/lac/la_parallel_vector.h>
# include <deal.II/lac/vector.h>
using namespace dealii;
namespace VTKUtils
{
template <int dim, int spacedim>
void
read_vtk(const std::string &vtk_filename,
Triangulation<dim, spacedim> &tria,
const bool cleanup = true);
void
read_cell_data(const std::string &vtk_filename,
const std::string &cell_data_name,
Vector<double> &output_vector);
void
read_vertex_data(const std::string &vtk_filename,
const std::string &vertex_data_name,
Vector<double> &output_vector);
void
read_data(const std::string &vtk_filename, Vector<double> &output_vector);
template <int dim, int spacedim>
std::pair<std::unique_ptr<FiniteElement<dim, spacedim>>,
std::vector<std::string>>
vtk_to_finite_element(const std::string &vtk_filename);
template <int dim, int spacedim>
BlockIndices
get_block_indices(const FiniteElement<dim, spacedim> &fe);
template <int dim, int spacedim, typename VectorType>
void
data_to_dealii_vector(const Triangulation<dim, spacedim> &serial_tria,
const Vector<double> &data,
const DoFHandler<dim, spacedim> &dh,
VectorType &output_vector);
template <int dim, int spacedim>
void
read_vtk(const std::string &vtk_filename,
DoFHandler<dim, spacedim> &dof_handler,
Vector<double> &output_vector,
std::vector<std::string> &data_names);
template <int dim, int spacedim>
void
serial_vector_to_distributed_vector(
const DoFHandler<dim, spacedim> &serial_dh,
const DoFHandler<dim, spacedim> ¶llel_dh,
const Vector<double> &serial_vec,
LinearAlgebra::distributed::Vector<double> ¶llel_vec);
template <int dim, int spacedim>
std::vector<types::global_vertex_index>
distributed_to_serial_vertex_indices(
const Triangulation<dim, spacedim> &serial_tria,
const Triangulation<dim, spacedim> ¶llel_tria);
# ifndef DOXYGEN
// Explicit implementation of template functions
template <int dim, int spacedim, typename VectorType>
void
data_to_dealii_vector(const Triangulation<dim, spacedim> &serial_tria,
const Vector<double> &data,
const DoFHandler<dim, spacedim> &dh,
VectorType &output_vector)
{
AssertDimension(dh.n_dofs(), output_vector.size());
const auto &fe = dh.get_fe();
const auto dist_to_serial_vertices =
distributed_to_serial_vertex_indices(serial_tria, dh.get_triangulation());
const auto &locally_owned_dofs = dh.locally_owned_dofs();
types::global_dof_index dofs_offset = 0;
unsigned int vertex_comp_offset = 0;
unsigned int cell_comp_offset = 0;
for (unsigned int field = 0; field < fe.n_blocks(); ++field)
{
const auto &field_fe = fe.base_element(field);
const unsigned int n_comps = field_fe.n_components();
if (field_fe.n_dofs_per_vertex() > 0)
{
// This is a vertex data field
const types::global_dof_index n_local_dofs =
n_comps * serial_tria.n_vertices();
for (const auto &cell : dh.active_cell_iterators())
if (cell->is_locally_owned())
for (const auto v : cell->vertex_indices())
{
const types::global_dof_index serial_vertex_index =
dist_to_serial_vertices[cell->vertex_index(v)];
if (serial_vertex_index != numbers::invalid_unsigned_int)
for (unsigned int c = 0; c < n_comps; ++c)
{
const types::global_dof_index dof_index =
cell->vertex_dof_index(v, vertex_comp_offset + c);
Assert(locally_owned_dofs.is_element(dof_index),
ExcInternalError());
output_vector[dof_index] =
data[dofs_offset + n_comps * serial_vertex_index +
c];
}
}
dofs_offset += n_local_dofs;
vertex_comp_offset += n_comps;
}
else if (field_fe.template n_dofs_per_object<dim>() > 0)
{
// this is a cell data field
const types::global_dof_index n_local_dofs =
n_comps * serial_tria.n_global_active_cells();
// Assumption: serial and parallel meshes have the same ordering of
// cells.
auto serial_cell = serial_tria.begin_active();
auto parallel_cell = dh.begin_active();
for (; parallel_cell != dh.end(); ++parallel_cell)
if (parallel_cell->is_locally_owned())
{
// Advanced serial cell until we reach the same cell index of
// the parallel cell
while (serial_cell->id() < parallel_cell->id())
++serial_cell;
const auto serial_cell_index =
serial_cell->global_active_cell_index();
for (unsigned int c = 0; c < n_comps; ++c)
{
const types::global_dof_index dof_index =
parallel_cell->dof_index(cell_comp_offset + c);
Assert(locally_owned_dofs.is_element(dof_index),
ExcInternalError());
output_vector[dof_index] =
data[dofs_offset + n_comps * serial_cell_index + c];
}
}
dofs_offset += n_local_dofs;
cell_comp_offset += n_comps;
}
}
}
template <int dim, int spacedim>
BlockIndices
get_block_indices(const FiniteElement<dim, spacedim> &fe)
{
BlockIndices block_indices;
for (unsigned int i = 0; i < fe.n_blocks(); ++i)
{
const auto &block_fe = fe.base_element(i);
block_indices.push_back(block_fe.n_components());
}
return block_indices;
};
# endif
} // namespace VTKUtils
#endif // DEAL_II_WITH_VTK
#endif