Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
Skinner.cpp
Go to the documentation of this file.
1 /**
2  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3  * storing and accessing finite element mesh data.
4  *
5  * Copyright 2004 Sandia Corporation. Under the terms of Contract
6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7  * retains certain rights in this software.
8  *
9  * This library is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 2.1 of the License, or (at your option) any later version.
13  *
14  */
15 
16 #ifdef __MFC_VER
17 #pragma warning( disable : 4786 )
18 #endif
19 
20 #ifdef WIN32 /* windows */
21 #define _USE_MATH_DEFINES // For M_PI
22 #endif
23 #include "moab/Skinner.hpp"
24 #include "moab/Range.hpp"
25 #include "moab/CN.hpp"
26 #include <vector>
27 #include <set>
28 #include <algorithm>
29 #include <cmath>
30 #include <cassert>
31 #include <iostream>
32 #include "moab/Util.hpp"
33 #include "Internals.hpp"
34 #include "MBTagConventions.hpp"
35 #include "moab/Core.hpp"
36 #include "AEntityFactory.hpp"
37 #include "moab/ScdInterface.hpp"
38 
39 #ifdef M_PI
40 #define SKINNER_PI M_PI
41 #else
42 #define SKINNER_PI 3.1415926535897932384626
43 #endif
44 
45 namespace moab
46 {
47 
49 {
50  // delete the adjacency tag
51 }
52 
54 {
55  // go through and mark all the target dimension entities
56  // that already exist as not deleteable
57  // also get the connectivity tags for each type
58  // also populate adjacency information
59  EntityType type;
60  DimensionPair target_ent_types = CN::TypeDimensionMap[mTargetDim];
61 
62  void* null_ptr = NULL;
63 
64  ErrorCode result = thisMB->tag_get_handle( "skinner adj", sizeof( void* ), MB_TYPE_OPAQUE, mAdjTag,
65  MB_TAG_DENSE | MB_TAG_CREAT, &null_ptr );MB_CHK_ERR( result );
66 
67  if( mDeletableMBTag == 0 )
68  {
69  result =
70  thisMB->tag_get_handle( "skinner deletable", 1, MB_TYPE_BIT, mDeletableMBTag, MB_TAG_BIT | MB_TAG_CREAT );MB_CHK_ERR( result );
71  }
72 
74 
75  // go through each type at this dimension
76  for( type = target_ent_types.first; type <= target_ent_types.second; ++type )
77  {
78  // get the entities of this type in the MB
80 
81  // go through each entity of this type in the MB
82  // and set its deletable tag to NO
83  Range::iterator iter, end_iter;
84  end_iter = entities.end();
85  for( iter = entities.begin(); iter != end_iter; ++iter )
86  {
87  unsigned char bit = 0x1;
88  result = thisMB->tag_set_data( mDeletableMBTag, &( *iter ), 1, &bit );
89  assert( MB_SUCCESS == result );
90  // add adjacency information too
91  if( TYPE_FROM_HANDLE( *iter ) != MBVERTEX ) add_adjacency( *iter );
92  }
93  }
94 
95  return MB_SUCCESS;
96 }
97 
99 {
100  ErrorCode result;
101  if( 0 != mDeletableMBTag )
102  {
103  result = thisMB->tag_delete( mDeletableMBTag );
104  mDeletableMBTag = 0;MB_CHK_ERR( result );
105  }
106 
107  // remove the adjacency tag
108  std::vector< std::vector< EntityHandle >* > adj_arr;
109  std::vector< std::vector< EntityHandle >* >::iterator i;
110  if( 0 != mAdjTag )
111  {
112  for( EntityType t = MBVERTEX; t != MBMAXTYPE; ++t )
113  {
114  Range entities;
115  result = thisMB->get_entities_by_type_and_tag( 0, t, &mAdjTag, 0, 1, entities );MB_CHK_ERR( result );
116  adj_arr.resize( entities.size() );
117  result = thisMB->tag_get_data( mAdjTag, entities, &adj_arr[0] );MB_CHK_ERR( result );
118  for( i = adj_arr.begin(); i != adj_arr.end(); ++i )
119  delete *i;
120  }
121 
122  result = thisMB->tag_delete( mAdjTag );
123  mAdjTag = 0;MB_CHK_ERR( result );
124  }
125 
126  return MB_SUCCESS;
127 }
128 
130 {
131  std::vector< EntityHandle >* adj = NULL;
132  const EntityHandle* nodes;
133  int num_nodes;
134  ErrorCode result = thisMB->get_connectivity( entity, nodes, num_nodes, true );MB_CHK_ERR( result );
135  const EntityHandle* iter = std::min_element( nodes, nodes + num_nodes );
136 
137  if( iter == nodes + num_nodes ) return MB_SUCCESS;
138 
139  // add this entity to the node
140  if( thisMB->tag_get_data( mAdjTag, iter, 1, &adj ) == MB_SUCCESS && adj != NULL )
141  {
142  adj->push_back( entity );
143  }
144  // create a new vector and add it
145  else
146  {
147  adj = new std::vector< EntityHandle >;
148  adj->push_back( entity );
149  result = thisMB->tag_set_data( mAdjTag, iter, 1, &adj );MB_CHK_ERR( result );
150  }
151 
152  return MB_SUCCESS;
153 }
154 
155 void Skinner::add_adjacency( EntityHandle entity, const EntityHandle* nodes, const int num_nodes )
156 {
157  std::vector< EntityHandle >* adj = NULL;
158  const EntityHandle* iter = std::min_element( nodes, nodes + num_nodes );
159 
160  if( iter == nodes + num_nodes ) return;
161 
162  // should not be setting adjacency lists in ho-nodes
163  assert( TYPE_FROM_HANDLE( entity ) == MBPOLYGON ||
164  num_nodes == CN::VerticesPerEntity( TYPE_FROM_HANDLE( entity ) ) );
165 
166  // add this entity to the node
167  if( thisMB->tag_get_data( mAdjTag, iter, 1, &adj ) == MB_SUCCESS && adj != NULL )
168  {
169  adj->push_back( entity );
170  }
171  // create a new vector and add it
172  else
173  {
174  adj = new std::vector< EntityHandle >;
175  adj->push_back( entity );
176  thisMB->tag_set_data( mAdjTag, iter, 1, &adj );
177  }
178 }
179 
180 ErrorCode Skinner::find_geometric_skin( const EntityHandle meshset, Range& forward_target_entities )
181 {
182  // attempts to find whole model skin, using geom topo sets first then
183  // normal find_skin function
184  bool debug = true;
185 
186  // look for geom topo sets
187  Tag geom_tag;
188  ErrorCode result =
190 
191  if( MB_SUCCESS != result ) return result;
192 
193  // get face sets (dimension = 2)
194  Range face_sets;
195  int two = 2;
196  const void* two_ptr = &two;
197  result = thisMB->get_entities_by_type_and_tag( meshset, MBENTITYSET, &geom_tag, &two_ptr, 1, face_sets );
198 
199  Range::iterator it;
200  if( MB_SUCCESS != result )
201  return result;
202  else if( face_sets.empty() )
203  return MB_ENTITY_NOT_FOUND;
204 
205  // ok, we have face sets; use those to determine skin
206  Range skin_sets;
207  if( debug ) std::cout << "Found " << face_sets.size() << " face sets total..." << std::endl;
208 
209  for( it = face_sets.begin(); it != face_sets.end(); ++it )
210  {
211  int num_parents;
212  result = thisMB->num_parent_meshsets( *it, &num_parents );
213  if( MB_SUCCESS != result )
214  return result;
215  else if( num_parents == 1 )
216  skin_sets.insert( *it );
217  }
218 
219  if( debug ) std::cout << "Found " << skin_sets.size() << " 1-parent face sets..." << std::endl;
220 
221  if( skin_sets.empty() ) return MB_FAILURE;
222 
223  // ok, we have the shell; gather up the elements, putting them all in forward for now
224  for( it = skin_sets.begin(); it != skin_sets.end(); ++it )
225  {
226  result = thisMB->get_entities_by_handle( *it, forward_target_entities, true );
227  if( MB_SUCCESS != result ) return result;
228  }
229 
230  return result;
231 }
232 
234  const Range& source_entities,
235  bool get_vertices,
236  Range& output_handles,
237  Range* output_reverse_handles,
238  bool create_vert_elem_adjs,
239  bool create_skin_elements,
240  bool look_for_scd )
241 {
242  if( source_entities.empty() ) return MB_SUCCESS;
243 
244  if( look_for_scd )
245  {
246  ErrorCode rval = find_skin_scd( source_entities, get_vertices, output_handles, create_skin_elements );
247  // if it returns success, it's all scd, and we don't need to do anything more
248  if( MB_SUCCESS == rval ) return rval;
249  }
250 
251  Core* this_core = dynamic_cast< Core* >( thisMB );
252  if( this_core && create_vert_elem_adjs && !this_core->a_entity_factory()->vert_elem_adjacencies() )
254 
255  return find_skin_vertices( meshset, source_entities, get_vertices ? &output_handles : 0,
256  get_vertices ? 0 : &output_handles, output_reverse_handles, create_skin_elements );
257 }
258 
259 ErrorCode Skinner::find_skin_scd( const Range& source_entities,
260  bool get_vertices,
261  Range& output_handles,
262  bool create_skin_elements )
263 {
264  // get the scd interface and check if it's been initialized
265  ScdInterface* scdi = NULL;
266  ErrorCode rval = thisMB->query_interface( scdi );
267  if( !scdi ) return MB_FAILURE;
268 
269  // ok, there's scd mesh; see if the entities passed in are all in a scd box
270  // a box needs to be wholly included in entities for this to work
271  std::vector< ScdBox* > boxes, myboxes;
272  Range myrange;
273  rval = scdi->find_boxes( boxes );
274  if( MB_SUCCESS != rval ) return rval;
275  for( std::vector< ScdBox* >::iterator bit = boxes.begin(); bit != boxes.end(); ++bit )
276  {
277  Range belems( ( *bit )->start_element(), ( *bit )->start_element() + ( *bit )->num_elements() - 1 );
278  if( source_entities.contains( belems ) )
279  {
280  myboxes.push_back( *bit );
281  myrange.merge( belems );
282  }
283  }
284  if( myboxes.empty() || myrange.size() != source_entities.size() ) return MB_FAILURE;
285 
286  // ok, we're all structured; get the skin for each box
287  for( std::vector< ScdBox* >::iterator bit = boxes.begin(); bit != boxes.end(); ++bit )
288  {
289  rval = skin_box( *bit, get_vertices, output_handles, create_skin_elements );
290  if( MB_SUCCESS != rval ) return rval;
291  }
292 
293  return MB_SUCCESS;
294 }
295 
296 ErrorCode Skinner::skin_box( ScdBox* box, bool get_vertices, Range& output_handles, bool create_skin_elements )
297 {
298  HomCoord bmin = box->box_min(), bmax = box->box_max();
299 
300  // don't support 1d boxes
301  if( bmin.j() == bmax.j() && bmin.k() == bmax.k() ) return MB_FAILURE;
302 
303  int dim = ( bmin.k() == bmax.k() ? 1 : 2 );
304 
305  ErrorCode rval;
306  EntityHandle ent;
307 
308  // i=min
309  for( int k = bmin.k(); k < bmax.k(); k++ )
310  {
311  for( int j = bmin.j(); j < bmax.j(); j++ )
312  {
313  ent = 0;
314  rval = box->get_adj_edge_or_face( dim, bmin.i(), j, k, 0, ent, create_skin_elements );
315  if( MB_SUCCESS != rval ) return rval;
316  if( ent ) output_handles.insert( ent );
317  }
318  }
319  // i=max
320  for( int k = bmin.k(); k < bmax.k(); k++ )
321  {
322  for( int j = bmin.j(); j < bmax.j(); j++ )
323  {
324  ent = 0;
325  rval = box->get_adj_edge_or_face( dim, bmax.i(), j, k, 0, ent, create_skin_elements );
326  if( MB_SUCCESS != rval ) return rval;
327  if( ent ) output_handles.insert( ent );
328  }
329  }
330  // j=min
331  for( int k = bmin.k(); k < bmax.k(); k++ )
332  {
333  for( int i = bmin.i(); i < bmax.i(); i++ )
334  {
335  ent = 0;
336  rval = box->get_adj_edge_or_face( dim, i, bmin.j(), k, 1, ent, create_skin_elements );
337  if( MB_SUCCESS != rval ) return rval;
338  if( ent ) output_handles.insert( ent );
339  }
340  }
341  // j=max
342  for( int k = bmin.k(); k < bmax.k(); k++ )
343  {
344  for( int i = bmin.i(); i < bmax.i(); i++ )
345  {
346  ent = 0;
347  rval = box->get_adj_edge_or_face( dim, i, bmax.j(), k, 1, ent, create_skin_elements );
348  if( MB_SUCCESS != rval ) return rval;
349  if( ent ) output_handles.insert( ent );
350  }
351  }
352  // k=min
353  for( int j = bmin.j(); j < bmax.j(); j++ )
354  {
355  for( int i = bmin.i(); i < bmax.i(); i++ )
356  {
357  ent = 0;
358  rval = box->get_adj_edge_or_face( dim, i, j, bmin.k(), 2, ent, create_skin_elements );
359  if( MB_SUCCESS != rval ) return rval;
360  if( ent ) output_handles.insert( ent );
361  }
362  }
363  // k=max
364  for( int j = bmin.j(); j < bmax.j(); j++ )
365  {
366  for( int i = bmin.i(); i < bmax.i(); i++ )
367  {
368  ent = 0;
369  rval = box->get_adj_edge_or_face( dim, i, j, bmax.k(), 2, ent, create_skin_elements );
370  if( MB_SUCCESS != rval ) return rval;
371  if( ent ) output_handles.insert( ent );
372  }
373  }
374 
375  if( get_vertices )
376  {
377  Range verts;
378  rval = thisMB->get_adjacencies( output_handles, 0, true, verts, Interface::UNION );
379  if( MB_SUCCESS != rval ) return rval;
380  output_handles.merge( verts );
381  }
382 
383  return MB_SUCCESS;
384 }
385 
386 ErrorCode Skinner::find_skin_noadj(const Range &source_entities,
387  Range &forward_target_entities,
388  Range &reverse_target_entities/*,
389  bool create_vert_elem_adjs*/)
390 {
391  if( source_entities.empty() ) return MB_FAILURE;
392 
393  // get our working dimensions
394  EntityType type = thisMB->type_from_handle( *( source_entities.begin() ) );
395  const int source_dim = CN::Dimension( type );
396  mTargetDim = source_dim - 1;
397 
398  // make sure we can handle the working dimensions
399  if( mTargetDim < 0 || source_dim > 3 ) return MB_FAILURE;
400 
401  initialize();
402 
403  Range::const_iterator iter, end_iter;
404  end_iter = source_entities.end();
405  const EntityHandle* conn;
406  EntityHandle match;
407 
408  direction direct;
409  ErrorCode result;
410  // assume we'll never have more than 32 vertices on a facet (checked
411  // with assert later)
412  EntityHandle sub_conn[32];
413  std::vector< EntityHandle > tmp_conn_vec;
414  int num_nodes, num_sub_nodes, num_sides;
415  int sub_indices[32]; // Also, assume that no polygon has more than 32 nodes
416  // we could increase that, but we will not display it right in visit moab h5m , anyway
417  EntityType sub_type;
418 
419  // for each source entity
420  for( iter = source_entities.begin(); iter != end_iter; ++iter )
421  {
422  // get the connectivity of this entity
423  int actual_num_nodes_polygon = 0;
424  result = thisMB->get_connectivity( *iter, conn, num_nodes, false, &tmp_conn_vec );
425  if( MB_SUCCESS != result ) return result;
426 
427  type = thisMB->type_from_handle( *iter );
428  Range::iterator seek_iter;
429 
430  // treat separately polygons (also, polyhedra will need special handling)
431  if( MBPOLYGON == type )
432  {
433  // treat padded polygons, if existing; count backwards, see how many of the last nodes
434  // are repeated assume connectivity is fine, otherwise we could be in trouble
435  actual_num_nodes_polygon = num_nodes;
436  while( actual_num_nodes_polygon >= 3 &&
437  conn[actual_num_nodes_polygon - 1] == conn[actual_num_nodes_polygon - 2] )
438  actual_num_nodes_polygon--;
439  num_sides = actual_num_nodes_polygon;
440  sub_type = MBEDGE;
441  num_sub_nodes = 2;
442  }
443  else // get connectivity of each n-1 dimension entity
444  num_sides = CN::NumSubEntities( type, mTargetDim );
445  for( int i = 0; i < num_sides; i++ )
446  {
447  if( MBPOLYGON == type )
448  {
449  sub_conn[0] = conn[i];
450  sub_conn[1] = conn[i + 1];
451  if( i + 1 == actual_num_nodes_polygon ) sub_conn[1] = conn[0];
452  }
453  else
454  {
455  CN::SubEntityNodeIndices( type, num_nodes, mTargetDim, i, sub_type, num_sub_nodes, sub_indices );
456  assert( (size_t)num_sub_nodes <= sizeof( sub_indices ) / sizeof( sub_indices[0] ) );
457  for( int j = 0; j < num_sub_nodes; j++ )
458  sub_conn[j] = conn[sub_indices[j]];
459  }
460 
461  // see if we can match this connectivity with
462  // an existing entity
463  find_match( sub_type, sub_conn, num_sub_nodes, match, direct );
464 
465  // if there is no match, create a new entity
466  if( match == 0 )
467  {
468  EntityHandle tmphndl = 0;
469  int indices[MAX_SUB_ENTITY_VERTICES];
470  EntityType new_type;
471  int num_new_nodes;
472  if( MBPOLYGON == type )
473  {
474  new_type = MBEDGE;
475  num_new_nodes = 2;
476  }
477  else
478  {
479  CN::SubEntityNodeIndices( type, num_nodes, mTargetDim, i, new_type, num_new_nodes, indices );
480  for( int j = 0; j < num_new_nodes; j++ )
481  sub_conn[j] = conn[indices[j]];
482  }
483  result = thisMB->create_element( new_type, sub_conn, num_new_nodes, tmphndl );
484  assert( MB_SUCCESS == result );
485  add_adjacency( tmphndl, sub_conn, CN::VerticesPerEntity( new_type ) );
486  forward_target_entities.insert( tmphndl );
487  }
488  // if there is a match, delete the matching entity
489  // if we can.
490  else
491  {
492  if( ( seek_iter = forward_target_entities.find( match ) ) != forward_target_entities.end() )
493  {
494  forward_target_entities.erase( seek_iter );
495  remove_adjacency( match );
496  if( /*!use_adjs &&*/ entity_deletable( match ) )
497  {
498  result = thisMB->delete_entities( &match, 1 );
499  assert( MB_SUCCESS == result );
500  }
501  }
502  else if( ( seek_iter = reverse_target_entities.find( match ) ) != reverse_target_entities.end() )
503  {
504  reverse_target_entities.erase( seek_iter );
505  remove_adjacency( match );
506  if( /*!use_adjs &&*/ entity_deletable( match ) )
507  {
508  result = thisMB->delete_entities( &match, 1 );
509  assert( MB_SUCCESS == result );
510  }
511  }
512  else
513  {
514  if( direct == FORWARD )
515  {
516  forward_target_entities.insert( match );
517  }
518  else
519  {
520  reverse_target_entities.insert( match );
521  }
522  }
523  }
524  }
525  }
526 
527  deinitialize();
528 
529  return MB_SUCCESS;
530 }
531 
532 void Skinner::find_match( EntityType type,
533  const EntityHandle* conn,
534  const int num_nodes,
535  EntityHandle& match,
536  Skinner::direction& direct )
537 {
538  match = 0;
539 
540  if( type == MBVERTEX )
541  {
542  match = *conn;
543  direct = FORWARD;
544  return;
545  }
546 
547  const EntityHandle* iter = std::min_element( conn, conn + num_nodes );
548 
549  std::vector< EntityHandle >* adj = NULL;
550 
551  ErrorCode result = thisMB->tag_get_data( mAdjTag, iter, 1, &adj );
552  if( result == MB_FAILURE || adj == NULL )
553  {
554  return;
555  }
556 
557  std::vector< EntityHandle >::iterator jter, end_jter;
558  end_jter = adj->end();
559 
560  const EntityHandle* tmp;
561  int num_verts;
562 
563  for( jter = adj->begin(); jter != end_jter; ++jter )
564  {
565  EntityType tmp_type;
566  tmp_type = thisMB->type_from_handle( *jter );
567 
568  if( type != tmp_type ) continue;
569 
570  result = thisMB->get_connectivity( *jter, tmp, num_verts, false );
571  assert( MB_SUCCESS == result && num_verts >= CN::VerticesPerEntity( type ) );
572  // FIXME: connectivity_match appears to work only for linear elements,
573  // so ignore higher-order nodes.
574  if( connectivity_match( conn, tmp, CN::VerticesPerEntity( type ), direct ) )
575  {
576  match = *jter;
577  break;
578  }
579  }
580 }
581 
583  const EntityHandle* conn2,
584  const int num_verts,
585  Skinner::direction& direct )
586 {
587  const EntityHandle* iter = std::find( conn2, conn2 + num_verts, conn1[0] );
588  if( iter == conn2 + num_verts ) return false;
589 
590  bool they_match = true;
591 
592  int i;
593  unsigned int j = iter - conn2;
594 
595  // first compare forward
596  for( i = 1; i < num_verts; ++i )
597  {
598  if( conn1[i] != conn2[( j + i ) % num_verts] )
599  {
600  they_match = false;
601  break;
602  }
603  }
604 
605  if( they_match == true )
606  {
607  // need to check for reversed edges here
608  direct = ( num_verts == 2 && j ) ? REVERSE : FORWARD;
609  return true;
610  }
611 
612  they_match = true;
613 
614  // then compare reverse
615  j += num_verts;
616  for( i = 1; i < num_verts; )
617  {
618  if( conn1[i] != conn2[( j - i ) % num_verts] )
619  {
620  they_match = false;
621  break;
622  }
623  ++i;
624  }
625  if( they_match )
626  {
627  direct = REVERSE;
628  }
629  return they_match;
630 }
631 
633 {
634  std::vector< EntityHandle > nodes, *adj = NULL;
635  ErrorCode result = thisMB->get_connectivity( &entity, 1, nodes );MB_CHK_ERR( result );
636  std::vector< EntityHandle >::iterator iter = std::min_element( nodes.begin(), nodes.end() );
637 
638  if( iter == nodes.end() ) return MB_FAILURE;
639 
640  // remove this entity from the node
641  if( thisMB->tag_get_data( mAdjTag, &( *iter ), 1, &adj ) == MB_SUCCESS && adj != NULL )
642  {
643  iter = std::find( adj->begin(), adj->end(), entity );
644  if( iter != adj->end() ) adj->erase( iter );
645  }
646 
647  return result;
648 }
649 
651 {
652  unsigned char deletable = 0;
653  ErrorCode result = thisMB->tag_get_data( mDeletableMBTag, &entity, 1, &deletable );
654  assert( MB_SUCCESS == result );
655  if( MB_SUCCESS == result && deletable == 1 ) return false;
656  return true;
657 }
658 
660  const Range& bar_elements,
661  EntityHandle boundary_edges,
662  EntityHandle inferred_edges,
663  EntityHandle non_manifold_edges,
664  EntityHandle other_edges,
665  int& number_boundary_nodes )
666 {
667  Range bedges, iedges, nmedges, oedges;
668  ErrorCode result =
669  classify_2d_boundary( boundary, bar_elements, bedges, iedges, nmedges, oedges, number_boundary_nodes );MB_CHK_ERR( result );
670 
671  // now set the input meshsets to the output ranges
672  result = thisMB->clear_meshset( &boundary_edges, 1 );MB_CHK_ERR( result );
673  result = thisMB->add_entities( boundary_edges, bedges );MB_CHK_ERR( result );
674 
675  result = thisMB->clear_meshset( &inferred_edges, 1 );MB_CHK_ERR( result );
676  result = thisMB->add_entities( inferred_edges, iedges );MB_CHK_ERR( result );
677 
678  result = thisMB->clear_meshset( &non_manifold_edges, 1 );MB_CHK_ERR( result );
679  result = thisMB->add_entities( non_manifold_edges, nmedges );MB_CHK_ERR( result );
680 
681  result = thisMB->clear_meshset( &other_edges, 1 );MB_CHK_ERR( result );
682  result = thisMB->add_entities( other_edges, oedges );MB_CHK_ERR( result );
683 
684  return MB_SUCCESS;
685 }
686 
688  const Range& bar_elements,
689  Range& boundary_edges,
690  Range& inferred_edges,
691  Range& non_manifold_edges,
692  Range& other_edges,
693  int& number_boundary_nodes )
694 {
695 
696  // clear out the edge lists
697 
698  boundary_edges.clear();
699  inferred_edges.clear();
700  non_manifold_edges.clear();
701  other_edges.clear();
702 
703  number_boundary_nodes = 0;
704 
705  // make sure we have something to work with
706  if( boundary.empty() )
707  {
708  return MB_FAILURE;
709  }
710 
711  // get our working dimensions
712  EntityType type = thisMB->type_from_handle( *( boundary.begin() ) );
713  const int source_dim = CN::Dimension( type );
714 
715  // make sure we can handle the working dimensions
716  if( source_dim != 2 )
717  {
718  return MB_FAILURE;
719  }
720  mTargetDim = source_dim - 1;
721 
722  // initialize
723  initialize();
724 
725  // additional initialization for this routine
726  // define a tag for MBEDGE which counts the occurrences of the edge below
727  // default should be 0 for existing edges, if any
728 
729  Tag count_tag;
730  int default_count = 0;
731  ErrorCode result =
732  thisMB->tag_get_handle( 0, 1, MB_TYPE_INTEGER, count_tag, MB_TAG_DENSE | MB_TAG_CREAT, &default_count );MB_CHK_ERR( result );
733 
734  Range::const_iterator iter, end_iter;
735  end_iter = boundary.end();
736 
737  std::vector< EntityHandle > conn;
738  EntityHandle sub_conn[2];
739  EntityHandle match;
740 
741  Range edge_list;
742  Range boundary_nodes;
743  Skinner::direction direct;
744 
745  EntityType sub_type;
746  int num_edge, num_sub_ent_vert;
747  const short* edge_verts;
748 
749  // now, process each entity in the boundary
750 
751  for( iter = boundary.begin(); iter != end_iter; ++iter )
752  {
753  // get the connectivity of this entity
754  conn.clear();
755  result = thisMB->get_connectivity( &( *iter ), 1, conn, false );
756  assert( MB_SUCCESS == result );
757 
758  // add node handles to boundary_node range
759  std::copy( conn.begin(), conn.begin() + CN::VerticesPerEntity( type ), range_inserter( boundary_nodes ) );
760 
761  type = thisMB->type_from_handle( *iter );
762 
763  // get connectivity of each n-1 dimension entity (edge in this case)
764  const struct CN::ConnMap* conn_map = &( CN::mConnectivityMap[type][0] );
765  num_edge = CN::NumSubEntities( type, 1 );
766  for( int i = 0; i < num_edge; i++ )
767  {
768  edge_verts = CN::SubEntityVertexIndices( type, 1, i, sub_type, num_sub_ent_vert );
769  assert( sub_type == MBEDGE && num_sub_ent_vert == 2 );
770  sub_conn[0] = conn[edge_verts[0]];
771  sub_conn[1] = conn[edge_verts[1]];
772  int num_sub_nodes = conn_map->num_corners_per_sub_element[i];
773 
774  // see if we can match this connectivity with
775  // an existing entity
776  find_match( MBEDGE, sub_conn, num_sub_nodes, match, direct );
777 
778  // if there is no match, create a new entity
779  if( match == 0 )
780  {
781  EntityHandle tmphndl = 0;
782  int indices[MAX_SUB_ENTITY_VERTICES];
783  EntityType new_type;
784  int num_new_nodes;
785  CN::SubEntityNodeIndices( type, conn.size(), 1, i, new_type, num_new_nodes, indices );
786  for( int j = 0; j < num_new_nodes; j++ )
787  sub_conn[j] = conn[indices[j]];
788 
789  result = thisMB->create_element( new_type, sub_conn, num_new_nodes, tmphndl );
790  assert( MB_SUCCESS == result );
791  add_adjacency( tmphndl, sub_conn, num_sub_nodes );
792  // target_entities.insert(tmphndl);
793  edge_list.insert( tmphndl );
794  int count;
795  result = thisMB->tag_get_data( count_tag, &tmphndl, 1, &count );
796  assert( MB_SUCCESS == result );
797  count++;
798  result = thisMB->tag_set_data( count_tag, &tmphndl, 1, &count );
799  assert( MB_SUCCESS == result );
800  }
801  else
802  {
803  // We found a match, we must increment the count on the match
804  int count;
805  result = thisMB->tag_get_data( count_tag, &match, 1, &count );
806  assert( MB_SUCCESS == result );
807  count++;
808  result = thisMB->tag_set_data( count_tag, &match, 1, &count );
809  assert( MB_SUCCESS == result );
810 
811  // if the entity is not deletable, it was pre-existing in
812  // the database. We therefore may need to add it to the
813  // edge_list. Since it will not hurt the range, we add
814  // whether it was added before or not
815  if( !entity_deletable( match ) )
816  {
817  edge_list.insert( match );
818  }
819  }
820  }
821  }
822 
823  // Any bar elements in the model should be classified separately
824  // If the element is in the skin edge_list, then it should be put in
825  // the non-manifold edge list. Edges not in the edge_list are stand-alone
826  // bars, and we make them simply boundary elements
827 
828  if( !bar_elements.empty() )
829  {
830  Range::iterator bar_iter;
831  for( iter = bar_elements.begin(); iter != bar_elements.end(); ++iter )
832  {
833  EntityHandle handle = *iter;
834  bar_iter = edge_list.find( handle );
835  if( bar_iter != edge_list.end() )
836  {
837  // it is in the list, erase it and put in non-manifold list
838  edge_list.erase( bar_iter );
839  non_manifold_edges.insert( handle );
840  }
841  else
842  {
843  // not in the edge list, make it a boundary edge
844  boundary_edges.insert( handle );
845  }
846  }
847  }
848 
849  // now all edges should be classified. Go through the edge_list,
850  // and put all in the appropriate lists
851 
852  Range::iterator edge_iter, edge_end_iter;
853  edge_end_iter = edge_list.end();
854  int count;
855  for( edge_iter = edge_list.begin(); edge_iter != edge_end_iter; ++edge_iter )
856  {
857  // check the count_tag
858  result = thisMB->tag_get_data( count_tag, &( *edge_iter ), 1, &count );
859  assert( MB_SUCCESS == result );
860  if( count == 1 )
861  {
862  boundary_edges.insert( *edge_iter );
863  }
864  else if( count == 2 )
865  {
866  other_edges.insert( *edge_iter );
867  }
868  else
869  {
870  non_manifold_edges.insert( *edge_iter );
871  }
872  }
873 
874  // find the inferred edges from the other_edge_list
875 
876  double min_angle_degrees = 20.0;
877  find_inferred_edges( const_cast< Range& >( boundary ), other_edges, inferred_edges, min_angle_degrees );
878 
879  // we now want to remove the inferred_edges from the other_edges
880 
881  Range temp_range;
882 
883  std::set_difference( other_edges.begin(), other_edges.end(), inferred_edges.begin(), inferred_edges.end(),
884  range_inserter( temp_range ), std::less< EntityHandle >() );
885 
886  other_edges = temp_range;
887 
888  // get rid of count tag and deinitialize
889 
890  result = thisMB->tag_delete( count_tag );
891  assert( MB_SUCCESS == result );
892  deinitialize();
893 
894  // set the node count
895  number_boundary_nodes = boundary_nodes.size();
896 
897  return MB_SUCCESS;
898 }
899 
900 void Skinner::find_inferred_edges( Range& skin_boundary,
901  Range& candidate_edges,
902  Range& inferred_edges,
903  double reference_angle_degrees )
904 {
905 
906  // mark all the entities in the skin boundary
907  Tag mark_tag;
908  ErrorCode result = thisMB->tag_get_handle( 0, 1, MB_TYPE_BIT, mark_tag, MB_TAG_CREAT );
909  assert( MB_SUCCESS == result );
910  unsigned char bit = true;
911  result = thisMB->tag_clear_data( mark_tag, skin_boundary, &bit );
912  assert( MB_SUCCESS == result );
913 
914  // find the cosine of the reference angle
915 
916  double reference_cosine = cos( reference_angle_degrees * SKINNER_PI / 180.0 );
917 
918  // check all candidate edges for an angle greater than the minimum
919 
920  Range::iterator iter, end_iter = candidate_edges.end();
921  std::vector< EntityHandle > adjacencies;
922  std::vector< EntityHandle >::iterator adj_iter;
923  EntityHandle face[2];
924 
925  for( iter = candidate_edges.begin(); iter != end_iter; ++iter )
926  {
927 
928  // get the 2D elements connected to this edge
929  adjacencies.clear();
930  result = thisMB->get_adjacencies( &( *iter ), 1, 2, false, adjacencies );
931  if( MB_SUCCESS != result ) continue;
932 
933  // there should be exactly two, that is why the edge is classified as nonBoundary
934  // and manifold
935 
936  int faces_found = 0;
937  for( adj_iter = adjacencies.begin(); adj_iter != adjacencies.end() && faces_found < 2; ++adj_iter )
938  {
939  // we need to find two of these which are in the skin
940  unsigned char is_marked = 0;
941  result = thisMB->tag_get_data( mark_tag, &( *adj_iter ), 1, &is_marked );
942  assert( MB_SUCCESS == result );
943  if( is_marked )
944  {
945  face[faces_found] = *adj_iter;
946  faces_found++;
947  }
948  }
949 
950  // assert(faces_found == 2 || faces_found == 0);
951  if( 2 != faces_found ) continue;
952 
953  // see if the two entities have a sufficient angle
954 
955  if( has_larger_angle( face[0], face[1], reference_cosine ) )
956  {
957  inferred_edges.insert( *iter );
958  }
959  }
960 
961  result = thisMB->tag_delete( mark_tag );
962  assert( MB_SUCCESS == result );
963 }
964 
965 bool Skinner::has_larger_angle( EntityHandle& entity1, EntityHandle& entity2, double reference_angle_cosine )
966 {
967  // compare normals to get angle. We assume that the surface quads
968  // which we test here will be approximately planar
969 
970  double norm[2][3];
971  Util::normal( thisMB, entity1, norm[0][0], norm[0][1], norm[0][2] );
972  Util::normal( thisMB, entity2, norm[1][0], norm[1][1], norm[1][2] );
973 
974  double cosine = norm[0][0] * norm[1][0] + norm[0][1] * norm[1][1] + norm[0][2] * norm[1][2];
975 
976  if( cosine < reference_angle_cosine )
977  {
978  return true;
979  }
980 
981  return false;
982 }
983 
984 // get skin entities of prescribed dimension
986  const Range& entities,
987  int dim,
988  Range& skin_entities,
989  bool create_vert_elem_adjs,
990  bool create_skin_elements )
991 {
992  Range tmp_skin;
993  ErrorCode result =
994  find_skin( this_set, entities, ( dim == 0 ), tmp_skin, 0, create_vert_elem_adjs, create_skin_elements );
995  if( MB_SUCCESS != result || tmp_skin.empty() ) return result;
996 
997  if( tmp_skin.all_of_dimension( dim ) )
998  {
999  if( skin_entities.empty() )
1000  skin_entities.swap( tmp_skin );
1001  else
1002  skin_entities.merge( tmp_skin );
1003  }
1004  else
1005  {
1006  result = thisMB->get_adjacencies( tmp_skin, dim, create_skin_elements, skin_entities, Interface::UNION );MB_CHK_ERR( result );
1007  if( this_set ) result = thisMB->add_entities( this_set, skin_entities );
1008  }
1009 
1010  return result;
1011 }
1012 
1014  const Range& entities,
1015  Range* skin_verts,
1016  Range* skin_elems,
1017  Range* skin_rev_elems,
1018  bool create_skin_elems,
1019  bool corners_only )
1020 {
1021  ErrorCode rval;
1022  if( entities.empty() ) return MB_SUCCESS;
1023 
1024  const int dim = CN::Dimension( TYPE_FROM_HANDLE( entities.front() ) );
1025  if( dim < 1 || dim > 3 || !entities.all_of_dimension( dim ) ) return MB_TYPE_OUT_OF_RANGE;
1026 
1027  // are we skinning all entities
1028  size_t count = entities.size();
1029  int num_total;
1030  rval = thisMB->get_number_entities_by_dimension( this_set, dim, num_total );
1031  if( MB_SUCCESS != rval ) return rval;
1032  bool all = ( count == (size_t)num_total );
1033 
1034  // Create a bit tag for fast intersection with input entities range.
1035  // If we're skinning all the entities in the mesh, we really don't
1036  // need the tag. To save memory, just create it with a default value
1037  // of one and don't set it. That way MOAB will return 1 for all
1038  // entities.
1039  Tag tag;
1040  char bit = all ? 1 : 0;
1041  rval = thisMB->tag_get_handle( NULL, 1, MB_TYPE_BIT, tag, MB_TAG_CREAT, &bit );
1042  if( MB_SUCCESS != rval ) return rval;
1043 
1044  // tag all entities in input range
1045  if( !all )
1046  {
1047  std::vector< unsigned char > vect( count, 1 );
1048  rval = thisMB->tag_set_data( tag, entities, &vect[0] );
1049  if( MB_SUCCESS != rval )
1050  {
1051  thisMB->tag_delete( tag );
1052  return rval;
1053  }
1054  }
1055 
1056  switch( dim )
1057  {
1058  case 1:
1059  if( skin_verts )
1060  rval = find_skin_vertices_1D( tag, entities, *skin_verts );
1061  else if( skin_elems )
1062  rval = find_skin_vertices_1D( tag, entities, *skin_elems );
1063  else
1064  rval = MB_SUCCESS;
1065  break;
1066  case 2:
1067  rval = find_skin_vertices_2D( this_set, tag, entities, skin_verts, skin_elems, skin_rev_elems,
1068  create_skin_elems, corners_only );
1069  break;
1070  case 3:
1071  rval = find_skin_vertices_3D( this_set, tag, entities, skin_verts, skin_elems, skin_rev_elems,
1072  create_skin_elems, corners_only );
1073  break;
1074  default:
1075  rval = MB_TYPE_OUT_OF_RANGE;
1076  break;
1077  }
1078 
1079  thisMB->tag_delete( tag );
1080  return rval;
1081 }
1082 
1083 ErrorCode Skinner::find_skin_vertices_1D( Tag tag, const Range& edges, Range& skin_verts )
1084 {
1085  // This rather simple algorithm is provided for completeness
1086  // (not sure how often one really wants the 'skin' of a chain
1087  // or tangle of edges.)
1088  //
1089  // A vertex is on the skin of the edges if it is contained in exactly
1090  // one of the edges *in the input range*.
1091  //
1092  // This function expects the caller to have tagged all edges in the
1093  // input range with a value of one for the passed bit tag, and all
1094  // other edges with a value of zero. This allows us to do a faster
1095  // intersection with the input range and the edges adjacent to a vertex.
1096 
1097  ErrorCode rval;
1098  Range::iterator hint = skin_verts.begin();
1099 
1100  // All input entities must be edges.
1101  if( !edges.all_of_dimension( 1 ) ) return MB_TYPE_OUT_OF_RANGE;
1102 
1103  // get all the vertices
1104  Range verts;
1105  rval = thisMB->get_connectivity( edges, verts, true );
1106  if( MB_SUCCESS != rval ) return rval;
1107 
1108  // Test how many edges each input vertex is adjacent to.
1109  std::vector< char > tag_vals;
1110  std::vector< EntityHandle > adj;
1111  int n;
1112  for( Range::const_iterator it = verts.begin(); it != verts.end(); ++it )
1113  {
1114  // get edges adjacent to vertex
1115  adj.clear();
1116  rval = thisMB->get_adjacencies( &*it, 1, 1, false, adj );
1117  if( MB_SUCCESS != rval ) return rval;
1118  if( adj.empty() ) continue;
1119 
1120  // Intersect adjacent edges with the input list of edges
1121  tag_vals.resize( adj.size() );
1122  rval = thisMB->tag_get_data( tag, &adj[0], adj.size(), &tag_vals[0] );
1123  if( MB_SUCCESS != rval ) return rval;
1124 #ifdef MOAB_OLD_STD_COUNT
1125  n = 0;
1126  std::count( tag_vals.begin(), tag_vals.end(), '\001' );
1127 #else
1128  n = std::count( tag_vals.begin(), tag_vals.end(), '\001' );
1129 #endif
1130  // If adjacent to only one input edge, then vertex is on skin
1131  if( n == 1 )
1132  {
1133  hint = skin_verts.insert( hint, *it );
1134  }
1135  }
1136 
1137  return MB_SUCCESS;
1138 }
1139 
1140 // A Container for storing a list of element sides adjacent
1141 // to a vertex. The template parameter is the number of
1142 // corners for the side.
1143 template < unsigned CORNERS >
1145 {
1146  public:
1147  /**
1148  * This struct is used to for a reduced representation of element
1149  * "sides" adjacent to a give vertex. As such, it
1150  * a) does not store the initial vertex that all sides are adjacent to
1151  * b) orders the remaining vertices in a specific way for fast comparison.
1152  *
1153  * For edge elements, only the opposite vertex is stored.
1154  * For triangle elements, only the other two vertices are stored,
1155  * and they are stored with the smaller of those two handles first.
1156  * For quad elements, only the other three vertices are stored.
1157  * They are stored such that the vertex opposite the implicit (not
1158  * stored) vertex is always in slot 1. The other two vertices
1159  * (in slots 0 and 2) are ordered such that the handle of the one in
1160  * slot 0 is smaller than the handle in slot 2.
1161  *
1162  * For each side, the adj_elem field is used to store the element that
1163  * it is a side of as long as the element is considered to be on the skin.
1164  * The adj_elem field is cleared (set to zero) to indicate that this
1165  * side is no longer considered to be on the skin (and is the side of
1166  * more than one element.)
1167  */
1168  struct Side
1169  {
1170  EntityHandle handles[CORNERS - 1]; //!< side vertices, except for implicit one
1171  EntityHandle adj_elem; //!< element that this is a side of, or zero
1172  bool skin() const
1173  {
1174  return 0 != adj_elem;
1175  }
1176 
1177  /** construct from connectivity of side
1178  *\param array The connectivity of the element side.
1179  *\param idx The index of the implicit vertex (contained
1180  * in all sides in the list.)
1181  *\param adj The element that this is a side of.
1182  */
1183  Side( const EntityHandle* array, int idx, EntityHandle adj, unsigned short ) : adj_elem( adj )
1184  {
1185  switch( CORNERS )
1186  {
1187  case 3:
1188  handles[1] = array[( idx + 2 ) % CORNERS];
1189  // fall through
1190  case 2:
1191  if( 3 == CORNERS ) handles[1] = array[( idx + 2 ) % CORNERS];
1192  if( 2 <= CORNERS ) handles[0] = array[( idx + 1 ) % CORNERS];
1193  break;
1194  default:
1195  assert( false );
1196  break;
1197  }
1198  if( CORNERS == 3 && handles[1] > handles[0] ) std::swap( handles[0], handles[1] );
1199  }
1200 
1201  /** construct from connectivity of parent element
1202  *\param array The connectivity of the parent element
1203  *\param idx The index of the implicit vertex (contained
1204  * in all sides in the list.) This is an index
1205  * into 'indices', not 'array'.
1206  *\param adj The element that this is a side of.
1207  *\param indices The indices into 'array' at which the vertices
1208  * representing the side occur.
1209  */
1210  Side( const EntityHandle* array, int idx, EntityHandle adj, unsigned short, const short* indices )
1211  : adj_elem( adj )
1212  {
1213  switch( CORNERS )
1214  {
1215  case 3:
1216  handles[1] = array[indices[( idx + 2 ) % CORNERS]];
1217  // fall through
1218  case 2:
1219  if( 3 == CORNERS ) handles[1] = array[indices[( idx + 2 ) % CORNERS]];
1220  if( 2 <= CORNERS ) handles[0] = array[indices[( idx + 1 ) % CORNERS]];
1221  break;
1222  default:
1223  assert( false );
1224  break;
1225  }
1226  if( CORNERS == 3 && handles[1] > handles[0] ) std::swap( handles[0], handles[1] );
1227  }
1228 
1229  // Compare two side instances. Relies in the ordering of
1230  // vertex handles as described above.
1231  bool operator==( const Side& other ) const
1232  {
1233  switch( CORNERS )
1234  {
1235  case 3:
1236  return handles[0] == other.handles[0] && handles[1] == other.handles[1];
1237  case 2:
1238  return handles[0] == other.handles[0];
1239  default:
1240  assert( false );
1241  return false;
1242  }
1243  }
1244  };
1245 
1246  private:
1247  std::vector< Side > data; //!< List of sides
1248  size_t skin_count; //!< Cached count of sides that are skin
1249 
1250  public:
1251  typedef typename std::vector< Side >::iterator iterator;
1252  typedef typename std::vector< Side >::const_iterator const_iterator;
1254  {
1255  return data.begin();
1256  }
1258  {
1259  return data.end();
1260  }
1261 
1262  void clear()
1263  {
1264  data.clear();
1265  skin_count = 0;
1266  }
1267  bool empty() const
1268  {
1269  return data.empty();
1270  }
1271 
1272  AdjSides() : skin_count( 0 ) {}
1273 
1274  size_t num_skin() const
1275  {
1276  return skin_count;
1277  }
1278 
1279  /** \brief insert side, specifying side connectivity
1280  *
1281  * Either insert a new side, or if the side is already in the
1282  * list, mark it as not on the skin.
1283  *
1284  *\param handles The connectivity of the element side.
1285  *\param skip_idx The index of the implicit vertex (contained
1286  * in all sides in the list.)
1287  *\param adj_elem The element that this is a side of.
1288  *\param elem_side Which side of adj_elem are we storing
1289  * (CN side number.)
1290  */
1291  void insert( const EntityHandle* handles, int skip_idx, EntityHandle adj_elem, unsigned short elem_side )
1292  {
1293  Side side( handles, skip_idx, adj_elem, elem_side );
1294  iterator p = std::find( data.begin(), data.end(), side );
1295  if( p == data.end() )
1296  {
1297  data.push_back( side );
1298  ++skin_count; // not in list yet, so skin side (so far)
1299  }
1300  else if( p->adj_elem )
1301  {
1302  p->adj_elem = 0; // mark as not on skin
1303  --skin_count; // decrement cached count of skin elements
1304  }
1305  }
1306 
1307  /** \brief insert side, specifying list of indices into parent element
1308  * connectivity.
1309  *
1310  * Either insert a new side, or if the side is already in the
1311  * list, mark it as not on the skin.
1312  *
1313  *\param handles The connectivity of the parent element
1314  *\param skip_idx The index of the implicit vertex (contained
1315  * in all sides in the list.) This is an index
1316  * into 'indices', not 'handles'.
1317  *\param adj_elem The element that this is a side of (parent handle).
1318  *\param indices The indices into 'handles' at which the vertices
1319  * representing the side occur.
1320  *\param elem_side Which side of adj_elem are we storing
1321  * (CN side number.)
1322  */
1323  void insert( const EntityHandle* handles,
1324  int skip_idx,
1325  EntityHandle adj_elem,
1326  unsigned short elem_side,
1327  const short* indices )
1328  {
1329  Side side( handles, skip_idx, adj_elem, elem_side, indices );
1330  iterator p = std::find( data.begin(), data.end(), side );
1331  if( p == data.end() )
1332  {
1333  data.push_back( side );
1334  ++skin_count; // not in list yet, so skin side (so far)
1335  }
1336  else if( p->adj_elem )
1337  {
1338  p->adj_elem = 0; // mark as not on skin
1339  --skin_count; // decrement cached count of skin elements
1340  }
1341  }
1342 
1343  /**\brief Search list for a given side, and if found, mark as not skin.
1344  *
1345  *\param other Connectivity of side
1346  *\param skip_index Index in 'other' at which implicit vertex occurs.
1347  *\param elem_out If return value is true, the element that the side is a
1348  * side of. If return value is false, not set.
1349  *\return true if found and marked not-skin, false if not found.
1350  *
1351  * Given the connectivity of some existing element, check if it occurs
1352  * in the list. If it does, clear the "is skin" state of the side so
1353  * that we know that we don't need to later create the side element.
1354  */
1355  bool find_and_unmark( const EntityHandle* other, int skip_index, EntityHandle& elem_out )
1356  {
1357  Side s( other, skip_index, 0, 0 );
1358  iterator p = std::find( data.begin(), data.end(), s );
1359  if( p == data.end() || !p->adj_elem )
1360  return false;
1361  else
1362  {
1363  elem_out = p->adj_elem;
1364  p->adj_elem = 0; // clear "is skin" state for side
1365  --skin_count; // decrement cached count of skin sides
1366  return true;
1367  }
1368  }
1369 };
1370 
1371 /** construct from connectivity of side
1372  *\param array The connectivity of the element side.
1373  *\param idx The index of the implicit vertex (contained
1374  * in all sides in the list.)
1375  *\param adj The element that this is a side of.
1376  */
1377 template <>
1378 AdjSides< 4 >::Side::Side( const EntityHandle* array, int idx, EntityHandle adj, unsigned short ) : adj_elem( adj )
1379 {
1380  const unsigned int CORNERS = 4;
1381  handles[2] = array[( idx + 3 ) % CORNERS];
1382  handles[1] = array[( idx + 2 ) % CORNERS];
1383  handles[0] = array[( idx + 1 ) % CORNERS];
1384  if( handles[2] > handles[0] ) std::swap( handles[0], handles[2] );
1385 }
1386 
1387 /** construct from connectivity of parent element
1388  *\param array The connectivity of the parent element
1389  *\param idx The index of the implicit vertex (contained
1390  * in all sides in the list.) This is an index
1391  * into 'indices', not 'array'.
1392  *\param adj The element that this is a side of.
1393  *\param indices The indices into 'array' at which the vertices
1394  * representing the side occur.
1395  */
1396 template <>
1397 AdjSides< 4 >::Side::Side( const EntityHandle* array, int idx, EntityHandle adj, unsigned short, const short* indices )
1398  : adj_elem( adj )
1399 {
1400  const unsigned int CORNERS = 4;
1401  handles[2] = array[indices[( idx + 3 ) % CORNERS]];
1402  handles[1] = array[indices[( idx + 2 ) % CORNERS]];
1403  handles[0] = array[indices[( idx + 1 ) % CORNERS]];
1404  if( handles[2] > handles[0] ) std::swap( handles[0], handles[2] );
1405 }
1406 
1407 // Compare two side instances. Relies in the ordering of
1408 // vertex handles as described above.
1409 template <>
1410 bool AdjSides< 4 >::Side::operator==( const Side& other ) const
1411 {
1412  return handles[0] == other.handles[0] && handles[1] == other.handles[1] && handles[2] == other.handles[2];
1413 }
1414 
1415 // Utility function used by find_skin_vertices_2D and
1416 // find_skin_vertices_3D to create elements representing
1417 // the skin side of a higher-dimension element if one
1418 // does not already exist.
1419 //
1420 // Some arguments may seem redundant, but they are used
1421 // to create the correct order of element when the input
1422 // element contains higher-order nodes.
1423 //
1424 // This function always creates elements that have a "forward"
1425 // orientation with respect to the parent element (have
1426 // nodes ordered the same as CN returns for the "side").
1427 //
1428 // elem - The higher-dimension element for which to create
1429 // a lower-dim element representing the side.
1430 // side_type - The EntityType of the desired side.
1431 // side_conn - The connectivity of the new side.
1433  EntityHandle elem,
1434  EntityType side_type,
1435  const EntityHandle* side_conn,
1436  EntityHandle& side_elem )
1437 {
1438  const int max_side = 9;
1439  const EntityHandle* conn;
1440  int len, side_len, side, sense, offset, indices[max_side];
1441  ErrorCode rval;
1442  EntityType type = TYPE_FROM_HANDLE( elem ), tmp_type;
1443  const int ncorner = CN::VerticesPerEntity( side_type );
1444  const int d = CN::Dimension( side_type );
1445  std::vector< EntityHandle > storage;
1446 
1447  // Get the connectivity of the parent element
1448  rval = thisMB->get_connectivity( elem, conn, len, false, &storage );
1449  if( MB_SUCCESS != rval ) return rval;
1450 
1451  // treat separately MBPOLYGON; we want to create the edge in the
1452  // forward sense always ; so figure out the sense first, then get out
1453  if( MBPOLYGON == type && 1 == d && MBEDGE == side_type )
1454  {
1455  // first find the first vertex in the conn list
1456  int i = 0;
1457  for( i = 0; i < len; i++ )
1458  {
1459  if( conn[i] == side_conn[0] ) break;
1460  }
1461  if( len == i ) return MB_FAILURE; // not found, big error
1462  // now, what if the polygon is padded?
1463  // the previous index is fine always. but the next one could be trouble :(
1464  int prevIndex = ( i + len - 1 ) % len;
1465  int nextIndex = ( i + 1 ) % len;
1466  // if the next index actually point to the same node, as current, it means it is padded
1467  if( conn[nextIndex] == conn[i] )
1468  {
1469  // it really means we are at the end of proper nodes, the last nodes are repeated, so it
1470  // should be the first node
1471  nextIndex = 0; // this is the first node!
1472  }
1473  EntityHandle conn2[2] = { side_conn[0], side_conn[1] };
1474  if( conn[prevIndex] == side_conn[1] )
1475  {
1476  // reverse, so the edge will be forward
1477  conn2[0] = side_conn[1];
1478  conn2[1] = side_conn[0];
1479  }
1480  else if( conn[nextIndex] != side_conn[1] )
1481  return MB_FAILURE; // it is not adjacent to the polygon
1482 
1483  rval = thisMB->create_element( MBEDGE, conn2, 2, side_elem );MB_CHK_ERR( rval );
1484  if( this_set )
1485  {
1486  rval = thisMB->add_entities( this_set, &side_elem, 1 );MB_CHK_ERR( rval );
1487  }
1488  return MB_SUCCESS;
1489  }
1490  // Find which side we are creating and get indices of all nodes
1491  // (including higher-order, if any.)
1492  CN::SideNumber( type, conn, side_conn, ncorner, d, side, sense, offset );
1493  CN::SubEntityNodeIndices( type, len, d, side, tmp_type, side_len, indices );
1494  assert( side_len <= max_side );
1495  assert( side_type == tmp_type );
1496 
1497  // NOTE: re-create conn array even when no higher-order nodes
1498  // because we want it to always be forward with respect
1499  // to the side ordering.
1500  EntityHandle side_conn_full[max_side];
1501  for( int i = 0; i < side_len; ++i )
1502  side_conn_full[i] = conn[indices[i]];
1503 
1504  rval = thisMB->create_element( side_type, side_conn_full, side_len, side_elem );MB_CHK_ERR( rval );
1505  if( this_set )
1506  {
1507  rval = thisMB->add_entities( this_set, &side_elem, 1 );MB_CHK_ERR( rval );
1508  }
1509  return MB_SUCCESS;
1510  ;
1511 }
1512 
1513 // Test if an edge is reversed with respect CN's ordering
1514 // for the "side" of a face.
1516 {
1517  const EntityHandle* conn;
1518  int len, idx;
1519  ErrorCode rval = thisMB->get_connectivity( face, conn, len, true );
1520  if( MB_SUCCESS != rval )
1521  {
1522  assert( false );
1523  return false;
1524  }
1525  idx = std::find( conn, conn + len, edge_ends[0] ) - conn;
1526  if( idx == len )
1527  {
1528  assert( false );
1529  return false;
1530  }
1531  return ( edge_ends[1] == conn[( idx + len - 1 ) % len] );
1532 }
1533 
1534 // Test if a 2D element representing the side or face of a
1535 // volume element is reversed with respect to the CN node
1536 // ordering for the corresponding region element side.
1537 bool Skinner::face_reversed( EntityHandle region, const EntityHandle* face_corners, EntityType face_type )
1538 {
1539  const EntityHandle* conn;
1540  int len, side, sense, offset;
1541  ErrorCode rval = thisMB->get_connectivity( region, conn, len, true );
1542  if( MB_SUCCESS != rval )
1543  {
1544  assert( false );
1545  return false;
1546  }
1547  short r = CN::SideNumber( TYPE_FROM_HANDLE( region ), conn, face_corners, CN::VerticesPerEntity( face_type ),
1548  CN::Dimension( face_type ), side, sense, offset );
1549  assert( 0 == r );
1550  return ( !r && sense == -1 );
1551 }
1552 
1554  Tag tag,
1555  const Range& faces,
1556  Range* skin_verts,
1557  Range* skin_edges,
1558  Range* reversed_edges,
1559  bool create_edges,
1560  bool corners_only )
1561 {
1562  // This function iterates over all the vertices contained in the
1563  // input face list. For each such vertex, it then iterates over
1564  // all of the sides of the face elements adjacent to the vertex.
1565  // If an adjacent side is the side of only one of the input
1566  // faces, then that side is on the skin.
1567  //
1568  // This algorithm will visit each skin vertex exactly once. It
1569  // will visit each skin side once for each vertex in the side.
1570  //
1571  // This function expects the caller to have created the passed bit
1572  // tag and set it to one only for the faces in the passed range. This
1573  // tag is used to do a fast intersection of the faces adjacent to a
1574  // vertex with the faces in the input range (discard any for which the
1575  // tag is not set to one.)
1576 
1577  ErrorCode rval;
1578  std::vector< EntityHandle >::iterator i, j;
1579  Range::iterator hint;
1580  if( skin_verts ) hint = skin_verts->begin();
1581  std::vector< EntityHandle > storage;
1582  const EntityHandle* conn;
1583  int len;
1584  bool find_edges = skin_edges || create_edges;
1585  bool printed_nonconformal_ho_warning = false;
1587 
1588  if( !faces.all_of_dimension( 2 ) ) return MB_TYPE_OUT_OF_RANGE;
1589 
1590  // get all the vertices
1591  Range verts;
1592  rval = thisMB->get_connectivity( faces, verts, true );
1593  if( MB_SUCCESS != rval ) return rval;
1594 
1595  std::vector< char > tag_vals;
1596  std::vector< EntityHandle > adj;
1597  AdjSides< 2 > adj_edges;
1598  for( Range::const_iterator it = verts.begin(); it != verts.end(); ++it )
1599  {
1600  bool higher_order = false;
1601 
1602  // get all adjacent faces
1603  adj.clear();
1604  rval = thisMB->get_adjacencies( &*it, 1, 2, false, adj );
1605  if( MB_SUCCESS != rval ) return rval;
1606  if( adj.empty() ) continue;
1607 
1608  // remove those not in the input list (intersect with input list)
1609  i = j = adj.begin();
1610  tag_vals.resize( adj.size() );
1611  rval = thisMB->tag_get_data( tag, &adj[0], adj.size(), &tag_vals[0] );
1612  if( MB_SUCCESS != rval ) return rval;
1613  // remove non-tagged entries
1614  i = j = adj.begin();
1615  for( ; i != adj.end(); ++i )
1616  if( tag_vals[i - adj.begin()] ) *( j++ ) = *i;
1617  adj.erase( j, adj.end() );
1618 
1619  // For each adjacent face, check the edges adjacent to the current vertex
1620  adj_edges.clear(); // other vertex for adjacent edges
1621  for( i = adj.begin(); i != adj.end(); ++i )
1622  {
1623  rval = thisMB->get_connectivity( *i, conn, len, false, &storage );
1624  if( MB_SUCCESS != rval ) return rval;
1625 
1626  // For a single face element adjacent to this vertex, there
1627  // will be exactly two sides (edges) adjacent to the vertex.
1628  // Find the other vertex for each of the two edges.
1629 
1630  EntityHandle prev, next; // vertices of two adjacent edge-sides
1631  const int idx = std::find( conn, conn + len, *it ) - conn;
1632  assert( idx != len );
1633 
1634  if( TYPE_FROM_HANDLE( *i ) == MBTRI && len > 3 )
1635  {
1636  len = 3;
1637  higher_order = true;
1638  if( idx > 2 )
1639  { // skip higher-order nodes for now
1640  if( !printed_nonconformal_ho_warning )
1641  {
1642  printed_nonconformal_ho_warning = true;
1643  std::cerr << "Non-conformal higher-order mesh detected in skinner: "
1644  << "vertex " << ID_FROM_HANDLE( *it ) << " is a corner in "
1645  << "some elements and a higher-order node in others" << std::endl;
1646  }
1647  continue;
1648  }
1649  }
1650  else if( TYPE_FROM_HANDLE( *i ) == MBQUAD && len > 4 )
1651  {
1652  len = 4;
1653  higher_order = true;
1654  if( idx > 3 )
1655  { // skip higher-order nodes for now
1656  if( !printed_nonconformal_ho_warning )
1657  {
1658  printed_nonconformal_ho_warning = true;
1659  std::cerr << "Non-conformal higher-order mesh detected in skinner: "
1660  << "vertex " << ID_FROM_HANDLE( *it ) << " is a corner in "
1661  << "some elements and a higher-order node in others" << std::endl;
1662  }
1663  continue;
1664  }
1665  }
1666 
1667  // so it must be a MBPOLYGON
1668  const int prev_idx = ( idx + len - 1 ) % len; // this should be fine, always, even for padded case
1669  prev = conn[prev_idx];
1670  next = conn[( idx + 1 ) % len];
1671  if( next == conn[idx] ) // it must be the padded case, so roll to the beginning
1672  next = conn[0];
1673 
1674  // Insert sides (edges) in our list of candidate skin sides
1675  adj_edges.insert( &prev, 1, *i, prev_idx );
1676  adj_edges.insert( &next, 1, *i, idx );
1677  }
1678 
1679  // If vertex is not on skin, advance to next vertex.
1680  // adj_edges handled checking for duplicates on insertion.
1681  // If every candidate skin edge occurred more than once (was
1682  // not in fact on the skin), then we're done with this vertex.
1683  if( 0 == adj_edges.num_skin() ) continue;
1684 
1685  // If user requested Range of *vertices* on the skin...
1686  if( skin_verts )
1687  {
1688  // Put skin vertex in output list
1689  hint = skin_verts->insert( hint, *it );
1690 
1691  // Add mid edge nodes to vertex list
1692  if( !corners_only && higher_order )
1693  {
1694  for( AdjSides< 2 >::const_iterator p = adj_edges.begin(); p != adj_edges.end(); ++p )
1695  {
1696  if( p->skin() )
1697  {
1698  face = p->adj_elem;
1699  EntityType type = TYPE_FROM_HANDLE( face );
1700 
1701  rval = thisMB->get_connectivity( face, conn, len, false );
1702  if( MB_SUCCESS != rval ) return rval;
1703  if( !CN::HasMidEdgeNodes( type, len ) ) continue;
1704 
1705  EntityHandle ec[2] = { *it, p->handles[0] };
1706  int side, sense, offset;
1707  CN::SideNumber( type, conn, ec, 2, 1, side, sense, offset );
1708  offset = CN::HONodeIndex( type, len, 1, side );
1709  assert( offset >= 0 && offset < len );
1710  skin_verts->insert( conn[offset] );
1711  }
1712  }
1713  }
1714  }
1715 
1716  // If user requested Range of *edges* on the skin...
1717  if( find_edges )
1718  {
1719  // Search list of existing adjacent edges for any that are on the skin
1720  adj.clear();
1721  rval = thisMB->get_adjacencies( &*it, 1, 1, false, adj );
1722  if( MB_SUCCESS != rval ) return rval;
1723  for( i = adj.begin(); i != adj.end(); ++i )
1724  {
1725  rval = thisMB->get_connectivity( *i, conn, len, true );
1726  if( MB_SUCCESS != rval ) return rval;
1727 
1728  // bool equality expression within find_and_unmark call
1729  // will be evaluate to the index of *it in the conn array.
1730  //
1731  // Note that the order of the terms in the if statement is important.
1732  // We want to unmark any existing skin edges even if we aren't
1733  // returning them. Otherwise we'll end up creating duplicates
1734  // if create_edges is true and skin_edges is not.
1735  if( adj_edges.find_and_unmark( conn, ( conn[1] == *it ), face ) && skin_edges )
1736  {
1737  if( reversed_edges && edge_reversed( face, conn ) )
1738  reversed_edges->insert( *i );
1739  else
1740  skin_edges->insert( *i );
1741  }
1742  }
1743  }
1744 
1745  // If the user requested that we create new edges for sides
1746  // on the skin for which there is no existing edge, and there
1747  // are still skin sides for which no corresponding edge was found...
1748  if( create_edges && adj_edges.num_skin() )
1749  {
1750  // Create any skin edges that don't exist
1751  for( AdjSides< 2 >::const_iterator p = adj_edges.begin(); p != adj_edges.end(); ++p )
1752  {
1753  if( p->skin() )
1754  {
1755  EntityHandle edge, ec[] = { *it, p->handles[0] };
1756  rval = create_side( this_set, p->adj_elem, MBEDGE, ec, edge );
1757  if( MB_SUCCESS != rval ) return rval;
1758  if( skin_edges ) skin_edges->insert( edge );
1759  }
1760  }
1761  }
1762 
1763  } // end for each vertex
1764 
1765  return MB_SUCCESS;
1766 }
1767 
1769  Tag tag,
1770  const Range& entities,
1771  Range* skin_verts,
1772  Range* skin_faces,
1773  Range* reversed_faces,
1774  bool create_faces,
1775  bool corners_only )
1776 {
1777  // This function iterates over all the vertices contained in the
1778  // input vol elem list. For each such vertex, it then iterates over
1779  // all of the sides of the vol elements adjacent to the vertex.
1780  // If an adjacent side is the side of only one of the input
1781  // elements, then that side is on the skin.
1782  //
1783  // This algorithm will visit each skin vertex exactly once. It
1784  // will visit each skin side once for each vertex in the side.
1785  //
1786  // This function expects the caller to have created the passed bit
1787  // tag and set it to one only for the elements in the passed range. This
1788  // tag is used to do a fast intersection of the elements adjacent to a
1789  // vertex with the elements in the input range (discard any for which the
1790  // tag is not set to one.)
1791  //
1792  // For each vertex, iterate over each adjacent element. Construct
1793  // lists of the sides of each adjacent element that contain the vertex.
1794  //
1795  // A list of three-vertex sides is kept for all triangular faces,
1796  // included three-vertex faces of type MBPOLYGON. Putting polygons
1797  // in the same list ensures that we find polyhedron and non-polyhedron
1798  // elements that are adjacent.
1799  //
1800  // A list of four-vertex sides is kept for all quadrilateral faces,
1801  // including four-vertex faces of type MBPOLYGON.
1802  //
1803  // Sides with more than four vertices must have an explicit MBPOLYGON
1804  // element representing them because MBPOLYHEDRON connectivity is a
1805  // list of faces rather than vertices. So the third list (vertices>=5),
1806  // need contain only the handle of the face rather than the vertex handles.
1807 
1808  ErrorCode rval;
1809  std::vector< EntityHandle >::iterator i, j;
1810  Range::iterator hint;
1811  if( skin_verts ) hint = skin_verts->begin();
1812  std::vector< EntityHandle > storage, storage2; // temp storage for conn lists
1813  const EntityHandle *conn, *conn2;
1814  int len, len2;
1815  bool find_faces = skin_faces || create_faces;
1816  int clen, side, sense, offset, indices[9];
1817  EntityType face_type;
1818  EntityHandle elem;
1819  bool printed_nonconformal_ho_warning = false;
1820 
1821  if( !entities.all_of_dimension( 3 ) ) return MB_TYPE_OUT_OF_RANGE;
1822 
1823  // get all the vertices
1824  Range verts;
1825  rval = thisMB->get_connectivity( entities, verts, true );
1826  if( MB_SUCCESS != rval ) return rval;
1827  // if there are polyhedra in the input list, need to make another
1828  // call to get vertices from faces
1829  if( !verts.all_of_dimension( 0 ) )
1830  {
1831  Range::iterator it = verts.upper_bound( MBVERTEX );
1832  Range pfaces;
1833  pfaces.merge( it, verts.end() );
1834  verts.erase( it, verts.end() );
1835  rval = thisMB->get_connectivity( pfaces, verts, true );
1836  if( MB_SUCCESS != rval ) return rval;
1837  assert( verts.all_of_dimension( 0 ) );
1838  }
1839 
1840  AdjSides< 4 > adj_quads; // 4-node sides adjacent to a vertex
1841  AdjSides< 3 > adj_tris; // 3-node sides adjacent to a vertex
1842  AdjSides< 2 > adj_poly; // n-node sides (n>5) adjacent to vertex
1843  // (must have an explicit polygon, so store
1844  // polygon handle rather than vertices.)
1845  std::vector< char > tag_vals;
1846  std::vector< EntityHandle > adj;
1847  for( Range::const_iterator it = verts.begin(); it != verts.end(); ++it )
1848  {
1849  bool higher_order = false;
1850 
1851  // get all adjacent elements
1852  adj.clear();
1853  rval = thisMB->get_adjacencies( &*it, 1, 3, false, adj );
1854  if( MB_SUCCESS != rval ) return rval;
1855  if( adj.empty() ) continue;
1856 
1857  // remove those not tagged (intersect with input range)
1858  i = j = adj.begin();
1859  tag_vals.resize( adj.size() );
1860  rval = thisMB->tag_get_data( tag, &adj[0], adj.size(), &tag_vals[0] );
1861  if( MB_SUCCESS != rval ) return rval;
1862  for( ; i != adj.end(); ++i )
1863  if( tag_vals[i - adj.begin()] ) *( j++ ) = *i;
1864  adj.erase( j, adj.end() );
1865 
1866  // Build lists of sides of 3D element adjacent to the current vertex
1867  adj_quads.clear(); // store three other vertices for each adjacent quad face
1868  adj_tris.clear(); // store two other vertices for each adjacent tri face
1869  adj_poly.clear(); // store handle of each adjacent polygonal face
1870  int idx;
1871  for( i = adj.begin(); i != adj.end(); ++i )
1872  {
1873  const EntityType type = TYPE_FROM_HANDLE( *i );
1874 
1875  // Special case for POLYHEDRA
1876  if( type == MBPOLYHEDRON )
1877  {
1878  rval = thisMB->get_connectivity( *i, conn, len );
1879  if( MB_SUCCESS != rval ) return rval;
1880  for( int k = 0; k < len; ++k )
1881  {
1882  rval = thisMB->get_connectivity( conn[k], conn2, len2, true, &storage2 );
1883  if( MB_SUCCESS != rval ) return rval;
1884  idx = std::find( conn2, conn2 + len2, *it ) - conn2;
1885  if( idx == len2 ) // vertex not in this face
1886  continue;
1887 
1888  // Treat 3- and 4-vertex faces specially, so that
1889  // if the mesh contains both elements and polyhedra,
1890  // we don't miss one type adjacent to the other.
1891  switch( len2 )
1892  {
1893  case 3:
1894  adj_tris.insert( conn2, idx, *i, k );
1895  break;
1896  case 4:
1897  adj_quads.insert( conn2, idx, *i, k );
1898  break;
1899  default:
1900  adj_poly.insert( conn + k, 1, *i, k );
1901  break;
1902  }
1903  }
1904  }
1905  else
1906  {
1907  rval = thisMB->get_connectivity( *i, conn, len, false, &storage );
1908  if( MB_SUCCESS != rval ) return rval;
1909 
1910  idx = std::find( conn, conn + len, *it ) - conn;
1911  assert( idx != len );
1912 
1913  if( len > CN::VerticesPerEntity( type ) )
1914  {
1915  higher_order = true;
1916  // skip higher-order nodes for now
1917  if( idx >= CN::VerticesPerEntity( type ) )
1918  {
1919  if( !printed_nonconformal_ho_warning )
1920  {
1921  printed_nonconformal_ho_warning = true;
1922  std::cerr << "Non-conformal higher-order mesh detected in skinner: "
1923  << "vertex " << ID_FROM_HANDLE( *it ) << " is a corner in "
1924  << "some elements and a higher-order node in others" << std::endl;
1925  }
1926  continue;
1927  }
1928  }
1929 
1930  // For each side of the element...
1931  const int num_faces = CN::NumSubEntities( type, 2 );
1932  for( int f = 0; f < num_faces; ++f )
1933  {
1934  int num_vtx;
1935  const short* face_indices = CN::SubEntityVertexIndices( type, 2, f, face_type, num_vtx );
1936  const short face_idx = std::find( face_indices, face_indices + num_vtx, (short)idx ) - face_indices;
1937  // skip sides that do not contain vertex from outer loop
1938  if( face_idx == num_vtx ) continue; // current vertex not in this face
1939 
1940  assert( num_vtx <= 4 ); // polyhedra handled above
1941  switch( face_type )
1942  {
1943  case MBTRI:
1944  adj_tris.insert( conn, face_idx, *i, f, face_indices );
1945  break;
1946  case MBQUAD:
1947  adj_quads.insert( conn, face_idx, *i, f, face_indices );
1948  break;
1949  default:
1950  return MB_TYPE_OUT_OF_RANGE;
1951  }
1952  }
1953  }
1954  } // end for (adj[3])
1955 
1956  // If vertex is not on skin, advance to next vertex
1957  if( 0 == ( adj_tris.num_skin() + adj_quads.num_skin() + adj_poly.num_skin() ) ) continue;
1958 
1959  // If user requested that skin *vertices* be passed back...
1960  if( skin_verts )
1961  {
1962  // Put skin vertex in output list
1963  hint = skin_verts->insert( hint, *it );
1964 
1965  // Add mid-edge and mid-face nodes to vertex list
1966  if( !corners_only && higher_order )
1967  {
1968  for( AdjSides< 3 >::const_iterator t = adj_tris.begin(); t != adj_tris.end(); ++t )
1969  {
1970  if( t->skin() )
1971  {
1972  elem = t->adj_elem;
1973  EntityType type = TYPE_FROM_HANDLE( elem );
1974 
1975  rval = thisMB->get_connectivity( elem, conn, len, false );
1976  if( MB_SUCCESS != rval ) return rval;
1977  if( !CN::HasMidNodes( type, len ) ) continue;
1978 
1979  EntityHandle ec[3] = { *it, t->handles[0], t->handles[1] };
1980  CN::SideNumber( type, conn, ec, 3, 2, side, sense, offset );
1981  CN::SubEntityNodeIndices( type, len, 2, side, face_type, clen, indices );
1982  assert( MBTRI == face_type );
1983  for( int k = 3; k < clen; ++k )
1984  skin_verts->insert( conn[indices[k]] );
1985  }
1986  }
1987  for( AdjSides< 4 >::const_iterator q = adj_quads.begin(); q != adj_quads.end(); ++q )
1988  {
1989  if( q->skin() )
1990  {
1991  elem = q->adj_elem;
1992  EntityType type = TYPE_FROM_HANDLE( elem );
1993 
1994  rval = thisMB->get_connectivity( elem, conn, len, false );
1995  if( MB_SUCCESS != rval ) return rval;
1996  if( !CN::HasMidNodes( type, len ) ) continue;
1997 
1998  EntityHandle ec[4] = { *it, q->handles[0], q->handles[1], q->handles[2] };
1999  CN::SideNumber( type, conn, ec, 4, 2, side, sense, offset );
2000  CN::SubEntityNodeIndices( type, len, 2, side, face_type, clen, indices );
2001  assert( MBQUAD == face_type );
2002  for( int k = 4; k < clen; ++k )
2003  skin_verts->insert( conn[indices[k]] );
2004  }
2005  }
2006  }
2007  }
2008 
2009  // If user requested that we pass back the list of 2D elements
2010  // representing the skin of the mesh...
2011  if( find_faces )
2012  {
2013  // Search list of adjacent faces for any that are on the skin
2014  adj.clear();
2015  rval = thisMB->get_adjacencies( &*it, 1, 2, false, adj );
2016  if( MB_SUCCESS != rval ) return rval;
2017 
2018  for( i = adj.begin(); i != adj.end(); ++i )
2019  {
2020  rval = thisMB->get_connectivity( *i, conn, len, true );
2021  if( MB_SUCCESS != rval ) return rval;
2022  const int idx2 = std::find( conn, conn + len, *it ) - conn;
2023  if( idx2 >= len )
2024  {
2025  assert( printed_nonconformal_ho_warning );
2026  continue;
2027  }
2028 
2029  // Note that the order of the terms in the if statements below
2030  // is important. We want to unmark any existing skin faces even
2031  // if we aren't returning them. Otherwise we'll end up creating
2032  // duplicates if create_faces is true.
2033  if( 3 == len )
2034  {
2035  if( adj_tris.find_and_unmark( conn, idx2, elem ) && skin_faces )
2036  {
2037  if( reversed_faces && face_reversed( elem, conn, MBTRI ) )
2038  reversed_faces->insert( *i );
2039  else
2040  skin_faces->insert( *i );
2041  }
2042  }
2043  else if( 4 == len )
2044  {
2045  if( adj_quads.find_and_unmark( conn, idx2, elem ) && skin_faces )
2046  {
2047  if( reversed_faces && face_reversed( elem, conn, MBQUAD ) )
2048  reversed_faces->insert( *i );
2049  else
2050  skin_faces->insert( *i );
2051  }
2052  }
2053  else
2054  {
2055  if( adj_poly.find_and_unmark( &*i, 1, elem ) && skin_faces ) skin_faces->insert( *i );
2056  }
2057  }
2058  }
2059 
2060  // If user does not want use to create new faces representing
2061  // sides for which there is currently no explicit element,
2062  // skip the remaining code and advance the outer loop to the
2063  // next vertex.
2064  if( !create_faces ) continue;
2065 
2066  // Polyhedra always have explicitly defined faces, so
2067  // there is no way we could need to create such a face.
2068  assert( 0 == adj_poly.num_skin() );
2069 
2070  // Create any skin tris that don't exist
2071  if( adj_tris.num_skin() )
2072  {
2073  for( AdjSides< 3 >::const_iterator t = adj_tris.begin(); t != adj_tris.end(); ++t )
2074  {
2075  if( t->skin() )
2076  {
2077  EntityHandle tri, c[3] = { *it, t->handles[0], t->handles[1] };
2078  rval = create_side( this_set, t->adj_elem, MBTRI, c, tri );
2079  if( MB_SUCCESS != rval ) return rval;
2080  if( skin_faces ) skin_faces->insert( tri );
2081  }
2082  }
2083  }
2084 
2085  // Create any skin quads that don't exist
2086  if( adj_quads.num_skin() )
2087  {
2088  for( AdjSides< 4 >::const_iterator q = adj_quads.begin(); q != adj_quads.end(); ++q )
2089  {
2090  if( q->skin() )
2091  {
2092  EntityHandle quad, c[4] = { *it, q->handles[0], q->handles[1], q->handles[2] };
2093  rval = create_side( this_set, q->adj_elem, MBQUAD, c, quad );
2094  if( MB_SUCCESS != rval ) return rval;
2095  if( skin_faces ) skin_faces->insert( quad );
2096  }
2097  }
2098  }
2099  } // end for each vertex
2100 
2101  return MB_SUCCESS;
2102 }
2103 
2104 } // namespace moab