From 36726bce8f8f380a6a2110629a266e93f46f4135 Mon Sep 17 00:00:00 2001 From: David Wells Date: Wed, 9 Sep 2020 13:03:14 -0400 Subject: [PATCH] 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 #include #include #include // 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 #include #include #include 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 #include #include #include 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 #include #include #include // 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 #include #include #include // 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); } --- doc/news/changes/major/20200909DavidWells | 4 + include/deal.II/base/exceptions.h | 65 +++ include/deal.II/base/undefine_macros.h | 4 + include/deal.II/grid/grid_in.h | 79 ++++ include/deal.II/grid/reference_cell.h | 154 +++++++ source/base/exceptions.cc | 30 +- source/grid/grid_in.cc | 380 +++++++++++++++- tests/grid/grid_in_exodusii.cc | 106 +++++ ...odusii.with_trilinos_with_seacas=on.output | 408 ++++++++++++++++++ tests/grid/grids/exodusii/five-3d.e | Bin 0 -> 4000 bytes tests/grid/grids/exodusii/five.e | Bin 0 -> 3636 bytes tests/grid/grids/exodusii/four-p.e | Bin 0 -> 4004 bytes tests/grid/grids/exodusii/four-sidesets.e | Bin 0 -> 5740 bytes tests/grid/grids/exodusii/four.e | Bin 0 -> 3532 bytes 14 files changed, 1227 insertions(+), 3 deletions(-) create mode 100644 doc/news/changes/major/20200909DavidWells create mode 100644 tests/grid/grid_in_exodusii.cc create mode 100644 tests/grid/grid_in_exodusii.with_trilinos_with_seacas=on.output create mode 100644 tests/grid/grids/exodusii/five-3d.e create mode 100644 tests/grid/grids/exodusii/five.e create mode 100644 tests/grid/grids/exodusii/four-p.e create mode 100644 tests/grid/grids/exodusii/four-sidesets.e create mode 100644 tests/grid/grids/exodusii/four.e 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. +
+(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 + * deal.II 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 \#undef + * this macro, as well as all other macros defined by deal.II that are not + * prefixed with either DEAL or deal, by including + * the header deal.II/base/undefine_macros.h 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> 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: + * + *
    + *
  1. Block ids: the block id of an element is loaded as its material + * id.
  2. + *
  3. 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 QUAD9 element in deal.II since all quadrilaterals have + * four vertices and additional geometric information is either stored in + * a Manifold or something like MappingQEulerian).
  4. + *
  5. 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.
  6. + *
+ * + * 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 false, 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 true 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 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 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 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 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 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 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 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 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 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 #endif +#ifdef DEAL_II_TRILINOS_WITH_SEACAS +# include +#endif + #ifdef DEAL_II_HAVE_GLIBC_STACKTRACE # include #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 // Output data structure #endif +#ifdef DEAL_II_TRILINOS_WITH_SEACAS +# include +#endif + DEAL_II_NAMESPACE_OPEN @@ -2820,6 +2824,364 @@ GridIn::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 + std::pair>> + read_exodusii_sidesets(const int ex_id, + const int n_side_sets, + const std::vector> &cells, + const bool apply_all_indicators_to_manifolds) + { + SubCellData subcelldata; + std::vector> 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 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::faces_per_cell; + std::map> 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 elements(n_sides); + std::vector 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>> + 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 &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 &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 +typename GridIn::ExodusIIData +GridIn::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 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 xs(n_nodes); + std::vector ys(n_nodes); + std::vector 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> 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 element_block_ids(n_element_blocks); + ierr = ex_get_ids(ex_id, EX_ELEM_BLOCK, element_block_ids.data()); + AssertThrowExodusII(ierr); + + std::vector> 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 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 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( + 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 void @@ -3053,6 +3415,10 @@ GridIn::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::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::default_suffix(const Format format) { case dbmesh: return ".dbmesh"; + case exodusii: + return ".e"; case msh: return ".msh"; case vtk: @@ -3161,6 +3536,9 @@ GridIn::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 std::string GridIn::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 +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +read_and_print(const std::string &filename, + const bool read_as_manifolds = false) +{ + Triangulation tria; + GridIn gi; + gi.attach_triangulation(tria); + const typename GridIn::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 0000000000000000000000000000000000000000..a6b3b21726728bc4e5b057b146ac15ce8d90a3c5 GIT binary patch literal 4000 zcmeH}OK#La5QZltnb?y_1}Q99u*em72(h$86u};`b)1ZoSop=*LmtWjSh8fnA}8QX z?6?O1*X^1NLD3M26`GQ|tE>9=U)9y^R=?icf6((XklI!Y^0bW0q)2a1 zmB61m{F?d^e#MFFbRupyf&X{>GlwVnTX;7Ces~kNfrOazK7iym?|UvA&Y{|d+GbK` z@uzgsWK}8Ct#Kv#`TY^|claM4=T*{XbB%nCoEHW%yDRFyI`(^rCe4Hj@VG#B{VxY6R zt1jsH?4J(hS<#%g2wM8deUibM@i)R8{L_` z$mTn|`AtXNV5V>3@B02<#kV|`dRc*G1(p?9R$y6yWd;803h<8M9m=DZw-Pp*_c!lm f=6IX&P6dOHIeWK)#gC2O_{@*y8}nNnc-MXbfd<$+ literal 0 HcmV?d00001 diff --git a/tests/grid/grids/exodusii/five.e b/tests/grid/grids/exodusii/five.e new file mode 100644 index 0000000000000000000000000000000000000000..ea7395846878d3e9ac01f0ffeb5b8e9f501af275 GIT binary patch literal 3636 zcmeH}OKuZE5Qd$Hof+FP2(nCIDlc;Wz3eC$w zvB;xYnxnBlL@&x&Qv3BlLx3&G5mq1U6c?2`kr?+{y^&cvTI?cP1vCoMKDfWLP5 zef1Bp%TH7#b8))`?7!pRI6T3x;oSoG@vYw)5`1=T0Ew?_dn+BzvD$&!D6Z1zOERz1 zvXJV|xHJ9n?TGa|>`yMTGH%l1GOEh7Xrk+KKGXTXD%U;Yf}h-l7ip%P_=V0I>yPvJ zCe6u%Y80thm(7*TtNqV3AB+6mv5FQq3({)bn)dol+GIGCaK<2Ct}G+zkJN^)=GRwc z3&bzYWS|LgN+!qZ!sP7MnRKw!+ywtEen9>L;j4LBk$b~pXTI;hJr~>FhZ9+R{I;L` zhIRdiKVmY$AAFL*Q+0vxU4LLV`Lq2mAUQ8psT<9FjC?O9a8123AQv9HbQaEY7q-zF zdEmIzh+0_AEVf~Enr7L@I9vFFfq@gz92FUSSH~H@3Q}S5zmA{sOOZ4u8c09{Z2zIaDV*^-$U5M z_fuWq{MYU-Aiw+$Ogi7vgn8iatBqLrjh8s#B|a0!<&}y5M%_JFZ^Ay;9^Ci|u z!b_7(UaGTrH(_51ON`dYr_C7=M(j?eL@k*fwWA-bP_5QkZ}F>zHT=|c<6GyyWB)qe zR^Mh|n}KZxwi(!F;6Kd(|HHhac(n0f&xp>OinkT-NPA=P4g?PlJG$*Ltnp$0+uUsU G*!=|@j?rHL literal 0 HcmV?d00001 diff --git a/tests/grid/grids/exodusii/four-p.e b/tests/grid/grids/exodusii/four-p.e new file mode 100644 index 0000000000000000000000000000000000000000..2c90ba402ce0c73add59473876ed57b230a7613c GIT binary patch literal 4004 zcmeHKOK;Oa5O!(aDK$_Y;U)FRk(P>UEKyZ(Pq>=MWv*RICVUfO}RCPvegUMBA>j>6GRJW9Op2ix=iE94&Zvz-@3p%60fjO=mhwIjITW3ci11&g~le;4IFf zLZ_23&vaUZm)Ue2&h!=fZu-SObimy+9i#u&vE-wFFNx+lVeheDAvSJOoWoqs-^X|d z!CxO?1Eb=C+rv7|#zHm=UBp}v#BGOA489QV4fiVD;`~S zdb1dGV%ZNNu7AfW&3(uY@cA^$xygd%x*Rv`IW)O#zdMG-<}dSKLT`;SbZ%~ZoWt`d z#N2*Zr5Qk4?7~>qU|TOr?A70}ZRWxnm@50rK1#h(M=w0lar_~QFKn}k;rb+u*%DIw z7al#VzdtDp4_z##t%9$df5%x{dN=FXPq1?&Pyq}6VoVizf% z;ZHjfd$8~Hs&w&(+<_$$y`1lGOAs>@zu?6$`taIX0_PAZc4sby#1mOM+(7IirB0E8 zU-7&8N_(Z=>O8@TUFs3L=%xSKj-;-@IYc_OR!H1h$uB-A5Lr1FVz;RJr{rBHYuA_c z+#3CEo`Ro&`5&3<@h-{xIPn}y4!pCPJMeGWbL;Z} literal 0 HcmV?d00001 diff --git a/tests/grid/grids/exodusii/four-sidesets.e b/tests/grid/grids/exodusii/four-sidesets.e new file mode 100644 index 0000000000000000000000000000000000000000..a259d50f6c6d42826d0b9691697cdb3a453f122a GIT binary patch literal 5740 zcmeHLJ#Q015Os)?+&N#6U@LSkC@A>|R7hCDiiAW#MTK;_Lg{f+_mT2*H zoT{MIs^HAZT6UHOQJiwNjOKHm<&i3_eG{EWok|ozZM2S^2ZzZ=1EOw|x-yQyue5u1 z1^O+AjWfcZ7(9d7rXGHRLXE&}4Bia>={&qK%d3a-3g>_@Pu9890 zXH~~Qn768q&gwc6vOZPZbIcjj2HG0hdRWB4Cp9YLJj2xZS{nTN%37|?(eE86d8p&; zBq;JY)4_Q@icp73)ZB10*UYPb94DxaHB@Z0@225JoU$I+D*-mmi9UsT@_ZlT4UoUs z$L0z38OzJK&89-XqGO%VA%tu6((y#Ij`;($j@f(Xr?~_2S$Q7A7IV!wA-x1YzB-&b z!<<;=UWESd+KI>KJm8B_Ua&@jWnMXM+cRKdTYfWy#>Owt?_utj<%~bZBtn0z1H^m8 zAK1CV%q$OW;1< z!+NrAl3&SNCpd_cs!OLBKhE^)m{t`2iEVY9K(H|3Um zlA;|G+fkln>PT~MkcVf@t8n>zw*`7!*ku|_`&60YaFKm9NVm$xb5AZGjqh8D>pCv4 z7>~j|bIXFLOv!h4&%G!n+9jm{=Dt^GeW>|Xe$aRG zT5q+}-wEFNom%XuZWlk&3tn(_z1Zu$;OF*kw-2PRnC$1|VpWZquiJnDFgz|Jv(L}x0CQ{CyD&=OjcYPGbwpyX+o ztx`2JDWxt(NIA=Adg<*Wuz)(9Xgc*UA}_1kHp#Bhe- z!2d;)zccg?{9iWtJHufh1b5=4S~0 z>Iw%K+ZB(8b-IHE?^|ONZV2K%2-3|~a|)UDMeh99i@b8=T3L^2hd3}#(C9lpp}yU1 z2ARIZLx}rdc%gp@9>BjS^MZ#gzTB7dp1+4Kx9+=ZXng*%{u$)KGzW6`&c-#ok3!7N z0fH1>WA1%7&V0=G^Vav_Z@vT8&Kj61@g$S1gm+eAg?yvRZi>Bo&|7E+C=(lFfJXbV_`n#8 zwTBqTVg<#fj8ij%n=6!W0yO%BkL1P&$awq>l6yyd{+OH97A6<#f)m90f{(=F<76!6 z#NVM~xeeBOHm~Ytt-Oz*(Kmd!e|<1N6w>_(?r1^;K*GJW`-)z;Sl;~dGtp1Rm0aUq zyZdhat9^GCKdakiAK~IJes#Ur>s