[MOAB-dev] r3373 - in MOAB/trunk: . tools tools/mbperf
kraftche at cae.wisc.edu
kraftche at cae.wisc.edu
Thu Nov 19 17:10:26 CST 2009
Author: kraftche
Date: 2009-11-19 17:10:26 -0600 (Thu, 19 Nov 2009)
New Revision: 3373
Modified:
MOAB/trunk/MBSkinner.cpp
MOAB/trunk/MBSkinner.hpp
MOAB/trunk/MBTest.cpp
MOAB/trunk/tools/mbperf/mbperf.cpp
MOAB/trunk/tools/skin.cpp
Log:
o New, faster algorithm for skinning using vertex-to-element adjacencies
o add tests to mbperf for measuring skinning performance
o add new skinning features:
- creation of skin entities is optional: can get skin vertices or
just existing skin entities w/out creating anything new
- splitting of results into forward and reverse groups such that
existing skin faces that are reversed are placed in the second
is now optional. If no MBRange is provided for reversed entities,
then everything is returned in the forward set.
o add lots of new skinning tests:
- test that entities are not created when they aren't supposed to be
- test that skinning of polyhedra works
- test that skinning of higher-order elements works
- test that returning of reversed sides in a separate MBRange works
- tests are done for both 2D and 3D cases for both algorithms
Modified: MOAB/trunk/MBSkinner.cpp
===================================================================
--- MOAB/trunk/MBSkinner.cpp 2009-11-19 23:08:51 UTC (rev 3372)
+++ MOAB/trunk/MBSkinner.cpp 2009-11-19 23:10:26 UTC (rev 3373)
@@ -13,7 +13,7 @@
*
*/
-#ifdef WIN32
+#ifdef __MFC_VER
#pragma warning(disable:4786)
#endif
@@ -29,23 +29,17 @@
#include "MBUtil.hpp"
#include "MBInternals.hpp"
#include "MBTagConventions.hpp"
-
-#ifndef IS_BUILDING_MB
-#define IS_BUILDING_MB
-#define I_CREATED_BUILDING
-#endif
#include "MBCore.hpp"
#include "AEntityFactory.hpp"
-#ifdef I_CREATED_BUILDING
-#undef IS_BUILDING_MB
-#undef I_CREATED_BUILDING
+
+#ifdef M_PI
+# define SKINNER_PI M_PI
+#else
+# define SKINNER_PI 3.1415926535897932384626
#endif
-#define SKINNER_PI 3.1415926535897932384626
-
-
MBSkinner::~MBSkinner()
{
// delete the adjacency tag
@@ -236,10 +230,86 @@
return result;
}
-MBErrorCode MBSkinner::find_skin(const MBRange &source_entities,
+MBErrorCode MBSkinner::find_skin( const MBRange& source_entities,
+ bool get_vertices,
+ MBRange& output_handles,
+ MBRange* output_reverse_handles,
+ bool create_vert_elem_adjs,
+ bool create_skin_elements )
+{
+ if (source_entities.empty())
+ return MB_SUCCESS;
+
+ MBCore* this_core = dynamic_cast<MBCore*>(thisMB);
+ if (this_core && create_vert_elem_adjs &&
+ !this_core->a_entity_factory()->vert_elem_adjacencies())
+ this_core->a_entity_factory()->create_vert_elem_adjacencies();
+
+ if (this_core && this_core->a_entity_factory()->vert_elem_adjacencies())
+ return find_skin_vertices( source_entities,
+ get_vertices ? &output_handles : 0,
+ get_vertices ? 0 : &output_handles,
+ output_reverse_handles,
+ create_skin_elements );
+
+ MBRange forward, reverse;
+ MBRange prev;
+ const int d = MBCN::Dimension(TYPE_FROM_HANDLE(source_entities.front()));
+ if (!source_entities.all_of_dimension(d))
+ return MB_TYPE_OUT_OF_RANGE;
+
+ MBErrorCode rval = thisMB->get_entities_by_dimension( 0, d-1, prev );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ rval = find_skin_noadj( source_entities, forward, reverse );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ if (get_vertices && !output_reverse_handles) {
+ forward.merge( reverse );
+ reverse.clear();
+ }
+
+ if (get_vertices) {
+ rval = thisMB->get_connectivity( forward, output_handles );
+ if (MB_SUCCESS != rval)
+ return rval;
+ }
+
+ if (!create_skin_elements) {
+ MBRange new_skin;
+ rval = thisMB->get_entities_by_dimension( 0, d-1, new_skin);
+ if (MB_SUCCESS != rval)
+ return rval;
+ new_skin = subtract( new_skin, prev );
+ forward = subtract( forward, new_skin );
+ reverse = subtract( reverse, new_skin );
+ rval = thisMB->delete_entities( new_skin );
+ if (MB_SUCCESS != rval)
+ return rval;
+ }
+
+ if (!get_vertices) {
+ if (output_handles.empty())
+ output_handles.swap( forward );
+ else
+ output_handles.merge( forward );
+ if (output_reverse_handles) {
+ if (output_reverse_handles->empty())
+ output_reverse_handles->swap( reverse );
+ else
+ output_reverse_handles->merge( reverse );
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode MBSkinner::find_skin_noadj(const MBRange &source_entities,
MBRange &forward_target_entities,
- MBRange &reverse_target_entities,
- bool create_vert_elem_adjs)
+ MBRange &reverse_target_entities/*,
+ bool create_vert_elem_adjs*/)
{
if(source_entities.empty())
return MB_FAILURE;
@@ -253,17 +323,18 @@
if(mTargetDim < 0 || source_dim > 3)
return MB_FAILURE;
- MBCore *this_core = dynamic_cast<MBCore*>(thisMB);
- bool use_adjs = false;
- if (!this_core->a_entity_factory()->vert_elem_adjacencies() &&
- create_vert_elem_adjs)
- this_core->a_entity_factory()->create_vert_elem_adjacencies();
+ //MBCore *this_core = dynamic_cast<MBCore*>(thisMB);
+ //bool use_adjs = false;
+ //if (!this_core->a_entity_factory()->vert_elem_adjacencies() &&
+ // create_vert_elem_adjs)
+ // this_core->a_entity_factory()->create_vert_elem_adjacencies();
+ //
+ //if (this_core->a_entity_factory()->vert_elem_adjacencies())
+ // use_adjs = true;
+ //
+ //else
+ initialize();
- if (this_core->a_entity_factory()->vert_elem_adjacencies())
- use_adjs = true;
-
- else initialize();
-
MBRange::const_iterator iter, end_iter;
end_iter = source_entities.end();
const MBEntityHandle *conn;
@@ -300,62 +371,62 @@
for(int j=0; j<num_sub_nodes; j++)
sub_conn[j] = conn[sub_indices[j]];
- if (use_adjs) {
- dum_elems.clear();
- result = thisMB->get_adjacencies(sub_conn, num_sub_nodes, source_dim, false,
- dum_elems);
- if (MB_SUCCESS != result) return result;
- dum_elems = intersect( dum_elems, source_entities );
- if (dum_elems.empty()) {
- assert(false); // should never happen
- return MB_FAILURE;
- }
- // if (dum_elems.size() > 2) {
- // source entities do not form a valid source-dim patch (t-junction).
- // do we care?
- // }
+// if (use_adjs) {
+// dum_elems.clear();
+// result = thisMB->get_adjacencies(sub_conn, num_sub_nodes, source_dim, false,
+// dum_elems);
+// if (MB_SUCCESS != result) return result;
+// dum_elems = intersect( dum_elems, source_entities );
+// if (dum_elems.empty()) {
+// assert(false); // should never happen
+// return MB_FAILURE;
+// }
+// // if (dum_elems.size() > 2) {
+// // source entities do not form a valid source-dim patch (t-junction).
+// // do we care?
+// // }
+//
+// if (1 == dum_elems.size()) {
+// // this sub_element is on the skin
+//
+// // check for existing entity
+// dum_sub_elems.clear();
+// result = thisMB->get_adjacencies(sub_conn, num_sub_nodes, mTargetDim, false,
+// dum_sub_elems);
+// if (MB_SUCCESS != result) return result;
+// if (dum_sub_elems.empty()) {
+// // need to create one
+// MBEntityHandle tmphndl=0;
+// int indices[MB_MAX_SUB_ENTITY_VERTICES];
+// MBEntityType new_type;
+// int num_new_nodes;
+// MBCN::SubEntityNodeIndices( type, num_nodes, mTargetDim, i, new_type, num_new_nodes, indices );
+// for(int j=0; j<num_new_nodes; j++)
+// sub_conn[j] = conn[indices[j]];
+//
+// result = thisMB->create_element(new_type, sub_conn,
+// num_new_nodes, tmphndl);
+// forward_target_entities.insert(tmphndl);
+// }
+// else {
+// // else find the relative sense of this entity to the source_entity in this set
+// int side_no, sense = 0, offset;
+// if (source_entities.find(*dum_elems.begin()) == source_entities.end()) {
+// result = thisMB->side_number(*dum_elems.rbegin(), *dum_sub_elems.begin(),
+// side_no, sense, offset);
+// }
+// else {
+// result = thisMB->side_number(*dum_elems.begin(), *dum_sub_elems.begin(),
+// side_no, sense, offset);
+// }
+// if (-1 == sense) reverse_target_entities.insert(*dum_sub_elems.begin());
+// else if (1 == sense) forward_target_entities.insert(*dum_sub_elems.begin());
+// else return MB_FAILURE;
+// }
+// }
+// }
+// else {
- if (1 == dum_elems.size()) {
- // this sub_element is on the skin
-
- // check for existing entity
- dum_sub_elems.clear();
- result = thisMB->get_adjacencies(sub_conn, num_sub_nodes, mTargetDim, false,
- dum_sub_elems);
- if (MB_SUCCESS != result) return result;
- if (dum_sub_elems.empty()) {
- // need to create one
- MBEntityHandle tmphndl=0;
- int indices[MB_MAX_SUB_ENTITY_VERTICES];
- MBEntityType new_type;
- int num_new_nodes;
- MBCN::SubEntityNodeIndices( type, num_nodes, mTargetDim, i, new_type, num_new_nodes, indices );
- for(int j=0; j<num_new_nodes; j++)
- sub_conn[j] = conn[indices[j]];
-
- result = thisMB->create_element(new_type, sub_conn,
- num_new_nodes, tmphndl);
- forward_target_entities.insert(tmphndl);
- }
- else {
- // else find the relative sense of this entity to the source_entity in this set
- int side_no, sense = 0, offset;
- if (source_entities.find(*dum_elems.begin()) == source_entities.end()) {
- result = thisMB->side_number(*dum_elems.rbegin(), *dum_sub_elems.begin(),
- side_no, sense, offset);
- }
- else {
- result = thisMB->side_number(*dum_elems.begin(), *dum_sub_elems.begin(),
- side_no, sense, offset);
- }
- if (-1 == sense) reverse_target_entities.insert(*dum_sub_elems.begin());
- else if (1 == sense) forward_target_entities.insert(*dum_sub_elems.begin());
- else return MB_FAILURE;
- }
- }
- }
- else {
-
// see if we can match this connectivity with
// an existing entity
find_match( sub_type, sub_conn, num_sub_nodes, match, direct );
@@ -384,7 +455,7 @@
{
forward_target_entities.erase(seek_iter);
remove_adjacency(match);
- if(!use_adjs && entity_deletable(match))
+ if(/*!use_adjs &&*/ entity_deletable(match))
{
result = thisMB->delete_entities(&match, 1);
assert(MB_SUCCESS == result);
@@ -394,7 +465,7 @@
{
reverse_target_entities.erase(seek_iter);
remove_adjacency(match);
- if(!use_adjs && entity_deletable(match))
+ if(/*!use_adjs &&*/ entity_deletable(match))
{
result = thisMB->delete_entities(&match, 1);
assert(MB_SUCCESS == result);
@@ -412,7 +483,7 @@
}
}
}
- }
+ //}
}
}
@@ -912,23 +983,786 @@
MBRange &skin_entities,
bool create_vert_elem_adjs)
{
- if (MBCN::Dimension(TYPE_FROM_HANDLE(*entities.begin())) !=
- MBCN::Dimension(TYPE_FROM_HANDLE(*entities.rbegin())))
- return MB_FAILURE;
+ MBRange tmp_skin;
+ MBErrorCode result = find_skin(entities, (dim==0), tmp_skin, 0,
+ create_vert_elem_adjs, true);
+ if (MB_SUCCESS != result || tmp_skin.empty()) return result;
- MBRange tmp_skin_for, tmp_skin_rev;
- MBErrorCode result = find_skin(entities, tmp_skin_for, tmp_skin_rev,
- create_vert_elem_adjs);
- if (MB_SUCCESS != result) return result;
+ if (tmp_skin.all_of_dimension(dim)) {
+ if (skin_entities.empty())
+ skin_entities.swap(tmp_skin);
+ else
+ skin_entities.merge(tmp_skin);
+ }
+ else {
+ result = thisMB->get_adjacencies( tmp_skin, dim, true, skin_entities,
+ MBInterface::UNION );
+ }
- tmp_skin_for.merge(tmp_skin_rev);
- if (MBCN::Dimension(TYPE_FROM_HANDLE(*tmp_skin_for.begin())) != dim ||
- MBCN::Dimension(TYPE_FROM_HANDLE(*tmp_skin_for.rbegin())) != dim)
- result = thisMB->get_adjacencies(tmp_skin_for, dim, true, skin_entities,
- MBInterface::UNION);
- else
- skin_entities.swap(tmp_skin_for);
-
return result;
}
+// Bit tag uses less memory, but is slower. Probably
+// because we have to query one entity at a time.
+// -- j.kraftcheck 2009-11-10
+const bool use_bit_tag = false;
+
+MBErrorCode MBSkinner::find_skin_vertices( const MBRange& entities,
+ MBRange* skin_verts,
+ MBRange* skin_elems,
+ MBRange* skin_rev_elems,
+ bool create_skin_elems,
+ bool corners_only )
+{
+ MBErrorCode rval;
+ if (entities.empty())
+ return MB_SUCCESS;
+
+ const int dim = MBCN::Dimension(TYPE_FROM_HANDLE(entities.front()));
+ if (dim < 1 || dim > 3 || !entities.all_of_dimension(dim))
+ return MB_TYPE_OUT_OF_RANGE;
+
+ // create a bit tag for use in a) tagging skin vertices, and
+ // b) fast intersection with input entities range.
+ MBTag tag;
+ char bit = 0;
+ rval = thisMB->tag_create( NULL, 1, use_bit_tag ? MB_TAG_BIT : MB_TAG_DENSE, tag, &bit );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ // tag all entities in input range
+ if (use_bit_tag) {
+ bit = '\001';
+ for (MBRange::const_iterator it = entities.begin(); it != entities.end(); ++it) {
+ rval = thisMB->tag_set_data( tag, &*it, 1, &bit );
+ if (MB_SUCCESS != rval) {
+ thisMB->tag_delete(tag);
+ return rval;
+ }
+ }
+ }
+ else {
+ size_t count = entities.size();
+ char* vect = new char[count];
+ memset( vect, 1, count );
+ rval = thisMB->tag_set_data( tag, entities, vect );
+ delete [] vect;
+ if (MB_SUCCESS != rval) {
+ thisMB->tag_delete(tag);
+ return rval;
+ }
+ }
+
+ switch (dim) {
+ case 1:
+ if (skin_verts)
+ rval = find_skin_vertices_1D( tag, entities, *skin_verts );
+ else if (skin_elems)
+ rval = find_skin_vertices_1D( tag, entities, *skin_elems );
+ else
+ rval = MB_SUCCESS;
+ break;
+ case 2:
+ rval = find_skin_vertices_2D( tag, entities, skin_verts,
+ skin_elems, skin_rev_elems,
+ create_skin_elems, corners_only );
+ break;
+ case 3:
+ rval = find_skin_vertices_3D( tag, entities, skin_verts,
+ skin_elems, skin_rev_elems,
+ create_skin_elems, corners_only );
+ break;
+ default:
+ rval = MB_TYPE_OUT_OF_RANGE;
+ break;
+ }
+
+ thisMB->tag_delete(tag);
+ return rval;
+}
+
+MBErrorCode MBSkinner::find_skin_vertices_1D( MBTag tag,
+ const MBRange& edges,
+ MBRange& skin_verts )
+{
+ MBErrorCode rval;
+ std::vector<MBEntityHandle>::iterator i;
+ MBRange::iterator hint = skin_verts.begin();
+
+ if (!edges.all_of_dimension(1))
+ return MB_TYPE_OUT_OF_RANGE;
+
+ // get all the vertices
+ MBRange verts;
+ rval = thisMB->get_connectivity( edges, verts, true );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ std::vector<char> tag_vals;
+ std::vector<MBEntityHandle> adj;
+ int count;
+ char bit;
+ for (MBRange::const_iterator it = verts.begin(); it != verts.end(); ++it) {
+ adj.clear();
+ rval = thisMB->get_adjacencies( &*it, 1, 1, false, adj );
+ if (MB_SUCCESS != rval) return rval;
+ if (adj.empty())
+ continue;
+
+ // remove those not in the input list
+ if (use_bit_tag) {
+ count = 0;
+ for (i = adj.begin(); i != adj.end(); ++i) {
+ rval = thisMB->tag_get_data( tag, &*i, 1, &bit );
+ if (MB_SUCCESS != rval) return rval;
+ if (bit)
+ ++count;
+ }
+ }
+ else {
+ tag_vals.resize( adj.size() );
+ rval = thisMB->tag_get_data( tag, &adj[0], adj.size(), &tag_vals[0] );
+ if (MB_SUCCESS != rval) return rval;
+ count = std::count( tag_vals.begin(), tag_vals.end(), '\001' );
+ }
+
+ if (count == 1) {
+ hint = skin_verts.insert( hint, *it );
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
+
+// A Container for storing a list of element sides adjacent
+// to a vertex. The template parameter is the number of
+// corners for the side.
+template <unsigned CORNERS>
+class AdjSides
+{
+public:
+ struct Side {
+ MBEntityHandle handles[CORNERS-1];
+ MBEntityHandle adj_elem;
+ bool skin() const { return 0 != adj_elem; }
+
+ Side( const MBEntityHandle* array, int idx,
+ MBEntityHandle adj, unsigned short )
+ : adj_elem(adj) /*, elem_side(side)*/
+ {
+ switch (CORNERS) {
+ default:
+ assert(false);
+ break;
+ case 4: handles[2] = array[(idx+3)%CORNERS];
+ case 3: handles[1] = array[(idx+2)%CORNERS];
+ case 2: handles[0] = array[(idx+1)%CORNERS];
+ }
+ if (CORNERS == 3 && handles[1] > handles[0])
+ std::swap( handles[0], handles[1] );
+ if (CORNERS == 4 && handles[2] > handles[0])
+ std::swap( handles[0], handles[2] );
+ }
+
+ Side( const MBEntityHandle* array, int idx,
+ MBEntityHandle adj, unsigned short ,
+ const short* indices )
+ : adj_elem(adj) /*, elem_side(side)*/
+ {
+ switch (CORNERS) {
+ default:
+ assert(false);
+ break;
+ case 4: handles[2] = array[indices[(idx+3)%CORNERS]];
+ case 3: handles[1] = array[indices[(idx+2)%CORNERS]];
+ case 2: handles[0] = array[indices[(idx+1)%CORNERS]];
+ }
+ if (CORNERS == 3 && handles[1] > handles[0])
+ std::swap( handles[0], handles[1] );
+ if (CORNERS == 4 && handles[2] > handles[0])
+ std::swap( handles[0], handles[2] );
+ }
+
+ bool operator==( const Side& other ) const
+ {
+ switch (CORNERS) {
+ default:
+ assert(false);
+ return false;
+ case 4:
+ return handles[0] == other.handles[0]
+ && handles[1] == other.handles[1]
+ && handles[2] == other.handles[2];
+ case 3:
+ return handles[0] == other.handles[0]
+ && handles[1] == other.handles[1];
+ case 2:
+ return handles[0] == other.handles[0];
+ }
+ }
+ };
+
+private:
+
+ std::vector<Side> data;
+ size_t skin_count;
+
+public:
+
+ typedef typename std::vector<Side>::iterator iterator;
+ typedef typename std::vector<Side>::const_iterator const_iterator;
+ const_iterator begin() const { return data.begin(); }
+ const_iterator end() const { return data.end(); }
+
+ void clear() { data.clear(); skin_count = 0; }
+ bool empty() const { return data.empty(); }
+
+ AdjSides() : skin_count(0) {}
+
+ size_t num_skin() const { return skin_count; }
+
+ void insert( const MBEntityHandle* handles, int skip_idx,
+ MBEntityHandle adj_elem, unsigned short elem_side )
+ {
+ Side side( handles, skip_idx, adj_elem, elem_side );
+ iterator p = std::find( data.begin(), data.end(), side );
+ if (p == data.end()) {
+ data.push_back( side );
+ ++skin_count;
+ }
+ else if (p->adj_elem) {
+ p->adj_elem = 0;
+ --skin_count;
+ }
+ }
+
+ void insert( const MBEntityHandle* handles, int skip_idx,
+ MBEntityHandle adj_elem, unsigned short elem_side,
+ const short* indices )
+ {
+ Side side( handles, skip_idx, adj_elem, elem_side, indices );
+ iterator p = std::find( data.begin(), data.end(), side );
+ if (p == data.end()) {
+ data.push_back( side );
+ ++skin_count;
+ }
+ else if (p->adj_elem) {
+ p->adj_elem = 0;
+ --skin_count;
+ }
+ }
+
+ bool find_and_unmark( const MBEntityHandle* other, int skip_index, MBEntityHandle& elem_out )
+ {
+ Side s( other, skip_index, 0, 0 );
+ iterator p = std::find( data.begin(), data.end(), s );
+ if (p == data.end() || !p->adj_elem)
+ return false;
+ else {
+ elem_out = p->adj_elem;
+ p->adj_elem = 0;
+ --skin_count;
+ return true;
+ }
+ }
+};
+
+MBErrorCode MBSkinner::create_side( MBEntityHandle elem,
+ MBEntityType side_type,
+ const MBEntityHandle* side_conn,
+ MBEntityHandle& side_elem )
+{
+ const int max_side = 9;
+ const MBEntityHandle* conn;
+ int len, side_len, side, sense, offset, indices[max_side];
+ MBErrorCode rval;
+ MBEntityType type = TYPE_FROM_HANDLE(elem), tmp_type;
+ const int ncorner = MBCN::VerticesPerEntity( side_type );
+ const int d = MBCN::Dimension(side_type);
+
+ rval = thisMB->get_connectivity( elem, conn, len, false );
+ if (MB_SUCCESS != rval) return rval;
+
+ MBCN::SideNumber( type, conn, side_conn, ncorner, d, side, sense, offset );
+ MBCN::SubEntityNodeIndices( type, len, d, side, tmp_type, side_len, indices );
+ assert(side_len <= max_side);
+ assert(side_type == tmp_type);
+
+ //NOTE: re-create conn array even when no higher-order nodes
+ // because we want it to always be forward with respect
+ // to the side ordering.
+ MBEntityHandle side_conn_full[max_side];
+ for (int i = 0; i < side_len; ++i)
+ side_conn_full[i] = conn[indices[i]];
+
+ return thisMB->create_element( side_type, side_conn_full, side_len, side_elem );
+}
+
+bool MBSkinner::edge_reversed( MBEntityHandle face,
+ const MBEntityHandle* edge_ends )
+{
+ const MBEntityHandle* conn;
+ int len, idx;
+ MBErrorCode rval = thisMB->get_connectivity( face, conn, len, true );
+ if (MB_SUCCESS != rval) {
+ assert(false);
+ return false;
+ }
+ idx = std::find( conn, conn+len, edge_ends[0] ) - conn;
+ if (idx == len) {
+ assert(false);
+ return false;
+ }
+ return (edge_ends[1] == conn[(idx+len-1)%len]);
+}
+
+bool MBSkinner::face_reversed( MBEntityHandle region,
+ const MBEntityHandle* face_corners,
+ MBEntityType face_type )
+{
+ const MBEntityHandle* conn;
+ int len, side, sense, offset;
+ MBErrorCode rval = thisMB->get_connectivity( region, conn, len, true );
+ if (MB_SUCCESS != rval) {
+ assert(false);
+ return false;
+ }
+ short r = MBCN::SideNumber( TYPE_FROM_HANDLE(region), conn, face_corners,
+ MBCN::VerticesPerEntity(face_type),
+ MBCN::Dimension(face_type),
+ side, sense, offset );
+ assert(0 == r);
+ return (!r && sense == -1);
+}
+
+MBErrorCode MBSkinner::find_skin_vertices_2D( MBTag tag,
+ const MBRange& faces,
+ MBRange* skin_verts,
+ MBRange* skin_edges,
+ MBRange* reversed_edges,
+ bool create_edges,
+ bool corners_only )
+{
+ MBErrorCode rval;
+ std::vector<MBEntityHandle>::iterator i, j;
+ MBRange::iterator hint;
+ if (skin_verts)
+ hint = skin_verts->begin();
+ std::vector<MBEntityHandle> storage;
+ const MBEntityHandle *conn;
+ int len;
+ char bit;
+ bool find_edges = skin_edges || create_edges;
+ MBEntityHandle face;
+
+ if (!faces.all_of_dimension(2))
+ return MB_TYPE_OUT_OF_RANGE;
+
+ // get all the vertices
+ MBRange verts;
+ rval = thisMB->get_connectivity( faces, verts, true );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ std::vector<char> tag_vals;
+ std::vector<MBEntityHandle> adj;
+ AdjSides<2> adj_edges;
+ for (MBRange::const_iterator it = verts.begin(); it != verts.end(); ++it) {
+ bool higher_order = false;
+
+ // get all adjacent faces
+ adj.clear();
+ rval = thisMB->get_adjacencies( &*it, 1, 2, false, adj );
+ if (MB_SUCCESS != rval) return rval;
+ if (adj.empty())
+ continue;
+
+ // remove those not in the input list
+ i = j = adj.begin();
+ if (use_bit_tag) {
+ for (; i != adj.end(); ++i) {
+ rval = thisMB->tag_get_data( tag, &*i, 1, &bit );
+ if (MB_SUCCESS != rval) return rval;
+ if (bit)
+ *(j++) = *i;
+ }
+ adj.erase( j, adj.end() );
+ }
+ else {
+ tag_vals.resize( adj.size() );
+ rval = thisMB->tag_get_data( tag, &adj[0], adj.size(), &tag_vals[0] );
+ if (MB_SUCCESS != rval) return rval;
+
+ i = j = adj.begin();
+ for (; i != adj.end(); ++i)
+ if (tag_vals[i - adj.begin()])
+ *(j++) = *i;
+ adj.erase( j, adj.end() );
+ }
+
+ // For each adjacent face, check the edges adjacent to the current vertex
+ adj_edges.clear(); // other vertex for adjacent edges
+ for (i = adj.begin(); i != adj.end(); ++i) {
+ rval = thisMB->get_connectivity( *i, conn, len, false, &storage );
+ if (MB_SUCCESS != rval) return rval;
+
+ MBEntityHandle prev, next; // vertices of two adjacent edge-sides
+ const int idx = std::find(conn, conn+len, *it) - conn;
+ assert(idx != len);
+
+ if (TYPE_FROM_HANDLE(*i) == MBTRI && len > 3) {
+ len = 3;
+ higher_order = true;
+ if (idx > 2) // skip higher-order nodes for now
+ continue;
+ }
+ else if (TYPE_FROM_HANDLE(*i) == MBQUAD && len > 4) {
+ len = 4;
+ higher_order = true;
+ if (idx > 3) // skip higher-order nodes for now
+ continue;
+ }
+
+ const int prev_idx = (idx + len - 1)%len;
+ prev = conn[prev_idx];
+ next = conn[(idx+1)%len];
+ adj_edges.insert( &prev, 1, *i, prev_idx );
+ adj_edges.insert( &next, 1, *i, idx );
+ }
+
+ // If vertex is not on skin, advance to next vertex
+ if (0 == adj_edges.num_skin())
+ continue;
+
+ // Put skin vertex in output list
+ if (skin_verts) {
+ hint = skin_verts->insert( hint, *it );
+
+ // Add mid edge nodes to vertex list
+ if (!corners_only && higher_order) {
+ for (AdjSides<2>::const_iterator p = adj_edges.begin(); p != adj_edges.end(); ++p) {
+ if (p->skin()) {
+ face = p->adj_elem;
+ MBEntityType type = TYPE_FROM_HANDLE(face);
+
+ rval = thisMB->get_connectivity( face, conn, len, false );
+ if (MB_SUCCESS != rval) return rval;
+ if (!MBCN::HasMidEdgeNodes( type, len ))
+ continue;
+
+ MBEntityHandle ec[2] = { *it, p->handles[0] };
+ int side, sense, offset;
+ MBCN::SideNumber( type, conn, ec, 2, 1, side, sense, offset );
+ offset = MBCN::HONodeIndex( type, len, 1, side );
+ assert(offset >= 0 && offset < len);
+ skin_verts->insert( conn[offset] );
+ }
+ }
+ }
+ }
+
+ if (find_edges) {
+ // Search list of adjacent edges for any that are on the skin
+ adj.clear();
+ rval = thisMB->get_adjacencies( &*it, 1, 1, false, adj );
+ if (MB_SUCCESS != rval) return rval;
+ for (i = adj.begin(); i != adj.end(); ++i) {
+ rval = thisMB->get_connectivity( *i, conn, len, true );
+ if (MB_SUCCESS != rval) return rval;
+
+ // bool equalality expression will be evaluate to the
+ // index of *it in the conn array. Note that the order
+ // of the terms in the if statement is important. We want
+ // to unmark any existing skin edges even if we aren't returning
+ // them. Otherwise we'll end up creating duplicates if create_edges
+ // is true.
+ if (adj_edges.find_and_unmark( conn, (conn[1] == *it), face ) && skin_edges) {
+ if (reversed_edges && edge_reversed( face, conn ))
+ reversed_edges->insert( *i );
+ else
+ skin_edges->insert( *i );
+ }
+ }
+ }
+
+ if (create_edges && adj_edges.num_skin()) {
+ // Create any skin edges that don't exist
+ for (AdjSides<2>::const_iterator p = adj_edges.begin(); p != adj_edges.end(); ++p) {
+ if (p->skin()) {
+ MBEntityHandle edge, ec[] = { *it, p->handles[0] };
+ rval = create_side( p->adj_elem, MBEDGE, ec, edge );
+ if (MB_SUCCESS != rval) return rval;
+ if (skin_edges)
+ skin_edges->insert( edge );
+ }
+ }
+ }
+
+ } // end for each vertex
+
+ return MB_SUCCESS;
+}
+
+
+MBErrorCode MBSkinner::find_skin_vertices_3D( MBTag tag,
+ const MBRange& entities,
+ MBRange* skin_verts,
+ MBRange* skin_faces,
+ MBRange* reversed_faces,
+ bool create_faces,
+ bool corners_only )
+{
+ MBErrorCode rval;
+ std::vector<MBEntityHandle>::iterator i, j;
+ MBRange::iterator hint;
+ if (skin_verts)
+ hint = skin_verts->begin();
+ std::vector<MBEntityHandle> storage, storage2;
+ const MBEntityHandle *conn, *conn2;
+ int len, len2;
+ char bit;
+ bool find_faces = skin_faces || create_faces;
+ int clen, side, sense, offset, indices[9];
+ MBEntityType face_type;
+ MBEntityHandle elem;
+
+ if (!entities.all_of_dimension(3))
+ return MB_TYPE_OUT_OF_RANGE;
+
+ MBRange verts;
+ rval = thisMB->get_adjacencies( entities, 0, false, verts, MBInterface::UNION );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ AdjSides<4> adj_quads;
+ AdjSides<3> adj_tris;
+ AdjSides<2> adj_poly;
+ std::vector<char> tag_vals;
+ std::vector<MBEntityHandle> adj;
+ for (MBRange::const_iterator it = verts.begin(); it != verts.end(); ++it) {
+ bool higher_order = false;
+
+ // get all adjacent elements
+ adj.clear();
+ rval = thisMB->get_adjacencies( &*it, 1, 3, false, adj );
+ if (MB_SUCCESS != rval) return rval;
+ if (adj.empty())
+ continue;
+
+ // remove those not in the input list
+ i = j = adj.begin();
+ if (use_bit_tag) {
+ for (; i != adj.end(); ++i) {
+ rval = thisMB->tag_get_data( tag, &*i, 1, &bit );
+ if (MB_SUCCESS != rval) return rval;
+ if (bit)
+ *(j++) = *i;
+ }
+ adj.erase( j, adj.end() );
+ }
+ else {
+ tag_vals.resize( adj.size() );
+ rval = thisMB->tag_get_data( tag, &adj[0], adj.size(), &tag_vals[0] );
+ if (MB_SUCCESS != rval) return rval;
+ for (; i != adj.end(); ++i)
+ if (tag_vals[i - adj.begin()])
+ *(j++) = *i;
+ adj.erase( j, adj.end() );
+ }
+
+ // Build lists of sides of 3D element adjacent to the current vertex
+ adj_quads.clear(); // store three other vertices for each adjacent quad face
+ adj_tris.clear(); // store two other vertices for each adjacent tri face
+ adj_poly.clear(); // store handle of each adjacent polygonal face
+ int idx;
+ for (i = adj.begin(); i != adj.end(); ++i) {
+ const MBEntityType type = TYPE_FROM_HANDLE(*i);
+
+ // Special case for POLYHEDRA
+ if (type == MBPOLYHEDRON) {
+ rval = thisMB->get_connectivity( *i, conn, len );
+ if (MB_SUCCESS != rval) return rval;
+ for (int k = 0; k < len; ++k) {
+ rval = thisMB->get_connectivity( conn[k], conn2, len2, true, &storage2 );
+ if (MB_SUCCESS != rval) return rval;
+ idx = std::find( conn2, conn2+len2, *it) - conn2;
+ if (idx == len2) // vertex not in this face
+ continue;
+
+ // Treat 3- and 4-vertex faces specially, so that
+ // if the mesh contains both elements and polyhedra,
+ // we don't miss one type adjacent to the other.
+ switch (len2) {
+ case 3:
+ adj_tris.insert( conn2, idx, *i, k );
+ break;
+ case 4:
+ adj_quads.insert( conn2, idx, *i, k );
+ break;
+ default:
+ adj_poly.insert( conn+k, 1, *i, k );
+ break;
+ }
+ }
+ }
+ else {
+ rval = thisMB->get_connectivity( *i, conn, len, false, &storage );
+ if (MB_SUCCESS != rval) return rval;
+
+ idx = std::find(conn, conn+len, *it) - conn;
+ assert(idx != len);
+
+ if (len > MBCN::VerticesPerEntity( type )) {
+ higher_order =true;
+ // skip higher-order nodes for now
+ if (idx >= MBCN::VerticesPerEntity( type ))
+ continue;
+ }
+
+ const int num_faces = MBCN::NumSubEntities( type, 2 );
+ for (int f = 0; f < num_faces; ++f) {
+ int num_vtx;
+ const short* face_indices = MBCN::SubEntityVertexIndices(type, 2, f, face_type, num_vtx );
+ const short face_idx = std::find(face_indices, face_indices+num_vtx, (short)idx) - face_indices;
+ if (face_idx == num_vtx)
+ continue; // current vertex not in this face
+
+ assert(num_vtx <= 4);
+ switch (face_type) {
+ case MBTRI:
+ adj_tris.insert( conn, face_idx, *i, f, face_indices );
+ break;
+ case MBQUAD:
+ adj_quads.insert( conn, face_idx, *i, f, face_indices );
+ break;
+ default:
+ return MB_TYPE_OUT_OF_RANGE;
+ }
+ }
+ }
+ } // end for (adj[3])
+
+ // If vertex is not on skin, advance to next vertex
+ if (0 == (adj_tris.num_skin() + adj_quads.num_skin() + adj_poly.num_skin()))
+ continue;
+
+ // Put skin vertex in output list
+ if (skin_verts) {
+ hint = skin_verts->insert( hint, *it );
+
+ // Add mid-edge and mid-face nodes to vertex list
+ if (!corners_only && higher_order) {
+ for (AdjSides<3>::const_iterator t = adj_tris.begin(); t != adj_tris.end(); ++t) {
+ if (t->skin()) {
+ elem = t->adj_elem;
+ MBEntityType type = TYPE_FROM_HANDLE(elem);
+
+ rval = thisMB->get_connectivity( elem, conn, len, false );
+ if (MB_SUCCESS != rval) return rval;
+ if (!MBCN::HasMidNodes( type, len ))
+ continue;
+
+ MBEntityHandle ec[3] = { *it, t->handles[0], t->handles[1] };
+ MBCN::SideNumber( type, conn, ec, 3, 2, side, sense, offset );
+ MBCN::SubEntityNodeIndices( type, len, 2, side, face_type, clen, indices );
+ assert(MBTRI == face_type);
+ for (int k = 3; k < clen; ++k)
+ skin_verts->insert( conn[indices[k]] );
+ }
+ }
+ for (AdjSides<4>::const_iterator q = adj_quads.begin(); q != adj_quads.end(); ++q) {
+ if (q->skin()) {
+ elem = q->adj_elem;
+ MBEntityType type = TYPE_FROM_HANDLE(elem);
+
+ rval = thisMB->get_connectivity( elem, conn, len, false );
+ if (MB_SUCCESS != rval) return rval;
+ if (!MBCN::HasMidNodes( type, len ))
+ continue;
+
+ MBEntityHandle ec[4] = { *it, q->handles[0], q->handles[1], q->handles[2] };
+ MBCN::SideNumber( type, conn, ec, 4, 2, side, sense, offset );
+ MBCN::SubEntityNodeIndices( type, len, 2, side, face_type, clen, indices );
+ assert(MBQUAD == face_type);
+ for (int k = 4; k < clen; ++k)
+ skin_verts->insert( conn[indices[k]] );
+ }
+ }
+ }
+ }
+
+ if (find_faces) {
+ // Search list of adjacent faces for any that are on the skin
+ adj.clear();
+ rval = thisMB->get_adjacencies( &*it, 1, 2, false, adj );
+ if (MB_SUCCESS != rval) return rval;
+
+ for (i = adj.begin(); i != adj.end(); ++i) {
+ rval = thisMB->get_connectivity( *i, conn, len, true );
+ if (MB_SUCCESS != rval) return rval;
+ const int idx = std::find( conn, conn+len, *it ) - conn;
+ assert(idx != len);
+ // Note that the order of the terms in the if statements below
+ // is important. We want to unmark any existing skin faces even
+ // if we aren't returning them. Otherwise we'll end up creating
+ // duplicates if create_faces is true.
+ if (3 == len) {
+ if (adj_tris.find_and_unmark( conn, idx, elem ) && skin_faces) {
+ if (reversed_faces && face_reversed( elem, conn, MBTRI ))
+ reversed_faces->insert( *i );
+ else
+ skin_faces->insert( *i );
+ }
+ }
+ else if (4 == len) {
+ if (adj_quads.find_and_unmark( conn, idx, elem ) && skin_faces) {
+ if (reversed_faces && face_reversed( elem, conn, MBQUAD ))
+ reversed_faces->insert( *i );
+ else
+ skin_faces->insert( *i );
+ }
+ }
+ else {
+ if (adj_poly.find_and_unmark( &*i, 1, elem ) && skin_faces)
+ skin_faces->insert( *i );
+ }
+ }
+ }
+
+ if (!create_faces)
+ continue;
+
+ // Polyhedra always have explictly defined faces, so
+ // there is no way we could need to create such a face.
+ assert(0 == adj_poly.num_skin());
+
+ // Create any skin tris that don't exist
+ if (adj_tris.num_skin()) {
+ for (AdjSides<3>::const_iterator t = adj_tris.begin(); t != adj_tris.end(); ++t) {
+ if (t->skin()) {
+ MBEntityHandle tri, c[3] = { *it, t->handles[0], t->handles[1] };
+ rval = create_side( t->adj_elem, MBTRI, c, tri );
+ if (MB_SUCCESS != rval) return rval;
+ if (skin_faces)
+ skin_faces->insert( tri );
+ }
+ }
+ }
+
+ // Create any skin quads that don't exist
+ if (adj_quads.num_skin()) {
+ for (AdjSides<4>::const_iterator q = adj_quads.begin(); q != adj_quads.end(); ++q) {
+ if (q->skin()) {
+ MBEntityHandle quad, c[4] = { *it, q->handles[0], q->handles[1], q->handles[2] };
+ rval = create_side( q->adj_elem, MBQUAD, c, quad );
+ if (MB_SUCCESS != rval) return rval;
+ if (skin_faces)
+ skin_faces->insert( quad );
+ }
+ }
+ }
+ } // end for each vertex
+
+ return MB_SUCCESS;
+}
Modified: MOAB/trunk/MBSkinner.hpp
===================================================================
--- MOAB/trunk/MBSkinner.hpp 2009-11-19 23:08:51 UTC (rev 3372)
+++ MOAB/trunk/MBSkinner.hpp 2009-11-19 23:10:26 UTC (rev 3373)
@@ -46,9 +46,11 @@
// will accept entities all of one dimension
// and return entities of n-1 dimension
MBErrorCode find_skin( const MBRange &entities,
- MBRange &forward_lower_entities,
- MBRange &reverse_lower_entities,
- bool create_vert_elem_adjs = false);
+ bool get_vertices,
+ MBRange &output_handles,
+ MBRange *output_reverse_handles = 0,
+ bool create_vert_elem_adjs = false,
+ bool create_skin_elements = true);
// get skin entities of prescribed dimension
MBErrorCode find_skin(const MBRange &entities,
@@ -56,6 +58,24 @@
MBRange &skin_entities,
bool create_vert_elem_adjs = false);
+ /**\brief Find vertices on the skin of a set of mesh entities.
+ *\param entities The elements for which to find the skin. Range
+ * may NOT contain vertices, polyhedra, or entity sets.
+ * All elements in range must be of the same dimension.
+ *\param skin_verts Output: the vertices on the skin.
+ *\param skin_elems Optional output: elements representing sides of entities
+ * that are on the skin
+ *\param create_if_missing If skin_elemts is non-null and this is true,
+ * create new elements representing the sides of
+ * entities on the skin. If this is false, skin_elems
+ * will contain only those skin elements that already
+ * exist.
+ */
+ MBErrorCode find_skin_vertices( const MBRange& entities,
+ MBRange& skin_verts,
+ MBRange* skin_elems = 0,
+ bool create_if_missing = true );
+
MBErrorCode classify_2d_boundary( const MBRange &boundary,
const MBRange &bar_elements,
MBEntityHandle boundary_edges,
@@ -80,6 +100,10 @@
void deinitialize();
+ MBErrorCode find_skin_noadj( const MBRange &source_entities,
+ MBRange &forward_target_entities,
+ MBRange &reverse_target_entities );
+
void add_adjacency(MBEntityHandle entity);
void add_adjacency(MBEntityHandle entity, const MBEntityHandle *conn,
@@ -109,6 +133,99 @@
MBEntityHandle &entity2,
double reference_angle_cosine);
+
+ /**\brief Find vertices on the skin of a set of mesh entities.
+ *\param entities The elements for which to find the skin. Range
+ * may NOT contain vertices, polyhedra, or entity sets.
+ * All elements in range must be of the same dimension.
+ *\param skin_verts Output: the vertices on the skin.
+ *\param skin_elems Optional output: elements representing sides of entities
+ * that are on the skin
+ *\param create_if_missing If skin_elemts is non-null and this is true,
+ * create new elements representing the sides of
+ * entities on the skin. If this is false, skin_elems
+ * will contain only those skin elements that already
+ * exist.
+ */
+ MBErrorCode find_skin_vertices( const MBRange& entities,
+ MBRange* skin_verts = 0,
+ MBRange* skin_elems = 0,
+ MBRange* rev_elems = 0,
+ bool create_if_missing = true,
+ bool corners_only = false );
+
+ /**\brief Skin edges
+ *
+ * Return any vertices adjacent to exactly one of the input edges.
+ */
+ MBErrorCode find_skin_vertices_1D( MBTag tag,
+ const MBRange& edges,
+ MBRange& skin_verts );
+
+ /**\brief Skin faces
+ *
+ * For the set of face sides (logical edges), return
+ * vertices on such sides and/or edges equivalent to such sides.
+ *\param faces Set of toplogically 2D entities to skin.
+ *\param skin_verts If non-NULL, skin vertices will be added to this container.
+ *\param skin_edges If non-NULL, skin edges will be added to this container
+ *\param reverse_edges If skin_edges is not NULL and this is not NULL, then
+ * any existing skin edges that are reversed with respect
+ * to the skin side will be placed in this range instead of
+ * skin_edges. Note: this argument is ignored if skin_edges
+ * is NULL.
+ *\param create_edges If true, edges equivalent to face sides on the skin
+ * that don't already exist will be created. Note: this
+ * parameter is honored regardless of whether or not skin
+ * edges or vertices are returned.
+ *\param corners_only If true, only skin vertices that correspond to the
+ * corners of sides will be returned (i.e. no higher-order
+ * nodes.) This argument is ignored if skin_verts is NULL.
+ */
+ MBErrorCode find_skin_vertices_2D( MBTag tag,
+ const MBRange& faces,
+ MBRange* skin_verts = 0,
+ MBRange* skin_edges = 0,
+ MBRange* reverse_edges = 0,
+ bool create_edges = false,
+ bool corners_only = false );
+
+ /**\brief Skin volume mesh
+ *
+ * For the set of element sides (logical faces), return
+ * vertices on such sides and/or faces equivalent to such sides.
+ *\param entities Set of toplogically 3D entities to skin.
+ *\param skin_verts If non-NULL, skin vertices will be added to this container.
+ *\param skin_faces If non-NULL, skin faces will be added to this container
+ *\param reverse_faces If skin_faces is not NULL and this is not NULL, then
+ * any existing skin faces that are reversed with respect
+ * to the skin side will be placed in this range instead of
+ * skin_faces. Note: this argument is ignored if skin_faces
+ * is NULL.
+ *\param create_faces If true, face equivalent to sides on the skin
+ * that don't already exist will be created. Note: this
+ * parameter is honored regardless of whether or not skin
+ * faces or vertices are returned.
+ *\param corners_only If true, only skin vertices that correspond to the
+ * corners of sides will be returned (i.e. no higher-order
+ * nodes.) This argument is ignored if skin_verts is NULL.
+ */
+ MBErrorCode find_skin_vertices_3D( MBTag tag,
+ const MBRange& entities,
+ MBRange* skin_verts = 0,
+ MBRange* skin_faces = 0,
+ MBRange* reverse_faces = 0,
+ bool create_faces = false,
+ bool corners_only = false );
+
+ MBErrorCode create_side( MBEntityHandle element,
+ MBEntityType side_type,
+ const MBEntityHandle* side_corners,
+ MBEntityHandle& side_elem_handle_out );
+
+ bool edge_reversed( MBEntityHandle face, const MBEntityHandle edge_ends[2] );
+ bool face_reversed( MBEntityHandle region, const MBEntityHandle* face_conn,
+ MBEntityType face_type );
};
Modified: MOAB/trunk/MBTest.cpp
===================================================================
--- MOAB/trunk/MBTest.cpp 2009-11-19 23:08:51 UTC (rev 3372)
+++ MOAB/trunk/MBTest.cpp 2009-11-19 23:10:26 UTC (rev 3373)
@@ -28,6 +28,7 @@
#endif
#include <iostream>
+#include <fstream>
#include <algorithm>
#include <cstdio>
#include <time.h>
@@ -4598,7 +4599,7 @@
//get Hexes from model
}
result = MB->get_entities_by_type(0, MBHEX, entities);
- MBSkinner_Obj.find_skin(entities,forward_lower,reverse_lower);
+ MBSkinner_Obj.find_skin(entities,false,forward_lower,&reverse_lower);
cout <<"num hexes = "<<entities.size()<<"\n";
cout <<"fl = "<<forward_lower.size()<<" rl = "<<reverse_lower.size()<<"\n";
@@ -6562,7 +6563,7 @@
MBErrorCode mb_skin_volume_adj_test( MBInterface* )
{ return mb_skin_volume_test_common( true ); }
-
+
MBErrorCode mb_skin_curve_test_common( bool use_adj )
{
MBErrorCode rval;
@@ -6884,6 +6885,774 @@
}
+// Test basic skinning using vert-to-elem adjacencies
+MBErrorCode mb_skin_verts_common( unsigned dim, bool skin_elems );
+
+MBErrorCode mb_skin_surf_verts_test( MBInterface* )
+ { return mb_skin_verts_common( 2, false ); }
+
+MBErrorCode mb_skin_vol_verts_test( MBInterface* )
+ { return mb_skin_verts_common( 3, false ); }
+
+MBErrorCode mb_skin_surf_verts_elems_test( MBInterface* )
+ { return mb_skin_verts_common( 2, true ); }
+
+MBErrorCode mb_skin_vol_verts_elems_test( MBInterface* )
+ { return mb_skin_verts_common( 3, true ); }
+
+MBErrorCode mb_skin_verts_common( unsigned dim, bool skin_elems )
+{
+ const int INT = 10; // intervals+1
+ const char* tmp_file = "structured.vtk";
+ std::ofstream str(tmp_file);
+ if (!str) {
+ std::cerr << tmp_file << ": filed to create temp file" << std::endl;
+ return MB_FAILURE;
+ }
+ str << "#vtk DataFile Version 2.0" << std::endl
+ << "mb_skin_verts_common temp file" << std::endl
+ << "ASCII" << std::endl
+ << "DATASET STRUCTURED_POINTS" << std::endl
+ << "DIMENSIONS " << INT << " " << (dim > 1 ? INT : 1) << " " << (dim > 2 ? INT : 1) << std::endl
+ << "ORIGIN 0 0 0" << std::endl
+ << "SPACING 1 1 1" << std::endl;
+ str.close();
+
+ MBCore moab;
+ MBInterface& mb = moab;
+ MBErrorCode rval;
+
+ rval = mb.load_file( tmp_file );
+ remove( tmp_file );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ MBRange ents;
+ rval = mb.get_entities_by_dimension( 0, dim, ents );
+ if (MB_SUCCESS != rval)
+ return rval;
+ if (ents.empty())
+ return MB_FAILURE;
+
+ MBSkinner tool( &mb );
+
+ // mesh is a structured quad/hex mesh, so we can
+ // determine skin vertices from the number of
+ // adjacent elements.
+ unsigned interior_adj = 1;
+ for (unsigned i = 0; i < dim; ++i)
+ interior_adj *= 2;
+ MBRange expected, verts;
+ rval = mb.get_entities_by_dimension( 0, 0, verts );
+ if (MB_SUCCESS != rval)
+ return rval;
+ MBRange::iterator h = expected.begin();
+ std::vector<MBEntityHandle> adj;
+ for (MBRange::iterator v = verts.begin(); v != verts.end(); ++v) {
+ adj.clear();
+ rval = mb.get_adjacencies( &*v, 1, dim, false, adj );
+ if (MB_SUCCESS != rval)
+ return rval;
+ if (adj.size() < interior_adj)
+ h = expected.insert( h, *v );
+ }
+
+ // Get skin vertices using skinner
+ MBRange actual;
+ rval = tool.find_skin( ents, !skin_elems, actual );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ MBRange extra, missing;
+ if (!skin_elems) {
+ // Check that we got expected result
+ extra = subtract( actual, expected );
+ missing = subtract( expected, actual );
+ if (!extra.empty() || !missing.empty()) {
+ std::cout << "Extra vertices returned: " << extra << std::endl
+ << "Missing vertices: " << missing << std::endl;
+ return MB_FAILURE;
+ }
+ return MB_SUCCESS;
+ }
+
+ // test that no extra elements we're created
+ extra.clear();
+ rval = mb.get_entities_by_dimension( 0, dim-1, extra );
+ if (MB_SUCCESS != rval)
+ return rval;
+ extra = subtract( extra, actual );
+ if (!extra.empty()) {
+ std::cout << "Extra/non-returned elements created: " << extra << std::endl;
+ return MB_FAILURE;
+ }
+
+ // check that each skin vertex has the correct number of adjacent quads
+ missing.clear(); extra.clear();
+ for (MBRange::iterator i = expected.begin(); i != expected.end(); ++i) {
+ std::vector<MBEntityHandle> elem, side;
+ rval = mb.get_adjacencies( &*i, 1, dim, false, elem );
+ if (MB_SUCCESS != rval) return rval;
+ rval = mb.get_adjacencies( &*i, 1, dim-1, false, side );
+ if (MB_SUCCESS != rval) return rval;
+ if (elem.size() == 1) {
+ if (side.size() < dim)
+ missing.insert( *i );
+ else if(side.size() > dim)
+ extra.insert( *i );
+ }
+ else if (elem.size() == interior_adj) {
+ if (!side.empty())
+ extra.insert( *i );
+ }
+ else {
+ if (side.size() < interior_adj/2)
+ missing.insert( *i );
+ else if (side.size() > interior_adj/2)
+ extra.insert( *i );
+ }
+ }
+ if (!missing.empty() || !extra.empty()) {
+ std::cout << "Skin elements missing at vertices: " << missing << std::endl
+ << "Extra skin elements at vertices: " << extra << std::endl;
+ return MB_FAILURE;
+ }
+
+ // check that all returned elements are actually on the skin
+ extra.clear();
+ for (MBRange::iterator i = actual.begin(); i != actual.end(); ++i) {
+ MBRange verts;
+ rval = mb.get_adjacencies( &*i, 1, 0, false, verts );
+ if (MB_SUCCESS != rval)
+ return rval;
+ verts = subtract( verts, expected );
+ if (!verts.empty())
+ extra.insert( *i );
+ }
+ if (!extra.empty()) {
+ std::cout << "Skinner returned elements not on skin: " << extra << std::endl;
+ return MB_FAILURE;
+ }
+
+ return MB_SUCCESS;
+}
+
+// Test that skinning of polyhedra works
+MBErrorCode mb_skin_poly_test( MBInterface* )
+{
+ /* Create a mesh composed of 8 hexagonal prisms and
+ two hexahedra by extruding the following cross section
+ two steps in the Z direction.
+
+ 0-----1
+ / (0) \
+ (11)/ \(1)
+ / \
+ 11 2
+ / \ Y / \
+ (10)/ \(12)^(13)/ \(2)
+ / \ | / \
+ 10 12-----13 3
+ | | | | |
+ (9)| | +--|-->X |(3)
+ | | | |
+ 9 15-----14 4
+ \ / \ /
+ (8)\ /(15) (14)\ /(4)
+ \ / \ /
+ 8 5
+ \ /
+ (7)\ /(5)
+ \ (6) /
+ 7-----6
+ */
+
+ const double coords2D[][2] = { {-1, 5}, // 0
+ { 1, 5},
+ { 3, 3},
+ { 5, 1},
+ { 5,-1},
+ { 3,-3}, // 5
+ { 1,-5},
+ {-1,-5},
+ {-3,-3},
+ {-5,-1},
+ {-5, 1}, // 10
+ {-3, 3},
+ {-1, 1},
+ { 1, 1},
+ { 1,-1},
+ {-1,-1} // 15
+ };
+ const int polyconn[4][6] = { { 0, 1, 2, 13, 12, 11 },
+ { 2, 3, 4, 5, 14, 13 },
+ { 5, 6, 7, 8, 15, 14 },
+ { 8, 9, 10, 11, 12, 15 } };
+ const int polyside[4][6] = { { 0, 1, 13, 16, 12, 11 },
+ { 2, 3, 4, 14, 17, 13 },
+ { 5, 6, 7, 15, 18, 14 },
+ { 8, 9, 10, 12, 19, 15 } };
+
+ MBErrorCode rval;
+ MBCore moab;
+ MBInterface& mb = moab;
+ MBRange regions, faces, interior_faces;
+
+ // create 48 vertices
+ MBEntityHandle verts[3][16];
+ for (int i = 0; i < 3; ++i) {
+ for (int j = 0; j < 16; ++j) {
+ double coords[3] = { coords2D[j][0], coords2D[j][1], 2*i };
+ rval = mb.create_vertex( coords, verts[i][j] );
+ if (MB_SUCCESS != rval) return rval;
+ }
+ }
+
+ // create two hexahedra
+ MBEntityHandle hexes[2];
+ for (int i = 0; i < 2; ++i) {
+ MBEntityHandle conn[8] = { verts[i ][15],
+ verts[i ][14],
+ verts[i ][13],
+ verts[i ][12],
+ verts[i+1][15],
+ verts[i+1][14],
+ verts[i+1][13],
+ verts[i+1][12] };
+ rval = mb.create_element( MBHEX, conn, 8, hexes[i] );
+ if (MB_SUCCESS != rval) return rval;
+ regions.insert(hexes[i]);
+ }
+
+ // create hexagonal faces
+ MBEntityHandle hexagons[3][4];
+ for (int i = 0; i < 3; ++i) {
+ for (int j = 0; j < 4; ++j) {
+ MBEntityHandle conn[6];
+ for (int k = 0; k < 6; ++k)
+ conn[k] = verts[i][polyconn[j][k]];
+ rval = mb.create_element( MBPOLYGON, conn, 6, hexagons[i][j] );
+ if (MB_SUCCESS != rval) return rval;
+ faces.insert( hexagons[i][j] );
+ if (i == 1)
+ interior_faces.insert( hexagons[i][j] );
+ }
+ }
+
+ // create quadrilateral faces
+ MBEntityHandle quads[2][20];
+ for (int i = 0; i < 2; ++i) {
+ for (int j = 0; j < 20; ++j) {
+ int c1, c2;
+ if (j < 12) {
+ c1 = j; c2 = (j+1)%12;
+ }
+ else if (j < 16) {
+ c1 = j; c2 = 2 + 3*((j-9)%4);
+ }
+ else {
+ c1 = j-4; c2 = 12 + (j-15)%4;
+ }
+ MBEntityHandle conn[4] = { verts[i ][c1],
+ verts[i ][c2],
+ verts[i+1][c2],
+ verts[i+1][c1] };
+ rval = mb.create_element( MBQUAD, conn, 4, quads[i][j] );
+ if (MB_SUCCESS != rval) return rval;
+ faces.insert( quads[i][j] );
+ if (j > 11)
+ interior_faces.insert( quads[i][j] );
+ }
+ }
+
+ // create polyhedra
+ MBEntityHandle poly[2][4];
+ for (int i = 0; i < 2; ++i) {
+ for (int j = 0; j < 4; ++j) {
+ MBEntityHandle conn[8];
+ for (int k = 0; k < 6; ++k)
+ conn[k] = quads[i][polyside[j][k]];
+ conn[6] = hexagons[ i][j];
+ conn[7] = hexagons[i+1][j];
+ rval = mb.create_element( MBPOLYHEDRON, conn, 8, poly[i][j] );
+ if (MB_SUCCESS != rval) return rval;
+ regions.insert( poly[i][j] );
+ }
+ }
+
+ MBRange interior_verts;
+ interior_verts.insert( verts[1][12] );
+ interior_verts.insert( verts[1][13] );
+ interior_verts.insert( verts[1][14] );
+ interior_verts.insert( verts[1][15] );
+
+ MBSkinner tool(&mb);
+ MBRange skin;
+ rval = tool.find_skin( regions, true, skin, 0, true, false );
+ if (MB_SUCCESS != rval) {
+ std::cout << "Vertex skinning failed with: " << mb.get_error_string(rval) << std::endl;
+ return rval;
+ }
+
+ MBRange all_verts, all_faces;
+ rval = mb.get_entities_by_dimension( 0, 0, all_verts );
+ if (MB_SUCCESS != rval) return rval;
+ rval = mb.get_entities_by_dimension( 0, 2, all_faces );
+ if (MB_SUCCESS != rval) return rval;
+
+ MBRange expected = subtract( all_verts, interior_verts );
+ if (expected != skin) {
+ std::cout << "Skinner returned incorrect vertices." << std::endl;
+ return MB_FAILURE;
+ }
+ if (all_faces != faces) {
+ std::cout << "Skinner created/deleted faces for vertex-only skinning" << std::endl;
+ return MB_FAILURE;
+ }
+
+ skin.clear();
+ rval = tool.find_skin( regions, false, skin, 0, true, false );
+ if (MB_SUCCESS != rval) {
+ std::cout << "Non-create face skinning failed with: " << mb.get_error_string(rval) << std::endl;
+ return rval;
+ }
+ expected = subtract( all_faces, interior_faces );
+ if (expected != skin) {
+ std::cout << "Skinner returned incorrect faces." << std::endl;
+ return MB_FAILURE;
+ }
+ if (all_faces != faces) {
+ std::cout << "Skinner created/deleted faces for no-create skinning" << std::endl;
+ return MB_FAILURE;
+ }
+
+ skin.clear();
+ rval = tool.find_skin( regions, false, skin, 0, true, true );
+ if (MB_SUCCESS != rval) {
+ std::cout << "Create face skinning failed with: " << mb.get_error_string(rval) << std::endl;
+ return rval;
+ }
+ MBRange all_faces2;
+ rval = mb.get_entities_by_dimension( 0, 2, all_faces2 );
+ if (MB_SUCCESS != rval) return rval;
+ MBRange difference = subtract( all_faces2, all_faces );
+ if (difference.size() != 2) { // should have created two quads for hex top/bottom
+ std::cout << "Skinner failed to create new quads or created to many." << std::endl;
+ return MB_FAILURE;
+ }
+ expected.merge(difference);
+ if (expected != skin) {
+ std::cout << "Skinner returned incorrect faces." << std::endl;
+ return MB_FAILURE;
+ }
+ // check that new faces are correct
+ MBEntityHandle expected_conn[2][4] = {
+ { verts[0][12],verts[0][13],verts[0][14],verts[0][15] },
+ { verts[2][12],verts[2][13],verts[2][14],verts[2][15] } };
+ MBEntityHandle nq[2] = { difference.front(), difference.back() };
+ for (int i = 0; i < 2; ++i) {
+ const MBEntityHandle* conn;
+ int len;
+ bool found = false;
+ for (int j = 0; !found && j < 2; ++j) {
+ rval = mb.get_connectivity( nq[j], conn, len );
+ if (MB_SUCCESS != rval) return rval;
+ int idx1 = std::find(conn,conn+len,expected_conn[i][0])-conn;
+ if (idx1 == len) continue;
+ found = true;
+ for (int k = 1; k < 4; ++k)
+ if (conn[(idx1+k)%4] != expected_conn[i][k])
+ found = false;
+ if (!found) {
+ found = true;
+ for (int k = 1; k < 4; ++k)
+ if (conn[(idx1+4-k)%4] != expected_conn[i][k])
+ found = false;
+ }
+ }
+ if (!found) {
+ std::cerr << "Skinner did not create & return expected quad " << i << std::endl;
+ return MB_FAILURE;
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
+// Test that skinning of higher-order elements works
+MBErrorCode mb_skin_higher_order_faces_common( bool use_adj )
+{
+ /* Create mesh:
+
+ 0---1---2---3---4
+ | | /
+ | | /
+ 5 6 7 8 9
+ | | /
+ | | /
+ 10--11--12
+ */
+
+ MBErrorCode rval;
+ MBCore moab;
+ MBInterface& mb = moab;
+
+ double coords[13][3] = {
+ {0,4,0}, {2,4,0}, {4,4,0}, {6,4,0}, {8,4,0},
+ {0,2,0}, {2,2,0}, {4,2,0}, {5,2,0}, {6,2,0},
+ {0,0,0}, {2,0,0}, {4,0,0} };
+ MBEntityHandle verts[13];
+ for (int i = 0; i < 13; ++i) {
+ rval = mb.create_vertex( coords[i], verts[i] );
+ if (MB_SUCCESS != rval) return rval;
+ }
+
+ MBEntityHandle qconn[9] = {
+ verts[0], verts[2], verts[12], verts[10],
+ verts[1], verts[7], verts[11], verts[5],
+ verts[6] };
+ MBEntityHandle tconn[7] = {
+ verts[2], verts[4], verts[12],
+ verts[3], verts[9], verts[7],
+ verts[8] };
+ MBEntityHandle quad, tri;
+ rval = mb.create_element( MBQUAD, qconn, 9, quad );
+ if (MB_SUCCESS != rval) return rval;
+ rval = mb.create_element( MBTRI, tconn, 7, tri );
+ if (MB_SUCCESS != rval) return rval;
+
+ MBRange faces;
+ faces.insert(quad);
+ faces.insert(tri);
+
+ MBRange skin_verts;
+ const int skin_vert_idx[] = { 0, 1, 2, 3, 4, 5, 9, 10, 11, 12 };
+ for (size_t i = 0; i < sizeof(skin_vert_idx)/sizeof(skin_vert_idx[0]); ++i)
+ skin_verts.insert( verts[skin_vert_idx[i]] );
+
+ MBSkinner tool(&mb);
+ MBRange skin;
+
+ rval = tool.find_skin( faces, true, skin, 0, use_adj, false );
+ if (MB_SUCCESS != rval) {
+ std::cout << "Vertex skinning failed with: " << mb.get_error_string(rval) << std::endl;
+ return rval;
+ }
+ if (skin != skin_verts) {
+ std::cout << "Skinner returned incorrect vertices." << std::endl;
+ return MB_FAILURE;
+ }
+
+ const int skin_edges[5][3] = {
+ {0,1,2}, {2,3,4}, {4,9,12}, {12,11,10}, {10,5,0} };
+ skin.clear();
+ rval = tool.find_skin( faces, false, skin, 0, use_adj, true );
+ if (MB_SUCCESS != rval) {
+ std::cout << "Edge skinning failed with: " << mb.get_error_string(rval) << std::endl;
+ return rval;
+ }
+ if (skin.size() != 5u) {
+ std::cout << "Skinner returned " << skin.size() << " vertices. Expected 5" << std::endl;
+ return MB_FAILURE;
+ }
+ int num_quadratic = 0;
+ const MBEntityHandle* conn;
+ int len;
+ for (MBRange::iterator i = skin.begin(); i != skin.end(); ++i) {
+ rval = mb.get_connectivity( *i, conn, len, false );
+ if (MB_SUCCESS != rval) return rval;
+ if (len == 3)
+ num_quadratic++;
+ else if(len != 2) {
+ std::cerr << "Skinner created edge with " << len << " vertices" << std::endl;
+ return MB_FAILURE;
+ }
+ }
+ if (num_quadratic != 5) {
+ std::cerr << num_quadratic << " of 5 created edges were quadratic" << std::endl;
+ return MB_FAILURE;
+ }
+
+ for (int i = 0; i < 5; ++i) {
+ bool found = false;
+ for (MBRange::iterator j = skin.begin(); j != skin.end(); ++j) {
+ mb.get_connectivity( *j, conn, len, false );
+ if (conn[2] == verts[skin_edges[i][1]]) {
+ found = true;
+ break;
+ }
+ }
+ if (!found) {
+ std::cerr << "One or more skin edges is incorrect" << std::endl;
+ return MB_FAILURE;
+ }
+ if ((conn[0] != verts[skin_edges[i][0]] || conn[1] != verts[skin_edges[i][2]])
+ && (conn[0] != verts[skin_edges[i][2]] || conn[1] != verts[skin_edges[i][0]])) {
+ std::cerr << "Invalid skin edge connectivity" << std::endl;
+ return MB_FAILURE;
+ }
+ }
+
+ return MB_SUCCESS;
+}
+MBErrorCode mb_skin_higher_order_faces_test( MBInterface* )
+ { return mb_skin_higher_order_faces_common( false ); }
+MBErrorCode mb_skin_adj_higher_order_faces_test( MBInterface* )
+ { return mb_skin_higher_order_faces_common( true ); }
+
+// Test that skinning of higher-order elements works
+MBErrorCode mb_skin_higher_order_regions_common( bool use_adj )
+{
+ // create mesh containing two 27-node hexes
+ /*
+ 0,2---1,2---2,2---3,2---4,2
+ | | |
+ | | |
+ 0,1 1,1 2,1 3,1 4,1
+ | | |
+ | | |
+ 0,0---1,0---2,0---3,0---4,0
+ */
+
+ MBErrorCode rval;
+ MBCore moab;
+ MBInterface& mb = moab;
+ MBRange hexes;
+
+
+ MBEntityHandle verts[5][3][3];
+ for (int i = 0; i < 5; ++i)
+ for (int j = 0; j < 3; ++j)
+ for (int k = 0; k < 3; ++k) {
+ double coords[] = { i, j, k };
+ rval = mb.create_vertex( coords, verts[i][j][k] );
+ if (MB_SUCCESS != rval) return rval;
+ }
+
+ int hex_conn[][3] = { // corners
+ {0,0,0}, {2,0,0}, {2,2,0}, {0,2,0},
+ {0,0,2}, {2,0,2}, {2,2,2}, {0,2,2},
+ // mid-edge
+ {1,0,0}, {2,1,0}, {1,2,0}, {0,1,0},
+ {0,0,1}, {2,0,1}, {2,2,1}, {0,2,1},
+ {1,0,2}, {2,1,2}, {1,2,2}, {0,1,2},
+ // mid-face
+ {1,0,1}, {2,1,1}, {1,2,1}, {0,1,1},
+ {1,1,0}, {1,1,2},
+ // mid-volume
+ {1,1,1} };
+
+ MBEntityHandle hexverts[2][27];
+ for (int i = 0; i < 2; ++i) {
+ MBEntityHandle h;
+ for (int j = 0; j < 27; ++j)
+ hexverts[i][j] = verts[ 2*i+hex_conn[j][0] ][ hex_conn[j][1] ][ hex_conn[j][2] ];
+ rval = mb.create_element( MBHEX, hexverts[i], 27, h );
+ if (MB_SUCCESS != rval)
+ return rval;
+ hexes.insert( h );
+ }
+
+ MBRange interior_verts;
+ interior_verts.insert( verts[1][1][1] ); // mid-node of hex 1
+ interior_verts.insert( verts[3][1][1] ); // mid-node of hex 2
+ interior_verts.insert( verts[2][1][1] ); // mid-node of shared face
+
+ MBSkinner tool(&mb);
+ MBRange skin;
+
+ rval = tool.find_skin( hexes, true, skin, 0, use_adj, false );
+ if (MB_SUCCESS != rval) {
+ std::cout << "Vertex skinning failed with: " << mb.get_error_string(rval)
+ << std::endl;
+ return rval;
+ }
+ MBRange extra = intersect( skin, interior_verts );
+ if (!extra.empty()) {
+ std::cout << "Skinner returned " << extra.size() << " interior vertices"
+ << std::endl;
+ std::cout << extra << std::endl;
+ return MB_FAILURE;
+ }
+ int num_vtx;
+ mb.get_number_entities_by_dimension( 0, 0, num_vtx );
+ size_t num_skin = num_vtx - interior_verts.size();
+ if (skin.size() != num_skin) {
+ std::cout << "Skinner returned " << skin.size() << " of "
+ << num_skin << " skin vertices" <<std::endl;
+ return MB_FAILURE;
+ }
+
+ skin.clear();
+ rval = tool.find_skin( hexes, false, skin, 0, use_adj, true );
+ if (MB_SUCCESS != rval) {
+ std::cout << "Element skinning failed with: " << mb.get_error_string(rval) << std::endl;
+ return rval;
+ }
+
+ if (skin.size() > 10u) {
+ std::cout << "Skinner created too many faces" << std::endl;
+ return MB_FAILURE;
+ }
+
+ bool all_okay = true;
+ bool faces[2][6] = { {0,0,0,0,0,0},{0,0,0,0,0,0} };
+ const MBEntityHandle *conn;
+ int len;
+ for (MBRange::iterator it = skin.begin(); it != skin.end(); ++it) {
+ rval = mb.get_connectivity( *it, conn, len );
+ if (MB_SUCCESS != rval) return rval;
+ if (len != 9) {
+ std::cout << "Skinner created face with " << len << " nodes" << std::endl;
+ all_okay = false;
+ continue;
+ }
+
+ int valid, side, sense, offset;
+ for (int h = 0; h < 2; ++h) {
+ valid = MBCN::SideNumber( MBHEX, hexverts[h], conn, 4, 2, side, sense, offset );
+ if (valid != 0)
+ continue;
+ if (sense != 1) {
+ std::cout << "Skinner created reversed face for hex "
+ << h << " side " << side << std::endl;
+ all_okay = false;
+ continue;
+ }
+
+ int idx[9], len2;
+ MBEntityType sidetype;
+ MBCN::SubEntityNodeIndices( MBHEX, 27, 2, side, sidetype, len2, idx );
+ assert(sidetype == MBQUAD);
+ assert(len2 == 9);
+ if ( conn[ offset ] != hexverts[h][idx[0]] ||
+ conn[ (offset+1)%4] != hexverts[h][idx[1]] ||
+ conn[ (offset+2)%4] != hexverts[h][idx[2]] ||
+ conn[ (offset+3)%4] != hexverts[h][idx[3]] ||
+ conn[4+ offset ] != hexverts[h][idx[4]] ||
+ conn[4+(offset+1)%4] != hexverts[h][idx[5]] ||
+ conn[4+(offset+2)%4] != hexverts[h][idx[6]] ||
+ conn[4+(offset+3)%4] != hexverts[h][idx[7]] ||
+ conn[ 8 ] != hexverts[h][idx[8]]) {
+ std::cout << "Side " << side << " of hex " << h
+ << " has invalid connectivity" << std::endl;
+ all_okay = false;
+ }
+
+
+ faces[h][side] = true;
+ }
+ }
+
+ for (int h = 0; h < 2; ++h) {
+ for (int i = 0; i < 6; ++i) {
+ if ((h == 0 && i == 1) || (h == 1 && i == 3)) {
+ if (faces[h][i]) {
+ std::cout << "Skinner created interior face for side "
+ << i << " of hex " << h << std::endl;
+ all_okay = false;
+ }
+ }
+ else if (!faces[h][i]) {
+ std::cout << "Skinner failed to create side "
+ << i << " of hex " << h << std::endl;
+ all_okay = false;
+ }
+ }
+ }
+
+ return all_okay ? MB_SUCCESS : MB_FAILURE;
+}
+
+MBErrorCode mb_skin_higher_order_regions_test( MBInterface* )
+ { return mb_skin_higher_order_regions_common(false); }
+MBErrorCode mb_skin_adj_higher_order_regions_test( MBInterface* )
+ { return mb_skin_higher_order_regions_common(true); }
+
+
+MBErrorCode mb_skin_reversed_common( int dim, bool use_adj )
+{
+ MBEntityType type, subtype;
+ switch (dim) {
+ case 2: type = MBTRI; subtype = MBEDGE; break;
+ case 3: type = MBTET; subtype = MBTRI; break;
+ default: assert(false); return MB_FAILURE;
+ }
+
+ /* 3
+ /|\
+ / | \
+ / | \
+ / | \
+ 0----1----2
+ */
+
+ MBErrorCode rval;
+ MBCore moab;
+ MBInterface& mb = moab;
+ MBRange hexes;
+
+ double coords[][3] = { { 0, 0, 0 },
+ { 1, 0, 0 },
+ { 2, 0, 0 },
+ { 1, 2, 0 },
+ { 1, 2, 2 } };
+ MBEntityHandle verts[5];
+ const int nvtx = 2+dim;
+ for (int i = 0; i < nvtx; ++i) {
+ rval = mb.create_vertex( coords[i], verts[i] );
+ if (MB_SUCCESS != rval) return rval;
+ }
+ // NOTE: order connectivity such that side 1 is shared!
+ MBEntityHandle conn[2][4] = {
+ { verts[0], verts[1], verts[3], verts[4] },
+ { verts[2], verts[3], verts[1], verts[4] } };
+ const int conn_len = dim+1;
+ MBRange elems;
+ for (int i = 0; i < 2; ++i) {
+ MBEntityHandle h;
+ rval = mb.create_element( type, conn[i], conn_len, h );
+ if (MB_SUCCESS != rval) return rval;
+ elems.insert(h);
+ }
+
+ // create one reversed side
+ MBEntityHandle side_conn[3];
+ int side_indices[3] = {0,0,0};
+ MBCN::SubEntityVertexIndices(type, dim-1, 0, side_indices );
+ side_conn[0] = conn[0][side_indices[1]];
+ side_conn[1] = conn[0][side_indices[0]];
+ side_conn[2] = conn[0][side_indices[2]];
+ MBEntityHandle side;
+ rval = mb.create_element( subtype, side_conn, dim, side );
+ if (MB_SUCCESS != rval) return rval;
+
+ MBRange forward, reverse;
+ MBSkinner tool(&mb);
+ rval = tool.find_skin( elems, false, forward, &reverse, use_adj, true );
+ if (MB_SUCCESS != rval) {
+ std::cout << "Skinner failed." << std::endl;
+ return rval;
+ }
+
+ // expect all faces created by the skinner to be forward,
+ // so the only reversed side should be the one created above.
+ if (reverse.size() != 1 || reverse.front() != side) {
+ std::cout << "Expected 1 reversed side, got: " << reverse.size() << std::endl;
+ return MB_FAILURE;
+ }
+
+ return MB_SUCCESS;
+}
+MBErrorCode mb_skin_faces_reversed_test( MBInterface* )
+ { return mb_skin_reversed_common( 2, false ); }
+MBErrorCode mb_skin_adj_faces_reversed_test( MBInterface* )
+ { return mb_skin_reversed_common( 2, true ); }
+MBErrorCode mb_skin_regions_reversed_test( MBInterface* )
+ { return mb_skin_reversed_common( 3, false ); }
+MBErrorCode mb_skin_adj_regions_reversed_test( MBInterface* )
+ { return mb_skin_reversed_common( 3, true ); }
+
static void usage(const char* exe) {
cerr << "Usage: " << exe << " [-nostress] [-d input_file_dir]\n";
exit (1);
@@ -6972,6 +7741,19 @@
RUN_TEST( mb_skin_surface_adj_test );
RUN_TEST( mb_skin_volume_test );
RUN_TEST( mb_skin_volume_adj_test );
+ RUN_TEST( mb_skin_surf_verts_test );
+ RUN_TEST( mb_skin_vol_verts_test );
+ RUN_TEST( mb_skin_surf_verts_elems_test );
+ RUN_TEST( mb_skin_vol_verts_elems_test );
+ RUN_TEST( mb_skin_poly_test );
+ RUN_TEST( mb_skin_higher_order_faces_test );
+ RUN_TEST( mb_skin_higher_order_regions_test );
+ RUN_TEST( mb_skin_adj_higher_order_faces_test );
+ RUN_TEST( mb_skin_adj_higher_order_regions_test );
+ RUN_TEST( mb_skin_faces_reversed_test );
+ RUN_TEST( mb_skin_adj_faces_reversed_test );
+ RUN_TEST( mb_skin_regions_reversed_test );
+ RUN_TEST( mb_skin_adj_regions_reversed_test );
RUN_TEST( mb_read_fail_test );
RUN_TEST( mb_enum_string_test );
RUN_TEST( mb_merge_update_test );
Modified: MOAB/trunk/tools/mbperf/mbperf.cpp
===================================================================
--- MOAB/trunk/tools/mbperf/mbperf.cpp 2009-11-19 23:08:51 UTC (rev 3372)
+++ MOAB/trunk/tools/mbperf/mbperf.cpp 2009-11-19 23:10:26 UTC (rev 3373)
@@ -51,12 +51,9 @@
const int DEFAULT_INTERVALS = 50;
-void testA(const int nelem);
-void testB(const int nelem);
-void testC(const int nelem);
-void skin_time(const int intervals);
-void skin_adj(const int intervals);
-void skin_vert(const int intervals);
+void testA(const int nelem, const int );
+void testB(const int nelem, const int );
+void testC(const int nelem, const int );
void print_time(const bool print_em, double &tot_time, double &utime, double &stime,
double &mem);
void query_vert_to_elem(MBInterface*);
@@ -65,6 +62,14 @@
void get_time_mem(double &tot_time, double &user_time,
double &sys_time, double &tot_mem);
void bulk_construct_mesh( MBInterface* gMB, const int nelem );
+void create_regular_mesh( int interval, int dimension );
+void skin_common( int interval, int dim, bool use_adj );
+void skin(const int intervals, const int dim)
+ { std::cout << "Skinning w/out adjacencies:" << std::endl;
+ skin_common( intervals, dim, false ); }
+void skin_adj(const int intervals, const int dim)
+ { std::cout << "Skinning with adjacencies:" << std::endl;
+ skin_common( intervals, dim, true ); }
void compute_edge(double *start, const int nelem, const double xint,
const int stride)
@@ -280,7 +285,7 @@
}
}
-typedef void (*test_func_t)( int );
+typedef void (*test_func_t)( const int, const int );
const struct {
std::string testName;
test_func_t testFunc;
@@ -289,20 +294,20 @@
{ "struct", &testA, "Conn. and adj. query time for structured mesh" },
{ "bulk", &testB, "Conn. and adj. query time for bulk-created mesh" },
{ "indiv", &testC, "Conn. and adj. query time for per-entity created mesh" },
- { "skin", &skin_time, "Test time to get skin quads" },
- { "skin_adj", &skin_adj, "Test skin quads using vert-to-elem adjacencies" },
- { "skin_vert", &skin_vert, "Test time to get skin verts" },
+ { "skin", &skin, "Test time to get skin mesh w/out adjacencies" },
+ { "skin_adj", &skin_adj, "Test time to get skin mesh with adjacencies" },
};
const int TestListSize = sizeof(TestList)/sizeof(TestList[0]);
void usage( const char* argv0, bool error = true )
{
std::ostream& str = error ? std::cerr : std::cout;
- str << "Usage: " << argv0 << " -i <ints_per_side> <test_name> [<test_name2> ...]" << std::endl;
+ str << "Usage: " << argv0 << " -i <ints_per_side> -d <dimension> <test_name> [<test_name2> ...]" << std::endl;
str << " " << argv0 << " [-h|-l]" << std::endl;
if (error)
return;
str << " -i : specify interverals per side (num hex = ints^3, default: " << DEFAULT_INTERVALS << std::endl;
+ str << " -d : specify element dimension (INGORED BY SOME TESTS), default: 3" << std::endl;
str << " -h : print this help text." << std::endl;
str << " -l : list available tests" << std::endl;
}
@@ -329,9 +334,10 @@
int main(int argc, char* argv[])
{
int intervals = DEFAULT_INTERVALS;
+ int dimension = 3;
std::vector<test_func_t> test_list;
bool did_help = false;
- bool expect_ints = false;
+ bool expect_ints = false, expect_dim = false;
for (int i = 1; i < argc; ++i) {
if (expect_ints) {
expect_ints = false;
@@ -343,10 +349,21 @@
return 1;
}
}
+ else if (expect_dim) {
+ expect_dim = false;
+ char* endptr;
+ dimension = strtol( argv[i], &endptr, 0 );
+ if (dimension < 1 || dimension > 3) {
+ usage(argv[0]);
+ std::cerr << "Invalid dimension: " << dimension << std::endl;
+ return 1;
+ }
+ }
else if (*argv[i] == '-') { // flag
for (int j = 1; argv[i][j]; ++j) {
switch (argv[i][j]) {
case 'i': expect_ints = true; break;
+ case 'd': expect_dim = true; break;
case 'h': did_help = true; usage(argv[0],false); break;
case 'l': did_help = true; list_tests(); break;
default:
@@ -380,7 +397,7 @@
// now run the tests
for (std::vector<test_func_t>::iterator i = test_list.begin(); i != test_list.end(); ++i) {
test_func_t fptr = *i;
- fptr(intervals);
+ fptr(intervals,dimension);
}
return 0;
@@ -512,7 +529,7 @@
}
#endif
-void testA( int nelem )
+void testA( const int nelem, const int )
{
MBCore moab;
MBInterface* gMB = &moab;
@@ -616,7 +633,7 @@
assert(MB_SUCCESS == result);
}
-void testB(const int nelem)
+void testB(const int nelem, const int )
{
MBCore moab;
MBInterface* gMB = &moab;
@@ -649,7 +666,7 @@
<< std::endl;
}
-void testC(const int nelem)
+void testC(const int nelem, const int )
{
MBCore moab;
MBInterface* gMB = &moab;
@@ -722,68 +739,166 @@
<< std::endl;
}
-void skin_time( const int nelem )
+
+void create_regular_mesh( MBInterface* gMB, int interval, int dim )
{
- MBCore moab;
- MBInterface* gMB = &moab;
- MBErrorCode rval;
+ if (dim < 1 || dim > 3 || interval < 1) {
+ std::cerr << "Invalid arguments" << std::endl;
+ exit(1);
+ }
+
+ const int nvi = interval+1;
+ const int dims[3] = { nvi, dim > 1 ? nvi : 1, dim > 2 ? nvi : 1 };
+ int num_vert = dims[0] * dims[1] * dims[2];
- bulk_construct_mesh( gMB, nelem );
-
- MBRange skin, hexes;
- rval = gMB->get_entities_by_dimension( 0, 3, hexes );
- assert(MB_SUCCESS == rval); assert(!hexes.empty());
-
- MBSkinner tool(gMB);
- clock_t t = clock();
- rval = tool.find_skin( hexes, 2, skin, false );
- double d = ((double)(clock() - t))/CLOCKS_PER_SEC;
- assert(MB_SUCCESS == rval);
+ void *ptr = 0;
+ // get the read interface
+ gMB->query_interface("MBReadUtilIface", &ptr);
+ MBReadUtilIface* readMeshIface = static_cast<MBReadUtilIface*>(ptr);
- std::cout << "Got 2D skin in " << d << " seconds w/out vertex-to-element adjacencies" << std::endl;
-}
+ MBEntityHandle vstart;
+ std::vector<double*> arrays;
+ MBErrorCode rval = readMeshIface->get_node_arrays(3, num_vert, 1, vstart, arrays);
+ if (MB_SUCCESS != rval || arrays.size() < 3) {
+ std::cerr << "Vertex creation failed" << std::endl;
+ exit(2);
+ }
+ double *x = arrays[0], *y = arrays[1], *z = arrays[2];
+
+ // Calculate vertex coordinates
+ for (int k = 0; k < dims[2]; ++k)
+ for (int j = 0; j < dims[1]; ++j)
+ for (int i = 0; i < dims[0]; ++i)
+ {
+ *x = i; ++x;
+ *y = j; ++y;
+ *z = k; ++z;
+ }
+
+ const long vert_per_elem = 1 << dim; // 2^dim
+ const long intervals[3] = { interval, dim>1?interval:1, dim>2?interval:1 };
+ const long num_elem = intervals[0]*intervals[1]*intervals[2];
+ const MBEntityType type = (dim == 1) ? MBEDGE : (dim == 2) ? MBQUAD : MBHEX;
+
+ MBEntityHandle estart, *conn = 0;
+ rval = readMeshIface->get_element_array( num_elem, vert_per_elem, type, 0, estart, conn );
+ if (MB_SUCCESS != rval || !conn) {
+ std::cerr << "Element creation failed" << std::endl;
+ exit(2);
+ }
+
+
+ // Offsets of element vertices in grid relative to corner closest to origin
+ long c = dims[0]*dims[1];
+ const long corners[8] = { 0, 1, 1+dims[0], dims[0], c, c+1, c+1+dims[0], c+dims[0] };
+
+ // Populate element list
+ MBEntityHandle* iter = conn;
+ for (long z = 0; z < intervals[2]; ++z)
+ for (long y = 0; y < intervals[1]; ++y)
+ for (long x = 0; x < intervals[0]; ++x)
+ {
+ const long index = x + y*dims[0] + z*(dims[0]*dims[1]);
+ for (long j = 0; j < vert_per_elem; ++j, ++iter)
+ *iter = index + corners[j] + vstart;
+ }
+
+ // notify MOAB of the new elements
+ rval = readMeshIface->update_adjacencies(estart, num_elem, vert_per_elem, conn);
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Element update failed" << std::endl;
+ exit(2);
+ }
+}
-void skin_adj( const int nelem )
+void skin_common( int interval, int dim, bool use_adj )
{
MBCore moab;
MBInterface* gMB = &moab;
MBErrorCode rval;
+ double d;
+ clock_t t, tt;
- bulk_construct_mesh( gMB, nelem );
+ create_regular_mesh( gMB, interval, dim );
- MBRange skin, verts, hexes;
- // force creation of adjacencies
- rval = gMB->get_entities_by_dimension( 0, 0, verts );
- assert(MB_SUCCESS == rval); assert(!verts.empty());
- rval = gMB->get_adjacencies( verts, 3, false, hexes, MBInterface::UNION );
- assert(MB_SUCCESS == rval); assert(!hexes.empty());
+ MBRange skin, verts, elems;
+ rval = gMB->get_entities_by_dimension( 0, dim, elems );
+ assert(MB_SUCCESS == rval); assert(!elems.empty());
MBSkinner tool(gMB);
- clock_t t = clock();
- rval = tool.find_skin( hexes, 2, skin, true );
- double d = ((double)(clock() - t))/CLOCKS_PER_SEC;
- assert(MB_SUCCESS == rval);
+
+ t = clock();
+ rval = tool.find_skin( elems, true, verts, 0, use_adj, false );
+ t = clock() - t;
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Search for skin vertices failed" << std::endl;
+ exit(2);
+ }
+ d = ((double)t)/CLOCKS_PER_SEC;
+ std::cout << "Got " << verts.size() << " skin vertices in " << d << " seconds." << std::endl;
+
+ t = 0;
+ long blocksize = elems.size() / 1000;
+ if (!blocksize) blocksize = 1;
+ long numblocks = elems.size()/blocksize;
+ MBRange::iterator it = elems.begin();
+ for (long i = 0; i < numblocks; ++i) {
+ verts.clear();
+ MBRange::iterator end = it + blocksize;
+ MBRange blockelems;
+ blockelems.merge( it, end );
+ it = end;
+ tt = clock();
+ rval = tool.find_skin( blockelems, true, verts, 0, use_adj, false );
+ t += clock() - tt;
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Search for skin vertices failed" << std::endl;
+ exit(2);
+ }
+ }
+ d = ((double)t)/CLOCKS_PER_SEC;
+ std::cout << "Got skin vertices for " << numblocks << " blocks of "
+ << blocksize << " elements in " << d << " seconds." << std::endl;
+
+ for (int e = 0; e < 2; ++e) { // do this twice
+ if (e == 1) {
+ // create all interior faces
+ skin.clear();
+ t = clock();
+ gMB->get_adjacencies( elems, dim-1, true, skin, MBInterface::UNION );
+ t = clock() - t;
+ d = ((double)t)/CLOCKS_PER_SEC;
+ std::cout << "Created " << skin.size() << " entities of dimension-1 in " << d << " seconds" << std::endl;
+ }
+
+ t = clock();
+ rval = tool.find_skin( elems, false, skin, 0, use_adj, true );
+ t = clock() - t;
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Search for skin vertices failed" << std::endl;
+ exit(2);
+ }
+ d = ((double)t)/CLOCKS_PER_SEC;
+ std::cout << "Got " << skin.size() << " skin elements in " << d << " seconds." << std::endl;
- std::cout << "Got 2D skin in " << d << " seconds with vertex-to-element adjacencies" << std::endl;
+ t = 0;
+ it = elems.begin();
+ for (long i = 0; i < numblocks; ++i) {
+ skin.clear();
+ MBRange::iterator end = it + blocksize;
+ MBRange blockelems;
+ blockelems.merge( it, end );
+ it = end;
+ tt = clock();
+ rval = tool.find_skin( blockelems, false, skin, 0, use_adj, true );
+ t += clock() - tt;
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Search for skin elements failed" << std::endl;
+ exit(2);
+ }
+ }
+ d = ((double)t)/CLOCKS_PER_SEC;
+ std::cout << "Got skin elements for " << numblocks << " blocks of "
+ << blocksize << " elements in " << d << " seconds." << std::endl;
+ }
}
-
-void skin_vert( const int nelem )
-{
- MBCore moab;
- MBInterface* gMB = &moab;
- MBErrorCode rval;
-
- bulk_construct_mesh( gMB, nelem );
-
- MBRange skin, hexes;
- rval = gMB->get_entities_by_dimension( 0, 3, hexes );
- assert(MB_SUCCESS == rval); assert(!hexes.empty());
-
- MBSkinner tool(gMB);
- clock_t t = clock();
- rval = tool.find_skin( hexes, 0, skin, false );
- double d = ((double)(clock() - t))/CLOCKS_PER_SEC;
- assert(MB_SUCCESS == rval);
-
- std::cout << "Got 0D skin in " << d << " seconds w/out vertex-to-element adjacencies" << std::endl;
-}
Modified: MOAB/trunk/tools/skin.cpp
===================================================================
--- MOAB/trunk/tools/skin.cpp 2009-11-19 23:08:51 UTC (rev 3372)
+++ MOAB/trunk/tools/skin.cpp 2009-11-19 23:10:26 UTC (rev 3373)
@@ -189,7 +189,7 @@
// skin the mesh
MBRange forward_lower, reverse_lower;
MBSkinner tool( iface );
- result = tool.find_skin( skin_ents, forward_lower, reverse_lower );
+ result = tool.find_skin( skin_ents, false, forward_lower, &reverse_lower );
MBRange boundary;
boundary.merge( forward_lower );
boundary.merge( reverse_lower );
More information about the moab-dev
mailing list