utilities.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2023 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_vtk_utilities_h
#define dealii_vtk_utilities_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/point.h>
 
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/fe/fe.h>
 
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
 
#include <deal.II/lac/vector.h>
 
#ifdef DEAL_II_WITH_VTK
 
#  include <vtkDataObject.h>
#  include <vtkDoubleArray.h>
#  include <vtkSmartPointer.h>
#  include <vtkUnstructuredGrid.h>
 
#  include <string>
 
#endif
 
DEAL_II_NAMESPACE_OPEN
 
#ifdef DEAL_II_WITH_VTK
namespace VTKWrappers
{
  namespace internal
  {
    vtkSmartPointer<vtkUnstructuredGrid>
    load_vtk_file(const std::string &vtk_filename,
                  const bool         cleanup            = true,
                  const double       relative_tolerance = 0.0);
 
    vtkSmartPointer<vtkUnstructuredGrid>
    convert_to_unstructured_grid(
      const vtkSmartPointer<vtkDataObject> &data_object);
 
    void
    write_vtk(const std::string                    &vtk_filename,
              const vtkSmartPointer<vtkDataObject> &data_object);
  } // namespace internal
 
  template <int dim>
  inline vtkSmartPointer<vtkDoubleArray>
  dealii_point_to_vtk_array(const ::Point<dim> &p);
 
  template <int dim, int spacedim>
  void
  unstructured_grid_to_dealii_triangulation(
    const vtkUnstructuredGrid    &unstructured_grid,
    Triangulation<dim, spacedim> &tria,
    const std::string            &material_id_field = "",
    const std::string            &boundary_id_field = "",
    const std::string            &manifold_id_field = "");
 
  template <int dim, int spacedim>
  vtkSmartPointer<vtkUnstructuredGrid>
  dealii_triangulation_to_unstructured_grid(
    const Triangulation<dim, spacedim> &tria,
    const std::string                  &material_id_field = "",
    const std::string                  &boundary_id_field = "",
    const std::string                  &manifold_id_field = "");
 
  template <int dim, int spacedim>
  void
  write_vtk(const std::string                  &vtk_filename,
            const Triangulation<dim, spacedim> &tria,
            const std::string                  &material_id_field = "",
            const std::string                  &boundary_id_field = "",
            const std::string                  &manifold_id_field = "");
 
  template <int dim, int spacedim>
  void
  read_tria(const std::string            &vtk_filename,
            Triangulation<dim, spacedim> &tria,
            const bool                    cleanup            = true,
            const double                  relative_tolerance = 0.0,
            const std::string            &material_id_field  = "",
            const std::string            &boundary_id_field  = "",
            const std::string            &manifold_id_field  = "");
 
  void
  read_cell_data(const std::string &vtk_filename,
                 const std::string &cell_data_name,
                 Vector<double>    &output_vector,
                 const bool         cleanup            = true,
                 const double       relative_tolerance = 0.0);
 
  void
  read_vertex_data(const std::string &vtk_filename,
                   const std::string &vertex_data_name,
                   Vector<double>    &output_vector,
                   const bool         cleanup            = true,
                   const double       relative_tolerance = 0.0);
 
  void
  read_all_data(const std::string &vtk_filename,
                Vector<double>    &output_vector,
                const bool         cleanup            = true,
                const double       relative_tolerance = 0.0);
 
  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, 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,
           const bool                 cleanup            = true,
           const double               relative_tolerance = 0.0,
           const std::string         &material_id_field  = "",
           const std::string         &boundary_id_field  = "",
           const std::string         &manifold_id_field  = "");
 
#  ifndef DOXYGEN
  // Template implementations
 
  template <int dim>
  inline vtkSmartPointer<vtkDoubleArray>
  dealii_point_to_vtk_array(const ::Point<dim> &p)
  {
    vtkSmartPointer<vtkDoubleArray> p_vtk =
      vtkSmartPointer<vtkDoubleArray>::New();
 
    p_vtk->SetNumberOfComponents(dim);
    p_vtk->SetNumberOfTuples(1);
 
    for (int d = 0; d < dim; ++d)
      p_vtk->FillComponent(d, p[d]);
 
    return p_vtk;
  }
 
 
 
  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 =
      GridTools::parallel_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())
                {
                  // Advance 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;
          }
      }
  }
#  endif
 
} // namespace VTKWrappers
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
#endif