Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
GeomTopoTool.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 #include "moab/GeomTopoTool.hpp"
17 #include "moab/GeomQueryTool.hpp"
19 #include "moab/Range.hpp"
20 #include "MBTagConventions.hpp"
21 #include "moab/Interface.hpp"
22 #include "moab/CN.hpp"
23 #include "moab/Skinner.hpp"
24 #include "Internals.hpp"
25 #include <iostream>
26 
27 namespace moab
28 {
29 
30 // Tag name used for saving sense of faces in volumes.
31 // We assume that the surface occurs in at most two volumes.
32 // Code will error out if more than two volumes per surface.
33 // The tag data is a pair of tag handles, representing the
34 // forward and reverse volumes, respectively. If a surface
35 // is non-manifold in a single volume, the same volume will
36 // be listed for both the forward and reverse slots.
37 const char GEOM_SENSE_2_TAG_NAME[] = "GEOM_SENSE_2";
38 
39 const char GEOM_SENSE_N_ENTS_TAG_NAME[] = "GEOM_SENSE_N_ENTS";
40 const char GEOM_SENSE_N_SENSES_TAG_NAME[] = "GEOM_SENSE_N_SENSES";
41 const char OBB_ROOT_TAG_NAME[] = "OBB_ROOT";
42 const char OBB_GSET_TAG_NAME[] = "OBB_GSET";
43 
44 const char IMPLICIT_COMPLEMENT_NAME[] = "impl_complement";
45 
47  bool find_geoments,
48  EntityHandle modelRootSet,
49  bool p_rootSets_vector,
50  bool restore_rootSets )
51  : mdbImpl( impl ), sense2Tag( 0 ), senseNEntsTag( 0 ), senseNSensesTag( 0 ), geomTag( 0 ), gidTag( 0 ),
52  obbRootTag( 0 ), obbGsetTag( 0 ), modelSet( modelRootSet ), updated( false ), setOffset( 0 ),
53  m_rootSets_vector( p_rootSets_vector ), oneVolRootSet( 0 )
54 {
55 
56  obbTree = new OrientedBoxTreeTool( impl, NULL, true );
57 
58  ErrorCode rval =
59  mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag, MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR_CONT( rval, "Error: Failed to create geometry dimension tag" );
60 
61  // global id tag is not really needed, but mbsize complains if we do not set it for
62  // geometry entities
64 
65  rval =
67 
68  rval = mdbImpl->tag_get_handle( OBB_ROOT_TAG_NAME, 1, MB_TYPE_HANDLE, obbRootTag, MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR_CONT( rval, "Error: Failed to create obb root tag" );
69 
70  rval = mdbImpl->tag_get_handle( OBB_GSET_TAG_NAME, 1, MB_TYPE_HANDLE, obbGsetTag, MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR_CONT( rval, "Error: Failed to create obb gset tag" );
71 
72  // set this value to zero for comparisons
74 
75  maxGlobalId[0] = maxGlobalId[1] = maxGlobalId[2] = maxGlobalId[3] = maxGlobalId[4] = 0;
76  if( find_geoments )
77  {
78  find_geomsets();
79  if( restore_rootSets )
80  {
81  rval = restore_obb_index();
82  if( MB_SUCCESS != rval )
83  {
84  rval = delete_all_obb_trees();MB_CHK_SET_ERR_CONT( rval, "Error: Failed to delete existing obb trees" );
85  rval = construct_obb_trees();MB_CHK_SET_ERR_CONT( rval, "Error: Failed to rebuild obb trees" );
86  }
87  }
88  }
89 }
90 
92 {
93  delete obbTree;
94 }
95 
97 {
98  ErrorCode result;
99  if( 0 == geomTag )
100  {
101  result = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag );MB_CHK_SET_ERR( result, "Failed to get the geometry dimension tag" );
102  }
103 
104  // check if the geo set belongs to this model
105  if( modelSet )
106  {
107  if( !mdbImpl->contains_entities( modelSet, &this_set, 1 ) )
108  {
109  // this g set does not belong to the current model
110  return -1;
111  }
112  }
113  // get the data for those tags
114  int dim;
115  result = mdbImpl->tag_get_data( geomTag, &this_set, 1, &dim );
116  if( MB_SUCCESS != result )
117  return -1;
118  else
119  return dim;
120 }
121 
123 {
124  ErrorCode result;
125  if( 0 == gidTag )
126  {
128  }
129 
130  // check if the geo set belongs to this model
131  if( modelSet )
132  {
133  if( !mdbImpl->contains_entities( modelSet, &this_set, 1 ) )
134  {
135  // this g set does not belong to the current model
136  return -1;
137  }
138  }
139 
140  // get the data for those tags
141  int id;
142  result = mdbImpl->tag_get_data( gidTag, &this_set, 1, &id );
143  if( MB_SUCCESS != result )
144  return -1;
145  else
146  return id;
147 }
148 
149 EntityHandle GeomTopoTool::entity_by_id( int dimension1, int id )
150 {
151  if( 0 > dimension1 || 3 < dimension1 )
152  {
153  MB_CHK_SET_ERR_CONT( MB_FAILURE, "Incorrect dimension provided" );
154  }
155  const Tag tags[] = { gidTag, geomTag };
156  const void* const vals[] = { &id, &dimension1 };
157  ErrorCode rval;
158 
159  Range results;
160  rval = mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, tags, vals, 2, results );
161 
162  if( MB_SUCCESS != rval )
163  return 0;
164  else
165  return results.front();
166 }
167 
169  EntityHandle not_this,
170  EntityHandle across,
171  EntityHandle& other )
172 {
173  other = 0;
174 
175  // get all children of bounded
176  Range bdy, tmpr;
177  ErrorCode rval = mdbImpl->get_child_meshsets( bounded, bdy );MB_CHK_SET_ERR( rval, "Failed to get the bounded entity's child meshsets" );
178 
179  // get all the parents of across
180  rval = mdbImpl->get_parent_meshsets( across, tmpr );
181 
182  // possible candidates is the intersection
183  bdy = intersect( bdy, tmpr );
184 
185  // if only two, choose the other
186  if( 1 == bdy.size() && *bdy.begin() == not_this )
187  {
188  return MB_SUCCESS;
189  }
190  else if( 2 == bdy.size() )
191  {
192  if( *bdy.begin() == not_this ) other = *bdy.rbegin();
193  if( *bdy.rbegin() == not_this )
194  other = *bdy.begin();
195  else
196  return MB_FAILURE;
197  }
198  else
199  {
200  return MB_FAILURE;
201  }
202 
203  return MB_SUCCESS;
204 }
205 
207 {
208 
210 
211  ErrorCode rval;
212  EntityHandle root;
213 
214  for( int dim = 2; dim <= 3; dim++ )
215  for( Range::iterator rit = geomRanges[dim].begin(); rit != geomRanges[dim].end(); ++rit )
216  {
217  rval = mdbImpl->tag_get_data( obbRootTag, &( *rit ), 1, &root );
218 
219  if( MB_SUCCESS == rval )
220  set_root_set( *rit, root );
221  else
222  {
223  return MB_TAG_NOT_FOUND;
224  }
225  }
226 
227  return MB_SUCCESS;
228 }
229 
231 {
232  ErrorCode rval;
233  // get all sets with this tag
234  Range geom_sets;
235 
236  if( 0 == geomTag )
237  {
238  rval = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag );MB_CHK_SET_ERR( rval, "Failed to get geom dimension tag handle" );
239  }
240 
241  rval = mdbImpl->get_entities_by_type_and_tag( modelSet, MBENTITYSET, &geomTag, NULL, 1, geom_sets );MB_CHK_SET_ERR( rval, "Failed to get the geometry entities" );
242 
243  rval = separate_by_dimension( geom_sets );MB_CHK_SET_ERR( rval, "Failed to separate geometry sets by dimension" );
244 
245  if( ranges )
246  {
247  for( int i = 0; i < 5; i++ )
248  {
249  ranges[i] = geomRanges[i];
250  }
251  }
252 
253  return MB_SUCCESS;
254 }
255 
257 {
258  ErrorCode rval;
259 
260  const int val = dim;
261  const void* const dim_val[] = { &val };
262  rval = mdbImpl->get_entities_by_type_and_tag( modelSet, MBENTITYSET, &geomTag, dim_val, 1, gset );MB_CHK_SET_ERR( rval, "Failed to get entity set by type and tag" );
263 
264  return MB_SUCCESS;
265 }
266 
268 {
269 
270  ErrorCode rval;
271 
272  // store original offset for later
273  EntityHandle orig_offset = setOffset;
274 
275  // get all surfaces and volumes
276  Range surfs, vols;
277  rval = get_gsets_by_dimension( 2, surfs );MB_CHK_SET_ERR( rval, "Could not get surface sets" );
278  rval = get_gsets_by_dimension( 3, vols );MB_CHK_SET_ERR( rval, "Could not get volume sets" );
279 
280  // check the vector size
281  Range surfs_and_vols;
282  surfs_and_vols = vols;
283  surfs_and_vols.merge( surfs );
284 
285  // update the setOffset
286  setOffset = surfs_and_vols.front();
287 
288  EntityHandle exp_size = surfs_and_vols.back() - setOffset + 1;
289 
290  // if new EnitytHandle(s) are lower than the original offset
291  if( setOffset < orig_offset )
292  {
293  // insert empty values at the beginning of the vector
294  rootSets.insert( rootSets.begin(), orig_offset - setOffset, 0 );
295  }
296 
297  if( exp_size != rootSets.size() )
298  {
299  // resize rootSets vector if necessary (new space will be added at the back)
300  rootSets.resize( exp_size );
301  }
302 
303  return MB_SUCCESS;
304 }
305 
307 {
308  // make sure entity set is part of the model
309  Range model_ents;
310  ErrorCode rval = mdbImpl->get_entities_by_handle( modelSet, model_ents );MB_CHK_SET_ERR( rval, "Failed to get entities" );
311  if( model_ents.find( eh ) == model_ents.end() )
312  {
313  MB_SET_ERR( MB_FAILURE, "Entity handle not in model set" );
314  }
315  return MB_SUCCESS;
316 }
317 
319 {
320 
321  ErrorCode rval;
322 
323  // Make sure this set is part of the model
324  rval = is_owned_set( gset );MB_CHK_SET_ERR( rval, "Entity set is not part of this model" );
325 
326  // Find the dimension of the entity
327  int dim;
328  rval = mdbImpl->tag_get_data( geomTag, &gset, 1, &dim );MB_CHK_SET_ERR( rval, "Failed to get dimension" );
329 
330  // Attempt to find a root for this set
331  EntityHandle root;
332  rval = get_root( gset, root );MB_CHK_SET_ERR( rval, "Failed to find an obb tree root for the entity set" );
333 
334  // Create range of tree nodes to delete
335  Range nodes_to_delete;
336  nodes_to_delete.insert( root );
337 
338  // If passed ent is a vol and 'vol_only' is true, delete vol root and all nodes between vol and
339  // surf root
340  if( dim == 3 && vol_only )
341  {
342  // Range of child nodes to check before adding to delete list
343  Range child_tree_nodes;
344  rval = mdbImpl->get_child_meshsets( root, child_tree_nodes );MB_CHK_SET_ERR( rval, "Problem getting child tree nodes" );
345 
346  // Traverse the tree, checking each child node until
347  // a surface root node is reached
348  while( child_tree_nodes.size() != 0 )
349  {
350  EntityHandle child = *child_tree_nodes.begin();
351  EntityHandle surf;
352  rval = mdbImpl->tag_get_data( obbGsetTag, &child, 1, &surf );
353  // If the node has a gset tag, it is a surf root. Stop here.
354  // If not, it is a tree node that needs to 1) have its children checked and
355  // 2) be added to delete range
356  if( MB_TAG_NOT_FOUND == rval )
357  {
358  Range new_child_tree_nodes;
359  rval = mdbImpl->get_child_meshsets( child, new_child_tree_nodes );MB_CHK_SET_ERR( rval, "Problem getting child nodes" );
360  child_tree_nodes.insert_list( new_child_tree_nodes.begin(), new_child_tree_nodes.end() );
361  nodes_to_delete.insert( child );
362  }
363  // We're done checking this node, so can erase from child list
364  child_tree_nodes.erase( child );
365  }
366  }
367  // If passed ent is a surf or a vol and 'vol_only' is false, recursively gather all child nodes
368  // and add them to delete list
369  else
370  {
371  Range all_tree_nodes;
372  rval = mdbImpl->get_child_meshsets( root, all_tree_nodes, 0 );MB_CHK_SET_ERR( rval, "Failed to get child tree node sets" );
373  nodes_to_delete.insert_list( all_tree_nodes.begin(), all_tree_nodes.end() );
374  }
375 
376  // Remove the root nodes from the GTT data structures
377  for( Range::iterator it = nodes_to_delete.begin(); it != nodes_to_delete.end(); ++it )
378  {
379  // Check to see if node is a root
380  EntityHandle vol_or_surf;
381  rval = mdbImpl->tag_get_data( obbGsetTag, &( *it ), 1, &vol_or_surf );
382  if( MB_SUCCESS == rval )
383  {
384  // Remove from set of all roots
385  rval = remove_root( vol_or_surf );MB_CHK_SET_ERR( rval, "Failed to remove node from GTT data structure" );
386  }
387  }
388 
389  // Delete the tree nodes from the database
390  rval = mdbImpl->delete_entities( nodes_to_delete );MB_CHK_SET_ERR( rval, "Failed to delete node set" );
391 
392  return MB_SUCCESS;
393 }
394 
396 {
397 
398  ErrorCode rval;
399 
400  for( Range::iterator rit = geomRanges[3].begin(); rit != geomRanges[3].end(); ++rit )
401  {
402  EntityHandle root;
403  rval = mdbImpl->tag_get_data( obbRootTag, &( *rit ), 1, &root );
404  if( MB_SUCCESS == rval )
405  {
406  rval = delete_obb_tree( *rit, false );MB_CHK_SET_ERR( rval, "Failed to delete obb tree" );
407  }
408  }
409 
410  return MB_SUCCESS;
411 }
412 
414 {
415  ErrorCode rval;
416  int dim;
417 
418  rval = is_owned_set( eh );MB_CHK_SET_ERR( rval, "Entity set is not part of this model" );
419 
420  // get the type
421  EntityType type = mdbImpl->type_from_handle( eh );
422 
423  // find the dimension of the entity
424  rval = mdbImpl->tag_get_data( geomTag, &eh, 1, &dim );MB_CHK_SET_ERR( rval, "Failed to get dimension" );
425 
426  // ensure that the rootSets vector is of the correct size
427  if( m_rootSets_vector && ( eh < setOffset || eh >= setOffset + rootSets.size() ) )
428  {
429  rval = resize_rootSets();MB_CHK_SET_ERR( rval, "Error setting offset and sizing rootSets vector." );
430  }
431 
432  EntityHandle root;
433  // if it's a surface
434  if( dim == 2 && type == MBENTITYSET )
435  {
436  rval = get_root( eh, root );
437  if( MB_SUCCESS == rval )
438  {
439  std::cerr << "Surface obb tree already exists" << std::endl;
440  return MB_SUCCESS;
441  }
442  else if( MB_INDEX_OUT_OF_RANGE != rval )
443  {
444  MB_CHK_SET_ERR( rval, "Failed to get surface obb tree root" );
445  }
446 
447  Range tris;
448  rval = mdbImpl->get_entities_by_dimension( eh, 2, tris );MB_CHK_SET_ERR( rval, "Failed to get entities by dimension" );
449 
450  if( tris.empty() )
451  {
452  std::cerr << "WARNING: Surface id " << global_id(eh) << " (handle: " << eh << ")" << "has no facets" << std::endl;
453  }
454 
455  rval = obbTree->build( tris, root );MB_CHK_SET_ERR( rval, "Failed to build obb Tree for surface" );
456 
457  rval = mdbImpl->add_entities( root, &eh, 1 );MB_CHK_SET_ERR( rval, "Failed to add entities to root set" );
458 
459  // add this root to the GeomTopoTool tree root indexing
460  set_root_set( eh, root );
461 
462  // if just building tree for surface, return here
463  return MB_SUCCESS;
464  }
465  // if it's a volume
466  else if( dim == 3 && type == MBENTITYSET )
467  {
468  // get its surfaces
469  Range tmp_surfs, surf_trees;
470  rval = mdbImpl->get_child_meshsets( eh, tmp_surfs );MB_CHK_SET_ERR( rval, "Failed to get surface meshsets" );
471 
472  // get OBB trees or create for each surface
473  for( Range::iterator j = tmp_surfs.begin(); j != tmp_surfs.end(); ++j )
474  {
475  rval = get_root( *j, root );
476  // if root doesn't exist, create obb tree
477  if( MB_INDEX_OUT_OF_RANGE == rval )
478  {
479  rval = construct_obb_tree( *j );MB_CHK_SET_ERR( rval, "Failed to get create surface obb tree" );
480  rval = get_root( *j, root );MB_CHK_SET_ERR( rval, "Failed to get surface obb tree root" );
481  }
482  else
483  {
484  MB_CHK_SET_ERR( rval, "Failed to get surface obb tree root" );
485  }
486 
487  surf_trees.insert( root );
488  }
489 
490  // build OBB tree for volume
491  rval = obbTree->join_trees( surf_trees, root );MB_CHK_SET_ERR( rval, "Failed to join the obb trees" );
492 
493  // add this root to the GeomTopoTool tree root indexing
494  set_root_set( eh, root );
495 
496  return MB_SUCCESS;
497  }
498  else
499  {
500  MB_SET_ERR( MB_FAILURE, "Improper dimension or type for constructing obb tree" );
501  }
502 }
503 
505 {
506 
507  // Tag the vol or surf with its obb root (obbRootTag)
508  ErrorCode rval;
509  rval = mdbImpl->tag_set_data( obbRootTag, &vol_or_surf, 1, &root );MB_CHK_SET_ERR( rval, "Failed to set the obb root tag" );
510 
511  // Tag obb root with corresponding gset (obbGsetTag)
512  rval = mdbImpl->tag_set_data( obbGsetTag, &root, 1, &vol_or_surf );MB_CHK_SET_ERR( rval, "Failed to set the obb gset tag" );
513 
514  // Add to the set of all roots
515  if( m_rootSets_vector )
516  rootSets[vol_or_surf - setOffset] = root;
517  else
518  mapRootSets[vol_or_surf] = root;
519 
520  return MB_SUCCESS;
521 }
522 
524 {
525 
526  // Find the root of the vol or surf
527  ErrorCode rval;
528  EntityHandle root;
529  rval = mdbImpl->tag_get_data( obbRootTag, &( vol_or_surf ), 1, &root );MB_CHK_SET_ERR( rval, "Failed to get obb root tag" );
530 
531  // If the ent is a vol, remove its root from obbtreetool
532  int dim;
533  rval = mdbImpl->tag_get_data( geomTag, &vol_or_surf, 1, &dim );MB_CHK_SET_ERR( rval, "Failed to get dimension" );
534  if( dim == 3 )
535  {
536  rval = obbTree->remove_root( root );MB_CHK_SET_ERR( rval, "Failed to remove root from obbTreeTool" );
537  }
538 
539  // Delete the obbGsetTag data from the root
540  rval = mdbImpl->tag_delete_data( obbGsetTag, &root, 1 );MB_CHK_SET_ERR( rval, "Failed to delete obb root tag" );
541 
542  // Delete the obbRootTag data from the vol or surf
543  rval = mdbImpl->tag_delete_data( obbRootTag, &vol_or_surf, 1 );MB_CHK_SET_ERR( rval, "Failed to delete obb root tag" );
544 
545  // Remove the root from set of all roots
546  if( m_rootSets_vector )
547  {
548  unsigned int index = vol_or_surf - setOffset;
549  if( index < rootSets.size() )
550  {
551  rootSets[index] = 0;
552  }
553  else
554  {
555  return MB_INDEX_OUT_OF_RANGE;
556  }
557  }
558  else
559  {
560  mapRootSets[vol_or_surf] = 0;
561  }
562 
563  return MB_SUCCESS;
564 }
565 
567 {
568  ErrorCode rval;
569  EntityHandle root;
570 
571  // get all surfaces and volumes
572  Range surfs, vols, vol_trees;
573  rval = get_gsets_by_dimension( 2, surfs );MB_CHK_SET_ERR( rval, "Could not get surface sets" );
574  rval = get_gsets_by_dimension( 3, vols );MB_CHK_SET_ERR( rval, "Could not get volume sets" );
575 
576  // for surface
577  Range one_vol_trees;
578  for( Range::iterator i = surfs.begin(); i != surfs.end(); ++i )
579  {
580  rval = construct_obb_tree( *i );MB_CHK_SET_ERR( rval, "Failed to construct obb tree for surface" );
581  // get the root set of this volume
582  rval = get_root( *i, root );MB_CHK_SET_ERR( rval, "Failed to get obb tree root for surface" );
583  // add to the Range of volume root sets
584  one_vol_trees.insert( root );
585  }
586 
587  // for volumes
588  for( Range::iterator i = vols.begin(); i != vols.end(); ++i )
589  {
590  // create tree for this volume
591  rval = construct_obb_tree( *i );MB_CHK_SET_ERR( rval, "Failed to construct obb tree for volume" );
592  }
593 
594  // build OBB tree for volume
595  if( make_one_vol )
596  {
597  rval = obbTree->join_trees( one_vol_trees, root );MB_CHK_SET_ERR( rval, "Failed to join surface trees into one volume" );
598  oneVolRootSet = root;
599  }
600 
601  return rval;
602 }
603 
604 //! Restore parent/child links between GEOM_TOPO mesh sets
606 {
607 
608  // look for geometric topology sets and restore parent/child links between them
609  // algorithm:
610  // - for each entity of dimension d=D-1..0:
611  // . get d-dimensional entity in entity
612  // . get all (d+1)-dim adjs to that entity
613  // . for each geom entity if dim d+1, if it contains any of the ents,
614  // add it to list of parents
615  // . make parent/child links with parents
616 
617  // the geomRanges are already known, separated by dimension
618 
619  std::vector< EntityHandle > parents;
620  Range tmp_parents;
621  ErrorCode result;
622 
623  // loop over dimensions
624  for( int dim = 2; dim >= 0; dim-- )
625  {
626  // mark entities of next higher dimension with their owners; regenerate tag
627  // each dimension so prev dim's tag data goes away
628  Tag owner_tag;
629  EntityHandle dum_val = 0;
630  result = mdbImpl->tag_get_handle( "__owner_tag", 1, MB_TYPE_HANDLE, owner_tag, MB_TAG_DENSE | MB_TAG_CREAT,
631  &dum_val );
632  if( MB_SUCCESS != result ) continue;
633  Range dp1ents;
634  std::vector< EntityHandle > owners;
635  for( Range::iterator rit = geomRanges[dim + 1].begin(); rit != geomRanges[dim + 1].end(); ++rit )
636  {
637  dp1ents.clear();
638  result = mdbImpl->get_entities_by_dimension( *rit, dim + 1, dp1ents );
639  if( MB_SUCCESS != result ) continue;
640  owners.resize( dp1ents.size() );
641  std::fill( owners.begin(), owners.end(), *rit );
642  result = mdbImpl->tag_set_data( owner_tag, dp1ents, &owners[0] );
643  if( MB_SUCCESS != result ) continue;
644  }
645 
646  for( Range::iterator d_it = geomRanges[dim].begin(); d_it != geomRanges[dim].end(); ++d_it )
647  {
648  Range dents;
649  result = mdbImpl->get_entities_by_dimension( *d_it, dim, dents );
650  if( MB_SUCCESS != result ) continue;
651  if( dents.empty() ) continue;
652 
653  // get (d+1)-dimensional adjs
654  dp1ents.clear();
655  result = mdbImpl->get_adjacencies( &( *dents.begin() ), 1, dim + 1, false, dp1ents );
656  if( MB_SUCCESS != result || dp1ents.empty() ) continue;
657 
658  // get owner tags
659  parents.resize( dp1ents.size() );
660  result = mdbImpl->tag_get_data( owner_tag, dp1ents, &parents[0] );
661  if( MB_TAG_NOT_FOUND == result )
662  {
663  MB_CHK_SET_ERR( result, "Could not find owner tag" );
664  }
665  if( MB_SUCCESS != result ) continue;
666 
667  // compress to a range to remove duplicates
668  tmp_parents.clear();
669  std::copy( parents.begin(), parents.end(), range_inserter( tmp_parents ) );
670  for( Range::iterator pit = tmp_parents.begin(); pit != tmp_parents.end(); ++pit )
671  {
672  result = mdbImpl->add_parent_child( *pit, *d_it );MB_CHK_SET_ERR( result, "Failed to create parent child relationship" );
673  }
674 
675  // store surface senses within regions, and edge senses within surfaces
676  if( dim == 0 ) continue;
677  const EntityHandle *conn3 = NULL, *conn2 = NULL;
678  int len3 = 0, len2 = 0, err = 0, num = 0, sense = 0, offset = 0;
679  for( size_t i = 0; i < parents.size(); ++i )
680  {
681  result = mdbImpl->get_connectivity( dp1ents[i], conn3, len3, true );MB_CHK_SET_ERR( result, "Failed to get the connectivity of the element" );
682  result = mdbImpl->get_connectivity( dents.front(), conn2, len2, true );MB_CHK_SET_ERR( result, "Failed to get the connectivity of the first element" );
683  if( len2 > 4 )
684  {
685  MB_SET_ERR( MB_FAILURE, "Connectivity of incorrect length" );
686  }
687  err = CN::SideNumber( TYPE_FROM_HANDLE( dp1ents[i] ), conn3, conn2, len2, dim, num, sense, offset );
688  if( err ) return MB_FAILURE;
689 
690  result = set_sense( *d_it, parents[i], sense );
691  if( MB_MULTIPLE_ENTITIES_FOUND == result )
692  {
693  if( 2 == dim )
694  std::cerr << "Warning: Multiple volumes use surface with same sense." << std::endl
695  << " Some geometric sense data lost." << std::endl;
696  }
697  else if( MB_SUCCESS != result )
698  {
699  return result;
700  }
701  }
702  }
703 
704  // now delete owner tag on this dimension, automatically removes tag data
705  result = mdbImpl->tag_delete( owner_tag );MB_CHK_SET_ERR( result, "Failed to delete the owner tag" );
706 
707  } // dim
708 
709  return result;
710 }
711 
713 {
714  ErrorCode result;
715 
716  result = check_geom_tag();MB_CHK_SET_ERR( result, "Could not verify geometry dimension tag" );
717 
718  // get the data for those tags
719  std::vector< int > tag_vals( geom_sets.size() );
720  result = mdbImpl->tag_get_data( geomTag, geom_sets, &tag_vals[0] );MB_CHK_SET_ERR( result, "Failed to get the geometry dimension tag" );
721 
723  std::vector< int >::iterator iit;
724 
725  for( int i = 0; i < 5; i++ )
726  this->geomRanges[i].clear();
727 
728  for( git = geom_sets.begin(), iit = tag_vals.begin(); git != geom_sets.end(); ++git, ++iit )
729  {
730  if( 0 <= *iit && 4 >= *iit )
731  geomRanges[*iit].insert( *git );
732  else
733  {
734  // assert(false);
735  // do nothing for now
736  }
737  }
738 
739  // establish the max global ids so far, per dimension
740  if( 0 == gidTag )
741  {
743  }
744 
745  for( int i = 0; i <= 4; i++ )
746  {
747  maxGlobalId[i] = 0;
748  for( Range::iterator it = geomRanges[i].begin(); it != geomRanges[i].end(); ++it )
749  {
750  EntityHandle set = *it;
751  int gid;
752 
753  result = mdbImpl->tag_get_data( gidTag, &set, 1, &gid );
754  if( MB_SUCCESS == result )
755  {
756  if( gid > maxGlobalId[i] ) maxGlobalId[i] = gid;
757  }
758  }
759  }
760 
761  return MB_SUCCESS;
762 }
763 
764 ErrorCode GeomTopoTool::construct_vertex_ranges( const Range& geom_sets, const Tag verts_tag )
765 {
766  // construct the vertex range for each entity and put on that tag
767  Range *temp_verts, temp_elems;
768  ErrorCode result = MB_SUCCESS;
769  for( Range::const_iterator it = geom_sets.begin(); it != geom_sets.end(); ++it )
770  {
771  temp_elems.clear();
772 
773  // get all the elements in the set, recursively
774  result = mdbImpl->get_entities_by_handle( *it, temp_elems, true );MB_CHK_SET_ERR( result, "Failed to get the geometry set entities" );
775 
776  // make the new range
777  temp_verts = new( std::nothrow ) Range();
778  if( NULL == temp_verts )
779  {
780  MB_SET_ERR( MB_FAILURE, "Could not construct Range object" );
781  }
782 
783  // get all the verts of those elements; use get_adjacencies 'cuz it handles ranges better
784  result = mdbImpl->get_adjacencies( temp_elems, 0, false, *temp_verts, Interface::UNION );
785  if( MB_SUCCESS != result )
786  {
787  delete temp_verts;
788  }
789  MB_CHK_SET_ERR( result, "Failed to get the element's adjacent vertices" );
790 
791  // store this range as a tag on the entity
792  result = mdbImpl->tag_set_data( verts_tag, &( *it ), 1, &temp_verts );
793  if( MB_SUCCESS != result )
794  {
795  delete temp_verts;
796  }
797  MB_CHK_SET_ERR( result, "Failed to get the adjacent vertex data" );
798 
799  delete temp_verts;
800  temp_verts = NULL;
801  }
802 
803  return result;
804 }
805 
806 //! Store sense of entity relative to wrt_entity.
807 //!\return MB_MULTIPLE_ENTITIES_FOUND if surface already has a forward volume.
808 //! MB_SUCCESS if successful
809 //! otherwise whatever internal error code occurred.
811 {
812  // entity is lower dim (edge or face), wrt_entity is face or volume
813  int edim = dimension( entity );
814  int wrtdim = dimension( wrt_entity );
815  if( -1 == edim || -1 == wrtdim ) MB_SET_ERR( MB_FAILURE, "Non-geometric entity provided" );
816  if( wrtdim - edim != 1 ) MB_SET_ERR( MB_FAILURE, "Entity dimension mismatch" );
817  if( sense < -1 || sense > 1 ) MB_SET_ERR( MB_FAILURE, "Invalid sense data provided" );
818 
819  ErrorCode rval;
820 
821  if( 1 == edim )
822  {
823  // this case is about setting the sense of an edge in a face
824  // it could be -1, 0 (rare, non manifold), or 1
825  rval = check_edge_sense_tags( true );MB_CHK_SET_ERR( rval, "Failed to check the curve to surface sense tag handles" );
826  std::vector< EntityHandle > higher_ents;
827  std::vector< int > senses;
828  rval = get_senses( entity, higher_ents, senses ); // the tags should be defined here
829  // if there are no higher_ents, we are fine, we will just set them
830  // if wrt_entity is not among higher_ents, we will add it to the list
831  // it is possible the entity (edge set) has no prior faces adjancent; in that case, the
832  // tag would not be set, and rval could be MB_TAG_NOT_FOUND; it is not a fatal error
833  if( MB_SUCCESS != rval && MB_TAG_NOT_FOUND != rval )
834  MB_CHK_SET_ERR( rval, "cannot determine sense tags for edge" );
835  bool append = true;
836  if( !higher_ents.empty() )
837  {
838  std::vector< EntityHandle >::iterator it = std::find( higher_ents.begin(), higher_ents.end(), wrt_entity );
839  if( it != higher_ents.end() )
840  {
841  // we should not reset the sense, if the sense is the same
842  // if the sense is different, put BOTH
843  unsigned int idx = it - higher_ents.begin();
844  int oldSense = senses[idx];
845  if( oldSense == sense ) return MB_SUCCESS; // sense already set fine, do not reset
846  if( 0 != oldSense && oldSense + sense != 0 ) return MB_MULTIPLE_ENTITIES_FOUND;
847  senses[idx] = SENSE_BOTH; // allow double senses
848  // do not need to add a new sense, but still need to reset the tag
849  // because of a new value
850  append = false;
851  }
852  }
853  if( append )
854  {
855  // what happens if a var tag data was already set before, and now it is
856  // reset with a different size??
857  higher_ents.push_back( wrt_entity );
858  senses.push_back( sense );
859  }
860  // finally, set the senses :
861  int dum_size = higher_ents.size();
862  void* dum_ptr = &higher_ents[0];
863  rval = mdbImpl->tag_set_by_ptr( senseNEntsTag, &entity, 1, &dum_ptr, &dum_size );MB_CHK_SET_ERR( rval, "Failed to set the sense data" );
864 
865  dum_ptr = &senses[0];
866  dum_size = higher_ents.size();
867  rval = mdbImpl->tag_set_by_ptr( senseNSensesTag, &entity, 1, &dum_ptr, &dum_size );MB_CHK_SET_ERR( rval, "Failed to set the sense data by pointer" );
868  }
869  else
870  {
871  // this case is about a face in the volume
872  // there could be only 2 volumes
873 
874  rval = check_face_sense_tag( true );MB_CHK_SET_ERR( rval, "Failed to verify the face sense tag" );
875 
876  EntityHandle sense_data[2] = { 0, 0 };
877  rval = mdbImpl->tag_get_data( sense2Tag, &entity, 1, sense_data );
878  if( MB_TAG_NOT_FOUND != rval && MB_SUCCESS != rval ) MB_CHK_SET_ERR( rval, "Failed to get the sense2Tag data" );
879 
880  if( 0 == sense )
881  {
882  if( 0 != sense_data[0] && wrt_entity != sense_data[0] ) return MB_MULTIPLE_ENTITIES_FOUND;
883  if( 0 != sense_data[1] && wrt_entity != sense_data[1] ) return MB_MULTIPLE_ENTITIES_FOUND;
884  sense_data[0] = sense_data[1] = wrt_entity;
885  }
886  else if( -1 == sense )
887  {
888  if( 0 != sense_data[1] && wrt_entity != sense_data[1] ) return MB_MULTIPLE_ENTITIES_FOUND;
889  if( sense_data[1] == wrt_entity ) return MB_SUCCESS; // already set as we want
890  sense_data[1] = wrt_entity;
891  }
892  else if( 1 == sense )
893  {
894  if( 0 != sense_data[0] && wrt_entity != sense_data[0] ) return MB_MULTIPLE_ENTITIES_FOUND;
895  if( sense_data[0] == wrt_entity ) return MB_SUCCESS; // already set as we want
896  sense_data[0] = wrt_entity;
897  }
898  return mdbImpl->tag_set_data( sense2Tag, &entity, 1, sense_data );
899  }
900  return MB_SUCCESS;
901 }
902 
903 //! Get the sense of entity with respect to wrt_entity
904 //! Returns MB_ENTITY_NOT_FOUND if no relationship found
906 {
907  // entity is lower dim (edge or face), wrt_entity is face or volume
908  int edim = dimension( entity );
909  int wrtdim = dimension( wrt_entity );
910  if( -1 == edim || -1 == wrtdim ) MB_SET_ERR( MB_FAILURE, "Non-geometric entity provided" );
911  if( wrtdim - edim != 1 ) MB_SET_ERR( MB_FAILURE, "Entity dimension mismatch" );
912 
913  ErrorCode rval;
914 
915  if( 1 == edim )
916  {
917  // edge in face
918  rval = check_edge_sense_tags( false );MB_CHK_SET_ERR( rval, "Failed to check the curve to surface sense tag handles" );
919 
920  std::vector< EntityHandle > faces;
921  std::vector< int > senses;
922  rval = get_senses( entity, faces, senses ); // the tags should be defined here
923  MB_CHK_SET_ERR( rval, "Failed to get the curve to surface sense data" );
924 
925  std::vector< EntityHandle >::iterator it = std::find( faces.begin(), faces.end(), wrt_entity );
926  if( it == faces.end() ) return MB_ENTITY_NOT_FOUND;
927  unsigned int index = it - faces.begin();
928  sense = senses[index];
929  }
930  else
931  {
932  // face in volume
933  rval = check_face_sense_tag( false );MB_CHK_SET_ERR( rval, "Failed to check the surface to volume sense tag handle" );
934  EntityHandle sense_data[2] = { 0, 0 };
935  rval = mdbImpl->tag_get_data( sense2Tag, &entity, 1, sense_data );
936  if( MB_TAG_NOT_FOUND != rval && MB_SUCCESS != rval )
937  MB_CHK_SET_ERR( rval, "Failed to get the surface to volume sense data" );
938  if( ( wrt_entity == sense_data[0] ) && ( wrt_entity == sense_data[1] ) )
939  sense = 0;
940  else if( wrt_entity == sense_data[0] )
941  sense = 1;
942  else if( wrt_entity == sense_data[1] )
943  sense = -1;
944  else
945  return MB_ENTITY_NOT_FOUND;
946  }
947  return MB_SUCCESS;
948 }
949 
951  EntityHandle& forward_vol,
952  EntityHandle& reverse_vol )
953 {
954  ErrorCode rval;
955  // this method should only be called to retrieve surface to volume
956  // sense relationships
957  int ent_dim = dimension( surface_ent );
958  // verify the incoming entity dimensions for this call
959  if( ent_dim != 2 )
960  {
961  MB_SET_ERR( MB_FAILURE, "Entity dimension is incorrect for surface meshset" );
962  }
963 
964  // get the sense information for this surface
965  EntityHandle parent_vols[2] = { 0, 0 };
966  rval = mdbImpl->tag_get_data( sense2Tag, &surface_ent, 1, parent_vols );MB_CHK_SET_ERR( rval, "Failed to get surface sense data" );
967 
968  // set the outgoing values
969  forward_vol = parent_vols[0];
970  reverse_vol = parent_vols[1];
971 
972  return MB_SUCCESS;
973 }
974 
976  EntityHandle forward_vol,
977  EntityHandle reverse_vol )
978 {
979  ErrorCode rval;
980  // this mthod should only be called to retrieve surface to volume
981  // sense relationships
982  int ent_dim = dimension( surface_ent );
983  // verify the incoming entity dimensions for this call
984  if( ent_dim != 2 )
985  {
986  MB_SET_ERR( MB_FAILURE, "Entity dimension is incorrect for surface meshset" );
987  }
988 
989  // set the sense information for this surface
990  EntityHandle parent_vols[2] = { forward_vol, reverse_vol };
991  rval = mdbImpl->tag_set_data( sense2Tag, &surface_ent, 1, parent_vols );MB_CHK_SET_ERR( rval, "Failed to set surface sense data" );
992 
993  return MB_SUCCESS;
994 }
995 
996 // get sense of surface(s) wrt volume
998  int num_surfaces,
999  const EntityHandle* surfaces,
1000  int* senses_out )
1001 {
1002 
1003  /* The sense tags do not reference the implicit complement handle.
1004  All surfaces that interact with the implicit complement should have
1005  a null handle in the direction of the implicit complement. */
1006  // if (volume == impl_compl_handle)
1007  // volume = (EntityHandle) 0;
1008 
1009  for( int surf_num = 0; surf_num < num_surfaces; surf_num++ )
1010  {
1011  get_sense( surfaces[surf_num], volume, senses_out[surf_num] );
1012  }
1013 
1014  return MB_SUCCESS;
1015 }
1016 
1018  std::vector< EntityHandle >& wrt_entities,
1019  std::vector< int >& senses )
1020 {
1021  // the question here is: the wrt_entities is supplied or not?
1022  // I assume not, we will obtain it !!
1023  int edim = dimension( entity );
1024 
1025  if( -1 == edim ) MB_SET_ERR( MB_FAILURE, "Non-geometric entity provided" );
1026 
1027  ErrorCode rval;
1028  wrt_entities.clear();
1029  senses.clear();
1030 
1031  if( 1 == edim ) // edge
1032  {
1033  rval = check_edge_sense_tags( false );MB_CHK_SET_ERR( rval, "Failed to check the curve to surface sense tag handles" );
1034  const void* dum_ptr;
1035  int num_ents;
1036  rval = mdbImpl->tag_get_by_ptr( senseNEntsTag, &entity, 1, &dum_ptr, &num_ents );MB_CHK_ERR( rval );
1037 
1038  const EntityHandle* ents_data = static_cast< const EntityHandle* >( dum_ptr );
1039  std::copy( ents_data, ents_data + num_ents, std::back_inserter( wrt_entities ) );
1040 
1041  rval = mdbImpl->tag_get_by_ptr( senseNSensesTag, &entity, 1, &dum_ptr, &num_ents );MB_CHK_ERR( rval );
1042 
1043  const int* senses_data = static_cast< const int* >( dum_ptr );
1044  std::copy( senses_data, senses_data + num_ents, std::back_inserter( senses ) );
1045  }
1046  else // face in volume, edim == 2
1047  {
1048  rval = check_face_sense_tag( false );MB_CHK_SET_ERR( rval, "Failed to check the surface to volume sense tag handle" );
1049  EntityHandle sense_data[2] = { 0, 0 };
1050  rval = mdbImpl->tag_get_data( sense2Tag, &entity, 1, sense_data );MB_CHK_SET_ERR( rval, "Failed to get the surface to volume sense data" );
1051  if( sense_data[0] != 0 && sense_data[1] == sense_data[0] )
1052  {
1053  wrt_entities.push_back( sense_data[0] );
1054  senses.push_back( 0 ); // both
1055  }
1056  else
1057  {
1058  if( sense_data[0] != 0 )
1059  {
1060  wrt_entities.push_back( sense_data[0] );
1061  senses.push_back( 1 );
1062  }
1063  if( sense_data[1] != 0 )
1064  {
1065  wrt_entities.push_back( sense_data[1] );
1066  senses.push_back( -1 );
1067  }
1068  }
1069  }
1070  // filter the results with the sets that are in the model at this time
1071  // this was introduced because extracting some sets (e.g. neumann set, with mbconvert)
1072  // from a model would leave some sense tags not defined correctly
1073  // also, the geom ent set really needs to be part of the current model set
1074  unsigned int currentSize = 0;
1075 
1076  for( unsigned int index = 0; index < wrt_entities.size(); index++ )
1077  {
1078  EntityHandle wrt_ent = wrt_entities[index];
1079  if( wrt_ent )
1080  {
1081  if( mdbImpl->contains_entities( modelSet, &wrt_ent, 1 ) )
1082  {
1083  wrt_entities[currentSize] = wrt_entities[index];
1084  senses[currentSize] = senses[index];
1085  currentSize++;
1086  }
1087  }
1088  }
1089  wrt_entities.resize( currentSize );
1090  senses.resize( currentSize );
1091  //
1092  return MB_SUCCESS;
1093 }
1094 
1096  std::vector< EntityHandle >& wrt_entities,
1097  std::vector< int >& senses )
1098 {
1099  // not efficient, and maybe wrong
1100  for( unsigned int i = 0; i < wrt_entities.size(); i++ )
1101  {
1102  ErrorCode rval = set_sense( entity, wrt_entities[i], senses[i] );MB_CHK_SET_ERR( rval, "Failed to set the sense" );
1103  }
1104 
1105  return MB_SUCCESS;
1106 }
1107 
1109 {
1110  std::vector< EntityHandle > parents;
1111  ErrorCode rval = mdbImpl->get_parent_meshsets( surface, parents );
1112 
1113  if( MB_SUCCESS == rval )
1114  {
1115  if( parents.size() != 2 )
1116  rval = MB_FAILURE;
1117  else if( parents.front() == old_volume )
1118  new_volume = parents.back();
1119  else if( parents.back() == old_volume )
1120  new_volume = parents.front();
1121  else
1122  rval = MB_FAILURE;
1123  }
1124 
1125  return rval;
1126 }
1127 
1129 {
1130  ErrorCode rval;
1131  unsigned flags = create ? MB_TAG_DENSE | MB_TAG_CREAT : MB_TAG_DENSE;
1132  if( !geomTag )
1133  {
1134  // get any kind of tag that already exists
1135  rval = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag, flags );MB_CHK_SET_ERR( rval, "Could not get/create the geometry dimension tag" );
1136  }
1137  return MB_SUCCESS;
1138 }
1139 
1141 {
1142  ErrorCode rval;
1143  unsigned flags = create ? MB_TAG_DENSE | MB_TAG_CREAT : MB_TAG_DENSE;
1144  if( !gidTag )
1145  {
1146  // get any kind of tag that already exists
1147  rval = mdbImpl->tag_get_handle( GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, gidTag, flags );MB_CHK_SET_ERR( rval, "Could not get/create the global id tag" );
1148  }
1149  return MB_SUCCESS;
1150 }
1151 
1152 // move the sense tag existence creation in private methods
1153 // verify sense face tag
1155 {
1156  ErrorCode rval;
1157  unsigned flags = create ? MB_TAG_SPARSE | MB_TAG_CREAT | MB_TAG_ANY : MB_TAG_SPARSE | MB_TAG_ANY;
1158  if( !sense2Tag )
1159  {
1160  EntityHandle def_val[2] = { 0, 0 };
1161  rval = mdbImpl->tag_get_handle( GEOM_SENSE_2_TAG_NAME, 2, MB_TYPE_HANDLE, sense2Tag, flags, def_val );MB_CHK_SET_ERR( rval, "Could not get/create the sense2Tag" );
1162  }
1163  return MB_SUCCESS;
1164 }
1165 
1166 // verify sense edge tags
1168 {
1169  ErrorCode rval;
1170  unsigned flags = MB_TAG_VARLEN | MB_TAG_SPARSE;
1171  if( create ) flags |= MB_TAG_CREAT;
1172  if( !senseNEntsTag )
1173  {
1174  rval = mdbImpl->tag_get_handle( GEOM_SENSE_N_ENTS_TAG_NAME, 0, MB_TYPE_HANDLE, senseNEntsTag, flags );MB_CHK_SET_ERR( rval, "Failed to get the curve to surface entity tag handle" );
1175  rval = mdbImpl->tag_get_handle( GEOM_SENSE_N_SENSES_TAG_NAME, 0, MB_TYPE_INTEGER, senseNSensesTag, flags );MB_CHK_SET_ERR( rval, "Failed to get the curve to surface sense tag handle" );
1176  }
1177  return MB_SUCCESS;
1178 }
1179 
1181 {
1182  if( dim < 0 || dim > 4 ) MB_SET_ERR( MB_FAILURE, "Invalid geometric dimension provided" );
1183 
1184  // see if it is not already set
1185  if( geomRanges[dim].find( set ) != geomRanges[dim].end() )
1186  {
1187  return MB_SUCCESS; // nothing to do, we already have it as a geo set of proper dimension
1188  }
1189  updated = false; // this should trigger at least an obb recomputation
1190  // get the geom topology tag
1191  ErrorCode result;
1192  if( 0 == geomTag )
1193  {
1194  result = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag );MB_CHK_SET_ERR( result, "Failed to get the geometry dimension tag handle" );
1195  }
1196 
1197  if( 0 == gidTag )
1198  {
1200  }
1201 
1202  // make sure the added set has the geom tag properly set
1203  result = mdbImpl->tag_set_data( geomTag, &set, 1, &dim );MB_CHK_SET_ERR( result, "Failed set the geometry dimension tag value" );
1204 
1205  geomRanges[dim].insert( set );
1206  // not only that, but also add it to the root model set
1207  if( modelSet )
1208  {
1209  result = mdbImpl->add_entities( modelSet, &set, 1 );MB_CHK_SET_ERR( result, "Failed to add new geometry set to the tool's modelSet" );
1210  }
1211 
1212  // set the global ID value
1213  // if passed 0, just increase the max id for the dimension
1214  if( 0 == gid )
1215  {
1216  gid = ++maxGlobalId[dim];
1217  }
1218 
1219  result = mdbImpl->tag_set_data( gidTag, &set, 1, &gid );MB_CHK_SET_ERR( result, "Failed to get the global id tag value for the geom entity" );
1220 
1221  return MB_SUCCESS;
1222 }
1223 
1224 // will assume no geo sets are defined for this surface
1225 // will output a mesh_set that contains everything (all sets of interest), for proper output
1227 {
1228  // usual scenario is to read a surface smf file, and "geometrize" it, and output it as a
1229  // h5m file with proper sets, tags defined for mesh-based geometry
1230 
1231  // get all triangles/quads from the surface, then build loops
1232  // we may even have to
1233  // proper care has to be given to the orientation, material to the left!!!
1234  // at some point we may have to reorient triangles, not only edges, for proper definition
1235  bool debugFlag = false;
1236 
1237  Range surface_ents, edge_ents, loop_range;
1238 
1239  // most of these should be triangles and quads
1240  ErrorCode rval = mdbImpl->get_entities_by_dimension( surface, 2, surface_ents );MB_CHK_SET_ERR( rval, "Failed to get the surface entities" );
1241 
1242  EntityHandle face = surface;
1243  if( !surface ) // in the case it is root set, create another set
1244  {
1245  rval = mdbImpl->create_meshset( MESHSET_SET, face );MB_CHK_SET_ERR( rval, "Failed to create a the new surface meshset" );
1246  }
1247  // set the geo tag
1248  rval = add_geo_set( face, 2 );MB_CHK_SET_ERR( rval, "Failed to add the geometry set to the tool" );
1249 
1250  // this will be our output set, will contain all our new geo sets
1251  rval = mdbImpl->create_meshset( MESHSET_SET, output );MB_CHK_SET_ERR( rval, "Failed to create the output meshset" );
1252 
1253  // add first geo set (face) to the output set
1254  rval = mdbImpl->add_entities( output, &face, 1 );MB_CHK_SET_ERR( rval, "Failed to add the new meshset to the output meshset" );
1255 
1256  // how many edges do we need to create?
1257  // depends on how many loops we have
1258  // also, we should avoid non-manifold topology
1259  if( !surface )
1260  { // in this case, surface is root, so we have to add entities
1261  rval = mdbImpl->add_entities( face, surface_ents );MB_CHK_SET_ERR( rval, "Failed to add surface entities to the surface meshset" );
1262  }
1263 
1264  Skinner tool( mdbImpl );
1265  rval = tool.find_skin( 0, surface_ents, 1, edge_ents );MB_CHK_SET_ERR( rval, "Failed to skin the surface entities" );
1266  if( debugFlag )
1267  {
1268  std::cout << "skinning edges: " << edge_ents.size() << "\n";
1269  for( Range::iterator it = edge_ents.begin(); it != edge_ents.end(); ++it )
1270  {
1271  EntityHandle ed = *it;
1272  std::cout << "edge: " << mdbImpl->id_from_handle( ed ) << " type:" << mdbImpl->type_from_handle( ed )
1273  << "\n";
1274  std::cout << mdbImpl->list_entity( ed );
1275  }
1276  }
1277 
1278  std::vector< EntityHandle > edges_loop;
1279 
1280  Range pool_of_edges = edge_ents;
1281  Range used_edges; // these edges are already used for some loops
1282  // get a new one
1283 
1284  while( !pool_of_edges.empty() )
1285  {
1286  // get the first edge, and start a loop with it
1287  EntityHandle current_edge = pool_of_edges[0];
1288  if( debugFlag )
1289  {
1290  std::cout << "Start current edge: " << mdbImpl->id_from_handle( current_edge ) << "\n ";
1291  std::cout << mdbImpl->list_entity( current_edge );
1292  }
1293  // get its triangle / quad and see its orientation
1294  std::vector< EntityHandle > tris;
1295  rval = mdbImpl->get_adjacencies( &current_edge, 1, 2, false, tris );MB_CHK_SET_ERR( rval, "Failed to get the adjacent triangles to the current edge" );
1296  if( tris.size() != 1 ) MB_SET_ERR( MB_FAILURE, "Edge not on boundary" );
1297 
1298  int side_n, sense, offset;
1299  rval = mdbImpl->side_number( tris[0], current_edge, side_n, sense, offset );MB_CHK_SET_ERR( rval, "Failed to get the current edge's side number" );
1300 
1301  const EntityHandle* conn2;
1302  int nnodes2;
1303  rval = mdbImpl->get_connectivity( current_edge, conn2, nnodes2 );MB_CHK_SET_ERR( rval, "Failed to get the current edge's connectivity" );
1304 
1305  if( nnodes2 != 2 ) MB_SET_ERR( MB_FAILURE, "Incorrect number of nodes found." );
1306 
1307  EntityHandle start_node = conn2[0];
1308  EntityHandle next_node = conn2[1];
1309 
1310  if( sense == -1 )
1311  {
1312  // revert the edge, and start well
1313  EntityHandle nn2[2] = { conn2[1], conn2[0] };
1314  rval = mdbImpl->set_connectivity( current_edge, nn2, 2 );MB_CHK_SET_ERR( rval, "Failed to set the connectivity of the current edge" );
1315 
1316  start_node = nn2[0]; // or conn2[0] !!! beware: conn2 is modified
1317  next_node = nn2[1]; // or conn2[1] !!!
1318  // reset connectivity of edge
1319  if( debugFlag ) std::cout << " current edge needs reversed\n";
1320  }
1321  // start a new loop of edges
1322  edges_loop.clear(); // every edge loop starts fresh
1323  edges_loop.push_back( current_edge );
1324  used_edges.insert( current_edge );
1325  pool_of_edges.erase( current_edge );
1326 
1327  if( debugFlag )
1328  {
1329  std::cout << " start node: " << start_node << "\n";
1330  std::cout << mdbImpl->list_entity( start_node );
1331  std::cout << " next node: " << next_node << "\n";
1332  std::cout << mdbImpl->list_entity( next_node );
1333  }
1334  while( next_node != start_node )
1335  {
1336  // find the next edge in the skin
1337  std::vector< EntityHandle > candidate_edges;
1338  rval = mdbImpl->get_adjacencies( &next_node, 1, 1, false, candidate_edges );MB_CHK_SET_ERR( rval, "Failed to get the adjacent edges to the next node" );
1339  // filter the edges that are used, or the edges not in the skin
1340  std::vector< EntityHandle > good_edges;
1341  for( int k = 0; k < (int)candidate_edges.size(); k++ )
1342  {
1343  EntityHandle edg = candidate_edges[k];
1344  if( used_edges.find( edg ) != used_edges.end() ) continue;
1345  if( pool_of_edges.find( edg ) != pool_of_edges.end() ) good_edges.push_back( edg );
1346  }
1347  if( good_edges.size() != 1 )
1348  {
1349  std::cout << " good_edges.size()=" << good_edges.size() << " STOP\n";
1350  MB_SET_ERR( MB_FAILURE, "Number of good edges is not one. Could not complete the loop" );
1351  }
1352  // see if the orientation is good; if not, revert it
1353 
1354  current_edge = good_edges[0];
1355  rval = mdbImpl->get_connectivity( current_edge, conn2, nnodes2 );MB_CHK_SET_ERR( rval, "Failed to get the connectivity of the current edge" );
1356  if( nnodes2 != 2 ) MB_SET_ERR( MB_FAILURE, "Incorrect number of nodes found" );
1357 
1358  if( conn2[0] != next_node )
1359  {
1360  if( conn2[1] != next_node )
1361  {
1362  // the edge is not connected then to current edge
1363  // bail out
1364  std::cout << "edge " << mdbImpl->id_from_handle( current_edge ) << " not connected to node "
1365  << next_node << "\n";
1366  MB_SET_ERR( MB_FAILURE, "Current edge is not connected to node" );
1367  ;
1368  }
1369  if( debugFlag )
1370  {
1371  std::cout << " revert edge " << mdbImpl->id_from_handle( current_edge ) << "\n";
1372  std::cout << mdbImpl->list_entity( current_edge );
1373  }
1374  // orientation should be reversed
1375  EntityHandle nn2[2] = { conn2[1], conn2[0] };
1376  rval = mdbImpl->set_connectivity( current_edge, nn2, 2 );MB_CHK_SET_ERR( rval, "Failed to set the connectivity of the current edge" );
1377 
1378  {
1379  std::cout << "after revert edge " << mdbImpl->id_from_handle( current_edge ) << "\n";
1380  std::cout << mdbImpl->list_entity( current_edge );
1381  std::cout << " conn2: " << conn2[0] << " " << conn2[1] << "\n";
1382  }
1383  }
1384  // before reversion, conn2 was something { n1, next_node}
1385  // after reversion, conn2 became {next_node, n1}, so the
1386  // new next node will be still conn2[1]; big surprise, as
1387  // I didn't expect the conn2 to change.
1388  // it seems that const EntityHandle * conn2 means conn2 cannot be
1389  // changed, but what is pointed to by it will change when we reset connectivity for edge
1390  next_node = conn2[1];
1391 
1392  if( debugFlag )
1393  {
1394  std::cout << " current edge: " << mdbImpl->id_from_handle( current_edge ) << "\n ";
1395  std::cout << mdbImpl->list_entity( current_edge );
1396  std::cout << "next node: " << next_node << "\n ";
1397  std::cout << mdbImpl->list_entity( next_node );
1398  }
1399  edges_loop.push_back( current_edge );
1400  used_edges.insert( current_edge );
1401  pool_of_edges.erase( current_edge );
1402  }
1403  // at this point, we have a loop formed;
1404  // create a geo edge, a vertex set, and add it to our sets
1405 
1407  rval = mdbImpl->create_meshset( MESHSET_ORDERED, edge );MB_CHK_SET_ERR( rval, "Failed to create the edge meshset" );
1408 
1409  rval = add_geo_set( edge, 1 );MB_CHK_SET_ERR( rval, "Failed to add the edge meshset to the tool's model set" );
1410  // add the mesh edges:
1411  // add loops edges to the edge set
1412  rval = mdbImpl->add_entities( edge, &edges_loop[0], edges_loop.size() ); //
1413  MB_CHK_SET_ERR( rval, "Failed to add entities to the edge meshset" );
1414  // create a vertex set
1416  rval = mdbImpl->create_meshset( MESHSET_SET, vertex );MB_CHK_SET_ERR( rval, "Failed to create the vertex meshset" );
1417  rval = add_geo_set( vertex, 0 );MB_CHK_SET_ERR( rval, "Failed to add the vertex meshset to the tool's model set" );
1418  // add one node to the vertex set
1419 
1420  rval = mdbImpl->add_entities( vertex, &start_node, 1 ); //
1421  MB_CHK_SET_ERR( rval, "Failed to add entities to the vertex meshset" );
1422 
1423  rval = mdbImpl->add_parent_child( face, edge );MB_CHK_SET_ERR( rval, "Failed to create the edge to face parent child relationship" );
1424 
1425  rval = mdbImpl->add_parent_child( edge, vertex );MB_CHK_SET_ERR( rval, "Failed to create the vertex to edge parent child relationship" );
1426 
1427  // the sense of the edge in face is for sure positive (forward)
1428  rval = set_sense( edge, face, 1 ); //
1429  MB_CHK_SET_ERR( rval, "Failed to set the edge to face sense" );
1430  // also add our sets to the output set, to be sure to be exported
1431 
1432  rval = mdbImpl->add_entities( output, &edge, 1 );MB_CHK_SET_ERR( rval, "Failed to add the edge meshset to the output set" );
1433  rval = mdbImpl->add_entities( output, &vertex, 1 );MB_CHK_SET_ERR( rval, "Failed to add the vertex meshset to the output set" );
1434 
1435  if( debugFlag )
1436  {
1437  std::cout << "add edge with start node " << start_node << " with " << edges_loop.size() << " edges\n";
1438  }
1439  }
1440 
1441  return MB_SUCCESS;
1442 }
1443 
1444 /*
1445  * This would create a deep copy of the geom topo model, into a new geom topo tool
1446  * sets will be duplicated, but entities not
1447  * modelSet will be a new one
1448  * if the original set was null (root), a new model set will be created for
1449  * original model, and its entities will be the original g sets
1450  * Scenario: split a face along a ground line, then write only one surface
1451  * the common line has 2 faces in var len tag for sense edge; if we write only one
1452  * surface to a new database, the var len tag of the edge will be extracted with 2 values, but
1453  * only one value will make sense, the other will be zero.
1454  *
1455  * There is no workaround; we need to create a duplicate model that has only that surface
1456  * and its children (and grand-children). Then the var len sense edge-face tags will have
1457  * the right size.
1458  *
1459  */
1460 ErrorCode GeomTopoTool::duplicate_model( GeomTopoTool*& duplicate, std::vector< EntityHandle >* pvGEnts )
1461 {
1462  // will
1463  EntityHandle rootModelSet;
1464  ErrorCode rval = mdbImpl->create_meshset( MESHSET_SET, rootModelSet );MB_CHK_SET_ERR( rval, "Failed to create the rootModelSet" );
1465 
1466  if( 0 == geomTag )
1467  {
1468  rval = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomTag );MB_CHK_SET_ERR( rval, "Failed to get the geometry dimension tag handle" );
1469  }
1470  if( 0 == gidTag )
1471  {
1473  }
1474  // extract from the geomSet the dimension, children, and grand-children
1475  Range depSets; // dependents of the geomSet, including the geomSet
1476  // add in this range all the dependents of this, to filter the ones that need to be deep copied
1477 
1478  if( pvGEnts != NULL )
1479  {
1480  unsigned int numGents = pvGEnts->size();
1481  for( unsigned int k = 0; k < numGents; k++ )
1482  {
1483  EntityHandle geomSet = ( *pvGEnts )[k];
1484  // will keep accumulating to the depSets range
1485  rval = mdbImpl->get_child_meshsets( geomSet, depSets, 0 ); // 0 for numHops means that all
1486  // dependents are returned, not only the direct children.
1487  MB_CHK_SET_ERR( rval, "Failed to get the geometry set's child meshsets" );
1488 
1489  depSets.insert( geomSet );
1490  }
1491  }
1492 
1493  // add to the root model set copies of the gsets, with correct sets
1494  // keep a map between sets to help in copying parent/child relations
1495  std::map< EntityHandle, EntityHandle > relate;
1496  // each set will get the same entities as the original
1497  for( int dim = 0; dim < 5; dim++ )
1498  {
1499  int gid = 0;
1500  unsigned int set_options = ( ( 1 != dim ) ? MESHSET_SET : MESHSET_ORDERED );
1501  for( Range::iterator it = geomRanges[dim].begin(); it != geomRanges[dim].end(); ++it )
1502  {
1503  EntityHandle set = *it;
1504  if( pvGEnts != NULL && depSets.find( set ) == depSets.end() )
1505  continue; // this means that this set is not of interest, skip it
1506  EntityHandle newSet;
1507  rval = mdbImpl->create_meshset( set_options, newSet );MB_CHK_SET_ERR( rval, "Failed to create new meshset" );
1508 
1509  relate[set] = newSet;
1510  rval = mdbImpl->add_entities( rootModelSet, &newSet, 1 );MB_CHK_SET_ERR( rval, "Failed to add the new meshset to the tool's modelSet" );
1511 
1512  // make it a geo set, and give also global id in order
1513  rval = mdbImpl->tag_set_data( geomTag, &newSet, 1, &dim );MB_CHK_SET_ERR( rval, "Failed to set the new meshset's geometry dimension data" );
1514 
1515  gid++; // increment global id, everything starts with 1 in the new model!
1516  rval = mdbImpl->tag_set_data( gidTag, &newSet, 1, &gid );MB_CHK_SET_ERR( rval, "Failed to get the new meshset's global id data" );
1517 
1518  if( dim == 1 )
1519  {
1520  // the entities are ordered, we need to retrieve them ordered, and set them ordered
1521  std::vector< EntityHandle > mesh_edges;
1522  rval = mdbImpl->get_entities_by_handle( set, mesh_edges );MB_CHK_SET_ERR( rval, "Failed to get the meshset entities by handle" );
1523 
1524  rval = mdbImpl->add_entities( newSet, &( mesh_edges[0] ), (int)mesh_edges.size() );MB_CHK_SET_ERR( rval, "Failed to add the new entities to the new meshset" );
1525  }
1526  else
1527  {
1528  Range ents;
1529  rval = mdbImpl->get_entities_by_handle( set, ents );MB_CHK_SET_ERR( rval, "Failed to add the entities to the existing meshset" );
1530 
1531  rval = mdbImpl->add_entities( newSet, ents );MB_CHK_SET_ERR( rval, "Failed to add the entities to the new meshset" );
1532  }
1533  // set parent/child relations if dim>=1
1534  if( dim >= 1 )
1535  {
1536  Range children;
1537  // the children of geo sets are only g sets
1538  rval = mdbImpl->get_child_meshsets( set, children ); // num_hops = 1 by default
1539  MB_CHK_SET_ERR( rval, "Failed to get the child meshsets of the existing set" );
1540 
1541  for( Range::iterator it2 = children.begin(); it2 != children.end(); ++it2 )
1542  {
1543  EntityHandle newChildSet = relate[*it2];
1544  rval = mdbImpl->add_parent_child( newSet, newChildSet );MB_CHK_SET_ERR( rval, "Failed to create parent child relationship to the new meshset" );
1545  }
1546  }
1547  }
1548  }
1549 
1550  duplicate = new GeomTopoTool( mdbImpl, true, rootModelSet ); // will retrieve the
1551  // sets and put them in ranges
1552 
1553  // this is the lazy way to it:
1554  // newgtt->restore_topology_from_adjacency(); // will reset the sense entities, and with this,
1555  // the model represented by this new gtt will be complete set senses by peeking at the old model
1556  // make sure we have the sense tags defined
1557  rval = check_face_sense_tag( true );MB_CHK_SET_ERR( rval, "Failed to check the face to volume sense tag handle" );
1558 
1559  rval = check_edge_sense_tags( true );MB_CHK_SET_ERR( rval, "Failed to check the curve to surface sense tag handles" );
1560 
1561  for( int dd = 1; dd <= 2; dd++ ) // do it for surfaces and edges
1562  {
1563  for( Range::iterator it = geomRanges[dd].begin(); it != geomRanges[dd].end(); ++it )
1564  {
1565  EntityHandle surf = *it;
1566  if( pvGEnts != NULL && depSets.find( surf ) == depSets.end() )
1567  continue; // this means that this set is not of interest, skip it
1568  EntityHandle newSurf = relate[surf];
1569  // we can actually look at the tag data, to be more efficient
1570  // or use the
1571  std::vector< EntityHandle > solids;
1572  std::vector< int > senses;
1573  rval = this->get_senses( surf, solids, senses );MB_CHK_SET_ERR( rval, "Failed to get the sense data for the surface with respect to its volumes" );
1574 
1575  std::vector< EntityHandle > newSolids;
1576  std::vector< int > newSenses;
1577  for( unsigned int i = 0; i < solids.size(); i++ )
1578  {
1579  if( pvGEnts != NULL && depSets.find( solids[i] ) == depSets.end() )
1580  continue; // this means that this set is not of interest, skip it
1581  EntityHandle newSolid = relate[solids[i]];
1582  // see which "solids" are in the new model
1583  newSolids.push_back( newSolid );
1584  newSenses.push_back( senses[i] );
1585  }
1586  rval = duplicate->set_senses( newSurf, newSolids, newSenses );MB_CHK_SET_ERR( rval, "Failed to set the sense data for the surface with respect to the new volumes" );
1587  }
1588  }
1589  // if the original root model set for this model is 0 (root set), then create
1590  // a new set and put all the old sets in the new model set
1591  // in this way, the existing gtt remains valid (otherwise, the modelSet would contain all the
1592  // gsets, the old ones and the new ones; the root set contains everything)
1593  if( modelSet == 0 )
1594  {
1595  rval = mdbImpl->create_meshset( MESHSET_SET, modelSet );MB_CHK_SET_ERR( rval, "Failed to create the modelSet meshset" );
1596 
1597  // add to this new set all previous sets (which are still in ranges)
1598  for( int dim = 0; dim < 5; dim++ )
1599  {
1600  rval = mdbImpl->add_entities( modelSet, geomRanges[dim] );MB_CHK_SET_ERR( rval, "Failed to add the geometric meshsets to the tool's modelSet" );
1601  }
1602  }
1603  return MB_SUCCESS;
1604 }
1605 
1607 {
1608  if( impl_compl_handle )
1609  {
1610  implicit_complement = impl_compl_handle;
1611  return MB_SUCCESS;
1612  }
1613  else
1614  {
1615  return MB_ENTITY_NOT_FOUND;
1616  }
1617 }
1618 
1620 {
1621 
1622  // if the implicit complement is already setup,
1623  // we're done
1624  if( impl_compl_handle != 0 )
1625  {
1626  std::cout << "IPC already exists!" << std::endl;
1627  return MB_SUCCESS;
1628  }
1629 
1630  // if not, then query for a set with it's name
1631  Range entities;
1632  const void* const tagdata[] = { IMPLICIT_COMPLEMENT_NAME };
1634  // query error
1635  MB_CHK_SET_ERR( rval, "Unable to query for implicit complement" );
1636 
1637  // if we found exactly one, set it as the implicit complement
1638  if( entities.size() == 1 )
1639  {
1640  impl_compl_handle = entities.front();
1641  return MB_SUCCESS;
1642  }
1643 
1644  // found too many
1645  if( entities.size() > 1 ) MB_CHK_SET_ERR( MB_MULTIPLE_ENTITIES_FOUND, "Too many implicit complement sets" );
1646 
1647  // found none
1648  if( entities.empty() )
1649  {
1650  // create implicit complement if requested
1651  rval = generate_implicit_complement( impl_compl_handle );MB_CHK_SET_ERR( rval, "Could not create implicit complement" );
1652 
1653  rval = mdbImpl->tag_set_data( nameTag, &impl_compl_handle, 1, &IMPLICIT_COMPLEMENT_NAME );MB_CHK_SET_ERR( rval, "Could not set the name tag for the implicit complement" );
1654 
1655  rval = add_geo_set( impl_compl_handle, 3 );MB_CHK_SET_ERR( rval, "Failed to add implicit complement to model" );
1656 
1657  // assign category tag - this is presumably for consistency so that the
1658  // implicit complement has all the appearance of being the same as any
1659  // other volume
1660  Tag category_tag;
1662  MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Could not get the category tag" );
1663 
1664  static const char volume_category[CATEGORY_TAG_SIZE] = "Volume\0";
1665  rval = mdbImpl->tag_set_data( category_tag, &impl_compl_handle, 1, volume_category );MB_CHK_SET_ERR( rval, "Could not set the category tag for the implicit complement" );
1666 
1667  return MB_SUCCESS;
1668  }
1669 
1670  return MB_FAILURE;
1671 }
1672 
1674 {
1675 
1676  ErrorCode rval;
1677  rval = mdbImpl->create_meshset( MESHSET_SET, implicit_complement_set );MB_CHK_SET_ERR( rval, "Failed to create mesh set for implicit complement" );
1678 
1679  // make sure the sense2Tag is set
1680  if( !sense2Tag )
1681  {
1682  check_face_sense_tag( true );
1683  }
1684 
1685  // get all geometric surface sets
1686  Range surfs;
1687  rval = get_gsets_by_dimension( 2, surfs );MB_CHK_SET_ERR( rval, "Could not get surface sets" );
1688 
1689  // search through all surfaces
1690  std::vector< EntityHandle > parent_vols;
1691  for( Range::iterator surf_i = surfs.begin(); surf_i != surfs.end(); ++surf_i )
1692  {
1693 
1694  parent_vols.clear();
1695  // get parents of each surface
1696  rval = mdbImpl->get_parent_meshsets( *surf_i, parent_vols );MB_CHK_SET_ERR( rval, "Failed to get volume meshsets" );
1697 
1698  // if only one parent, get the OBB root for this surface
1699  if( parent_vols.size() == 1 )
1700  {
1701 
1702  // add this surf to the topology of the implicit complement volume
1703  rval = mdbImpl->add_parent_child( implicit_complement_set, *surf_i );MB_CHK_SET_ERR( rval, "Could not add surface to implicit complement set" );
1704 
1705  // get the surface sense wrt original volume
1706  EntityHandle sense_data[2] = { 0, 0 };
1707  rval = get_surface_senses( *surf_i, sense_data[0], sense_data[1] );MB_CHK_SET_ERR( rval, "Could not get surface sense data" );
1708 
1709  // set the surface sense wrt implicit complement volume
1710  if( 0 == sense_data[0] && 0 == sense_data[1] )
1711  MB_SET_ERR( MB_FAILURE, "No sense data for current surface" );
1712  if( 0 == sense_data[0] )
1713  sense_data[0] = implicit_complement_set;
1714  else if( 0 == sense_data[1] )
1715  sense_data[1] = implicit_complement_set;
1716  else
1717  MB_SET_ERR( MB_FAILURE, "Could not insert implicit complement into surface sense data" );
1718 
1719  // set the new sense data for this surface
1720  rval = set_surface_senses( *surf_i, sense_data[0], sense_data[1] );MB_CHK_SET_ERR( rval, "Failed to set sense tag data" );
1721  }
1722  } // end surface loop
1723 
1724  return MB_SUCCESS;
1725 }
1726 
1727 #define RETFALSE( a, b ) \
1728  { \
1729  std::cout << ( a ) << "\n"; \
1730  mdbImpl->list_entity( b ); \
1731  return false; \
1732  }
1734 {
1735  // vertex sets should have one node
1736  Range::iterator rit;
1737  ErrorCode rval;
1738  for( rit = geomRanges[0].begin(); rit != geomRanges[0].end(); ++rit )
1739  {
1740  EntityHandle vSet = *rit;
1741  Range nodes;
1742  rval = mdbImpl->get_entities_by_handle( vSet, nodes );
1743  if( MB_SUCCESS != rval ) RETFALSE( " failed to get nodes from vertex set ", vSet );
1744  if( nodes.size() != 1 ) RETFALSE( " number of nodes is different from 1 ", vSet )
1745  EntityType type = mdbImpl->type_from_handle( nodes[0] );
1746  if( type != MBVERTEX ) RETFALSE( " entity in vertex set is not a node ", nodes[0] )
1747  // get all parents, and see if they belong to geomRanges[1]
1748  Range edges;
1749  rval = mdbImpl->get_parent_meshsets( vSet, edges );
1750  if( MB_SUCCESS != rval ) RETFALSE( " can't get parent edges for a node set ", vSet )
1751  Range notEdges = subtract( edges, geomRanges[1] );
1752  if( !notEdges.empty() ) RETFALSE( " some parents of a node set are not geo edges ", notEdges[0] )
1753  }
1754 
1755  // edges to be formed by continuous chain of mesh edges, oriented correctly
1756  for( rit = geomRanges[1].begin(); rit != geomRanges[1].end(); ++rit )
1757  {
1758  EntityHandle edge = *rit;
1759  std::vector< EntityHandle > mesh_edges;
1760  rval = mdbImpl->get_entities_by_type( edge, MBEDGE, mesh_edges );
1761  if( MB_SUCCESS != rval ) RETFALSE( " can't get mesh edges from edge set", edge )
1762  int num_edges = (int)mesh_edges.size();
1763  if( num_edges == 0 ) RETFALSE( " no mesh edges in edge set ", edge )
1764  EntityHandle firstNode;
1765  EntityHandle currentNode; // will also hold the last node in chain of edges
1766  const EntityHandle* conn2;
1767  int nnodes2;
1768  // get all parents, and see if they belong to geomRanges[1]
1769  for( int i = 0; i < num_edges; i++ )
1770  {
1771  rval = mdbImpl->get_connectivity( mesh_edges[i], conn2, nnodes2 );
1772  if( MB_SUCCESS != rval || nnodes2 != 2 ) RETFALSE( " mesh edge connectivity is wrong ", mesh_edges[i] )
1773  if( i == 0 )
1774  {
1775  firstNode = conn2[0];
1776  currentNode = conn2[1];
1777  }
1778 
1779  else // if ( (i>0) )
1780  {
1781  // check the current node is conn[0]
1782  if( conn2[0] != currentNode )
1783  {
1784  std::cout << "i=" << i << " conn2:" << conn2[0] << " " << conn2[1] << " currentNode:" << currentNode
1785  << "\n";
1786  mdbImpl->list_entity( mesh_edges[i] );
1787  RETFALSE( " edges are not contiguous in edge set ", edge )
1788  }
1789  currentNode = conn2[1];
1790  }
1791  }
1792  // check the children of the edge set; do they have the first and last node?
1793  Range vertSets;
1794  rval = mdbImpl->get_child_meshsets( edge, vertSets );
1795  if( MB_SUCCESS != rval ) RETFALSE( " can't get vertex children ", edge )
1796  Range notVertices = subtract( vertSets, geomRanges[0] );
1797  if( !notVertices.empty() ) RETFALSE( " children sets that are not vertices ", notVertices[0] )
1798  for( Range::iterator it = vertSets.begin(); it != vertSets.end(); ++it )
1799  {
1800  if( !mdbImpl->contains_entities( *it, &firstNode, 1 ) &&
1801  !mdbImpl->contains_entities( *it, &currentNode, 1 ) )
1802  RETFALSE( " a vertex set is not containing the first and last nodes ", *it )
1803  }
1804  // check now the faces / parents
1805  Range faceSets;
1806  rval = mdbImpl->get_parent_meshsets( edge, faceSets );
1807  if( MB_SUCCESS != rval ) RETFALSE( " can't get edge parents ", edge )
1808  Range notFaces = subtract( faceSets, geomRanges[2] );
1809  if( !notFaces.empty() ) RETFALSE( " parent sets that are not faces ", notFaces[0] )
1810 
1811  // for a geo edge, check the sense tags with respect to the adjacent faces
1812  // in general, it is sufficient to check one mesh edge (the first one)
1813  // edge/face senses
1814  EntityHandle firstMeshEdge = mesh_edges[0];
1815  // get all entities/elements adjacent to it
1816  Range adjElem;
1817  rval = mdbImpl->get_adjacencies( &firstMeshEdge, 1, 2, false, adjElem );
1818  if( MB_SUCCESS != rval ) RETFALSE( " can't get adjacent elements to the edge ", firstMeshEdge )
1819  for( Range::iterator it2 = adjElem.begin(); it2 != adjElem.end(); ++it2 )
1820  {
1821  EntityHandle elem = *it2;
1822  // find the geo face it belongs to
1823  EntityHandle gFace = 0;
1824  for( Range::iterator fit = faceSets.begin(); fit != faceSets.end(); ++fit )
1825  {
1826  EntityHandle possibleFace = *fit;
1827  if( mdbImpl->contains_entities( possibleFace, &elem, 1 ) )
1828  {
1829  gFace = possibleFace;
1830  break;
1831  }
1832  }
1833  if( 0 == gFace )
1834  RETFALSE( " can't find adjacent surface that contains the adjacent element to the edge ",
1835  firstMeshEdge )
1836 
1837  // now, check the sense of mesh_edge in element, and the sense of gedge in gface
1838  // side_number
1839  int side_n, sense, offset;
1840  rval = mdbImpl->side_number( elem, firstMeshEdge, side_n, sense, offset );
1841  if( MB_SUCCESS != rval ) RETFALSE( " can't get sense and side number of an element ", elem )
1842  // now get the sense
1843  int topoSense;
1844  rval = this->get_sense( edge, gFace, topoSense );
1845  if( topoSense != sense ) RETFALSE( " geometric topo sense and element sense do not agree ", edge )
1846  }
1847  }
1848  // surfaces to be true meshes
1849  // for surfaces, check that the skinner will find the correct boundary
1850 
1851  // use the skinner for boundary check
1852  Skinner tool( mdbImpl );
1853 
1854  for( rit = geomRanges[2].begin(); rit != geomRanges[2].end(); ++rit )
1855  {
1856  EntityHandle faceSet = *rit;
1857  // get all boundary edges (adjacent edges)
1858 
1859  Range edges;
1860  rval = mdbImpl->get_child_meshsets( faceSet, edges );
1861  if( MB_SUCCESS != rval ) RETFALSE( " can't get children edges for a face set ", faceSet )
1862  Range notEdges = subtract( edges, geomRanges[1] );
1863  if( !notEdges.empty() ) RETFALSE( " some children of a face set are not geo edges ", notEdges[0] )
1864 
1865  Range boundary_mesh_edges;
1866  for( Range::iterator it = edges.begin(); it != edges.end(); ++it )
1867  {
1868  rval = mdbImpl->get_entities_by_type( *it, MBEDGE, boundary_mesh_edges );
1869  if( MB_SUCCESS != rval ) RETFALSE( " can't get edge elements from the edge set ", *it )
1870  }
1871  // skin the elements of the surface
1872  // most of these should be triangles and quads
1873  Range surface_ents, edge_ents;
1874  rval = mdbImpl->get_entities_by_dimension( faceSet, 2, surface_ents );
1875  if( MB_SUCCESS != rval ) RETFALSE( " can't get surface elements from the face set ", faceSet )
1876 
1877  rval = tool.find_skin( 0, surface_ents, 1, edge_ents );
1878  if( MB_SUCCESS != rval ) RETFALSE( "can't skin a surface ", surface_ents[0] )
1879 
1880  // those 2 ranges for boundary edges now must be equal
1881  if( boundary_mesh_edges != edge_ents ) RETFALSE( "boundary ranges are different", boundary_mesh_edges[0] )
1882  }
1883 
1884  // solids to be filled correctly, maybe a skin procedure too.
1885  // (maybe the solids are empty)
1886 
1887  return true;
1888 }
1889 
1891 {
1892  return rootSets.size() != 0 || mapRootSets.size() != 0;
1893 }
1894 
1895 // This function gets coordinates of the minimum and maxmiumum points
1896 // from an OBB/AABB, ie. such that these points represent
1897 // the maximum and minium extents of an AABB
1898 ErrorCode GeomTopoTool::get_bounding_coords( EntityHandle volume, double minPt[3], double maxPt[3] )
1899 {
1900  double center[3], axis1[3], axis2[3], axis3[3];
1901 
1902  // get center point and vectors to OBB faces
1903  ErrorCode rval = get_obb( volume, center, axis1, axis2, axis3 );MB_CHK_SET_ERR( rval, "Failed to get the oriented bounding box of the volume" );
1904 
1905  // compute min and max vertices
1906  for( int i = 0; i < 3; i++ )
1907  {
1908  double sum = fabs( axis1[i] ) + fabs( axis2[i] ) + fabs( axis3[i] );
1909  minPt[i] = center[i] - sum;
1910  maxPt[i] = center[i] + sum;
1911  }
1912  return MB_SUCCESS;
1913 }
1914 
1916  double center[3],
1917  double axis1[3],
1918  double axis2[3],
1919  double axis3[3] )
1920 {
1921  // find EntityHandle node_set for use in box
1922  EntityHandle root;
1923  ErrorCode rval = get_root( volume, root );MB_CHK_SET_ERR( rval, "Failed to get volume's obb tree root" );
1924 
1925  // call box to get center and vectors to faces
1926  return obbTree->box( root, center, axis1, axis2, axis3 );
1927 }
1928 
1930 {
1931  Range all_children, desired_children;
1932  Range::iterator it;
1933  int actual_dimension;
1934 
1935  desired_children.clear();
1936  all_children.clear();
1937  mdbImpl->get_child_meshsets( parent, all_children );
1938 
1939  for( it = all_children.begin(); it != all_children.end(); ++it )
1940  {
1941  mdbImpl->tag_get_data( geomTag, &( *it ), 1, &actual_dimension );
1942  if( actual_dimension == desired_dimension ) desired_children.insert( *it );
1943  }
1944 
1945  return desired_children;
1946 }
1947 
1948 // runs GeomQueryTool point_in_vol and to test if vol A is inside vol B
1949 // returns true or false
1951 {
1952  ErrorCode rval;
1953 
1954  Range child_surfaces, triangles, vertices;
1955  double coord[3]; // coord[0] = x, etc.
1956  int result; // point in vol result; 0=F, 1=T
1957 
1958  // find coordinates of any point on surface of A
1959  // get surface corresponding to volume, then get the triangles
1960  child_surfaces = get_ct_children_by_dimension( volume_A, 2 );
1961  rval = mdbImpl->get_entities_by_type( *child_surfaces.begin(), MBTRI, triangles );MB_CHK_ERR( rval );
1962 
1963  // now get 1st triangle vertices
1964  rval = mdbImpl->get_connectivity( &( *triangles.begin() ), 1, vertices );MB_CHK_ERR( rval );
1965 
1966  // now get coordinates of first vertex
1967  rval = mdbImpl->get_coords( &( *vertices.begin() ), 1, &( coord[0] ) );MB_CHK_ERR( rval );
1968 
1969  // if point on A is inside vol B, return T; o.w. return F
1970  rval = GQT->point_in_volume( volume_B, coord, result );MB_CHK_SET_ERR( rval, "Failed to complete point in volume query." );
1971 
1972  return ( result != 0 );
1973 }
1974 
1976 {
1977  ErrorCode rval;
1978 
1979  bool inserted = false;
1980  EntityHandle current_volume = volume; // volume to be inserted
1981  EntityHandle tree_volume = ct_root; // volume already existing in the tree
1982  EntityHandle parent = ct_root;
1983  Range child_volumes;
1984 
1985  // while not inserted in tree
1986  while( !inserted )
1987  {
1988  // if current volume is insde of tree volume-- always true if tree volume
1989  // is the root of the tree
1990  if( tree_volume == ct_root || ( tree_volume != ct_root && A_is_in_B( current_volume, tree_volume, GQT ) ) )
1991  {
1992 
1993  parent = tree_volume;
1994 
1995  // if tree_volume has children then we must test them,
1996  // (tree_volume will change)
1997  child_volumes = get_ct_children_by_dimension( tree_volume, 3 );
1998  if( child_volumes.size() > 0 ) tree_volume = child_volumes.pop_front();
1999  // otherwise current_volume is the only child of the tree volume
2000  else
2001  {
2002  rval = mdbImpl->add_parent_child( parent, current_volume );MB_CHK_SET_ERR( rval, "Failed to add parent-child relationship." );
2003 
2004  inserted = true;
2005  }
2006  // if current volume is not in the tree volume, the converse may be true
2007  }
2008  else
2009  {
2010  // if the tree volume is inside the current volume
2011  if( A_is_in_B( tree_volume, current_volume, GQT ) )
2012  {
2013  // reverse their parentage
2014  rval = mdbImpl->remove_parent_child( parent, tree_volume );MB_CHK_SET_ERR( rval, "Failed to remove parent-child relationship." );
2015  rval = mdbImpl->add_parent_child( current_volume, tree_volume );MB_CHK_SET_ERR( rval, "Failed to add parent-child relationship." );
2016  }
2017 
2018  if( child_volumes.size() == 0 )
2019  {
2020  rval = mdbImpl->add_parent_child( parent, current_volume );MB_CHK_SET_ERR( rval, "Failed to add parent-child relationship." );
2021  inserted = true;
2022  }
2023  else
2024  tree_volume = child_volumes.pop_front();
2025  }
2026  }
2027  return MB_SUCCESS;
2028 }
2029 
2031 {
2032 
2033  ErrorCode rval;
2034  // local var will go out of scope if errors appear, no need to free it also
2035  GeomQueryTool GQT( this );
2036  std::map< EntityHandle, EntityHandle > volume_surface; // map of volume
2037  // to its surface
2038 
2039  EntityHandle ct_root;
2040  // create root meshset-- this will be top of tree
2041  std::string meshset_name = "build_hierarchy_root";
2042  rval = mdbImpl->create_meshset( MESHSET_SET, ct_root );MB_CHK_ERR( rval );
2043  rval = mdbImpl->tag_set_data( nameTag, &ct_root, 1, meshset_name.c_str() );MB_CHK_ERR( rval );
2044 
2045  for( Range::iterator vol = flat_volumes.begin(); vol != flat_volumes.end(); vol++ )
2046  {
2047  // get the surface corresponding to each volume
2048  // at this point, each volume meshset only has one 'child' surface
2049  // which exactly corresponds to that volume
2050  Range child_surfaces = get_ct_children_by_dimension( *vol, 2 );
2051  volume_surface[*vol] = *child_surfaces.begin();
2052 
2053  rval = insert_in_tree( ct_root, *vol, &GQT );MB_CHK_SET_ERR( rval, "Failed to insert volume into tree." );
2054  }
2055 
2056  // for each original volume, get its child volumes
2057  for( Range::iterator parent_it = flat_volumes.begin(); parent_it != flat_volumes.end(); parent_it++ )
2058  {
2059  Range volume_children = get_ct_children_by_dimension( *parent_it, 3 );
2060 
2061  if( volume_children.size() != 0 )
2062  {
2063  // loop over all of original volume's child volumes
2064  for( Range::iterator child_it = volume_children.begin(); child_it != volume_children.end(); ++child_it )
2065  {
2066  // set the sense of the surface mapped to the child volume to REVERSE
2067  // wrt the parent volume
2068  rval = set_sense( volume_surface[*child_it], *parent_it, SENSE_REVERSE );MB_CHK_SET_ERR( rval, "Failed to set sense." );
2069 
2070  // add the child volume's surface as a child of the original volume
2071  // and delete the child volume as a child of original volume
2072  rval = mdbImpl->add_parent_child( *parent_it, volume_surface[*child_it] );MB_CHK_SET_ERR( rval, "Failed to add parent-child relationship." );
2073  rval = mdbImpl->remove_parent_child( *parent_it, *child_it );MB_CHK_SET_ERR( rval, "Failed to remove parent-child relationship." );
2074  }
2075  }
2076  }
2077 
2078  return MB_SUCCESS;
2079 }
2080 
2081 } // namespace moab