vtk_utils.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// 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>            &parallel_dh,
    const Vector<double>                       &serial_vec,
    LinearAlgebra::distributed::Vector<double> &parallel_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> &parallel_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