From: David Wells <drwells@email.unc.edu>
Date: Wed, 9 Sep 2020 17:03:14 +0000 (-0400)
Subject: Implement ExodusII in GridIn.
X-Git-Tag: v9.3.0-rc1~917^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10899%2Fhead;p=dealii.git

Implement ExodusII in GridIn.

The tests were all generated by hand. For future reference I included the code here.

Here is the basic 2D test (five quadrilaterals):

#include <exodusII.h>

#include <algorithm>
#include <cassert>
#include <cstring>

// Create a disk mesh with five quadrilaterals.
//
//  ____________________________
// |\                          /|
// | \                        / |
// |  \           3          /  |
// |   \                    /   |
// |    \                  /    |
// |     \                /     |
// |      \______________/      |
// |      |              |      |
// |      |              |      |
// |      |              |      |
// |   4  |      5       |   2  |
// |      |              |      |
// |      |______________|      |
// |      /              \      |
// |     /                \     |
// |    /                  \    |
// |   /                    \   |
// |  /          1           \  |
// | /                        \ |
// |/                          \|
// -----------------------------
//
// element 1 has external side sets 0, 1
// element 2 has external side sets 0, 1, 10
// element 3 has no external side sets
// element 4 has no external side set 0
//
// The inner element is in block 10 - the rest are in block 1.
int main()
{
    const int n_elems = 5;
    const int n_nodes = 8;
    int elements_1[n_elems - 1][4] = {
    {1, 2, 3, 4},
    {3, 2, 8, 7},
    {5, 7, 8, 6},
    {1, 4, 5, 6}};
    int elements_10[1][4] = {{4, 3, 7, 5}};

    double xs[n_nodes] = {0, 7, 5, 2, 2, 0, 5, 7};
    double ys[n_nodes] = {0, 0, 2, 2, 5, 7, 5, 7};

    int sizeof_data = sizeof(double);
    int sizeof_geom = sizeof(double);
    const int exid = ex_create("five.e", EX_CLOBBER, &sizeof_data, &sizeof_geom);

    const int dim = 2;
    const int n_blocks = 2;
    const int n_node_sets = 0;
    const int n_side_sets = 3;
    int ierr = ex_put_init(exid, "mytitle", dim, n_nodes, n_elems, n_blocks,
                           n_node_sets, n_side_sets);
    assert(ierr == 0);

    // nodes:
    ierr = ex_put_coord(exid, xs, ys, nullptr);
    assert(ierr == 0);

    // element blocks:
    ierr = ex_put_elem_block(exid, 1, "quad", 4, 4, 0);
    assert(ierr == 0);
    ierr = ex_put_elem_block(exid, 10, "quad", 1, 4, 0);
    assert(ierr == 0);

    ierr = ex_put_elem_conn(exid, 10, elements_10);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 1, elements_1);
    assert(ierr == 0);

    // side sets:
    ierr = ex_put_side_set_param(exid, 0, 3, 0);
    assert(ierr == 0);
    int sideset_0_elems[] = {1, 2, 4};
    int sideset_0_sides[] = {1, 2, 4};
    ierr = ex_put_side_set(exid, 0, sideset_0_elems, sideset_0_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 1, 2, 0);
    assert(ierr == 0);
    int sideset_1_elems[] = {1, 2};
    int sideset_1_sides[] = {1, 2};
    ierr = ex_put_side_set(exid, 1, sideset_1_elems, sideset_1_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 10, 1, 0);
    assert(ierr == 0);
    int sideset_10_elems[] = {2};
    int sideset_10_sides[] = {2};
    ierr = ex_put_side_set(exid, 10, sideset_10_elems, sideset_10_sides);
    assert(ierr == 0);

    ierr = ex_close(exid);
    assert(ierr == 0);
}

A 3D mixed mesh:

#include <exodusII.h>

#include <algorithm>
#include <cassert>
#include <cstring>

int main()
{
    const int n_elems = 4;
    const int n_nodes = 13;
    int elements_1[1][8] = {{1, 2, 3, 4, 5, 6, 7, 8}};
    int elements_2[1][5] = {{5, 6, 7, 8, 9}};
    int elements_3[1][4] = {{6, 10, 7, 9}};
    int elements_4[1][6] = {{10, 6, 7, 13, 11, 12}};

    double xs[n_nodes] = {0, 3, 3, 0, 0, 3, 3, 0, 1, 3, 8, 8, 8};
    double ys[n_nodes] = {0, 0, 3, 3, 0, 0, 3, 3, 1, 1, 0, 3, 1};
    double zs[n_nodes] = {0, 0, 0, 0, 3, 3, 3, 3, 5, 5, 3, 3, 5};

    int sizeof_data = sizeof(double);
    int sizeof_geom = sizeof(double);
    const int exid = ex_create("four.e", EX_CLOBBER, &sizeof_data, &sizeof_geom);

    const int dim = 3;
    const int n_blocks = n_elems;
    const int n_node_sets = 0;
    const int n_side_sets = 0;
    int ierr = ex_put_init(exid, "mytitle", dim, n_nodes, n_elems, n_blocks,
                           n_node_sets, n_side_sets);
    assert(ierr == 0);

    // nodes:
    ierr = ex_put_coord(exid, xs, ys, zs);
    assert(ierr == 0);

    // block 0 (hex):
    ierr = ex_put_elem_block(exid, 1, "hex", 1, 8, 0);
    assert(ierr == 0);
    // block 1 (pyramid):
    ierr = ex_put_elem_block(exid, 10, "pyramid", 1, 5, 0);
    assert(ierr == 0);
    // block 2 (tet):
    ierr = ex_put_elem_block(exid, 100, "tet", 1, 4, 0);
    assert(ierr == 0);
    // block 3 (wedge):
    ierr = ex_put_elem_block(exid, 1000, "wedge", 1, 6, 0);
    assert(ierr == 0);

    // connectivity:
    ierr = ex_put_elem_conn(exid, 1, elements_1);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 10, elements_2);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 100, elements_3);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 1000, elements_4);
    assert(ierr == 0);

    ierr = ex_close(exid);
    assert(ierr == 0);
}

Check that we can truncate higher-order elements correctly (like PYRAMID14 ->
our five-point pyramid)

#include <exodusII.h>

#include <algorithm>
#include <cassert>
#include <cstring>

int main()
{
    const int n_elems = 4;
    const int n_nodes = 29;
    int elements_1[1][20] = {{1, 2, 3, 4, 5, 6, 7, 8, 17, 14, 15, 16, 21, 18, 19, 20, 25, 22, 23, 24}};
    // int elements_1[1][8] = {{1, 2, 3, 4, 5, 6, 7, 8}};
    int elements_2[1][13] = {{5, 6, 7, 8, 9, 25, 22, 23, 24, 26, 27, 28, 29}};
    int elements_3[1][4] = {{6, 10, 7, 9}};
    int elements_4[1][6] = {{10, 6, 7, 13, 11, 12}};

    double xs[n_nodes] = {0, 3, 3, 0, 0, 3, 3, 0, 1, 3, 8, 8, 8, 3, 1.5, 0, 1.5, 3, 3, 0, 0, 3, 1.5, 0, 1.5, 0.5, 2, 2, 0.5};
    double ys[n_nodes] = {0, 0, 3, 3, 0, 0, 3, 3, 1, 1, 0, 3, 1, 1.5, 3, 1.5, 0, 0, 3, 3, 0, 1.5, 3, 1.5, 0, 0.5, 0.5, 2, 2};
    double zs[n_nodes] = {0, 0, 0, 0, 3, 3, 3, 3, 5, 5, 3, 3, 5, 0, 0, 0, 0, 1.5, 1.5, 1.5, 1.5, 3, 3, 3, 3, 4, 4, 4, 4};

    int sizeof_data = sizeof(double);
    int sizeof_geom = sizeof(double);
    const int exid = ex_create("four-p.e", EX_CLOBBER, &sizeof_data, &sizeof_geom);

    const int dim = 3;
    const int n_blocks = n_elems;
    const int n_node_sets = 0;
    const int n_side_sets = 0;
    int ierr = ex_put_init(exid, "mytitle", dim, n_nodes, n_elems, n_blocks,
                           n_node_sets, n_side_sets);
    assert(ierr == 0);

    // nodes:
    ierr = ex_put_coord(exid, xs, ys, zs);
    assert(ierr == 0);

    // block 0 (hex):
    ierr = ex_put_elem_block(exid, 1, "hex20", 1, 20, 0);
    assert(ierr == 0);
    // block 1 (pyramid):
    ierr = ex_put_elem_block(exid, 10, "pyramid13", 1, 13, 0);
    assert(ierr == 0);
    // block 2 (tet):
    ierr = ex_put_elem_block(exid, 100, "tet", 1, 4, 0);
    assert(ierr == 0);
    // block 3 (wedge):
    ierr = ex_put_elem_block(exid, 1000, "wedge", 1, 6, 0);
    assert(ierr == 0);

    // connectivity:
    ierr = ex_put_elem_conn(exid, 1, elements_1);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 10, elements_2);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 100, elements_3);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 1000, elements_4);
    assert(ierr == 0);

    ierr = ex_close(exid);
    assert(ierr == 0);
}

Get things working with a nonzero codimension:

#include <exodusII.h>

#include <algorithm>
#include <cassert>
#include <cstring>

// Create a disk mesh with five quadrilaterals.
//
//  ____________________________
// |\                          /|
// | \                        / |
// |  \           3          /  |
// |   \                    /   |
// |    \                  /    |
// |     \                /     |
// |      \______________/      |
// |      |              |      |
// |      |              |      |
// |      |              |      |
// |   4  |      5       |   2  |
// |      |              |      |
// |      |______________|      |
// |      /              \      |
// |     /                \     |
// |    /                  \    |
// |   /                    \   |
// |  /          1           \  |
// | /                        \ |
// |/                          \|
// -----------------------------
//
// element 1 has external side sets 0, 1
// element 2 has external side sets 0, 1, 10
// element 3 has no external side sets
// element 4 has no external side set 0
//
// The inner element is in block 10 - the rest are in block 1.
int main()
{
    const int n_elems = 5;
    const int n_nodes = 8;
    int elements_1[n_elems - 1][4] = {
    {1, 2, 3, 4},
    {1, 4, 5, 6},
    {2, 8, 7, 3},
    {6, 5, 7, 8}};
    int elements_10[1][4] = {{4, 3, 7, 5}};

    double xs[n_nodes] = {0, 7, 5, 2, 2, 0, 5, 7};
    double ys[n_nodes] = {0, 0, 2, 2, 5, 7, 5, 7};

    int sizeof_data = sizeof(double);
    int sizeof_geom = sizeof(double);
    const int exid = ex_create("five-3d.e", EX_CLOBBER, &sizeof_data, &sizeof_geom);

    const int dim = 3;
    const int n_blocks = 2;
    const int n_node_sets = 0;
    const int n_side_sets = 3;
    int ierr = ex_put_init(exid, "mytitle", dim, n_nodes, n_elems, n_blocks,
                           n_node_sets, n_side_sets);
    assert(ierr == 0);

    // nodes:
    ierr = ex_put_coord(exid, xs, ys, xs);
    assert(ierr == 0);

    // element blocks:
    ierr = ex_put_elem_block(exid, 1, "quad", 4, 4, 0);
    assert(ierr == 0);
    ierr = ex_put_elem_block(exid, 10, "quad", 1, 4, 0);
    assert(ierr == 0);

    ierr = ex_put_elem_conn(exid, 10, elements_10);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 1, elements_1);
    assert(ierr == 0);

    // side sets:
    ierr = ex_put_side_set_param(exid, 0, 3, 0);
    assert(ierr == 0);
    int sideset_0_elems[] = {1, 2, 4};
    int sideset_0_sides[] = {3, 4, 6};
    ierr = ex_put_side_set(exid, 0, sideset_0_elems, sideset_0_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 1, 2, 0);
    assert(ierr == 0);
    int sideset_1_elems[] = {1, 2};
    int sideset_1_sides[] = {3, 4};
    ierr = ex_put_side_set(exid, 1, sideset_1_elems, sideset_1_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 10, 1, 0);
    assert(ierr == 0);
    int sideset_10_elems[] = {2};
    int sideset_10_sides[] = {4};
    ierr = ex_put_side_set(exid, 10, sideset_10_elems, sideset_10_sides);
    assert(ierr == 0);

    ierr = ex_close(exid);
    assert(ierr == 0);
}

Check sidesets in 3D:

#include <exodusII.h>

#include <algorithm>
#include <cassert>
#include <cstring>

// Like writer-2 but with side sets.

int main()
{
    const int n_elems = 4;
    const int n_nodes = 13;
    int elements_1[1][8] = {{1, 2, 3, 4, 5, 6, 7, 8}};
    int elements_2[1][5] = {{5, 6, 7, 8, 9}};
    int elements_3[1][4] = {{6, 10, 7, 9}};
    int elements_4[1][6] = {{10, 6, 7, 13, 11, 12}};

    double xs[n_nodes] = {0, 3, 3, 0, 0, 3, 3, 0, 1, 3, 8, 8, 8};
    double ys[n_nodes] = {0, 0, 3, 3, 0, 0, 3, 3, 1, 1, 0, 3, 1};
    double zs[n_nodes] = {0, 0, 0, 0, 3, 3, 3, 3, 5, 5, 3, 3, 5};

    int sizeof_data = sizeof(double);
    int sizeof_geom = sizeof(double);
    const int exid = ex_create("four-sidesets.e", EX_CLOBBER, &sizeof_data, &sizeof_geom);

    const int dim = 3;
    const int n_blocks = n_elems;
    const int n_node_sets = 0;
    const int n_side_sets = 5;
    int ierr = ex_put_init(exid, "mytitle", dim, n_nodes, n_elems, n_blocks,
                           n_node_sets, n_side_sets);
    assert(ierr == 0);

    // nodes:
    ierr = ex_put_coord(exid, xs, ys, zs);
    assert(ierr == 0);

    // block 0 (hex):
    ierr = ex_put_elem_block(exid, 1, "hex", 1, 8, 0);
    assert(ierr == 0);
    // block 1 (pyramid):
    ierr = ex_put_elem_block(exid, 10, "pyramid", 1, 5, 0);
    assert(ierr == 0);
    // block 2 (tet):
    ierr = ex_put_elem_block(exid, 100, "tet", 1, 4, 0);
    assert(ierr == 0);
    // block 3 (wedge):
    ierr = ex_put_elem_block(exid, 1000, "wedge", 1, 6, 0);
    assert(ierr == 0);

    // connectivity:
    ierr = ex_put_elem_conn(exid, 1, elements_1);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 10, elements_2);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 100, elements_3);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 1000, elements_4);
    assert(ierr == 0);

    // side sets:
    ierr = ex_put_side_set_param(exid, 1, 5, 0);
    assert(ierr == 0);
    int sideset_1_elems[] = {1, 1, 2, 3, 4};
    int sideset_1_sides[] = {1, 5, 1, 1, 1};
    ierr = ex_put_side_set(exid, 1, sideset_1_elems, sideset_1_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 2, 3, 0);
    assert(ierr == 0);
    int sideset_2_elems[] = {1, 3, 4};
    int sideset_2_sides[] = {2, 2, 2};
    ierr = ex_put_side_set(exid, 2, sideset_2_elems, sideset_2_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 3, 3, 0);
    assert(ierr == 0);
    int sideset_3_elems[] = {1, 2, 4};
    int sideset_3_sides[] = {3, 3, 3};
    ierr = ex_put_side_set(exid, 3, sideset_3_elems, sideset_3_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 4, 2, 0);
    assert(ierr == 0);
    int sideset_4_elems[] = {1, 2};
    int sideset_4_sides[] = {4, 4};
    ierr = ex_put_side_set(exid, 4, sideset_4_elems, sideset_4_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 5, 4, 0);
    assert(ierr == 0);
    int sideset_5_elems[] = {1, 1, 2, 4};
    int sideset_5_sides[] = {4, 5, 4, 5};
    ierr = ex_put_side_set(exid, 5, sideset_5_elems, sideset_5_sides);
    assert(ierr == 0);

    ierr = ex_close(exid);
    assert(ierr == 0);
}
---

diff --git a/doc/news/changes/major/20200909DavidWells b/doc/news/changes/major/20200909DavidWells
new file mode 100644
index 0000000000..5664f10da5
--- /dev/null
+++ b/doc/news/changes/major/20200909DavidWells
@@ -0,0 +1,4 @@
+New: The GridIn class can now read ExodusII files when deal.II is configured
+with Trilinos and SEACAS.
+<br>
+(David Wells, 2020/09/09)
diff --git a/include/deal.II/base/exceptions.h b/include/deal.II/base/exceptions.h
index 588fa60d9d..73d7f8b139 100644
--- a/include/deal.II/base/exceptions.h
+++ b/include/deal.II/base/exceptions.h
@@ -1148,6 +1148,15 @@ namespace StandardExceptions
 #endif
   //@}
 
+  /**
+   * This function requires support for the Exodus II library.
+   */
+  DeclExceptionMsg(
+    ExcNeedsExodusII,
+    "You are attempting to use functionality that is only available if deal.II "
+    "was configured to use Trilinos' SEACAS library (which provides ExodusII), "
+    "but cmake did not find find a valid SEACAS library.");
+
 #ifdef DEAL_II_WITH_MPI
   /**
    * Exception for MPI errors. This exception is only defined if
@@ -1181,6 +1190,40 @@ namespace StandardExceptions
     const int error_code;
   };
 #endif // DEAL_II_WITH_MPI
+
+
+
+#ifdef DEAL_II_TRILINOS_WITH_SEACAS
+  /**
+   * Exception for ExodusII errors. This exception is only defined if
+   * <code>deal.II</code> is compiled with SEACAS support, which is available
+   * through Trilinos. This function should be used with the convenience macro
+   * AssertThrowExodusII.
+   *
+   * @ingroup Exceptions
+   */
+  class ExcExodusII : public ExceptionBase
+  {
+  public:
+    /**
+     * Constructor.
+     *
+     * @param error_code The error code returned by an ExodusII function.
+     */
+    ExcExodusII(const int error_code);
+
+    /**
+     * Print a description of the error to the given stream.
+     */
+    virtual void
+    print_info(std::ostream &out) const override;
+
+    /**
+     * Store the error code.
+     */
+    const int error_code;
+  };
+#endif // DEAL_II_TRILINOS_WITH_SEACAS
 } /*namespace StandardExceptions*/
 
 
@@ -1888,6 +1931,28 @@ namespace internal
 
 #endif
 
+#ifdef DEAL_II_TRILINOS_WITH_SEACAS
+/**
+ * Assertion that checks that the error code produced by calling an ExodusII
+ * routine is equal to EX_NOERR (which is zero).
+ *
+ * @note This and similar macro names are examples of preprocessor definitions
+ * in the deal.II library that are not prefixed by a string that likely makes
+ * them unique to deal.II. As a consequence, it is possible that other
+ * libraries your code interfaces with define the same name, and the result
+ * will be name collisions (see
+ * https://en.wikipedia.org/wiki/Name_collision). One can <code>\#undef</code>
+ * this macro, as well as all other macros defined by deal.II that are not
+ * prefixed with either <code>DEAL</code> or <code>deal</code>, by including
+ * the header <code>deal.II/base/undefine_macros.h</code> after all other
+ * deal.II headers have been included.
+ *
+ * @ingroup Exceptions
+ */
+#  define AssertThrowExodusII(error_code) \
+    AssertThrow(error_code == 0, ExcExodusII(error_code));
+#endif // DEAL_II_TRILINOS_WITH_SEACAS
+
 using namespace StandardExceptions;
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/include/deal.II/base/undefine_macros.h b/include/deal.II/base/undefine_macros.h
index 9f0dd01986..3d5b9a73c8 100644
--- a/include/deal.II/base/undefine_macros.h
+++ b/include/deal.II/base/undefine_macros.h
@@ -82,6 +82,10 @@
 #  undef AssertThrowMPI
 #endif // #ifdef AssertThrowMPI
 
+#ifdef AssertThrowExodusII
+#  undef AssertThrowExodusII
+#endif // #ifdef AssertThrowExodusII
+
 #ifdef AssertVectorVectorDimension
 #  undef AssertVectorVectorDimension
 #endif // #ifdef AssertVectorVectorDimension
diff --git a/include/deal.II/grid/grid_in.h b/include/deal.II/grid/grid_in.h
index 56bacfe3fd..03f7a9e5dc 100644
--- a/include/deal.II/grid/grid_in.h
+++ b/include/deal.II/grid/grid_in.h
@@ -334,6 +334,8 @@ public:
     vtu,
     /// Use read_assimp()
     assimp,
+    /// Use read_exodusii()
+    exodusii,
   };
 
   /**
@@ -567,6 +569,83 @@ public:
               const double       tol                              = 1e-12,
               const bool         ignore_unsupported_element_types = true);
 
+  /**
+   * A structure containing some of the information provided by ExodusII that
+   * doesn't have a direct representation in the Triangulation object.
+   *
+   * @note This struct exists to enable forward compatibility with future
+   * versions of read_exodusii that may provide additional output data, but for
+   * now it has a single field.
+   */
+  struct ExodusIIData
+  {
+    /**
+     * A vector containing a mapping from deal.II boundary ids (or manifold ids)
+     * to the provided ExodusII sideset ids.
+     */
+    std::vector<std::vector<int>> id_to_sideset_ids;
+  };
+
+  /**
+   * Read in a mesh stored in the ExodusII file format.
+   *
+   * ExodusII is a feature-rich file format that supports many more features
+   * (like node sets, finite element fields, quality assurance data, and more)
+   * than most other grid formats supported by this class. Many of these
+   * features do not have equivalent representations in deal.II and are
+   * therefore not supported (for example, deal.II does not assign degrees of
+   * freedom directly to nodes, so data stored in a nodal format is not loaded
+   * by this function). At the current time only the following information is
+   * extracted from the input file:
+   *
+   * <ol>
+   *   <li>Block ids: the block id of an element is loaded as its material
+   *     id.</li>
+   *   <li>Elements and vertices: the core geometric information stored in the
+   *     ExodusII file populates the attached Triangulation object. Higher-order
+   *     elements are automatically truncated to lower-order elements since
+   *     deal.II does not support this feature (e.g., there is no equivalent to
+   *     the <code>QUAD9</code> element in deal.II since all quadrilaterals have
+   *     four vertices and additional geometric information is either stored in
+   *     a Manifold or something like MappingQEulerian).</li>
+   *   <li>Sideset ids: these are interpreted as boundary ids or manifold ids
+   *     (see the note on the output value below). An error will occur if you
+   *     attempt to read an ExodusII file that assigns a sideset id to an
+   *     internal face boundary id.</li>
+   * </ol>
+   *
+   * Sideset ids are not translated for Triangulations with nonzero codimension
+   * since those Triangulations do not support the setting of boundary ids.
+   *
+   * @param filename The name of the file to read from.
+   *
+   * @param apply_all_indicators_to_manifolds Boolean determining if the sideset
+   * ids should be interpreted as manifold ids or boundary ids. The default
+   * value is <tt>false</tt>, i.e., treat all sideset ids as boundary ids. If
+   * your mesh sets sideset ids on internal faces then it will be necessary to
+   * set this argument to <code>true</code> and then do some postprocessing to
+   * set the boundary ids correctly.
+   *
+   * @return This function returns a struct containing some extra data stored by
+   * the ExodusII file that cannot be loaded into a Triangulation - see
+   * ExodusIIData for more information.
+   *
+   * A cell face in ExodusII can be in an arbitrary number of sidesets (i.e., it
+   * can have an arbitrary number of sideset ids) - however, a boundary cell
+   * face in deal.II has exactly one boundary id. All boundary faces that are
+   * not in a sideset are given the (default) boundary id of $0$. This function
+   * then groups sidesets together into unique sets and gives each one a
+   * boundary id. For example: Consider a single-quadrilateral mesh whose left
+   * side has no sideset id, right side has sideset ids $0$ and $1$, and whose
+   * bottom and top sides have sideset ids of $0$. The left face will have a
+   * boundary id of $0$, the top and bottom faces boundary ids of $1$, and the
+   * right face a boundary id of $2$. Hence the vector returned by this function
+   * in that case will be $\{\{\}, \{0\}, \{0, 1\}\}$.
+   */
+  ExodusIIData
+  read_exodusii(const std::string &filename,
+                const bool         apply_all_indicators_to_manifolds = false);
+
   /**
    * Return the standard suffix for a file in this format.
    */
diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index 017607e8f2..32d1c49b72 100644
--- a/include/deal.II/grid/reference_cell.h
+++ b/include/deal.II/grid/reference_cell.h
@@ -111,6 +111,10 @@ namespace ReferenceCell
        {{X, X, X, Type::Tri, Type::Quad, X, X, X, X}},
        // dim 3
        {{X, X, X, X, Type::Tet, Type::Pyramid, Type::Wedge, X, Type::Hex}}}};
+    Assert(table[dim][n_vertices] != Type::Invalid,
+           ExcMessage("The combination of dim = " + std::to_string(dim) +
+                      " and n_vertices = " + std::to_string(n_vertices) +
+                      " does not correspond to a known reference cell type."));
     return table[dim][n_vertices];
   }
 
@@ -347,6 +351,30 @@ namespace ReferenceCell
 
           return 0;
         }
+
+        /**
+         * Map an ExodusII vertex number to a deal.II vertex number.
+         */
+        virtual unsigned int
+        exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const
+        {
+          Assert(false, ExcNotImplemented());
+          (void)vertex_n;
+
+          return 0;
+        }
+
+        /**
+         * Map an ExodusII face number to a deal.II face number.
+         */
+        virtual unsigned int
+        exodusii_face_to_deal_face(const unsigned int face_n) const
+        {
+          Assert(false, ExcNotImplemented());
+          (void)face_n;
+
+          return 0;
+        }
       };
 
 
@@ -415,6 +443,15 @@ namespace ReferenceCell
           (void)face_no;
           return ReferenceCell::Type::Invalid;
         }
+
+        virtual unsigned int
+        exodusii_face_to_deal_face(const unsigned int face_n) const override
+        {
+          (void)face_n;
+          AssertIndexRange(face_n, n_faces());
+
+          return 0;
+        }
       };
 
 
@@ -430,6 +467,21 @@ namespace ReferenceCell
           (void)face_no;
           return ReferenceCell::Type::Vertex;
         }
+
+        virtual unsigned int
+        exodusii_vertex_to_deal_vertex(
+          const unsigned int vertex_n) const override
+        {
+          AssertIndexRange(vertex_n, n_vertices());
+          return vertex_n;
+        }
+
+        virtual unsigned int
+        exodusii_face_to_deal_face(const unsigned int face_n) const override
+        {
+          AssertIndexRange(face_n, n_faces());
+          return face_n;
+        }
       };
 
 
@@ -518,6 +570,21 @@ namespace ReferenceCell
 
           return table[face][face_orientation ? vertex : (1 - vertex)];
         }
+
+        virtual unsigned int
+        exodusii_vertex_to_deal_vertex(
+          const unsigned int vertex_n) const override
+        {
+          AssertIndexRange(vertex_n, n_vertices());
+          return vertex_n;
+        }
+
+        virtual unsigned int
+        exodusii_face_to_deal_face(const unsigned int face_n) const override
+        {
+          AssertIndexRange(face_n, n_faces());
+          return face_n;
+        }
       };
 
 
@@ -553,6 +620,23 @@ namespace ReferenceCell
           (void)face_no;
           return ReferenceCell::Type::Line;
         }
+
+        virtual unsigned int
+        exodusii_vertex_to_deal_vertex(
+          const unsigned int vertex_n) const override
+        {
+          AssertIndexRange(vertex_n, n_vertices());
+          constexpr std::array<unsigned int, 4> exodus_to_deal{{0, 1, 3, 2}};
+          return exodus_to_deal[vertex_n];
+        }
+
+        virtual unsigned int
+        exodusii_face_to_deal_face(const unsigned int face_n) const override
+        {
+          AssertIndexRange(face_n, n_faces());
+          constexpr std::array<unsigned int, 4> exodus_to_deal{{2, 1, 3, 0}};
+          return exodus_to_deal[face_n];
+        }
       };
 
 
@@ -691,6 +775,22 @@ namespace ReferenceCell
           return table[face][standard_to_real_face_vertex(
             vertex, face, face_orientation)];
         }
+
+        virtual unsigned int
+        exodusii_vertex_to_deal_vertex(
+          const unsigned int vertex_n) const override
+        {
+          AssertIndexRange(vertex_n, n_vertices());
+          return vertex_n;
+        }
+
+        virtual unsigned int
+        exodusii_face_to_deal_face(const unsigned int face_n) const override
+        {
+          AssertIndexRange(face_n, n_faces());
+          constexpr std::array<unsigned int, 4> exodus_to_deal{{1, 3, 2, 0}};
+          return exodus_to_deal[face_n];
+        }
       };
 
 
@@ -838,6 +938,23 @@ namespace ReferenceCell
           return table[face][standard_to_real_face_vertex(
             vertex, face, face_orientation)];
         }
+
+        virtual unsigned int
+        exodusii_vertex_to_deal_vertex(
+          const unsigned int vertex_n) const override
+        {
+          AssertIndexRange(vertex_n, n_vertices());
+          constexpr std::array<unsigned int, 5> exodus_to_deal{{0, 1, 3, 2, 4}};
+          return exodus_to_deal[vertex_n];
+        }
+
+        virtual unsigned int
+        exodusii_face_to_deal_face(const unsigned int face_n) const override
+        {
+          AssertIndexRange(face_n, n_faces());
+          constexpr std::array<unsigned int, 5> exodus_to_deal{{3, 2, 4, 1, 0}};
+          return exodus_to_deal[face_n];
+        }
       };
 
 
@@ -985,6 +1102,24 @@ namespace ReferenceCell
           return table[face][standard_to_real_face_vertex(
             vertex, face, face_orientation)];
         }
+
+        virtual unsigned int
+        exodusii_vertex_to_deal_vertex(
+          const unsigned int vertex_n) const override
+        {
+          AssertIndexRange(vertex_n, n_vertices());
+          constexpr std::array<unsigned int, 6> exodus_to_deal{
+            {2, 1, 0, 5, 4, 3}};
+          return exodus_to_deal[vertex_n];
+        }
+
+        virtual unsigned int
+        exodusii_face_to_deal_face(const unsigned int face_n) const override
+        {
+          AssertIndexRange(face_n, n_faces());
+          constexpr std::array<unsigned int, 6> exodus_to_deal{{3, 4, 2, 0, 1}};
+          return exodus_to_deal[face_n];
+        }
       };
 
 
@@ -1074,6 +1209,25 @@ namespace ReferenceCell
           (void)face_no;
           return ReferenceCell::Type::Quad;
         }
+
+        virtual unsigned int
+        exodusii_vertex_to_deal_vertex(
+          const unsigned int vertex_n) const override
+        {
+          AssertIndexRange(vertex_n, n_vertices());
+          constexpr std::array<unsigned int, 8> exodus_to_deal{
+            {0, 1, 3, 2, 4, 5, 7, 6}};
+          return exodus_to_deal[vertex_n];
+        }
+
+        virtual unsigned int
+        exodusii_face_to_deal_face(const unsigned int face_n) const override
+        {
+          AssertIndexRange(face_n, n_faces());
+          constexpr std::array<unsigned int, 6> exodus_to_deal{
+            {2, 1, 3, 0, 4, 5}};
+          return exodus_to_deal[face_n];
+        }
       };
 
       /**
diff --git a/source/base/exceptions.cc b/source/base/exceptions.cc
index c5a23250e4..e3d1e2bb5a 100644
--- a/source/base/exceptions.cc
+++ b/source/base/exceptions.cc
@@ -28,6 +28,10 @@
 #  include <mpi.h>
 #endif
 
+#ifdef DEAL_II_TRILINOS_WITH_SEACAS
+#  include <exodusII.h>
+#endif
+
 #ifdef DEAL_II_HAVE_GLIBC_STACKTRACE
 #  include <execinfo.h>
 #endif
@@ -347,9 +351,9 @@ ExceptionBase::generate_message() const
 
 
 
-#ifdef DEAL_II_WITH_MPI
 namespace StandardExceptions
 {
+#ifdef DEAL_II_WITH_MPI
   ExcMPI::ExcMPI(const int error_code)
     : error_code(error_code)
   {}
@@ -397,9 +401,31 @@ namespace StandardExceptions
     out << "The numerical value of the original error code is " << error_code
         << "." << std::endl;
   }
-} // namespace StandardExceptions
 #endif // DEAL_II_WITH_MPI
 
+
+
+#ifdef DEAL_II_TRILINOS_WITH_SEACAS
+  ExcExodusII::ExcExodusII(const int error_code)
+    : error_code(error_code)
+  {
+    // To avoid including yet another header in exceptions.h we assume that
+    // EX_NOERR is zero. Check that here:
+    static_assert(EX_NOERR == 0,
+                  "EX_NOERR is assumed to be zero in all versions of ExodusII");
+  }
+
+
+
+  void
+  ExcExodusII::print_info(std::ostream &out) const
+  {
+    out << "Error code is " << error_code << '\n';
+    out << "String description: " << ex_strerror(error_code) << '\n';
+  }
+#endif // DEAL_II_TRILINOS_WITH_SEACAS
+} // namespace StandardExceptions
+
 namespace deal_II_exceptions
 {
   namespace internals
diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc
index 84eec6bef3..50d432be0f 100644
--- a/source/grid/grid_in.cc
+++ b/source/grid/grid_in.cc
@@ -41,6 +41,10 @@
 #  include <assimp/scene.h>       // Output data structure
 #endif
 
+#ifdef DEAL_II_TRILINOS_WITH_SEACAS
+#  include <exodusII.h>
+#endif
+
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -2820,6 +2824,364 @@ GridIn<dim, spacedim>::read_assimp(const std::string &filename,
 #endif
 }
 
+#ifdef DEAL_II_TRILINOS_WITH_SEACAS
+// Namespace containing some extra functions for reading ExodusII files
+namespace
+{
+  // Convert ExodusII strings to cell types. Use the number of nodes per element
+  // to disambiguate some cases.
+  ReferenceCell::Type
+  exodusii_name_to_type(const std::string &type_name,
+                        const int          n_nodes_per_element)
+  {
+    Assert(type_name.size() > 0, ExcInternalError());
+    // Try to canonify the name by switching to upper case and removing trailing
+    // numbers. This makes, e.g., pyramid, PYRAMID, PYRAMID5, and PYRAMID13 all
+    // equal.
+    std::string type_name_2 = type_name;
+    std::transform(type_name_2.begin(),
+                   type_name_2.end(),
+                   type_name_2.begin(),
+                   [](unsigned char c) { return std::toupper(c); });
+    const std::string numbers = "0123456789";
+    type_name_2.erase(std::find_first_of(type_name_2.begin(),
+                                         type_name_2.end(),
+                                         numbers.begin(),
+                                         numbers.end()),
+                      type_name_2.end());
+
+    if (type_name_2 == "TRI" || type_name_2 == "TRIANGLE")
+      return ReferenceCell::Type::Tri;
+    else if (type_name_2 == "QUAD" || type_name_2 == "QUADRILATERAL")
+      return ReferenceCell::Type::Quad;
+    else if (type_name_2 == "SHELL")
+      {
+        if (n_nodes_per_element == 3)
+          return ReferenceCell::Type::Tri;
+        else
+          return ReferenceCell::Type::Quad;
+      }
+    else if (type_name_2 == "TET" || type_name_2 == "TETRA" ||
+             type_name_2 == "TETRAHEDRON")
+      return ReferenceCell::Type::Tet;
+    else if (type_name_2 == "PYRA" || type_name_2 == "PYRAMID")
+      return ReferenceCell::Type::Pyramid;
+    else if (type_name_2 == "WEDGE")
+      return ReferenceCell::Type::Wedge;
+    else if (type_name_2 == "HEX" || type_name_2 == "HEXAHEDRON")
+      return ReferenceCell::Type::Hex;
+
+    Assert(false, ExcNotImplemented());
+    return ReferenceCell::Type::Invalid;
+  }
+
+  // Associate deal.II boundary ids with sidesets (a face can be in multiple
+  // sidesets - to translate we assign each set of side set ids to a
+  // boundary_id or manifold_id)
+  template <int dim, int spacedim = dim>
+  std::pair<SubCellData, std::vector<std::vector<int>>>
+  read_exodusii_sidesets(const int                         ex_id,
+                         const int                         n_side_sets,
+                         const std::vector<CellData<dim>> &cells,
+                         const bool apply_all_indicators_to_manifolds)
+  {
+    SubCellData                   subcelldata;
+    std::vector<std::vector<int>> b_or_m_id_to_sideset_ids;
+    // boundary id 0 is the default
+    b_or_m_id_to_sideset_ids.emplace_back();
+    // deal.II does not support assigning boundary ids with nonzero codimension
+    // meshes so completely skip this information in that case.
+    //
+    // Exodus prints warnings if we try to get empty sets so always check first
+    if (dim == spacedim && n_side_sets > 0)
+      {
+        std::vector<int> side_set_ids(n_side_sets);
+        int ierr = ex_get_ids(ex_id, EX_SIDE_SET, side_set_ids.data());
+        AssertThrowExodusII(ierr);
+
+        // First collect all side sets on all boundary faces (indexed here as
+        // max_faces_per_cell * cell_n + face_n). We then sort and uniquify the
+        // side sets so that we can convert a set of side set indices into a
+        // single deal.II boundary or manifold id (and save the correspondence).
+        constexpr auto max_faces_per_cell = GeometryInfo<dim>::faces_per_cell;
+        std::map<std::size_t, std::vector<int>> face_side_sets;
+        for (const int side_set_id : side_set_ids)
+          {
+            int n_sides                = -1;
+            int n_distribution_factors = -1;
+
+            ierr = ex_get_set_param(ex_id,
+                                    EX_SIDE_SET,
+                                    side_set_id,
+                                    &n_sides,
+                                    &n_distribution_factors);
+            AssertThrowExodusII(ierr);
+            if (n_sides > 0)
+              {
+                std::vector<int> elements(n_sides);
+                std::vector<int> faces(n_sides);
+                ierr = ex_get_set(ex_id,
+                                  EX_SIDE_SET,
+                                  side_set_id,
+                                  elements.data(),
+                                  faces.data());
+                AssertThrowExodusII(ierr);
+
+                // According to the manual (subsection 4.8): "The internal
+                // number of an element numbering is defined implicitly by the
+                // order in which it appears in the file. Elements are numbered
+                // internally (beginning with 1) consecutively across all
+                // element blocks." Hence element i in Exodus numbering is entry
+                // i - 1 in the cells array.
+                for (int side_n = 0; side_n < n_sides; ++side_n)
+                  {
+                    const long        element_n = elements[side_n] - 1;
+                    const long        face_n    = faces[side_n] - 1;
+                    const std::size_t face_id =
+                      element_n * max_faces_per_cell + face_n;
+                    face_side_sets[face_id].push_back(side_set_id);
+                  }
+              }
+          }
+
+        // Collect into a sortable data structure:
+        std::vector<std::pair<std::size_t, std::vector<int>>>
+          face_id_to_side_sets;
+        for (auto &pair : face_side_sets)
+          {
+            Assert(pair.second.size() > 0, ExcInternalError());
+            face_id_to_side_sets.push_back(std::move(pair));
+          }
+
+        // sort by side sets:
+        std::sort(face_id_to_side_sets.begin(),
+                  face_id_to_side_sets.end(),
+                  [](const auto &a, const auto &b) {
+                    return std::lexicographical_compare(a.second.begin(),
+                                                        a.second.end(),
+                                                        b.second.begin(),
+                                                        b.second.end());
+                  });
+
+        types::boundary_id current_b_or_m_id = 0;
+        for (const auto &pair : face_id_to_side_sets)
+          {
+            const std::size_t       face_id          = pair.first;
+            const std::vector<int> &face_sideset_ids = pair.second;
+            if (face_sideset_ids != b_or_m_id_to_sideset_ids.back())
+              {
+                // Since we sorted by sideset ids we are guaranteed that if this
+                // doesn't match the last set then it has not yet been seen
+                ++current_b_or_m_id;
+                b_or_m_id_to_sideset_ids.push_back(face_sideset_ids);
+                Assert(current_b_or_m_id == b_or_m_id_to_sideset_ids.size() - 1,
+                       ExcInternalError());
+              }
+            // Record the b_or_m_id of the current face.
+            const unsigned int   local_face_n = face_id % max_faces_per_cell;
+            const CellData<dim> &cell = cells[face_id / max_faces_per_cell];
+            const ReferenceCell::Type cell_type =
+              ReferenceCell::n_vertices_to_type(dim, cell.vertices.size());
+            const ReferenceCell::internal::Info::Base &info =
+              ReferenceCell::internal::Info::get_cell(cell_type);
+            const unsigned int deal_face_n =
+              info.exodusii_face_to_deal_face(local_face_n);
+            const ReferenceCell::internal::Info::Base &face_info =
+              ReferenceCell::internal::Info::get_face(cell_type, deal_face_n);
+
+            // The orientation we pick doesn't matter here since when we create
+            // the Triangulation we will sort the vertices for each CellData
+            // object created here.
+            if (dim == 2)
+              {
+                CellData<1> boundary_line(face_info.n_vertices());
+                if (apply_all_indicators_to_manifolds)
+                  boundary_line.manifold_id = current_b_or_m_id;
+                else
+                  boundary_line.boundary_id = current_b_or_m_id;
+                for (unsigned int j = 0; j < face_info.n_vertices(); ++j)
+                  boundary_line.vertices[j] =
+                    cell
+                      .vertices[info.face_to_cell_vertices(deal_face_n, j, 0)];
+
+                subcelldata.boundary_lines.push_back(std::move(boundary_line));
+              }
+            else if (dim == 3)
+              {
+                CellData<2> boundary_quad(face_info.n_vertices());
+                if (apply_all_indicators_to_manifolds)
+                  boundary_quad.manifold_id = current_b_or_m_id;
+                else
+                  boundary_quad.boundary_id = current_b_or_m_id;
+                for (unsigned int j = 0; j < face_info.n_vertices(); ++j)
+                  boundary_quad.vertices[j] =
+                    cell
+                      .vertices[info.face_to_cell_vertices(deal_face_n, j, 0)];
+
+                subcelldata.boundary_quads.push_back(std::move(boundary_quad));
+              }
+          }
+      }
+
+    return std::make_pair(std::move(subcelldata),
+                          std::move(b_or_m_id_to_sideset_ids));
+  }
+} // namespace
+#endif
+
+template <int dim, int spacedim>
+typename GridIn<dim, spacedim>::ExodusIIData
+GridIn<dim, spacedim>::read_exodusii(
+  const std::string &filename,
+  const bool         apply_all_indicators_to_manifolds)
+{
+#ifdef DEAL_II_TRILINOS_WITH_SEACAS
+  // deal.II always uses double precision numbers for geometry
+  int component_word_size = sizeof(double);
+  // setting to zero uses the stored word size
+  int   floating_point_word_size = 0;
+  float ex_version               = 0.0;
+
+  const int ex_id = ex_open(filename.c_str(),
+                            EX_READ,
+                            &component_word_size,
+                            &floating_point_word_size,
+                            &ex_version);
+  AssertThrow(ex_id > 0,
+              ExcMessage("ExodusII failed to open the specified input file."));
+
+  // Read basic mesh information:
+  std::vector<char> string_temp(MAX_LINE_LENGTH + 1, '\0');
+  int               mesh_dimension   = 0;
+  int               n_nodes          = 0;
+  int               n_elements       = 0;
+  int               n_element_blocks = 0;
+  int               n_node_sets      = 0;
+  int               n_side_sets      = 0;
+
+  int ierr = ex_get_init(ex_id,
+                         string_temp.data(),
+                         &mesh_dimension,
+                         &n_nodes,
+                         &n_elements,
+                         &n_element_blocks,
+                         &n_node_sets,
+                         &n_side_sets);
+  AssertThrowExodusII(ierr);
+  AssertDimension(mesh_dimension, spacedim);
+
+  // Read nodes:
+  std::vector<double> xs(n_nodes);
+  std::vector<double> ys(n_nodes);
+  std::vector<double> zs(n_nodes);
+
+  ierr = ex_get_coord(ex_id, xs.data(), ys.data(), zs.data());
+  AssertThrowExodusII(ierr);
+
+  // Even if there is a node numbering array the values stored inside the
+  // ExodusII file must use the contiguous, internal ordering (see Section 4.5
+  // of the manual - "Internal (contiguously numbered) node and element IDs must
+  // be used for all data structures that contain node or element numbers (IDs),
+  // including node set node lists, side set element lists, and element
+  // connectivity.")
+  std::vector<Point<spacedim>> vertices;
+  vertices.reserve(n_nodes);
+  for (int vertex_n = 0; vertex_n < n_nodes; ++vertex_n)
+    {
+      switch (spacedim)
+        {
+          case 1:
+            vertices.emplace_back(xs[vertex_n]);
+            break;
+          case 2:
+            vertices.emplace_back(xs[vertex_n], ys[vertex_n]);
+            break;
+          case 3:
+            vertices.emplace_back(xs[vertex_n], ys[vertex_n], zs[vertex_n]);
+            break;
+          default:
+            Assert(spacedim <= 3, ExcNotImplemented());
+        }
+    }
+
+  std::vector<int> element_block_ids(n_element_blocks);
+  ierr = ex_get_ids(ex_id, EX_ELEM_BLOCK, element_block_ids.data());
+  AssertThrowExodusII(ierr);
+
+  std::vector<CellData<dim>> cells;
+  // Elements are grouped together by same reference cell type in element
+  // blocks. There may be multiple blocks for a single reference cell type,
+  // but "each element block may contain only one element type".
+  for (const int element_block_id : element_block_ids)
+    {
+      std::fill(string_temp.begin(), string_temp.end(), '\0');
+      int n_block_elements         = 0;
+      int n_nodes_per_element      = 0;
+      int n_edges_per_element      = 0;
+      int n_faces_per_element      = 0;
+      int n_attributes_per_element = 0;
+
+      // Extract element data.
+      ierr = ex_get_block(ex_id,
+                          EX_ELEM_BLOCK,
+                          element_block_id,
+                          string_temp.data(),
+                          &n_block_elements,
+                          &n_nodes_per_element,
+                          &n_edges_per_element,
+                          &n_faces_per_element,
+                          &n_attributes_per_element);
+      AssertThrowExodusII(ierr);
+      const ReferenceCell::Type type =
+        exodusii_name_to_type(string_temp.data(), n_nodes_per_element);
+      const ReferenceCell::internal::Info::Base &info =
+        ReferenceCell::internal::Info::get_cell(type);
+      // The number of nodes per element may be larger than what we want to
+      // read - for example, if the Exodus file contains a QUAD9 element, we
+      // only want to read the first four values and ignore the rest.
+      Assert(int(info.n_vertices()) <= n_nodes_per_element, ExcInternalError());
+
+      std::vector<int> connection(n_nodes_per_element * n_block_elements);
+      ierr = ex_get_conn(ex_id,
+                         EX_ELEM_BLOCK,
+                         element_block_id,
+                         connection.data(),
+                         nullptr,
+                         nullptr);
+      AssertThrowExodusII(ierr);
+
+      for (unsigned int elem_n = 0; elem_n < connection.size();
+           elem_n += n_nodes_per_element)
+        {
+          CellData<dim> cell(info.n_vertices());
+          for (unsigned int i : info.vertex_indices())
+            {
+              cell.vertices[info.exodusii_vertex_to_deal_vertex(i)] =
+                connection[elem_n + i] - 1;
+            }
+          cell.material_id = element_block_id;
+          cells.push_back(cell);
+        }
+    }
+
+  // Extract boundary data.
+  auto pair = read_exodusii_sidesets<dim, spacedim>(
+    ex_id, n_side_sets, cells, apply_all_indicators_to_manifolds);
+  ierr = ex_close(ex_id);
+  AssertThrowExodusII(ierr);
+
+  tria->create_triangulation(vertices, cells, pair.first);
+  ExodusIIData out;
+  out.id_to_sideset_ids = std::move(pair.second);
+  return out;
+#else
+  (void)filename;
+  (void)apply_all_indicators_to_manifolds;
+  AssertThrow(false, ExcNeedsExodusII());
+  return {};
+#endif
+}
+
 
 template <int dim, int spacedim>
 void
@@ -3053,6 +3415,10 @@ GridIn<dim, spacedim>::read(const std::string &filename, Format format)
     {
       read_assimp(name);
     }
+  else if (format == exodusii)
+    {
+      read_exodusii(name);
+    }
   else
     {
       std::ifstream in(name.c_str());
@@ -3113,6 +3479,13 @@ GridIn<dim, spacedim>::read(std::istream &in, Format format)
                           "functions, instead."));
         return;
 
+      case exodusii:
+        Assert(false,
+               ExcMessage("There is no read_exodusii(istream &) function. "
+                          "Use the read_exodusii(string &filename, ...) "
+                          "function, instead."));
+        return;
+
       case Default:
         break;
     }
@@ -3129,6 +3502,8 @@ GridIn<dim, spacedim>::default_suffix(const Format format)
     {
       case dbmesh:
         return ".dbmesh";
+      case exodusii:
+        return ".e";
       case msh:
         return ".msh";
       case vtk:
@@ -3161,6 +3536,9 @@ GridIn<dim, spacedim>::parse_format(const std::string &format_name)
   if (format_name == "dbmesh")
     return dbmesh;
 
+  if (format_name == "exodusii")
+    return exodusii;
+
   if (format_name == "msh")
     return msh;
 
@@ -3212,7 +3590,7 @@ template <int dim, int spacedim>
 std::string
 GridIn<dim, spacedim>::get_format_names()
 {
-  return "dbmesh|msh|unv|vtk|vtu|ucd|abaqus|xda|tecplot|assimp";
+  return "dbmesh|exodusii|msh|unv|vtk|vtu|ucd|abaqus|xda|tecplot|assimp";
 }
 
 
diff --git a/tests/grid/grid_in_exodusii.cc b/tests/grid/grid_in_exodusii.cc
new file mode 100644
index 0000000000..b260e2bb3b
--- /dev/null
+++ b/tests/grid/grid_in_exodusii.cc
@@ -0,0 +1,106 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2020 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.
+//
+// ---------------------------------------------------------------------
+
+// Read in several Exodus meshes and save the results to output files.
+
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <string>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim, int spacedim = dim>
+void
+read_and_print(const std::string &filename,
+               const bool         read_as_manifolds = false)
+{
+  Triangulation<dim, spacedim> tria;
+  GridIn<dim, spacedim>        gi;
+  gi.attach_triangulation(tria);
+  const typename GridIn<dim, spacedim>::ExodusIIData exodus_data =
+    gi.read_exodusii(filename, read_as_manifolds);
+  const auto &boundary_id_to_sideset_ids = exodus_data.id_to_sideset_ids;
+
+  for (unsigned int i = 0; i < boundary_id_to_sideset_ids.size(); ++i)
+    {
+      deallog << "boundary id = " << i << " sideset ids = ";
+      for (const int v : boundary_id_to_sideset_ids[i])
+        deallog << v << ", ";
+      deallog << std::endl;
+    }
+
+  for (const auto face : tria.active_face_iterators())
+    {
+      if (face->at_boundary())
+        {
+          deallog << "face center = " << face->center();
+          if (read_as_manifolds)
+            deallog << "   face manifold id = " << face->manifold_id();
+          else
+            deallog << "   face boundary id = " << face->boundary_id();
+          deallog << std::endl;
+        }
+    }
+
+  deallog << "Number of vertices: " << tria.get_vertices().size() << std::endl;
+  deallog << "Number of cells: " << tria.n_cells() << std::endl;
+
+  GridOut go;
+  go.write_vtk(tria, deallog.get_file_stream());
+}
+
+int
+main()
+{
+  initlog();
+  deallog.get_file_stream() << std::setprecision(2);
+
+  deallog << "-----------" << std::endl;
+  deallog << "five-square" << std::endl;
+  deallog << "-----------" << std::endl;
+  read_and_print<2>(SOURCE_DIR "/grids/exodusii/five.e");
+
+  deallog << "-----------------" << std::endl;
+  deallog << "five-square in 3D" << std::endl;
+  deallog << "-----------------" << std::endl;
+  read_and_print<2, 3>(SOURCE_DIR "/grids/exodusii/five-3d.e");
+
+  deallog << "-----------------------------" << std::endl;
+  deallog << "Straightforward 3D mixed mesh" << std::endl;
+  deallog << "-----------------------------" << std::endl;
+  read_and_print<3>(SOURCE_DIR "/grids/exodusii/four.e");
+
+  deallog << "--------------------------" << std::endl;
+  deallog << "Higher-order 3D mixed mesh" << std::endl;
+  deallog << "--------------------------" << std::endl;
+  read_and_print<3>(SOURCE_DIR "/grids/exodusii/four-p.e");
+
+  deallog << "--------------------------------------" << std::endl;
+  deallog << "3D mixed mesh with boundary indicators" << std::endl;
+  deallog << "--------------------------------------" << std::endl;
+  read_and_print<3>(SOURCE_DIR "/grids/exodusii/four-sidesets.e");
+
+  deallog << "--------------------------------------" << std::endl;
+  deallog << "3D mixed mesh with manifold indicators" << std::endl;
+  deallog << "--------------------------------------" << std::endl;
+  read_and_print<3>(SOURCE_DIR "/grids/exodusii/four-sidesets.e", true);
+
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/grid/grid_in_exodusii.with_trilinos_with_seacas=on.output b/tests/grid/grid_in_exodusii.with_trilinos_with_seacas=on.output
new file mode 100644
index 0000000000..20404599dd
--- /dev/null
+++ b/tests/grid/grid_in_exodusii.with_trilinos_with_seacas=on.output
@@ -0,0 +1,408 @@
+
+DEAL::-----------
+DEAL::five-square
+DEAL::-----------
+DEAL::boundary id = 0 sideset ids = 
+DEAL::boundary id = 1 sideset ids = 0, 
+DEAL::boundary id = 2 sideset ids = 0, 1, 
+DEAL::boundary id = 3 sideset ids = 0, 1, 10, 
+DEAL::face center = 3.50000 0.00000   face boundary id = 2
+DEAL::face center = 0.00000 3.50000   face boundary id = 1
+DEAL::face center = 7.00000 3.50000   face boundary id = 3
+DEAL::face center = 3.50000 7.00000   face boundary id = 0
+DEAL::Number of vertices: 8
+DEAL::Number of cells: 5
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 8 double
+0.0 0.0 0
+7.0 0.0 0
+5.0 2.0 0
+2.0 2.0 0
+2.0 5.0 0
+0.0 7.0 0
+5.0 5.0 0
+7.0 7.0 0
+
+CELLS 8 34
+4 0 1 2 3
+4 2 1 7 6
+4 4 6 7 5
+4 0 3 4 5
+4 3 2 6 4
+2 0 1
+2 0 5
+2 1 7
+
+CELL_TYPES 8
+9 9 9 9 9 
+3 3 3 
+
+
+CELL_DATA 8
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+1 1 1 1 10 
+2 1 3 
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 -1 
+-1 -1 -1 
+
+DEAL::-----------------
+DEAL::five-square in 3D
+DEAL::-----------------
+DEAL::boundary id = 0 sideset ids = 
+DEAL::face center = 3.50000 0.00000 3.50000   face boundary id = 0
+DEAL::face center = 0.00000 3.50000 0.00000   face boundary id = 0
+DEAL::face center = 7.00000 3.50000 7.00000   face boundary id = 0
+DEAL::face center = 3.50000 7.00000 3.50000   face boundary id = 0
+DEAL::Number of vertices: 8
+DEAL::Number of cells: 5
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 8 double
+0.0 0.0 0.0
+7.0 0.0 7.0
+5.0 2.0 5.0
+2.0 2.0 2.0
+2.0 5.0 2.0
+0.0 7.0 0.0
+5.0 5.0 5.0
+7.0 7.0 7.0
+
+CELLS 5 25
+4 0 1 2 3
+4 0 3 4 5
+4 1 7 6 2
+4 5 4 6 7
+4 3 2 6 4
+
+CELL_TYPES 5
+9 9 9 9 9 
+
+
+
+CELL_DATA 5
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+1 1 1 1 10 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 -1 
+
+
+DEAL::-----------------------------
+DEAL::Straightforward 3D mixed mesh
+DEAL::-----------------------------
+DEAL::boundary id = 0 sideset ids = 
+DEAL::face center = 1.50000 0.00000 1.50000   face boundary id = 0
+DEAL::face center = 1.50000 1.50000 0.00000   face boundary id = 0
+DEAL::face center = 0.00000 1.50000 1.50000   face boundary id = 0
+DEAL::face center = 3.00000 1.50000 1.50000   face boundary id = 0
+DEAL::face center = 1.50000 3.00000 1.50000   face boundary id = 0
+DEAL::face center = 1.33333 0.333333 3.66667   face boundary id = 0
+DEAL::face center = 0.333333 1.33333 3.66667   face boundary id = 0
+DEAL::face center = 2.33333 0.666667 4.33333   face boundary id = 0
+DEAL::face center = 5.50000 0.500000 4.00000   face boundary id = 0
+DEAL::face center = 5.50000 1.50000 3.00000   face boundary id = 0
+DEAL::face center = 1.33333 2.33333 3.66667   face boundary id = 0
+DEAL::face center = 2.33333 1.66667 4.33333   face boundary id = 0
+DEAL::face center = 5.50000 2.00000 4.00000   face boundary id = 0
+DEAL::face center = 8.00000 1.33333 3.66667   face boundary id = 0
+DEAL::Number of vertices: 13
+DEAL::Number of cells: 4
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 13 double
+0.0 0.0 0.0
+3.0 0.0 0.0
+3.0 3.0 0.0
+0.0 3.0 0.0
+0.0 0.0 3.0
+3.0 0.0 3.0
+3.0 3.0 3.0
+0.0 3.0 3.0
+1.0 1.0 5.0
+3.0 1.0 5.0
+8.0 0.0 3.0
+8.0 3.0 3.0
+8.0 1.0 5.0
+
+CELLS 4 27
+8 0 1 5 4 3 2 6 7
+5 4 5 7 6 8
+4 5 9 6 8
+6 6 5 9 11 10 12
+
+CELL_TYPES 4
+12 14 10 13 
+
+
+
+CELL_DATA 4
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+1 10 100 1000 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 
+
+
+DEAL::--------------------------
+DEAL::Higher-order 3D mixed mesh
+DEAL::--------------------------
+DEAL::boundary id = 0 sideset ids = 
+DEAL::face center = 1.50000 0.00000 1.50000   face boundary id = 0
+DEAL::face center = 1.50000 1.50000 0.00000   face boundary id = 0
+DEAL::face center = 0.00000 1.50000 1.50000   face boundary id = 0
+DEAL::face center = 3.00000 1.50000 1.50000   face boundary id = 0
+DEAL::face center = 1.50000 3.00000 1.50000   face boundary id = 0
+DEAL::face center = 1.33333 0.333333 3.66667   face boundary id = 0
+DEAL::face center = 0.333333 1.33333 3.66667   face boundary id = 0
+DEAL::face center = 2.33333 0.666667 4.33333   face boundary id = 0
+DEAL::face center = 5.50000 0.500000 4.00000   face boundary id = 0
+DEAL::face center = 5.50000 1.50000 3.00000   face boundary id = 0
+DEAL::face center = 1.33333 2.33333 3.66667   face boundary id = 0
+DEAL::face center = 2.33333 1.66667 4.33333   face boundary id = 0
+DEAL::face center = 5.50000 2.00000 4.00000   face boundary id = 0
+DEAL::face center = 8.00000 1.33333 3.66667   face boundary id = 0
+DEAL::Number of vertices: 29
+DEAL::Number of cells: 4
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 29 double
+0.0 0.0 0.0
+3.0 0.0 0.0
+3.0 3.0 0.0
+0.0 3.0 0.0
+0.0 0.0 3.0
+3.0 0.0 3.0
+3.0 3.0 3.0
+0.0 3.0 3.0
+1.0 1.0 5.0
+3.0 1.0 5.0
+8.0 0.0 3.0
+8.0 3.0 3.0
+8.0 1.0 5.0
+3.0 1.5 0.0
+1.5 3.0 0.0
+0.0 1.5 0.0
+1.5 0.0 0.0
+3.0 0.0 1.5
+3.0 3.0 1.5
+0.0 3.0 1.5
+0.0 0.0 1.5
+3.0 1.5 3.0
+1.5 3.0 3.0
+0.0 1.5 3.0
+1.5 0.0 3.0
+0.50 0.50 4.0
+2.0 0.50 4.0
+2.0 2.0 4.0
+0.50 2.0 4.0
+
+CELLS 4 27
+8 0 1 5 4 3 2 6 7
+5 4 5 7 6 8
+4 5 9 6 8
+6 6 5 9 11 10 12
+
+CELL_TYPES 4
+12 14 10 13 
+
+
+
+CELL_DATA 4
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+1 10 100 1000 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 
+
+
+DEAL::--------------------------------------
+DEAL::3D mixed mesh with boundary indicators
+DEAL::--------------------------------------
+DEAL::boundary id = 0 sideset ids = 
+DEAL::boundary id = 1 sideset ids = 1, 
+DEAL::boundary id = 2 sideset ids = 1, 5, 
+DEAL::boundary id = 3 sideset ids = 2, 
+DEAL::boundary id = 4 sideset ids = 3, 
+DEAL::boundary id = 5 sideset ids = 4, 5, 
+DEAL::boundary id = 6 sideset ids = 5, 
+DEAL::face center = 1.50000 0.00000 1.50000   face boundary id = 1
+DEAL::face center = 1.50000 1.50000 0.00000   face boundary id = 2
+DEAL::face center = 0.00000 1.50000 1.50000   face boundary id = 5
+DEAL::face center = 3.00000 1.50000 1.50000   face boundary id = 3
+DEAL::face center = 1.50000 3.00000 1.50000   face boundary id = 4
+DEAL::face center = 1.33333 0.333333 3.66667   face boundary id = 1
+DEAL::face center = 0.333333 1.33333 3.66667   face boundary id = 5
+DEAL::face center = 2.33333 0.666667 4.33333   face boundary id = 1
+DEAL::face center = 5.50000 0.500000 4.00000   face boundary id = 1
+DEAL::face center = 5.50000 1.50000 3.00000   face boundary id = 4
+DEAL::face center = 1.33333 2.33333 3.66667   face boundary id = 4
+DEAL::face center = 2.33333 1.66667 4.33333   face boundary id = 3
+DEAL::face center = 5.50000 2.00000 4.00000   face boundary id = 3
+DEAL::face center = 8.00000 1.33333 3.66667   face boundary id = 6
+DEAL::Number of vertices: 13
+DEAL::Number of cells: 4
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 13 double
+0.0 0.0 0.0
+3.0 0.0 0.0
+3.0 3.0 0.0
+0.0 3.0 0.0
+0.0 0.0 3.0
+3.0 0.0 3.0
+3.0 3.0 3.0
+0.0 3.0 3.0
+1.0 1.0 5.0
+3.0 1.0 5.0
+8.0 0.0 3.0
+8.0 3.0 3.0
+8.0 1.0 5.0
+
+CELLS 18 91
+8 0 1 5 4 3 2 6 7
+5 4 5 7 6 8
+4 5 9 6 8
+6 6 5 9 11 10 12
+4 0 4 5 1
+4 0 1 2 3
+4 0 3 7 4
+4 1 2 6 5
+4 3 7 6 2
+3 5 4 8
+3 4 7 8
+3 9 5 8
+4 5 9 12 10
+4 6 5 10 11
+3 7 6 8
+3 6 9 8
+4 9 6 11 12
+3 11 10 12
+
+CELL_TYPES 18
+12 14 10 13 
+9 9 9 9 9 5 5 5 9 9 5 5 9 5 
+
+
+CELL_DATA 18
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+1 10 100 1000 
+1 2 5 3 4 1 5 1 1 4 4 3 3 6 
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 
+-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 
+
+DEAL::--------------------------------------
+DEAL::3D mixed mesh with manifold indicators
+DEAL::--------------------------------------
+DEAL::boundary id = 0 sideset ids = 
+DEAL::boundary id = 1 sideset ids = 1, 
+DEAL::boundary id = 2 sideset ids = 1, 5, 
+DEAL::boundary id = 3 sideset ids = 2, 
+DEAL::boundary id = 4 sideset ids = 3, 
+DEAL::boundary id = 5 sideset ids = 4, 5, 
+DEAL::boundary id = 6 sideset ids = 5, 
+DEAL::face center = 1.50000 0.00000 1.50000   face manifold id = 1
+DEAL::face center = 1.50000 1.50000 0.00000   face manifold id = 2
+DEAL::face center = 0.00000 1.50000 1.50000   face manifold id = 5
+DEAL::face center = 3.00000 1.50000 1.50000   face manifold id = 3
+DEAL::face center = 1.50000 3.00000 1.50000   face manifold id = 4
+DEAL::face center = 1.33333 0.333333 3.66667   face manifold id = 1
+DEAL::face center = 0.333333 1.33333 3.66667   face manifold id = 5
+DEAL::face center = 2.33333 0.666667 4.33333   face manifold id = 1
+DEAL::face center = 5.50000 0.500000 4.00000   face manifold id = 1
+DEAL::face center = 5.50000 1.50000 3.00000   face manifold id = 4
+DEAL::face center = 1.33333 2.33333 3.66667   face manifold id = 4
+DEAL::face center = 2.33333 1.66667 4.33333   face manifold id = 3
+DEAL::face center = 5.50000 2.00000 4.00000   face manifold id = 3
+DEAL::face center = 8.00000 1.33333 3.66667   face manifold id = 6
+DEAL::Number of vertices: 13
+DEAL::Number of cells: 4
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 13 double
+0.0 0.0 0.0
+3.0 0.0 0.0
+3.0 3.0 0.0
+0.0 3.0 0.0
+0.0 0.0 3.0
+3.0 0.0 3.0
+3.0 3.0 3.0
+0.0 3.0 3.0
+1.0 1.0 5.0
+3.0 1.0 5.0
+8.0 0.0 3.0
+8.0 3.0 3.0
+8.0 1.0 5.0
+
+CELLS 18 91
+8 0 1 5 4 3 2 6 7
+5 4 5 7 6 8
+4 5 9 6 8
+6 6 5 9 11 10 12
+4 0 4 5 1
+4 0 1 2 3
+4 0 3 7 4
+4 1 2 6 5
+4 3 7 6 2
+3 5 4 8
+3 4 7 8
+3 9 5 8
+4 5 9 12 10
+4 6 5 10 11
+3 7 6 8
+3 6 9 8
+4 9 6 11 12
+3 11 10 12
+
+CELL_TYPES 18
+12 14 10 13 
+9 9 9 9 9 5 5 5 9 9 5 5 9 5 
+
+
+CELL_DATA 18
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+1 10 100 1000 
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 
+1 2 5 3 4 1 5 1 1 4 4 3 3 6 
+
+DEAL::OK
diff --git a/tests/grid/grids/exodusii/five-3d.e b/tests/grid/grids/exodusii/five-3d.e
new file mode 100644
index 0000000000..a6b3b21726
Binary files /dev/null and b/tests/grid/grids/exodusii/five-3d.e differ
diff --git a/tests/grid/grids/exodusii/five.e b/tests/grid/grids/exodusii/five.e
new file mode 100644
index 0000000000..ea73958468
Binary files /dev/null and b/tests/grid/grids/exodusii/five.e differ
diff --git a/tests/grid/grids/exodusii/four-p.e b/tests/grid/grids/exodusii/four-p.e
new file mode 100644
index 0000000000..2c90ba402c
Binary files /dev/null and b/tests/grid/grids/exodusii/four-p.e differ
diff --git a/tests/grid/grids/exodusii/four-sidesets.e b/tests/grid/grids/exodusii/four-sidesets.e
new file mode 100644
index 0000000000..a259d50f6c
Binary files /dev/null and b/tests/grid/grids/exodusii/four-sidesets.e differ
diff --git a/tests/grid/grids/exodusii/four.e b/tests/grid/grids/exodusii/four.e
new file mode 100644
index 0000000000..0a33d08f04
Binary files /dev/null and b/tests/grid/grids/exodusii/four.e differ