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);
}
--- /dev/null
+New: The GridIn class can now read ExodusII files when deal.II is configured
+with Trilinos and SEACAS.
+<br>
+(David Wells, 2020/09/09)
#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
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*/
#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
# undef AssertThrowMPI
#endif // #ifdef AssertThrowMPI
+#ifdef AssertThrowExodusII
+# undef AssertThrowExodusII
+#endif // #ifdef AssertThrowExodusII
+
#ifdef AssertVectorVectorDimension
# undef AssertVectorVectorDimension
#endif // #ifdef AssertVectorVectorDimension
vtu,
/// Use read_assimp()
assimp,
+ /// Use read_exodusii()
+ exodusii,
};
/**
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.
*/
{{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];
}
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;
+ }
};
(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;
+ }
};
(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;
+ }
};
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;
+ }
};
(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];
+ }
};
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];
+ }
};
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];
+ }
};
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];
+ }
};
(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];
+ }
};
/**
# include <mpi.h>
#endif
+#ifdef DEAL_II_TRILINOS_WITH_SEACAS
+# include <exodusII.h>
+#endif
+
#ifdef DEAL_II_HAVE_GLIBC_STACKTRACE
# include <execinfo.h>
#endif
-#ifdef DEAL_II_WITH_MPI
namespace StandardExceptions
{
+#ifdef DEAL_II_WITH_MPI
ExcMPI::ExcMPI(const int error_code)
: error_code(error_code)
{}
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
# include <assimp/scene.h> // Output data structure
#endif
+#ifdef DEAL_II_TRILINOS_WITH_SEACAS
+# include <exodusII.h>
+#endif
+
DEAL_II_NAMESPACE_OPEN
#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
{
read_assimp(name);
}
+ else if (format == exodusii)
+ {
+ read_exodusii(name);
+ }
else
{
std::ifstream in(name.c_str());
"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;
}
{
case dbmesh:
return ".dbmesh";
+ case exodusii:
+ return ".e";
case msh:
return ".msh";
case vtk:
if (format_name == "dbmesh")
return dbmesh;
+ if (format_name == "exodusii")
+ return exodusii;
+
if (format_name == "msh")
return msh;
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";
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff