 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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
48 Skinner::~Skinner()
49 {
50 // delete the adjacency tag
51 }
52
53 ErrorCode Skinner::initialize()
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
73 Range entities;
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
79 thisMB->get_entities_by_type( 0, type, entities );
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
98 ErrorCode Skinner::deinitialize()
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
129 ErrorCode Skinner::add_adjacency( EntityHandle entity )
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 =
189 thisMB->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom_tag, MB_TAG_SPARSE | MB_TAG_CREAT );
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
233 ErrorCode Skinner::find_skin( const EntityHandle meshset,
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() )
253 this_core->a_entity_factory()->create_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
582 bool Skinner::connectivity_match( const EntityHandle* conn1,
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
632 ErrorCode Skinner::remove_adjacency( EntityHandle entity )
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
650 bool Skinner::entity_deletable( EntityHandle entity )
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
659 ErrorCode Skinner::classify_2d_boundary( const Range& boundary,
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
687 ErrorCode Skinner::classify_2d_boundary( const Range& boundary,
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
985 ErrorCode Skinner::find_skin( const EntityHandle this_set,
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
1013 ErrorCode Skinner::find_skin_vertices( const EntityHandle this_set,
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 >
1144 class AdjSides
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;
1253 const_iterator begin() const
1254 {
1255 return data.begin();
1256 }
1257 const_iterator end() const
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.
1432 ErrorCode Skinner::create_side( const EntityHandle this_set,
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.
1515 bool Skinner::edge_reversed( EntityHandle face, const EntityHandle* edge_ends )
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
1553 ErrorCode Skinner::find_skin_vertices_2D( const EntityHandle this_set,
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;
1586 EntityHandle face;
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
1768 ErrorCode Skinner::find_skin_vertices_3D( const EntityHandle this_set,
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
1.9.1.