mint/mesh/UnstructuredMesh.hpp

// Copyright (c) 2017-2021, Lawrence Livermore National Security, LLC and
// other Axom Project Developers. See the top-level LICENSE file for details.
//
// SPDX-License-Identifier: (BSD-3-Clause)

#ifndef MINT_UNSTRUCTUREDMESH_HPP_
#define MINT_UNSTRUCTUREDMESH_HPP_

// Axom includes
#include "axom/core/Macros.hpp"
#include "axom/core/Types.hpp"

// Mint includes
#include "axom/mint/config.hpp"
#include "axom/mint/mesh/CellTypes.hpp"
#include "axom/mint/mesh/ConnectivityArray.hpp"
#include "axom/mint/mesh/Mesh.hpp"
#include "axom/mint/mesh/MeshCoordinates.hpp"
#include "axom/mint/mesh/MeshTypes.hpp"
#include "axom/mint/mesh/internal/MeshHelpers.hpp"

// Slic includes
#include "axom/slic/interface/slic.hpp"

// C/C++ includes
#include <cstring>  // for std::memcpy

namespace axom
{
namespace mint
{
enum Topology
{
  SINGLE_SHAPE,
  MIXED_SHAPE
};

template <Topology TOPO>
struct topology_traits
{ };

template <>
struct topology_traits<SINGLE_SHAPE>
{
  constexpr static ConnectivityType cell_to_nodes = NO_INDIRECTION;
  constexpr static ConnectivityType cell_to_faces = NO_INDIRECTION;
  constexpr static ConnectivityType face_to_nodes = NO_INDIRECTION;
};

template <>
struct topology_traits<MIXED_SHAPE>
{
  constexpr static ConnectivityType cell_to_nodes = TYPED_INDIRECTION;
  constexpr static ConnectivityType cell_to_faces = INDIRECTION;
  constexpr static ConnectivityType face_to_nodes = TYPED_INDIRECTION;
};

/*!
 * \class UnstructuredMesh
 *
 * \brief Provides the ability to store and operate on Unstructured meshes.
 *
 *  An <em> unstructured mesh </em> stores both node and topology information
 *  explicitly. This allows the flexibility of discretizing the solution
 *  domain using a variety of cell types not just quadrilateral (in 2D) or
 *  hexahedral (in 3D) cells. Due to this added flexibility, the use of
 *  <em> unstructured meshes </em> is more common when dealing with complex
 *  geometries. However, <em> unstructured meshes </em> require additional
 *  storage and generally incur some performance penalty to store, create and
 *  access mesh topology information respectively.
 *
 *  Mint classifies <em> unstructured meshes </em> in two basic types based on
 *  the underlying mesh topology:
 *
 *  * <b> Single Shape Topology </b>
 *
 *   In this case, the <em> unstructured mesh </em> consists of a single cell
 *   type, e.g., a quad or triangle mesh in 2D, or, a hex or tet mesh in 3D.
 *   In this case  the underlying implementation is optimized for the
 *   specified cell type (specified in the constructor).
 *
 *  * <b> Mixed Shape Topology </b>
 *
 *    When <em> mixed cell topology </em> is specified, the <em> unstructured
 *    mesh </em> can be composed of any of the supported cell types, e.g.,
 *    a mesh consisting of both quads and triangles. This mode incurs
 *    additional overhead for storage and access to mesh topology information,
 *    since it requires indirection.
 *
 *  The list of supported cell types for an <em> unstructured mesh </em> is
 *  available in CellTypes.hpp
 *
 * An UnstructuredMesh object may be constructed using (a) native storage,
 * (b) external storage, or, (c) with Sidre storage when Mint is compiled
 *  with Sidre support:
 *
 *  * <b> Native Storage </b> <br />
 *
 *    When using native storage, the Unstructured object owns all memory
 *    associated with the particle data. The storage can grow dynamically as
 *    needed by the application, i.e., adding more nodes/cells. Once the
 *    mesh object goes out-of-scope, all memory associated with it is
 *    returned to the system.
 *
 *  * <b> External Storage </b> <br />
 *
 *    An UnstructuredMesh may also be constructed from external, user-supplied
 *    buffers. In this case, all memory associated with the mesh data
 *    is owned by the caller. Consequently, the mesh cannot grow dynamically.
 *    All calls to resize(), append(), etc. will fail with an error.
 *
 *  * <b> Sidre </b> <br />
 *
 *    An UnstructuredMesh may also be constructed using Sidre as the back-end
 *    data-store manager. The mesh data is laid out in Sidre according to the
 *    <a href="http://llnl-conduit.readthedocs.io/en/latest/"> computational
 *    mesh blueprint </a> conventions. In this case, all operations are
 *    supported, including dynamically resizing the mesh and growing the
 *    associated storage as needed. However, Sidre owns all the memory. Once the
 *    mesh object goes out-of-scope, the data remains persistent in Sidre.
 *
 * \see mint::Mesh
 * \see mint::CellTypes
 */
template <Topology TOPO>
class UnstructuredMesh : public Mesh
{
  AXOM_STATIC_ASSERT(TOPO == SINGLE_SHAPE || TOPO == MIXED_SHAPE);

public:
  /*! \brief The types for face-cell and cell-face connectivity.
   *
   * Usually, each face will connect two cells.  The exceptions are
   * - edge faces, which will connect to one cell.  The other will be -1.
   * - some faces in a malformed mesh, which may join more than two cells.
   *   The concensus is that storing such a non-manifold mesh is not useful
   *   so we only store the first two incident cells and return an error.
   *
   * Face types are stored in the face-node connectivity.
   */
  using CellToNodeConnectivity =
    ConnectivityArray<topology_traits<TOPO>::cell_to_nodes>;
  using CellToFaceConnectivity =
    ConnectivityArray<topology_traits<TOPO>::cell_to_faces>;
  using FaceToCellConnectivity = ConnectivityArray<NO_INDIRECTION>;
  using FaceToNodeConnectivity =
    ConnectivityArray<topology_traits<TOPO>::face_to_nodes>;

  /*!
   * \brief Default constructor. Disabled.
   */
  UnstructuredMesh() = delete;

  /// \name Native Storage Constructors
  /// @{

  /*!
   * \brief Constructs an Unstructured single topology mesh.
   *
   * \param [in] ndims the number of dimensions.
   * \param [in] cell_type the cell type of the mesh.
   * \param [in] node_capacity the number of nodes to allocate space for.
   * \param [in] cell_capacity the number of cells to allocate space for.
   *
   * \note This constructor is only active when TOPO == SINGLE_SHAPE.
   *
   * \post getCellType() == cell_type
   * \post getNumberOfNodes() == 0
   * \post getNumberOfCells() == 0
   */
  UnstructuredMesh(int ndims,
                   CellType cell_type,
                   IndexType node_capacity = USE_DEFAULT,
                   IndexType cell_capacity = USE_DEFAULT)
    : Mesh(ndims, UNSTRUCTURED_MESH)
    , m_coordinates(new MeshCoordinates(ndims, 0, node_capacity))
    , m_cell_to_node(new CellToNodeConnectivity(cell_type, cell_capacity))
    , m_cell_to_face(initializeCellToFace(cell_type))
    , m_face_to_cell(new FaceToCellConnectivity(2, 0))
    , m_face_to_node(initializeFaceToNode(cell_type))
  {
    AXOM_STATIC_ASSERT_MSG(
      TOPO == SINGLE_SHAPE,
      "This constructor is only active for single topology meshes.");

    SLIC_ERROR_IF(
      cell_type == PRISM || cell_type == PYRAMID,
      "Single shape unstructured meshes do not support prisms or pyramids");

    initialize();
  }

  /*!
   * \brief Constructs an Unstructured mixed topology mesh.
   *
   * \param [in] ndims the number of dimensions.
   * \param [in] node_capacity the number of nodes to allocate space for.
   * \param [in] cell_capacity the number of cells to allocate space for.
   * \param [in] connectivity_capacity the number of vertices to allocate space
   *  for in the cell connectivity array.
   *
   * \note This constructor is only active when TOPO == MIXED_SHAPE.
   *
   * \post getCellType() == UNDEFINED_CELL
   * \post getNumberOfNodes() == 0
   * \post getNumberOfCells() == 0
   */
  UnstructuredMesh(int ndims,
                   IndexType node_capacity = USE_DEFAULT,
                   IndexType cell_capacity = USE_DEFAULT,
                   IndexType connectivity_capacity = USE_DEFAULT)
    : Mesh(ndims, UNSTRUCTURED_MESH)
    , m_coordinates(new MeshCoordinates(ndims, 0, node_capacity))
    , m_cell_to_node(
        new CellToNodeConnectivity(cell_capacity, connectivity_capacity))
    , m_cell_to_face(initializeCellToFace())
    , m_face_to_cell(new FaceToCellConnectivity(2, 0))
    , m_face_to_node(initializeFaceToNode())
  {
    AXOM_STATIC_ASSERT_MSG(
      TOPO == MIXED_SHAPE,
      "This constructor is only active for mixed topology meshes.");

    m_has_mixed_topology = true;
    initialize();
  }

  /// @}

  /// \name External Storage Constructors
  /// @{

  /*!
   * \brief Constructs an Unstructured single topology mesh using the provided
   *  external buffers.
   *
   * \param [in] cell_type the cell type of the mesh.
   * \param [in] n_cells the number of cells in the mesh.
   * \param [in] cell_capacity max number of cells the mesh is able to hold.
   * \param [in] connectivity the cell connectivity array.
   * \param [in] n_nodes the number of nodes in the mesh.
   * \param [in] node_capacity max number of nodes the mesh is able to hold
   * \param [in] x pointer to the x-coordinates
   * \param [in] y pointer to the y-coordinates (required only for 2D and 3D)
   * \param [in] z pointer to the z-coordinates (required only for 3D)
   *
   * \note The length of the connectivity array must be at least
   *  cell_capacity * getCellInfo( cell_type ).num_nodes.
   * \note The provided coordinate arrays have to be of length at least
   *  node_capacity.
   *
   * \note This constructor is only supported when TOPO == SINGLE_SHAPE.
   *
   * \post getCellType() == cell_type
   * \post getNumberOfCells() == n_cells
   * \post getCellCapacity == cell_capacity
   * \post getNumberOfNodes() == n_nodes
   * \post getNodeCapacity() == node_capacity
   * \post isExternal() == true
   */
  UnstructuredMesh(CellType cell_type,
                   IndexType n_cells,
                   IndexType cell_capacity,
                   IndexType* connectivity,
                   IndexType n_nodes,
                   IndexType node_capacity,
                   double* x,
                   double* y = nullptr,
                   double* z = nullptr)
    : Mesh(internal::dim(x, y, z), UNSTRUCTURED_MESH)
    , m_coordinates(new MeshCoordinates(n_nodes, node_capacity, x, y, z))
    , m_cell_to_node(
        new CellToNodeConnectivity(cell_type, n_cells, connectivity, cell_capacity))
    , m_cell_to_face(initializeCellToFace(cell_type))
    , m_face_to_cell(new FaceToCellConnectivity(2, 0))
    , m_face_to_node(initializeFaceToNode(cell_type))
  {
    AXOM_STATIC_ASSERT_MSG(
      TOPO == SINGLE_SHAPE,
      "This constructor is only active for single topology meshes.");

    SLIC_ERROR_IF(
      cell_type == PRISM || cell_type == PYRAMID,
      "Single shape unstructured meshes do not support prisms or pyramids");

    SLIC_ASSERT(x != nullptr);
    SLIC_ASSERT(m_ndims < 2 || y != nullptr);
    SLIC_ASSERT(m_ndims < 3 || z != nullptr);

    initialize();
  }

  /*!
   * \brief Constructs an Unstructured single topology mesh using the provided
   *  external buffers.
   *
   * \param [in] cell_type the cell type of the mesh
   * \param [in] n_cells the number of cells in the mesh
   * \param [in] connectivity the cell connectivity array.
   * \param [in] n_nodes the number of nodes in the mesh
   * \param [in] x pointer to the x-coordinates
   * \param [in] y pointer to the y-coordinates (required only for 2D or 3D)
   * \param [in] z pointer to the z-coordinates (required only for 3D)
   *
   * \note This constructor is only supported when TOPO == SINGLE_SHAPE.
   *
   * \post getCellType() == cell_type
   * \post getNumberOfCells() == n_cells
   * \post getCellCapacity == n_cells
   * \post getNumberOfNodes() == n_nodes
   * \post getNodeCapacity() == n_nodes
   * \post isExternal() == true
   */
  UnstructuredMesh(CellType cell_type,
                   IndexType n_cells,
                   IndexType* connectivity,
                   IndexType n_nodes,
                   double* x,
                   double* y = nullptr,
                   double* z = nullptr)
    : UnstructuredMesh(cell_type,
                       n_cells,
                       n_cells,
                       connectivity,
                       n_nodes,
                       n_nodes,
                       x,
                       y,
                       z)
  { }

  /*!
   * \brief Constructs an Unstructured mixed topology mesh using the provided
   *  external buffers.
   *
   * \param [in] n_cells the number of cells in the mesh.
   * \param [in] cell_capacity max number of cells the mesh is able to hold
   * \param [in] connectivity_capacity max capacity of the connectivity array
   * \param [in] connectivity the connectivity array.
   * \param [in] offsets array of cell offsets (of length cell_capacity + 1)
   * \param [in] types array of cell types (of length cell_capacity)
   * \param [in] n_nodes the number of nodes in the mesh.
   * \param [in] node_capacity max number of nodes the mesh can hold
   * \param [in] x pointer to the x-coordinates
   * \param [in] y pointer to the y-coordinates (required only for 2D and 3D)
   * \param [in] z pointer to the z-coordinates (required only 3D)
   *
   * \note The supplied connectivity array must have a length that is at least
   *  equal to the specified connectivity_capacity.
   * \note the provided coordinate arrays are to be of length at least
   *  node_capacity.
   *
   * \note This constructor is only supported when TOPO == MIXED_SHAPE.
   *
   * \post getCellType() == UNDEFINED_CELL
   * \post getNumberOfCells() == n_cells
   * \post getCellCapacity == cell_capacity
   * \post getNumberOfNodes() == n_nodes
   * \post getNodeCapacity() == node_capacity
   * \post isExternal() == true
   */
  UnstructuredMesh(IndexType n_cells,
                   IndexType cell_capacity,
                   IndexType connectivity_capacity,
                   IndexType* connectivity,
                   IndexType* offsets,
                   CellType* types,
                   IndexType n_nodes,
                   IndexType node_capacity,
                   double* x,
                   double* y = nullptr,
                   double* z = nullptr)
    : Mesh(internal::dim(x, y, z), UNSTRUCTURED_MESH)
    , m_coordinates(new MeshCoordinates(n_nodes, node_capacity, x, y, z))
    , m_cell_to_node(new CellToNodeConnectivity(n_cells,
                                                connectivity,
                                                offsets,
                                                types,
                                                cell_capacity,
                                                connectivity_capacity))
    , m_cell_to_face(initializeCellToFace())
    , m_face_to_cell(new FaceToCellConnectivity(2, 0))
    , m_face_to_node(initializeFaceToNode())
  {
    AXOM_STATIC_ASSERT_MSG(
      TOPO == MIXED_SHAPE,
      "This constructor is only active for mixed topology meshes.");

    SLIC_ASSERT(x != nullptr);
    SLIC_ASSERT(m_ndims < 2 || y != nullptr);
    SLIC_ASSERT(m_ndims < 3 || z != nullptr);

    m_has_mixed_topology = true;
    initialize();
  }

  /*!
   * \brief Constructs an Unstructured mixed topology mesh using the provided
   *  external mesh buffers.
   *
   * \param [in] n_cells the number of cells in the mesh.
   * \param [in] connectivity_size the size of the connectivity array
   * \param [in] connectivity the connectivity array
   * \param [in] offsets array of cell offsets (of length n_cells+1)
   * \param [in] types array of cell types (of length n_cells)
   * \param [in] n_nodes the number of nodes in the mesh
   * \param [in] x pointer to the x-coordinates
   * \param [in] y pointer to the y-coordinates (required only for 2D and 3D)
   * \param [in] z pointer to the z-coordinates (required only for 3D)
   *
   * \note The supplied connectivity array must have a length that is at least
   *  equal to the specified connectivity_capacity.
   * \note the provided coordinate arrays are to be of length at least
   *  node_capacity.
   *
   * \post getCellType() == UNDEFINED_CELL
   * \post getNumberOfCells() == n_cells
   * \post getCellCapacity == n_cells
   * \post getNumberOfNodes() == n_nodes
   * \post getNodeCapacity() == n_cells
   * \post isExternal() == true
   */
  UnstructuredMesh(IndexType n_cells,
                   IndexType connectivity_size,
                   IndexType* connectivity,
                   IndexType* offsets,
                   CellType* types,
                   IndexType n_nodes,
                   double* x,
                   double* y = nullptr,
                   double* z = nullptr)
    : UnstructuredMesh(n_cells,
                       n_cells,
                       connectivity_size,
                       connectivity,
                       offsets,
                       types,
                       n_nodes,
                       n_nodes,
                       x,
                       y,
                       z)
  { }

  /// @}

  /// \name Sidre Storage Constructors
  /// @{

#ifdef AXOM_MINT_USE_SIDRE

  /*!
   * \brief Creates an UnstructuredMesh instance from a given Sidre group that
   *  holds mesh data for an unstructured mesh according to the conventions
   *  described in the computational mesh blueprint.
   *
   * \param [in] group the sidre::Group to use.
   * \param [in] topo optional argument specifying the name of the topology
   *  associated with this Mesh instance.
   *
   * \note If a topology name is not provided, the implementation will construct
   *  a mesh based on the 1st topology group under the parent "topologies"
   *  group.
   *
   * \pre group != nullptr.
   * \pre blueprint::isValidRootGroup( group ) == true
   * \post isInSidre() == true
   */
  UnstructuredMesh(sidre::Group* group, const std::string& topo = "")
    : Mesh(group, topo)
    , m_coordinates(new MeshCoordinates(getCoordsetGroup()))
    , m_cell_to_node(new CellToNodeConnectivity(getTopologyGroup()))
    , m_cell_to_face(nullptr)
    , m_face_to_cell(new FaceToCellConnectivity(2, 0))
    , m_face_to_node(nullptr)
  {
    SLIC_ERROR_IF(
      m_type != UNSTRUCTURED_MESH,
      "Supplied sidre::Group does not correspond to a UnstructuredMesh.");

    if(TOPO == MIXED_SHAPE)
    {
      m_has_mixed_topology = true;
      m_cell_to_face = initializeCellToFace();
      m_face_to_node = initializeFaceToNode();
    }
    else
    {
      SLIC_ERROR_IF(
        getCellType() == PRISM || getCellType() == PYRAMID,
        "Single shape unstructured meshes do not support prisms or pyramids");
      m_cell_to_face = initializeCellToFace(getCellType());
      m_face_to_node = initializeFaceToNode(getCellType());
    }

    initialize();
  }

  /*!
   * \brief Creates an UnstructuredMesh instance on an empty Sidre group.
   *
   * \param [in] ndims the number of dimensions.
   * \param [in] cell_type the cell type of the mesh.
   * \param [in] group the sidre::Group to use.
   * \param [in] topo the name of the associated topology group.
   * \param [in] coordset the name of the associated coordset group.
   * \param [in] node_capacity the number of nodes to allocate space for.
   * \param [in] cell_capacity the number of cells to allocate space for.
   * \param [in] connectivity_capacity the number of vertices to allocate space
   *  for in the cell connectivity array.
   *
   * \note If a topology and coordset name are not provided a default name is
   *  used by the implementation.
   * \note The first two constructors are only active when
   *  TOPO == SINGLE_SHAPE and the last two are active only when
   *  TOPO == MIXED_SHAPE.
   *
   * \pre group != nullptr.
   * \pre group->getNumGroups() == 0
   * \pre group->getNumViews() == 0
   * \post blueprint::isValidRootGroup( group )
   * \post getNumberOfNodes() == 0
   * \post getNumberOfCells() == 0
   * \post isInSidre() == true
   */
  /// @{

  UnstructuredMesh(int ndims,
                   CellType cell_type,
                   sidre::Group* group,
                   const std::string& topo,
                   const std::string& coordset,
                   IndexType node_capacity = USE_DEFAULT,
                   IndexType cell_capacity = USE_DEFAULT)
    : Mesh(ndims, UNSTRUCTURED_MESH, group, topo, coordset)
    , m_coordinates(
        new MeshCoordinates(getCoordsetGroup(), ndims, 0, node_capacity))
    , m_cell_to_node(new CellToNodeConnectivity(cell_type,
                                                getTopologyGroup(),
                                                getCoordsetName(),
                                                cell_capacity))
    , m_cell_to_face(initializeCellToFace(cell_type))
    , m_face_to_cell(new FaceToCellConnectivity(2, 0))
    , m_face_to_node(initializeFaceToNode(cell_type))
  {
    AXOM_STATIC_ASSERT_MSG(
      TOPO == SINGLE_SHAPE,
      "This constructor is only active for single topology meshes.");

    SLIC_ERROR_IF(
      cell_type == PRISM || cell_type == PYRAMID,
      "Single shape unstructured meshes do not support prisms or pyramids");

    initialize();
  }

  UnstructuredMesh(int ndims,
                   CellType cell_type,
                   sidre::Group* group,
                   IndexType node_capacity = USE_DEFAULT,
                   IndexType cell_capacity = USE_DEFAULT)
    : UnstructuredMesh(ndims, cell_type, group, "", "", node_capacity, cell_capacity)
  { }

  UnstructuredMesh(int ndims,
                   sidre::Group* group,
                   const std::string& topo,
                   const std::string& coordset,
                   IndexType node_capacity = USE_DEFAULT,
                   IndexType cell_capacity = USE_DEFAULT,
                   IndexType connectivity_capacity = USE_DEFAULT)
    : Mesh(ndims, UNSTRUCTURED_MESH, group, topo, coordset)
    , m_coordinates(
        new MeshCoordinates(getCoordsetGroup(), ndims, 0, node_capacity))
    , m_cell_to_node(new CellToNodeConnectivity(getTopologyGroup(),
                                                getCoordsetName(),
                                                cell_capacity,
                                                connectivity_capacity))
    , m_cell_to_face(initializeCellToFace())
    , m_face_to_cell(new FaceToCellConnectivity(2, 0))
    , m_face_to_node(initializeFaceToNode())
  {
    AXOM_STATIC_ASSERT_MSG(
      TOPO == MIXED_SHAPE,
      "This constructor is only active for mixed topology meshes.");

    m_has_mixed_topology = true;
    initialize();
  }

  UnstructuredMesh(int ndims,
                   sidre::Group* group,
                   IndexType node_capacity = USE_DEFAULT,
                   IndexType cell_capacity = USE_DEFAULT,
                   IndexType connectivity_capacity = USE_DEFAULT)
    : UnstructuredMesh(ndims,
                       group,
                       "",
                       "",
                       node_capacity,
                       cell_capacity,
                       connectivity_capacity)
  { }

    /// @}

#endif /* AXOM_MINT_USE_SIDRE */

  /// @}

  /// \name Virtual methods
  /// @{

  /*!
   * \brief Destructor, deletes the MeshCoordinates and ConnectivityArray.
   */
  virtual ~UnstructuredMesh()
  {
    delete m_coordinates;
    m_coordinates = nullptr;

    delete m_cell_to_node;
    m_cell_to_node = nullptr;

    delete m_cell_to_face;
    m_cell_to_face = nullptr;

    delete m_face_to_cell;
    m_face_to_cell = nullptr;

    delete m_face_to_node;
    m_face_to_node = nullptr;
  }

  /// \name Cells
  /// @{

  /*!
   * \brief Return the number of cells in the mesh.
   */
  virtual IndexType getNumberOfCells() const final override
  {
    return m_cell_to_node->getNumberOfIDs();
  }

  /*!
   * \brief Return the capacity for cells.
   */
  virtual IndexType getCellCapacity() const final override
  {
    return m_cell_to_node->getIDCapacity();
  }

  /*!
   * \brief Return the type of the given cell.
   *
   * \param [in] cellID the ID of the cell in question, this parameter is
   *  ignored if TOPO == SINGLE_SHAPE. If TOPO == MIXED_SHAPE and no
   *  cellID is provided the returned type is UNDEFINED_CELL.
   *
   * \pre 0 <= cellID < getNumberOfCells()
   */
  virtual CellType getCellType(IndexType cellID = -1) const final override
  {
    return m_cell_to_node->getIDType(cellID);
  }

  /*!
   * \brief Return the number of nodes associated with the given cell.
   *
   * \param [in] cellID the ID of the cell in question, this parameter is
   *  ignored if TOPO == SINGLE_SHAPE.
   *
   * \pre 0 <= cellID < getNumberOfCells()
   */
  virtual IndexType getNumberOfCellNodes(IndexType cellID = 0) const final override
  {
    return m_cell_to_node->getNumberOfValuesForID(cellID);
  }

  /*!
   * \brief Copy the connectivity of the given cell into the provided buffer.
   *  The buffer must be of length at least getNumberOfCellNodes( cellID ).
   *
   * \param [in] cellID the ID of the cell in question.
   * \param [out] nodes the buffer into which the connectivity is copied.
   *
   * \return The number of nodes for the given cell.
   *
   * \pre nodes != nullptr
   * \pre 0 <= cellID < getNumberOfCells()
   */
  virtual IndexType getCellNodeIDs(IndexType cellID,
                                   IndexType* nodes) const final override
  {
    SLIC_ASSERT(nodes != nullptr);
    const IndexType n_nodes = getNumberOfCellNodes(cellID);
    std::memcpy(nodes, getCellNodeIDs(cellID), n_nodes * sizeof(IndexType));
    return n_nodes;
  }

  /*!
   * \brief Return the number of faces associated with the given cell.
   *
   * \param [in] cellID the ID of the cell in question.
   *
   * \note Codes must call initializeFaceConnectivity() before calling
   *       this method.
   */
  virtual IndexType getNumberOfCellFaces(IndexType cellID = 0) const final override
  {
    return getCellInfo(getCellType(cellID)).num_faces;
  }

  /*!
   * \brief Copy the face IDs of the given cell into the provided buffer.
   * The buffer must be of length at least getNumberOfCellFaces( cellID ).
   *
   * \param [in] cellID the ID of the cell in question
   * \param [out] faces the buffer into which the face IDs are copied.
   *
   * \return The number of faces for the given cell.
   *
   * \note Codes must call initializeFaceConnectivity() before calling
   *       this method.
   *
   * \pre faces != nullptr
   * \pre 0 <= cellID < getNumberOfCells()
   */
  virtual IndexType getCellFaceIDs(IndexType cellID,
                                   IndexType* faces) const final override
  {
    SLIC_ASSERT(faces != nullptr);
    const IndexType n_faces = getNumberOfCellFaces(cellID);
    std::memcpy(faces, getCellFaceIDs(cellID), n_faces * sizeof(IndexType));
    return n_faces;
  }

  /// @}

  /// \name Nodes
  /// @{

  /*!
   * \brief Return the number of nodes in the mesh.
   */
  virtual IndexType getNumberOfNodes() const final override
  {
    return m_coordinates->numNodes();
  }

  /*!
   * \brief Return the capacity for nodes.
   */
  virtual IndexType getNodeCapacity() const final override
  {
    return m_coordinates->capacity();
  }

  /*!
   * \brief Copy the coordinates of the given node into the provided buffer.
   *
   * \param [in] nodeID the ID of the node in question.
   * \param [in] coords the buffer to copy the coordinates into, of length at
   *  least getDimension().
   *
   * \pre 0 <= nodeID < getNumberOfNodes()
   * \pre coords != nullptr
   */
  virtual void getNode(IndexType nodeID, double* coords) const final override
  {
    m_coordinates->getCoordinates(nodeID, coords);
  }

  /*!
   * \brief Return a pointer to the array of nodal coordinates of the
   *  given dimension.
   *
   * \param [in] dim the dimension to return.
   *
   * \pre 0 <= dim < getDimension()
   */
  /// @{

  virtual double* getCoordinateArray(int dim) final override
  {
    return m_coordinates->getCoordinateArray(dim);
  }

  virtual const double* getCoordinateArray(int dim) const final override
  {
    return m_coordinates->getCoordinateArray(dim);
  }

  /// @}

  /// @}

  /// \name Faces
  /// @{

  /*!
   * \brief Return the number of faces in the mesh.
   *
   * \note Codes must call initializeFaceConnectivity() before calling
   *       this method.
   */
  virtual IndexType getNumberOfFaces() const final override
  {
    return m_face_to_cell->getNumberOfIDs();
  }

  /*!
   * \brief Return the capacity for faces.
   *
   * \note Codes must call initializeFaceConnectivity() before calling
   *       this method.
   */
  virtual IndexType getFaceCapacity() const final override
  {
    return m_face_to_cell->getIDCapacity();
  }

  /*!
   * \brief Return the type of the given face.
   *
   * \param [in] faceID the ID of the face in question.
   *
   * \note Codes must call initializeFaceConnectivity() before calling
   *       this method.
   */
  virtual CellType getFaceType(IndexType faceID) const final override
  {
    return m_face_to_node->getIDType(faceID);
  }

  /*!
   * \brief Return the number of nodes associated with the given face.
   *
   * \param [in] faceID the ID of the face in question.
   *
   * \note Codes must call initializeFaceConnectivity() before calling
   *       this method.
   */
  virtual IndexType getNumberOfFaceNodes(IndexType faceID = 0) const final override
  {
    return m_face_to_node->getNumberOfValuesForID(faceID);
  }

  /*!
   * \brief Copy the IDs of the nodes that compose the given face into the
   *  provided buffer.
   *
   * \param [in] faceID the ID of the face in question.
   * \param [out] nodes the buffer into which the node IDs are copied, must
   *  be of length at least getNumberOfFaceNodes().
   *
   * \return The number of nodes for the given face.
   *
   * A face with ID faceID will have a normal pointing outward from the
   * first cell it is adjacent to, as returned by
   * getFaceCellIDs(faceID, cellID1, cellID2).
   *
   * \pre nodes != nullptr
   * \pre 0 <= faceID < getNumberOfCells()
   */
  virtual IndexType getFaceNodeIDs(IndexType faceID,
                                   IndexType* nodes) const final override
  {
    SLIC_ASSERT(nodes != nullptr);
    const IndexType n_nodes = getNumberOfFaceNodes(faceID);
    std::memcpy(nodes, getFaceNodeIDs(faceID), n_nodes * sizeof(IndexType));
    return n_nodes;
  }

  /*!
   * \brief Copy the IDs of the cells adjacent to the given face into the
   *  provided indices.
   *
   * \param [in] faceID the ID of the face in question.
   * \param [out] cellIDOne the ID of the first cell.
   * \param [out] cellIDTwo the ID of the second cell.
   *
   * \note A face can be associated with one or two cells, depending on whether
   *  it is an external boundary face or interior face. By convention, if a face
   *  is an external boundary face, then only cellIDOne exists and cellIDTwo
   *  will be set to -1.
   *
   * \note Codes must call initializeFaceConnectivity() before calling
   *       this method.
   *
   * \pre 0 <= faceID < getNumberOfFaces()
   */
  virtual void getFaceCellIDs(IndexType faceID,
                              IndexType& cellIDOne,
                              IndexType& cellIDTwo) const final override
  {
    IndexType* faces = (*m_face_to_cell)[faceID];
    cellIDOne = faces[0];
    cellIDTwo = faces[1];
  }

  /// @}

  /// \name Edges
  /// @{

  /*!
   * \brief Return the number of edges in the mesh.
   */
  virtual IndexType getNumberOfEdges() const final override
  {
    SLIC_ERROR("NOT IMPLEMENTED!!!");
    return 0;
  }

  /*!
   * \brief Return the capacity for edges.
   */
  virtual IndexType getEdgeCapacity() const final override
  {
    SLIC_ERROR("NOT IMPLEMENTED!!!");
    return 0;
  }

  /// @}

  /*!
   * \brief Return true iff both the connectivity and coordinates are stored in
   *  external arrays.
   */
  virtual bool isExternal() const final override
  {
    bool connec_external = m_cell_to_node->isExternal();
    bool coords_external = m_coordinates->isExternal();

    if(connec_external != coords_external)
    {
      SLIC_WARNING("External state not consistent.");
      return false;
    }

    return connec_external;
  }

  /// @}

  /// \name Attribute get/set Methods
  /// @{

  /// \name Cells
  /// @{

  /*!
   * \brief Return the cell resize ratio.
   */
  double getCellResizeRatio() const { return m_cell_to_node->getResizeRatio(); }

  /*!
   * \brief Set the cell resize ratio.
   *
   * \param [in] ratio the new cell resize ratio.
   *
   * \post getCellResizeRatio() == ratio
   */
  void setCellResizeRatio(double ratio)
  {
    m_cell_to_node->setResizeRatio(ratio);
    m_mesh_fields[CELL_CENTERED]->setResizeRatio(ratio);
  }

  /*!
   * \brief Return the size of the connectivity array.
   */
  IndexType getCellNodesSize() const
  {
    return m_cell_to_node->getNumberOfValues();
  }

  /*!
   * \brief Return the capacity of the connectivity array.
   */
  IndexType getCellNodesCapacity() const
  {
    return m_cell_to_node->getValueCapacity();
  }

  /*!
   * \brief Resizes the cell connectivity array and cell-centered fields of this
   *  mesh instance to hold the specified number of cells.
   *
   * \param [in] cell_size the number of cells to resize to.
   *
   * \post getNumberOfCells() == cell_size
   */
  void resizeCells(IndexType cell_size)
  {
    IndexType connectivity_size =
      (hasMixedCellTypes()) ? USE_DEFAULT : getNumberOfCellNodes() * cell_size;
    m_cell_to_node->resize(cell_size, connectivity_size);
    m_mesh_fields[CELL_CENTERED]->resize(cell_size);
  }

  /*!
   * \brief Reserve space for the given number of cells.
   *
   * \param [in] cell_capacity the number of cells to reserve space for.
   * \param [in] connectivity_capacity the ammount of space to reserve in the
   *  connectivity array. Ignored if TOPO == SINGLE_SHAPE.
   *
   * \post getCellCapacity() >= cell_capacity
   */
  void reserveCells(IndexType cell_capacity,
                    IndexType connectivity_capacity = USE_DEFAULT)
  {
    m_cell_to_node->reserve(cell_capacity, connectivity_capacity);
    m_mesh_fields[CELL_CENTERED]->reserve(cell_capacity);
  }

  /*!
   * \brief Shrink the cell capacity to be equal to the number of cells.
   *
   * \post getCellCapacity() == getNumberOfCells()
   * \post getCellNodesCapacity() == getCellNodesSize()
   */
  void shrinkCells()
  {
    m_cell_to_node->shrink();
    m_mesh_fields[CELL_CENTERED]->shrink();
  }

  /// @}

  /// \name Nodes
  /// @{

  /*!
   * \brief Return the node resize ratio.
   */
  double getNodeResizeRatio() const { return m_coordinates->getResizeRatio(); }

  /*!
   * \brief Set the node resize ratio.
   *
   * \param [in] ratio the new node resize ratio.
   *
   * \post getNodeResizeRatio() == ratio
   */
  void setNodeResizeRatio(double ratio)
  {
    m_coordinates->setResizeRatio(ratio);
    m_mesh_fields[NODE_CENTERED]->setResizeRatio(ratio);
  }

  /*!
   * \brief Resizes the nodal coordinates and fields of this mesh instance to
   *  the specified number of nodes.
   *
   * \param [in] nodes_size the number of nodes to resize to.
   *
   * \post getNumberOfNodes() == nodes_size
   */
  void resizeNodes(IndexType nodes_size)
  {
    m_coordinates->resize(nodes_size);
    m_mesh_fields[NODE_CENTERED]->resize(nodes_size);
  }

  /*!
   * \brief Reserve space for the given number of nodes.
   *
   * \param [in] node_capacity the number of nodes to reserve space for.
   *
   * \post getNodeCapacity() >= node_capacity
   */
  void reserveNodes(IndexType node_capacity)
  {
    m_coordinates->reserve(node_capacity);
    m_mesh_fields[NODE_CENTERED]->reserve(node_capacity);
  }

  /*!
   * \brief Shrink the node capacity to be equal to the number of nodes.
   *
   * \post getNodeCapacity() == getNumberOfNodes()
   */
  void shrinkNodes()
  {
    m_coordinates->shrink();
    m_mesh_fields[NODE_CENTERED]->shrink();
  }

  /// @}

  /// \name Faces
  /// @{

  /*!
   * \brief Return the face resize ratio.
   */
  double getFaceResizeRatio() const
  {
    const double ratio = m_face_to_node->getResizeRatio();
    SLIC_WARNING_IF(m_face_to_cell->getResizeRatio() != ratio,
                    "Resize ratios are inconsistent");
    return ratio;
  }

  /*!
   * \brief Return the size of the connectivity array.
   */
  IndexType getFaceNodesSize() const
  {
    return m_face_to_node->getNumberOfValues();
  }

  /*!
   * \brief Return the capacity of the connectivity array.
   */
  IndexType getFaceNodesCapacity() const
  {
    return m_face_to_node->getValueCapacity();
  }

  /// @}

  /// \name Edges
  /// @{

  /*!
   * \brief Return the edge resize ratio.
   */
  double getEdgeResizeRatio() const
  {
    SLIC_ERROR("NOT IMPLEMENTED!!!");
    return 0.0;
  }

  /// @}

  /*!
   * \brief Resizes this mesh instance to the specified number of nodes & cells.
   *
   * \param [in] node_size the desired number of nodes
   * \param [in] cell_size the desired number of cells
   *
   * \note This method will also resize the node-centered and cell-centered
   *  fields accordingly.
   *
   * \post getNumberOfNodes() == nodes_size
   * \post getNumberOfCells() == cell_size
   *
   * \see resizeNodes()
   * \see resizeCells()
   */
  void resize(IndexType node_size, IndexType cell_size)
  {
    resizeNodes(node_size);
    resizeCells(cell_size);
  }

  /*!
   * \brief Reserve space for the given number of nodes and cells.
   *
   * \param [in] node_capacity the number of nodes to reserve space for.
   * \param [in] cell_capacity the number of cells to reserve space for.
   * \param [in] connectivity_capacity the ammount of space to reserve in the
   *  connectivity array. Ignored if TOPO == SINGLE_SHAPE.
   *
   * \post getNodeCapacity() >= node_capacity
   * \post getCellCapacity() >= cell_capacity
   */
  void reserve(IndexType node_capacity,
               IndexType cell_capacity,
               IndexType connectivity_capacity = USE_DEFAULT)
  {
    reserveNodes(node_capacity);
    reserveCells(cell_capacity, connectivity_capacity);
  }

  /*!
   * \brief Shrink the node capacity to be equal to the number of nodes and the
   *  cell capacity to be equal to the number of cells.
   *
   * \post getNodeCapacity() == getNumberOfNodes()
   * \post getCellCapacity() == getNumberOfCells()
   * \post getCellNodesCapacity() == getCellNodesSize()
   */
  void shrink()
  {
    shrinkNodes();
    shrinkCells();
  }

  /*!
   * \brief Return true iff the mesh holds no nodes and no cells.
   */
  bool empty() const
  {
    return m_coordinates->empty() && m_cell_to_node->empty();
  }

  /*!
   * \brief Return true iff both the connectivity and coordinates are stored in
   *  sidre.
   */
  bool isInSidre() const
  {
    bool connec_sidre = m_cell_to_node->isInSidre();
    bool coords_sidre = m_coordinates->isInSidre();

    if(connec_sidre != coords_sidre)
    {
      SLIC_WARNING("Sidre state not consistent.");
      return false;
    }

    return connec_sidre;
  }

  /// @}

  /// \name Data Access Methods
  /// @{

  /// \name Cells
  /// @{

  /*!
   * \brief Return a pointer to the connectivity of the given cell. The
   *  buffer is guarenteed to be of length at least
   *  getNumberOfCellNodes( cellID ).
   *
   * \param [in] cellID the ID of the cell in question.
   *
   * \pre 0 <= cellID < getNumberOfCells()
   */
  /// @{

  IndexType* getCellNodeIDs(IndexType cellID)
  {
    return (*m_cell_to_node)[cellID];
  }

  const IndexType* getCellNodeIDs(IndexType cellID) const
  {
    return (*m_cell_to_node)[cellID];
  }

  /*!
   * \brief Return a pointer to the faces of the given cell. The
   *  buffer is guarenteed to be of length at least
   *  getNumberOfCellFaces( cellID ).
   *
   * \param [in] cellID the ID of the cell in question.
   *
   * \note Codes must call initializeFaceConnectivity() before calling
   *       this method.
   *
   * \pre 0 <= cellID < getNumberOfCells()
   */
  /// @{

  IndexType* getCellFaceIDs(IndexType cellID)
  {
    return (*m_cell_to_face)[cellID];
  }

  const IndexType* getCellFaceIDs(IndexType cellID) const
  {
    return (*m_cell_to_face)[cellID];
  }

  /// @}

  /*!
   * \brief Return a pointer to the cell nodes array, of length
   *  getCellNodesSize().
   */
  /// @{

  IndexType* getCellNodesArray() { return m_cell_to_node->getValuePtr(); }

  const IndexType* getCellNodesArray() const
  {
    return m_cell_to_node->getValuePtr();
  }

  /// @}

  /*!
   * \brief Return a pointer to the cell nodes offset array, of length
   *  getNumberOfCells() + 1. Returns nullptr if
   *  TOPO == SINGLE_SHAPE.
   */
  /// @{

  IndexType* getCellNodesOffsetsArray()
  {
    return m_cell_to_node->getOffsetPtr();
  }

  const IndexType* getCellNodesOffsetsArray() const
  {
    return m_cell_to_node->getOffsetPtr();
  }

  /// @}

  /*!
   * \brief Return a pointer to the cell types array, of length
   *  getNumberOfCells(). Returns nullptr if
   *  TOPO == SINGLE_SHAPE.
   */
  /// @{

  CellType* getCellTypesArray() { return m_cell_to_node->getTypePtr(); }

  const CellType* getCellTypesArray() const
  {
    return m_cell_to_node->getTypePtr();
  }

  /// @}

  /*!
   * \brief Return a pointer to the cell faces array, of length
   *  getCellNodesSize().
   */
  /// @{

  IndexType* getCellFacesArray() { return m_cell_to_face->getValuePtr(); }

  const IndexType* getCellFacesArray() const
  {
    return m_cell_to_face->getValuePtr();
  }

  /// @}

  /*!
   * \brief Return a pointer to the cell faces offset array, of length
   *  getNumberOfCells() + 1.
   */
  /// @{

  IndexType* getCellFacesOffsetsArray()
  {
    return m_cell_to_face->getOffsetPtr();
  }

  const IndexType* getCellFacesOffsetsArray() const
  {
    return m_cell_to_face->getOffsetPtr();
  }

  /// @}

  /*!
   * \brief Append a cell to the mesh.
   *
   * \param [in] connec the connectivity of the new cell.
   * \param [in] type the type of the new cell, ignored if
   *  TOPO == SINGLE_SHAPE.
   *
   * \pre connec != nullptr
   */
  void appendCell(const IndexType* connec, CellType type = UNDEFINED_CELL)
  {
    IndexType n_values =
      (type == UNDEFINED_CELL) ? 0 : getCellInfo(type).num_nodes;
    m_cell_to_node->append(connec, n_values, type);
    m_mesh_fields[CELL_CENTERED]->resize(getNumberOfCells());
  }

  /*!
   * \brief Append multiple cells to the mesh.
   *
   * \param [in] connec the connectivity of the new cells.
   * \param [in] n_cells the number of cells to append.
   * \param [in] offsets the offsets array of the cells to append, ignored
   *  if TOPO == SINGLE_SHAPE.
   * \param [in] types the types array of the new cells, ignored if
   *  TOPO == SINGLE_SHAPE.
   *
   * \pre connec != nullptr
   * \pre n_cells >= 0
   */
  void appendCells(const IndexType* connec,
                   IndexType n_cells,
                   const IndexType* offsets = nullptr,
                   const CellType* types = nullptr)
  {
    m_cell_to_node->appendM(connec, n_cells, offsets, types);
    m_mesh_fields[CELL_CENTERED]->resize(getNumberOfCells());
  }

  /*!
   * \brief Insert a cell in to the mesh at the given position.
   *
   * \param [in] connec the connectivity of the new cell.
   * \param [in] ID the position to insert at.
   * \param [in] n_values the number of values in the connectivity, ignored
   *  if TOPO == SINGLE_SHAPE.
   * \param [in] type the type of the new cells, ignored if
   *  TOPO == SINGLE_SHAPE.
   *
   * \pre connec != nullptr
   * \pre 0 <= ID <= getNumberOfCells()
   */
  void insertCell(const IndexType* connec,
                  IndexType ID,
                  CellType type = UNDEFINED_CELL)
  {
    IndexType n_values =
      (type == UNDEFINED_CELL) ? 0 : getCellInfo(type).num_nodes;
    m_cell_to_node->insert(connec, ID, n_values, type);
    m_mesh_fields[CELL_CENTERED]->emplace(ID, 1);
  }

  /*!
   * \brief Insert multiple cells in to the mesh at the given position.
   *
   * \param [in] connec the connectivity of the new cells.
   * \param [in] start_ID the position to insert at.
   * \param [in] n_cells the number of cells to insert
   * \param [in] offsets the offsets array of the cells to append, ignored
   *  if TOPO == SINGLE_SHAPE.
   * \param [in] types the types array of the new cells, ignored if
   *  TOPO == SINGLE_SHAPE.
   *
   * \pre connec != nullptr
   * \pre 0 <= start_ID <= getNumberOfCells()
   */
  void insertCells(const IndexType* connec,
                   IndexType start_ID,
                   IndexType n_cells,
                   const IndexType* offsets = nullptr,
                   const CellType* types = nullptr)
  {
    m_cell_to_node->insertM(connec, start_ID, n_cells, offsets, types);
    m_mesh_fields[CELL_CENTERED]->emplace(start_ID, n_cells);
  }

  /// @}

  /// \name Nodes
  /// @{

  /*!
   * \brief Return the coordinate of the given dimension of the given node.
   *
   * \param [in] nodeID the ID of the node in question.
   * \param [in] dim the dimension to return.
   *
   * \pre 0 <= nodeID < getNumberOfNodes()
   * \pre 0 <= dim < getDimension()
   */
  double getNodeCoordinate(IndexType nodeID, int dim) const
  {
    return m_coordinates->getCoordinate(nodeID, dim);
  }

  /*!
   * \brief Appends a new node to the mesh.
   *
   * \param [in] x the first coordinate to append.
   * \param [in] y the second coordinate to append.
   * \param [in] z the third coordinate to append.
   *
   * \note Each method is valid only for the appropriate dimension of the mesh.
   */
  /// @{

  IndexType appendNode(double x)
  {
    IndexType n_index = m_coordinates->append(x);
    m_mesh_fields[NODE_CENTERED]->resize(getNumberOfNodes());
    return n_index;
  }

  IndexType appendNode(double x, double y)
  {
    IndexType n_index = m_coordinates->append(x, y);
    m_mesh_fields[NODE_CENTERED]->resize(getNumberOfNodes());
    return n_index;
  }

  IndexType appendNode(double x, double y, double z)
  {
    IndexType n_index = m_coordinates->append(x, y, z);
    m_mesh_fields[NODE_CENTERED]->resize(getNumberOfNodes());
    return n_index;
  }

  /// @}

  /*!
   * \brief Appends multiple nodes to the mesh.
   *
   * \param [in] coords pointer to the nodes to append, of length
   *  n * getDimension().
   * \param [in] n the number of nodes to append.
   *
   * \note coords is assumed to be in the array of structs format, ie
   *  coords = {x0, y0, z0, x1, y1, z1, ..., xn, yn, zn}.
   *
   * \pre coords != nullptr
   * \pre n >= 0
   */
  void appendNodes(const double* coords, IndexType n = 1)
  {
    m_coordinates->append(coords, n);
    m_mesh_fields[NODE_CENTERED]->resize(getNumberOfNodes());
  }

  /*!
   * \brief Appends new nodes to the mesh.
   *
   * \param [in] x array of the first coordinates to append, of length n.
   * \param [in] y array of the second coordinates to append, of length n.
   * \param [in] z array of the third coordinates to append, of length n.
   * \param [in] n the number of coordinates to append.
   *
   * \note The first method is only valid for 2D meshes while the second
   *  is only for 3D.
   * \pre x != nullptr
   * \pre y != nullptr
   * \pre z != nullptr
   * \pre n >= 0
   */
  /// @{

  void appendNodes(const double* x, const double* y, IndexType n)
  {
    m_coordinates->append(x, y, n);
    m_mesh_fields[NODE_CENTERED]->resize(getNumberOfNodes());
  }

  void appendNodes(const double* x, const double* y, const double* z, IndexType n)
  {
    m_coordinates->append(x, y, z, n);
    m_mesh_fields[NODE_CENTERED]->resize(getNumberOfNodes());
  }

  /// @}

  /*!
   * \brief Insert a node to the mesh.
   *
   * \param [in] nodeID the position to insert at.
   * \param [in] x the value of the first coordinate to insert.
   * \param [in] y the value of the second coordinate to insert.
   * \param [in] z the value of the third coordinate to insert.
   * \param [in] update_connectivity if true will update the connectivity so
   *  that all elements remain connected to the same coordinates as before.
   *
   * \note Each method is valid only for the appropriate dimension of the mesh.
   * \pre 0 <= nodeID <= getNumberOfNodes
   */
  /// @{

  void insertNode(IndexType nodeID, double x, bool update_connectivity = true)
  {
    m_coordinates->insert(nodeID, x);
    m_mesh_fields[NODE_CENTERED]->emplace(nodeID, 1);
    if(update_connectivity)
    {
      cellConnectivityUpdateInsert(nodeID, 1);
    }
  }

  void insertNode(IndexType nodeID,
                  double x,
                  double y,
                  bool update_connectivity = true)
  {
    m_coordinates->insert(nodeID, x, y);
    m_mesh_fields[NODE_CENTERED]->emplace(nodeID, 1);
    if(update_connectivity)
    {
      cellConnectivityUpdateInsert(nodeID, 1);
    }
  }

  void insertNode(IndexType nodeID,
                  double x,
                  double y,
                  double z,
                  bool update_connectivity = true)
  {
    m_coordinates->insert(nodeID, x, y, z);
    m_mesh_fields[NODE_CENTERED]->emplace(nodeID, 1);
    if(update_connectivity)
    {
      cellConnectivityUpdateInsert(nodeID, 1);
    }
  }

  /// @}

  /*!
   * \brief Inserts multiple nodes to the mesh.
   *
   * \param [in] coords pointer to the nodes to insert, of length
   *  n * getDimension().
   * \param [in] n the number of nodes to append.
   * \param [in] update_connectivity if true will update the connectivity so
   *  that all elements remain connected to the same coordinates as before.
   *
   * \note coords is assumed to be in the array of structs format, ie
   *  coords = {x0, y0, z0, x1, y1, z1, ..., xn, yn, zn}.
   *
   * \pre 0 <= nodeID <= getNumberOfNodes
   * \pre coords != nullptr
   * \pre n >= 0
   */
  void insertNodes(IndexType nodeID,
                   const double* coords,
                   IndexType n = 1,
                   bool update_connectivity = true)
  {
    m_coordinates->insert(nodeID, coords, n);
    m_mesh_fields[NODE_CENTERED]->emplace(nodeID, n);
    if(update_connectivity)
    {
      cellConnectivityUpdateInsert(nodeID, n);
    }
  }

  /*!
   * \brief Insert multiple nodes to the mesh.
   *
   * \param [in] nodeID the position to insert at.
   * \param [in] x the array of the first coordinates to insert.
   * \param [in] y the array of the second coordinates to insert.
   * \param [in] z the array of the third coordinates to insert.
   * \param [in] n the number of nodes to insert.
   * \param [in] update_connectivity if true will update the connectivity so
   *  that all elements remain connected to the same coordinates as before.
   *
   * \note The first method is only valid for 2D meshes while the second
   *  is only for 3D.
   * \pre 0 <= nodeID <= getNumberOfNodes
   * \pre x != nullptr
   * \pre y != nullptr if 2-D or 3-D
   * \pre z != nullptr if 3-D
   * \pre n >= 0
   */
  /// @{

  void insertNodes(IndexType nodeID,
                   const double* x,
                   const double* y,
                   IndexType n,
                   bool update_connectivity = true)
  {
    m_coordinates->insert(nodeID, x, y, n);
    m_mesh_fields[NODE_CENTERED]->emplace(nodeID, n);
    if(update_connectivity)
    {
      cellConnectivityUpdateInsert(nodeID, n);
    }
  }

  void insertNodes(IndexType nodeID,
                   const double* x,
                   const double* y,
                   const double* z,
                   IndexType n,
                   bool update_connectivity = true)
  {
    m_coordinates->insert(nodeID, x, y, z, n);
    m_mesh_fields[NODE_CENTERED]->emplace(nodeID, n);
    if(update_connectivity)
    {
      cellConnectivityUpdateInsert(nodeID, n);
    }
  }

  /// @}

  /// @}

  /// \name Faces
  /// @{

  /*!
   * \brief Sets up cell-face, face-cell, and face-node connectivity.
   *
   * \param [in] force re-initialize face-related connectivity, even if it
   *             has already been done.
   */
  bool initializeFaceConnectivity(bool force = false)
  {
    if(getDimension() == 1)
    {
      return true;
    }

    if(!force && getNumberOfFaces() > 0)
    {
      return true;
    }

    if(force)
    {
      m_cell_to_face->resize(0, 0);
      m_face_to_cell->resize(0, 0);
      m_face_to_node->resize(0, 0);
    }

    IndexType facecount = 0;
    IndexType* f2cdata = nullptr;
    IndexType* c2fdata = nullptr;
    IndexType* c2ndata = nullptr;
    IndexType* c2foffsets = nullptr;
    IndexType* f2ndata = nullptr;
    IndexType* f2noffsets = nullptr;
    CellType* f2ntypes = nullptr;

    bool retval = internal::initFaces(this,
                                      facecount,
                                      f2cdata,
                                      c2fdata,
                                      c2ndata,
                                      c2foffsets,
                                      f2ndata,
                                      f2noffsets,
                                      f2ntypes);

    if(retval)
    {
      // Copy in the face connectivity data.
      const IndexType cellCount = getNumberOfCells();
      m_cell_to_face->reserve(cellCount, c2foffsets[cellCount]);
      m_cell_to_face->appendM(c2fdata, cellCount, c2foffsets);

      m_face_to_cell->reserve(facecount);
      m_face_to_cell->appendM(f2cdata, facecount);

      m_face_to_node->reserve(facecount, f2noffsets[facecount]);
      m_face_to_node->appendM(f2ndata, facecount, f2noffsets, f2ntypes);
    }

    m_mesh_fields[FACE_CENTERED]->resize(getNumberOfFaces());

    delete[] f2cdata;
    delete[] c2fdata;
    delete[] c2ndata;
    delete[] c2foffsets;
    delete[] f2ndata;
    delete[] f2noffsets;
    delete[] f2ntypes;

    return retval;
  }

  /*!
   * \brief Return a pointer to the nodes of the given face. The
   *  buffer is guaranteed to be of length at least
   *  getNumberOfFaceNodes( faceID ).
   *
   * \param [in] faceID the ID of the face in question.
   *
   * \note Codes must call initializeFaceConnectivity() before calling
   *       this method.
   *
   * \pre 0 <= faceID < getNumberOfFaces()
   */
  /// @{

  IndexType* getFaceNodeIDs(IndexType faceID)
  {
    return (*m_face_to_node)[faceID];
  }

  const IndexType* getFaceNodeIDs(IndexType faceID) const
  {
    return (*m_face_to_node)[faceID];
  }

  /// @}

  /*!
   * \brief Return a pointer to the face nodes array, of length
   *  getFaceNodesSize().
   */
  /// @{

  IndexType* getFaceNodesArray() { return m_face_to_node->getValuePtr(); }

  const IndexType* getFaceNodesArray() const
  {
    return m_face_to_node->getValuePtr();
  }

  /// @}

  /*!
   * \brief Return a pointer to the face nodes offset array, of length
   *  getNumberOfFaces() + 1.
   */
  /// @{

  IndexType* getFaceNodesOffsetsArray()
  {
    return m_face_to_node->getOffsetPtr();
  }

  const IndexType* getFaceNodesOffsetsArray() const
  {
    return m_face_to_node->getOffsetPtr();
  }

  /// @}

  /*!
   * \brief Return a pointer to the face cells array, of length
   *  2 * getNumberOfFaces().
   */
  /// @{

  IndexType* getFaceCellsArray() { return m_face_to_cell->getValuePtr(); }

  const IndexType* getFaceCellsArray() const
  {
    return m_face_to_cell->getValuePtr();
  }

  /// @}

  /// @}

  /// @}

private:
  /*! \brief Construct and fill the cell-to-face connectivity. */
  void buildCellFaceConnectivity(IndexType* c2fdata, IndexType* c2foffsets);

  /*! \brief Return a new empty CellToFaceConnectivity instance. */
  CellToFaceConnectivity* initializeCellToFace(
    CellType cell_type = UNDEFINED_CELL) const;

  /*! \brief Return a new empty FaceToNodeConnectivity instance. */
  FaceToNodeConnectivity* initializeFaceToNode(
    CellType cell_type = UNDEFINED_CELL) const;

  /*!
   * \brief Update the connectivity given an nodal insert at position pos of
   *  length n.
   *
   * \param [in] pos the position of the insert.
   * \param [in] n the length of the insert.
   */
  void cellConnectivityUpdateInsert(IndexType pos, IndexType n)
  {
    SLIC_ASSERT(0 <= pos && pos < getNumberOfNodes());

    const IndexType n_values = getCellNodesSize();
    IndexType* values = getCellNodesArray();
    SLIC_ASSERT(n_values == 0 || values != nullptr);

    for(IndexType i = 0; i < n_values; ++i)
    {
      if(values[i] >= pos)
      {
        values[i] += n;
      }
    }
  }

  /*!
   * \brief Performs common initialization.
   */
  void initialize()
  {
    m_explicit_coords = true;
    m_explicit_connectivity = true;
    m_mesh_fields[NODE_CENTERED]->setResizeRatio(getNodeResizeRatio());
    m_mesh_fields[CELL_CENTERED]->setResizeRatio(getCellResizeRatio());
    m_mesh_fields[FACE_CENTERED]->setResizeRatio(getFaceResizeRatio());

    m_mesh_fields[NODE_CENTERED]->reserve(getNodeCapacity());
    m_mesh_fields[CELL_CENTERED]->reserve(getCellCapacity());
    m_mesh_fields[FACE_CENTERED]->reserve(getFaceCapacity());

    m_mesh_fields[NODE_CENTERED]->resize(getNumberOfNodes());
    m_mesh_fields[CELL_CENTERED]->resize(getNumberOfCells());
    m_mesh_fields[FACE_CENTERED]->resize(getNumberOfFaces());
  }

  MeshCoordinates* m_coordinates;

  /*! \brief The nodes for each cell */
  CellToNodeConnectivity* m_cell_to_node;

  /*! \brief Each cell's faces */
  CellToFaceConnectivity* m_cell_to_face;

  /*! \brief Each face's cells */
  FaceToCellConnectivity* m_face_to_cell;

  /*! \brief Each face's nodes */
  FaceToNodeConnectivity* m_face_to_node;

  DISABLE_COPY_AND_ASSIGNMENT(UnstructuredMesh);
  DISABLE_MOVE_AND_ASSIGNMENT(UnstructuredMesh);
};

/*!
 * \brief Return a new CellToFaceConnectivity. This instance is for
 *  SINGLE_SHAPE meshes. If the mesh has cell type VERTEX or SEGMENT then the
 *  stride of the returned connectivity is 1 instead of 0 which is the proper
 *  number of cell faces, this avoids an error with 0 stride.
 *
 * \param [in] cell_type the cell_type of the mesh.
 */
template <>
inline UnstructuredMesh<SINGLE_SHAPE>::CellToFaceConnectivity*
UnstructuredMesh<SINGLE_SHAPE>::initializeCellToFace(CellType cell_type) const
{
  IndexType num_faces = getCellInfo(cell_type).num_faces;
  if(num_faces == 0)
  {
    num_faces = 1;
  }

  return new UnstructuredMesh<SINGLE_SHAPE>::CellToFaceConnectivity(num_faces, 0);
}

/*!
 * \brief Return a new CellToFaceConnectivity. This instance is for
 *  MIXED_SHAPE meshes.
 */
template <>
inline UnstructuredMesh<MIXED_SHAPE>::CellToFaceConnectivity*
  UnstructuredMesh<MIXED_SHAPE>::initializeCellToFace(CellType) const
{
  return new UnstructuredMesh<MIXED_SHAPE>::CellToFaceConnectivity(UNDEFINED_CELL,
                                                                   0,
                                                                   0);
}

/*!
 * \brief Return a new FaceToNodeConnectivity. This instance is for
 *  SINGLE_SHAPE meshes. If the mesh has cell type VERTEX or SEGMENT then the
 *  cell type of the returned ConnectivityArray is UNDEFINED_CELL with stride 1.
 *
 * \param [in] cell_type the cell_type of the mesh.
 */
template <>
inline UnstructuredMesh<SINGLE_SHAPE>::FaceToNodeConnectivity*
UnstructuredMesh<SINGLE_SHAPE>::initializeFaceToNode(CellType cell_type) const
{
  const CellType face_type = getCellInfo(cell_type).face_types[0];
  if(face_type == UNDEFINED_CELL)
  {
    return new UnstructuredMesh<SINGLE_SHAPE>::FaceToNodeConnectivity(1, 0);
  }

  return new UnstructuredMesh<SINGLE_SHAPE>::FaceToNodeConnectivity(face_type, 0);
}

/*!
 * \brief Return a new FaceToNodeConnectivity. This instance is for
 *  MIXED_SHAPE meshes.
 */
template <>
inline UnstructuredMesh<MIXED_SHAPE>::FaceToNodeConnectivity*
  UnstructuredMesh<MIXED_SHAPE>::initializeFaceToNode(CellType) const
{
  return new UnstructuredMesh<MIXED_SHAPE>::FaceToNodeConnectivity(0, 0);
}

} /* namespace mint */
} /* namespace axom */

#endif /* MINT_UNSTRUCTUREDMESH_HPP_ */