grid_in.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 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 directory of deal.II.
//
// ---------------------------------------------------------------------


#include <deal.II/base/exceptions.h>
#include <deal.II/base/path_search.h>
#include <deal.II/base/utilities.h>

#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_reordering.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>

#include <boost/io/ios_state.hpp>

#include <algorithm>
#include <cctype>
#include <fstream>
#include <functional>
#include <map>


#ifdef DEAL_II_WITH_NETCDF
#  include <netcdfcpp.h>
#endif

#ifdef DEAL_II_WITH_ASSIMP
#  include <assimp/Importer.hpp>  // C++ importer interface
#  include <assimp/postprocess.h> // Post processing flags
#  include <assimp/scene.h>       // Output data structure
#endif


DEAL_II_NAMESPACE_OPEN


namespace
{
  template <int spacedim>
  void
  assign_1d_boundary_ids(
    const std::map<unsigned int, types::boundary_id> &boundary_ids,
    Triangulation<1, spacedim> &                      triangulation)
  {
    if (boundary_ids.size() > 0)
      for (typename Triangulation<1, spacedim>::active_cell_iterator cell =
             triangulation.begin_active();
           cell != triangulation.end();
           ++cell)
        for (unsigned int f = 0; f < GeometryInfo<1>::faces_per_cell; ++f)
          if (boundary_ids.find(cell->vertex_index(f)) != boundary_ids.end())
            {
              AssertThrow(
                cell->at_boundary(f),
                ExcMessage(
                  "You are trying to prescribe boundary ids on the face "
                  "of a 1d cell (i.e., on a vertex), but this face is not actually at "
                  "the boundary of the mesh. This is not allowed."));
              cell->face(f)->set_boundary_id(
                boundary_ids.find(cell->vertex_index(f))->second);
            }
  }


  template <int dim, int spacedim>
  void
  assign_1d_boundary_ids(const std::map<unsigned int, types::boundary_id> &,
                         Triangulation<dim, spacedim> &)
  {
    // we shouldn't get here since boundary ids are not assigned to
    // vertices except in 1d
    Assert(dim != 1, ExcInternalError());
  }
} // namespace

template <int dim, int spacedim>
GridIn<dim, spacedim>::GridIn()
  : tria(nullptr, typeid(*this).name())
  , default_format(ucd)
{}


template <int dim, int spacedim>
void
GridIn<dim, spacedim>::attach_triangulation(Triangulation<dim, spacedim> &t)
{
  tria = &t;
}



template <int dim, int spacedim>
void
GridIn<dim, spacedim>::read_vtk(std::istream &in)
{
  Assert((dim == 2) || (dim == 3), ExcNotImplemented());
  std::string line;

  // verify that the first, third and fourth lines match
  // expectations. the second line of the file may essentially be
  // anything the author of the file chose to identify what's in
  // there, so we just ensure that we can read it
  {
    std::string text[4];
    text[0] = "# vtk DataFile Version 3.0";
    text[1] = "****";
    text[2] = "ASCII";
    text[3] = "DATASET UNSTRUCTURED_GRID";

    for (unsigned int i = 0; i < 4; ++i)
      {
        getline(in, line);
        if (i != 1)
          AssertThrow(
            line.compare(text[i]) == 0,
            ExcMessage(
              std::string(
                "While reading VTK file, failed to find a header line with text <") +
              text[i] + ">"));
      }
  }


  std::vector<Point<spacedim>> vertices;
  std::vector<CellData<dim>>   cells;
  SubCellData                  subcelldata;

  std::string keyword;

  in >> keyword;


  if (keyword == "POINTS")
    {
      unsigned int n_vertices;
      in >> n_vertices;

      in.ignore(256,
                '\n'); // ignoring the number beside the total no. of points.

      for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
        {
          // VTK format always specifies vertex coordinates with 3 components
          Point<3> x;
          in >> x(0) >> x(1) >> x(2);

          vertices.emplace_back();
          for (unsigned int d = 0; d < spacedim; ++d)
            vertices.back()(d) = x(d);
        }
    }

  else
    AssertThrow(false,
                ExcMessage(
                  "While reading VTK file, failed to find POINTS section"));


  std::string checkline;
  int         no;
  in.ignore(256,
            '\n'); // this move pointer to the next line ignoring unwanted no.
  no = in.tellg();
  getline(in, checkline);
  if (checkline.compare("") != 0)
    {
      in.seekg(no);
    }

  in >> keyword;


  if (keyword == "CELLS")
    {
      unsigned int n_geometric_objects;
      in >> n_geometric_objects;
      in.ignore(256, '\n');

      if (dim == 3)
        {
          for (unsigned int count = 0; count < n_geometric_objects; count++)
            {
              unsigned int type;
              in >> type;

              if (type == 8)
                {
                  // we assume that the file contains first all cells,
                  // and only then any faces or lines
                  AssertThrow(subcelldata.boundary_quads.size() == 0 &&
                                subcelldata.boundary_lines.size() == 0,
                              ExcNotImplemented());

                  cells.emplace_back();

                  for (unsigned int j = 0; j < type; j++) // loop to feed data
                    in >> cells.back().vertices[j];

                  cells.back().material_id = 0;
                }

              else if (type == 4)
                {
                  subcelldata.boundary_quads.emplace_back();

                  for (unsigned int j = 0; j < type;
                       j++) // loop to feed the data to the boundary
                    in >> subcelldata.boundary_quads.back().vertices[j];

                  subcelldata.boundary_quads.back().material_id = 0;
                }

              else
                AssertThrow(
                  false,
                  ExcMessage(
                    "While reading VTK file, unknown file type encountered"));
            }
        }

      else if (dim == 2)
        {
          for (unsigned int count = 0; count < n_geometric_objects; count++)
            {
              unsigned int type;
              in >> type;

              if (type == 4)
                {
                  // we assume that the file contains first all cells,
                  // and only then any faces
                  AssertThrow(subcelldata.boundary_lines.size() == 0,
                              ExcNotImplemented());

                  cells.emplace_back();

                  for (unsigned int j = 0; j < type; j++) // loop to feed data
                    in >> cells.back().vertices[j];

                  cells.back().material_id = 0;
                }

              else if (type == 2)
                {
                  // If this is encountered, the pointer comes out of the loop
                  // and starts processing boundaries.
                  subcelldata.boundary_lines.emplace_back();

                  for (unsigned int j = 0; j < type;
                       j++) // loop to feed the data to the boundary
                    {
                      in >> subcelldata.boundary_lines.back().vertices[j];
                    }

                  subcelldata.boundary_lines.back().material_id = 0;
                }

              else
                AssertThrow(
                  false,
                  ExcMessage(
                    "While reading VTK file, unknown file type encountered"));
            }
        }
      else
        AssertThrow(false,
                    ExcMessage(
                      "While reading VTK file, failed to find CELLS section"));


      in >> keyword;

      if (keyword ==
          "CELL_TYPES") // Entering the cell_types section and ignoring data.
        {
          in.ignore(256, '\n');

          while (!in.fail() && !in.eof())
            {
              in >> keyword;
              if (keyword != "12" && keyword != "9")
                {
                  break;
                }
            }
        }


      if (keyword == "CELL_DATA")
        {
          unsigned int n_ids;
          in >> n_ids;

          AssertThrow(
            n_ids == cells.size() +
                       (dim == 3 ?
                          subcelldata.boundary_quads.size() :
                          (dim == 2 ? subcelldata.boundary_lines.size() : 0)),
            ExcMessage(
              "The VTK reader found a CELL_DATA statement "
              "that lists a total of " +
              Utilities::int_to_string(n_ids) +
              " cell data objects, but this needs to "
              "equal the number of cells (which is " +
              Utilities::int_to_string(cells.size()) +
              ") plus the number of quads (" +
              Utilities::int_to_string(subcelldata.boundary_quads.size()) +
              " in 3d or the number of lines (" +
              Utilities::int_to_string(subcelldata.boundary_lines.size()) +
              ") in 2d."));


          std::string linenew;
          std::string textnew[2];
          textnew[0] = "SCALARS MaterialID double";
          textnew[1] = "LOOKUP_TABLE default";

          in.ignore(256, '\n');

          for (unsigned int i = 0; i < 2; i++)
            {
              getline(in, linenew);
              if (i == 0)
                if (linenew.size() > textnew[0].size())
                  linenew.resize(textnew[0].size());

              AssertThrow(linenew.compare(textnew[i]) == 0,
                          ExcMessage(
                            std::string(
                              "While reading VTK file, failed to find <") +
                            textnew[i] + "> section"));
            }

          // read material ids first for all cells, then for all
          // faces. the assumption that cells come before all faces
          // has been verified above via an assertion, so the order
          // used in the following blocks makes sense
          for (unsigned int i = 0; i < cells.size(); i++)
            {
              double id;
              in >> id;
              cells[i].material_id = id;
            }

          if (dim == 3)
            {
              for (unsigned int i = 0; i < subcelldata.boundary_quads.size();
                   i++)
                {
                  double id;
                  in >> id;
                  subcelldata.boundary_quads[i].material_id = id;
                }
            }
          else if (dim == 2)
            {
              for (unsigned int i = 0; i < subcelldata.boundary_lines.size();
                   i++)
                {
                  double id;
                  in >> id;
                  subcelldata.boundary_lines[i].material_id = id;
                }
            }
        }

      Assert(subcelldata.check_consistency(dim), ExcInternalError());

      GridTools::delete_unused_vertices(vertices, cells, subcelldata);

      if (dim == spacedim)
        GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(
          vertices, cells);

      GridReordering<dim, spacedim>::reorder_cells(cells);
      tria->create_triangulation_compatibility(vertices, cells, subcelldata);

      return;
    }
  else
    AssertThrow(false,
                ExcMessage(
                  "While reading VTK file, failed to find CELLS section"));
}



template <int dim, int spacedim>
void
GridIn<dim, spacedim>::read_unv(std::istream &in)
{
  Assert(tria != nullptr, ExcNoTriangulationSelected());
  Assert((dim == 2) || (dim == 3), ExcNotImplemented());

  AssertThrow(in, ExcIO());
  skip_comment_lines(in, '#'); // skip comments (if any) at beginning of file

  int tmp;

  AssertThrow(in, ExcIO());
  in >> tmp;
  AssertThrow(in, ExcIO());
  in >> tmp;

  // section 2411 describes vertices: see
  // http://www.sdrl.uc.edu/sdrl/referenceinfo/universalfileformats/file-format-storehouse/universal-dataset-number-2411
  AssertThrow(tmp == 2411, ExcUnknownSectionType(tmp));

  std::vector<Point<spacedim>> vertices; // vector of vertex coordinates
  std::map<int, int>
    vertex_indices; // # vert in unv (key) ---> # vert in deal.II (value)

  int no_vertex = 0; // deal.II

  while (tmp != -1) // we do until reach end of 2411
    {
      int    no; // unv
      int    dummy;
      double x[3];

      AssertThrow(in, ExcIO());
      in >> no;

      tmp = no;
      if (tmp == -1)
        break;

      in >> dummy >> dummy >> dummy;

      AssertThrow(in, ExcIO());
      in >> x[0] >> x[1] >> x[2];

      vertices.emplace_back();

      for (unsigned int d = 0; d < spacedim; d++)
        vertices.back()(d) = x[d];

      vertex_indices[no] = no_vertex;

      no_vertex++;
    }

  AssertThrow(in, ExcIO());
  in >> tmp;
  AssertThrow(in, ExcIO());
  in >> tmp;

  // section 2412 describes elements: see
  // http://www.sdrl.uc.edu/sdrl/referenceinfo/universalfileformats/file-format-storehouse/universal-dataset-number-2412
  AssertThrow(tmp == 2412, ExcUnknownSectionType(tmp));

  std::vector<CellData<dim>> cells; // vector of cells
  SubCellData                subcelldata;

  std::map<int, int>
    cell_indices; // # cell in unv (key) ---> # cell in deal.II (value)
  std::map<int, int>
    line_indices; // # line in unv (key) ---> # line in deal.II (value)
  std::map<int, int>
    quad_indices; // # quad in unv (key) ---> # quad in deal.II (value)

  int no_cell = 0; // deal.II
  int no_line = 0; // deal.II
  int no_quad = 0; // deal.II

  while (tmp != -1) // we do until reach end of 2412
    {
      int no; // unv
      int type;
      int dummy;

      AssertThrow(in, ExcIO());
      in >> no;

      tmp = no;
      if (tmp == -1)
        break;

      in >> type >> dummy >> dummy >> dummy >> dummy;

      AssertThrow((type == 11) || (type == 44) || (type == 94) || (type == 115),
                  ExcUnknownElementType(type));

      if ((((type == 44) || (type == 94)) && (dim == 2)) ||
          ((type == 115) && (dim == 3))) // cell
        {
          cells.emplace_back();

          AssertThrow(in, ExcIO());
          for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
               v++)
            in >> cells.back().vertices[v];

          cells.back().material_id = 0;

          for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
               v++)
            cells.back().vertices[v] = vertex_indices[cells.back().vertices[v]];

          cell_indices[no] = no_cell;

          no_cell++;
        }
      else if (((type == 11) && (dim == 2)) ||
               ((type == 11) && (dim == 3))) // boundary line
        {
          AssertThrow(in, ExcIO());
          in >> dummy >> dummy >> dummy;

          subcelldata.boundary_lines.emplace_back();

          AssertThrow(in, ExcIO());
          for (unsigned int v = 0; v < 2; v++)
            in >> subcelldata.boundary_lines.back().vertices[v];

          subcelldata.boundary_lines.back().material_id = 0;

          for (unsigned int v = 0; v < 2; v++)
            subcelldata.boundary_lines.back().vertices[v] =
              vertex_indices[subcelldata.boundary_lines.back().vertices[v]];

          line_indices[no] = no_line;

          no_line++;
        }
      else if (((type == 44) || (type == 94)) && (dim == 3)) // boundary quad
        {
          subcelldata.boundary_quads.emplace_back();

          AssertThrow(in, ExcIO());
          for (unsigned int v = 0; v < 4; v++)
            in >> subcelldata.boundary_quads.back().vertices[v];

          subcelldata.boundary_quads.back().material_id = 0;

          for (unsigned int v = 0; v < 4; v++)
            subcelldata.boundary_quads.back().vertices[v] =
              vertex_indices[subcelldata.boundary_quads.back().vertices[v]];

          quad_indices[no] = no_quad;

          no_quad++;
        }
      else
        AssertThrow(false,
                    ExcMessage("Unknown element label <" +
                               Utilities::int_to_string(type) +
                               "> when running in dim=" +
                               Utilities::int_to_string(dim)));
    }

  // note that so far all materials and bcs are explicitly set to 0
  // if we do not need more info on materials and bcs - this is end of file
  // if we do - section 2467 or 2477 comes

  in >> tmp; // tmp can be either -1 or end-of-file

  if (!in.eof())
    {
      AssertThrow(in, ExcIO());
      in >> tmp;

      // section 2467 (2477) describes (materials - first and bcs - second) or
      // (bcs - first and materials - second) - sequence depends on which
      // group is created first: see
      // http://www.sdrl.uc.edu/sdrl/referenceinfo/universalfileformats/file-format-storehouse/universal-dataset-number-2467
      AssertThrow((tmp == 2467) || (tmp == 2477), ExcUnknownSectionType(tmp));

      while (tmp != -1) // we do until reach end of 2467 or 2477
        {
          int n_entities; // number of entities in group
          int id;         // id is either material or bc
          int no;         // unv
          int dummy;

          AssertThrow(in, ExcIO());
          in >> dummy;

          tmp = dummy;
          if (tmp == -1)
            break;

          in >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >>
            n_entities;

          AssertThrow(in, ExcIO());
          in >> id;

          const unsigned int n_lines =
            (n_entities % 2 == 0) ? (n_entities / 2) : ((n_entities + 1) / 2);

          for (unsigned int line = 0; line < n_lines; line++)
            {
              unsigned int n_fragments;

              if (line == n_lines - 1)
                n_fragments = (n_entities % 2 == 0) ? (2) : (1);
              else
                n_fragments = 2;

              for (unsigned int no_fragment = 0; no_fragment < n_fragments;
                   no_fragment++)
                {
                  AssertThrow(in, ExcIO());
                  in >> dummy >> no >> dummy >> dummy;

                  if (cell_indices.count(no) > 0) // cell - material
                    cells[cell_indices[no]].material_id = id;

                  if (line_indices.count(no) > 0) // boundary line - bc
                    subcelldata.boundary_lines[line_indices[no]].material_id =
                      id;

                  if (quad_indices.count(no) > 0) // boundary quad - bc
                    subcelldata.boundary_quads[quad_indices[no]].material_id =
                      id;
                }
            }
        }
    }

  Assert(subcelldata.check_consistency(dim), ExcInternalError());

  GridTools::delete_unused_vertices(vertices, cells, subcelldata);

  if (dim == spacedim)
    GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(vertices,
                                                                     cells);

  GridReordering<dim, spacedim>::reorder_cells(cells);

  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
}



template <int dim, int spacedim>
void
GridIn<dim, spacedim>::read_ucd(std::istream &in,
                                const bool    apply_all_indicators_to_manifolds)
{
  Assert(tria != nullptr, ExcNoTriangulationSelected());
  AssertThrow(in, ExcIO());

  // skip comments at start of file
  skip_comment_lines(in, '#');


  unsigned int n_vertices;
  unsigned int n_cells;
  int          dummy;

  in >> n_vertices >> n_cells >> dummy // number of data vectors
    >> dummy                           // cell data
    >> dummy;                          // model data
  AssertThrow(in, ExcIO());

  // set up array of vertices
  std::vector<Point<spacedim>> vertices(n_vertices);
  // set up mapping between numbering
  // in ucd-file (key) and in the
  // vertices vector
  std::map<int, int> vertex_indices;

  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
    {
      int    vertex_number;
      double x[3];

      // read vertex
      AssertThrow(in, ExcIO());
      in >> vertex_number >> x[0] >> x[1] >> x[2];

      // store vertex
      for (unsigned int d = 0; d < spacedim; ++d)
        vertices[vertex](d) = x[d];
      // store mapping; note that
      // vertices_indices[i] is automatically
      // created upon first usage
      vertex_indices[vertex_number] = vertex;
    };

  // set up array of cells
  std::vector<CellData<dim>> cells;
  SubCellData                subcelldata;

  for (unsigned int cell = 0; cell < n_cells; ++cell)
    {
      // note that since in the input
      // file we found the number of
      // cells at the top, there
      // should still be input here,
      // so check this:
      AssertThrow(in, ExcIO());

      std::string cell_type;

      // we use an unsigned int because we
      // fill this variable through an read-in process
      unsigned int material_id;

      in >> dummy // cell number
        >> material_id;
      in >> cell_type;

      if (((cell_type == "line") && (dim == 1)) ||
          ((cell_type == "quad") && (dim == 2)) ||
          ((cell_type == "hex") && (dim == 3)))
        // found a cell
        {
          // allocate and read indices
          cells.emplace_back();
          for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell;
               ++i)
            in >> cells.back().vertices[i];

          // to make sure that the cast won't fail
          Assert(material_id <= std::numeric_limits<types::material_id>::max(),
                 ExcIndexRange(material_id,
                               0,
                               std::numeric_limits<types::material_id>::max()));
          // we use only material_ids in the range from 0 to
          // numbers::invalid_material_id-1
          Assert(material_id < numbers::invalid_material_id,
                 ExcIndexRange(material_id, 0, numbers::invalid_material_id));

          cells.back().material_id =
            static_cast<types::material_id>(material_id);

          // transform from ucd to
          // consecutive numbering
          for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell;
               ++i)
            if (vertex_indices.find(cells.back().vertices[i]) !=
                vertex_indices.end())
              // vertex with this index exists
              cells.back().vertices[i] =
                vertex_indices[cells.back().vertices[i]];
            else
              {
                // no such vertex index
                AssertThrow(false,
                            ExcInvalidVertexIndex(cell,
                                                  cells.back().vertices[i]));

                cells.back().vertices[i] = numbers::invalid_unsigned_int;
              }
        }
      else if ((cell_type == "line") && ((dim == 2) || (dim == 3)))
        // boundary info
        {
          subcelldata.boundary_lines.emplace_back();
          in >> subcelldata.boundary_lines.back().vertices[0] >>
            subcelldata.boundary_lines.back().vertices[1];

          // to make sure that the cast won't fail
          Assert(material_id <= std::numeric_limits<types::boundary_id>::max(),
                 ExcIndexRange(material_id,
                               0,
                               std::numeric_limits<types::boundary_id>::max()));
          // we use only boundary_ids in the range from 0 to
          // numbers::internal_face_boundary_id-1
          Assert(material_id < numbers::internal_face_boundary_id,
                 ExcIndexRange(material_id,
                               0,
                               numbers::internal_face_boundary_id));

          if (apply_all_indicators_to_manifolds)
            subcelldata.boundary_lines.back().manifold_id =
              static_cast<types::manifold_id>(material_id);
          else
            subcelldata.boundary_lines.back().boundary_id =
              static_cast<types::boundary_id>(material_id);

          // transform from ucd to
          // consecutive numbering
          for (unsigned int i = 0; i < 2; ++i)
            if (vertex_indices.find(
                  subcelldata.boundary_lines.back().vertices[i]) !=
                vertex_indices.end())
              // vertex with this index exists
              subcelldata.boundary_lines.back().vertices[i] =
                vertex_indices[subcelldata.boundary_lines.back().vertices[i]];
            else
              {
                // no such vertex index
                AssertThrow(false,
                            ExcInvalidVertexIndex(
                              cell,
                              subcelldata.boundary_lines.back().vertices[i]));
                subcelldata.boundary_lines.back().vertices[i] =
                  numbers::invalid_unsigned_int;
              };
        }
      else if ((cell_type == "quad") && (dim == 3))
        // boundary info
        {
          subcelldata.boundary_quads.emplace_back();
          in >> subcelldata.boundary_quads.back().vertices[0] >>
            subcelldata.boundary_quads.back().vertices[1] >>
            subcelldata.boundary_quads.back().vertices[2] >>
            subcelldata.boundary_quads.back().vertices[3];

          // to make sure that the cast won't fail
          Assert(material_id <= std::numeric_limits<types::boundary_id>::max(),
                 ExcIndexRange(material_id,
                               0,
                               std::numeric_limits<types::boundary_id>::max()));
          // we use only boundary_ids in the range from 0 to
          // numbers::internal_face_boundary_id-1
          Assert(material_id < numbers::internal_face_boundary_id,
                 ExcIndexRange(material_id,
                               0,
                               numbers::internal_face_boundary_id));

          if (apply_all_indicators_to_manifolds)
            subcelldata.boundary_quads.back().manifold_id =
              static_cast<types::manifold_id>(material_id);
          else
            subcelldata.boundary_quads.back().boundary_id =
              static_cast<types::boundary_id>(material_id);

          // transform from ucd to
          // consecutive numbering
          for (unsigned int i = 0; i < 4; ++i)
            if (vertex_indices.find(
                  subcelldata.boundary_quads.back().vertices[i]) !=
                vertex_indices.end())
              // vertex with this index exists
              subcelldata.boundary_quads.back().vertices[i] =
                vertex_indices[subcelldata.boundary_quads.back().vertices[i]];
            else
              {
                // no such vertex index
                Assert(false,
                       ExcInvalidVertexIndex(
                         cell, subcelldata.boundary_quads.back().vertices[i]));
                subcelldata.boundary_quads.back().vertices[i] =
                  numbers::invalid_unsigned_int;
              };
        }
      else
        // cannot read this
        AssertThrow(false, ExcUnknownIdentifier(cell_type));
    };


  // check that no forbidden arrays are used
  Assert(subcelldata.check_consistency(dim), ExcInternalError());

  AssertThrow(in, ExcIO());

  // do some clean-up on vertices...
  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
  // ... and cells
  if (dim == spacedim)
    GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(vertices,
                                                                     cells);
  GridReordering<dim, spacedim>::reorder_cells(cells);
  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
}

namespace
{
  template <int dim>
  class Abaqus_to_UCD
  {
  public:
    Abaqus_to_UCD();

    void
    read_in_abaqus(std::istream &in);
    void
    write_out_avs_ucd(std::ostream &out) const;

  private:
    const double tolerance;

    std::vector<double>
    get_global_node_numbers(const int face_cell_no,
                            const int face_cell_face_no) const;

    // NL: Stored as [ global node-id (int), x-coord, y-coord, z-coord ]
    std::vector<std::vector<double>> node_list;
    // CL: Stored as [ material-id (int), node1, node2, node3, node4, node5,
    // node6, node7, node8 ]
    std::vector<std::vector<double>> cell_list;
    // FL: Stored as [ sideset-id (int), node1, node2, node3, node4 ]
    std::vector<std::vector<double>> face_list;
    // ELSET: Stored as [ (std::string) elset_name = (std::vector) of cells
    // numbers]
    std::map<std::string, std::vector<int>> elsets_list;
  };
} // namespace

template <int dim, int spacedim>
void
GridIn<dim, spacedim>::read_abaqus(std::istream &in,
                                   const bool apply_all_indicators_to_manifolds)
{
  Assert(tria != nullptr, ExcNoTriangulationSelected());
  Assert(dim == 2 || dim == 3, ExcNotImplemented());
  AssertThrow(in, ExcIO());

  // Read in the Abaqus file into an intermediate object
  // that is to be passed along to the UCD reader
  Abaqus_to_UCD<dim> abaqus_to_ucd;
  abaqus_to_ucd.read_in_abaqus(in);

  std::stringstream in_ucd;
  abaqus_to_ucd.write_out_avs_ucd(in_ucd);

  // This next call is wrapped in a try-catch for the following reason:
  // It ensures that if the Abaqus mesh is read in correctly but produces
  // an erroneous result then the user is alerted to the source of the problem
  // and doesn't think that they've somehow called the wrong function.
  try
    {
      read_ucd(in_ucd, apply_all_indicators_to_manifolds);
    }
  catch (std::exception &exc)
    {
      std::cerr << "Exception on processing internal UCD data: " << std::endl
                << exc.what() << std::endl;

      AssertThrow(
        false,
        ExcMessage(
          "Internal conversion from ABAQUS file to UCD format was unsuccessful. "
          "More information is provided in an error message printed above. "
          "Are you sure that your ABAQUS mesh file conforms with the requirements "
          "listed in the documentation?"));
    }
  catch (...)
    {
      AssertThrow(
        false,
        ExcMessage(
          "Internal conversion from ABAQUS file to UCD format was unsuccessful. "
          "Are you sure that your ABAQUS mesh file conforms with the requirements "
          "listed in the documentation?"));
    }
}


template <int dim, int spacedim>
void
GridIn<dim, spacedim>::read_dbmesh(std::istream &in)
{
  Assert(tria != nullptr, ExcNoTriangulationSelected());
  Assert(dim == 2, ExcNotImplemented());

  AssertThrow(in, ExcIO());

  // skip comments at start of file
  skip_comment_lines(in, '#');

  // first read in identifier string
  std::string line;
  getline(in, line);

  AssertThrow(line == "MeshVersionFormatted 0", ExcInvalidDBMESHInput(line));

  skip_empty_lines(in);

  // next read dimension
  getline(in, line);
  AssertThrow(line == "Dimension", ExcInvalidDBMESHInput(line));
  unsigned int dimension;
  in >> dimension;
  AssertThrow(dimension == dim, ExcDBMESHWrongDimension(dimension));
  skip_empty_lines(in);

  // now there are a lot of fields of
  // which we don't know the exact
  // meaning and which are far from
  // being properly documented in the
  // manual. we skip everything until
  // we find a comment line with the
  // string "# END". at some point in
  // the future, someone may have the
  // knowledge to parse and interpret
  // the other fields in between as
  // well...
  while (getline(in, line), line.find("# END") == std::string::npos)
    ;
  skip_empty_lines(in);


  // now read vertices
  getline(in, line);
  AssertThrow(line == "Vertices", ExcInvalidDBMESHInput(line));

  unsigned int n_vertices;
  double       dummy;

  in >> n_vertices;
  std::vector<Point<spacedim>> vertices(n_vertices);
  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
    {
      // read vertex coordinates
      for (unsigned int d = 0; d < dim; ++d)
        in >> vertices[vertex][d];
      // read Ref phi_i, whatever that may be
      in >> dummy;
    };
  AssertThrow(in, ExcInvalidDBMeshFormat());

  skip_empty_lines(in);

  // read edges. we ignore them at
  // present, so just read them and
  // discard the input
  getline(in, line);
  AssertThrow(line == "Edges", ExcInvalidDBMESHInput(line));

  unsigned int n_edges;
  in >> n_edges;
  for (unsigned int edge = 0; edge < n_edges; ++edge)
    {
      // read vertex indices
      in >> dummy >> dummy;
      // read Ref phi_i, whatever that may be
      in >> dummy;
    };
  AssertThrow(in, ExcInvalidDBMeshFormat());

  skip_empty_lines(in);



  // read cracked edges (whatever
  // that may be). we ignore them at
  // present, so just read them and
  // discard the input
  getline(in, line);
  AssertThrow(line == "CrackedEdges", ExcInvalidDBMESHInput(line));

  in >> n_edges;
  for (unsigned int edge = 0; edge < n_edges; ++edge)
    {
      // read vertex indices
      in >> dummy >> dummy;
      // read Ref phi_i, whatever that may be
      in >> dummy;
    };
  AssertThrow(in, ExcInvalidDBMeshFormat());

  skip_empty_lines(in);


  // now read cells.
  // set up array of cells
  getline(in, line);
  AssertThrow(line == "Quadrilaterals", ExcInvalidDBMESHInput(line));

  std::vector<CellData<dim>> cells;
  SubCellData                subcelldata;
  unsigned int               n_cells;
  in >> n_cells;
  for (unsigned int cell = 0; cell < n_cells; ++cell)
    {
      // read in vertex numbers. they
      // are 1-based, so subtract one
      cells.emplace_back();
      for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
        {
          in >> cells.back().vertices[i];

          AssertThrow((cells.back().vertices[i] >= 1) &&
                        (static_cast<unsigned int>(cells.back().vertices[i]) <=
                         vertices.size()),
                      ExcInvalidVertexIndex(cell, cells.back().vertices[i]));

          --cells.back().vertices[i];
        };

      // read and discard Ref phi_i
      in >> dummy;
    };
  AssertThrow(in, ExcInvalidDBMeshFormat());

  skip_empty_lines(in);


  // then there are again a whole lot
  // of fields of which I have no
  // clue what they mean. skip them
  // all and leave the interpretation
  // to other implementors...
  while (getline(in, line), ((line.find("End") == std::string::npos) && (in)))
    ;
  // ok, so we are not at the end of
  // the file, that's it, mostly


  // check that no forbidden arrays are used
  Assert(subcelldata.check_consistency(dim), ExcInternalError());

  AssertThrow(in, ExcIO());

  // do some clean-up on vertices...
  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
  // ...and cells
  GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(vertices,
                                                                   cells);
  GridReordering<dim, spacedim>::reorder_cells(cells);
  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
}



template <int dim, int spacedim>
void
GridIn<dim, spacedim>::read_xda(std::istream &)
{
  Assert(false, ExcNotImplemented());
}



template <>
void
GridIn<2>::read_xda(std::istream &in)
{
  Assert(tria != nullptr, ExcNoTriangulationSelected());
  AssertThrow(in, ExcIO());

  std::string line;
  // skip comments at start of file
  getline(in, line);


  unsigned int n_vertices;
  unsigned int n_cells;

  // read cells, throw away rest of line
  in >> n_cells;
  getline(in, line);

  in >> n_vertices;
  getline(in, line);

  // ignore following 8 lines
  for (unsigned int i = 0; i < 8; ++i)
    getline(in, line);

  // set up array of cells
  std::vector<CellData<2>> cells(n_cells);
  SubCellData              subcelldata;

  for (unsigned int cell = 0; cell < n_cells; ++cell)
    {
      // note that since in the input
      // file we found the number of
      // cells at the top, there
      // should still be input here,
      // so check this:
      AssertThrow(in, ExcIO());
      Assert(GeometryInfo<2>::vertices_per_cell == 4, ExcInternalError());

      for (unsigned int i = 0; i < 4; ++i)
        in >> cells[cell].vertices[i];
    };



  // set up array of vertices
  std::vector<Point<2>> vertices(n_vertices);
  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
    {
      double x[3];

      // read vertex
      in >> x[0] >> x[1] >> x[2];

      // store vertex
      for (unsigned int d = 0; d < 2; ++d)
        vertices[vertex](d) = x[d];
    };
  AssertThrow(in, ExcIO());

  // do some clean-up on vertices...
  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
  // ... and cells
  GridReordering<2>::invert_all_cells_of_negative_grid(vertices, cells);
  GridReordering<2>::reorder_cells(cells);
  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
}



template <>
void
GridIn<3>::read_xda(std::istream &in)
{
  Assert(tria != nullptr, ExcNoTriangulationSelected());
  AssertThrow(in, ExcIO());

  static const unsigned int xda_to_dealII_map[] = {0, 1, 5, 4, 3, 2, 6, 7};

  std::string line;
  // skip comments at start of file
  getline(in, line);


  unsigned int n_vertices;
  unsigned int n_cells;

  // read cells, throw away rest of line
  in >> n_cells;
  getline(in, line);

  in >> n_vertices;
  getline(in, line);

  // ignore following 8 lines
  for (unsigned int i = 0; i < 8; ++i)
    getline(in, line);

  // set up array of cells
  std::vector<CellData<3>> cells(n_cells);
  SubCellData              subcelldata;

  for (unsigned int cell = 0; cell < n_cells; ++cell)
    {
      // note that since in the input
      // file we found the number of
      // cells at the top, there
      // should still be input here,
      // so check this:
      AssertThrow(in, ExcIO());
      Assert(GeometryInfo<3>::vertices_per_cell == 8, ExcInternalError());

      unsigned int xda_ordered_nodes[8];

      for (unsigned int i = 0; i < 8; ++i)
        in >> xda_ordered_nodes[i];

      for (unsigned int i = 0; i < 8; i++)
        cells[cell].vertices[i] = xda_ordered_nodes[xda_to_dealII_map[i]];
    };



  // set up array of vertices
  std::vector<Point<3>> vertices(n_vertices);
  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
    {
      double x[3];

      // read vertex
      in >> x[0] >> x[1] >> x[2];

      // store vertex
      for (unsigned int d = 0; d < 3; ++d)
        vertices[vertex](d) = x[d];
    };
  AssertThrow(in, ExcIO());

  // do some clean-up on vertices...
  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
  // ... and cells
  GridReordering<3>::invert_all_cells_of_negative_grid(vertices, cells);
  GridReordering<3>::reorder_cells(cells);
  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
}



template <int dim, int spacedim>
void
GridIn<dim, spacedim>::read_msh(std::istream &in)
{
  Assert(tria != nullptr, ExcNoTriangulationSelected());
  AssertThrow(in, ExcIO());

  unsigned int n_vertices;
  unsigned int n_cells;
  unsigned int dummy;
  std::string  line;

  in >> line;

  // first determine file format
  unsigned int gmsh_file_format = 0;
  if (line == "@f$NOD")
    gmsh_file_format = 1;
  else if (line == "@f$MeshFormat")
    gmsh_file_format = 2;
  else
    AssertThrow(false, ExcInvalidGMSHInput(line));

  // if file format is 2 or greater
  // then we also have to read the
  // rest of the header
  if (gmsh_file_format == 2)
    {
      double       version;
      unsigned int file_type, data_size;

      in >> version >> file_type >> data_size;

      Assert((version >= 2.0) && (version <= 2.2), ExcNotImplemented());
      Assert(file_type == 0, ExcNotImplemented());
      Assert(data_size == sizeof(double), ExcNotImplemented());

      // read the end of the header
      // and the first line of the
      // nodes description to synch
      // ourselves with the format 1
      // handling above
      in >> line;
      AssertThrow(line == "@f$EndMeshFormat", ExcInvalidGMSHInput(line));

      in >> line;
      // if the next block is of kind
      // @f$PhysicalNames, ignore it
      if (line == "@f$PhysicalNames")
        {
          do
            {
              in >> line;
            }
          while (line != "@f$EndPhysicalNames");
          in >> line;
        }

      // but the next thing should,
      // in any case, be the list of
      // nodes:
      AssertThrow(line == "@f$Nodes", ExcInvalidGMSHInput(line));
    }

  // now read the nodes list
  in >> n_vertices;
  std::vector<Point<spacedim>> vertices(n_vertices);
  // set up mapping between numbering
  // in msh-file (nod) and in the
  // vertices vector
  std::map<int, int> vertex_indices;

  for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
    {
      int    vertex_number;
      double x[3];

      // read vertex
      in >> vertex_number >> x[0] >> x[1] >> x[2];

      for (unsigned int d = 0; d < spacedim; ++d)
        vertices[vertex](d) = x[d];
      // store mapping
      vertex_indices[vertex_number] = vertex;
    }

  // Assert we reached the end of the block
  in >> line;
  static const std::string end_nodes_marker[] = {"@f$ENDNOD", "@f$EndNodes"};
  AssertThrow(line == end_nodes_marker[gmsh_file_format - 1],
              ExcInvalidGMSHInput(line));

  // Now read in next bit
  in >> line;
  static const std::string begin_elements_marker[] = {"@f$ELM", "@f$Elements"};
  AssertThrow(line == begin_elements_marker[gmsh_file_format - 1],
              ExcInvalidGMSHInput(line));

  in >> n_cells;

  // set up array of cells and subcells (faces). In 1d, there is currently no
  // standard way in deal.II to pass boundary indicators attached to individual
  // vertices, so do this by hand via the boundary_ids_1d array
  std::vector<CellData<dim>>                 cells;
  SubCellData                                subcelldata;
  std::map<unsigned int, types::boundary_id> boundary_ids_1d;

  for (unsigned int cell = 0; cell < n_cells; ++cell)
    {
      // note that since in the input
      // file we found the number of
      // cells at the top, there
      // should still be input here,
      // so check this:
      AssertThrow(in, ExcIO());

      unsigned int cell_type;
      unsigned int material_id;
      unsigned int nod_num;

      /*
        For file format version 1, the format of each cell is as follows:
          elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list

        However, for version 2, the format reads like this:
          elm-number elm-type number-of-tags < tag > ... node-number-list

        In the following, we will ignore the element number (we simply enumerate
        them in the order in which we read them, and we will take reg-phys
        (version 1) or the first tag (version 2, if any tag is given at all) as
        material id.
      */

      unsigned int elm_number;
      in >> elm_number // ELM-NUMBER
        >> cell_type;  // ELM-TYPE

      switch (gmsh_file_format)
        {
          case 1:
            {
              in >> material_id // REG-PHYS
                >> dummy        // reg_elm
                >> nod_num;
              break;
            }

          case 2:
            {
              // read the tags; ignore all but the first one which we will
              // interpret as the material_id (for cells) or boundary_id
              // (for faces)
              unsigned int n_tags;
              in >> n_tags;
              if (n_tags > 0)
                in >> material_id;
              else
                material_id = 0;

              for (unsigned int i = 1; i < n_tags; ++i)
                in >> dummy;

              nod_num = GeometryInfo<dim>::vertices_per_cell;

              break;
            }

          default:
            AssertThrow(false, ExcNotImplemented());
        }


      /*       `ELM-TYPE'
               defines the geometrical type of the N-th element:
               `1'
               Line (2 nodes, 1 edge).

               `3'
               Quadrangle (4 nodes, 4 edges).

               `5'
               Hexahedron (8 nodes, 12 edges, 6 faces).

               `15'
               Point (1 node).
      */

      if (((cell_type == 1) && (dim == 1)) ||
          ((cell_type == 3) && (dim == 2)) || ((cell_type == 5) && (dim == 3)))
        // found a cell
        {
          AssertThrow(nod_num == GeometryInfo<dim>::vertices_per_cell,
                      ExcMessage("Number of nodes does not coincide with the "
                                 "number required for this object"));

          // allocate and read indices
          cells.emplace_back();
          for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell;
               ++i)
            in >> cells.back().vertices[i];

          // to make sure that the cast won't fail
          Assert(material_id <= std::numeric_limits<types::material_id>::max(),
                 ExcIndexRange(material_id,
                               0,
                               std::numeric_limits<types::material_id>::max()));
          // we use only material_ids in the range from 0 to
          // numbers::invalid_material_id-1
          Assert(material_id < numbers::invalid_material_id,
                 ExcIndexRange(material_id, 0, numbers::invalid_material_id));

          cells.back().material_id =
            static_cast<types::material_id>(material_id);

          // transform from ucd to
          // consecutive numbering
          for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell;
               ++i)
            {
              AssertThrow(vertex_indices.find(cells.back().vertices[i]) !=
                            vertex_indices.end(),
                          ExcInvalidVertexIndexGmsh(cell,
                                                    elm_number,
                                                    cells.back().vertices[i]));

              // vertex with this index exists
              cells.back().vertices[i] =
                vertex_indices[cells.back().vertices[i]];
            }
        }
      else if ((cell_type == 1) && ((dim == 2) || (dim == 3)))
        // boundary info
        {
          subcelldata.boundary_lines.emplace_back();
          in >> subcelldata.boundary_lines.back().vertices[0] >>
            subcelldata.boundary_lines.back().vertices[1];

          // to make sure that the cast won't fail
          Assert(material_id <= std::numeric_limits<types::boundary_id>::max(),
                 ExcIndexRange(material_id,
                               0,
                               std::numeric_limits<types::boundary_id>::max()));
          // we use only boundary_ids in the range from 0 to
          // numbers::internal_face_boundary_id-1
          Assert(material_id < numbers::internal_face_boundary_id,
                 ExcIndexRange(material_id,
                               0,
                               numbers::internal_face_boundary_id));

          subcelldata.boundary_lines.back().boundary_id =
            static_cast<types::boundary_id>(material_id);

          // transform from ucd to
          // consecutive numbering
          for (unsigned int i = 0; i < 2; ++i)
            if (vertex_indices.find(
                  subcelldata.boundary_lines.back().vertices[i]) !=
                vertex_indices.end())
              // vertex with this index exists
              subcelldata.boundary_lines.back().vertices[i] =
                vertex_indices[subcelldata.boundary_lines.back().vertices[i]];
            else
              {
                // no such vertex index
                AssertThrow(false,
                            ExcInvalidVertexIndex(
                              cell,
                              subcelldata.boundary_lines.back().vertices[i]));
                subcelldata.boundary_lines.back().vertices[i] =
                  numbers::invalid_unsigned_int;
              };
        }
      else if ((cell_type == 3) && (dim == 3))
        // boundary info
        {
          subcelldata.boundary_quads.emplace_back();
          in >> subcelldata.boundary_quads.back().vertices[0] >>
            subcelldata.boundary_quads.back().vertices[1] >>
            subcelldata.boundary_quads.back().vertices[2] >>
            subcelldata.boundary_quads.back().vertices[3];

          // to make sure that the cast won't fail
          Assert(material_id <= std::numeric_limits<types::boundary_id>::max(),
                 ExcIndexRange(material_id,
                               0,
                               std::numeric_limits<types::boundary_id>::max()));
          // we use only boundary_ids in the range from 0 to
          // numbers::internal_face_boundary_id-1
          Assert(material_id < numbers::internal_face_boundary_id,
                 ExcIndexRange(material_id,
                               0,
                               numbers::internal_face_boundary_id));

          subcelldata.boundary_quads.back().boundary_id =
            static_cast<types::boundary_id>(material_id);

          // transform from gmsh to
          // consecutive numbering
          for (unsigned int i = 0; i < 4; ++i)
            if (vertex_indices.find(
                  subcelldata.boundary_quads.back().vertices[i]) !=
                vertex_indices.end())
              // vertex with this index exists
              subcelldata.boundary_quads.back().vertices[i] =
                vertex_indices[subcelldata.boundary_quads.back().vertices[i]];
            else
              {
                // no such vertex index
                Assert(false,
                       ExcInvalidVertexIndex(
                         cell, subcelldata.boundary_quads.back().vertices[i]));
                subcelldata.boundary_quads.back().vertices[i] =
                  numbers::invalid_unsigned_int;
              }
        }
      else if (cell_type == 15)
        {
          // read the indices of nodes given
          unsigned int node_index = 0;
          switch (gmsh_file_format)
            {
              case 1:
                {
                  for (unsigned int i = 0; i < nod_num; ++i)
                    in >> node_index;
                  break;
                }
              case 2:
                {
                  in >> node_index;
                  break;
                }
              default:
                Assert(false, ExcInternalError());
            }

          // we only care about boundary indicators assigned to individual
          // vertices in 1d (because otherwise the vertices are not faces)
          if (dim == 1)
            boundary_ids_1d[vertex_indices[node_index]] = material_id;
        }
      else
        // cannot read this, so throw
        // an exception. treat
        // triangles and tetrahedra
        // specially since this
        // deserves a more explicit
        // error message
        {
          AssertThrow(cell_type != 2,
                      ExcMessage("Found triangles while reading a file "
                                 "in gmsh format. deal.II does not "
                                 "support triangles"));
          AssertThrow(cell_type != 11,
                      ExcMessage("Found tetrahedra while reading a file "
                                 "in gmsh format. deal.II does not "
                                 "support tetrahedra"));

          AssertThrow(false, ExcGmshUnsupportedGeometry(cell_type));
        }
    }

  // Assert we reached the end of the block
  in >> line;
  static const std::string end_elements_marker[] = {"@f$ENDELM", "@f$EndElements"};
  AssertThrow(line == end_elements_marker[gmsh_file_format - 1],
              ExcInvalidGMSHInput(line));

  // check that no forbidden arrays are used
  Assert(subcelldata.check_consistency(dim), ExcInternalError());

  AssertThrow(in, ExcIO());

  // check that we actually read some
  // cells.
  AssertThrow(cells.size() > 0, ExcGmshNoCellInformation());

  // do some clean-up on
  // vertices...
  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
  // ... and cells
  if (dim == spacedim)
    GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(vertices,
                                                                     cells);
  GridReordering<dim, spacedim>::reorder_cells(cells);
  tria->create_triangulation_compatibility(vertices, cells, subcelldata);

  // in 1d, we also have to attach boundary ids to vertices, which does not
  // currently work through the call above
  if (dim == 1)
    assign_1d_boundary_ids(boundary_ids_1d, *tria);
}


template <>
void
GridIn<1>::read_netcdf(const std::string &)
{
  AssertThrow(false, ExcImpossibleInDim(1));
}

template <>
void
GridIn<1, 2>::read_netcdf(const std::string &)
{
  AssertThrow(false, ExcImpossibleInDim(1));
}


template <>
void
GridIn<1, 3>::read_netcdf(const std::string &)
{
  AssertThrow(false, ExcImpossibleInDim(1));
}


template <>
void
GridIn<2, 3>::read_netcdf(const std::string &)
{
  Assert(false, ExcNotImplemented());
}

template <>
void
GridIn<2>::read_netcdf(const std::string &filename)
{
#ifndef DEAL_II_WITH_NETCDF
  (void)filename;
  AssertThrow(false, ExcNeedsNetCDF());
#else
  const unsigned int dim      = 2;
  const unsigned int spacedim = 2;
  Assert(tria != nullptr, ExcNoTriangulationSelected());
  // this function assumes the TAU
  // grid format.
  //
  // This format stores 2d grids as
  // 3d grids. In particular, a 2d
  // grid of n_cells quadrilaterals
  // in the y=0 plane is duplicated
  // to y=1 to build n_cells
  // hexaeders.  The surface
  // quadrilaterals of this 3d grid
  // are marked with boundary
  // marker. In the following we read
  // in all data required, find the
  // boundary marker associated with
  // the plane y=0, and extract the
  // corresponding 2d data to build a
  // Triangulation<2>.

  // In the following, we assume that
  // the 2d grid lies in the x-z
  // plane (y=0). I.e. we choose:
  // point[coord]=0, with coord=1
  const unsigned int coord = 1;
  // Also x-y-z (0-1-2) point
  // coordinates will be transformed
  // to x-y (x2d-y2d) coordinates.
  // With coord=1 as above, we have
  // x-z (0-2) -> (x2d-y2d)
  const unsigned int x2d = 0;
  const unsigned int y2d = 2;
  // For the case, the 2d grid lies
  // in x-y or y-z plane instead, the
  // following code must be extended
  // to find the right value for
  // coord, and setting x2d and y2d
  // accordingly.

  // First, open the file
  NcFile nc(filename.c_str());
  AssertThrow(nc.is_valid(), ExcIO());

  // then read n_cells
  NcDim *elements_dim = nc.get_dim("no_of_elements");
  AssertThrow(elements_dim->is_valid(), ExcIO());
  const unsigned int n_cells = elements_dim->size();

  // then we read
  //   int marker(no_of_markers)
  NcDim *marker_dim = nc.get_dim("no_of_markers");
  AssertThrow(marker_dim->is_valid(), ExcIO());
  const unsigned int n_markers = marker_dim->size();

  NcVar *marker_var = nc.get_var("marker");
  AssertThrow(marker_var->is_valid(), ExcIO());
  AssertThrow(marker_var->num_dims() == 1, ExcIO());
  AssertThrow(static_cast<unsigned int>(marker_var->get_dim(0)->size()) ==
                n_markers,
              ExcIO());

  std::vector<int> marker(n_markers);
  // use &* to convert
  // vector<int>::iterator to int *
  marker_var->get(&*marker.begin(), n_markers);

  // next we read
  // int boundarymarker_of_surfaces(
  //   no_of_surfaceelements)
  NcDim *bquads_dim = nc.get_dim("no_of_surfacequadrilaterals");
  AssertThrow(bquads_dim->is_valid(), ExcIO());
  const unsigned int n_bquads = bquads_dim->size();

  NcVar *bmarker_var = nc.get_var("boundarymarker_of_surfaces");
  AssertThrow(bmarker_var->is_valid(), ExcIO());
  AssertThrow(bmarker_var->num_dims() == 1, ExcIO());
  AssertThrow(static_cast<unsigned int>(bmarker_var->get_dim(0)->size()) ==
                n_bquads,
              ExcIO());

  std::vector<int> bmarker(n_bquads);
  bmarker_var->get(&*bmarker.begin(), n_bquads);

  // for each marker count the
  // number of boundary quads
  // which carry this marker
  std::map<int, unsigned int> n_bquads_per_bmarker;
  for (unsigned int i = 0; i < n_markers; ++i)
    {
      // the markers should all be
      // different
      AssertThrow(n_bquads_per_bmarker.find(marker[i]) ==
                    n_bquads_per_bmarker.end(),
                  ExcIO());

      n_bquads_per_bmarker[marker[i]] =
        count(bmarker.begin(), bmarker.end(), marker[i]);
    }
  // Note: the n_bquads_per_bmarker
  // map could be used to find the
  // right coord by finding the
  // marker0 such that
  // a/ n_bquads_per_bmarker[marker0]==n_cells
  // b/ point[coord]==0,
  // Condition a/ would hold for at
  // least two markers, marker0 and
  // marker1, whereas b/ would hold
  // for marker0 only. For marker1 we
  // then had point[coord]=constant
  // with e.g. constant=1 or -1

  // next we read
  // int points_of_surfacequadrilaterals(
  //   no_of_surfacequadrilaterals,
  //   points_per_surfacequadrilateral)
  NcDim *quad_vertices_dim = nc.get_dim("points_per_surfacequadrilateral");
  AssertThrow(quad_vertices_dim->is_valid(), ExcIO());
  const unsigned int vertices_per_quad = quad_vertices_dim->size();
  AssertThrow(vertices_per_quad == GeometryInfo<dim>::vertices_per_cell,
              ExcIO());

  NcVar *vertex_indices_var = nc.get_var("points_of_surfacequadrilaterals");
  AssertThrow(vertex_indices_var->is_valid(), ExcIO());
  AssertThrow(vertex_indices_var->num_dims() == 2, ExcIO());
  AssertThrow(static_cast<unsigned int>(
                vertex_indices_var->get_dim(0)->size()) == n_bquads,
              ExcIO());
  AssertThrow(static_cast<unsigned int>(
                vertex_indices_var->get_dim(1)->size()) == vertices_per_quad,
              ExcIO());

  std::vector<int> vertex_indices(n_bquads * vertices_per_quad);
  vertex_indices_var->get(&*vertex_indices.begin(),
                          n_bquads,
                          vertices_per_quad);

  for (unsigned int i = 0; i < vertex_indices.size(); ++i)
    AssertThrow(vertex_indices[i] >= 0, ExcIO());

  // next we read
  //   double points_xc(no_of_points)
  //   double points_yc(no_of_points)
  //   double points_zc(no_of_points)
  NcDim *vertices_dim = nc.get_dim("no_of_points");
  AssertThrow(vertices_dim->is_valid(), ExcIO());
  const unsigned int n_vertices = vertices_dim->size();

  NcVar *points_xc = nc.get_var("points_xc");
  NcVar *points_yc = nc.get_var("points_yc");
  NcVar *points_zc = nc.get_var("points_zc");
  AssertThrow(points_xc->is_valid(), ExcIO());
  AssertThrow(points_yc->is_valid(), ExcIO());
  AssertThrow(points_zc->is_valid(), ExcIO());
  AssertThrow(points_xc->num_dims() == 1, ExcIO());
  AssertThrow(points_yc->num_dims() == 1, ExcIO());
  AssertThrow(points_zc->num_dims() == 1, ExcIO());
  AssertThrow(points_yc->get_dim(0)->size() == static_cast<int>(n_vertices),
              ExcIO());
  AssertThrow(points_zc->get_dim(0)->size() == static_cast<int>(n_vertices),
              ExcIO());
  AssertThrow(points_xc->get_dim(0)->size() == static_cast<int>(n_vertices),
              ExcIO());
  std::vector<std::vector<double>> point_values(
    3, std::vector<double>(n_vertices));
  points_xc->get(&*point_values[0].begin(), n_vertices);
  points_yc->get(&*point_values[1].begin(), n_vertices);
  points_zc->get(&*point_values[2].begin(), n_vertices);

  // and fill the vertices
  std::vector<Point<spacedim>> vertices(n_vertices);
  for (unsigned int i = 0; i < n_vertices; ++i)
    {
      vertices[i](0) = point_values[x2d][i];
      vertices[i](1) = point_values[y2d][i];
    }

  // For all boundary quads in the
  // point[coord]=0 plane add the
  // bmarker to zero_plane_markers
  std::map<int, bool> zero_plane_markers;
  for (unsigned int quad = 0; quad < n_bquads; ++quad)
    {
      bool zero_plane = true;
      for (unsigned int i = 0; i < vertices_per_quad; ++i)
        if (point_values[coord][vertex_indices[quad * vertices_per_quad + i]] !=
            0)
          {
            zero_plane = false;
            break;
          }

      if (zero_plane)
        zero_plane_markers[bmarker[quad]] = true;
    }
  unsigned int sum_of_zero_plane_cells = 0;
  for (std::map<int, bool>::const_iterator iter = zero_plane_markers.begin();
       iter != zero_plane_markers.end();
       ++iter)
    sum_of_zero_plane_cells += n_bquads_per_bmarker[iter->first];
  AssertThrow(sum_of_zero_plane_cells == n_cells, ExcIO());

  // fill cells with all quads
  // associated with
  // zero_plane_markers
  std::vector<CellData<dim>> cells(n_cells);
  for (unsigned int quad = 0, cell = 0; quad < n_bquads; ++quad)
    {
      bool zero_plane = false;
      for (std::map<int, bool>::const_iterator iter =
             zero_plane_markers.begin();
           iter != zero_plane_markers.end();
           ++iter)
        if (bmarker[quad] == iter->first)
          {
            zero_plane = true;
            break;
          }

      if (zero_plane)
        {
          for (unsigned int i = 0; i < vertices_per_quad; ++i)
            {
              Assert(
                point_values[coord]
                            [vertex_indices[quad * vertices_per_quad + i]] == 0,
                ExcNotImplemented());
              cells[cell].vertices[i] =
                vertex_indices[quad * vertices_per_quad + i];
            }
          ++cell;
        }
    }

  SubCellData subcelldata;
  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
  GridReordering<dim, spacedim>::reorder_cells(cells);
  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
#endif
}


template <>
void
GridIn<3>::read_netcdf(const std::string &filename)
{
#ifndef DEAL_II_WITH_NETCDF
  // do something with the function argument
  // to make sure it at least looks used,
  // even if it is not
  (void)filename;
  AssertThrow(false, ExcNeedsNetCDF());
#else
  const unsigned int dim      = 3;
  const unsigned int spacedim = 3;
  Assert(tria != nullptr, ExcNoTriangulationSelected());
  // this function assumes the TAU
  // grid format.

  // First, open the file
  NcFile nc(filename.c_str());
  AssertThrow(nc.is_valid(), ExcIO());

  // then read n_cells
  NcDim *elements_dim = nc.get_dim("no_of_elements");
  AssertThrow(elements_dim->is_valid(), ExcIO());
  const unsigned int n_cells = elements_dim->size();
  // and n_hexes
  NcDim *hexes_dim = nc.get_dim("no_of_hexaeders");
  AssertThrow(hexes_dim->is_valid(), ExcIO());
  const unsigned int n_hexes = hexes_dim->size();
  AssertThrow(n_hexes == n_cells,
              ExcMessage("deal.II can handle purely hexaedral grids, only."));

  // next we read
  // int points_of_hexaeders(
  //   no_of_hexaeders,
  //   points_per_hexaeder)
  NcDim *hex_vertices_dim = nc.get_dim("points_per_hexaeder");
  AssertThrow(hex_vertices_dim->is_valid(), ExcIO());
  const unsigned int vertices_per_hex = hex_vertices_dim->size();
  AssertThrow(vertices_per_hex == GeometryInfo<dim>::vertices_per_cell,
              ExcIO());

  NcVar *vertex_indices_var = nc.get_var("points_of_hexaeders");
  AssertThrow(vertex_indices_var->is_valid(), ExcIO());
  AssertThrow(vertex_indices_var->num_dims() == 2, ExcIO());
  AssertThrow(static_cast<unsigned int>(
                vertex_indices_var->get_dim(0)->size()) == n_cells,
              ExcIO());
  AssertThrow(static_cast<unsigned int>(
                vertex_indices_var->get_dim(1)->size()) == vertices_per_hex,
              ExcIO());

  std::vector<int> vertex_indices(n_cells * vertices_per_hex);
  // use &* to convert
  // vector<int>::iterator to int *
  vertex_indices_var->get(&*vertex_indices.begin(), n_cells, vertices_per_hex);

  for (unsigned int i = 0; i < vertex_indices.size(); ++i)
    AssertThrow(vertex_indices[i] >= 0, ExcIO());

  // next we read
  //   double points_xc(no_of_points)
  //   double points_yc(no_of_points)
  //   double points_zc(no_of_points)
  NcDim *vertices_dim = nc.get_dim("no_of_points");
  AssertThrow(vertices_dim->is_valid(), ExcIO());
  const unsigned int n_vertices = vertices_dim->size();

  NcVar *points_xc = nc.get_var("points_xc");
  NcVar *points_yc = nc.get_var("points_yc");
  NcVar *points_zc = nc.get_var("points_zc");
  AssertThrow(points_xc->is_valid(), ExcIO());
  AssertThrow(points_yc->is_valid(), ExcIO());
  AssertThrow(points_zc->is_valid(), ExcIO());
  AssertThrow(points_xc->num_dims() == 1, ExcIO());
  AssertThrow(points_yc->num_dims() == 1, ExcIO());
  AssertThrow(points_zc->num_dims() == 1, ExcIO());
  AssertThrow(points_yc->get_dim(0)->size() == static_cast<int>(n_vertices),
              ExcIO());
  AssertThrow(points_zc->get_dim(0)->size() == static_cast<int>(n_vertices),
              ExcIO());
  AssertThrow(points_xc->get_dim(0)->size() == static_cast<int>(n_vertices),
              ExcIO());
  std::vector<std::vector<double>> point_values(
    3, std::vector<double>(n_vertices));
  points_xc->get(&*point_values[0].begin(), n_vertices);
  points_yc->get(&*point_values[1].begin(), n_vertices);
  points_zc->get(&*point_values[2].begin(), n_vertices);

  // and fill the vertices
  std::vector<Point<spacedim>> vertices(n_vertices);
  for (unsigned int i = 0; i < n_vertices; ++i)
    {
      vertices[i](0) = point_values[0][i];
      vertices[i](1) = point_values[1][i];
      vertices[i](2) = point_values[2][i];
    }

  // and cells
  std::vector<CellData<dim>> cells(n_cells);
  for (unsigned int cell = 0; cell < n_cells; ++cell)
    for (unsigned int i = 0; i < vertices_per_hex; ++i)
      cells[cell].vertices[i] = vertex_indices[cell * vertices_per_hex + i];

  // for setting up the SubCellData
  // we read the vertex indices of
  // the boundary quadrilaterals and
  // their boundary markers

  // first we read
  // int points_of_surfacequadrilaterals(
  //   no_of_surfacequadrilaterals,
  //   points_per_surfacequadrilateral)
  NcDim *quad_vertices_dim = nc.get_dim("points_per_surfacequadrilateral");
  AssertThrow(quad_vertices_dim->is_valid(), ExcIO());
  const unsigned int vertices_per_quad = quad_vertices_dim->size();
  AssertThrow(vertices_per_quad == GeometryInfo<dim>::vertices_per_face,
              ExcIO());

  NcVar *bvertex_indices_var = nc.get_var("points_of_surfacequadrilaterals");
  AssertThrow(bvertex_indices_var->is_valid(), ExcIO());
  AssertThrow(bvertex_indices_var->num_dims() == 2, ExcIO());
  const unsigned int n_bquads = bvertex_indices_var->get_dim(0)->size();
  AssertThrow(static_cast<unsigned int>(
                bvertex_indices_var->get_dim(1)->size()) ==
                GeometryInfo<dim>::vertices_per_face,
              ExcIO());

  std::vector<int> bvertex_indices(n_bquads * vertices_per_quad);
  bvertex_indices_var->get(&*bvertex_indices.begin(),
                           n_bquads,
                           vertices_per_quad);

  // next we read
  // int boundarymarker_of_surfaces(
  //   no_of_surfaceelements)
  NcDim *bquads_dim = nc.get_dim("no_of_surfacequadrilaterals");
  AssertThrow(bquads_dim->is_valid(), ExcIO());
  AssertThrow(static_cast<unsigned int>(bquads_dim->size()) == n_bquads,
              ExcIO());

  NcVar *bmarker_var = nc.get_var("boundarymarker_of_surfaces");
  AssertThrow(bmarker_var->is_valid(), ExcIO());
  AssertThrow(bmarker_var->num_dims() == 1, ExcIO());
  AssertThrow(static_cast<unsigned int>(bmarker_var->get_dim(0)->size()) ==
                n_bquads,
              ExcIO());

  std::vector<int> bmarker(n_bquads);
  bmarker_var->get(&*bmarker.begin(), n_bquads);
  // we only handle boundary
  // indicators that fit into an
  // types::boundary_id. Also, we don't
  // take numbers::internal_face_boundary_id
  // as it denotes an internal face
  for (unsigned int i = 0; i < bmarker.size(); ++i)
    Assert(0 <= bmarker[i] && static_cast<types::boundary_id>(bmarker[i]) !=
                                numbers::internal_face_boundary_id,
           ExcIO());

  // finally we setup the boundary
  // information
  SubCellData subcelldata;
  subcelldata.boundary_quads.resize(n_bquads);
  for (unsigned int i = 0; i < n_bquads; ++i)
    {
      for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
        subcelldata.boundary_quads[i].vertices[v] =
          bvertex_indices[i * GeometryInfo<dim>::vertices_per_face + v];
      subcelldata.boundary_quads[i].boundary_id =
        static_cast<types::boundary_id>(bmarker[i]);
    }

  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
  GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(vertices,
                                                                   cells);
  GridReordering<dim, spacedim>::reorder_cells(cells);
  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
#endif
}


template <int dim, int spacedim>
void
GridIn<dim, spacedim>::parse_tecplot_header(
  std::string &              header,
  std::vector<unsigned int> &tecplot2deal,
  unsigned int &             n_vars,
  unsigned int &             n_vertices,
  unsigned int &             n_cells,
  std::vector<unsigned int> &IJK,
  bool &                     structured,
  bool &                     blocked)
{
  Assert(tecplot2deal.size() == dim, ExcInternalError());
  Assert(IJK.size() == dim, ExcInternalError());
  // initialize the output variables
  n_vars     = 0;
  n_vertices = 0;
  n_cells    = 0;
  switch (dim)
    {
      case 3:
        IJK[2] = 0;
        DEAL_II_FALLTHROUGH;
      case 2:
        IJK[1] = 0;
        DEAL_II_FALLTHROUGH;
      case 1:
        IJK[0] = 0;
    }
  structured = true;
  blocked    = false;

  // convert the string to upper case
  std::transform(header.begin(), header.end(), header.begin(), ::toupper);

  // replace all tabs, commas, newlines by
  // whitespaces
  std::replace(header.begin(), header.end(), '\t', ' ');
  std::replace(header.begin(), header.end(), ',', ' ');
  std::replace(header.begin(), header.end(), '\n', ' ');

  // now remove whitespace in front of and
  // after '='
  std::string::size_type pos = header.find('=');

  while (pos != static_cast<std::string::size_type>(std::string::npos))
    if (header[pos + 1] == ' ')
      header.erase(pos + 1, 1);
    else if (header[pos - 1] == ' ')
      {
        header.erase(pos - 1, 1);
        --pos;
      }
    else
      pos = header.find('=', ++pos);

  // split the string into individual entries
  std::vector<std::string> entries =
    Utilities::break_text_into_lines(header, 1, ' ');

  // now go through the list and try to extract
  for (unsigned int i = 0; i < entries.size(); ++i)
    {
      if (Utilities::match_at_string_start(entries[i], "VARIABLES=\""))
        {
          ++n_vars;
          // we assume, that the first variable
          // is x or no coordinate at all (not y or z)
          if (Utilities::match_at_string_start(entries[i], "VARIABLES=\"X\""))
            {
              tecplot2deal[0] = 0;
            }
          ++i;
          while (entries[i][0] == '"')
            {
              if (entries[i] == "\"X\"")
                tecplot2deal[0] = n_vars;
              else if (entries[i] == "\"Y\"")
                {
                  // we assume, that y contains
                  // zero data in 1d, so do
                  // nothing
                  if (dim > 1)
                    tecplot2deal[1] = n_vars;
                }
              else if (entries[i] == "\"Z\"")
                {
                  // we assume, that z contains
                  // zero data in 1d and 2d, so
                  // do nothing
                  if (dim > 2)
                    tecplot2deal[2] = n_vars;
                }
              ++n_vars;
              ++i;
            }
          // set i back, so that the next
          // string is treated correctly
          --i;

          AssertThrow(
            n_vars >= dim,
            ExcMessage(
              "Tecplot file must contain at least one variable for each dimension"));
          for (unsigned int d = 1; d < dim; ++d)
            AssertThrow(
              tecplot2deal[d] > 0,
              ExcMessage(
                "Tecplot file must contain at least one variable for each dimension."));
        }
      else if (Utilities::match_at_string_start(entries[i], "ZONETYPE=ORDERED"))
        structured = true;
      else if (Utilities::match_at_string_start(entries[i],
                                                "ZONETYPE=FELINESEG") &&
               dim == 1)
        structured = false;
      else if (Utilities::match_at_string_start(entries[i],
                                                "ZONETYPE=FEQUADRILATERAL") &&
               dim == 2)
        structured = false;
      else if (Utilities::match_at_string_start(entries[i],
                                                "ZONETYPE=FEBRICK") &&
               dim == 3)
        structured = false;
      else if (Utilities::match_at_string_start(entries[i], "ZONETYPE="))
        // unsupported ZONETYPE
        {
          AssertThrow(false,
                      ExcMessage(
                        "The tecplot file contains an unsupported ZONETYPE."));
        }
      else if (Utilities::match_at_string_start(entries[i],
                                                "DATAPACKING=POINT"))
        blocked = false;
      else if (Utilities::match_at_string_start(entries[i],
                                                "DATAPACKING=BLOCK"))
        blocked = true;
      else if (Utilities::match_at_string_start(entries[i], "F=POINT"))
        {
          structured = true;
          blocked    = false;
        }
      else if (Utilities::match_at_string_start(entries[i], "F=BLOCK"))
        {
          structured = true;
          blocked    = true;
        }
      else if (Utilities::match_at_string_start(entries[i], "F=FEPOINT"))
        {
          structured = false;
          blocked    = false;
        }
      else if (Utilities::match_at_string_start(entries[i], "F=FEBLOCK"))
        {
          structured = false;
          blocked    = true;
        }
      else if (Utilities::match_at_string_start(entries[i],
                                                "ET=QUADRILATERAL") &&
               dim == 2)
        structured = false;
      else if (Utilities::match_at_string_start(entries[i], "ET=BRICK") &&
               dim == 3)
        structured = false;
      else if (Utilities::match_at_string_start(entries[i], "ET="))
        // unsupported ElementType
        {
          AssertThrow(
            false,
            ExcMessage(
              "The tecplot file contains an unsupported ElementType."));
        }
      else if (Utilities::match_at_string_start(entries[i], "I="))
        IJK[0] = Utilities::get_integer_at_position(entries[i], 2).first;
      else if (Utilities::match_at_string_start(entries[i], "J="))
        {
          IJK[1] = Utilities::get_integer_at_position(entries[i], 2).first;
          AssertThrow(
            dim > 1 || IJK[1] == 1,
            ExcMessage(
              "Parameter 'J=' found in tecplot, although this is only possible for dimensions greater than 1."));
        }
      else if (Utilities::match_at_string_start(entries[i], "K="))
        {
          IJK[2] = Utilities::get_integer_at_position(entries[i], 2).first;
          AssertThrow(
            dim > 2 || IJK[2] == 1,
            ExcMessage(
              "Parameter 'K=' found in tecplot, although this is only possible for dimensions greater than 2."));
        }
      else if (Utilities::match_at_string_start(entries[i], "N="))
        n_vertices = Utilities::get_integer_at_position(entries[i], 2).first;
      else if (Utilities::match_at_string_start(entries[i], "E="))
        n_cells = Utilities::get_integer_at_position(entries[i], 2).first;
    }

  // now we have read all the fields we are
  // interested in. do some checks and
  // calculate the variables
  if (structured)
    {
      n_vertices = 1;
      n_cells    = 1;
      for (unsigned int d = 0; d < dim; ++d)
        {
          AssertThrow(
            IJK[d] > 0,
            ExcMessage(
              "Tecplot file does not contain a complete and consistent set of parameters"));
          n_vertices *= IJK[d];
          n_cells *= (IJK[d] - 1);
        }
    }
  else
    {
      AssertThrow(
        n_vertices > 0,
        ExcMessage(
          "Tecplot file does not contain a complete and consistent set of parameters"));
      if (n_cells == 0)
        // this means an error, although
        // tecplot itself accepts entries like
        // 'J=20' instead of 'E=20'. therefore,
        // take the max of IJK
        n_cells = *std::max_element(IJK.begin(), IJK.end());
      AssertThrow(
        n_cells > 0,
        ExcMessage(
          "Tecplot file does not contain a complete and consistent set of parameters"));
    }
}



template <>
void
GridIn<2>::read_tecplot(std::istream &in)
{
  const unsigned int dim      = 2;
  const unsigned int spacedim = 2;
  Assert(tria != nullptr, ExcNoTriangulationSelected());
  AssertThrow(in, ExcIO());

  // skip comments at start of file
  skip_comment_lines(in, '#');

  // some strings for parsing the header
  std::string line, header;

  // first, concatenate all header lines
  // create a searchstring with almost all
  // letters. exclude e and E from the letters
  // to search, as they might appear in
  // exponential notation
  std::string letters = "abcdfghijklmnopqrstuvwxyzABCDFGHIJKLMNOPQRSTUVWXYZ";

  getline(in, line);
  while (line.find_first_of(letters) != std::string::npos)
    {
      header += " " + line;
      getline(in, line);
    }

  // now create some variables holding
  // important information on the mesh, get
  // this information from the header string
  std::vector<unsigned int> tecplot2deal(dim);
  std::vector<unsigned int> IJK(dim);
  unsigned int              n_vars, n_vertices, n_cells;
  bool                      structured, blocked;

  parse_tecplot_header(header,
                       tecplot2deal,
                       n_vars,
                       n_vertices,
                       n_cells,
                       IJK,
                       structured,
                       blocked);

  // reserve space for vertices. note, that in
  // tecplot vertices are ordered beginning
  // with 1, whereas in deal all indices start
  // with 0. in order not to use -1 for all the
  // connectivity information, a 0th vertex
  // (unused) is inserted at the origin.
  std::vector<Point<spacedim>> vertices(n_vertices + 1);
  vertices[0] = Point<spacedim>();
  // reserve space for cells
  std::vector<CellData<dim>> cells(n_cells);
  SubCellData                subcelldata;

  if (blocked)
    {
      // blocked data format. first we get all
      // the values of the first variable for
      // all points, after that we get all
      // values for the second variable and so
      // on.

      // dummy variable to read in all the info
      // we do not want to use
      double dummy;
      // which is the first index to read in
      // the loop (see below)
      unsigned int next_index = 0;

      // note, that we have already read the
      // first line containing the first variable
      if (tecplot2deal[0] == 0)
        {
          // we need the information in this
          // line, so extract it
          std::vector<std::string> first_var =
            Utilities::break_text_into_lines(line, 1);
          char *endptr;
          for (unsigned int i = 1; i < first_var.size() + 1; ++i)
            vertices[i](0) = std::strtod(first_var[i - 1].c_str(), &endptr);

          // if there are many points, the data
          // for this var might continue in the
          // next line(s)
          for (unsigned int j = first_var.size() + 1; j < n_vertices + 1; ++j)
            in >> vertices[j](next_index);
          // now we got all values of the first
          // variable, so increase the counter
          next_index = 1;
        }

      // main loop over all variables
      for (unsigned int i = 1; i < n_vars; ++i)
        {
          // if we read all the important
          // variables and do not want to
          // read further, because we are
          // using a structured grid, we can
          // stop here (and skip, for
          // example, a whole lot of solution
          // variables)
          if (next_index == dim && structured)
            break;

          if ((next_index < dim) && (i == tecplot2deal[next_index]))
            {
              // we need this line, read it in
              for (unsigned int j = 1; j < n_vertices + 1; ++j)
                in >> vertices[j](next_index);
              ++next_index;
            }
          else
            {
              // we do not need this line, read
              // it in and discard it
              for (unsigned int j = 1; j < n_vertices + 1; ++j)
                in >> dummy;
            }
        }
      Assert(next_index == dim, ExcInternalError());
    }
  else
    {
      // the data is not blocked, so we get all
      // the variables for one point, then the
      // next and so on. create a vector to
      // hold these components
      std::vector<double> vars(n_vars);

      // now fill the first vertex. note, that we
      // have already read the first line
      // containing the first vertex
      std::vector<std::string> first_vertex =
        Utilities::break_text_into_lines(line, 1);
      char *endptr;
      for (unsigned int d = 0; d < dim; ++d)
        vertices[1](d) =
          std::strtod(first_vertex[tecplot2deal[d]].c_str(), &endptr);

      // read the remaining vertices from the
      // list
      for (unsigned int v = 2; v < n_vertices + 1; ++v)
        {
          for (unsigned int i = 0; i < n_vars; ++i)
            in >> vars[i];
          // fill the vertex
          // coordinates. respect the position
          // of coordinates in the list of
          // variables
          for (unsigned int i = 0; i < dim; ++i)
            vertices[v](i) = vars[tecplot2deal[i]];
        }
    }

  if (structured)
    {
      // this is the part of the code that only
      // works in 2d
      unsigned int I = IJK[0], J = IJK[1];

      unsigned int cell = 0;
      // set up array of cells
      for (unsigned int j = 0; j < J - 1; ++j)
        for (unsigned int i = 1; i < I; ++i)
          {
            cells[cell].vertices[0] = i + j * I;
            cells[cell].vertices[1] = i + 1 + j * I;
            cells[cell].vertices[2] = i + 1 + (j + 1) * I;
            cells[cell].vertices[3] = i + (j + 1) * I;
            ++cell;
          }
      Assert(cell == n_cells, ExcInternalError());
      std::vector<unsigned int> boundary_vertices(2 * I + 2 * J - 4);
      unsigned int              k = 0;
      for (unsigned int i = 1; i < I + 1; ++i)
        {
          boundary_vertices[k] = i;
          ++k;
          boundary_vertices[k] = i + (J - 1) * I;
          ++k;
        }
      for (unsigned int j = 1; j < J - 1; ++j)
        {
          boundary_vertices[k] = 1 + j * I;
          ++k;
          boundary_vertices[k] = I + j * I;
          ++k;
        }
      Assert(k == boundary_vertices.size(), ExcInternalError());
      // delete the duplicated vertices at the
      // boundary, which occur, e.g. in c-type
      // or o-type grids around a body
      // (airfoil). this automatically deletes
      // unused vertices as well.
      GridTools::delete_duplicated_vertices(vertices,
                                            cells,
                                            subcelldata,
                                            boundary_vertices);
    }
  else
    {
      // set up array of cells, unstructured
      // mode, so the connectivity is
      // explicitly given
      for (unsigned int i = 0; i < n_cells; ++i)
        {
          // note that since in the input file
          // we found the number of cells at
          // the top, there should still be
          // input here, so check this:
          AssertThrow(in, ExcIO());

          // get the connectivity from the
          // input file. the vertices are
          // ordered like in the ucd format
          for (unsigned int j = 0; j < GeometryInfo<dim>::vertices_per_cell;
               ++j)
            in >> cells[i].vertices[j];
        }
      // do some clean-up on vertices
      GridTools::delete_unused_vertices(vertices, cells, subcelldata);
    }

  // check that no forbidden arrays are
  // used. as we do not read in any
  // subcelldata, nothing should happen here.
  Assert(subcelldata.check_consistency(dim), ExcInternalError());
  AssertThrow(in, ExcIO());

  // do some cleanup on cells
  GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(vertices,
                                                                   cells);
  GridReordering<dim, spacedim>::reorder_cells(cells);
  tria->create_triangulation_compatibility(vertices, cells, subcelldata);
}



template <int dim, int spacedim>
void
GridIn<dim, spacedim>::read_tecplot(std::istream &)
{
  Assert(false, ExcNotImplemented());
}



template <int dim, int spacedim>
void
GridIn<dim, spacedim>::read_assimp(const std::string &filename,
                                   const unsigned int mesh_index,
                                   const bool         remove_duplicates,
                                   const double       tol,
                                   const bool         ignore_unsupported_types)
{
#ifdef DEAL_II_WITH_ASSIMP
  // Only good for surface grids.
  AssertThrow(dim < 3, ExcImpossibleInDim(dim));

  // Create an instance of the Importer class
  Assimp::Importer importer;

  // And have it read the given file with some  postprocessing
  const aiScene *scene =
    importer.ReadFile(filename.c_str(),
                      aiProcess_RemoveComponent |
                        aiProcess_JoinIdenticalVertices |
                        aiProcess_ImproveCacheLocality | aiProcess_SortByPType |
                        aiProcess_OptimizeGraph | aiProcess_OptimizeMeshes);

  // If the import failed, report it
  AssertThrow(scene, ExcMessage(importer.GetErrorString()))

    AssertThrow(scene->mNumMeshes,
                ExcMessage("Input file contains no meshes."));

  AssertThrow((mesh_index == numbers::invalid_unsigned_int) ||
                (mesh_index < scene->mNumMeshes),
              ExcMessage("Too few meshes in the file."));

  unsigned int start_mesh =
    (mesh_index == numbers::invalid_unsigned_int ? 0 : mesh_index);
  unsigned int end_mesh =
    (mesh_index == numbers::invalid_unsigned_int ? scene->mNumMeshes :
                                                   mesh_index + 1);

  // Deal.II objects are created empty, and then filled with imported file.
  std::vector<Point<spacedim>> vertices;
  std::vector<CellData<dim>>   cells;
  SubCellData                  subcelldata;

  // A series of counters to merge cells.
  unsigned int v_offset = 0;
  unsigned int c_offset = 0;

  // The index of the mesh will be used as a material index.
  for (unsigned int m = start_mesh; m < end_mesh; ++m)
    {
      const aiMesh *mesh = scene->mMeshes[m];

      // Check that we know what to do with this mesh, otherwise just
      // ignore it
      if ((dim == 2) && mesh->mPrimitiveTypes != aiPrimitiveType_POLYGON)
        {
          AssertThrow(ignore_unsupported_types,
                      ExcMessage("Incompatible mesh " + std::to_string(m) +
                                 "/" + std::to_string(scene->mNumMeshes)));
          continue;
        }
      else if ((dim == 1) && mesh->mPrimitiveTypes != aiPrimitiveType_LINE)
        {
          AssertThrow(ignore_unsupported_types,
                      ExcMessage("Incompatible mesh " + std::to_string(m) +
                                 "/" + std::to_string(scene->mNumMeshes)));
          continue;
        }
      // Vertices
      const unsigned int n_vertices = mesh->mNumVertices;
      const aiVector3D * mVertices  = mesh->mVertices;

      // Faces
      const unsigned int n_faces = mesh->mNumFaces;
      const aiFace *     mFaces  = mesh->mFaces;

      vertices.resize(v_offset + n_vertices);
      cells.resize(c_offset + n_faces);

      for (unsigned int i = 0; i < n_vertices; ++i)
        for (unsigned int d = 0; d < spacedim; ++d)
          vertices[i + v_offset][d] = mVertices[i][d];

      unsigned int valid_cell = c_offset;
      for (unsigned int i = 0; i < n_faces; ++i)
        {
          if (mFaces[i].mNumIndices == GeometryInfo<dim>::vertices_per_cell)
            {
              for (unsigned int f = 0; f < GeometryInfo<dim>::vertices_per_cell;
                   ++f)
                {
                  cells[valid_cell].vertices[f] =
                    mFaces[i].mIndices[f] + v_offset;
                }
              cells[valid_cell].material_id = (types::material_id)m;
              ++valid_cell;
            }
          else
            {
              AssertThrow(ignore_unsupported_types,
                          ExcMessage("Face " + std::to_string(i) + " of mesh " +
                                     std::to_string(m) + " has " +
                                     std::to_string(mFaces[i].mNumIndices) +
                                     " vertices. We expected only " +
                                     std::to_string(
                                       GeometryInfo<dim>::vertices_per_cell)));
            }
        }
      cells.resize(valid_cell);

      // The vertices are added all at once. Cells are checked for
      // validity, so only valid_cells are now present in the deal.II
      // list of cells.
      v_offset += n_vertices;
      c_offset = valid_cell;
    }

  // No cells were read
  if (cells.size() == 0)
    return;

  if (remove_duplicates)
    {
      // The function delete_duplicated_vertices() needs to be called more
      // than once if a vertex is duplicated more than once. So we keep
      // calling it until the number of vertices does not change any more.
      unsigned int n_verts = 0;
      while (n_verts != vertices.size())
        {
          n_verts = vertices.size();
          std::vector<unsigned int> considered_vertices;
          GridTools::delete_duplicated_vertices(
            vertices, cells, subcelldata, considered_vertices, tol);
        }
    }

  GridTools::delete_unused_vertices(vertices, cells, subcelldata);
  if (dim == spacedim)
    GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(vertices,
                                                                     cells);

  GridReordering<dim, spacedim>::reorder_cells(cells);
  if (dim == 2)
    tria->create_triangulation_compatibility(vertices, cells, subcelldata);
  else
    tria->create_triangulation(vertices, cells, subcelldata);
#else
  (void)filename;
  (void)mesh_index;
  (void)remove_duplicates;
  (void)tol;
  (void)ignore_unsupported_types;
  Assert(false, ExcNeedsAssimp());
#endif
}


template <int dim, int spacedim>
void
GridIn<dim, spacedim>::skip_empty_lines(std::istream &in)
{
  std::string line;
  while (in)
    {
      // get line
      getline(in, line);

      // check if this is a line that
      // consists only of spaces, and
      // if not put the whole thing
      // back and return
      if (std::find_if(line.begin(),
                       line.end(),
                       std::bind(std::not_equal_to<char>(),
                                 std::placeholders::_1,
                                 ' ')) != line.end())
        {
          in.putback('\n');
          for (int i = line.length() - 1; i >= 0; --i)
            in.putback(line[i]);
          return;
        }

      // else: go on with next line
    }
}



template <int dim, int spacedim>
void
GridIn<dim, spacedim>::skip_comment_lines(std::istream &in,
                                          const char    comment_start)
{
  char c;
  // loop over the following comment
  // lines
  while (in.get(c) && c == comment_start)
    // loop over the characters after
    // the comment starter
    while (in.get() != '\n')
      ;


  // put back first character of
  // first non-comment line
  if (in)
    in.putback(c);

  // at last: skip additional empty lines, if present
  skip_empty_lines(in);
}



template <int dim, int spacedim>
void
GridIn<dim, spacedim>::debug_output_grid(
  const std::vector<CellData<dim>> & /*cells*/,
  const std::vector<Point<spacedim>> & /*vertices*/,
  std::ostream & /*out*/)
{
  Assert(false, ExcNotImplemented());
}



template <>
void
GridIn<2>::debug_output_grid(const std::vector<CellData<2>> &cells,
                             const std::vector<Point<2>> &   vertices,
                             std::ostream &                  out)
{
  double min_x = vertices[cells[0].vertices[0]](0),
         max_x = vertices[cells[0].vertices[0]](0),
         min_y = vertices[cells[0].vertices[0]](1),
         max_y = vertices[cells[0].vertices[0]](1);

  for (unsigned int i = 0; i < cells.size(); ++i)
    {
      for (unsigned int v = 0; v < 4; ++v)
        {
          const Point<2> &p = vertices[cells[i].vertices[v]];

          if (p(0) < min_x)
            min_x = p(0);
          if (p(0) > max_x)
            max_x = p(0);
          if (p(1) < min_y)
            min_y = p(1);
          if (p(1) > max_y)
            max_y = p(1);
        };

      out << "# cell " << i << std::endl;
      Point<2> center;
      for (unsigned int f = 0; f < 4; ++f)
        center += vertices[cells[i].vertices[f]];
      center /= 4;

      out << "set label \"" << i << "\" at " << center(0) << ',' << center(1)
          << " center" << std::endl;

      // first two line right direction
      for (unsigned int f = 0; f < 2; ++f)
        out << "set arrow from " << vertices[cells[i].vertices[f]](0) << ','
            << vertices[cells[i].vertices[f]](1) << " to "
            << vertices[cells[i].vertices[(f + 1) % 4]](0) << ','
            << vertices[cells[i].vertices[(f + 1) % 4]](1) << std::endl;
      // other two lines reverse direction
      for (unsigned int f = 2; f < 4; ++f)
        out << "set arrow from " << vertices[cells[i].vertices[(f + 1) % 4]](0)
            << ',' << vertices[cells[i].vertices[(f + 1) % 4]](1) << " to "
            << vertices[cells[i].vertices[f]](0) << ','
            << vertices[cells[i].vertices[f]](1) << std::endl;
      out << std::endl;
    };


  out << std::endl
      << "set nokey" << std::endl
      << "pl [" << min_x << ':' << max_x << "][" << min_y << ':' << max_y
      << "] " << min_y << std::endl
      << "pause -1" << std::endl;
}



template <>
void
GridIn<3>::debug_output_grid(const std::vector<CellData<3>> &cells,
                             const std::vector<Point<3>> &   vertices,
                             std::ostream &                  out)
{
  for (unsigned int cell = 0; cell < cells.size(); ++cell)
    {
      // line 0
      out << vertices[cells[cell].vertices[0]] << std::endl
          << vertices[cells[cell].vertices[1]] << std::endl
          << std::endl
          << std::endl;
      // line 1
      out << vertices[cells[cell].vertices[1]] << std::endl
          << vertices[cells[cell].vertices[2]] << std::endl
          << std::endl
          << std::endl;
      // line 2
      out << vertices[cells[cell].vertices[3]] << std::endl
          << vertices[cells[cell].vertices[2]] << std::endl
          << std::endl
          << std::endl;
      // line 3
      out << vertices[cells[cell].vertices[0]] << std::endl
          << vertices[cells[cell].vertices[3]] << std::endl
          << std::endl
          << std::endl;
      // line 4
      out << vertices[cells[cell].vertices[4]] << std::endl
          << vertices[cells[cell].vertices[5]] << std::endl
          << std::endl
          << std::endl;
      // line 5
      out << vertices[cells[cell].vertices[5]] << std::endl
          << vertices[cells[cell].vertices[6]] << std::endl
          << std::endl
          << std::endl;
      // line 6
      out << vertices[cells[cell].vertices[7]] << std::endl
          << vertices[cells[cell].vertices[6]] << std::endl
          << std::endl
          << std::endl;
      // line 7
      out << vertices[cells[cell].vertices[4]] << std::endl
          << vertices[cells[cell].vertices[7]] << std::endl
          << std::endl
          << std::endl;
      // line 8
      out << vertices[cells[cell].vertices[0]] << std::endl
          << vertices[cells[cell].vertices[4]] << std::endl
          << std::endl
          << std::endl;
      // line 9
      out << vertices[cells[cell].vertices[1]] << std::endl
          << vertices[cells[cell].vertices[5]] << std::endl
          << std::endl
          << std::endl;
      // line 10
      out << vertices[cells[cell].vertices[2]] << std::endl
          << vertices[cells[cell].vertices[6]] << std::endl
          << std::endl
          << std::endl;
      // line 11
      out << vertices[cells[cell].vertices[3]] << std::endl
          << vertices[cells[cell].vertices[7]] << std::endl
          << std::endl
          << std::endl;
    };
}



template <int dim, int spacedim>
void
GridIn<dim, spacedim>::read(const std::string &filename, Format format)
{
  // Search file class for meshes
  PathSearch  search("MESH");
  std::string name;
  // Open the file and remember its name
  if (format == Default)
    name = search.find(filename);
  else
    name = search.find(filename, default_suffix(format));

  std::ifstream in(name.c_str());

  if (format == Default)
    {
      const std::string::size_type slashpos = name.find_last_of('/');
      const std::string::size_type dotpos   = name.find_last_of('.');
      if (dotpos < name.length() &&
          (dotpos > slashpos || slashpos == std::string::npos))
        {
          std::string ext = name.substr(dotpos + 1);
          format          = parse_format(ext);
        }
    }
  if (format == netcdf)
    read_netcdf(filename);
  else
    read(in, format);
}


template <int dim, int spacedim>
void
GridIn<dim, spacedim>::read(std::istream &in, Format format)
{
  if (format == Default)
    format = default_format;

  switch (format)
    {
      case dbmesh:
        read_dbmesh(in);
        return;

      case msh:
        read_msh(in);
        return;

      case vtk:
        read_vtk(in);
        return;

      case unv:
        read_unv(in);
        return;

      case ucd:
        read_ucd(in);
        return;

      case abaqus:
        read_abaqus(in);
        return;

      case xda:
        read_xda(in);
        return;

      case netcdf:
        Assert(false,
               ExcMessage("There is no read_netcdf(istream &) function. "
                          "Use the read_netcdf(string &filename) "
                          "functions, instead."));
        return;

      case tecplot:
        read_tecplot(in);
        return;

      case assimp:
        Assert(false,
               ExcMessage("There is no read_assimp(istream &) function. "
                          "Use the read_assimp(string &filename, ...) "
                          "functions, instead."));
        return;

      case Default:
        break;
    }
  Assert(false, ExcInternalError());
}



template <int dim, int spacedim>
std::string
GridIn<dim, spacedim>::default_suffix(const Format format)
{
  switch (format)
    {
      case dbmesh:
        return ".dbmesh";
      case msh:
        return ".msh";
      case vtk:
        return ".vtk";
      case unv:
        return ".unv";
      case ucd:
        return ".inp";
      case abaqus:
        return ".inp"; // Typical suffix for Abaqus mesh files conflicts with
                       // UCD.
      case xda:
        return ".xda";
      case netcdf:
        return ".nc";
      case tecplot:
        return ".dat";
      default:
        Assert(false, ExcNotImplemented());
        return ".unknown_format";
    }
}



template <int dim, int spacedim>
typename GridIn<dim, spacedim>::Format
GridIn<dim, spacedim>::parse_format(const std::string &format_name)
{
  if (format_name == "dbmesh")
    return dbmesh;

  if (format_name == "msh")
    return msh;

  if (format_name == "unv")
    return unv;

  if (format_name == "vtk")
    return vtk;

  // This is also the typical extension of Abaqus input files.
  if (format_name == "inp")
    return ucd;

  if (format_name == "ucd")
    return ucd;

  if (format_name == "xda")
    return xda;

  if (format_name == "netcdf")
    return netcdf;

  if (format_name == "nc")
    return netcdf;

  if (format_name == "tecplot")
    return tecplot;

  if (format_name == "dat")
    return tecplot;

  if (format_name == "plt")
    // Actually, this is the extension for the
    // tecplot binary format, which we do not
    // support right now. However, some people
    // tend to create tecplot ascii files with
    // the extension 'plt' instead of
    // 'dat'. Thus, include this extension
    // here. If it actually is a binary file,
    // the read_tecplot() function will fail
    // and throw an exception, anyway.
    return tecplot;

  AssertThrow(false, ExcInvalidState());
  // return something weird
  return Format(Default);
}



template <int dim, int spacedim>
std::string
GridIn<dim, spacedim>::get_format_names()
{
  return "dbmesh|msh|unv|vtk|ucd|abaqus|xda|netcdf|tecplot|assimp";
}

namespace
{
  template <int dim>
  Abaqus_to_UCD<dim>::Abaqus_to_UCD()
    : tolerance(5e-16) // Used to offset Cubit tolerance error when outputting
                       // value close to zero
  {
    AssertThrow(dim == 2 || dim == 3, ExcNotImplemented());
  }

  // Convert from a string to some other data type
  // Reference: http://www.codeguru.com/forum/showthread.php?t=231054
  template <class T>
  bool
  from_string(T &t, const std::string &s, std::ios_base &(*f)(std::ios_base &))
  {
    std::istringstream iss(s);
    return !(iss >> f >> t).fail();
  }

  // Extract an integer from a string
  int
  extract_int(const std::string &s)
  {
    std::string tmp;
    for (unsigned int i = 0; i < s.size(); ++i)
      {
        if (isdigit(s[i]))
          {
            tmp += s[i];
          }
      }

    int number = 0;
    from_string(number, tmp, std::dec);
    return number;
  }

  template <int dim>
  void
  Abaqus_to_UCD<dim>::read_in_abaqus(std::istream &input_stream)
  {
    // References:
    // http://www.egr.msu.edu/software/abaqus/Documentation/docs/v6.7/books/usb/default.htm?startat=pt01ch02.html
    // http://www.cprogramming.com/tutorial/string.html

    AssertThrow(input_stream, ExcIO());
    std::string line;
    std::getline(input_stream, line);

    while (!input_stream.fail() && !input_stream.eof())
      {
        std::transform(line.begin(), line.end(), line.begin(), ::toupper);

        if (line.compare("*HEADING") == 0 || line.compare(0, 2, "**") == 0 ||
            line.compare(0, 5, "*PART") == 0)
          {
            // Skip header and comments
            while (!input_stream.fail() && !input_stream.eof())
              {
                std::getline(input_stream, line);
                if (line[0] == '*')
                  goto cont; // My eyes, they burn!
              }
          }
        else if (line.compare(0, 5, "*NODE") == 0)
          {
            // Extract list of vertices
            // Header line might be:
            // *NODE, NSET=ALLNODES
            // *NODE

            // Contains lines in the form:
            // Index, x, y, z
            while (!input_stream.fail() && !input_stream.eof())
              {
                std::getline(input_stream, line);
                if (line[0] == '*')
                  goto cont;

                std::vector<double> node(dim + 1);

                std::istringstream iss(line);
                char               comma;
                for (unsigned int i = 0; i < dim + 1; ++i)
                  iss >> node[i] >> comma;

                node_list.push_back(node);
              }
          }
        else if (line.compare(0, 8, "*ELEMENT") == 0)
          {
            // Element construction.
            // There are different header formats, the details
            // of which we're not particularly interested in except
            // whether they represent quads or hexahedrals.
            // *ELEMENT, TYPE=S4R, ELSET=EB<material id>
            // *ELEMENT, TYPE=C3D8R, ELSET=EB<material id>
            // *ELEMENT, TYPE=C3D8
            // Elements itself (n=4 or n=8):
            // Index, i[0], ..., i[n]

            int material = 0;
            // Scan for material id
            {
              const std::string before_material = "ELSET=EB";
              const std::size_t idx             = line.find(before_material);
              if (idx != std::string::npos)
                {
                  from_string(material,
                              line.substr(idx + before_material.size()),
                              std::dec);
                }
            }

            // Read ELEMENT definition
            std::getline(input_stream, line);
            while (!input_stream.fail() && !input_stream.eof())
              {
                if (line[0] == '*')
                  goto cont;

                std::istringstream iss(line);
                char               comma;

                // We will store the material id in the zeroth entry of the
                // vector and the rest of the elements represent the global
                // node numbers
                const unsigned int n_data_per_cell =
                  1 + GeometryInfo<dim>::vertices_per_cell;
                std::vector<double> cell(n_data_per_cell);
                for (unsigned int i = 0; i < n_data_per_cell; ++i)
                  iss >> cell[i] >> comma;

                // Overwrite cell index from file by material
                cell[0] = static_cast<double>(material);
                cell_list.push_back(cell);

                std::getline(input_stream, line);
              }
          }
        else if (line.compare(0, 8, "*SURFACE") == 0)
          {
            // Extract the definitions of boundary surfaces
            // Old format from Cubit:
            // *SURFACE, NAME=SS<boundary indicator>
            //    <element index>,     S<face number>
            // Abaqus default format:
            // *SURFACE, TYPE=ELEMENT, NAME=SURF-<indicator>

            // Get name of the surface and extract id from it;
            // this will be the boundary indicator
            const std::string name_key = "NAME=";
            const std::size_t name_idx_start =
              line.find(name_key) + name_key.size();
            std::size_t name_idx_end = line.find(',', name_idx_start);
            if (name_idx_end == std::string::npos)
              {
                name_idx_end = line.size();
              }
            const int b_indicator = extract_int(
              line.substr(name_idx_start, name_idx_end - name_idx_start));

            // Read SURFACE definition
            // Note that the orientation of the faces is embedded within the
            // definition of each "set" of faces that comprise the surface
            // These are either marked by an "S" or "E" in 3d or 2d
            // respectively.
            std::getline(input_stream, line);
            while (!input_stream.fail() && !input_stream.eof())
              {
                if (line[0] == '*')
                  goto cont;

                // Change all characters to upper case
                std::transform(line.begin(),
                               line.end(),
                               line.begin(),
                               ::toupper);

                // Surface can be created from ELSET, or directly from cells
                // If elsets_list contains a key with specific name - refers to
                // that ELSET, otherwise refers to cell
                std::istringstream iss(line);
                int                el_idx;
                int                face_number;
                char               temp;

                // Get relevant faces, taking into account the element
                // orientation
                std::vector<double> quad_node_list;
                const std::string   elset_name = line.substr(0, line.find(','));
                if (elsets_list.count(elset_name) != 0)
                  {
                    // Surface refers to ELSET
                    std::string stmp;
                    iss >> stmp >> temp >> face_number;

                    const std::vector<int> cells = elsets_list[elset_name];
                    for (unsigned int i = 0; i < cells.size(); ++i)
                      {
                        el_idx = cells[i];
                        quad_node_list =
                          get_global_node_numbers(el_idx, face_number);
                        quad_node_list.insert(quad_node_list.begin(),
                                              b_indicator);

                        face_list.push_back(quad_node_list);
                      }
                  }
                else
                  {
                    // Surface refers directly to elements
                    char comma;
                    iss >> el_idx >> comma >> temp >> face_number;
                    quad_node_list =
                      get_global_node_numbers(el_idx, face_number);
                    quad_node_list.insert(quad_node_list.begin(), b_indicator);

                    face_list.push_back(quad_node_list);
                  }

                std::getline(input_stream, line);
              }
          }
        else if (line.compare(0, 6, "*ELSET") == 0)
          {
            // Get ELSET name.
            // Materials are attached to elsets with specific name
            std::string elset_name;
            {
              const std::string elset_key = "*ELSET, ELSET=";
              const std::size_t idx       = line.find(elset_key);
              if (idx != std::string::npos)
                {
                  const std::string comma       = ",";
                  const std::size_t first_comma = line.find(comma);
                  const std::size_t second_comma =
                    line.find(comma, first_comma + 1);
                  const std::size_t elset_name_start =
                    line.find(elset_key) + elset_key.size();
                  elset_name = line.substr(elset_name_start,
                                           second_comma - elset_name_start);
                }
            }

            // There are two possibilities of storing cells numbers in ELSET:
            // 1. If the header contains the 'GENERATE' keyword, then the next
            // line describes range of cells as:
            //    cell_id_start, cell_id_end, cell_step
            // 2. If the header does not contain the 'GENERATE' keyword, then
            // the next lines contain cells numbers
            std::vector<int>  elements;
            const std::size_t generate_idx = line.find("GENERATE");
            if (generate_idx != std::string::npos)
              {
                // Option (1)
                std::getline(input_stream, line);
                std::istringstream iss(line);
                char               comma;
                int                elid_start;
                int                elid_end;
                int                elis_step =
                  1; // Default value set in case stride not provided
                // Some files don't have the stride size
                // Compare mesh test cases ./grids/abaqus/3d/other_simple.inp to
                // ./grids/abaqus/2d/2d_test_abaqus.inp
                iss >> elid_start >> comma >> elid_end;
                // https://stackoverflow.com/questions/8046357/how-do-i-check-if-a-stringstream-variable-is-empty-null
                if (iss.rdbuf()->in_avail() != 0)
                  iss >> comma >> elis_step;
                for (int i = elid_start; i <= elid_end; i += elis_step)
                  {
                    elements.push_back(i);
                  }
                elsets_list[elset_name] = elements;

                std::getline(input_stream, line);
              }
            else
              {
                // Option (2)
                std::getline(input_stream, line);
                while (!input_stream.fail() && !input_stream.eof())
                  {
                    if (line[0] == '*')
                      break;

                    std::istringstream iss(line);
                    char               comma;
                    int                elid;
                    while (!iss.eof())
                      {
                        iss >> elid >> comma;
                        elements.push_back(elid);
                      }

                    std::getline(input_stream, line);
                  }

                elsets_list[elset_name] = elements;
              }

            goto cont;
          }
        else if (line.compare(0, 5, "*NSET") == 0)
          {
            // Skip nodesets; we have no use for them
            while (!input_stream.fail() && !input_stream.eof())
              {
                std::getline(input_stream, line);
                if (line[0] == '*')
                  goto cont;
              }
          }
        else if (line.compare(0, 14, "*SOLID SECTION") == 0)
          {
            // The ELSET name, which describes a section for particular material
            const std::string elset_key = "ELSET=";
            const std::size_t elset_start =
              line.find("ELSET=") + elset_key.size();
            const std::size_t elset_end = line.find(',', elset_start + 1);
            const std::string elset_name =
              line.substr(elset_start, elset_end - elset_start);

            // Solid material definition.
            // We assume that material id is taken from material name,
            // eg. "Material-1" -> ID=1
            const std::string material_key = "MATERIAL=";
            const std::size_t last_equal =
              line.find("MATERIAL=") + material_key.size();
            const std::size_t material_id_start = line.find('-', last_equal);
            int               material_id       = 0;
            from_string(material_id,
                        line.substr(material_id_start + 1),
                        std::dec);

            // Assign material id to cells
            const std::vector<int> &elset_cells = elsets_list[elset_name];
            for (unsigned int i = 0; i < elset_cells.size(); ++i)
              {
                const int cell_id     = elset_cells[i] - 1;
                cell_list[cell_id][0] = material_id;
              }
          }
        // Note: All other lines / entries are ignored

        std::getline(input_stream, line);

      cont:
        (void)0;
      }
  }

  template <int dim>
  std::vector<double>
  Abaqus_to_UCD<dim>::get_global_node_numbers(const int face_cell_no,
                                              const int face_cell_face_no) const
  {
    std::vector<double> quad_node_list(GeometryInfo<dim>::vertices_per_face);

    // These orderings were reverse engineered by hand and may
    // conceivably be erroneous.
    // TODO: Currently one test (2d unstructured mesh) in the test
    // suite fails, presumably because of an ordering issue.
    if (dim == 2)
      {
        if (face_cell_face_no == 1)
          {
            quad_node_list[0] = cell_list[face_cell_no - 1][1];
            quad_node_list[1] = cell_list[face_cell_no - 1][2];
          }
        else if (face_cell_face_no == 2)
          {
            quad_node_list[0] = cell_list[face_cell_no - 1][2];
            quad_node_list[1] = cell_list[face_cell_no - 1][3];
          }
        else if (face_cell_face_no == 3)
          {
            quad_node_list[0] = cell_list[face_cell_no - 1][3];
            quad_node_list[1] = cell_list[face_cell_no - 1][4];
          }
        else if (face_cell_face_no == 4)
          {
            quad_node_list[0] = cell_list[face_cell_no - 1][4];
            quad_node_list[1] = cell_list[face_cell_no - 1][1];
          }
        else
          {
            AssertThrow(face_cell_face_no <= 4,
                        ExcMessage("Invalid face number in 2d"));
          }
      }
    else if (dim == 3)
      {
        if (face_cell_face_no == 1)
          {
            quad_node_list[0] = cell_list[face_cell_no - 1][1];
            quad_node_list[1] = cell_list[face_cell_no - 1][4];
            quad_node_list[2] = cell_list[face_cell_no - 1][3];
            quad_node_list[3] = cell_list[face_cell_no - 1][2];
          }
        else if (face_cell_face_no == 2)
          {
            quad_node_list[0] = cell_list[face_cell_no - 1][5];
            quad_node_list[1] = cell_list[face_cell_no - 1][8];
            quad_node_list[2] = cell_list[face_cell_no - 1][7];
            quad_node_list[3] = cell_list[face_cell_no - 1][6];
          }
        else if (face_cell_face_no == 3)
          {
            quad_node_list[0] = cell_list[face_cell_no - 1][1];
            quad_node_list[1] = cell_list[face_cell_no - 1][2];
            quad_node_list[2] = cell_list[face_cell_no - 1][6];
            quad_node_list[3] = cell_list[face_cell_no - 1][5];
          }
        else if (face_cell_face_no == 4)
          {
            quad_node_list[0] = cell_list[face_cell_no - 1][2];
            quad_node_list[1] = cell_list[face_cell_no - 1][3];
            quad_node_list[2] = cell_list[face_cell_no - 1][7];
            quad_node_list[3] = cell_list[face_cell_no - 1][6];
          }
        else if (face_cell_face_no == 5)
          {
            quad_node_list[0] = cell_list[face_cell_no - 1][3];
            quad_node_list[1] = cell_list[face_cell_no - 1][4];
            quad_node_list[2] = cell_list[face_cell_no - 1][8];
            quad_node_list[3] = cell_list[face_cell_no - 1][7];
          }
        else if (face_cell_face_no == 6)
          {
            quad_node_list[0] = cell_list[face_cell_no - 1][1];
            quad_node_list[1] = cell_list[face_cell_no - 1][5];
            quad_node_list[2] = cell_list[face_cell_no - 1][8];
            quad_node_list[3] = cell_list[face_cell_no - 1][4];
          }
        else
          {
            AssertThrow(face_cell_no <= 6,
                        ExcMessage("Invalid face number in 3d"));
          }
      }
    else
      {
        AssertThrow(dim == 2 || dim == 3, ExcNotImplemented());
      }

    return quad_node_list;
  }

  template <int dim>
  void
  Abaqus_to_UCD<dim>::write_out_avs_ucd(std::ostream &output) const
  {
    // References:
    // http://www.dealii.org/developer/doxygen/deal.II/structGeometryInfo.html
    // http://people.scs.fsu.edu/~burkardt/data/ucd/ucd.html

    AssertThrow(output, ExcIO());

    // save old formatting options
    const boost::io::ios_base_all_saver formatting_saver(output);

    // Write out title - Note: No other commented text can be inserted below the
    // title in a UCD file
    output << "# Abaqus to UCD mesh conversion" << std::endl;
    output << "# Mesh type: AVS UCD" << std::endl;

    // ========================================================
    // ASCII UCD File Format
    // The input file cannot contain blank lines or lines with leading blanks.
    // Comments, if present, must precede all data in the file.
    // Comments within the data will cause read errors.
    // The general order of the data is as follows:
    // 1. Numbers defining the overall structure, including the number of nodes,
    //    the number of cells, and the length of the vector of data associated
    //    with the nodes, cells, and the model.
    //     e.g. 1:
    //        <num_nodes> <num_cells> <num_ndata> <num_cdata> <num_mdata>
    //     e.g. 2:
    //        n_elements = n_hex_cells + n_bc_quads + n_quad_cells + n_bc_edges
    //        outfile.write(str(n_nodes) + " " + str(n_elements) + " 0 0 0\n")
    // 2. For each node, its node id and the coordinates of that node in space.
    //    Node-ids must be integers, but any number including non sequential
    //    numbers can be used. Mid-edge nodes are treated like any other node.
    // 3. For each cell: its cell-id, material, cell type (hexahedral, pyramid,
    //    etc.), and the list of node-ids that correspond to each of the cell's
    //    vertices. The below table specifies the different cell types and the
    //    keyword used to represent them in the file.

    // Write out header
    output << node_list.size() << "\t" << (cell_list.size() + face_list.size())
           << "\t0\t0\t0" << std::endl;

    output.width(16);
    output.precision(8);

    // Write out node numbers
    // Loop over all nodes
    for (unsigned int ii = 0; ii < node_list.size(); ++ii)
      {
        // Node number
        output << node_list[ii][0] << "\t";

        // Node coordinates
        output.setf(std::ios::scientific, std::ios::floatfield);
        for (unsigned int jj = 1; jj < dim + 1; ++jj)
          {
            // invoke tolerance -> set points close to zero equal to zero
            if (std::abs(node_list[ii][jj]) > tolerance)
              output << static_cast<double>(node_list[ii][jj]) << "\t";
            else
              output << 0.0 << "\t";
          }
        if (dim == 2)
          output << 0.0 << "\t";

        output << std::endl;
        output.unsetf(std::ios::floatfield);
      }

    // Write out cell node numbers
    for (unsigned int ii = 0; ii < cell_list.size(); ++ii)
      {
        output << ii + 1 << "\t" << cell_list[ii][0] << "\t"
               << (dim == 2 ? "quad" : "hex") << "\t";
        for (unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_cell + 1;
             ++jj)
          output << cell_list[ii][jj] << "\t";

        output << std::endl;
      }

    // Write out quad node numbers
    for (unsigned int ii = 0; ii < face_list.size(); ++ii)
      {
        output << ii + 1 << "\t" << face_list[ii][0] << "\t"
               << (dim == 2 ? "line" : "quad") << "\t";
        for (unsigned int jj = 1; jj < GeometryInfo<dim>::vertices_per_face + 1;
             ++jj)
          output << face_list[ii][jj] << "\t";

        output << std::endl;
      }
  }
} // namespace


// explicit instantiations
#include "grid_in.inst"

DEAL_II_NAMESPACE_CLOSE