diff --git a/include/mesh/exodusII_io_helper.h b/include/mesh/exodusII_io_helper.h index 30aee769ce..4e0566decd 100644 --- a/include/mesh/exodusII_io_helper.h +++ b/include/mesh/exodusII_io_helper.h @@ -225,6 +225,11 @@ class ExodusII_IO_Helper : public ParallelObject */ void read_elem_in_block(int block); + /** + * Reads NSIDED face blocks used by NFACED element blocks. + */ + void read_face_blocks(); + /** * Read in edge blocks, storing information in the BoundaryInfo object. */ @@ -703,6 +708,13 @@ class ExodusII_IO_Helper : public ParallelObject // each block must equal num_edge. int & num_edge_blk; + // Total number of faces + int & num_face; + + // Total number of face blocks. The sum of the number of faces in + // each block must equal num_face. + int & num_face_blk; + // Total number of node sets int & num_node_sets; @@ -748,6 +760,15 @@ class ExodusII_IO_Helper : public ParallelObject // Vector of nodes in an element std::vector connect; + // For NSIDED element blocks, the number of nodes in each element + std::vector elem_node_counts; + + // For NFACED element blocks, the number of faces in each element + std::vector elem_face_counts; + + // For NFACED element blocks, the nodes in each Exodus face + std::vector> c0polyhedron_face_connect; + // Vector of the sideset IDs std::vector ss_ids; diff --git a/include/mesh/exodus_header_info.h b/include/mesh/exodus_header_info.h index bc7d1f7b97..4bc43a6429 100644 --- a/include/mesh/exodus_header_info.h +++ b/include/mesh/exodus_header_info.h @@ -57,7 +57,8 @@ class ExodusHeaderInfo // Pack individual integers into vector std::vector buffer = {num_dim, num_elem, num_elem_blk, num_node_sets, - num_side_sets, num_elem_sets, num_edge_blk, num_edge}; + num_side_sets, num_elem_sets, num_edge_blk, num_edge, + num_face_blk, num_face}; // broadcast integers comm.broadcast(buffer); @@ -72,6 +73,8 @@ class ExodusHeaderInfo num_elem_sets = buffer[ctr++]; num_edge_blk = buffer[ctr++]; num_edge = buffer[ctr++]; + num_face_blk = buffer[ctr++]; + num_face = buffer[ctr++]; } std::vector title; @@ -84,6 +87,8 @@ class ExodusHeaderInfo int num_elem_sets; int num_edge_blk; int num_edge; + int num_face_blk; + int num_face; }; } // namespace libMesh diff --git a/src/geom/cell_c0polyhedron.C b/src/geom/cell_c0polyhedron.C index ce5a207729..ba1f87f0f9 100644 --- a/src/geom/cell_c0polyhedron.C +++ b/src/geom/cell_c0polyhedron.C @@ -18,6 +18,7 @@ // Local includes #include "libmesh/cell_c0polyhedron.h" +#include "libmesh/enum_io_package.h" #include "libmesh/enum_order.h" #include "libmesh/face_polygon.h" #include "libmesh/mesh_tools.h" @@ -188,11 +189,44 @@ C0Polyhedron::edges_adjacent_to_node(const unsigned int n) const -void C0Polyhedron::connectivity(const unsigned int /*sf*/, - const IOPackage /*iop*/, - std::vector & /*conn*/) const +void C0Polyhedron::connectivity(const unsigned int sf, + const IOPackage iop, + std::vector & conn) const { - libmesh_not_implemented(); + libmesh_assert_less (sf, this->n_sub_elem()); + libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE); + + const auto & subtet = this->_triangulation[sf]; + + switch (iop) + { + case TECPLOT: + { + conn.resize(8); + conn[0] = this->node_id(subtet[0])+1; + conn[1] = this->node_id(subtet[1])+1; + conn[2] = this->node_id(subtet[2])+1; + conn[3] = this->node_id(subtet[2])+1; + conn[4] = this->node_id(subtet[3])+1; + conn[5] = this->node_id(subtet[3])+1; + conn[6] = this->node_id(subtet[3])+1; + conn[7] = this->node_id(subtet[3])+1; + return; + } + + case VTK: + { + conn.resize(4); + conn[0] = this->node_id(subtet[0]); + conn[1] = this->node_id(subtet[1]); + conn[2] = this->node_id(subtet[2]); + conn[3] = this->node_id(subtet[3]); + return; + } + + default: + libmesh_error_msg("Unsupported IO package " << iop); + } } diff --git a/src/mesh/exodusII_io.C b/src/mesh/exodusII_io.C index 4bf7bcf164..be706be80b 100644 --- a/src/mesh/exodusII_io.C +++ b/src/mesh/exodusII_io.C @@ -20,6 +20,7 @@ #include "libmesh/exodusII_io.h" #include "libmesh/boundary_info.h" +#include "libmesh/cell_c0polyhedron.h" #include "libmesh/dof_map.h" #include "libmesh/dyna_io.h" // ElementDefinition for BEX #include "libmesh/enum_elem_type.h" @@ -27,6 +28,7 @@ #include "libmesh/enum_to_string.h" #include "libmesh/equation_systems.h" #include "libmesh/exodusII_io_helper.h" +#include "libmesh/face_c0polygon.h" #include "libmesh/fpe_disabler.h" #include "libmesh/int_range.h" #include "libmesh/libmesh_logging.h" @@ -452,20 +454,89 @@ void ExodusII_IO::read (const std::string & fname) // Set any relevant node/edge maps for this element const std::string type_str (exio_helper->get_elem_type()); const auto & conv = exio_helper->get_conversion(type_str); + const bool is_c0polygon = (conv.libmesh_elem_type() == C0POLYGON); + const bool is_c0polyhedron = (conv.libmesh_elem_type() == C0POLYHEDRON); + std::size_t c0polygon_connect_offset = 0; + std::size_t c0polyhedron_connect_offset = 0; // Loop over all the faces in this block int jmax = nelem_last_block+exio_helper->num_elem_this_blk; for (int j=nelem_last_block; j uelem; + + if (is_c0polygon) + { + const int n_nodes = exio_helper->elem_node_counts[elem_num]; + libmesh_error_msg_if(n_nodes < 3, + "Error: Exodus NSIDED block element " + << elem_num + << " has only " << n_nodes << " nodes."); + + uelem = std::make_unique(cast_int(n_nodes)); + } + else if (is_c0polyhedron) + { + const int n_faces = exio_helper->elem_face_counts[elem_num]; + libmesh_error_msg_if(n_faces < 4, + "Error: Exodus NFACED block element " + << elem_num + << " has only " << n_faces << " faces."); + + std::vector> + sides(cast_int(n_faces)); + + for (auto s : index_range(sides)) + { + libmesh_assert_less(c0polyhedron_connect_offset, + exio_helper->connect.size()); + const int exodus_face_id = + exio_helper->connect[c0polyhedron_connect_offset++]; + + libmesh_error_msg_if + (exodus_face_id <= 0 || + exodus_face_id > + cast_int(exio_helper->c0polyhedron_face_connect.size()), + "Error: Exodus NFACED block element " + << elem_num + << " references face ID " << exodus_face_id + << ", but the file contains " + << exio_helper->c0polyhedron_face_connect.size() + << " faces."); + + const auto & face_nodes = + exio_helper->c0polyhedron_face_connect[exodus_face_id - 1]; + auto side = std::make_shared + (cast_int(face_nodes.size())); + + for (auto n : index_range(face_nodes)) + { + const auto libmesh_node_id = + exio_helper->get_libmesh_node_id(face_nodes[n]); + side->set_node(n, mesh.node_ptr(libmesh_node_id)); + } + + sides[s] = std::move(side); + } + + std::unique_ptr mid_elem_node; + uelem = std::make_unique(sides, mid_elem_node); + if (mid_elem_node) + { + Node * added_node = mesh.add_node(std::move(mid_elem_node)); + if (added_node->id() >= n_nodes) + n_nodes = added_node->id() + 1; + } + } + else + uelem = Elem::build(conv.libmesh_elem_type()); // Make sure that Exodus's number of nodes per Elem matches // the number of Nodes for this type of Elem. We only check // this for the first Elem in each block, since these values // are the same for every Elem in the block. - if (!elem_num) + if (!is_c0polygon && !is_c0polyhedron && !elem_num) libmesh_error_msg_if(exio_helper->num_nodes_per_elem != static_cast(uelem->n_nodes()), "Error: Exodus file says " << exio_helper->num_nodes_per_elem @@ -549,12 +620,26 @@ void ExodusII_IO::read (const std::string & fname) // If we don't have any Bezier extraction operators, this // is easy: we've already built all our nodes and just need // to link to them. - if (exio_helper->bex_cv_conn.empty()) + if (is_c0polyhedron) + { + // Node pointers were already assigned while constructing + // the polygonal sides above. + } + else if (exio_helper->bex_cv_conn.empty()) { - for (int k=0; knum_nodes_per_elem; k++) + const int n_nodes_this_elem = + is_c0polygon ? + exio_helper->elem_node_counts[elem_num] : + exio_helper->num_nodes_per_elem; + + for (int k=0; knum_nodes_per_elem + conv.get_node_map(k); + std::size_t gi = + is_c0polygon ? + c0polygon_connect_offset++ : + elem_num * exio_helper->num_nodes_per_elem + conv.get_node_map(k); + libmesh_assert_less(gi, exio_helper->connect.size()); // Get the 1-based Exodus node id from the "connect" array auto exodus_node_id = exio_helper->connect[gi]; @@ -689,6 +774,11 @@ void ExodusII_IO::read (const std::string & fname) } } + libmesh_assert(!is_c0polygon || + c0polygon_connect_offset == exio_helper->connect.size()); + libmesh_assert(!is_c0polyhedron || + c0polyhedron_connect_offset == exio_helper->connect.size()); + // running sum of # of elements per block, // (should equal total number of elements in the end) nelem_last_block += exio_helper->num_elem_this_blk; diff --git a/src/mesh/exodusII_io_helper.C b/src/mesh/exodusII_io_helper.C index 13f74055cf..e57b80e6ff 100644 --- a/src/mesh/exodusII_io_helper.C +++ b/src/mesh/exodusII_io_helper.C @@ -53,6 +53,7 @@ extern "C" { #include // workaround for HDF5 bug #include // std::strtol #include +#include #include // Anonymous namespace for file local data and helper functions @@ -66,6 +67,8 @@ static constexpr int libmesh_max_str_length = MAX_LINE_LENGTH; using namespace libMesh; +constexpr int c0polyhedron_face_block_id = 1; + // File scope constant node/edge/face mapping arrays. // 2D inverse face map definitions. // These take a libMesh ID and turn it into an Exodus ID @@ -274,6 +277,8 @@ ExodusII_IO_Helper::ExodusII_IO_Helper(const ParallelObject & parent, num_elem_blk(header_info.num_elem_blk), num_edge(header_info.num_edge), num_edge_blk(header_info.num_edge_blk), + num_face(header_info.num_face), + num_face_blk(header_info.num_face_blk), num_node_sets(header_info.num_node_sets), num_side_sets(header_info.num_side_sets), num_elem_sets(header_info.num_elem_sets), @@ -361,6 +366,14 @@ void ExodusII_IO_Helper::init_conversion_map() convert_type(QUAD4, "QUAD4"); convert_type(QUAD8, "QUAD8"); convert_type(QUAD9, "QUAD9"); + convert_type(C0POLYGON, "NSIDED"); + { + auto & conv = conversion_map[3][C0POLYHEDRON]; + conv.libmesh_type = C0POLYHEDRON; + conv.exodus_type = "NFACED"; + conv.dim = 3; + conv.n_nodes = 0; + } convert_type(QUADSHELL4, "SHELL4", nullptr, nullptr, nullptr, /* inverse_side_map = */ &quadshell4_inverse_edge_map, nullptr, nullptr, /* shellface_index_offset = */ 2); @@ -477,6 +490,10 @@ void ExodusII_IO_Helper::init_element_equivalence_map() // QUADSHELL9 equivalences element_equivalence_map["SHELL9"] = QUADSHELL9; + // Runtime-topology polytope equivalences + element_equivalence_map["NSIDED"] = C0POLYGON; + element_equivalence_map["NFACED"] = C0POLYHEDRON; + // TRI3 equivalences element_equivalence_map["TRI"] = TRI3; element_equivalence_map["TRI3"] = TRI3; @@ -706,6 +723,9 @@ void ExodusII_IO_Helper::open(const char * filename, bool read_only) if (read_only) { opened_for_reading = true; + elem_node_counts.clear(); + elem_face_counts.clear(); + c0polyhedron_face_connect.clear(); // ExodusII reads truncate to 32-char strings by default; we'd // like to support whatever's in the file, so as early as possible @@ -756,6 +776,8 @@ ExodusII_IO_Helper::read_header() const h.num_elem_sets = params.num_elem_sets; h.num_edge_blk = params.num_edge_blk; h.num_edge = params.num_edge; + h.num_face_blk = params.num_face_blk; + h.num_face = params.num_face; // And return it return h; @@ -1171,6 +1193,124 @@ std::string ExodusII_IO_Helper::get_node_set_name(int index) } +void ExodusII_IO_Helper::read_face_blocks() +{ + LOG_SCOPE("read_face_blocks()", "ExodusII_IO_Helper"); + + if (!c0polyhedron_face_connect.empty()) + return; + + libmesh_error_msg_if(num_face_blk == 0, + "Error: Exodus NFACED element block found, " + "but the file has no face blocks."); + + std::vector face_block_ids(num_face_blk); + ex_err = exII::ex_get_ids(ex_id, + exII::EX_FACE_BLOCK, + face_block_ids.data()); + EX_CHECK_ERR(ex_err, "Error getting face block IDs."); + + c0polyhedron_face_connect.clear(); + c0polyhedron_face_connect.reserve(num_face); + + for (auto block : index_range(face_block_ids)) + { + std::vector face_type(libmesh_max_str_length+1); + int num_face_this_blk = 0; + int num_node_data_this_blk = 0; + int num_edges_per_face = 0; + int num_faces_per_face = 0; + int num_attr_face = 0; + + ex_err = exII::ex_get_block(ex_id, + exII::EX_FACE_BLOCK, + face_block_ids[block], + face_type.data(), + &num_face_this_blk, + &num_node_data_this_blk, + &num_edges_per_face, + &num_faces_per_face, + &num_attr_face); + EX_CHECK_ERR(ex_err, "Error getting face block info."); + + const auto & conv = get_conversion(std::string(face_type.data())); + libmesh_error_msg_if(conv.libmesh_elem_type() != C0POLYGON, + "Error: NFACED polyhedron input currently expects " + "NSIDED face blocks, but face block " + << face_block_ids[block] << " has Exodus type " + << face_type.data() << "."); + + libmesh_error_msg_if + (!(num_edges_per_face == 0) && !(num_edges_per_face == -1), + "Error: Exodus NSIDED face block " + << face_block_ids[block] + << " has edge connectivity, which NFACED polyhedron input " + << "does not currently support."); + libmesh_error_msg_if + (!(num_faces_per_face == 0) && !(num_faces_per_face == -1), + "Error: Exodus NSIDED face block " + << face_block_ids[block] + << " has face-in-face connectivity, which NFACED polyhedron " + << "input does not currently support."); + + std::vector face_node_counts(num_face_this_blk); + if (!face_node_counts.empty()) + { + ex_err = exII::ex_get_entity_count_per_polyhedra + (ex_id, + exII::EX_FACE_BLOCK, + face_block_ids[block], + face_node_counts.data()); + EX_CHECK_ERR(ex_err, "Error reading polyhedron face node counts"); + } + + int counted_nodes = 0; + for (const auto count : face_node_counts) + counted_nodes += count; + + libmesh_error_msg_if(counted_nodes != num_node_data_this_blk, + "Error: Exodus NSIDED face block " + << face_block_ids[block] + << " says it has " << num_node_data_this_blk + << " total node entries, but its per-face " + << "node counts sum to " << counted_nodes << "."); + + std::vector face_connect(num_node_data_this_blk); + if (!face_connect.empty()) + { + ex_err = exII::ex_get_conn(ex_id, + exII::EX_FACE_BLOCK, + face_block_ids[block], + face_connect.data(), + nullptr, + nullptr); + EX_CHECK_ERR(ex_err, "Error reading polyhedron face connectivity."); + } + + std::size_t offset = 0; + for (const auto count : face_node_counts) + { + libmesh_error_msg_if(count < 3, + "Error: Exodus NSIDED face block " + << face_block_ids[block] + << " has a face with only " + << count << " nodes."); + + c0polyhedron_face_connect.emplace_back + (face_connect.begin() + offset, + face_connect.begin() + offset + count); + offset += count; + } + } + + libmesh_error_msg_if(c0polyhedron_face_connect.size() != + cast_int(num_face), + "Error: Exodus file says it has " + << num_face << " faces, but its face blocks contain " + << c0polyhedron_face_connect.size() << " faces."); +} + + void ExodusII_IO_Helper::read_elem_in_block(int block) @@ -1178,6 +1318,8 @@ void ExodusII_IO_Helper::read_elem_in_block(int block) LOG_SCOPE("read_elem_in_block()", "ExodusII_IO_Helper"); libmesh_assert_less (block, block_ids.size()); + elem_node_counts.clear(); + elem_face_counts.clear(); // Unlike the other "extended" APIs, this one does not use a parameter struct. int num_edges_per_elem = 0; @@ -1196,12 +1338,19 @@ void ExodusII_IO_Helper::read_elem_in_block(int block) EX_CHECK_ERR(ex_err, "Error getting block info."); message("Info retrieved successfully for block: ", block); - // Warn that we don't currently support reading blocks with extended info. + const auto & conv = get_conversion(std::string(elem_type.data())); + const bool is_c0polygon = (conv.libmesh_elem_type() == C0POLYGON); + const bool is_c0polyhedron = (conv.libmesh_elem_type() == C0POLYHEDRON); + + // Warn or error when we don't currently support reading blocks with extended info. // Note: the docs say -1 will be returned for this but I found that it was // actually 0, so not sure which it will be in general. - if (!(num_edges_per_elem == 0) && !(num_edges_per_elem == -1)) + if (is_c0polyhedron && !(num_edges_per_elem == 0) && !(num_edges_per_elem == -1)) + libmesh_error_msg("Error: Exodus NFACED element blocks with edge " + "connectivity are not currently supported."); + else if (!(num_edges_per_elem == 0) && !(num_edges_per_elem == -1)) libmesh_warning("Exodus files with extended edge connectivity not currently supported."); - if (!(num_faces_per_elem == 0) && !(num_faces_per_elem == -1)) + if (!is_c0polyhedron && !(num_faces_per_elem == 0) && !(num_faces_per_elem == -1)) libmesh_warning("Exodus files with extended face connectivity not currently supported."); // If we have a Bezier element here, then we've packed constraint @@ -1209,32 +1358,100 @@ void ExodusII_IO_Helper::read_elem_in_block(int block) // num_nodes_per_elem reflected both. const bool is_bezier = is_bezier_elem(elem_type.data()); if (is_bezier) + num_nodes_per_elem = conv.n_nodes; + else if (is_c0polygon) { - const auto & conv = get_conversion(std::string(elem_type.data())); - num_nodes_per_elem = conv.n_nodes; + elem_node_counts.resize(num_elem_this_blk); + + if (!elem_node_counts.empty()) + { + ex_err = exII::ex_get_entity_count_per_polyhedra + (ex_id, + exII::EX_ELEM_BLOCK, + block_ids[block], + elem_node_counts.data()); + EX_CHECK_ERR(ex_err, "Error reading polygon node counts"); + } + + int counted_nodes = 0; + for (const auto count : elem_node_counts) + counted_nodes += count; + + libmesh_error_msg_if(counted_nodes != num_node_data_per_elem, + "Error: Exodus NSIDED block " + << block_ids[block] + << " says it has " << num_node_data_per_elem + << " total node entries, but its per-element " + << "node counts sum to " << counted_nodes << "."); + + num_nodes_per_elem = 0; + } + else if (is_c0polyhedron) + { + if (c0polyhedron_face_connect.empty()) + this->read_face_blocks(); + + elem_face_counts.resize(num_elem_this_blk); + + if (!elem_face_counts.empty()) + { + ex_err = exII::ex_get_entity_count_per_polyhedra + (ex_id, + exII::EX_ELEM_BLOCK, + block_ids[block], + elem_face_counts.data()); + EX_CHECK_ERR(ex_err, "Error reading polyhedron face counts"); + } + + int counted_faces = 0; + for (const auto count : elem_face_counts) + counted_faces += count; + + libmesh_error_msg_if(counted_faces != num_faces_per_elem, + "Error: Exodus NFACED block " + << block_ids[block] + << " says it has " << num_faces_per_elem + << " total face entries, but its per-element " + << "face counts sum to " << counted_faces << "."); + + num_nodes_per_elem = 0; } else num_nodes_per_elem = num_node_data_per_elem; if (verbose) - libMesh::out << "Read a block of " << num_elem_this_blk - << " " << elem_type.data() << "(s)" - << " having " << num_nodes_per_elem - << " nodes per element." << std::endl; + { + libMesh::out << "Read a block of " << num_elem_this_blk + << " " << elem_type.data() << "(s)"; + if (is_c0polygon) + libMesh::out << " having " << num_node_data_per_elem + << " total node entries."; + else if (is_c0polyhedron) + libMesh::out << " having " << num_faces_per_elem + << " total face entries."; + else + libMesh::out << " having " << num_nodes_per_elem + << " nodes per element."; + libMesh::out << std::endl; + } // Read in the connectivity of the elements of this block, // watching out for the case where we actually have no // elements in this block (possible with parallel files) - connect.resize(num_node_data_per_elem*num_elem_this_blk); + connect.resize(is_c0polygon ? + num_node_data_per_elem : + is_c0polyhedron ? + num_faces_per_elem : + num_node_data_per_elem*num_elem_this_blk); if (!connect.empty()) { ex_err = exII::ex_get_conn(ex_id, exII::EX_ELEM_BLOCK, block_ids[block], - connect.data(), // node_conn - nullptr, // elem_edge_conn (unused) - nullptr); // elem_face_conn (unused) + is_c0polyhedron ? nullptr : connect.data(), + nullptr, // elem_edge_conn (unused) + is_c0polyhedron ? connect.data() : nullptr); EX_CHECK_ERR(ex_err, "Error reading block connectivity."); message("Connectivity retrieved successfully for block: ", block); @@ -2361,6 +2578,25 @@ void ExodusII_IO_Helper::initialize(std::string str_title, const MeshBase & mesh num_elem = n_active_elem; num_nodes = 0; + num_face = 0; + num_face_blk = 0; + + dof_id_type local_num_c0polyhedron_faces = 0; + int has_c0polyhedron = 0; + for (const auto & elem : mesh.active_local_element_ptr_range()) + if (elem->type() == C0POLYHEDRON) + { + has_c0polyhedron = 1; + local_num_c0polyhedron_faces += elem->n_sides(); + } + + mesh.comm().sum(local_num_c0polyhedron_faces); + mesh.comm().max(has_c0polyhedron); + if (has_c0polyhedron) + { + num_face = cast_int(local_num_c0polyhedron_faces); + num_face_blk = 1; + } // If we're adding face elements they'll need copies of their nodes. // We also have to count of how many nodes (and gaps between nodes!) @@ -2534,6 +2770,8 @@ void ExodusII_IO_Helper::initialize(std::string str_title, const MeshBase & mesh params.num_elem_sets = num_elem_sets; params.num_edge_blk = num_edge_blk; params.num_edge = num_edge; + params.num_face_blk = num_face_blk; + params.num_face = num_face; ex_err = exII::ex_put_init_ext(ex_id, ¶ms); EX_CHECK_ERR(ex_err, "Error initializing new Exodus file."); @@ -2729,6 +2967,10 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti NamesData names_table(num_elem_blk, _max_name_length); num_elem = 0; + bool has_c0polygon_blocks = false; + bool has_c0polyhedron_blocks = false; + int c0polyhedron_total_faces = 0; + int c0polyhedron_total_face_nodes = 0; // counter indexes into the block_ids vector unsigned int counter = 0; @@ -2754,8 +2996,13 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti libmesh_assert(!element_id_vec.empty()); num_elem_this_blk_vec.push_back (cast_int(element_id_vec.size())); - names_table.push_back_entry - (mesh.subdomain_name(subdomain_id)); + + std::string block_name = mesh.subdomain_name(subdomain_id); + if (block_name.empty() && elem_t == C0POLYGON) + block_name = "NSIDED_" + std::to_string(counter + 1); + if (block_name.empty() && elem_t == C0POLYHEDRON) + block_name = "NFACED_" + std::to_string(counter + 1); + names_table.push_back_entry(block_name); } num_elem += num_elem_this_blk_vec.back(); @@ -2764,19 +3011,79 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti // Note that Exodus assumes all elements in a block are of the same type! // We are using that same assumption here! const auto & conv = get_conversion(elem_t); - num_nodes_per_elem = Elem::type_to_n_nodes_map[elem_t]; - if (Elem::type_to_n_nodes_map[elem_t] == invalid_uint) - libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented"); + int num_edges_per_elem = 0; + int num_faces_per_elem = 0; + if (elem_t == C0POLYGON) + { + if (subdomain_id >= subdomain_id_end) + libmesh_not_implemented_msg("Support for C0POLYGON side blocks not yet implemented"); + + has_c0polygon_blocks = true; + num_nodes_per_elem = 0; + + for (auto elem_id : element_id_vec) + { + const Elem & elem = mesh.elem_ref(elem_id); + + libmesh_error_msg_if(elem.type() != C0POLYGON, + "Error: Exodus requires all elements with a given subdomain ID " + "to be the same type.\n" + << "Can't write both " + << Utility::enum_to_string(elem.type()) + << " and C0POLYGON in the same block!"); + + num_nodes_per_elem += cast_int(elem.n_nodes()); + } + } + else if (elem_t == C0POLYHEDRON) + { + if (subdomain_id >= subdomain_id_end) + libmesh_not_implemented_msg("Support for C0POLYHEDRON side blocks not yet implemented"); + + has_c0polyhedron_blocks = true; + num_nodes_per_elem = 0; + + for (auto elem_id : element_id_vec) + { + const Elem & elem = mesh.elem_ref(elem_id); + + libmesh_error_msg_if(elem.type() != C0POLYHEDRON, + "Error: Exodus requires all elements with a given subdomain ID " + "to be the same type.\n" + << "Can't write both " + << Utility::enum_to_string(elem.type()) + << " and C0POLYHEDRON in the same block!"); + + const int elem_n_sides = cast_int(elem.n_sides()); + num_faces_per_elem += elem_n_sides; + c0polyhedron_total_faces += elem_n_sides; + + for (auto s : elem.side_index_range()) + c0polyhedron_total_face_nodes += cast_int(elem.nodes_on_side(s).size()); + } + } + else + { + num_nodes_per_elem = Elem::type_to_n_nodes_map[elem_t]; + if (Elem::type_to_n_nodes_map[elem_t] == invalid_uint) + libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented"); + } elem_blk_id.push_back(subdomain_id); elem_type_table.push_back_entry(conv.exodus_elem_type().c_str()); num_nodes_per_elem_vec.push_back(num_nodes_per_elem); num_attr_vec.push_back(0); // we don't currently use elem block attributes. - num_edges_per_elem_vec.push_back(0); // We don't currently store any edge blocks - num_faces_per_elem_vec.push_back(0); // We don't currently store any face blocks + num_edges_per_elem_vec.push_back(num_edges_per_elem); // We don't currently store any edge blocks + num_faces_per_elem_vec.push_back(num_faces_per_elem); ++counter; } + if (has_c0polyhedron_blocks) + { + libmesh_assert_equal_to(num_face_blk, 1); + libmesh_assert_equal_to(num_face, c0polyhedron_total_faces); + } + // Here we reserve() space so that we can push_back() onto the // elem_num_map in the loops below. this->elem_num_map.reserve(num_elem); @@ -2915,6 +3222,9 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti // We also build a data structure of edge block names which can // later be passed to exII::ex_put_names(). NamesData edge_block_names_table(num_edge_blk, _max_name_length); + NamesData face_block_names_table(num_face_blk, _max_name_length); + if (has_c0polyhedron_blocks) + face_block_names_table.push_back_entry("NSIDED_FACES"); // Note: We are going to use the edge **boundary** ids as **block** ids. for (const auto & pr : edge_id_to_conn) @@ -2940,41 +3250,91 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti edge_block_names_table.push_back_entry(bi.get_edgeset_name(id)); } - // Zero-initialize and then fill in an exII::ex_block_params struct - // with the data we have collected. This new API replaces the old - // exII::ex_put_concat_elem_block() API, and will eventually allow - // us to also allocate space for edge/face blocks if desired. - // - // TODO: It seems like we should be able to take advantage of the - // optimization where you set define_maps==1, but when I tried this - // I got the error: "failed to find node map size". I think the - // problem is that we need to first specify a nonzero number of - // node/elem maps during the call to ex_put_init_ext() in order for - // this to work correctly. - exII::ex_block_params params = {}; - - // Set pointers for information about elem blocks. - params.elem_blk_id = elem_blk_id.data(); - params.elem_type = elem_type_table.get_char_star_star(); - params.num_elem_this_blk = num_elem_this_blk_vec.data(); - params.num_nodes_per_elem = num_nodes_per_elem_vec.data(); - params.num_edges_per_elem = num_edges_per_elem_vec.data(); - params.num_faces_per_elem = num_faces_per_elem_vec.data(); - params.num_attr_elem = num_attr_vec.data(); - params.define_maps = 0; - - // Set pointers to edge block information only if we actually have some. - if (num_edge_blk) + if (has_c0polygon_blocks || has_c0polyhedron_blocks) { - params.edge_blk_id = edge_blk_id.data(); - params.edge_type = edge_type_table.get_char_star_star(); - params.num_edge_this_blk = num_edge_this_blk_vec.data(); - params.num_nodes_per_edge = num_nodes_per_edge_vec.data(); - params.num_attr_edge = num_attr_edge_vec.data(); + // ex_put_concat_all_blocks() does not define the per-polytope + // entity count arrays required by NSIDED/NFACED blocks in all supported + // Exodus versions. Define blocks individually through ex_put_block(). + if (has_c0polyhedron_blocks) + { + ex_err = exII::ex_put_block(ex_id, + exII::EX_FACE_BLOCK, + c0polyhedron_face_block_id, + "NSIDED", + c0polyhedron_total_faces, + c0polyhedron_total_face_nodes, + 0, + 0, + 0); + EX_CHECK_ERR(ex_err, "Error writing polyhedron face block."); + } + + for (auto i : index_range(elem_blk_id)) + { + ex_err = exII::ex_put_block(ex_id, + exII::EX_ELEM_BLOCK, + elem_blk_id[i], + elem_type_table.get_char_star(cast_int(i)), + num_elem_this_blk_vec[i], + num_nodes_per_elem_vec[i], + num_edges_per_elem_vec[i], + num_faces_per_elem_vec[i], + num_attr_vec[i]); + EX_CHECK_ERR(ex_err, "Error writing element block."); + } + + for (auto i : index_range(edge_blk_id)) + { + ex_err = exII::ex_put_block(ex_id, + exII::EX_EDGE_BLOCK, + edge_blk_id[i], + edge_type_table.get_char_star(cast_int(i)), + num_edge_this_blk_vec[i], + num_nodes_per_edge_vec[i], + 0, + 0, + num_attr_edge_vec[i]); + EX_CHECK_ERR(ex_err, "Error writing edge block."); + } } + else + { + // Zero-initialize and then fill in an exII::ex_block_params struct + // with the data we have collected. This new API replaces the old + // exII::ex_put_concat_elem_block() API, and will eventually allow + // us to also allocate space for edge/face blocks if desired. + // + // TODO: It seems like we should be able to take advantage of the + // optimization where you set define_maps==1, but when I tried this + // I got the error: "failed to find node map size". I think the + // problem is that we need to first specify a nonzero number of + // node/elem maps during the call to ex_put_init_ext() in order for + // this to work correctly. + exII::ex_block_params params = {}; + + // Set pointers for information about elem blocks. + params.elem_blk_id = elem_blk_id.data(); + params.elem_type = elem_type_table.get_char_star_star(); + params.num_elem_this_blk = num_elem_this_blk_vec.data(); + params.num_nodes_per_elem = num_nodes_per_elem_vec.data(); + params.num_edges_per_elem = num_edges_per_elem_vec.data(); + params.num_faces_per_elem = num_faces_per_elem_vec.data(); + params.num_attr_elem = num_attr_vec.data(); + params.define_maps = 0; + + // Set pointers to edge block information only if we actually have some. + if (num_edge_blk) + { + params.edge_blk_id = edge_blk_id.data(); + params.edge_type = edge_type_table.get_char_star_star(); + params.num_edge_this_blk = num_edge_this_blk_vec.data(); + params.num_nodes_per_edge = num_nodes_per_edge_vec.data(); + params.num_attr_edge = num_attr_edge_vec.data(); + } - ex_err = exII::ex_put_concat_all_blocks(ex_id, ¶ms); - EX_CHECK_ERR(ex_err, "Error writing element blocks."); + ex_err = exII::ex_put_concat_all_blocks(ex_id, ¶ms); + EX_CHECK_ERR(ex_err, "Error writing element blocks."); + } // This counter is used to fill up the libmesh_elem_num_to_exodus map in the loop below. unsigned libmesh_elem_num_to_exodus_counter = 0; @@ -2993,6 +3353,12 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti next_fake_id = mesh.next_unique_id(); #endif + std::vector c0polyhedron_face_connect; + std::vector c0polyhedron_face_node_counts; + c0polyhedron_face_connect.reserve(c0polyhedron_total_face_nodes); + c0polyhedron_face_node_counts.reserve(c0polyhedron_total_faces); + int next_c0polyhedron_face_id = 1; + for (auto & [subdomain_id, element_id_vec] : subdomain_map) { // Use the first element in the block to get representative @@ -3004,15 +3370,39 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti mesh.elem_ref(element_id_vec[0]).type(); const auto & conv = get_conversion(elem_t); - num_nodes_per_elem = Elem::type_to_n_nodes_map[elem_t]; - if (Elem::type_to_n_nodes_map[elem_t] == invalid_uint) - libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented"); + const bool is_c0polygon_block = (elem_t == C0POLYGON); + const bool is_c0polyhedron_block = (elem_t == C0POLYHEDRON); + std::vector c0polygon_node_counts; + std::vector c0polyhedron_face_counts; + + if (is_c0polygon_block && subdomain_id >= subdomain_id_end) + libmesh_not_implemented_msg("Support for C0POLYGON side blocks not yet implemented"); + if (is_c0polyhedron_block && subdomain_id >= subdomain_id_end) + libmesh_not_implemented_msg("Support for C0POLYHEDRON side blocks not yet implemented"); + + if (!is_c0polygon_block && !is_c0polyhedron_block) + { + num_nodes_per_elem = Elem::type_to_n_nodes_map[elem_t]; + if (Elem::type_to_n_nodes_map[elem_t] == invalid_uint) + libmesh_not_implemented_msg("Support for Polygons/Polyhedra not yet implemented"); + } // If this is a *real* block, we just loop over vectors of // element ids to add. if (subdomain_id < subdomain_id_end) { - connect.resize(element_id_vec.size()*num_nodes_per_elem); + if (is_c0polygon_block) + { + connect.clear(); + c0polygon_node_counts.reserve(element_id_vec.size()); + } + else if (is_c0polyhedron_block) + { + connect.clear(); + c0polyhedron_face_counts.reserve(element_id_vec.size()); + } + else + connect.resize(element_id_vec.size()*num_nodes_per_elem); for (auto i : index_range(element_id_vec)) { @@ -3039,30 +3429,92 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti << Utility::enum_to_string(conv.libmesh_elem_type()) << " in the same block!"); - for (unsigned int j=0; j(num_nodes_per_elem); ++j) + if (is_c0polygon_block) { - unsigned int connect_index = cast_int((i*num_nodes_per_elem)+j); - unsigned elem_node_index = conv.get_inverse_node_map(j); // inverse node map is for writing. - if (!use_discontinuous) + c0polygon_node_counts.push_back(cast_int(elem.n_nodes())); + + for (auto elem_node_index : elem.node_index_range()) { - // The global id for the current node in libmesh. - dof_id_type libmesh_node_id = elem.node_id(elem_node_index); - - // Write the Exodus global node id associated with - // this libmesh node number to the connectivity - // array, or throw an error if it's not found. - connect[connect_index] = - libmesh_map_find(libmesh_node_num_to_exodus, - cast_int(libmesh_node_id)); + if (!use_discontinuous) + { + // The global id for the current node in libmesh. + dof_id_type libmesh_node_id = elem.node_id(elem_node_index); + + // Write the Exodus global node id associated with + // this libmesh node number to the connectivity + // array, or throw an error if it's not found. + connect.push_back + (libmesh_map_find(libmesh_node_num_to_exodus, + cast_int(libmesh_node_id))); + } + else + { + // Look up the (elem_id, elem_node_index) pair in the map. + connect.push_back + (libmesh_map_find(discontinuous_node_indices, + std::make_pair(elem_id, elem_node_index))); + } } - else + } + else if (is_c0polyhedron_block) + { + c0polyhedron_face_counts.push_back(cast_int(elem.n_sides())); + + for (auto s : elem.side_index_range()) { - // Look up the (elem_id, elem_node_index) pair in the map. - connect[connect_index] = - libmesh_map_find(discontinuous_node_indices, - std::make_pair(elem_id, elem_node_index)); + connect.push_back(next_c0polyhedron_face_id++); + + const std::vector side_nodes = + elem.nodes_on_side(s); + c0polyhedron_face_node_counts.push_back + (cast_int(side_nodes.size())); + + for (const auto elem_node_index : side_nodes) + { + if (!use_discontinuous) + { + const dof_id_type libmesh_node_id = + elem.node_id(elem_node_index); + c0polyhedron_face_connect.push_back + (libmesh_map_find(libmesh_node_num_to_exodus, + cast_int(libmesh_node_id))); + } + else + { + c0polyhedron_face_connect.push_back + (libmesh_map_find(discontinuous_node_indices, + std::make_pair(elem_id, elem_node_index))); + } + } } - } // end for(j) + } + else + { + for (unsigned int j=0; j(num_nodes_per_elem); ++j) + { + unsigned int connect_index = cast_int((i*num_nodes_per_elem)+j); + unsigned elem_node_index = conv.get_inverse_node_map(j); // inverse node map is for writing. + if (!use_discontinuous) + { + // The global id for the current node in libmesh. + dof_id_type libmesh_node_id = elem.node_id(elem_node_index); + + // Write the Exodus global node id associated with + // this libmesh node number to the connectivity + // array, or throw an error if it's not found. + connect[connect_index] = + libmesh_map_find(libmesh_node_num_to_exodus, + cast_int(libmesh_node_id)); + } + else + { + // Look up the (elem_id, elem_node_index) pair in the map. + connect[connect_index] = + libmesh_map_find(discontinuous_node_indices, + std::make_pair(elem_id, elem_node_index)); + } + } // end for(j) + } // push_back() either elem_id+1 or the current Elem's // unique_id+1 into the elem_num_map, depending on the value @@ -3126,12 +3578,54 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti (ex_id, exII::EX_ELEM_BLOCK, subdomain_id, - connect.data(), // node_conn - nullptr, // elem_edge_conn (unused) - nullptr); // elem_face_conn (unused) + is_c0polyhedron_block ? nullptr : connect.data(), // node_conn + nullptr, // elem_edge_conn (unused) + is_c0polyhedron_block ? connect.data() : nullptr); EX_CHECK_ERR(ex_err, "Error writing element connectivities"); + + if (is_c0polygon_block) + { + ex_err = exII::ex_put_entity_count_per_polyhedra + (ex_id, + exII::EX_ELEM_BLOCK, + subdomain_id, + c0polygon_node_counts.data()); + EX_CHECK_ERR(ex_err, "Error writing polygon node counts"); + } + if (is_c0polyhedron_block) + { + ex_err = exII::ex_put_entity_count_per_polyhedra + (ex_id, + exII::EX_ELEM_BLOCK, + subdomain_id, + c0polyhedron_face_counts.data()); + EX_CHECK_ERR(ex_err, "Error writing polyhedron face counts"); + } } // end for (auto & [subdomain_id, element_id_vec] : subdomain_map) + if (has_c0polyhedron_blocks) + { + libmesh_assert_equal_to(c0polyhedron_face_node_counts.size(), + cast_int(num_face)); + libmesh_assert_equal_to(next_c0polyhedron_face_id, num_face + 1); + + ex_err = exII::ex_put_conn + (ex_id, + exII::EX_FACE_BLOCK, + c0polyhedron_face_block_id, + c0polyhedron_face_connect.data(), // node_conn + nullptr, // elem_edge_conn (unused) + nullptr); // elem_face_conn (unused) + EX_CHECK_ERR(ex_err, "Error writing polyhedron face connectivities"); + + ex_err = exII::ex_put_entity_count_per_polyhedra + (ex_id, + exII::EX_FACE_BLOCK, + c0polyhedron_face_block_id, + c0polyhedron_face_node_counts.data()); + EX_CHECK_ERR(ex_err, "Error writing polyhedron face node counts"); + } + // write out the element number map that we created ex_err = exII::ex_put_elem_num_map(ex_id, elem_num_map.data()); EX_CHECK_ERR(ex_err, "Error writing element map"); @@ -3143,6 +3637,15 @@ void ExodusII_IO_Helper::write_elements(const MeshBase & mesh, bool use_disconti EX_CHECK_ERR(ex_err, "Error writing element block names"); } + if (num_face_blk > 0) + { + ex_err = exII::ex_put_names + (ex_id, + exII::EX_FACE_BLOCK, + face_block_names_table.get_char_star_star()); + EX_CHECK_ERR(ex_err, "Error writing face block names"); + } + // Write out edge blocks if we have any for (const auto & pr : edge_id_to_conn) { diff --git a/src/mesh/mesh_generation.C b/src/mesh/mesh_generation.C index 900b359b32..3848052eee 100644 --- a/src/mesh/mesh_generation.C +++ b/src/mesh/mesh_generation.C @@ -31,6 +31,7 @@ #include "libmesh/face_quad8.h" #include "libmesh/face_quad9.h" #include "libmesh/face_c0polygon.h" +#include "libmesh/cell_c0polyhedron.h" #include "libmesh/cell_hex8.h" #include "libmesh/cell_hex20.h" #include "libmesh/cell_hex27.h" @@ -57,6 +58,7 @@ #include "libmesh/enum_to_string.h" // C++ includes +#include #include // *must* precede for proper std:abs() on PGI, Sun Studio CC #include // for std::sqrt #include @@ -125,6 +127,7 @@ unsigned int idx(const ElemType type, case INVALID_ELEM: case HEX8: case PRISM6: + case C0POLYHEDRON: { return i + (nx+1)*(j + k*(ny+1)); } @@ -921,8 +924,8 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh, std::vector node_list; // Start with a layer of triangles on the boundary - const auto dx_tri = (xmax - xmin) / (nx); - const auto dy_tri = (ymax - ymin) / (ny + 1); + const auto dx_tri = Real(1) / nx; + const auto dy_tri = Real(1) / (ny + 1); std::unique_ptr new_elem; for (const auto i : make_range(nx + 1)) { @@ -973,7 +976,9 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh, // Build layers of hexagons const auto hex_side = - (ymax - ymin - (ny == 1 ? dy_tri : (1. + (ny - 1) / 2.) * dy_tri)) / ny; + (Real(1) - (ny == 1 ? + dy_tri : + (Real(1) + (ny - 1) / 2.) * dy_tri)) / ny; for (const auto j : make_range(ny)) { for (const auto i : make_range(nx + (j % 2))) @@ -1066,7 +1071,7 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh, Node *node0, *node1, *node2; if (i == 0 && ny_odd) { - node0 = mesh.add_point(Point(0., ymax, 0.)); + node0 = mesh.add_point(Point(0., 1., 0.)); node1 = node_list[running_index++]; node2 = node_list[running_index]; } @@ -1081,7 +1086,7 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh, { node0 = node_list[running_index++]; node1 = node_list[running_index]; - node2 = mesh.add_point(Point(xmax, ymax, 0.)); + node2 = mesh.add_point(Point(1., 1., 0.)); } new_elem = std::make_unique(3); @@ -1171,6 +1176,7 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh, case HEX8: case HEX20: case HEX27: + case C0POLYHEDRON: case TET4: // TET4's are created from an initial HEX27 discretization case TET10: // TET10's are created from an initial HEX27 discretization case TET14: // TET14's are created from an initial HEX27 discretization @@ -1207,8 +1213,14 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh, case INVALID_ELEM: case HEX8: case PRISM6: + case C0POLYHEDRON: { - mesh.reserve_nodes( (nx+1)*(ny+1)*(nz+1) ); + const dof_id_type grid_nodes = + cast_int((nx+1)*(ny+1)*(nz+1)); + mesh.reserve_nodes(grid_nodes + + (type == C0POLYHEDRON ? + cast_int(nx*ny*nz) : + 0)); break; } @@ -1267,6 +1279,7 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh, case INVALID_ELEM: case HEX8: case PRISM6: + case C0POLYHEDRON: { for (unsigned int k=0; k<=nz; k++) for (unsigned int j=0; j<=ny; j++) @@ -1415,6 +1428,69 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh, } + case C0POLYHEDRON: + { + const std::array, 6> side_nodes = + {{{0, 1, 2, 3}, // min z + {0, 1, 5, 4}, // min y + {2, 6, 5, 1}, // max x + {2, 3, 7, 6}, // max y + {0, 4, 7, 3}, // min x + {5, 6, 7, 4}}}; // max z + + for (unsigned int k=0; k elem_nodes = + {{mesh.node_ptr(idx(type,nx,ny,i,j,k) ), + mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ), + mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ), + mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ), + mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ), + mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ), + mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)), + mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) )}}; + + std::vector> sides(side_nodes.size()); + for (auto s : index_range(side_nodes)) + { + sides[s] = std::make_shared(side_nodes[s].size()); + for (auto n : index_range(side_nodes[s])) + sides[s]->set_node(n, elem_nodes[side_nodes[s][n]]); + } + + std::unique_ptr mid_elem_node; + std::unique_ptr new_elem = + std::make_unique(sides, mid_elem_node); + if (mid_elem_node) + mesh.add_node(std::move(mid_elem_node)); + + new_elem->set_id() = elem_id++; + Elem * elem = mesh.add_elem(std::move(new_elem)); + + if (k == 0) + boundary_info.add_side(elem, 0, 0); + + if (k == (nz-1)) + boundary_info.add_side(elem, 5, 5); + + if (j == 0) + boundary_info.add_side(elem, 1, 1); + + if (j == (ny-1)) + boundary_info.add_side(elem, 3, 3); + + if (i == 0) + boundary_info.add_side(elem, 4, 4); + + if (i == (nx-1)) + boundary_info.add_side(elem, 2, 2); + } + break; + } + + case PRISM6: diff --git a/tests/Makefile.am b/tests/Makefile.am index c1965e31c6..6427362120 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -415,6 +415,10 @@ CLEANFILES = cube_mesh.xda \ side_discontinuous_QUAD9_disc.e \ side_discontinuous_TRI6.e \ side_discontinuous_TRI6_disc.e \ + write_exodus_C0POLYGON.e \ + write_exodus_C0POLYHEDRON.e \ + write_exodus_C0POLYHEDRON_HEXPRISM.e \ + write_exodus_C0POLYHEDRON_HEXPRISM_READ.e \ write_exodus_EDGE2.e \ write_exodus_EDGE3.e \ write_exodus_EDGE4.e \ diff --git a/tests/Makefile.in b/tests/Makefile.in index 679d9c5524..2fc6d9bd86 100644 --- a/tests/Makefile.in +++ b/tests/Makefile.in @@ -2521,9 +2521,11 @@ CLEANFILES = cube_mesh.xda slit_mesh.xda slit_solution.xda out.e \ side_discontinuous_EDGE3_disc.e side_discontinuous_HEX27.e \ side_discontinuous_HEX27_disc.e side_discontinuous_QUAD9.e \ side_discontinuous_QUAD9_disc.e side_discontinuous_TRI6.e \ - side_discontinuous_TRI6_disc.e write_exodus_EDGE2.e \ - write_exodus_EDGE3.e write_exodus_EDGE4.e write_exodus_HEX20.e \ - write_exodus_HEX27.e write_exodus_HEX8.e \ + side_discontinuous_TRI6_disc.e write_exodus_C0POLYGON.e \ + write_exodus_C0POLYHEDRON.e write_exodus_C0POLYHEDRON_HEXPRISM.e \ + write_exodus_C0POLYHEDRON_HEXPRISM_READ.e \ + write_exodus_EDGE2.e write_exodus_EDGE3.e write_exodus_EDGE4.e \ + write_exodus_HEX20.e write_exodus_HEX27.e write_exodus_HEX8.e \ write_exodus_PRISM15.e write_exodus_PRISM18.e \ write_exodus_PRISM20.e write_exodus_PRISM21.e \ write_exodus_PRISM6.e write_exodus_PYRAMID13.e \ diff --git a/tests/mesh/exodus_test.C b/tests/mesh/exodus_test.C index 6effe982fa..2773e72a8e 100644 --- a/tests/mesh/exodus_test.C +++ b/tests/mesh/exodus_test.C @@ -2,10 +2,19 @@ #ifdef LIBMESH_HAVE_EXODUS_API +#include "libmesh/cell_c0polyhedron.h" #include "libmesh/enum_to_string.h" #include "libmesh/exodusII_io.h" +#include "libmesh/face_c0polygon.h" +#include "libmesh/int_range.h" #include "libmesh/mesh_communication.h" +#include "libmesh/mesh_generation.h" #include "libmesh/mesh_serializer.h" +#include "libmesh/node.h" + +#include +#include +#include using namespace libMesh; @@ -13,7 +22,6 @@ template class ExodusTest : public MeshPerElemTest { public: - void test_read_gold() { LOG_UNIT_TEST; @@ -64,26 +72,27 @@ public: } }; -#define EXODUSTEST \ - CPPUNIT_TEST( test_read_gold ); \ - CPPUNIT_TEST( test_write ); - -#define INSTANTIATE_EXODUSTEST(elemtype) \ - class ExodusTest_##elemtype : public ExodusTest { \ - public: \ - ExodusTest_##elemtype() : \ - ExodusTest() { \ - if (unitlog->summarized_logs_enabled()) \ - this->libmesh_suite_name = "ExodusTest"; \ - else \ - this->libmesh_suite_name = "ExodusTest_" #elemtype; \ - } \ - CPPUNIT_TEST_SUITE( ExodusTest_##elemtype ); \ - EXODUSTEST; \ - CPPUNIT_TEST_SUITE_END(); \ - }; \ - \ - CPPUNIT_TEST_SUITE_REGISTRATION( ExodusTest_##elemtype ) +#define EXODUSTEST \ + CPPUNIT_TEST(test_read_gold); \ + CPPUNIT_TEST(test_write); + +#define INSTANTIATE_EXODUSTEST(elemtype) \ + class ExodusTest_##elemtype : public ExodusTest \ + { \ + public: \ + ExodusTest_##elemtype() : ExodusTest() \ + { \ + if (unitlog->summarized_logs_enabled()) \ + this->libmesh_suite_name = "ExodusTest"; \ + else \ + this->libmesh_suite_name = "ExodusTest_" #elemtype; \ + } \ + CPPUNIT_TEST_SUITE(ExodusTest_##elemtype); \ + EXODUSTEST; \ + CPPUNIT_TEST_SUITE_END(); \ + }; \ + \ + CPPUNIT_TEST_SUITE_REGISTRATION(ExodusTest_##elemtype) INSTANTIATE_EXODUSTEST(EDGE2); INSTANTIATE_EXODUSTEST(EDGE3); @@ -101,9 +110,246 @@ INSTANTIATE_EXODUSTEST(QUAD8); INSTANTIATE_EXODUSTEST(QUADSHELL8); INSTANTIATE_EXODUSTEST(QUAD9); INSTANTIATE_EXODUSTEST(QUADSHELL9); + +class ExodusC0PolygonTest : public CppUnit::TestCase +{ +public: + LIBMESH_CPPUNIT_TEST_SUITE(ExodusC0PolygonTest); + + CPPUNIT_TEST(test_write_and_read_pentagon); + + CPPUNIT_TEST_SUITE_END(); + + void build_pentagon(Mesh &mesh) + { + const std::vector points = + { {0, 0}, {1, 0}, {1.5, 0.5}, {1, 1}, {0, 1} }; + + for (auto p : index_range(points)) + mesh.add_point(points[p], /*id=*/p); + + std::unique_ptr polygon = + std::make_unique(cast_int(points.size())); + for (auto i : index_range(points)) + polygon->set_node(i, mesh.node_ptr(i)); + + polygon->set_id() = 0; + Elem *elem = mesh.add_elem(std::move(polygon)); + elem->subdomain_id() = 1; + mesh.prepare_for_use(); + } + + void test_write_and_read_pentagon() + { + LOG_UNIT_TEST; + + Mesh mesh(*TestCommWorld); + this->build_pentagon(mesh); + + { + ExodusII_IO exii(mesh); + exii.write("write_exodus_C0POLYGON.e"); + } + + Mesh input_mesh(*TestCommWorld); + ExodusII_IO exii_input(input_mesh); + if (input_mesh.processor_id() == 0) + exii_input.read("write_exodus_C0POLYGON.e"); + + MeshCommunication().broadcast(input_mesh); + input_mesh.prepare_for_use(); + + CPPUNIT_ASSERT_EQUAL(cast_int(1), input_mesh.n_elem()); + + const Elem *elem = input_mesh.elem_ptr(0); + CPPUNIT_ASSERT(elem); + CPPUNIT_ASSERT_EQUAL(C0POLYGON, elem->type()); + CPPUNIT_ASSERT_EQUAL(5u, elem->n_nodes()); + + for (auto i : make_range(5)) + CPPUNIT_ASSERT_EQUAL(cast_int(i), elem->node_id(i)); + } +}; + +CPPUNIT_TEST_SUITE_REGISTRATION(ExodusC0PolygonTest); #endif // LIBMESH_DIM > 1 #if LIBMESH_DIM > 2 +class ExodusC0PolyhedronTest : public CppUnit::TestCase +{ +public: + LIBMESH_CPPUNIT_TEST_SUITE(ExodusC0PolyhedronTest); + + CPPUNIT_TEST(test_write_cube_header); + CPPUNIT_TEST(test_write_hexagonal_prism_header); + CPPUNIT_TEST(test_write_and_read_hexagonal_prism); + + CPPUNIT_TEST_SUITE_END(); + + ExodusHeaderInfo write_and_read_header(Mesh &mesh, + const std::string &filename) + { + { + ExodusII_IO exii(mesh); + exii.write(filename); + } + + TestCommWorld->barrier(); + + Mesh header_mesh(*TestCommWorld); + ExodusII_IO exii(header_mesh); + return exii.read_header(filename); + } + + void build_c0polyhedron(Mesh &mesh, + const std::vector &points, + const std::vector> &nodes_on_side) + { + for (auto p : index_range(points)) + mesh.add_point(points[p], /*id=*/p); + + std::vector> sides(nodes_on_side.size()); + for (auto s : index_range(nodes_on_side)) + { + const auto &nodes_on_s = nodes_on_side[s]; + sides[s] = std::make_shared(nodes_on_s.size()); + for (auto i : index_range(nodes_on_s)) + sides[s]->set_node(i, mesh.node_ptr(nodes_on_s[i])); + } + + std::unique_ptr mid_elem_node; + std::unique_ptr polyhedron = + std::make_unique(sides, mid_elem_node); + if (mid_elem_node) + mesh.add_node(std::move(mid_elem_node)); + + polyhedron->set_id() = 0; + Elem *elem = mesh.add_elem(std::move(polyhedron)); + elem->subdomain_id() = 1; + mesh.cache_elem_data(); + mesh.prepare_for_use(); + } + + void test_write_cube_header() + { + LOG_UNIT_TEST; + + Mesh mesh(*TestCommWorld); + MeshTools::Generation::build_cube(mesh, 2, 2, 2, + -1., 1., + -1., 1., + -1., 1., + C0POLYHEDRON); + mesh.get_boundary_info().clear(); + + for (auto &elem : mesh.element_ptr_range()) + elem->subdomain_id() = 1; + mesh.cache_elem_data(); + + ExodusHeaderInfo header_info = + this->write_and_read_header(mesh, "write_exodus_C0POLYHEDRON.e"); + + CPPUNIT_ASSERT_EQUAL(header_info.num_dim, 3); + CPPUNIT_ASSERT_EQUAL(header_info.num_elem, 8); + CPPUNIT_ASSERT_EQUAL(header_info.num_elem_blk, 1); + CPPUNIT_ASSERT_EQUAL(header_info.num_face, 48); + CPPUNIT_ASSERT_EQUAL(header_info.num_face_blk, 1); + CPPUNIT_ASSERT_EQUAL(header_info.num_node_sets, 0); + CPPUNIT_ASSERT_EQUAL(header_info.num_side_sets, 0); + } + + void test_write_hexagonal_prism_header() + { + LOG_UNIT_TEST; + + Mesh mesh(*TestCommWorld); + const std::vector points = + { { 0, -2, 0}, {-1, -1, 0}, {-1, 1, 0}, + { 0, 2, 0}, { 1, 1, 0}, { 1, -1, 0}, + { 0, -2, 1}, {-1, -1, 1}, {-1, 1, 1}, + { 0, 2, 1}, { 1, 1, 1}, { 1, -1, 1} }; + + const std::vector> nodes_on_side = + { {0, 1, 2, 3, 4, 5}, + {0, 1, 7, 6}, + {1, 2, 8, 7}, + {2, 3, 9, 8}, + {3, 4, 10, 9}, + {4, 5, 11, 10}, + {5, 0, 6, 11}, + {6, 7, 8, 9, 10, 11} }; + + this->build_c0polyhedron(mesh, points, nodes_on_side); + + ExodusHeaderInfo header_info = + this->write_and_read_header(mesh, "write_exodus_C0POLYHEDRON_HEXPRISM.e"); + + CPPUNIT_ASSERT_EQUAL(header_info.num_dim, 3); + CPPUNIT_ASSERT_EQUAL(header_info.num_elem, 1); + CPPUNIT_ASSERT_EQUAL(header_info.num_elem_blk, 1); + CPPUNIT_ASSERT_EQUAL(header_info.num_face, 8); + CPPUNIT_ASSERT_EQUAL(header_info.num_face_blk, 1); + CPPUNIT_ASSERT_EQUAL(header_info.num_node_sets, 0); + CPPUNIT_ASSERT_EQUAL(header_info.num_side_sets, 0); + } + + void test_write_and_read_hexagonal_prism() + { + LOG_UNIT_TEST; + + Mesh mesh(*TestCommWorld); + const std::vector points = + { { 0, -2, 0}, {-1, -1, 0}, {-1, 1, 0}, + { 0, 2, 0}, { 1, 1, 0}, { 1, -1, 0}, + { 0, -2, 1}, {-1, -1, 1}, {-1, 1, 1}, + { 0, 2, 1}, { 1, 1, 1}, { 1, -1, 1} }; + + const std::vector> nodes_on_side = + { {0, 1, 2, 3, 4, 5}, + {0, 1, 7, 6}, + {1, 2, 8, 7}, + {2, 3, 9, 8}, + {3, 4, 10, 9}, + {4, 5, 11, 10}, + {5, 0, 6, 11}, + {6, 7, 8, 9, 10, 11} }; + + this->build_c0polyhedron(mesh, points, nodes_on_side); + + { + ExodusII_IO exii(mesh); + exii.write("write_exodus_C0POLYHEDRON_HEXPRISM_READ.e"); + } + + Mesh input_mesh(*TestCommWorld); + ExodusII_IO exii_input(input_mesh); + if (input_mesh.processor_id() == 0) + exii_input.read("write_exodus_C0POLYHEDRON_HEXPRISM_READ.e"); + + MeshCommunication().broadcast(input_mesh); + input_mesh.prepare_for_use(); + + CPPUNIT_ASSERT_EQUAL(cast_int(1), input_mesh.n_elem()); + + const Elem *elem = input_mesh.elem_ptr(0); + CPPUNIT_ASSERT(elem); + CPPUNIT_ASSERT_EQUAL(C0POLYHEDRON, elem->type()); + CPPUNIT_ASSERT_EQUAL(12u, elem->n_vertices()); + CPPUNIT_ASSERT_EQUAL(8u, elem->n_sides()); + + for (auto s : index_range(nodes_on_side)) + { + const auto side_nodes = elem->nodes_on_side(s); + CPPUNIT_ASSERT_EQUAL(nodes_on_side[s].size(), side_nodes.size()); + for (auto n : index_range(nodes_on_side[s])) + CPPUNIT_ASSERT_EQUAL(cast_int(nodes_on_side[s][n]), + elem->node_id(side_nodes[n])); + } + } +}; + +CPPUNIT_TEST_SUITE_REGISTRATION(ExodusC0PolyhedronTest); + INSTANTIATE_EXODUSTEST(TET4); INSTANTIATE_EXODUSTEST(TET10); INSTANTIATE_EXODUSTEST(TET14); diff --git a/tests/mesh/mesh_generation_test.C b/tests/mesh/mesh_generation_test.C index 56c232f755..476c18192a 100644 --- a/tests/mesh/mesh_generation_test.C +++ b/tests/mesh/mesh_generation_test.C @@ -51,6 +51,7 @@ public: CPPUNIT_TEST( buildCubeHex8 ); CPPUNIT_TEST( buildCubeHex20 ); CPPUNIT_TEST( buildCubeHex27 ); + CPPUNIT_TEST( buildCubeC0Polyhedron ); CPPUNIT_TEST( buildCubePrism6 ); CPPUNIT_TEST( buildCubePrism15 ); CPPUNIT_TEST( buildCubePrism18 ); @@ -152,6 +153,23 @@ public: CPPUNIT_ASSERT(bbox.min()(1) <= Real(-4.0)); CPPUNIT_ASSERT(bbox.max()(1) >= Real(5.0)); + if (type == C0POLYGON) + { + BoundingBox nodal_bbox = MeshTools::create_nodal_bounding_box(mesh); + LIBMESH_ASSERT_FP_EQUAL(nodal_bbox.min()(0), + Real(-2.0), + TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(nodal_bbox.max()(0), + Real(3.0), + TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(nodal_bbox.min()(1), + Real(-4.0), + TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(nodal_bbox.max()(1), + Real(5.0), + TOLERANCE*TOLERANCE); + } + // Do serial assertions *after* all parallel assertions, so we // stay in sync after failure on only some processor(s) if (type != C0POLYGON) @@ -162,81 +180,92 @@ public: void testBuildCube(UnstructuredMesh & mesh, unsigned int n, ElemType type) { MeshTools::Generation::build_cube (mesh, n, n, n, -2.0, 3.0, -4.0, 5.0, -6.0, 7.0, type); - switch (Elem::type_to_n_sides_map[type]) + if (type == C0POLYHEDRON) { - case 4: // tets - CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), cast_int(n*n*n*24)); - break; - case 5: // prisms, pyramids - if (type == PRISM6 || type == PRISM15 || type == PRISM18 || - type == PRISM20 || type == PRISM21) - CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), cast_int(n*n*n*2)); - else - CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), cast_int(n*n*n*6)); - break; - case 6: // hexes + const dof_id_type grid_nodes = cast_int((n+1)*(n+1)*(n+1)); + CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), cast_int(n*n*n)); - break; - default: - libmesh_error(); + CPPUNIT_ASSERT(mesh.n_nodes() >= grid_nodes); + CPPUNIT_ASSERT(mesh.n_nodes() <= grid_nodes + cast_int(n*n*n)); } + else + { + switch (Elem::type_to_n_sides_map[type]) + { + case 4: // tets + CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), cast_int(n*n*n*24)); + break; + case 5: // prisms, pyramids + if (type == PRISM6 || type == PRISM15 || type == PRISM18 || + type == PRISM20 || type == PRISM21) + CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), cast_int(n*n*n*2)); + else + CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), cast_int(n*n*n*6)); + break; + case 6: // hexes + CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), cast_int(n*n*n)); + break; + default: + libmesh_error(); + } - switch (Elem::type_to_n_nodes_map[type]) - { - case 4: // First-order tets - CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), - cast_int((n+1)*(n+1)*(n+1) + n*n*n + 3*(n+1)*n*n)); - break; - case 6: // First-order prisms and hexes use the same nodes - case 8: - CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), - cast_int((n+1)*(n+1)*(n+1))); - break; - case 10: // Second-order tets - CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), - cast_int((2*n+1)*(2*n+1)*(2*n+1) + 14*n*n*n + 4*3*(n+1)*n*n)); - break; - case 18: - case 27: // Second-order prisms and hexes use the same nodes - CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), - cast_int((2*n+1)*(2*n+1)*(2*n+1))); - break; - case 20: - if (type == HEX20) - CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), - cast_int((2*n+1)*(2*n+1)*(2*n+1) - n*n*n - 3*(n+1)*n*n)); - if (type == PRISM20) - CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), - cast_int((2*n+1)*(2*n+1)*(2*n+1) + 2*(n+1)*n*n)); - break; - case 21: // Prisms based on full Tri7 cross sections - CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), - cast_int((2*n+1)*(2*n+1)*(2*n+1) + 2*(2*n+1)*n*n)); - break; - case 15: // weird partial order prism - CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), - cast_int((2*n+1)*(2*n+1)*(2*n+1) - n*n*n - 2*(n+1)*n*n)); - break; - case 5: // pyramids - CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), - cast_int((n+1)*(n+1)*(n+1) + n*n*n)); - break; - case 13: - CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), - cast_int((2*n+1)*(2*n+1)*(2*n+1) + 8*n*n*n - 3*(n+1)*n*n)); - break; - case 14: // pyramids, tets - if (type == PYRAMID14) - CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), - cast_int((2*n+1)*(2*n+1)*(2*n+1) + 8*n*n*n)); - else // TET14 - CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), - cast_int((2*n+1)*(2*n+1)*(2*n+1) + 14*n*n*n + 4*3*(n+1)*n*n + - 36*n*n*n + 4*3*(n+1)*n*n)); - break; - default: - libmesh_error(); + switch (Elem::type_to_n_nodes_map[type]) + { + case 4: // First-order tets + CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), + cast_int((n+1)*(n+1)*(n+1) + n*n*n + 3*(n+1)*n*n)); + break; + case 6: // First-order prisms and hexes use the same nodes + case 8: + CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), + cast_int((n+1)*(n+1)*(n+1))); + break; + case 10: // Second-order tets + CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), + cast_int((2*n+1)*(2*n+1)*(2*n+1) + 14*n*n*n + 4*3*(n+1)*n*n)); + break; + case 18: + case 27: // Second-order prisms and hexes use the same nodes + CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), + cast_int((2*n+1)*(2*n+1)*(2*n+1))); + break; + case 20: + if (type == HEX20) + CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), + cast_int((2*n+1)*(2*n+1)*(2*n+1) - n*n*n - 3*(n+1)*n*n)); + if (type == PRISM20) + CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), + cast_int((2*n+1)*(2*n+1)*(2*n+1) + 2*(n+1)*n*n)); + break; + case 21: // Prisms based on full Tri7 cross sections + CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), + cast_int((2*n+1)*(2*n+1)*(2*n+1) + 2*(2*n+1)*n*n)); + break; + case 15: // weird partial order prism + CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), + cast_int((2*n+1)*(2*n+1)*(2*n+1) - n*n*n - 2*(n+1)*n*n)); + break; + case 5: // pyramids + CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), + cast_int((n+1)*(n+1)*(n+1) + n*n*n)); + break; + case 13: + CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), + cast_int((2*n+1)*(2*n+1)*(2*n+1) + 8*n*n*n - 3*(n+1)*n*n)); + break; + case 14: // pyramids, tets + if (type == PYRAMID14) + CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), + cast_int((2*n+1)*(2*n+1)*(2*n+1) + 8*n*n*n)); + else // TET14 + CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), + cast_int((2*n+1)*(2*n+1)*(2*n+1) + 14*n*n*n + 4*3*(n+1)*n*n + + 36*n*n*n + 4*3*(n+1)*n*n)); + break; + default: + libmesh_error(); + } } // Our bounding boxes can be loose on higher order elements, but @@ -249,6 +278,44 @@ public: CPPUNIT_ASSERT(bbox.min()(2) <= Real(-6.0)); CPPUNIT_ASSERT(bbox.max()(2) >= Real(7.0)); + if (type == C0POLYHEDRON) + { + BoundingBox nodal_bbox = MeshTools::create_nodal_bounding_box(mesh); + const Real expected_volume = Real(5) * Real(9) * Real(13); + + LIBMESH_ASSERT_FP_EQUAL(nodal_bbox.min()(0), + Real(-2.0), + TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(nodal_bbox.max()(0), + Real(3.0), + TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(nodal_bbox.min()(1), + Real(-4.0), + TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(nodal_bbox.max()(1), + Real(5.0), + TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(nodal_bbox.min()(2), + Real(-6.0), + TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(nodal_bbox.max()(2), + Real(7.0), + TOLERANCE*TOLERANCE); + LIBMESH_ASSERT_FP_EQUAL(MeshTools::volume(mesh), + expected_volume, + TOLERANCE*TOLERANCE); + + for (auto & elem : mesh.element_ptr_range()) + { + CPPUNIT_ASSERT_EQUAL(elem->type(), C0POLYHEDRON); + CPPUNIT_ASSERT_EQUAL(elem->n_sides(), 6u); + CPPUNIT_ASSERT(elem->n_sub_elem() > 0); + CPPUNIT_ASSERT(elem->volume() > 0); + } + + return; + } + // We don't yet try to do affine map optimizations on pyramids if (type == PYRAMID5 || type == PYRAMID13 || @@ -314,6 +381,8 @@ public: void buildCubeHex8 () { LOG_UNIT_TEST; tester(&MeshGenerationTest::testBuildCube, 2, HEX8); } void buildCubeHex20 () { LOG_UNIT_TEST; tester(&MeshGenerationTest::testBuildCube, 2, HEX20); } void buildCubeHex27 () { LOG_UNIT_TEST; tester(&MeshGenerationTest::testBuildCube, 2, HEX27); } + void buildCubeC0Polyhedron () + { LOG_UNIT_TEST; tester(&MeshGenerationTest::testBuildCube, 2, C0POLYHEDRON); } void buildCubePrism6 () { LOG_UNIT_TEST; tester(&MeshGenerationTest::testBuildCube, 2, PRISM6); } void buildCubePrism15 () { LOG_UNIT_TEST; tester(&MeshGenerationTest::testBuildCube, 2, PRISM15); } void buildCubePrism18 () { LOG_UNIT_TEST; tester(&MeshGenerationTest::testBuildCube, 2, PRISM18); }