Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
MeshTopoUtil.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 WIN32
17 #pragma warning( disable : 4786 )
18 #endif
19 
20 #include "moab/MeshTopoUtil.hpp"
21 #include "moab/Range.hpp"
22 #include "Internals.hpp"
23 #include "moab/Interface.hpp"
24 #include "moab/CN.hpp"
25 
26 #include <cassert>
27 
28 #define RR \
29  { \
30  if( MB_SUCCESS != result ) return result; \
31  }
32 
33 namespace moab
34 {
35 
36 //! generate all the AEntities bounding the vertices
38 {
39  Range out_range;
40  ErrorCode result;
41  result = mbImpl->get_adjacencies( vertices, 1, true, out_range, Interface::UNION );
42  if( MB_SUCCESS != result ) return result;
43  out_range.clear();
44  result = mbImpl->get_adjacencies( vertices, 2, true, out_range, Interface::UNION );
45  if( MB_SUCCESS != result ) return result;
46  out_range.clear();
47  result = mbImpl->get_adjacencies( vertices, 3, true, out_range, Interface::UNION );
48 
49  return result;
50 }
51 
52 //! given an entity, get its average position (avg vertex locations)
54 {
55  std::vector< EntityHandle > ent_vec;
56  std::copy( entities.begin(), entities.end(), std::back_inserter( ent_vec ) );
57  return get_average_position( &ent_vec[0], ent_vec.size(), avg_position );
58 }
59 
60 //! given an entity, get its average position (avg vertex locations)
62  const int num_entities,
63  double* avg_position )
64 {
65  double dum_pos[3];
66  avg_position[0] = avg_position[1] = avg_position[2] = 0.0;
67 
68  Range connect;
69  ErrorCode result = mbImpl->get_adjacencies( entities, num_entities, 0, false, connect, Interface::UNION );
70  if( MB_SUCCESS != result ) return result;
71 
72  if( connect.empty() ) return MB_FAILURE;
73 
74  for( Range::iterator rit = connect.begin(); rit != connect.end(); ++rit )
75  {
76  result = mbImpl->get_coords( &( *rit ), 1, dum_pos );
77  if( MB_SUCCESS != result ) return result;
78  avg_position[0] += dum_pos[0];
79  avg_position[1] += dum_pos[1];
80  avg_position[2] += dum_pos[2];
81  }
82  avg_position[0] /= (double)connect.size();
83  avg_position[1] /= (double)connect.size();
84  avg_position[2] /= (double)connect.size();
85 
86  return MB_SUCCESS;
87 }
88 
89 //! given an entity, get its average position (avg vertex locations)
90 ErrorCode MeshTopoUtil::get_average_position( const EntityHandle entity, double* avg_position )
91 {
92  const EntityHandle* connect = NULL;
93  int num_connect = 0;
94  if( MBVERTEX == mbImpl->type_from_handle( entity ) ) return mbImpl->get_coords( &entity, 1, avg_position );
95 
96  ErrorCode result = mbImpl->get_connectivity( entity, connect, num_connect );
97  if( MB_SUCCESS != result ) return result;
98 
99  return get_average_position( connect, num_connect, avg_position );
100 }
101 
102 // given an entity, find the entities of next higher dimension around
103 // that entity, ordered by connection through next higher dimension entities;
104 // if any of the star entities is in only one entity of next higher dimension,
105 // on_boundary is returned true
107  std::vector< EntityHandle >& star_ents,
108  bool& bdy_entity,
109  const EntityHandle starting_star_entity,
110  std::vector< EntityHandle >* star_entities_dp2,
111  Range* star_candidates_dp2 )
112 {
113  // now start the traversal
114  bdy_entity = false;
115  EntityHandle last_entity = starting_star_entity, last_dp2 = 0, next_entity, next_dp2;
116  std::vector< EntityHandle > star_dp2;
117  ErrorCode result;
118  int center_dim = mbImpl->dimension_from_handle( star_center );
119 
120  Range tmp_candidates_dp2;
121  if( NULL != star_candidates_dp2 )
122  tmp_candidates_dp2 = *star_candidates_dp2;
123  else
124  {
125  result = mbImpl->get_adjacencies( &star_center, 1, center_dim + 2, false, tmp_candidates_dp2 );
126  if( MB_SUCCESS != result ) return result;
127  }
128 
129  do
130  {
131  // get the next star entity
132  result = star_next_entity( star_center, last_entity, last_dp2, &tmp_candidates_dp2, next_entity, next_dp2 );
133  if( MB_SUCCESS != result ) return result;
134 
135  // special case: if starting_star_entity isn't connected to any entities of next
136  // higher dimension, it's the only entity in the star; put it on the list and return
137  if( star_ents.empty() && next_entity == 0 && next_dp2 == 0 )
138  {
139  star_ents.push_back( last_entity );
140  bdy_entity = true;
141  return MB_SUCCESS;
142  }
143 
144  // if we're at a bdy and bdy_entity hasn't been set yet, we're at the
145  // first bdy; reverse the lists and start traversing in the other direction; but,
146  // pop the last star entity off the list and find it again, so that we properly
147  // check for next_dp2
148  if( 0 == next_dp2 && !bdy_entity )
149  {
150  star_ents.push_back( next_entity );
151  bdy_entity = true;
152  std::reverse( star_ents.begin(), star_ents.end() );
153  star_ents.pop_back();
154  last_entity = star_ents.back();
155  if( !star_dp2.empty() )
156  {
157  std::reverse( star_dp2.begin(), star_dp2.end() );
158  last_dp2 = star_dp2.back();
159  }
160  }
161  // else if we're not on the bdy and next_entity is already in star, that means
162  // we've come all the way around; don't put next_entity on list again, and
163  // zero out last_dp2 to terminate while loop
164  else if( !bdy_entity && std::find( star_ents.begin(), star_ents.end(), next_entity ) != star_ents.end() &&
165  ( std::find( star_dp2.begin(), star_dp2.end(), next_dp2 ) != star_dp2.end() || !next_dp2 ) )
166  {
167  last_dp2 = 0;
168  }
169 
170  // else, just assign last entities seen and go on to next iteration
171  else
172  {
173  if( std::find( star_ents.begin(), star_ents.end(), next_entity ) == star_ents.end() )
174  star_ents.push_back( next_entity );
175  if( 0 != next_dp2 )
176  {
177  star_dp2.push_back( next_dp2 );
178  tmp_candidates_dp2.erase( next_dp2 );
179  }
180  last_entity = next_entity;
181  last_dp2 = next_dp2;
182  }
183  } while( 0 != last_dp2 );
184 
185  // copy over the star_dp2 list, if requested
186  if( NULL != star_entities_dp2 ) ( *star_entities_dp2 ).swap( star_dp2 );
187 
188  return MB_SUCCESS;
189 }
190 
192  const EntityHandle last_entity,
193  const EntityHandle last_dp1,
194  Range* star_candidates_dp1,
195  EntityHandle& next_entity,
196  EntityHandle& next_dp1 )
197 {
198  // given a star_center, a last_entity (whose dimension should be 1 greater than center)
199  // and last_dp1 (dimension 2 higher than center), returns the next star entity across
200  // last_dp1, and the next dp1 entity sharing next_entity; if star_candidates is non-empty,
201  // star must come from those
202  Range from_ents, to_ents;
203  from_ents.insert( star_center );
204  if( 0 != last_dp1 ) from_ents.insert( last_dp1 );
205 
206  int dim = mbImpl->dimension_from_handle( star_center );
207 
208  ErrorCode result = mbImpl->get_adjacencies( from_ents, dim + 1, true, to_ents );
209  if( MB_SUCCESS != result ) return result;
210 
211  // remove last_entity from result, and should only have 1 left, if any
212  if( 0 != last_entity ) to_ents.erase( last_entity );
213 
214  // if no last_dp1, contents of to_ents should share dp1-dimensional entity with last_entity
215  if( 0 != last_entity && 0 == last_dp1 )
216  {
217  Range tmp_to_ents;
218  for( Range::iterator rit = to_ents.begin(); rit != to_ents.end(); ++rit )
219  {
220  if( 0 != common_entity( last_entity, *rit, dim + 2 ) ) tmp_to_ents.insert( *rit );
221  }
222  to_ents = tmp_to_ents;
223  }
224 
225  if( 0 == last_dp1 && to_ents.size() > 1 && NULL != star_candidates_dp1 && !star_candidates_dp1->empty() )
226  {
227  // if we have a choice of to_ents and no previous dp1 and there are dp1 candidates,
228  // the one we choose needs to be adjacent to one of the candidates
229  result = mbImpl->get_adjacencies( *star_candidates_dp1, dim + 1, true, from_ents, Interface::UNION );
230  if( MB_SUCCESS != result ) return result;
231  to_ents = intersect( to_ents, from_ents );
232  }
233 
234  if( !to_ents.empty() )
235  next_entity = *to_ents.begin();
236  else
237  {
238  next_entity = 0;
239  next_dp1 = 0;
240  return MB_SUCCESS;
241  }
242 
243  // get next_dp1
244  if( 0 != star_candidates_dp1 )
245  to_ents = *star_candidates_dp1;
246  else
247  to_ents.clear();
248 
249  result = mbImpl->get_adjacencies( &next_entity, 1, dim + 2, true, to_ents );
250  if( MB_SUCCESS != result ) return result;
251 
252  // can't be last one
253  if( 0 != last_dp1 ) to_ents.erase( last_dp1 );
254 
255  if( !to_ents.empty() ) next_dp1 = *to_ents.begin();
256 
257  // could be zero, means we're at bdy
258  else
259  next_dp1 = 0;
260 
261  return MB_SUCCESS;
262 }
263 
265  std::vector< std::vector< EntityHandle > >& stars,
266  std::vector< bool >* bdy_flags,
267  std::vector< std::vector< EntityHandle > >* dp2_stars )
268 {
269  // Get a series of (d+1)-dimensional stars around a d-dimensional entity, such that
270  // each star is on a (d+2)-manifold containing the d-dimensional entity; each star
271  // is either open or closed, and also defines a (d+2)-star whose entities are bounded by
272  // (d+1)-entities on the star and on the (d+2)-manifold
273  //
274  // Algorithm:
275  // get the (d+2)-manifold entities; for d=1 / d+2=3, just assume all connected elements, since
276  // we don't do 4d yet
277  // get intersection of (d+1)-entities adjacent to star entity and union of (d+1)-entities
278  // adjacent to (d+2)-manifold entities; these will be the entities in the star
279  // while (d+1)-entities
280  // remove (d+1)-entity from (d+1)-entities
281  // get the (d+1)-star and (d+2)-star around that (d+1)-entity (using star_entities)
282  // save that star to the star list, and the bdy flag and (d+2)-star if requested
283  // remove (d+2)-entities from the (d+2)-manifold entities
284  // remove (d+1)-entities from the (d+1)-entities
285  // (end while)
286 
287  int this_dim = mbImpl->dimension_from_handle( star_entity );
288  if( 3 <= this_dim || 0 > this_dim ) return MB_FAILURE;
289 
290  // get the (d+2)-manifold entities; for d=1 / d+2=3, just assume all connected elements, since
291  // we don't do 4d yet
292  Range dp2_manifold;
293  ErrorCode result = get_manifold( star_entity, this_dim + 2, dp2_manifold );
294  if( MB_SUCCESS != result ) return result;
295 
296  // get intersection of (d+1)-entities adjacent to star and union of (d+1)-entities
297  // adjacent to (d+2)-manifold entities; also add manifold (d+1)-entities, to catch
298  // any not connected to (d+2)-entities
299  Range dp1_manifold;
300  result = mbImpl->get_adjacencies( dp2_manifold, this_dim + 1, false, dp1_manifold, Interface::UNION );
301  if( MB_SUCCESS != result ) return result;
302 
303  result = mbImpl->get_adjacencies( &star_entity, 1, this_dim + 1, false, dp1_manifold );
304  if( MB_SUCCESS != result ) return result;
305 
306  result = get_manifold( star_entity, this_dim + 1, dp1_manifold );
307  if( MB_SUCCESS != result ) return result;
308 
309  // while (d+1)-entities
310  while( !dp1_manifold.empty() )
311  {
312 
313  // get (d+1)-entity from (d+1)-entities (don't remove it until after star,
314  // since the star entities must come from dp1_manifold)
315  EntityHandle this_ent = *dp1_manifold.begin();
316 
317  // get the (d+1)-star and (d+2)-star around that (d+1)-entity (using star_entities)
318  std::vector< EntityHandle > this_star_dp1, this_star_dp2;
319  bool on_bdy;
320  result = star_entities( star_entity, this_star_dp1, on_bdy, this_ent, &this_star_dp2, &dp2_manifold );
321  if( MB_SUCCESS != result ) return result;
322 
323  // if there's no star entities, it must mean this_ent isn't bounded by any dp2
324  // entities (wasn't put into star in star_entities 'cuz we're passing in non-null
325  // dp2_manifold above); put it in
326  if( this_star_dp1.empty() )
327  {
328  Range dum_range;
329  result = mbImpl->get_adjacencies( &this_ent, 1, this_dim + 2, false, dum_range );
330  if( MB_SUCCESS != result ) return result;
331  if( dum_range.empty() ) this_star_dp1.push_back( this_ent );
332  }
333 
334  // now we can remove it
335  dp1_manifold.erase( dp1_manifold.begin() );
336 
337  // save that star to the star list, and the bdy flag and (d+2)-star if requested
338  if( !this_star_dp1.empty() )
339  {
340  stars.push_back( this_star_dp1 );
341  if( NULL != bdy_flags ) bdy_flags->push_back( on_bdy );
342  if( NULL != dp2_stars ) dp2_stars->push_back( this_star_dp2 );
343  }
344 
345  // remove (d+2)-entities from the (d+2)-manifold entities
346  for( std::vector< EntityHandle >::iterator vit = this_star_dp2.begin(); vit != this_star_dp2.end(); ++vit )
347  dp2_manifold.erase( *vit );
348 
349  // remove (d+1)-entities from the (d+1)-entities
350  for( std::vector< EntityHandle >::iterator vit = this_star_dp1.begin(); vit != this_star_dp1.end(); ++vit )
351  dp1_manifold.erase( *vit );
352 
353  // (end while)
354  }
355 
356  // check for leftover dp2 manifold entities, these should be in one of the
357  // stars
358  if( !dp2_manifold.empty() )
359  {
360  for( Range::iterator rit = dp2_manifold.begin(); rit != dp2_manifold.end(); ++rit )
361  {
362  }
363  }
364 
365  return MB_SUCCESS;
366 }
367 
368 //! get (target_dim)-dimensional manifold entities connected to star_entity; that is,
369 //! the entities with <= 1 connected (target_dim+2)-dimensional adjacent entities;
370 //! for target_dim=3, just return all of them
371 //! just insert into the list, w/o clearing manifold list first
372 ErrorCode MeshTopoUtil::get_manifold( const EntityHandle star_entity, const int target_dim, Range& manifold )
373 {
374  // get all the entities of target dimension connected to star
375  Range tmp_range;
376  ErrorCode result = mbImpl->get_adjacencies( &star_entity, 1, target_dim, false, tmp_range );
377  if( MB_SUCCESS != result ) return result;
378 
379  // now save the ones which are (target_dim+1)-dimensional manifold;
380  // for target_dim=3, just return whole range, since we don't do 4d
381  if( target_dim == 3 )
382  {
383  manifold.merge( tmp_range );
384  return MB_SUCCESS;
385  }
386 
387  for( Range::iterator rit = tmp_range.begin(); rit != tmp_range.end(); ++rit )
388  {
389  Range dum_range;
390  // get (target_dim+1)-dimensional entities
391  result = mbImpl->get_adjacencies( &( *rit ), 1, target_dim + 1, false, dum_range );
392  if( MB_SUCCESS != result ) return result;
393 
394  // if there are only 1 or zero, add to manifold list
395  if( 1 >= dum_range.size() ) manifold.insert( *rit );
396  }
397 
398  return MB_SUCCESS;
399 }
400 
401 //! get "bridge" or "2nd order" adjacencies, going through dimension bridge_dim
403  int bridge_dim,
404  int to_dim,
405  Range& to_ents,
406  int num_layers )
407 {
408  Range bridge_ents, accum_layers, new_toents( from_entities );
409  ErrorCode result;
410  if( 0 == num_layers || from_entities.empty() ) return MB_FAILURE;
411 
412  // for each layer, get bridge-adj entities and accumulate
413  for( int nl = 0; nl < num_layers; nl++ )
414  {
415  Range new_bridges;
416  // get bridge ents
417  result = mbImpl->get_adjacencies( new_toents, bridge_dim, true, new_bridges, Interface::UNION );
418  if( MB_SUCCESS != result ) return result;
419 
420  // get to_dim adjacencies, merge into to_ents
421  Range new_layer;
422  if( -1 == to_dim )
423  {
424  result = mbImpl->get_adjacencies( new_bridges, 3, false, new_layer, Interface::UNION );
425  if( MB_SUCCESS != result ) return result;
426  for( int d = 2; d >= 1; d-- )
427  {
428  result = mbImpl->get_adjacencies( to_ents, d, true, new_layer, Interface::UNION );
429  if( MB_SUCCESS != result ) return result;
430  }
431  }
432  else
433  {
434  result = mbImpl->get_adjacencies( new_bridges, to_dim, false, new_layer, Interface::UNION );
435  if( MB_SUCCESS != result ) return result;
436  }
437 
438  // subtract last_toents to get new_toents
439  accum_layers.merge( new_layer );
440  if( nl < num_layers - 1 ) new_toents = subtract( new_layer, new_toents );
441  }
442  to_ents.merge( accum_layers );
443 
444  return MB_SUCCESS;
445 }
446 
447 //! get "bridge" or "2nd order" adjacencies, going through dimension bridge_dim
449  const int bridge_dim,
450  const int to_dim,
451  Range& to_adjs )
452 {
453  // get pointer to connectivity for this entity
454  const EntityHandle* connect;
455  int num_connect;
456  ErrorCode result = MB_SUCCESS;
457  EntityType from_type = TYPE_FROM_HANDLE( from_entity );
458  if( from_type == MBVERTEX )
459  {
460  connect = &from_entity;
461  num_connect = 1;
462  }
463  else
464  {
465  result = mbImpl->get_connectivity( from_entity, connect, num_connect );
466  if( MB_SUCCESS != result ) return result;
467  }
468 
469  if( from_type >= MBENTITYSET ) return MB_FAILURE;
470 
471  int from_dim = CN::Dimension( from_type );
472 
473  Range to_ents;
474 
475  if( bridge_dim < from_dim )
476  {
477  // looping over each sub-entity of dimension bridge_dim...
478  if( MBPOLYGON == from_type )
479  {
480  for( int i = 0; i < num_connect; i++ )
481  {
482  // loop over edges, and get the vertices
483  EntityHandle verts_on_edge[2] = { connect[i], connect[( i + 1 ) % num_connect] };
484  to_ents.clear();
485  ErrorCode tmp_result =
486  mbImpl->get_adjacencies( verts_on_edge, 2, to_dim, false, to_ents, Interface::INTERSECT );
487  if( MB_SUCCESS != tmp_result ) result = tmp_result;
488  to_adjs.merge( to_ents );
489  }
490  }
491  else
492  {
493  EntityHandle bridge_verts[MAX_SUB_ENTITIES];
494  int bridge_indices[MAX_SUB_ENTITIES];
495  for( int i = 0; i < CN::NumSubEntities( from_type, bridge_dim ); i++ )
496  {
497 
498  // get the vertices making up this sub-entity
499  int num_bridge_verts = CN::VerticesPerEntity( CN::SubEntityType( from_type, bridge_dim, i ) );
500  assert( num_bridge_verts >= 0 && num_bridge_verts <= MAX_SUB_ENTITIES );
501  CN::SubEntityVertexIndices( from_type, bridge_dim, i, bridge_indices );
502  for( int j = 0; j < num_bridge_verts; ++j )
503  {
504  if( bridge_indices[j] >= 0 && bridge_indices[j] < num_connect )
505  bridge_verts[j] = connect[bridge_indices[j]];
506  else
507  bridge_verts[j] = 0;
508  }
509  // CN::SubEntityConn(connect, from_type, bridge_dim, i, &bridge_verts[0],
510  // num_bridge_verts);
511 
512  // get the to_dim entities adjacent
513  to_ents.clear();
514  ErrorCode tmp_result = mbImpl->get_adjacencies( bridge_verts, num_bridge_verts, to_dim, false, to_ents,
516  if( MB_SUCCESS != tmp_result ) result = tmp_result;
517 
518  to_adjs.merge( to_ents );
519  }
520  }
521  }
522 
523  // now get the direct ones too, or only in the case where we're
524  // going to higher dimension for bridge
525  Range bridge_ents, tmp_ents;
526  tmp_ents.insert( from_entity );
527  ErrorCode tmp_result = mbImpl->get_adjacencies( tmp_ents, bridge_dim, false, bridge_ents, Interface::UNION );
528  if( MB_SUCCESS != tmp_result ) return tmp_result;
529 
530  tmp_result = mbImpl->get_adjacencies( bridge_ents, to_dim, false, to_adjs, Interface::UNION );
531  if( MB_SUCCESS != tmp_result ) return tmp_result;
532 
533  // if to_dimension is same as that of from_entity, make sure from_entity isn't
534  // in list
535  if( to_dim == from_dim ) to_adjs.erase( from_entity );
536 
537  return result;
538 }
539 
540 //! return a common entity of the specified dimension, or 0 if there isn't one
542 {
543  Range tmp_range, tmp_range2;
544  tmp_range.insert( ent1 );
545  tmp_range.insert( ent2 );
546  ErrorCode result = mbImpl->get_adjacencies( tmp_range, dim, false, tmp_range2 );
547  if( MB_SUCCESS != result || tmp_range2.empty() )
548  return 0;
549  else
550  return *tmp_range2.begin();
551 }
552 
553 //! return the opposite side entity given a parent and bounding entity.
554 //! This function is only defined for certain types of parent/child types;
555 //! See CN.hpp::OppositeSide for details.
556 //!
557 //! \param parent The parent element
558 //! \param child The child element
559 //! \param opposite_element The index of the opposite element
561  const EntityHandle child,
562  EntityHandle& opposite_element )
563 {
564  // get the side no.
565  int side_no, offset, sense;
566  ErrorCode result = mbImpl->side_number( parent, child, side_no, offset, sense );
567  if( MB_SUCCESS != result ) return result;
568 
569  // get the child index from CN
570  int opposite_index, opposite_dim;
571  int status = CN::OppositeSide( mbImpl->type_from_handle( parent ), side_no, mbImpl->dimension_from_handle( child ),
572  opposite_index, opposite_dim );
573  if( 0 != status ) return MB_FAILURE;
574 
575  // now get the side element from MOAB
576  result = mbImpl->side_element( parent, opposite_dim, opposite_index, opposite_element );
577  if( MB_SUCCESS != result ) return result;
578 
579  return MB_SUCCESS;
580 }
581 
583 {
584  Range tmp_range, *tmp_ptr_fill_entity;
585  if( NULL != fill_entities )
586  tmp_ptr_fill_entity = &tmp_range;
587  else
588  tmp_ptr_fill_entity = NULL;
589 
590  for( Range::iterator rit = entities.begin(); rit != entities.end(); ++rit )
591  {
592  EntityHandle new_entity;
593  if( NULL != tmp_ptr_fill_entity ) tmp_ptr_fill_entity->clear();
594 
595  EntityHandle this_ent = *rit;
596  ErrorCode result = split_entities_manifold( &this_ent, 1, &new_entity, tmp_ptr_fill_entity );
597  if( MB_SUCCESS != result ) return result;
598 
599  new_entities.insert( new_entity );
600  if( NULL != fill_entities ) fill_entities->merge( *tmp_ptr_fill_entity );
601  }
602 
603  return MB_SUCCESS;
604 }
605 
607  const int num_entities,
608  EntityHandle* new_entities,
609  Range* fill_entities,
610  EntityHandle* gowith_ents )
611 {
612  // split entities by duplicating them; splitting manifold means that there is at
613  // most two higher-dimension entities bounded by a given entity; after split, the
614  // new entity bounds one and the original entity bounds the other
615 
616 #define ITERATE_RANGE( range, it ) for( Range::iterator it = ( range ).begin(); ( it ) != ( range ).end(); ++( it ) )
617 #define GET_CONNECT_DECL( ent, connect, num_connect ) \
618  const EntityHandle* connect = NULL; \
619  int num_connect = 0; \
620  { \
621  ErrorCode connect_result = mbImpl->get_connectivity( ent, connect, num_connect ); \
622  if( MB_SUCCESS != connect_result ) return connect_result; \
623  }
624 #define GET_CONNECT( ent, connect, num_connect ) \
625  { \
626  ErrorCode connect_result = mbImpl->get_connectivity( ent, connect, num_connect ); \
627  if( MB_SUCCESS != connect_result ) return connect_result; \
628  }
629 #define TC \
630  if( MB_SUCCESS != tmp_result ) \
631  { \
632  result = tmp_result; \
633  continue; \
634  }
635 
636  ErrorCode result = MB_SUCCESS;
637  for( int i = 0; i < num_entities; i++ )
638  {
639  ErrorCode tmp_result;
640 
641  // get original higher-dimensional bounding entities
642  Range up_adjs[4];
643  // can only do a split_manifold if there are at most 2 entities of each
644  // higher dimension; otherwise it's a split non-manifold
645  bool valid_up_adjs = true;
646  for( int dim = 1; dim <= 3; dim++ )
647  {
648  tmp_result = mbImpl->get_adjacencies( entities + i, 1, dim, false, up_adjs[dim] );
649  TC;
650  if( dim > CN::Dimension( TYPE_FROM_HANDLE( entities[i] ) ) && up_adjs[dim].size() > 2 )
651  {
652  valid_up_adjs = false;
653  break;
654  }
655  }
656  if( !valid_up_adjs ) return MB_FAILURE;
657 
658  // ok to split; create the new entity, with connectivity of the original
659  GET_CONNECT_DECL( entities[i], connect, num_connect );
660  EntityHandle new_entity;
661  result = mbImpl->create_element( mbImpl->type_from_handle( entities[i] ), connect, num_connect, new_entity );
662  TC;
663 
664  // by definition, new entity and original will be equivalent; need to add explicit
665  // adjs to distinguish them; don't need to check if there's already one there,
666  // 'cuz add_adjacency does that for us
667  for( int dim = 1; dim <= 3; dim++ )
668  {
669  if( up_adjs[dim].empty() || dim == CN::Dimension( TYPE_FROM_HANDLE( entities[i] ) ) ) continue;
670 
671  if( dim < CN::Dimension( TYPE_FROM_HANDLE( entities[i] ) ) )
672  {
673  // adjacencies from other entities to this one; if any of those are equivalent
674  // entities, need to make explicit adjacency to new entity too
675  for( Range::iterator rit = up_adjs[dim].begin(); rit != up_adjs[dim].end(); ++rit )
676  {
677  if( equivalent_entities( *rit ) ) result = mbImpl->add_adjacencies( *rit, &new_entity, 1, false );
678  }
679  }
680  else
681  {
682 
683  // get the two up-elements
684  EntityHandle up_elem1 = *( up_adjs[dim].begin() ),
685  up_elem2 = ( up_adjs[dim].size() > 1 ? *( up_adjs[dim].rbegin() ) : 0 );
686 
687  // if two, and a gowith entity was input, make sure the new entity goes with
688  // that one
689  if( gowith_ents && up_elem2 && gowith_ents[i] != up_elem1 && gowith_ents[i] == up_elem2 )
690  {
691  EntityHandle tmp_elem = up_elem1;
692  up_elem1 = up_elem2;
693  up_elem2 = tmp_elem;
694  }
695 
696  mbImpl->remove_adjacencies( entities[i], &up_elem1, 1 );
697  // (ok if there's an error, that just means there wasn't an explicit adj)
698 
699  tmp_result = mbImpl->add_adjacencies( new_entity, &up_elem1, 1, false );
700  TC;
701  if( !up_elem2 ) continue;
702 
703  // add adj to other up_adj
704  tmp_result = mbImpl->add_adjacencies( entities[i], &up_elem2, 1, false );
705  TC;
706  }
707  }
708 
709  // if we're asked to build a next-higher-dimension object, do so
710  EntityHandle fill_entity = 0;
711  EntityHandle tmp_ents[2];
712  if( NULL != fill_entities )
713  {
714  // how to do this depends on dimension
715  switch( CN::Dimension( TYPE_FROM_HANDLE( entities[i] ) ) )
716  {
717  case 0:
718  tmp_ents[0] = entities[i];
719  tmp_ents[1] = new_entity;
720  tmp_result = mbImpl->create_element( MBEDGE, tmp_ents, 2, fill_entity );
721  TC;
722  break;
723  case 1:
724  tmp_result = mbImpl->create_element( MBPOLYGON, connect, 2, fill_entity );
725  TC;
726  // need to create explicit adj in this case
727  tmp_result = mbImpl->add_adjacencies( entities[i], &fill_entity, 1, false );
728  TC;
729  tmp_result = mbImpl->add_adjacencies( new_entity, &fill_entity, 1, false );
730  TC;
731  break;
732  case 2:
733  tmp_ents[0] = entities[i];
734  tmp_ents[1] = new_entity;
735  tmp_result = mbImpl->create_element( MBPOLYHEDRON, tmp_ents, 2, fill_entity );
736  TC;
737  break;
738  }
739  if( 0 == fill_entity )
740  {
741  result = MB_FAILURE;
742  continue;
743  }
744  fill_entities->insert( fill_entity );
745  }
746 
747  new_entities[i] = new_entity;
748 
749  } // end for over input entities
750 
751  return result;
752 }
753 
755  Range& old_adjs,
756  Range& new_adjs,
757  EntityHandle& new_entity )
758 {
759  // split an entity into two entities; new entity gets explicit adj to new_adjs,
760  // old to old_adjs
761 
762  // make new entities and add adjacencies
763  // create the new entity
764  EntityType split_type = mbImpl->type_from_handle( split_ent );
765 
766  ErrorCode result;
767  if( MBVERTEX == split_type )
768  {
769  double coords[3];
770  result = mbImpl->get_coords( &split_ent, 1, coords );RR;
771  result = mbImpl->create_vertex( coords, new_entity );RR;
772  }
773  else
774  {
775  const EntityHandle* connect;
776  int num_connect;
777  result = mbImpl->get_connectivity( split_ent, connect, num_connect );RR;
778  result = mbImpl->create_element( split_type, connect, num_connect, new_entity );RR;
779 
780  // remove any explicit adjacencies between new_adjs and split entity
781  for( Range::iterator rit = new_adjs.begin(); rit != new_adjs.end(); ++rit )
782  mbImpl->remove_adjacencies( split_ent, &( *rit ), 1 );
783  }
784 
785  if( MBVERTEX != split_type )
786  {
787  // add adj's between new_adjs & new entity, old_adjs & split_entity
788  for( Range::iterator rit = new_adjs.begin(); rit != new_adjs.end(); ++rit )
789  mbImpl->add_adjacencies( new_entity, &( *rit ), 1, true );
790  for( Range::iterator rit = old_adjs.begin(); rit != old_adjs.end(); ++rit )
791  mbImpl->add_adjacencies( split_ent, &( *rit ), 1, true );
792  }
793  else if( split_ent != new_entity )
794  {
795  // in addition to explicit adjs, need to check if vertex is part of any
796  // other entities, and check those entities against ents in old and new adjs
797  Range other_adjs;
798  for( int i = 1; i < 4; i++ )
799  {
800  result = mbImpl->get_adjacencies( &split_ent, 1, i, false, other_adjs, Interface::UNION );RR;
801  }
802  other_adjs = subtract( other_adjs, old_adjs );
803  other_adjs = subtract( other_adjs, new_adjs );
804  for( Range::iterator rit1 = other_adjs.begin(); rit1 != other_adjs.end(); ++rit1 )
805  {
806  // find an adjacent lower-dimensional entity in old_ or new_ adjs
807  bool found = false;
808  for( Range::iterator rit2 = old_adjs.begin(); rit2 != old_adjs.end(); ++rit2 )
809  {
810  if( mbImpl->dimension_from_handle( *rit1 ) != mbImpl->dimension_from_handle( *rit2 ) &&
811  common_entity( *rit1, *rit2, mbImpl->dimension_from_handle( *rit1 ) ) )
812  {
813  found = true;
814  old_adjs.insert( *rit1 );
815  break;
816  }
817  }
818  if( found ) continue;
819  for( Range::iterator rit2 = new_adjs.begin(); rit2 != new_adjs.end(); ++rit2 )
820  {
821  if( mbImpl->dimension_from_handle( *rit1 ) != mbImpl->dimension_from_handle( *rit2 ) &&
822  common_entity( *rit1, *rit2, mbImpl->dimension_from_handle( *rit1 ) ) )
823  {
824  found = true;
825  new_adjs.insert( *rit1 );
826  break;
827  }
828  }
829  if( !found ) return MB_FAILURE;
830  }
831 
832  // instead of adjs replace in connectivity
833  std::vector< EntityHandle > connect;
834  for( Range::iterator rit = new_adjs.begin(); rit != new_adjs.end(); ++rit )
835  {
836  connect.clear();
837  result = mbImpl->get_connectivity( &( *rit ), 1, connect );RR;
838  std::replace( connect.begin(), connect.end(), split_ent, new_entity );
839  result = mbImpl->set_connectivity( *rit, &connect[0], connect.size() );RR;
840  }
841  }
842 
843  return result;
844 
845  /*
846 
847  Commented out for now, because I decided to do a different implementation
848  for the sake of brevity. However, I still think this function is the right
849  way to do it, if I ever get the time. Sigh.
850 
851  // split entity d, producing entity nd; generates various new entities,
852  // see algorithm description in notes from 2/25/05
853  const EntityHandle split_types = {MBEDGE, MBPOLYGON, MBPOLYHEDRON};
854  ErrorCode result = MB_SUCCESS;
855  const int dim = CN::Dimension(TYPE_FROM_HANDLE(d));
856  MeshTopoUtil mtu(this);
857 
858  // get all (d+2)-, (d+1)-cells connected to d
859  Range dp2s, dp1s, dp1s_manif, dp2s_manif;
860  result = get_adjacencies(&d, 1, dim+2, false, dp2s); RR;
861  result = get_adjacencies(&d, 1, dim+1, false, dp1s); RR;
862 
863  // also get (d+1)-cells connected to d which are manifold
864  get_manifold_dp1s(d, dp1s_manif);
865  get_manifold_dp2s(d, dp2s_manif);
866 
867  // make new cell nd, then ndp1
868  result = copy_entity(d, nd); RR;
869  EntityHandle tmp_connect[] = {d, nd};
870  EntityHandle ndp1;
871  result = create_element(split_types[dim],
872  tmp_connect, 2, ndp1); RR;
873 
874  // modify (d+2)-cells, depending on what type they are
875  ITERATE_RANGE(dp2s, dp2) {
876  // first, get number of connected manifold (d+1)-entities
877  Range tmp_range, tmp_range2(dp1s_manif);
878  tmp_range.insert(*dp2);
879  tmp_range.insert(d);
880  tmp_result = get_adjacencies(tmp_range, 1, false, tmp_range2); TC;
881  EntityHandle ndp2;
882 
883  // a. manif (d+1)-cells is zero
884  if (tmp_range2.empty()) {
885  // construct new (d+1)-cell
886  EntityHandle ndp1a;
887  EntityHandle tmp_result = create_element(split_types[dim],
888  tmp_connect, 2, ndp1a); TC;
889  // now make new (d+2)-cell
890  EntityHandle tmp_connect2[] = {ndp1, ndp1a};
891  tmp_result = create_element(split_types[dim+1],
892  tmp_connect2, 2, ndp2); TC;
893  // need to add explicit adjacencies, since by definition ndp1, ndp1a will be equivalent
894  tmp_result = add_adjacencies(ndp1a, &dp2, 1, false); TC;
895  tmp_result = add_adjacencies(ndp1a, &ndp2, 1, false); TC;
896  tmp_result = add_adjacencies(ndp1, &ndp2, 1, false); TC;
897 
898  // now insert nd into connectivity of dp2, right after d if dim < 1
899  std::vector<EntityHandle> connect;
900  tmp_result = get_connectivity(&dp2, 1, connect); TC;
901  if (dim < 1) {
902  std::vector<EntityHandle>::iterator vit = std::find(connect.begin(), connect.end(), d);
903  if (vit == connect.end()) {
904  result = MB_FAILURE;
905  continue;
906  }
907  connect.insert(vit, nd);
908  }
909  else
910  connect.push_back(nd);
911  tmp_result = set_connectivity(dp2, connect); TC;
912 
913  // if dim < 1, need to add explicit adj from ndp2 to higher-dim ents, since it'll
914  // be equiv to other dp2 entities
915  if (dim < 1) {
916  Range tmp_dp3s;
917  tmp_result = get_adjacencies(&dp2, 1, dim+3, false, tmp_dp3s); TC;
918  tmp_result = add_adjacencies(ndp2, tmp_dp3s, false); TC;
919  }
920  } // end if (tmp_range2.empty())
921 
922  // b. single manifold (d+1)-cell, which isn't adjacent to manifold (d+2)-cell
923  else if (tmp_range2.size() == 1) {
924  // b1. check validity, and skip if not valid
925 
926  // only change if not dp1-adjacent to manifold dp2cell; check that...
927  Range tmp_adjs(dp2s_manif);
928  tmp_result = get_adjacencies(&(*tmp_range2.begin()), 1, dim+2, false, tmp_adjs); TC;
929  if (!tmp_adjs.empty()) continue;
930 
931  EntityHandle dp1 = *tmp_range2.begin();
932 
933  // b2. make new (d+1)- and (d+2)-cell next to dp2
934 
935  // get the (d+2)-cell on the other side of dp1
936  tmp_result = get_adjacencies(&dp1, 1, dim+2, false, tmp_adjs); TC;
937  EntityHandle odp2 = *tmp_adjs.begin();
938  if (odp2 == dp2) odp2 = *tmp_adjs.rbegin();
939 
940  // get od, the d-cell on dp1_manif which isn't d
941  tmp_result = get_adjacencies(&dp1_manif, 1, dim, false, tmp_adjs); TC;
942  tmp_adjs.erase(d);
943  if (tmp_adjs.size() != 1) {
944  result = MB_FAILURE;
945  continue;
946  }
947  EntityHandle od = *tmp_adjs.begin();
948 
949  // make a new (d+1)-cell from od and nd
950  tmp_adjs.insert(nd);
951  tmp_result = create_element(split_types[1], tmp_adjs, ndp1a); TC;
952 
953  // construct new (d+2)-cell from dp1, ndp1, ndp1a
954  tmp_adjs.clear();
955  tmp_adjs.insert(dp1); tmp_adjs.insert(ndp1); tmp_adjs.insert(ndp1a);
956  tmp_result = create_element(split_types[2], tmp_adjs, ndp2); TC;
957 
958  // b3. replace d, dp1 in connect/adjs of odp2
959  std::vector<EntityHandle> connect;
960  tmp_result = get_connectivity(&odp2, 1, connect); TC;
961  if (dim == 0) {
962  *(std::find(connect.begin(), connect.end(), d)) = nd;
963  remove_adjacency(dp1, odp2);
964 
965 
966 
967  // if dp1 was explicitly adj to odp2, remove it
968  remove_adjacency
969 
970  ...
971 
972  */
973 }
974 
975 //! return whether entity is equivalent to any other of same type and same vertices;
976 //! if equivalent entity is found, it's returned in equiv_ents and return value is true,
977 //! false otherwise.
978 bool MeshTopoUtil::equivalent_entities( const EntityHandle entity, Range* equiv_ents )
979 {
980  const EntityHandle* connect = NULL;
981  int num_connect = 0;
982  ErrorCode result = mbImpl->get_connectivity( entity, connect, num_connect );
983  if( MB_SUCCESS != result ) return false;
984 
985  Range dum;
986  result = mbImpl->get_adjacencies( connect, num_connect, mbImpl->dimension_from_handle( entity ), false, dum );
987  dum.erase( entity );
988 
989  if( NULL != equiv_ents )
990  {
991  equiv_ents->swap( dum );
992  }
993 
994  if( !dum.empty() )
995  return true;
996  else
997  return false;
998 }
999 
1000 } // namespace moab