 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /**
2 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3 * storing and accessing finite element mesh data.
4 *
5 * Copyright 2004 Sandia Corporation. Under the terms of Contract
6 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7 * retains certain rights in this software.
8 *
9 * This library is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public
11 * License as published by the Free Software Foundation; either
12 * version 2.1 of the License, or (at your option) any later version.
13 *
14 */
15
16 #include "moab/GeomTopoTool.hpp"
17 #include "moab/GeomQueryTool.hpp"
18 #include "moab/OrientedBoxTreeTool.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
46 GeomTopoTool::GeomTopoTool( Interface* impl,
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
63 gidTag = mdbImpl->globalId_tag();
64
65 rval =
66 mdbImpl->tag_get_handle( NAME_TAG_NAME, NAME_TAG_SIZE, MB_TYPE_OPAQUE, nameTag, MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR_CONT( rval, "Error: Failed to create name tag" );
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
73 impl_compl_handle = 0;
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
91 GeomTopoTool::~GeomTopoTool()
92 {
93 delete obbTree;
94 }
95
96 int GeomTopoTool::dimension( EntityHandle this_set )
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
122 int GeomTopoTool::global_id( EntityHandle this_set )
123 {
124 ErrorCode result;
125 if( 0 == gidTag )
126 {
127 gidTag = mdbImpl->globalId_tag();
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
168 ErrorCode GeomTopoTool::other_entity( EntityHandle bounded,
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
206 ErrorCode GeomTopoTool::restore_obb_index()
207 {
208
209 if( m_rootSets_vector ) resize_rootSets();
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
230 ErrorCode GeomTopoTool::find_geomsets( Range* ranges )
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
256 ErrorCode GeomTopoTool::get_gsets_by_dimension( int dim, Range& gset )
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
267 ErrorCode GeomTopoTool::resize_rootSets()
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
306 ErrorCode GeomTopoTool::is_owned_set( EntityHandle eh )
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
318 ErrorCode GeomTopoTool::delete_obb_tree( EntityHandle gset, bool vol_only )
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
395 ErrorCode GeomTopoTool::delete_all_obb_trees()
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
413 ErrorCode GeomTopoTool::construct_obb_tree( EntityHandle eh )
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
504 ErrorCode GeomTopoTool::set_root_set( EntityHandle vol_or_surf, EntityHandle root )
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
523 ErrorCode GeomTopoTool::remove_root( EntityHandle vol_or_surf )
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
566 ErrorCode GeomTopoTool::construct_obb_trees( bool make_one_vol )
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
605 ErrorCode GeomTopoTool::restore_topology_from_adjacency()
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
712 ErrorCode GeomTopoTool::separate_by_dimension( const Range& geom_sets )
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
722 Range::const_iterator git;
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 {
742 gidTag = mdbImpl->globalId_tag();
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.
810 ErrorCode GeomTopoTool::set_sense( EntityHandle entity, EntityHandle wrt_entity, int sense )
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
905 ErrorCode GeomTopoTool::get_sense( EntityHandle entity, EntityHandle wrt_entity, int& sense )
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
950 ErrorCode GeomTopoTool::get_surface_senses( EntityHandle surface_ent,
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
975 ErrorCode GeomTopoTool::set_surface_senses( EntityHandle surface_ent,
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
997 ErrorCode GeomTopoTool::get_surface_senses( EntityHandle 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
1017 ErrorCode GeomTopoTool::get_senses( EntityHandle entity,
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
1095 ErrorCode GeomTopoTool::set_senses( EntityHandle entity,
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
1108 ErrorCode GeomTopoTool::next_vol( EntityHandle surface, EntityHandle old_volume, EntityHandle& new_volume )
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
1128 ErrorCode GeomTopoTool::check_geom_tag( bool create )
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
1140 ErrorCode GeomTopoTool::check_gid_tag( bool create )
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
1154 ErrorCode GeomTopoTool::check_face_sense_tag( bool create )
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
1167 ErrorCode GeomTopoTool::check_edge_sense_tags( bool create )
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
1180 ErrorCode GeomTopoTool::add_geo_set( EntityHandle set, int dim, int gid )
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 {
1199 gidTag = mdbImpl->globalId_tag();
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
1226 ErrorCode GeomTopoTool::geometrize_surface_set( EntityHandle surface, EntityHandle& 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( ¤t_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
1406 EntityHandle edge;
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
1415 EntityHandle vertex;
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 {
1472 gidTag = mdbImpl->globalId_tag();
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
1606 ErrorCode GeomTopoTool::get_implicit_complement( EntityHandle& implicit_complement )
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
1619 ErrorCode GeomTopoTool::setup_implicit_complement()
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 };
1633 ErrorCode rval = mdbImpl->get_entities_by_type_and_tag( modelSet, MBENTITYSET, &nameTag, tagdata, 1, entities );
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;
1661 rval = mdbImpl->tag_get_handle( CATEGORY_TAG_NAME, CATEGORY_TAG_SIZE, MB_TYPE_OPAQUE, 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
1673 ErrorCode GeomTopoTool::generate_implicit_complement( EntityHandle& implicit_complement_set )
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 }
1733 bool GeomTopoTool::check_model()
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, ¤tNode, 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
1890 bool GeomTopoTool::have_obb_tree()
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
1915 ErrorCode GeomTopoTool::get_obb( EntityHandle volume,
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
1929 Range GeomTopoTool::get_ct_children_by_dimension( EntityHandle parent, int desired_dimension )
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
1950 bool GeomTopoTool::A_is_in_B( EntityHandle volume_A, EntityHandle volume_B, GeomQueryTool* GQT )
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
1975 ErrorCode GeomTopoTool::insert_in_tree( EntityHandle ct_root, EntityHandle volume, GeomQueryTool* GQT )
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
2030 ErrorCode GeomTopoTool::restore_topology_from_geometric_inclusion( const Range& flat_volumes )
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
1.9.1.