[MOAB-dev] r2874 - in MOAB/trunk: . test/h5file
kraftche at cae.wisc.edu
kraftche at cae.wisc.edu
Mon May 4 12:14:30 CDT 2009
Author: kraftche
Date: 2009-05-04 12:14:29 -0500 (Mon, 04 May 2009)
New Revision: 2874
Modified:
MOAB/trunk/README.IO
MOAB/trunk/ReadHDF5.cpp
MOAB/trunk/ReadHDF5.hpp
MOAB/trunk/test/h5file/h5partial.cpp
Log:
during partial read of H5M file, read any explicitly defined element sides
Modified: MOAB/trunk/README.IO
===================================================================
--- MOAB/trunk/README.IO 2009-05-01 22:15:22 UTC (rev 2873)
+++ MOAB/trunk/README.IO 2009-05-04 17:14:29 UTC (rev 2874)
@@ -134,13 +134,15 @@
MOAB Native (HDF5-based MHDF) format
------------------------------------
- ELEMENTS={EXPLICIT|NODES}
+ ELEMENTS={EXPLICIT|NODES|SIDES}
If reading only part of a file, specify which additional elements to read
in addition to those contained in the specified set. The options are:
EXPLICIT - read only explicitly designated elements
NODES - read any element for which all the nodes are being read.
-Default is EXPLICIT unless only nodesa are explicitly specified, in which
+ SIDES - read explicilty specified elements and any elements that
+ are sides of those elements.
+Default is SIDES unless only nodes are explicitly specified, in which
case NODES will be used.
CHILDREN={NONE|SETS|CONTENTS}
Modified: MOAB/trunk/ReadHDF5.cpp
===================================================================
--- MOAB/trunk/ReadHDF5.cpp 2009-05-01 22:15:22 UTC (rev 2873)
+++ MOAB/trunk/ReadHDF5.cpp 2009-05-04 17:14:29 UTC (rev 2874)
@@ -472,59 +472,63 @@
DEBUGOUT( "READING ELEMENTS" );
// decide if we need to read additional elements
- int mode;
- const char* const options[] = { "EXPLICIT", "ADJACENT", 0 };
- rval = opts.match_option( "ELEMENTS", options, mode );
- bool read_adjacent_elements;
- if (MB_SUCCESS == rval) {
- read_adjacent_elements = (mode == 1);
- }
+ int side_mode;
+ const char* const options[] = { "EXPLICIT", "NODES", "SIDES", 0 };
+ rval = opts.match_option( "ELEMENTS", options, side_mode );
if (MB_ENTITY_NOT_FOUND == rval) {
// Chose default based on whether or not any elements have been
// specified.
MBRange tmp;
intersect( fileInfo->nodes, file_ids, tmp );
- if (tmp.empty()) // can't do 'ADJACENT' mode if no nodes
- read_adjacent_elements = false;
- else { // check if we have any elements, default to 'ADJACENT' only if no
- // explicitly specified elements
- tmp.merge( sets );
- MBRange elems = file_ids.subtract( tmp );
- read_adjacent_elements = elems.empty();
- }
+ // If only nodes were specified, then default to "NODES", otherwise
+ // default to "SIDES".
+ if (!tmp.empty() && (tmp.size() + sets.size()) == file_ids.size())
+ side_mode = 1;
+ else
+ side_mode = 2;
}
- else {
+ else if (MB_SUCCESS != rval) {
readUtil->report_error( "Invalid value for 'ELEMENTS' option" );
return error(rval);
}
- if (read_adjacent_elements) {
+ switch (side_mode) {
+ case 0: // ELEMENTS=EXPLICIT : read only specified element IDS
for (int dim = 1; dim <= 3; ++dim) {
for (int i = 0; i < fileInfo->num_elem_desc; ++i) {
MBEntityType type = MBCN::EntityTypeFromName( fileInfo->elems[i].type );
if (MBCN::Dimension(type) == dim) {
- rval = read_node_adj_elems( fileInfo->elems[i] );
+ MBRange subset;
+ intersect( fileInfo->elems[i].desc, file_ids, subset );
+ rval = read_elems( fileInfo->elems[i], subset );
if (MB_SUCCESS != rval)
return error(rval);
}
}
}
- }
- else {
+ break;
+
+ case 1: // ELEMENTS=NODES : read all elements for which all nodes have been read
for (int dim = 1; dim <= 3; ++dim) {
for (int i = 0; i < fileInfo->num_elem_desc; ++i) {
MBEntityType type = MBCN::EntityTypeFromName( fileInfo->elems[i].type );
if (MBCN::Dimension(type) == dim) {
- MBRange subset;
- intersect( fileInfo->elems[i].desc, file_ids, subset );
- rval = read_elems( fileInfo->elems[i], subset );
+ rval = read_node_adj_elems( fileInfo->elems[i] );
if (MB_SUCCESS != rval)
return error(rval);
}
}
}
+ break;
+
+ case 2: // ELEMENTS=SIDES : read explicitly specified elems and any sides of those elems
+ rval = read_elements_and_sides( file_ids );
+ if (MB_SUCCESS != rval)
+ return error(rval);
+ break;
}
+
DEBUGOUT( "READING SETS" );
// If reading contained/child sets but not their contents then find
@@ -1155,6 +1159,163 @@
}
*/
+MBErrorCode ReadHDF5::read_elements_and_sides( const MBRange& file_ids )
+{
+ MBErrorCode rval;
+ MBEntityType type;
+
+ // determine largest element dimension
+ int max_dim = 0;
+ for (int i = 0; i < fileInfo->num_elem_desc; ++i) {
+ type = MBCN::EntityTypeFromName( fileInfo->elems[i].type );
+ int dim = MBCN::Dimension(type);
+ if (dim > max_dim) {
+ MBEntityHandle start = (MBEntityHandle)fileInfo->elems[i].desc.start_id;
+ MBRange::iterator it = file_ids.lower_bound( start );
+ if (it != file_ids.end() && (long)(*it - start ) < fileInfo->elems[i].desc.count)
+ max_dim = dim;
+ }
+ }
+
+ // Get MBRange of element IDs only
+ MBRange elem_ids( file_ids );
+ MBRange::iterator s, e;
+ if (fileInfo->nodes.count) {
+ MBEntityHandle first = (MBEntityHandle)fileInfo->nodes.start_id;
+ s = elem_ids.lower_bound( first );
+ e = MBRange::lower_bound( s, elem_ids.end(), first + fileInfo->nodes.count );
+ elem_ids.erase( s, e );
+ }
+ if (fileInfo->sets.count) {
+ MBEntityHandle first = (MBEntityHandle)fileInfo->sets.start_id;
+ s = elem_ids.lower_bound( first );
+ e = MBRange::lower_bound( s, elem_ids.end(), first + fileInfo->sets.count );
+ elem_ids.erase( s, e );
+ }
+
+ // read all node-adjacent elements of smaller dimensions
+ for (int i = 0; i < fileInfo->num_elem_desc; ++i) {
+ MBEntityType type = MBCN::EntityTypeFromName( fileInfo->elems[i].type );
+ if (MBCN::Dimension(type) < max_dim) {
+ rval = read_node_adj_elems( fileInfo->elems[i] );
+ if (MB_SUCCESS != rval)
+ return error(rval);
+ }
+ }
+
+ // Read the subset of the explicitly sepecified elements that are of the
+ // largest dimension. We could read all explicitly specified elements,
+ // but then the lower-dimension elements will be read a second when reading
+ // node-adjacent elements.
+ for (int i = 0; i < fileInfo->num_elem_desc; ++i) {
+ type = MBCN::EntityTypeFromName( fileInfo->elems[i].type );
+ if (MBCN::Dimension(type) == max_dim) {
+ MBRange subset;
+ intersect( fileInfo->elems[i].desc, elem_ids, subset );
+ if (!subset.empty()) {
+ elem_ids.subtract( subset );
+ rval = read_elems( fileInfo->elems[i], subset );
+ if (MB_SUCCESS != rval)
+ return error(rval);
+ }
+ }
+ }
+
+ // delete anything we read in that we should not have
+ // (e.g. an edge spanning two disjoint blocks of elements)
+
+ // get MBRange of all explicitly specified elements, grouped by dimension
+ MBRange explicit_elems[4];
+ MBRange::iterator hints[4] = { explicit_elems[0].begin(),
+ explicit_elems[1].begin(),
+ explicit_elems[2].begin(),
+ explicit_elems[3].begin() };
+ RangeMap<long, MBEntityHandle>::iterator rit;
+ while (!elem_ids.empty()) {
+ long start = elem_ids.front();
+ long count = elem_ids.const_pair_begin()->second - start + 1;
+ rit = idMap.lower_bound( start );
+ assert( rit->begin <= start && rit->begin + rit->count > start );
+ long offset = start - rit->begin;
+ long avail = rit->count - offset;
+ if (avail > count)
+ count = avail;
+ MBEntityHandle start_h = rit->value + offset;
+ int d = MBCN::Dimension( TYPE_FROM_HANDLE( start_h ) );
+ if (MBCN::Dimension( TYPE_FROM_HANDLE( start_h + count - 1 ) ) != d)
+ count = start - LAST_HANDLE( TYPE_FROM_HANDLE(start_h) ) + 1;
+
+ elem_ids.erase( elem_ids.begin(), elem_ids.begin() + count );
+ hints[d] = explicit_elems[d].insert( hints[d], rit->value + offset, rit->value + offset + count - 1 );
+ }
+
+ // get MBRange of all read elements, and remove any that were explicity specified
+ // get handles for everything we've read
+ MBRange all_elems;
+ MBRange::iterator hint = all_elems.begin();
+ for (rit = idMap.begin(); rit != idMap.end(); ++rit)
+ hint = all_elems.insert( hint, rit->value, rit->value + rit->count - 1 );
+ // remove any vertex handles
+ all_elems.erase( all_elems.begin(), all_elems.upper_bound( MBVERTEX ) );
+ // remove any set handles and any handles for elements >= max_dim
+ all_elems.erase( all_elems.lower_bound( MBCN::TypeDimensionMap[max_dim].first ), all_elems.end() );
+ // remove explicit elements < max_dim
+ if (!explicit_elems[1].empty() && max_dim > 1)
+ all_elems = all_elems.subtract( explicit_elems[1] );
+ if (!explicit_elems[2].empty() && max_dim > 2)
+ all_elems = all_elems.subtract( explicit_elems[2] );
+
+ // remove any elements that are adjacent to some explicity specified element.
+ if (max_dim > 1 && !explicit_elems[2].empty()) {
+ MBRange adj;
+ rval = iFace->get_adjacencies( explicit_elems[2], 1, false, adj, MBInterface::UNION );
+ if (MB_SUCCESS != rval)
+ return error(rval);
+ if (!adj.empty())
+ all_elems = all_elems.subtract( adj );
+ }
+ if (max_dim == 3) {
+ MBRange adj;
+ rval = iFace->get_adjacencies( explicit_elems[3], 1, false, adj, MBInterface::UNION );
+ if (MB_SUCCESS != rval)
+ return error(rval);
+ if (!adj.empty())
+ all_elems = all_elems.subtract( adj );
+ adj.clear();
+ rval = iFace->get_adjacencies( explicit_elems[3], 2, false, adj, MBInterface::UNION );
+ if (MB_SUCCESS != rval)
+ return error(rval);
+ if (!adj.empty())
+ all_elems = all_elems.subtract( adj );
+ }
+
+ // now delete anything remaining in all_elems
+ rval = iFace->delete_entities( all_elems );
+ if (MB_SUCCESS != rval)
+ return error(rval);
+
+ // remove dead entities from ID map
+ while (!all_elems.empty()) {
+ MBEntityHandle start = all_elems.front();
+ MBEntityID count = all_elems.const_pair_begin()->second - start + 1;
+ for (rit = idMap.begin(); rit != idMap.end(); ++rit)
+ if (rit->value <= start && (long)(start - rit->value) < rit->count)
+ break;
+ if (rit == idMap.end())
+ return error(MB_FAILURE);
+
+ MBEntityID offset = start - rit->value;
+ MBEntityID avail = rit->count - offset;
+ if (avail < count)
+ count = avail;
+
+ all_elems.erase( all_elems.begin(), all_elems.begin() + count );
+ idMap.erase( rit->begin + offset, count );
+ }
+
+ return MB_SUCCESS;
+}
+
MBErrorCode ReadHDF5::read_sets( const MBRange& file_ids )
{
mhdf_Status status;
Modified: MOAB/trunk/ReadHDF5.hpp
===================================================================
--- MOAB/trunk/ReadHDF5.hpp 2009-05-01 22:15:22 UTC (rev 2873)
+++ MOAB/trunk/ReadHDF5.hpp 2009-05-04 17:14:29 UTC (rev 2874)
@@ -132,6 +132,9 @@
//! Read poly(gons|hedra)
MBErrorCode read_poly( const mhdf_ElemDesc& elems, const MBRange& file_ids );
+ //! Read specified elements and any adjacent elemnets of lower dimension
+ MBErrorCode read_elements_and_sides( const MBRange& file_ids );
+
//! Read sets
MBErrorCode read_sets( const MBRange& set_file_ids );
Modified: MOAB/trunk/test/h5file/h5partial.cpp
===================================================================
--- MOAB/trunk/test/h5file/h5partial.cpp 2009-05-01 22:15:22 UTC (rev 2873)
+++ MOAB/trunk/test/h5file/h5partial.cpp 2009-05-04 17:14:29 UTC (rev 2874)
@@ -132,6 +132,7 @@
void test_read_tagged_nodes();
+void test_read_sides();
int main( int argc, char* argv[] )
@@ -167,6 +168,7 @@
result += RUN_TEST(test_read_adjacencies);
result += RUN_TEST(test_read_tagged_elems);
result += RUN_TEST(test_read_tagged_nodes);
+ result += RUN_TEST(test_read_sides);
if (argc == 1)
remove( TEST_FILE );
@@ -1396,3 +1398,85 @@
}
+void test_read_sides()
+{
+ MBErrorCode rval;
+ MBCore instance;
+ MBInterface& mb = instance;
+
+ // create 4x4 grid of quads with edges
+ const int INT = 4;
+ MBEntityHandle verts[INT+1][INT+1];
+ for (int j = 0; j <= INT; ++j) {
+ for (int i = 0; i <= INT; ++i) {
+ double coords[3] = { i, j, 0 };
+ rval = mb.create_vertex( coords, verts[INT-j][i] );
+ CHECK_ERR(rval);
+ }
+ }
+ MBEntityHandle quads[INT][INT];
+ for (int j = 0; j < INT; ++j) {
+ for (int i = 0; i < INT; ++i) {
+ MBEntityHandle conn[4] = { verts[INT-j][i],
+ verts[INT-j][i+1],
+ verts[INT-j-1][i+1],
+ verts[INT-j-1][i] };
+ rval = mb.create_element( MBQUAD, conn, 4, quads[INT-j-1][i] );
+ CHECK_ERR(rval);
+ }
+ }
+ MBRange edges;
+ rval = mb.get_adjacencies( &quads[0][0], INT*INT, 1, true, edges, MBInterface::UNION );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 40, (int)edges.size() );
+
+ // group quads into two sets
+ MBEntityHandle sets[2];
+ rval = mb.create_meshset( MESHSET_SET, sets[0] ); CHECK_ERR(rval);
+ rval = mb.create_meshset( MESHSET_SET, sets[1] ); CHECK_ERR(rval);
+ rval = mb.add_entities( sets[0], quads[0], INT ); CHECK_ERR(rval);
+ rval = mb.add_entities( sets[1], quads[1], INT ); CHECK_ERR(rval);
+ rval = mb.add_entities( sets[0], quads[2], INT ); CHECK_ERR(rval);
+ rval = mb.add_entities( sets[1], quads[3], INT ); CHECK_ERR(rval);
+
+ // assign IDS
+ MBTag id_tag;
+ rval = mb.tag_create( ID_TAG_NAME, sizeof(int), MB_TAG_SPARSE, MB_TYPE_INTEGER, id_tag, 0 );
+ CHECK_ERR(rval);
+ int ids[2] = { 4, 5 };
+ rval = mb.tag_set_data( id_tag, sets, 2, ids );
+ CHECK_ERR(rval);
+
+ // write mesh
+ rval = mb.write_file( TEST_FILE, "MOAB" );
+ CHECK_ERR(rval);
+
+ // read first set back in
+ rval = mb.delete_mesh(); CHECK_ERR(rval);
+ MBEntityHandle file;
+ rval = mb.load_file( TEST_FILE, file, READ_OPTS, ID_TAG_NAME, ids, 1 );
+ CHECK_ERR(rval);
+
+ // check expected counts
+ int count;
+ rval = mb.get_number_entities_by_type( 0, MBVERTEX, count );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( (INT+1)*INT, count );
+ rval = mb.get_number_entities_by_type( 0, MBQUAD, count );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( INT*INT/2, count );
+ rval = mb.get_number_entities_by_type( 0, MBEDGE, count );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 2*(INT+1) + INT*INT, count );
+
+ // check edges adjacent to each quad
+ MBRange elems;
+ rval = mb.get_entities_by_type( 0, MBQUAD, elems );
+ CHECK_ERR(rval);
+ for (MBRange::iterator it = elems.begin(); it != elems.end(); ++it) {
+ edges.clear();
+ rval = mb.get_adjacencies( &*it, 1, 1, false, edges );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 4, (int)edges.size() );
+ }
+}
More information about the moab-dev
mailing list