 |
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/DualTool.hpp"
17 #include "moab/Range.hpp"
18 // using Core for call to check_adjacencies
19 #include "moab/Core.hpp"
20 #include "Internals.hpp"
21 #include "MBTagConventions.hpp"
22 #include "moab/Skinner.hpp"
23 #include "moab/Core.hpp"
24 #include "moab/MeshTopoUtil.hpp"
25 #include "AEntityFactory.hpp"
26 #include "moab/CN.hpp"
27 #include <string>
28 #include <algorithm>
29 #include <iostream>
30 #include <sstream>
31 #include <cassert>
32
33 #define RR \
34 if( MB_SUCCESS != result ) return result
35 #define SWAP( a, b ) \
36 { \
37 EntityHandle tmp_ent = a; \
38 ( a ) = b; \
39 ( b ) = tmp_ent; \
40 }
41
42 namespace moab
43 {
44
45 bool debug = false;
46 bool debug_ap = false;
47
48 //! tag name for dual surfaces
49 const char* DualTool::DUAL_SURFACE_TAG_NAME = "DUAL_SURFACE";
50
51 //! tag name for dual curves
52 const char* DualTool::DUAL_CURVE_TAG_NAME = "DUAL_CURVE";
53
54 //! tag name for dual cells
55 const char* DualTool::IS_DUAL_CELL_TAG_NAME = "__IS_DUAL_CELL";
56
57 //! tag name for dual entities
58 const char* DualTool::DUAL_ENTITY_TAG_NAME = "__DUAL_ENTITY";
59
60 //! tag name for extra dual entities
61 const char* DualTool::EXTRA_DUAL_ENTITY_TAG_NAME = "__EXTRA_DUAL_ENTITY";
62
63 //! tag name for graphics point
64 const char* DualTool::DUAL_GRAPHICS_POINT_TAG_NAME = "__DUAL_GRAPHICS_POINT";
65
66 // const int DualTool::GP_SIZE = 20;
67
68 DualTool::DualTool( Interface* impl ) : mbImpl( impl )
69 {
70 EntityHandle dum_handle = 0;
71 ErrorCode result;
72
73 result = mbImpl->tag_get_handle( DUAL_SURFACE_TAG_NAME, 1, MB_TYPE_HANDLE, dualSurfaceTag,
74 MB_TAG_SPARSE | MB_TAG_CREAT, &dum_handle );
75 assert( MB_SUCCESS == result );
76
77 result = mbImpl->tag_get_handle( DUAL_CURVE_TAG_NAME, 1, MB_TYPE_HANDLE, dualCurveTag, MB_TAG_SPARSE | MB_TAG_CREAT,
78 &dum_handle );
79 assert( MB_SUCCESS == result );
80
81 unsigned int dummy = 0;
82 result = mbImpl->tag_get_handle( IS_DUAL_CELL_TAG_NAME, 1, MB_TYPE_INTEGER, isDualCellTag,
83 MB_TAG_SPARSE | MB_TAG_CREAT, &dummy );
84 assert( MB_SUCCESS == result );
85
86 result = mbImpl->tag_get_handle( DUAL_ENTITY_TAG_NAME, 1, MB_TYPE_HANDLE, dualEntityTag,
87 MB_TAG_DENSE | MB_TAG_CREAT, &dum_handle );
88 assert( MB_SUCCESS == result );
89
90 result = mbImpl->tag_get_handle( EXTRA_DUAL_ENTITY_TAG_NAME, 1, MB_TYPE_HANDLE, extraDualEntityTag,
91 MB_TAG_SPARSE | MB_TAG_CREAT, &dum_handle );
92 assert( MB_SUCCESS == result );
93
94 static const char dum_name[CATEGORY_TAG_SIZE] = { 0 };
95 result = mbImpl->tag_get_handle( CATEGORY_TAG_NAME, CATEGORY_TAG_SIZE, MB_TYPE_OPAQUE, categoryTag,
96 MB_TAG_SPARSE | MB_TAG_CREAT, dum_name );
97 assert( MB_SUCCESS == result );
98
99 DualTool::GraphicsPoint dum_pt( 0.0, 0.0, 0.0, -1 );
100 result = mbImpl->tag_get_handle( DUAL_GRAPHICS_POINT_TAG_NAME, sizeof( DualTool::GraphicsPoint ), MB_TYPE_DOUBLE,
101 dualGraphicsPointTag, MB_TAG_DENSE | MB_TAG_CREAT | MB_TAG_BYTES, &dum_pt );
102 assert( MB_SUCCESS == result );
103
104 globalIdTag = mbImpl->globalId_tag();
105
106 if( MB_SUCCESS == result )
107 {
108 } // empty statement to get rid of warning.
109
110 maxHexId = -1;
111 }
112
113 DualTool::~DualTool() {}
114
115 //! construct the dual entities for the entire mesh
116 ErrorCode DualTool::construct_dual( EntityHandle* entities, const int num_entities )
117 {
118 // allocate a dual entity for each primal entity in the mesh, starting
119 // with highest dimension and working downward; do each dimension in a separate code
120 // block, since they're all handled slightly differently
121
122 Range regions, faces, edges, vertices;
123 ErrorCode result;
124
125 if( NULL == entities || 0 == num_entities )
126 {
127
128 // first, construct all the aentities, since they're currently needed to
129 // compute the dual
130 result = mbImpl->get_entities_by_dimension( 0, 0, vertices );
131 if( MB_SUCCESS != result ) return result;
132
133 result = MeshTopoUtil( mbImpl ).construct_aentities( vertices );
134 if( MB_SUCCESS != result ) return result;
135
136 // get all edges, faces and regions now, so we don't need to filter out dual
137 // entities later
138
139 result = mbImpl->get_entities_by_dimension( 0, 1, edges );
140 if( MB_SUCCESS != result ) return result;
141 result = mbImpl->get_entities_by_dimension( 0, 2, faces );
142 if( MB_SUCCESS != result ) return result;
143 result = mbImpl->get_entities_by_dimension( 0, 3, regions );
144 if( MB_SUCCESS != result ) return result;
145
146 // get the max global id for hexes, we'll need for modification ops
147 std::vector< int > gid_vec( regions.size() );
148 result = mbImpl->tag_get_data( globalId_tag(), regions, &gid_vec[0] );
149 if( MB_SUCCESS != result ) return result;
150 maxHexId = -1;
151 Range::iterator rit;
152 unsigned int i;
153 for( rit = regions.begin(), i = 0; rit != regions.end(); ++rit, i++ )
154 {
155 if( gid_vec[i] > maxHexId && mbImpl->type_from_handle( *rit ) == MBHEX ) maxHexId = gid_vec[i];
156 }
157 }
158 else
159 {
160 // get entities of various dimensions adjacent to these
161 result = mbImpl->get_adjacencies( entities, num_entities, 0, true, vertices, Interface::UNION );
162 if( MB_SUCCESS != result ) return result;
163 result = mbImpl->get_adjacencies( entities, num_entities, 1, true, edges, Interface::UNION );
164 if( MB_SUCCESS != result ) return result;
165 result = mbImpl->get_adjacencies( entities, num_entities, 2, true, faces, Interface::UNION );
166 if( MB_SUCCESS != result ) return result;
167 result = mbImpl->get_adjacencies( entities, num_entities, 3, true, regions, Interface::UNION );
168 if( MB_SUCCESS != result ) return result;
169 }
170
171 Range dual_verts;
172 result = construct_dual_vertices( regions, dual_verts );
173 if( MB_SUCCESS != result || dual_verts.size() != regions.size() ) return result;
174 if( debug ) std::cout << "Constructed " << dual_verts.size() << " dual vertices." << std::endl;
175
176 // don't really need dual edges, but construct 'em anyway
177 Range dual_edges;
178 result = construct_dual_edges( faces, dual_edges );
179 if( MB_SUCCESS != result || dual_edges.size() != faces.size() ) return result;
180 if( debug ) std::cout << "Constructed " << dual_edges.size() << " dual edges." << std::endl;
181
182 // construct dual faces
183 Range dual_faces;
184 result = construct_dual_faces( edges, dual_faces );
185 if( MB_SUCCESS != result || dual_faces.size() != edges.size() ) return result;
186 if( debug ) std::cout << "Constructed " << dual_faces.size() << " dual faces." << std::endl;
187
188 // construct dual cells
189 Range dual_cells;
190 result = construct_dual_cells( vertices, dual_cells );
191 if( MB_SUCCESS != result || dual_cells.size() != vertices.size() ) return result;
192 if( debug ) std::cout << "Constructed " << dual_cells.size() << " dual cells." << std::endl;
193
194 return MB_SUCCESS;
195 }
196
197 ErrorCode DualTool::construct_dual_vertices( const Range& all_regions, Range& dual_ents )
198 {
199 if( all_regions.empty() ) return MB_SUCCESS;
200
201 // make sure they're all regions
202 assert( 3 == CN::Dimension( TYPE_FROM_HANDLE( *all_regions.begin() ) ) &&
203 3 == CN::Dimension( TYPE_FROM_HANDLE( *all_regions.rbegin() ) ) );
204
205 Range::const_iterator rit;
206 EntityHandle dual_ent;
207 ErrorCode tmp_result = MB_SUCCESS;
208 ErrorCode result = MB_SUCCESS;
209
210 for( rit = all_regions.begin(); rit != all_regions.end(); ++rit )
211 {
212 if( tmp_result != MB_SUCCESS ) result = tmp_result;
213
214 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
215 if( MB_SUCCESS == tmp_result && 0 != dual_ent )
216 {
217 dual_ents.insert( dual_ent );
218 continue;
219 }
220 else if( MB_SUCCESS != tmp_result )
221 continue;
222
223 tmp_result = construct_dual_vertex( *rit, dual_ent, false, true );
224 if( MB_SUCCESS != tmp_result ) continue;
225
226 // save it in the list of new dual ents
227 dual_ents.insert( dual_ent );
228 }
229
230 return result;
231 }
232
233 ErrorCode DualTool::construct_dual_vertex( EntityHandle entity,
234 EntityHandle& dual_ent,
235 const bool extra,
236 const bool add_graphics_pt )
237 {
238 // no dual entity; construct one; first need the avg coordinates
239 unsigned int is_dual = 0x1;
240 double avg_pos[3];
241 ErrorCode result = MeshTopoUtil( mbImpl ).get_average_position( entity, avg_pos );
242 if( MB_SUCCESS != result ) return result;
243
244 // now construct the new dual entity
245 result = mbImpl->create_vertex( avg_pos, dual_ent );
246 if( MB_SUCCESS != result ) return result;
247
248 // tag it indicating it's a dual entity
249 result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
250 if( MB_SUCCESS != result ) return result;
251
252 // tag the primal entity with its dual entity and vica versa
253 if( extra )
254 result = mbImpl->tag_set_data( extraDualEntity_tag(), &( entity ), 1, &dual_ent );
255 else
256 result = mbImpl->tag_set_data( dualEntity_tag(), &( entity ), 1, &dual_ent );
257 if( MB_SUCCESS != result ) return result;
258
259 result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( entity ) );
260 if( MB_SUCCESS != result ) return result;
261
262 if( add_graphics_pt )
263 // put a graphics point on that vertex too
264 result = add_graphics_point( dual_ent, avg_pos );
265
266 return result;
267 }
268
269 ErrorCode DualTool::add_graphics_point( EntityHandle entity, double* avg_pos )
270 {
271 // add a graphics pt, placed at the same position as the vertex
272 double my_pos[3];
273 ErrorCode result;
274
275 if( NULL == avg_pos )
276 {
277 result = MeshTopoUtil( mbImpl ).get_average_position( entity, my_pos );
278 if( MB_SUCCESS != result ) return result;
279 }
280 else
281 for( int i = 0; i < 3; i++ )
282 my_pos[i] = avg_pos[i];
283
284 DualTool::GraphicsPoint dum_pt( my_pos, -1 );
285 result = mbImpl->tag_set_data( dualGraphicsPoint_tag(), &entity, 1, &dum_pt );
286 return result;
287 }
288
289 ErrorCode DualTool::construct_dual_edges( const Range& all_faces, Range& dual_ents )
290 {
291 if( all_faces.empty() ) return MB_SUCCESS;
292
293 // make sure they're all faces
294 assert( 2 == CN::Dimension( TYPE_FROM_HANDLE( *all_faces.begin() ) ) &&
295 2 == CN::Dimension( TYPE_FROM_HANDLE( *all_faces.rbegin() ) ) );
296
297 Range::const_iterator rit;
298 EntityHandle dual_ent;
299 unsigned int is_dual = 0x1;
300 ErrorCode tmp_result = MB_SUCCESS;
301 ErrorCode result = MB_SUCCESS;
302
303 for( rit = all_faces.begin(); rit != all_faces.end(); ++rit )
304 {
305 if( tmp_result != MB_SUCCESS ) result = tmp_result;
306
307 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
308 if( MB_SUCCESS == tmp_result && 0 != dual_ent )
309 {
310 dual_ents.insert( dual_ent );
311 continue;
312 }
313
314 // no dual entity; construct one; get the bounding regions
315 std::vector< EntityHandle > out_ents;
316 tmp_result = mbImpl->get_adjacencies( &( *rit ), 1, 3, false, out_ents );
317 if( MB_SUCCESS != tmp_result || out_ents.empty() ) continue;
318
319 // get the dual vertices
320 std::vector< EntityHandle > dual_verts( out_ents.size() );
321 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &out_ents[0], out_ents.size(), &dual_verts[0] );
322 if( MB_SUCCESS != tmp_result ) continue;
323 assert( dual_verts.size() <= 2 );
324
325 double avg_pos[3];
326 bool bdy_face = ( dual_verts.size() == 1 ? true : false );
327 if( bdy_face )
328 {
329 // boundary face - make a dual vertex at the face center and put in list
330 tmp_result = construct_dual_vertex( *rit, dual_ent, true, true );
331
332 // put it on vertex list
333 dual_verts.push_back( dual_ent );
334 }
335
336 assert( dual_verts.size() == 2 );
337
338 // now create the dual edge
339 tmp_result = mbImpl->create_element( MBEDGE, &dual_verts[0], 2, dual_ent );
340 if( MB_SUCCESS != tmp_result || 0 == dual_ent ) continue;
341
342 // save it in the list of new dual ents
343 dual_ents.insert( dual_ent );
344
345 // tag the primal entity with its dual entity and vica versa
346 tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
347 if( MB_SUCCESS != tmp_result ) continue;
348
349 tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) );
350 if( MB_SUCCESS != tmp_result ) continue;
351
352 // tag the edge indicating it's a dual entity
353 tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
354 if( MB_SUCCESS != tmp_result ) continue;
355
356 // add a graphics point to the edge; position depends on whether it's a
357 // bdy face (mid-pt of dual edge) or not (mid-pt of primal face)
358 if( bdy_face )
359 tmp_result = add_graphics_point( dual_ent );
360 else
361 {
362 // get the face's position
363 tmp_result = MeshTopoUtil( mbImpl ).get_average_position( *rit, avg_pos );
364 if( MB_SUCCESS != tmp_result ) continue;
365 tmp_result = add_graphics_point( dual_ent, avg_pos );
366 }
367 if( MB_SUCCESS != tmp_result ) continue;
368 }
369
370 return result;
371 }
372
373 ErrorCode DualTool::construct_dual_faces( const Range& all_edges, Range& dual_ents )
374 {
375 if( all_edges.empty() ) return MB_SUCCESS;
376
377 // make sure they're all edges
378 assert( 1 == CN::Dimension( TYPE_FROM_HANDLE( *all_edges.begin() ) ) &&
379 1 == CN::Dimension( TYPE_FROM_HANDLE( *all_edges.rbegin() ) ) );
380
381 Range::const_iterator rit;
382 EntityHandle dual_ent;
383 unsigned int is_dual = 0x1;
384 ErrorCode tmp_result = MB_SUCCESS;
385 ErrorCode result = MB_SUCCESS;
386 Range equiv_edges;
387 #define TRC \
388 if( MB_SUCCESS != tmp_result ) \
389 { \
390 result = tmp_result; \
391 continue; \
392 }
393 for( rit = all_edges.begin(); rit != all_edges.end(); ++rit )
394 {
395
396 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
397 if( MB_SUCCESS == tmp_result && 0 != dual_ent )
398 {
399 dual_ents.insert( dual_ent );
400 continue;
401 }
402
403 // no dual entity; construct one; get the dual vertices bounding the edge in radial order,
404 // then construct the dual face
405 std::vector< EntityHandle > rad_dverts;
406 bool bdy_edge;
407 tmp_result = get_radial_dverts( *rit, rad_dverts, bdy_edge );
408 TRC if( rad_dverts.empty() ) continue;
409
410 tmp_result = mbImpl->create_element( MBPOLYGON, &rad_dverts[0], rad_dverts.size(), dual_ent );
411 TRC
412
413 // tag it indicating it's a dual entity, and tag primal/dual with dual/primal
414 tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
415 TRC tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
416 TRC tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) );
417 TRC
418
419 // save it in the list of new dual ents
420 dual_ents.insert( dual_ent );
421
422 // add a graphics point to the cell; position depends on whether it's a
423 // bdy cell (mid-pt of cell's vertices) or not (mid-pt of primal edge)
424 double avg_pos[3];
425 tmp_result = MeshTopoUtil( mbImpl ).get_average_position( *rit, avg_pos );
426 TRC if( bdy_edge )
427 {
428
429 // add a new dual edge betw last 2 verts
430 EntityHandle new_edge;
431 tmp_result = mbImpl->create_element( MBEDGE, &rad_dverts[rad_dverts.size() - 2], 2, new_edge );
432 TRC tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &new_edge, 1, &is_dual );
433 TRC
434
435 // tag the new dual edge with the primal edge as it's dual entity; primal
436 // edge IS NOT likewise tagged, since it's already tagged with the 2cell
437 tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &new_edge, 1, &( *rit ) );
438 TRC
439
440 // add a graphics pt, position is center of primal edge
441 tmp_result = add_graphics_point( dual_ent );
442 TRC tmp_result = add_graphics_point( new_edge, avg_pos );
443 TRC
444 }
445
446 else
447 {
448 // if inside, point goes on the 2cell, at primal edge mid-pt
449 tmp_result = add_graphics_point( dual_ent, avg_pos );
450 TRC
451 }
452
453 // check to see whether we have equiv entities; if we find any, save for later fixup
454 Range dum_edges, dum_poly( dual_ent, dual_ent );
455 tmp_result = mbImpl->get_adjacencies( dum_poly, 1, false, dum_edges );
456 if( MB_MULTIPLE_ENTITIES_FOUND == tmp_result )
457 {
458 // we do - need to add adjacencies to disambiguate; use the primal
459 equiv_edges.merge( dum_edges );
460 }
461 }
462
463 if( !equiv_edges.empty() ) result = check_dual_equiv_edges( equiv_edges );
464
465 return result;
466 }
467
468 ErrorCode DualTool::check_dual_equiv_edges( Range& dual_edges )
469 {
470 // fix equivalent dual edges (i.e. edges whose vertices define multiple edges)
471 // by explicitly adding adjacencies to containing polygons; adjacent polygons
472 // found by going through primal
473 ErrorCode tmp_result, result = MB_SUCCESS;
474
475 Range all_dedges( dual_edges );
476 // first, go through all dual edges and find equivalent edges (by looking for
477 // up-adjacent edges on the vertices of each edge)
478 for( Range::iterator rit = dual_edges.begin(); rit != dual_edges.end(); ++rit )
479 {
480 Range connect, dum_range( *rit, *rit );
481 tmp_result = mbImpl->get_adjacencies( dum_range, 0, false, connect );
482 if( MB_SUCCESS != tmp_result ) continue;
483 tmp_result = mbImpl->get_adjacencies( connect, 1, false, all_dedges, Interface::UNION );
484 if( MB_SUCCESS != tmp_result ) continue;
485 }
486
487 // save a copy for checking later
488 Range save_all_2cells;
489
490 // go through each edge
491 while( !all_dedges.empty() )
492 {
493 EntityHandle this_edge = *all_dedges.begin();
494 all_dedges.erase( all_dedges.begin() );
495
496 const EntityHandle* connect;
497 int num_connect;
498 result = mbImpl->get_connectivity( this_edge, connect, num_connect );
499 if( MB_SUCCESS != result ) continue;
500
501 Range dum_edges, verts;
502 verts.insert( connect[0] );
503 verts.insert( connect[1] );
504 tmp_result = mbImpl->get_adjacencies( verts, 1, false, dum_edges );
505 if( MB_SUCCESS != tmp_result )
506 {
507 result = tmp_result;
508 continue;
509 }
510 if( dum_edges.size() == 1 )
511 {
512 // not an equiv edge - already removed from list, so just continue
513 continue;
514 }
515
516 // ok, have an equiv entity - fix by looking through primal
517 // pre-get the primal of these
518 EntityHandle dedge_quad;
519 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &this_edge, 1, &dedge_quad );
520 if( MB_SUCCESS != tmp_result )
521 {
522 result = tmp_result;
523 continue;
524 }
525
526 if( MBQUAD == mbImpl->type_from_handle( dedge_quad ) )
527 {
528
529 // get the primal edges adj to quad
530 Range dum_quad_range( dedge_quad, dedge_quad ), adj_pedges;
531 tmp_result = mbImpl->get_adjacencies( dum_quad_range, 1, false, adj_pedges );
532 if( MB_SUCCESS != tmp_result )
533 {
534 result = tmp_result;
535 continue;
536 }
537 // get the dual 2cells corresponding to those pedges
538 std::vector< EntityHandle > dcells;
539 dcells.resize( adj_pedges.size() );
540 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), adj_pedges, &dcells[0] );
541 if( MB_SUCCESS != tmp_result )
542 {
543 result = tmp_result;
544 continue;
545 }
546 // now add explicit adjacencies from the dedge to those dcells
547 std::vector< EntityHandle >::iterator vit;
548 for( vit = dcells.begin(); vit != dcells.end(); ++vit )
549 {
550 save_all_2cells.insert( *vit );
551
552 assert( MBPOLYGON == mbImpl->type_from_handle( *vit ) );
553 tmp_result = mbImpl->add_adjacencies( this_edge, &( *vit ), 1, false );
554 if( MB_SUCCESS != tmp_result )
555 {
556 result = tmp_result;
557 continue;
558 }
559 // check that there are really adjacencies and *vit is in them
560 const EntityHandle* adjs;
561 int num_adjs;
562 tmp_result = reinterpret_cast< Core* >( mbImpl )->a_entity_factory()->get_adjacencies( this_edge, adjs,
563 num_adjs );
564 if( NULL == adjs || std::find( adjs, adjs + num_adjs, *vit ) == adjs + num_adjs )
565 std::cout << "Add_adjacencies failed in construct_dual_faces." << std::endl;
566 }
567 }
568 else
569 {
570 // else, have a dual edge representing a bdy edge - tie directly to
571 // dual entity if its dual entity
572 EntityHandle bdy_dcell;
573 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &dedge_quad, 1, &bdy_dcell );
574 TRC assert( MBPOLYGON == mbImpl->type_from_handle( bdy_dcell ) );
575
576 tmp_result = mbImpl->add_adjacencies( this_edge, &bdy_dcell, 1, false );
577 if( MB_SUCCESS != tmp_result )
578 {
579 result = tmp_result;
580 continue;
581 }
582 }
583 }
584
585 // sanity check - look for adj edges again, and check for equiv entities
586 for( Range::iterator vit = save_all_2cells.begin(); vit != save_all_2cells.end(); ++vit )
587 {
588 Range adj_edges, dum_quad_range;
589 dum_quad_range.insert( *vit );
590 assert( MBPOLYGON == mbImpl->type_from_handle( *vit ) );
591 tmp_result = mbImpl->get_adjacencies( dum_quad_range, 1, false, adj_edges );
592 if( MB_MULTIPLE_ENTITIES_FOUND == tmp_result )
593 {
594 std::cout << "Multiple entities returned for polygon " << mbImpl->id_from_handle( *vit ) << "."
595 << std::endl;
596 continue;
597 }
598 }
599 // success!
600 return result;
601 }
602
603 ErrorCode DualTool::construct_dual_cells( const Range& all_verts, Range& dual_ents )
604 {
605 if( all_verts.empty() ) return MB_SUCCESS;
606
607 // make sure they're all edges
608 assert( 0 == CN::Dimension( TYPE_FROM_HANDLE( *all_verts.begin() ) ) &&
609 0 == CN::Dimension( TYPE_FROM_HANDLE( *all_verts.rbegin() ) ) );
610
611 Range::const_iterator rit;
612 EntityHandle dual_ent;
613 unsigned int is_dual = 0x1;
614 ErrorCode tmp_result = MB_SUCCESS;
615 ErrorCode result = MB_SUCCESS;
616 std::vector< EntityHandle > edges, dfaces;
617
618 for( rit = all_verts.begin(); rit != all_verts.end(); ++rit )
619 {
620 if( tmp_result != MB_SUCCESS ) result = tmp_result;
621
622 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
623 if( MB_SUCCESS == tmp_result && 0 != dual_ent )
624 {
625 dual_ents.insert( dual_ent );
626 continue;
627 }
628
629 // no dual entity; construct one; get the edges bounding the vertex
630 edges.clear();
631 dfaces.clear();
632 tmp_result = mbImpl->get_adjacencies( &( *rit ), 1, 1, false, edges );
633 if( MB_SUCCESS != tmp_result ) continue;
634
635 // get the dual faces corresponding to the edges
636 dfaces.resize( edges.size() );
637 tmp_result = mbImpl->tag_get_data( dualEntity_tag(), &edges[0], edges.size(), &dfaces[0] );
638 if( MB_SUCCESS != tmp_result ) continue;
639
640 // create the dual cell from those faces
641 tmp_result = mbImpl->create_element( MBPOLYHEDRON, &dfaces[0], dfaces.size(), dual_ent );
642 if( MB_SUCCESS != tmp_result || 0 == dual_ent ) continue;
643
644 // save it in the list of new dual ents
645 dual_ents.insert( dual_ent );
646
647 // tag it indicating it's a dual entity
648 tmp_result = mbImpl->tag_set_data( isDualCell_tag(), &dual_ent, 1, &is_dual );
649 if( MB_SUCCESS != tmp_result ) continue;
650
651 // tag the primal entity with its dual entity and vica versa
652 tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &( *rit ), 1, &dual_ent );
653 if( MB_SUCCESS != tmp_result ) continue;
654 tmp_result = mbImpl->tag_set_data( dualEntity_tag(), &dual_ent, 1, &( *rit ) );
655 if( MB_SUCCESS != tmp_result ) continue;
656 }
657
658 return result;
659 }
660
661 //! given an edge handle, return a list of dual vertices in radial order
662 //! around the edge
663 ErrorCode DualTool::get_radial_dverts( const EntityHandle edge,
664 std::vector< EntityHandle >& rad_dverts,
665 bool& bdy_edge )
666 {
667 rad_dverts.clear();
668
669 std::vector< EntityHandle > rad_faces, rad_ents;
670 ErrorCode result = MeshTopoUtil( mbImpl ).star_entities( edge, rad_faces, bdy_edge, 0, &rad_ents );
671 if( MB_SUCCESS != result ) return result;
672
673 if( bdy_edge )
674 {
675 // if we're a bdy edge, change the order back to what DualTool expects
676 rad_ents.push_back( *rad_faces.rbegin() );
677 rad_ents.push_back( *rad_faces.begin() );
678 }
679
680 rad_dverts.resize( rad_ents.size() );
681 for( unsigned int i = 0; i < rad_ents.size(); i++ )
682 {
683 EntityHandle dual_ent;
684 result = mbImpl->tag_get_data( dualEntity_tag(), &rad_ents[i], 1, &dual_ent );
685 if( !bdy_edge || i < rad_ents.size() - 2 )
686 rad_dverts[i] = dual_ent;
687 else
688 {
689 // fix up this entry
690 assert( mbImpl->type_from_handle( dual_ent ) == MBEDGE );
691
692 // get connectivity of that edge
693 const EntityHandle* connect;
694 int num_connect;
695 result = mbImpl->get_connectivity( dual_ent, connect, num_connect );
696 if( MB_SUCCESS != result ) return result;
697
698 // we want the one that's not already on the list; reuse last_face
699 int last_hex = ( i == rad_ents.size() - 1 ? 0 : i - 1 );
700 EntityHandle last_face = ( connect[0] == rad_dverts[last_hex] ? connect[1] : connect[0] );
701 rad_dverts[i] = last_face;
702 }
703 }
704
705 return result;
706 }
707
708 //! construct the dual entities for a hex mesh, including dual surfaces & curves
709 ErrorCode DualTool::construct_hex_dual( Range& entities )
710 {
711 std::vector< EntityHandle > evec;
712 std::copy( entities.begin(), entities.end(), std::back_inserter( evec ) );
713 return construct_hex_dual( &evec[0], evec.size() );
714 }
715
716 //! construct the dual entities for a hex mesh, including dual surfaces & curves
717 ErrorCode DualTool::construct_hex_dual( EntityHandle* entities, const int num_entities )
718 {
719 // really quite simple:
720
721 // construct the dual...
722 ErrorCode result = construct_dual( entities, num_entities );
723 if( MB_SUCCESS != result )
724 {
725 std::cerr << "Error constructing dual entities for primal entities." << std::endl;
726 return result;
727 }
728
729 // now traverse to build 1d and 2d hyperplanes
730 result = construct_dual_hyperplanes( 1, entities, num_entities );
731 if( MB_SUCCESS != result )
732 {
733 std::cerr << "Problem traversing 1d hyperplanes." << std::endl;
734 return result;
735 }
736
737 result = construct_dual_hyperplanes( 2, entities, num_entities );
738 if( MB_SUCCESS != result )
739 {
740 std::cerr << "Problem traversing 2d hyperplanes." << std::endl;
741 return result;
742 }
743
744 result = construct_hp_parent_child();
745 if( MB_SUCCESS != result )
746 {
747 std::cerr << "Problem constructing parent/child relations between hyperplanes." << std::endl;
748 return result;
749 }
750
751 // see? simple, just like I said
752 return MB_SUCCESS;
753 }
754
755 //! get the cells of the dual
756 ErrorCode DualTool::get_dual_entities( const int dim, EntityHandle* entities, const int num_entities, Range& dual_ents )
757 {
758 if( 0 == isDualCell_tag() ) return MB_SUCCESS;
759 if( 0 > dim || 3 < dim ) return MB_INDEX_OUT_OF_RANGE;
760
761 unsigned int dum = 0x1;
762 const void* dum_ptr = &dum;
763 static EntityType dual_type[] = { MBVERTEX, MBEDGE, MBPOLYGON, MBPOLYHEDRON };
764
765 Range dim_ents;
766
767 ErrorCode result;
768
769 if( 0 == entities || 0 == num_entities )
770 {
771 // just get all the dual entities of this dimension
772 result = mbImpl->get_entities_by_type_and_tag( 0, dual_type[dim], &isDualCellTag, &dum_ptr, 1, dual_ents );
773 }
774 else
775 {
776 // else look for specific dual entities
777 result = mbImpl->get_adjacencies( entities, num_entities, 3 - dim, false, dim_ents, Interface::UNION );
778 if( MB_SUCCESS != result ) return result;
779 std::vector< EntityHandle > dual_ents_vec( dim_ents.size() );
780 result = mbImpl->tag_get_data( dualEntity_tag(), dim_ents, &dual_ents_vec[0] );
781 if( MB_SUCCESS != result ) return result;
782 std::copy( dual_ents_vec.begin(), dual_ents_vec.end(), range_inserter( dual_ents ) );
783 }
784
785 return result;
786 }
787
788 //! get the faces of the dual
789 ErrorCode DualTool::get_dual_entities( const int dim,
790 EntityHandle* entities,
791 const int num_entities,
792 std::vector< EntityHandle >& dual_ents )
793 {
794 Range tmp_range;
795 ErrorCode result = get_dual_entities( dim, entities, num_entities, tmp_range );
796 if( MB_SUCCESS != result ) return result;
797
798 // dual_ents.insert(dual_ents.end(), tmp_range.begin(), tmp_range.end());
799 dual_ents.reserve( dual_ents.size() + tmp_range.size() );
800 for( Range::const_iterator it = tmp_range.begin(); it != tmp_range.end(); ++it )
801 {
802 dual_ents.push_back( *it );
803 }
804 return MB_SUCCESS;
805 }
806
807 ErrorCode DualTool::get_dual_hyperplanes( const Interface* impl, const int dim, Range& dual_ents )
808 {
809 if( dim != 1 && dim != 2 ) return MB_INDEX_OUT_OF_RANGE;
810
811 Tag dual_tag;
812 ErrorCode result;
813
814 if( dim == 1 )
815 result = impl->tag_get_handle( DUAL_CURVE_TAG_NAME, 1, MB_TYPE_HANDLE, dual_tag );
816 else
817 result = impl->tag_get_handle( DUAL_SURFACE_TAG_NAME, 1, MB_TYPE_HANDLE, dual_tag );
818
819 if( MB_SUCCESS == result )
820 result = impl->get_entities_by_type_and_tag( 0, MBENTITYSET, &dual_tag, NULL, 1, dual_ents, Interface::UNION );
821
822 return result;
823 }
824
825 ErrorCode DualTool::construct_dual_hyperplanes( const int dim, EntityHandle* entities, const int num_entities )
826 {
827 // this function traverses dual faces of input dimension, constructing
828 // dual hyperplanes of them in sets as it goes
829
830 // check various inputs
831 int num_quads, num_hexes;
832 if(
833 // this function only makes sense for dim == 1 or dim == 2
834 ( dim != 1 && dim != 2 ) ||
835 // should either be quads or hexes around
836 mbImpl->get_number_entities_by_type( 0, MBQUAD, num_quads ) != MB_SUCCESS ||
837 mbImpl->get_number_entities_by_type( 0, MBHEX, num_hexes ) != MB_SUCCESS ||
838 // if we're asking for 1d dual ents, should be quads around
839 ( num_quads == 0 && dim == 1 ) ||
840 // if we're asking for 2d dual ents, should be hexes around
841 ( num_hexes == 0 && dim == 2 ) )
842 return MB_FAILURE;
843
844 Tag hp_tag = ( 1 == dim ? dualCurve_tag() : dualSurface_tag() );
845
846 // two stacks: one completely untreated entities, and the other untreated
847 // entities on the current dual hyperplane
848 std::vector< EntityHandle > tot_untreated;
849
850 // put dual entities of this dimension on the untreated list
851 ErrorCode result = get_dual_entities( dim, entities, num_entities, tot_untreated );
852 if( MB_SUCCESS != result ) return result;
853
854 // main part of traversal loop
855 EntityHandle this_ent;
856 EntityHandle this_hp;
857
858 while( !tot_untreated.empty() )
859 {
860 if( debug && dim == 2 /*(tot_untreated.size()%report == 0)*/ )
861 std::cout << "Untreated list size " << tot_untreated.size() << "." << std::endl;
862
863 this_ent = tot_untreated.back();
864 tot_untreated.pop_back();
865 result = mbImpl->tag_get_data( hp_tag, &this_ent, 1, &this_hp );
866 if( MB_SUCCESS != result && MB_TAG_NOT_FOUND != result ) return result;
867
868 // d for this entity having a hyperplane assignment already
869 else if( this_hp != 0 )
870 continue;
871
872 if( 1 == dim && check_1d_loop_edge( this_ent ) ) continue;
873
874 // inner loop: traverse the hyperplane 'till we don't have any more
875 result = traverse_hyperplane( hp_tag, this_hp, this_ent );
876 if( MB_SUCCESS != result )
877 {
878 std::cout << "Failed to traverse hyperplane ";
879 if( this_hp )
880 std::cout << mbImpl->id_from_handle( this_hp ) << "." << std::endl;
881 else
882 std::cout << "0." << std::endl;
883 return result;
884 }
885
886 // ok, now order the edges if it's a chord
887 if( 1 == dim ) order_chord( this_hp );
888 }
889
890 return MB_SUCCESS;
891 }
892
893 ErrorCode DualTool::traverse_hyperplane( const Tag hp_tag, EntityHandle& this_hp, EntityHandle this_ent )
894 {
895 Range tmp_star, star, tmp_range, new_hyperplane_ents;
896 std::vector< EntityHandle > hp_untreated;
897 int dim = mbImpl->dimension_from_handle( this_ent );
898 MeshTopoUtil mtu( mbImpl );
899 this_hp = 0;
900 ErrorCode result;
901
902 unsigned short mark_val = 0x0;
903 Tag mark_tag;
904 result = mbImpl->tag_get_handle( "__hyperplane_mark", 1, MB_TYPE_BIT, mark_tag, MB_TAG_CREAT | MB_TAG_BIT );
905 if( MB_SUCCESS != result ) return result;
906 mark_val = 0x1;
907
908 while( 0 != this_ent )
909 {
910 EntityHandle tmp_hp = get_dual_hyperplane( this_ent );
911 if( 0 == this_hp && 0 != tmp_hp ) this_hp = tmp_hp;
912
913 if( 0 == tmp_hp ) new_hyperplane_ents.insert( this_ent );
914
915 if( debug && hp_untreated.size() % 10 == 0 )
916 std::cout << "Dual surface " << this_hp << ", hp_untreated list size = " << hp_untreated.size() << "."
917 << std::endl;
918
919 // get the 2nd order adjacencies through lower dimension
920 tmp_range.clear();
921 tmp_star.clear();
922 star.clear();
923 result = mtu.get_bridge_adjacencies( this_ent, dim - 1, dim, star );RR;
924
925 // get the bridge adjacencies through higher dimension
926 result = mtu.get_bridge_adjacencies( this_ent, dim + 1, dim, tmp_star );RR;
927 tmp_range = subtract( star, tmp_star );
928
929 for( Range::iterator rit = tmp_range.begin(); rit != tmp_range.end(); ++rit )
930 {
931 if( new_hyperplane_ents.find( *rit ) != new_hyperplane_ents.end() ) continue;
932
933 // check for tag first, 'cuz it's probably faster than checking adjacencies
934 // assign to avoid valgrind warning
935 unsigned short tmp_mark = 0x0;
936 result = mbImpl->tag_get_data( mark_tag, &( *rit ), 1, &tmp_mark );
937 if( MB_SUCCESS == result && mark_val == tmp_mark ) continue;
938
939 // if it's on the loop, it's not eligible
940 if( 1 == dim && check_1d_loop_edge( *rit ) ) continue;
941
942 // have one on this hp; just put it on the hp_untreated list for now,
943 // will get tagged and put in the hp set later
944 hp_untreated.push_back( *rit );
945 result = mbImpl->tag_set_data( mark_tag, &( *rit ), 1, &mark_val );
946 if( MB_SUCCESS != result ) return result;
947 }
948
949 // end of inner loop; get the next this_ent, or set to zero
950 if( hp_untreated.empty() )
951 this_ent = 0;
952 else
953 {
954 this_ent = hp_untreated.back();
955 hp_untreated.pop_back();
956 }
957 }
958
959 if( debug_ap )
960 {
961 std::string hp_name;
962 if( 2 == dim )
963 hp_name = "sheet";
964 else
965 hp_name = "chord";
966
967 if( 0 == this_hp )
968 std::cout << "Constructed new " << hp_name << " with ";
969 else
970 {
971 int this_id;
972 result = mbImpl->tag_get_data( globalId_tag(), &this_hp, 1, &this_id );RR;
973 std::cout << "Added to " << hp_name << " " << this_id << " ";
974 }
975 if( dim == 2 )
976 std::cout << "edges:" << std::endl;
977 else
978 std::cout << "quads:" << std::endl;
979 std::vector< EntityHandle > pents( new_hyperplane_ents.size() );
980 result = mbImpl->tag_get_data( dualEntity_tag(), new_hyperplane_ents, &pents[0] );RR;
981 for( std::vector< EntityHandle >::iterator vit = pents.begin(); vit != pents.end(); ++vit )
982 {
983 if( vit != pents.begin() ) std::cout << ", ";
984 std::cout << mbImpl->id_from_handle( *vit );
985 }
986 std::cout << std::endl;
987 }
988
989 if( 0 == this_hp )
990 {
991 // ok, doesn't have one; make a new hyperplane
992 int new_id = -1;
993 result = construct_new_hyperplane( dim, this_hp, new_id );
994 if( MB_SUCCESS != result ) return result;
995
996 if( debug_ap )
997 {
998 std::cout << "New ";
999 if( 2 == dim )
1000 std::cout << " sheet ";
1001 else
1002 std::cout << " chord ";
1003 std::cout << new_id << " constructed." << std::endl;
1004 }
1005 }
1006
1007 // set the hp_val for entities which didn't have one before
1008 std::vector< EntityHandle > hp_tags( new_hyperplane_ents.size() );
1009 std::fill( hp_tags.begin(), hp_tags.end(), this_hp );
1010 result = mbImpl->tag_set_data( hp_tag, new_hyperplane_ents, &hp_tags[0] );
1011 if( MB_SUCCESS != result ) return result;
1012 result = mbImpl->add_entities( this_hp, new_hyperplane_ents );
1013 if( MB_SUCCESS != result ) return result;
1014
1015 // unmark the entities by removing the tag
1016 result = mbImpl->tag_delete( mark_tag );
1017 if( MB_SUCCESS != result ) return result;
1018
1019 return MB_SUCCESS;
1020 }
1021
1022 ErrorCode DualTool::order_chord( EntityHandle chord_set )
1023 {
1024 // re-order the 1cells in the set so they are in order along the chord
1025 // start by finding the vertex dual to a quad
1026 Range verts, one_cells;
1027 ErrorCode result = mbImpl->get_entities_by_dimension( chord_set, 1, one_cells );
1028 if( MB_SUCCESS != result || one_cells.empty() ) return MB_FAILURE;
1029
1030 result = mbImpl->get_adjacencies( one_cells, 0, false, verts, Interface::UNION );
1031 if( MB_SUCCESS != result || verts.empty() ) return MB_FAILURE;
1032
1033 EntityHandle last_vert = 0;
1034 for( Range::iterator rit = verts.begin(); rit != verts.end(); ++rit )
1035 {
1036 if( TYPE_FROM_HANDLE( get_dual_entity( *rit ) ) == MBQUAD )
1037 {
1038 last_vert = *rit;
1039 break;
1040 }
1041 }
1042 // if there's no vertex owned by a quad, just start with 1st one
1043 if( 0 == last_vert ) last_vert = *verts.begin();
1044
1045 // now, skip from vertex to vertex, building a list of 1cells
1046 std::vector< EntityHandle > ordered_1cells;
1047 EntityHandle last_1cell = 0;
1048 Range dum1, dum2;
1049 const EntityHandle* connect;
1050 int num_connect;
1051 ErrorCode tmp_result = MB_SUCCESS;
1052 while( ordered_1cells.size() != one_cells.size() )
1053 {
1054 dum1 = one_cells;
1055 result = mbImpl->get_adjacencies( &last_vert, 1, 1, false, dum1 );
1056 if( 0 != last_1cell ) dum1.erase( last_1cell );
1057 // assert(1 == dum1.size());
1058 if( 0 != last_1cell && 1 != dum1.size() )
1059 {
1060 std::cerr << "unexpected size traversing chord." << std::endl;
1061 tmp_result = MB_FAILURE;
1062 }
1063
1064 last_1cell = *dum1.begin();
1065 ordered_1cells.push_back( last_1cell );
1066 result = mbImpl->get_connectivity( last_1cell, connect, num_connect );RR;
1067 if( last_vert == connect[0] )
1068 last_vert = connect[1];
1069 else
1070 last_vert = connect[0];
1071 }
1072
1073 // now have the 1cells in order, replace them in the set
1074 if( MB_SUCCESS == tmp_result )
1075 {
1076 result = mbImpl->remove_entities( chord_set, one_cells );RR;
1077 result = mbImpl->add_entities( chord_set, &ordered_1cells[0], ordered_1cells.size() );RR;
1078 }
1079
1080 return MB_SUCCESS;
1081 }
1082
1083 ErrorCode DualTool::construct_new_hyperplane( const int dim, EntityHandle& new_hyperplane, int& id )
1084 {
1085 ErrorCode result;
1086 if( 1 == dim )
1087 result = mbImpl->create_meshset( ( MESHSET_ORDERED | MESHSET_TRACK_OWNER ), new_hyperplane );
1088 else
1089 result = mbImpl->create_meshset( ( MESHSET_SET | MESHSET_TRACK_OWNER ), new_hyperplane );
1090 if( MB_SUCCESS != result ) return result;
1091
1092 if( -1 == id )
1093 {
1094 Range all_hyperplanes;
1095 result = get_dual_hyperplanes( mbImpl, dim, all_hyperplanes );RR;
1096 std::vector< int > gids( all_hyperplanes.size() );
1097 result = mbImpl->tag_get_data( globalIdTag, all_hyperplanes, ( gids.empty() ) ? NULL : &gids[0] );RR;
1098 for( unsigned int i = 0; i < gids.size(); i++ )
1099 if( gids[i] > id ) id = gids[i];
1100 id++;
1101 if( 0 == id ) id++;
1102 }
1103
1104 result = mbImpl->tag_set_data( globalId_tag(), &new_hyperplane, 1, &id );RR;
1105 Tag hp_tag = ( 1 == dim ? dualCurve_tag() : dualSurface_tag() );
1106 result = mbImpl->tag_set_data( hp_tag, &new_hyperplane, 1, &new_hyperplane );RR;
1107
1108 // assign a category name to these sets
1109 static const char dual_category_names[2][CATEGORY_TAG_SIZE] = { "Chord\0", "Sheet\0" };
1110
1111 result = mbImpl->tag_set_data( categoryTag, &new_hyperplane, 1, dual_category_names[dim - 1] );
1112
1113 return result;
1114 }
1115
1116 bool DualTool::check_1d_loop_edge( EntityHandle this_ent )
1117 {
1118 // make sure it's an edge
1119 if( MBEDGE != mbImpl->type_from_handle( this_ent ) ) return false;
1120
1121 // also has to be a dual entity
1122 unsigned int dum;
1123 ErrorCode result = mbImpl->tag_get_data( isDualCell_tag(), &this_ent, 1, &dum );
1124 if( MB_SUCCESS != result || dum != 0x1 ) return false;
1125
1126 const EntityHandle* verts;
1127 EntityHandle vert_tags[2];
1128 int num_verts;
1129 result = mbImpl->get_connectivity( this_ent, verts, num_verts );
1130 if( MB_SUCCESS != result ) return false;
1131
1132 result = mbImpl->tag_get_data( dualEntity_tag(), verts, 2, vert_tags );
1133 if( MB_SUCCESS != result || mbImpl->type_from_handle( vert_tags[0] ) != MBQUAD ||
1134 mbImpl->type_from_handle( vert_tags[1] ) != MBQUAD )
1135 return false;
1136
1137 else
1138 return true;
1139 }
1140
1141 ErrorCode DualTool::construct_hp_parent_child()
1142 {
1143 Range dual_surfs, dual_cells, dual_edges;
1144 ErrorCode result = this->get_dual_hyperplanes( mbImpl, 2, dual_surfs );
1145 if( MB_SUCCESS != result || dual_surfs.empty() ) return result;
1146 std::vector< EntityHandle > dual_curve_sets;
1147
1148 for( Range::iterator surf_it = dual_surfs.begin(); surf_it != dual_surfs.end(); ++surf_it )
1149 {
1150 // get all the cells, edges in those cells, and chords for those edges
1151 dual_cells.clear();
1152 result = mbImpl->get_entities_by_handle( *surf_it, dual_cells );
1153 if( MB_SUCCESS != result ) return result;
1154 dual_edges.clear();
1155 result = mbImpl->get_adjacencies( dual_cells, 1, false, dual_edges, Interface::UNION );
1156 if( MB_SUCCESS != result ) return result;
1157 dual_curve_sets.resize( dual_edges.size() );
1158 result = mbImpl->tag_get_data( dualCurve_tag(), dual_edges, &dual_curve_sets[0] );
1159 if( MB_SUCCESS != result ) return result;
1160
1161 // reuse dual_cells to get unique list of chord sets
1162 dual_cells.clear();
1163 for( unsigned int i = 0; i < dual_edges.size(); i++ )
1164 if( dual_curve_sets[i] != 0 ) dual_cells.insert( dual_curve_sets[i] );
1165
1166 // now connect up this dual surf with all the 1d ones
1167 for( Range::iterator rit = dual_cells.begin(); rit != dual_cells.end(); ++rit )
1168 {
1169 result = mbImpl->add_parent_child( *surf_it, *rit );
1170 if( MB_SUCCESS != result ) return result;
1171 }
1172 }
1173
1174 return MB_SUCCESS;
1175 }
1176
1177 ErrorCode DualTool::get_graphics_points( EntityHandle dual_ent,
1178 std::vector< int >& npts,
1179 std::vector< GraphicsPoint >& points )
1180 {
1181 // shouldn't be a set
1182 assert( MBENTITYSET != mbImpl->type_from_handle( dual_ent ) );
1183
1184 // get the graphics points comprising the given entity
1185 GraphicsPoint gp_array[DualTool::GP_SIZE];
1186
1187 ErrorCode result = MB_SUCCESS;
1188
1189 // switch based on topological dimension
1190 switch( mbImpl->dimension_from_handle( dual_ent ) )
1191 {
1192 case 0:
1193 // just return the vertex point
1194 result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, gp_array );
1195 if( MB_SUCCESS == result ) points.push_back( gp_array[0] );
1196
1197 break;
1198
1199 case 1:
1200 // get my graphics point then those of my vertices
1201 const EntityHandle* connect;
1202 int num_connect;
1203 result = mbImpl->get_connectivity( dual_ent, connect, num_connect );
1204 if( MB_SUCCESS != result ) break;
1205
1206 result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), connect, 2, gp_array );
1207 if( MB_SUCCESS == result )
1208 {
1209 points.push_back( gp_array[0] );
1210 points.push_back( gp_array[0] );
1211 points.push_back( gp_array[1] );
1212 result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, gp_array );
1213 if( MB_SUCCESS == result ) points[1] = gp_array[0];
1214 }
1215
1216 npts.push_back( 3 );
1217
1218 break;
1219
1220 case 2:
1221 result = get_cell_points( dual_ent, npts, points );
1222 break;
1223 }
1224
1225 return result;
1226 }
1227
1228 ErrorCode DualTool::get_cell_points( EntityHandle dual_ent,
1229 std::vector< int >& npts,
1230 std::vector< GraphicsPoint >& points )
1231 {
1232 assert( MBPOLYGON == mbImpl->type_from_handle( dual_ent ) );
1233
1234 // get the 1cells in this 2cell
1235 Range one_cells;
1236
1237 Range tc_range;
1238 tc_range.insert( dual_ent );
1239 ErrorCode result = mbImpl->get_adjacencies( tc_range, 1, false, one_cells, Interface::UNION );RR;
1240
1241 int num_edges = one_cells.size();
1242 std::vector< GraphicsPoint > dum_gps( num_edges + 1 );
1243
1244 // get graphics points for 0cells and for this cell
1245 result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), one_cells, &dum_gps[0] );RR;
1246 result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), &dual_ent, 1, &( dum_gps[num_edges] ) );RR;
1247
1248 Range::iterator eit;
1249 const EntityHandle* connect;
1250 int num_connect;
1251 GraphicsPoint vert_gps[2];
1252 int i;
1253 for( i = 0, eit = one_cells.begin(); i < num_edges; i++, ++eit )
1254 {
1255 // get the vertices and the graphics points for them
1256 result = mbImpl->get_connectivity( *eit, connect, num_connect );RR;
1257 result = mbImpl->tag_get_data( dualGraphicsPoint_tag(), connect, 2, vert_gps );RR;
1258
1259 // make the 2 tris corresponding to this edge; don't worry about order
1260 // for now
1261 npts.push_back( 3 );
1262 points.push_back( dum_gps[num_edges] );
1263 points.push_back( vert_gps[0] );
1264 points.push_back( dum_gps[i] );
1265
1266 npts.push_back( 3 );
1267 points.push_back( dum_gps[num_edges] );
1268 points.push_back( dum_gps[i] );
1269 points.push_back( vert_gps[1] );
1270 }
1271
1272 return result;
1273 }
1274
1275 ErrorCode DualTool::get_graphics_points( const Range& in_range,
1276 std::vector< GraphicsPoint >& points,
1277 const bool assign_ids,
1278 const int start_id )
1279 {
1280 // return graphics points on dual entities in in_range or in entities
1281 // in sets in in_range
1282 ErrorCode result;
1283
1284 // for each dual hyperplane set:
1285 Range::const_iterator rit;
1286
1287 Range two_cells, all_cells;
1288 for( rit = in_range.begin(); rit != in_range.end(); ++rit )
1289 {
1290 // for each entity:
1291 two_cells.clear();
1292 EntityType this_type = mbImpl->type_from_handle( *rit );
1293 if( MBENTITYSET == this_type )
1294 {
1295 result = mbImpl->get_entities_by_handle( *rit, two_cells );RR;
1296
1297 std::copy( two_cells.begin(), two_cells.end(), range_inserter( all_cells ) );
1298 }
1299
1300 else
1301 {
1302 two_cells.insert( *rit );
1303 assert( this_type == MBVERTEX || this_type == MBEDGE || this_type == MBPOLYGON ||
1304 this_type == MBPOLYHEDRON );
1305 }
1306
1307 result = mbImpl->get_adjacencies( two_cells, 0, false, all_cells, Interface::UNION );RR;
1308 result = mbImpl->get_adjacencies( two_cells, 1, false, all_cells, Interface::UNION );RR;
1309 }
1310
1311 // get graphics points
1312 points.resize( all_cells.size() );
1313
1314 result = mbImpl->tag_get_data( dualGraphicsPointTag, all_cells, &points[0] );RR;
1315
1316 if( assign_ids )
1317 {
1318 int i = start_id;
1319
1320 for( std::vector< GraphicsPoint >::iterator vit = points.begin(); vit != points.end(); ++vit )
1321 vit->id = i++;
1322
1323 result = mbImpl->tag_set_data( dualGraphicsPoint_tag(), all_cells, &points[0] );RR;
1324 }
1325
1326 return result;
1327 }
1328
1329 EntityHandle DualTool::next_loop_vertex( const EntityHandle last_v,
1330 const EntityHandle this_v,
1331 const EntityHandle dual_surf )
1332 {
1333 // given two vertices, find the next one on the loop; if one is a dual
1334 // surface, then just choose either one for that surface
1335 assert( ( 0 == last_v || mbImpl->type_from_handle( last_v ) == MBVERTEX ) &&
1336 mbImpl->type_from_handle( this_v ) == MBVERTEX && mbImpl->type_from_handle( dual_surf ) == MBENTITYSET );
1337
1338 // get the connected vertices
1339 MeshTopoUtil tpu( mbImpl );
1340 Range other_verts;
1341 ErrorCode result = tpu.get_bridge_adjacencies( this_v, 1, 0, other_verts );
1342 if( MB_SUCCESS != result || other_verts.empty() ) return 0;
1343
1344 // if (mbImpl->type_from_handle(last_v) == MBENTITYSET) {
1345 // dual surface, choose either; first get a 2cell on this surface
1346 Range tcells, tcells2, verts;
1347 result = mbImpl->get_entities_by_type( dual_surf, MBPOLYGON, tcells );
1348 if( MB_SUCCESS != result || tcells.empty() ) return 0;
1349
1350 // ok, pay attention here: first get 2cells common to dual surface and this_v
1351 verts.insert( this_v );
1352 result = mbImpl->get_adjacencies( verts, 2, false, tcells );
1353 if( MB_SUCCESS != result || tcells.empty() ) return 0;
1354
1355 // next get vertices common to both 2cells and subtract from other_verts; also
1356 // remove last_v if it's non-zero
1357 verts.clear();
1358 result = mbImpl->get_adjacencies( tcells, 0, false, verts );
1359 if( MB_SUCCESS != result || verts.empty() ) return 0;
1360
1361 Range tmp_verts = subtract( other_verts, verts );
1362 other_verts.swap( tmp_verts );
1363 if( 0 != last_v ) other_verts.erase( last_v );
1364
1365 // now get intersection of remaining vertices and 2 2cells vertices
1366 // look at each one successively; should find one, maybe not on both
1367 tmp_verts = other_verts;
1368 Range tmp_faces( *tcells.begin(), *tcells.begin() );
1369 result = mbImpl->get_adjacencies( tmp_faces, 0, false, tmp_verts );
1370 if( MB_SUCCESS == result && !tmp_verts.empty() ) return *tmp_verts.begin();
1371 tmp_faces.clear();
1372 tmp_faces.insert( *tcells.rbegin() );
1373 result = mbImpl->get_adjacencies( tmp_faces, 0, false, other_verts );
1374 if( MB_SUCCESS == result && !other_verts.empty() ) return *other_verts.begin();
1375
1376 // if we got here, there isn't any
1377 return 0;
1378 }
1379
1380 EntityHandle DualTool::get_dual_hyperplane( const EntityHandle ncell )
1381 {
1382 // get the sheet or chord it's in
1383 std::vector< EntityHandle > adj_sets;
1384 ErrorCode result = mbImpl->get_adjacencies( &ncell, 1, 4, false, adj_sets );
1385 if( MB_SUCCESS != result ) return 0;
1386
1387 EntityHandle dum_set;
1388 for( std::vector< EntityHandle >::iterator vit = adj_sets.begin(); vit != adj_sets.end(); ++vit )
1389 {
1390 if( mbImpl->tag_get_data( dualCurve_tag(), &( *vit ), 1, &dum_set ) != MB_TAG_NOT_FOUND ||
1391 mbImpl->tag_get_data( dualSurface_tag(), &( *vit ), 1, &dum_set ) != MB_TAG_NOT_FOUND )
1392 return *vit;
1393 }
1394
1395 return 0;
1396 }
1397
1398 //! set the dual surface or curve for an entity
1399 ErrorCode DualTool::set_dual_surface_or_curve( EntityHandle entity,
1400 const EntityHandle dual_hyperplane,
1401 const int dual_entity_dimension )
1402 {
1403 if( 1 == dual_entity_dimension )
1404 mbImpl->tag_set_data( dualCurve_tag(), &entity, 1, &dual_hyperplane );
1405 else if( 2 == dual_entity_dimension )
1406 mbImpl->tag_set_data( dualSurface_tag(), &entity, 1, &dual_hyperplane );
1407 else
1408 return MB_INDEX_OUT_OF_RANGE;
1409
1410 return MB_SUCCESS;
1411 }
1412
1413 //! return the corresponding dual entity
1414 EntityHandle DualTool::get_dual_entity( const EntityHandle this_ent ) const
1415 {
1416 EntityHandle dual_ent;
1417 ErrorCode result = mbImpl->tag_get_data( dualEntity_tag(), &this_ent, 1, &dual_ent );
1418 if( MB_SUCCESS != result || MB_TAG_NOT_FOUND == result )
1419 return 0;
1420 else
1421 return dual_ent;
1422 }
1423
1424 //! return the corresponding dual entity
1425 EntityHandle DualTool::get_extra_dual_entity( const EntityHandle this_ent )
1426 {
1427 EntityHandle dual_ent;
1428 ErrorCode result = mbImpl->tag_get_data( extraDualEntity_tag(), &this_ent, 1, &dual_ent );
1429 if( MB_SUCCESS != result || MB_TAG_NOT_FOUND == result )
1430 return 0;
1431 else
1432 return dual_ent;
1433 }
1434
1435 ErrorCode DualTool::atomic_pillow( EntityHandle odedge, EntityHandle& quad1, EntityHandle& quad2 )
1436 {
1437 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
1438
1439 if( debug_ap )
1440 {
1441 Range sets;
1442 Tag ms_tag;
1443
1444 ErrorCode result = mbImpl->tag_get_handle( "MATERIAL_SET", 1, MB_TYPE_INTEGER, ms_tag );
1445 if( MB_SUCCESS == result )
1446 {
1447 result = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &ms_tag, NULL, 1, sets );
1448 if( MB_SUCCESS == result ) result = mbImpl->delete_entities( sets );
1449 }
1450 }
1451
1452 std::cout << "-AP(";
1453 print_cell( odedge );
1454 std::cout << ")" << std::endl;
1455
1456 // perform an atomic pillow operation around dedge
1457
1458 // grab the quad before deleting the odedge
1459 quad1 = get_dual_entity( odedge );
1460 assert( 0 != quad1 );
1461
1462 // 0. get star 2cells around odedge (before odedge changes) and 3cells around
1463 // those 2cells (going to delete the 2cells, therefore need to delete the 3cells
1464 // that depend on those too)
1465 MeshTopoUtil mtu( mbImpl );
1466 Range star_cells, tmp_cells;
1467 ErrorCode result = mbImpl->get_adjacencies( &odedge, 1, 2, false, star_cells );RR;
1468 result = mbImpl->get_adjacencies( star_cells, 3, false, tmp_cells, Interface::UNION );RR;
1469 star_cells.merge( tmp_cells );
1470 star_cells.insert( odedge );
1471
1472 // tear down the dual entities which will be modified by the ap first
1473 result = delete_dual_entities( star_cells );RR;
1474
1475 // now change the quad to an ap
1476 std::vector< EntityHandle > verts;
1477 result = mbImpl->get_connectivity( &quad1, 1, verts );RR;
1478
1479 // get average position of vertices
1480 double coords[12], avg[3] = { 0.0, 0.0, 0.0 };
1481 result = mbImpl->get_coords( &verts[0], verts.size(), coords );RR;
1482 for( int i = 0; i < 4; i++ )
1483 {
1484 avg[0] += coords[3 * i];
1485 avg[1] += coords[3 * i + 1];
1486 avg[2] += coords[3 * i + 2];
1487 }
1488 for( int i = 0; i < 3; i++ )
1489 avg[i] *= 0.25;
1490
1491 // for each position, get a corresponding position 1/2 way to avg
1492 double new_coords[12];
1493 for( int i = 0; i < 4; i++ )
1494 {
1495 new_coords[3 * i] = avg[0] + .5 * ( coords[3 * i] - avg[0] );
1496 new_coords[3 * i + 1] = avg[1] + .5 * ( coords[3 * i + 1] - avg[1] );
1497 new_coords[3 * i + 2] = avg[2] + .5 * ( coords[3 * i + 2] - avg[2] );
1498 }
1499
1500 // make the 4 new vertices; store in vector long enough for hex connectivity
1501 for( int i = 0; i < 4; i++ )
1502 {
1503 verts.push_back( 0 );
1504 result = mbImpl->create_vertex( &new_coords[3 * i], verts[4 + i] );RR;
1505 }
1506
1507 // get the hexes connected to the quad
1508 Range hexes;
1509 result = mbImpl->get_adjacencies( &quad1, 1, 3, false, hexes );RR;
1510 assert( hexes.size() <= 2 );
1511
1512 // remove any explicit adjacency from the first hex, since that'll get connected
1513 // to the new quad; add adjacency between quad and second hex, if there is a 2nd
1514 result = mbImpl->remove_adjacencies( quad1, &( *hexes.begin() ), 1 );RR;
1515 if( hexes.size() == 2 )
1516 {
1517 result = mbImpl->add_adjacencies( quad1, &( *hexes.rbegin() ), 1, false );RR;
1518 }
1519
1520 // create the new, inner quad, and make it explicitly adjacent to 1st hex;
1521 // make the connectivity of this quad same as the original one
1522 std::vector< EntityHandle > tmp_verts;
1523 std::copy( verts.begin(), verts.end(), std::back_inserter( tmp_verts ) );
1524
1525 result = mbImpl->create_element( MBQUAD, &tmp_verts[0], 4, quad2 );RR;
1526 result = mbImpl->add_adjacencies( quad2, &( *hexes.begin() ), 1, false );RR;
1527
1528 // reverse the connectivity of the 1st hex
1529 std::reverse( verts.begin(), verts.begin() + 4 );
1530 std::reverse( verts.begin() + 4, verts.end() );
1531
1532 // now make two inner hexes; note connectivity array is flipped for the two hexes
1533 EntityHandle new_hexes[2];
1534 result = mbImpl->create_element( MBHEX, &verts[0], 8, new_hexes[0] );RR;
1535 result = mbImpl->create_element( MBHEX, &tmp_verts[0], 8, new_hexes[1] );RR;
1536
1537 // set the global id tag on the new hexes
1538 int new_hex_ids[2] = { maxHexId + 1, maxHexId + 2 };
1539 maxHexId += 2;
1540 result = mbImpl->tag_set_data( globalId_tag(), new_hexes, 2, new_hex_ids );
1541 if( MB_SUCCESS != result ) return result;
1542
1543 // by definition, quad1 is adj to new_hexes[0]
1544 result = mbImpl->add_adjacencies( quad1, &new_hexes[0], 1, false );RR;
1545 result = mbImpl->add_adjacencies( quad2, &new_hexes[1], 1, false );RR;
1546
1547 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
1548
1549 // now update the dual
1550 result = construct_hex_dual( &new_hexes[0], 2 );RR;
1551
1552 // get the new dual surface, by getting one of the edges between the center
1553 // and outer vertex rings
1554 Range new_edge;
1555 verts[1] = verts[4];
1556 result = mbImpl->get_adjacencies( &verts[0], 2, 1, false, new_edge );
1557 if( MB_SUCCESS != result || new_edge.size() != 1 ) return result;
1558
1559 return MB_SUCCESS;
1560 }
1561
1562 //! effect reverse atomic pillow operation
1563 ErrorCode DualTool::rev_atomic_pillow( EntityHandle pillow, Range& chords )
1564 {
1565 // get the dual entities associated with elements in the pillow; go through
1566 // the elements instead of the pillow sheet so you get all of them, not just
1567 // the ones on the sheet
1568 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
1569
1570 std::cout << "-AP(";
1571 print_cell( pillow );
1572 std::cout << ")" << std::endl;
1573
1574 Range dverts;
1575 ErrorCode result = get_dual_entities( pillow, NULL, NULL, &dverts, NULL, NULL );
1576 if( MB_SUCCESS != result ) return result;
1577 assert( 2 == dverts.size() );
1578
1579 EntityHandle hexes[2];
1580 result = mbImpl->tag_get_data( dualEntity_tag(), dverts, hexes );RR;
1581 assert( hexes[0] != 0 && hexes[1] != 0 );
1582
1583 std::vector< EntityHandle > dcells[4];
1584 Range pcells[4];
1585 std::copy( hexes, hexes + 2, range_inserter( pcells[3] ) );
1586 std::copy( dverts.begin(), dverts.end(), std::back_inserter( dcells[0] ) );
1587 for( int dim = 0; dim <= 2; dim++ )
1588 {
1589 result = mbImpl->get_adjacencies( hexes, 2, dim, false, pcells[dim], Interface::UNION );RR;
1590 dcells[3 - dim].resize( pcells[dim].size() );
1591 result = mbImpl->tag_get_data( dualEntity_tag(), pcells[dim], &dcells[3 - dim][0] );RR;
1592 }
1593
1594 // delete the dual entities which are part of the original pillow
1595 result = mbImpl->delete_entities( &pillow, 1 );
1596 if( MB_SUCCESS != result ) return result;
1597
1598 result = mbImpl->delete_entities( chords );
1599 if( MB_SUCCESS != result ) return result;
1600
1601 for( int i = 3; i >= 0; i-- )
1602 {
1603 result = delete_dual_entities( &dcells[i][0], dcells[i].size() );RR;
1604 }
1605
1606 // delete the primal entities inserted by the ap; be careful to get the right
1607 // faces, edges and vertices
1608 Range del_faces, del_edges, del_verts, tmp_faces, tmp_verts;
1609 // faces are the shared 5 and the 1 other one with greater handle (which
1610 // we know will be later in the range)
1611 result = mbImpl->get_adjacencies( hexes, 2, 2, false, del_faces );RR;
1612 assert( 5 == del_faces.size() );
1613 std::copy( pcells[2].begin(), pcells[2].end(), range_inserter( tmp_faces ) );
1614 tmp_faces = subtract( tmp_faces, del_faces );
1615 del_faces.insert( *tmp_faces.rbegin() );
1616 result = mbImpl->get_adjacencies( tmp_faces, 0, false, tmp_verts );RR;
1617 std::copy( pcells[0].begin(), pcells[0].end(), range_inserter( del_verts ) );
1618 del_verts = subtract( del_verts, tmp_verts );
1619 assert( 4 == del_verts.size() );
1620 result = mbImpl->get_adjacencies( del_verts, 1, false, del_edges, Interface::UNION );RR;
1621 assert( 8 == del_edges.size() );
1622
1623 result = mbImpl->delete_entities( hexes, 2 );RR;
1624 result = mbImpl->delete_entities( del_faces );RR;
1625 result = mbImpl->delete_entities( del_edges );RR;
1626 result = mbImpl->delete_entities( del_verts );RR;
1627
1628 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
1629
1630 // recompute the dual for the hexes on either side of the quad affected
1631 // by the ap removal
1632 Range tmp_hexes;
1633 result = mbImpl->get_adjacencies( tmp_verts, 3, false, tmp_hexes, Interface::UNION );RR;
1634 result = construct_hex_dual( tmp_hexes );RR;
1635
1636 return MB_SUCCESS;
1637 }
1638
1639 ErrorCode DualTool::delete_dual_entities( EntityHandle* entities, int num_entities )
1640 {
1641 Range tmp_ents;
1642 std::copy( entities, entities + num_entities, range_inserter( tmp_ents ) );
1643 return delete_dual_entities( tmp_ents );
1644 }
1645
1646 ErrorCode DualTool::delete_dual_entities( Range& entities )
1647 {
1648 if( entities.empty() ) return delete_whole_dual();
1649
1650 EntityHandle null_entity = 0;
1651 ErrorCode result;
1652 Range ents_to_delete;
1653
1654 while( !entities.empty() )
1655 {
1656 EntityHandle this_entity = entities.pop_back();
1657
1658 // reset the primal's dual entity
1659 EntityHandle primal = get_dual_entity( this_entity );
1660 if( get_dual_entity( primal ) == this_entity )
1661 {
1662 result = mbImpl->tag_set_data( dualEntity_tag(), &primal, 1, &null_entity );RR;
1663 }
1664 EntityHandle extra = get_extra_dual_entity( primal );
1665 if( 0 != extra )
1666 {
1667 result = mbImpl->tag_set_data( extraDualEntity_tag(), &primal, 1, &null_entity );RR;
1668 }
1669
1670 ents_to_delete.insert( this_entity );
1671
1672 // check for extra dual entities
1673 if( mbImpl->type_from_handle( this_entity ) == MBPOLYGON )
1674 {
1675 // for 2cell, might be a loop edge
1676 Range loop_edges;
1677 result = mbImpl->get_adjacencies( &this_entity, 1, 1, false, loop_edges );
1678 for( Range::iterator rit = loop_edges.begin(); rit != loop_edges.end(); ++rit )
1679 if( check_1d_loop_edge( *rit ) ) entities.insert( *rit );
1680 }
1681 else if( extra && extra != this_entity )
1682 // just put it on the list; primal for which we're extra has already been
1683 // reset to not point to extra entity
1684 ents_to_delete.insert( extra );
1685 }
1686
1687 // now delete the entities (sheets and chords will be updated automatically)
1688 return mbImpl->delete_entities( ents_to_delete );
1689 }
1690
1691 void DualTool::print_cell( EntityHandle cell )
1692 {
1693 const EntityHandle* connect;
1694 int num_connect;
1695 ErrorCode result = mbImpl->get_connectivity( cell, connect, num_connect );
1696 if( MB_SUCCESS != result ) return;
1697 bool first = true;
1698 EntityHandle primals[20];
1699 std::vector< int > ids;
1700
1701 assert( num_connect < 20 );
1702 result = mbImpl->tag_get_data( dualEntityTag, connect, num_connect, primals );
1703 if( MB_SUCCESS != result ) return;
1704 ids.resize( num_connect );
1705 result = mbImpl->tag_get_data( globalIdTag, primals, num_connect, &ids[0] );
1706 if( MB_SUCCESS != result ) return;
1707 for( int i = 0; i < num_connect; i++ )
1708 {
1709 if( !first ) std::cout << "-";
1710 EntityType this_type = mbImpl->type_from_handle( primals[i] );
1711 if( this_type == MBHEX )
1712 std::cout << "h";
1713 else if( this_type == MBQUAD )
1714 std::cout << "f";
1715 else
1716 std::cout << "u";
1717
1718 if( ids[i] != 0 )
1719 std::cout << ids[i];
1720 else
1721 std::cout << mbImpl->id_from_handle( primals[i] );
1722
1723 first = false;
1724 }
1725 }
1726
1727 ErrorCode DualTool::face_open_collapse( EntityHandle ocl, EntityHandle ocr )
1728 {
1729 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
1730
1731 MeshTopoUtil mtu( mbImpl );
1732
1733 std::cout << "OC(";
1734 print_cell( ocl );
1735 std::cout << ")-(";
1736 print_cell( ocr );
1737 std::cout << ")" << std::endl;
1738
1739 // get the primal entities we're dealing with
1740 EntityHandle split_quads[2] = { 0 }, split_edges[3] = { 0 }, split_nodes[2] = { 0 }, other_edges[6] = { 0 },
1741 other_nodes[6] = { 0 };
1742 Range hexes;
1743 ErrorCode result = foc_get_ents( ocl, ocr, split_quads, split_edges, split_nodes, hexes, other_edges, other_nodes );RR;
1744
1745 // get star entities around edges, separated into halves
1746 std::vector< EntityHandle > star_dp1[2], star_dp2[2];
1747 result = foc_get_stars( split_quads, split_edges, star_dp1, star_dp2 );RR;
1748
1749 if( MBQUAD != mbImpl->type_from_handle( split_quads[0] ) || MBQUAD != mbImpl->type_from_handle( split_quads[1] ) )
1750 return MB_TYPE_OUT_OF_RANGE;
1751
1752 result = foc_delete_dual( split_quads, split_edges, hexes );
1753 if( MB_SUCCESS != result ) return result;
1754
1755 EntityHandle new_quads[2], new_edges[3], new_nodes[2];
1756 result = split_pair_nonmanifold( split_quads, split_edges, split_nodes, star_dp1, star_dp2, other_edges,
1757 other_nodes, new_quads, new_edges, new_nodes );
1758 if( MB_SUCCESS != result ) return result;
1759
1760 // now merge entities, the C of foc
1761 EntityHandle keepit, deleteit;
1762 #define MIN( a, b ) ( ( a ) < ( b ) ? ( a ) : ( b ) )
1763 #define MAX( a, b ) ( ( a ) > ( b ) ? ( a ) : ( b ) )
1764 #define KEEP_DELETE( a, b, c, d ) \
1765 { \
1766 ( c ) = MIN( a, b ); \
1767 ( d ) = MAX( a, b ); \
1768 }
1769
1770 // find how many shared edges there were
1771 int num_shared_edges = ( split_edges[2] ? 3 : ( split_edges[1] ? 2 : 1 ) );
1772
1773 // first the node(s)
1774 for( int i = 0; i < 3 - num_shared_edges; i++ )
1775 {
1776 KEEP_DELETE( other_nodes[2 + 2 * i], other_nodes[3 + 2 * i], keepit, deleteit );
1777 result = mbImpl->merge_entities( keepit, deleteit, false, true );RR;
1778 }
1779
1780 // now the edges
1781 for( int i = 0; i < 4 - num_shared_edges; i++ )
1782 {
1783 KEEP_DELETE( other_edges[2 * i], other_edges[2 * i + 1], keepit, deleteit );
1784 result = mbImpl->merge_entities( keepit, deleteit, false, true );RR;
1785 }
1786
1787 // now the faces
1788 KEEP_DELETE( split_quads[0], split_quads[1], keepit, deleteit );
1789 result = mbImpl->merge_entities( keepit, deleteit, false, true );RR;
1790
1791 result = mbImpl->merge_entities( new_quads[0], new_quads[1], false, true );RR;
1792
1793 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
1794
1795 // reconstruct dual
1796 result = construct_hex_dual( hexes );
1797 if( MB_SUCCESS != result ) return result;
1798
1799 return check_dual_adjs();
1800 }
1801
1802 ErrorCode DualTool::foc_get_ents( EntityHandle ocl,
1803 EntityHandle ocr,
1804 EntityHandle* split_quads,
1805 EntityHandle* split_edges,
1806 EntityHandle* split_nodes,
1807 Range& hexes,
1808 EntityHandle* other_edges,
1809 EntityHandle* other_nodes )
1810 {
1811 // get the entities used for foc; ocl and ocr are dual 1-cells
1812 // representing quads to be split; returned from this function:
1813 // quads[2] - 2 quads to be split
1814 // split_edges[2] - edge(s) to be split (2nd is 0 if only one)
1815 // split_node - node to be split, if any (otherwise 0)
1816 // hexes - connected hexes to split_edges
1817 // other_edges[0], [1] - edges in quads[0] and [1] sharing node with
1818 // one end of split_edges[0]
1819 // other_edges[2], [3] - other end of split_edges[0] (or [1] if 2
1820 // split_edges)
1821 // other_edges[4], [5] - edges in quads[0], [1] opposite to split_edges[0]
1822 // other_nodes[0], [1] - nodes on other_edges[0], [1] not shared with
1823 // split_edges[0]
1824 // other_nodes[2], [3] - nodes on other_edges[2], [3] not shared with
1825 // split_edges[0] (if 2 split_edges, there's only 1 opposite node
1826 // in each split quad)
1827 // (for diagram, see Tim's notes from 11/12/07)
1828
1829 split_quads[0] = get_dual_entity( ocl );
1830 split_quads[1] = get_dual_entity( ocr );
1831 if( MBQUAD != mbImpl->type_from_handle( split_quads[0] ) || MBQUAD != mbImpl->type_from_handle( split_quads[1] ) )
1832 return MB_TYPE_OUT_OF_RANGE;
1833
1834 Range common_edges;
1835 ErrorCode result = mbImpl->get_adjacencies( split_quads, 2, 1, false, common_edges );
1836 if( MB_SUCCESS != result ) return result;
1837
1838 if( common_edges.empty() ) return MB_FAILURE;
1839 for( unsigned int i = 0; i < common_edges.size(); i++ )
1840 split_edges[i] = common_edges[i];
1841
1842 MeshTopoUtil mtu( mbImpl );
1843
1844 if( common_edges.size() == 3 )
1845 {
1846 // find other (non-shared) edges
1847 for( int i = 0; i < 2; i++ )
1848 {
1849 Range tmp_edges;
1850 result = mbImpl->get_adjacencies( &split_quads[i], 1, 1, false, tmp_edges );
1851 if( MB_SUCCESS != result ) return result;
1852 tmp_edges = subtract( tmp_edges, common_edges );
1853 assert( tmp_edges.size() == 1 );
1854 other_edges[i] = *tmp_edges.begin();
1855 }
1856 assert( other_edges[0] && other_edges[1] && other_edges[0] != other_edges[1] );
1857
1858 // arrange common edges so middle is in middle
1859 result = mtu.opposite_entity( split_quads[0], other_edges[0], split_edges[1] );RR;
1860 common_edges.erase( split_edges[1] );
1861 split_edges[0] = *common_edges.begin();
1862 split_edges[2] = *common_edges.rbegin();
1863 common_edges.insert( split_edges[1] );
1864
1865 // get split nodes and other nodes
1866 split_nodes[0] = mtu.common_entity( split_edges[0], split_edges[1], 0 );
1867 split_nodes[1] = mtu.common_entity( split_edges[2], split_edges[1], 0 );
1868 other_nodes[0] = mtu.common_entity( split_edges[0], other_edges[0], 0 );
1869 other_nodes[1] = mtu.common_entity( split_edges[2], other_edges[1], 0 );
1870
1871 assert( other_nodes[0] && other_nodes[1] && split_nodes[0] && split_nodes[1] );
1872 assert( split_edges[0] && split_edges[1] && split_edges[2] && split_edges[0] != split_edges[1] &&
1873 split_edges[1] != split_edges[2] && split_edges[0] != split_edges[2] );
1874 }
1875 else if( common_edges.size() == 2 )
1876 {
1877 // split node is shared by split edges
1878 split_nodes[0] = mtu.common_entity( split_edges[0], split_edges[1], 0 );
1879 if( 0 == split_nodes[0] ) return MB_FAILURE;
1880 // first two other nodes are on split edges opposite split node
1881 result = mtu.opposite_entity( split_edges[0], split_nodes[0], other_nodes[0] );RR;
1882 result = mtu.opposite_entity( split_edges[1], split_nodes[0], other_nodes[1] );RR;
1883 // over split quads:
1884 for( int i = 0; i < 2; i++ )
1885 {
1886 // 1st other edge is opposite second split edge
1887 result = mtu.opposite_entity( split_quads[i], split_edges[1], other_edges[i] );RR;
1888 // 2nd other edge is opposite first split edge
1889 result = mtu.opposite_entity( split_quads[i], split_edges[0], other_edges[2 + i] );RR;
1890 // last other node is opposite split node on split quad
1891 result = mtu.opposite_entity( split_quads[i], split_nodes[0], other_nodes[2 + i] );RR;
1892 }
1893 }
1894 else
1895 {
1896 const EntityHandle* connect;
1897 int num_connect;
1898 result = mbImpl->get_connectivity( split_edges[0], connect, num_connect );
1899 if( MB_SUCCESS != result ) return result;
1900 // other_nodes[0], [1] are on split edge
1901 other_nodes[0] = connect[0];
1902 other_nodes[1] = connect[1];
1903
1904 // for each of the split quads
1905 for( int i = 0; i < 2; i++ )
1906 {
1907 // get the other edge on the split quad adj to node 0 on the split edge, by getting
1908 // edges adj to split quad and node and removing split edge; that's other_edge[i]
1909 Range tmp_range1, tmp_range2;
1910 tmp_range1.insert( connect[0] );
1911 tmp_range1.insert( split_quads[i] );
1912 result = mbImpl->get_adjacencies( tmp_range1, 1, false, tmp_range2 );
1913 if( MB_SUCCESS != result ) return result;
1914 tmp_range2.erase( split_edges[0] );
1915 assert( tmp_range2.size() == 1 );
1916 other_edges[i] = *tmp_range2.begin();
1917 // get edge connected to other node on split edge & split quad; that's
1918 // opposite prev other_edges on the split quad; that's other_edges[4+i]
1919 result = mtu.opposite_entity( split_quads[i], other_edges[i], other_edges[4 + i] );RR;
1920 // get the edge on the split quad opposite the split edge; that's other_edges[2+i]
1921 result = mtu.opposite_entity( split_quads[i], split_edges[0], other_edges[2 + i] );RR;
1922 // get nodes on other side of split quad from split edge, by getting common
1923 // node between top/bottom edge and opposite edge
1924 other_nodes[2 + i] = mtu.common_entity( other_edges[i], other_edges[2 + i], 0 );
1925 other_nodes[4 + i] = mtu.common_entity( other_edges[4 + i], other_edges[2 + i], 0 );
1926 if( 0 == other_nodes[2 + i] || 0 == other_nodes[4 + i] ) return MB_FAILURE;
1927 }
1928 }
1929
1930 result = mbImpl->get_adjacencies( split_edges, common_edges.size(), 3, false, hexes, Interface::UNION );
1931 if( MB_SUCCESS != result ) return result;
1932
1933 assert( "split node not in other_nodes" && other_nodes[0] != split_nodes[0] && other_nodes[0] != split_nodes[1] &&
1934 other_nodes[1] != split_nodes[0] && other_nodes[1] != split_nodes[1] );
1935 assert( "each split node on an end of a split edge" && mtu.common_entity( other_nodes[0], split_edges[0], 0 ) &&
1936 ( ( ( split_edges[2] && mtu.common_entity( other_nodes[1], split_edges[2], 0 ) ) ||
1937 ( split_edges[1] && mtu.common_entity( other_nodes[1], split_edges[1], 0 ) ) ||
1938 mtu.common_entity( other_nodes[1], split_edges[0], 0 ) ) ) );
1939 assert( "opposite other edges meet at an other node" &&
1940 ( mtu.common_entity( other_edges[0], other_edges[1], 0 ) == other_nodes[0] ||
1941 ( split_edges[2] && mtu.common_entity( other_edges[0], other_edges[1], 0 ) == other_nodes[1] ) ) &&
1942 ( split_edges[2] ||
1943 ( split_edges[1] && mtu.common_entity( other_edges[2], other_edges[3], 0 ) == other_nodes[1] ) ||
1944 mtu.common_entity( other_edges[4], other_edges[5], 0 ) == other_nodes[1] ) );
1945
1946 return MB_SUCCESS;
1947 }
1948
1949 ErrorCode DualTool::split_pair_nonmanifold( EntityHandle* split_quads,
1950 EntityHandle* split_edges,
1951 EntityHandle* split_nodes,
1952 std::vector< EntityHandle >* star_dp1,
1953 std::vector< EntityHandle >* star_dp2,
1954 EntityHandle* /*other_edges*/,
1955 EntityHandle* /*other_nodes*/,
1956 EntityHandle* new_quads,
1957 EntityHandle* new_edges,
1958 EntityHandle* new_nodes )
1959 {
1960
1961 // if there's a bdy in the star around the shared edge(s), get the quads on that
1962 // bdy so we know which edges to merge after the split-nonmanifold
1963 MeshTopoUtil mtu( mbImpl );
1964 ErrorCode result;
1965
1966 // get which star the split faces are in, and choose the other one
1967 int new_side = -1;
1968 if( std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) != star_dp1[0].end() )
1969 new_side = 1;
1970 else if( std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) != star_dp1[1].end() )
1971 new_side = 0;
1972 assert( -1 != new_side );
1973 if( -1 == new_side ) return MB_FAILURE;
1974
1975 //=============== split faces
1976
1977 for( int i = 0; i < 2; i++ )
1978 {
1979 // get a hex in star_dp2[new_side] that's adj to this split quad, to tell
1980 // mtu which one the new quad should go with; there should be at least one,
1981 // if we have any hexes connected to the split quad
1982 EntityHandle gowith_hex = 0;
1983 for( std::vector< EntityHandle >::iterator vit = star_dp2[new_side].begin(); vit != star_dp2[new_side].end();
1984 ++vit )
1985 {
1986 if( mtu.common_entity( *vit, split_quads[i], 2 ) )
1987 {
1988 gowith_hex = *vit;
1989 break;
1990 }
1991 }
1992 assert( 0 != gowith_hex );
1993
1994 // split manifold each of the split_quads, and put the results on the merge list
1995 result =
1996 mtu.split_entities_manifold( split_quads + i, 1, new_quads + i, NULL, ( gowith_hex ? &gowith_hex : NULL ) );RR;
1997 }
1998
1999 // make ranges of faces which need to be explicitly adj to old, new
2000 // edge; faces come from stars and new_quads (which weren't in the stars);
2001 // new_quads go with side j, which does not have split quads
2002 Range tmp_addl_faces[2], addl_faces[2];
2003 for( int i = 0; i < 2; i++ )
2004 {
2005 std::copy( star_dp1[i].begin(), star_dp1[i].end(), range_inserter( tmp_addl_faces[i] ) );
2006 tmp_addl_faces[new_side].insert( new_quads[i] );
2007 }
2008 #ifndef NDEBUG
2009 bool cond1 = ( "split_quads on 1, new_quads on 0" &&
2010 tmp_addl_faces[0].find( split_quads[0] ) == tmp_addl_faces[0].end() &&
2011 tmp_addl_faces[0].find( split_quads[1] ) == tmp_addl_faces[0].end() &&
2012 tmp_addl_faces[1].find( split_quads[0] ) != tmp_addl_faces[1].end() &&
2013 tmp_addl_faces[1].find( split_quads[1] ) != tmp_addl_faces[1].end() &&
2014 tmp_addl_faces[0].find( new_quads[0] ) != tmp_addl_faces[0].end() &&
2015 tmp_addl_faces[0].find( new_quads[1] ) != tmp_addl_faces[0].end() &&
2016 tmp_addl_faces[1].find( new_quads[0] ) == tmp_addl_faces[1].end() &&
2017 tmp_addl_faces[1].find( new_quads[1] ) == tmp_addl_faces[1].end() ),
2018 cond2 = ( "split_quads on 0, new_quads on 1" &&
2019 tmp_addl_faces[0].find( split_quads[0] ) != tmp_addl_faces[0].end() &&
2020 tmp_addl_faces[0].find( split_quads[1] ) != tmp_addl_faces[0].end() &&
2021 tmp_addl_faces[1].find( split_quads[0] ) == tmp_addl_faces[1].end() &&
2022 tmp_addl_faces[1].find( split_quads[1] ) == tmp_addl_faces[1].end() &&
2023 tmp_addl_faces[0].find( new_quads[0] ) == tmp_addl_faces[0].end() &&
2024 tmp_addl_faces[0].find( new_quads[1] ) == tmp_addl_faces[0].end() &&
2025 tmp_addl_faces[1].find( new_quads[0] ) != tmp_addl_faces[1].end() &&
2026 tmp_addl_faces[1].find( new_quads[1] ) != tmp_addl_faces[1].end() );
2027
2028 assert( cond1 || cond2 );
2029 #endif
2030
2031 //=============== split edge(s)
2032 for( int j = 0; j < 3; j++ )
2033 {
2034 if( !split_edges[j] ) break;
2035
2036 // filter add'l faces to only those adj to split_edges[j]
2037 addl_faces[0] = tmp_addl_faces[0];
2038 addl_faces[1] = tmp_addl_faces[1];
2039 for( int i = 0; i < 2; i++ )
2040 {
2041 result = mbImpl->get_adjacencies( &split_edges[j], 1, 2, false, addl_faces[i] );RR;
2042 }
2043
2044 // split edge
2045 result = mtu.split_entity_nonmanifold( split_edges[j], addl_faces[1 - new_side], addl_faces[new_side],
2046 new_edges[j] );RR;
2047 }
2048
2049 //=============== split node(s)
2050
2051 for( int j = 0; j < 2; j++ )
2052 {
2053 if( !split_nodes[j] ) break;
2054
2055 // if we're splitting multiple edges, there might be other edges that have the split
2056 // node; also need to know which side they're on
2057 Range tmp_addl_edges[2];
2058 result = foc_get_addl_ents( star_dp1, star_dp2, split_edges, split_nodes[j], tmp_addl_edges );RR;
2059
2060 // also, we need to know which of the split/new edges go
2061 // with the split/new node; new edges go with side 0, split with 1
2062 for( int i = 0; i < 3; i++ )
2063 {
2064 if( !split_edges[i] ) break;
2065 tmp_addl_edges[new_side].insert( new_edges[i] );
2066 tmp_addl_edges[1 - new_side].insert( split_edges[i] );
2067 }
2068
2069 // same for star faces and hexes
2070 for( int i = 0; i < 2; i++ )
2071 {
2072 std::copy( star_dp1[i].begin(), star_dp1[i].end(), range_inserter( tmp_addl_edges[i] ) );
2073 std::copy( star_dp2[i].begin(), star_dp2[i].end(), range_inserter( tmp_addl_edges[i] ) );
2074 }
2075
2076 // finally, new quads
2077 for( int i = 0; i < 2; i++ )
2078 tmp_addl_edges[new_side].insert( new_quads[i] );
2079
2080 // filter the entities, keeping only the ones adjacent to this node
2081 Range addl_edges[2];
2082 for( int i = 0; i < 2; i++ )
2083 {
2084 for( Range::reverse_iterator rit = tmp_addl_edges[i].rbegin(); rit != tmp_addl_edges[i].rend(); ++rit )
2085 {
2086 if( mtu.common_entity( *rit, split_nodes[j], 0 ) ) addl_edges[i].insert( *rit );
2087 }
2088 }
2089
2090 // now split the node too
2091 result = mtu.split_entity_nonmanifold( split_nodes[j], addl_edges[1 - new_side], addl_edges[new_side],
2092 new_nodes[j] );RR;
2093 }
2094
2095 return MB_SUCCESS;
2096 }
2097
2098 ErrorCode DualTool::foc_get_addl_ents( std::vector< EntityHandle >* star_dp1,
2099 std::vector< EntityHandle >* /*star_dp2*/,
2100 EntityHandle* split_edges,
2101 EntityHandle split_node,
2102 Range* addl_ents )
2103 {
2104 // if we're splitting 2 edges, there might be other edges that have the split
2105 // node; also need to know which side they're on
2106
2107 // algorithm: for a given star_dp1 (faces) on a side:
2108 // - get all edges adj to all faces -> R1
2109 // - get all edges adj to split_node -> R2
2110 // - R3 = R1 & R2 (edges connected to split_node & adj to a star face)
2111 // - R3 -= split_edges (take split edges off addl_ents)
2112
2113 Range R2;
2114 MeshTopoUtil mtu( mbImpl );
2115 ErrorCode result = mbImpl->get_adjacencies( &split_node, 1, 1, false, R2 );RR;
2116 Range::iterator rit;
2117
2118 for( int i = 0; i < 2; i++ )
2119 {
2120 Range R1, R3;
2121 result = mbImpl->get_adjacencies( &star_dp1[i][0], star_dp1[i].size(), 1, false, R1, Interface::UNION );RR;
2122 R3 = intersect( R1, R2 );
2123 for( int j = 0; j < 3; j++ )
2124 if( split_edges[j] ) R3.erase( split_edges[j] );
2125 addl_ents[i].merge( R3 );
2126 }
2127
2128 return MB_SUCCESS;
2129 }
2130
2131 ErrorCode DualTool::foc_get_stars( EntityHandle* split_quads,
2132 EntityHandle* split_edges,
2133 std::vector< EntityHandle >* star_dp1,
2134 std::vector< EntityHandle >* star_dp2 )
2135 {
2136 bool on_bdy = false, on_bdy_tmp;
2137 ErrorCode result;
2138 MeshTopoUtil mtu( mbImpl );
2139
2140 // get the star around the split_edge
2141 std::vector< EntityHandle > qstar, hstar;
2142 unsigned int qpos = 0;
2143
2144 for( int i = 0; i < 3; i++ )
2145 {
2146 if( !split_edges[i] ) break;
2147
2148 // get the star around this split edge
2149 unsigned int qpos_tmp = 0;
2150 std::vector< EntityHandle > qstar_tmp, hstar_tmp;
2151 result = mtu.star_entities( split_edges[i], qstar_tmp, on_bdy_tmp, 0, &hstar_tmp );RR;
2152 // if we're on the bdy, add a null to the hex star too
2153 if( on_bdy_tmp )
2154 {
2155 assert( hstar_tmp.size() == qstar_tmp.size() - 1 );
2156 hstar_tmp.push_back( 0 );
2157 on_bdy = true;
2158 }
2159
2160 // get the position of first split quad in star
2161 while( qpos_tmp < qstar_tmp.size() && qstar_tmp[qpos_tmp] != split_quads[0] )
2162 qpos_tmp++;
2163 if( qpos_tmp == qstar_tmp.size() ) return MB_FAILURE;
2164
2165 bool forward;
2166 // 1st iteration is forward by definition
2167 if( 0 == i ) forward = true;
2168
2169 // need to be careful about direction on later iters
2170 else if( hstar[qpos] == hstar_tmp[qpos_tmp] )
2171 forward = true;
2172 else if( hstar[qpos] == hstar_tmp[( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size()] &&
2173 hstar_tmp[qpos_tmp] == hstar[( qpos + qstar.size() - 1 ) % qstar.size()] )
2174 forward = false;
2175 else
2176 return MB_FAILURE;
2177
2178 if( forward )
2179 {
2180 // 1st half of star
2181 // save hex right after split_quad[0] first
2182 star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
2183 qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
2184 while( qstar_tmp[qpos_tmp] != split_quads[1] )
2185 {
2186 star_dp1[0].push_back( qstar_tmp[qpos_tmp] );
2187 star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
2188 qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
2189 }
2190 // 2nd half of star
2191 // save hex right after split_quad[1] first
2192 star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
2193 qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
2194 while( qstar_tmp[qpos_tmp] != split_quads[0] )
2195 {
2196 star_dp1[1].push_back( qstar_tmp[qpos_tmp] );
2197 star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
2198 qpos_tmp = ( qpos_tmp + 1 ) % qstar_tmp.size();
2199 }
2200 }
2201 else
2202 {
2203 // go in reverse - take prev hex instead of current
2204 // one, and step in reverse
2205
2206 // save hex right after split_quad[0] first
2207 qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
2208 star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
2209 while( qstar_tmp[qpos_tmp] != split_quads[1] )
2210 {
2211 star_dp1[0].push_back( qstar_tmp[qpos_tmp] );
2212 qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
2213 star_dp2[0].push_back( hstar_tmp[qpos_tmp] );
2214 }
2215 // 2nd half of star
2216 // save hex right after split_quad[1] first
2217 qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
2218 star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
2219 while( qstar_tmp[qpos_tmp] != split_quads[0] )
2220 {
2221 star_dp1[1].push_back( qstar_tmp[qpos_tmp] );
2222 qpos_tmp = ( qpos_tmp + qstar_tmp.size() - 1 ) % qstar_tmp.size();
2223 star_dp2[1].push_back( hstar_tmp[qpos_tmp] );
2224 }
2225 }
2226
2227 if( 0 == i )
2228 {
2229 // if we're on the first iteration, save results and continue, other iters
2230 // get compared to this one
2231 qstar.swap( qstar_tmp );
2232 hstar.swap( hstar_tmp );
2233 on_bdy = on_bdy_tmp;
2234 qpos = qpos_tmp;
2235 }
2236 }
2237
2238 // split quads go on list with NULLs, if any, otherwise on 2nd
2239 if( on_bdy )
2240 {
2241 if( std::find( star_dp2[0].begin(), star_dp2[0].end(), 0 ) != star_dp2[0].end() )
2242 {
2243 // remove *all* the zeros
2244 star_dp2[0].erase( std::remove( star_dp2[0].begin(), star_dp2[0].end(), 0 ), star_dp2[0].end() );
2245 // put the split quads on this half
2246 star_dp1[0].push_back( split_quads[0] );
2247 star_dp1[0].push_back( split_quads[1] );
2248 }
2249 else
2250 {
2251 star_dp2[1].erase( std::remove( star_dp2[1].begin(), star_dp2[1].end(), 0 ), star_dp2[1].end() );
2252 // put the split quads on this half
2253 star_dp1[1].push_back( split_quads[0] );
2254 star_dp1[1].push_back( split_quads[1] );
2255 }
2256 }
2257 else
2258 {
2259 star_dp1[1].push_back( split_quads[0] );
2260 star_dp1[1].push_back( split_quads[1] );
2261 }
2262
2263 // some error checking
2264 if( !( ( ( std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) == star_dp1[0].end() &&
2265 std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[1] ) == star_dp1[0].end() &&
2266 std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) != star_dp1[1].end() &&
2267 std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[1] ) != star_dp1[1].end() ) ||
2268
2269 ( std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[0] ) == star_dp1[1].end() &&
2270 std::find( star_dp1[1].begin(), star_dp1[1].end(), split_quads[1] ) == star_dp1[1].end() &&
2271 std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[0] ) != star_dp1[0].end() &&
2272 std::find( star_dp1[0].begin(), star_dp1[0].end(), split_quads[1] ) != star_dp1[0].end() ) ) ) )
2273 {
2274 std::cerr << "foc_get_stars: both split quads should be on the same star list half and not "
2275 << "on the other, failed" << std::endl;
2276 return MB_FAILURE;
2277 }
2278
2279 if( !( std::find( star_dp2[0].begin(), star_dp2[0].end(), 0 ) == star_dp2[0].end() &&
2280 std::find( star_dp2[1].begin(), star_dp2[1].end(), 0 ) == star_dp2[1].end() ) )
2281 {
2282 std::cerr << "foc_get_stars: no NULLs on the hstar lists, failed";
2283 return MB_FAILURE;
2284 }
2285
2286 return MB_SUCCESS;
2287 }
2288
2289 ErrorCode DualTool::foc_delete_dual( EntityHandle* split_quads, EntityHandle* split_edges, Range& hexes )
2290 {
2291 // special delete dual procedure, because in some cases we need to delete
2292 // a sheet too since it'll get merged into another
2293
2294 // figure out whether we'll need to delete a sheet
2295 EntityHandle sheet1 = get_dual_hyperplane( get_dual_entity( split_edges[0] ) );
2296 if( split_edges[1] ) sheet1 = get_dual_hyperplane( get_dual_entity( split_edges[1] ) );
2297 EntityHandle chordl = get_dual_hyperplane( get_dual_entity( split_quads[0] ) );
2298 EntityHandle chordr = get_dual_hyperplane( get_dual_entity( split_quads[1] ) );
2299 assert( 0 != sheet1 && 0 != chordl && 0 != chordr );
2300 Range parentsl, parentsr;
2301 ErrorCode result = mbImpl->get_parent_meshsets( chordl, parentsl );
2302 if( MB_SUCCESS != result ) return result;
2303 result = mbImpl->get_parent_meshsets( chordr, parentsr );
2304 if( MB_SUCCESS != result ) return result;
2305 parentsl.erase( sheet1 );
2306 parentsr.erase( sheet1 );
2307
2308 // before deciding which one to delete, collect the other cells which must
2309 // be deleted, and all the chords/sheets they're on
2310 Range adj_ents, dual_ents, cells1or2;
2311 for( int i = 0; i < 3; i++ )
2312 {
2313 result = mbImpl->get_adjacencies( hexes, i, false, adj_ents, Interface::UNION );
2314 if( MB_SUCCESS != result ) return result;
2315 }
2316
2317 // cache any adjacent hexes, for rebuilding the dual later
2318 result = mbImpl->get_adjacencies( adj_ents, 3, false, hexes, Interface::UNION );
2319 if( MB_SUCCESS != result ) return result;
2320
2321 for( Range::iterator rit = adj_ents.begin(); rit != adj_ents.end(); ++rit )
2322 {
2323 EntityHandle this_ent = get_dual_entity( *rit );
2324 dual_ents.insert( this_ent );
2325 int dim = mbImpl->dimension_from_handle( this_ent );
2326 if( 1 == dim || 2 == dim ) cells1or2.insert( this_ent );
2327 }
2328
2329 Range dual_hps;
2330 for( Range::iterator rit = cells1or2.begin(); rit != cells1or2.end(); ++rit )
2331 dual_hps.insert( get_dual_hyperplane( *rit ) );
2332
2333 result = delete_dual_entities( dual_ents );
2334 if( MB_SUCCESS != result ) return result;
2335
2336 // now decide which sheet to delete (to be merged into the other)
2337 EntityHandle sheet_delete = 0;
2338 if( is_blind( *parentsl.begin() ) )
2339 sheet_delete = *parentsl.begin();
2340 else if( is_blind( *parentsr.begin() ) )
2341 sheet_delete = *parentsr.begin();
2342 else
2343 {
2344 // neither is blind, take the one with fewer cells
2345 Range tmp_ents;
2346 int numl, numr;
2347 result = mbImpl->get_number_entities_by_handle( *parentsl.begin(), numl );
2348 if( MB_SUCCESS != result ) return result;
2349 result = mbImpl->get_number_entities_by_handle( *parentsr.begin(), numr );
2350 if( MB_SUCCESS != result ) return result;
2351 sheet_delete = ( numl > numr ? *parentsr.begin() : *parentsl.begin() );
2352 }
2353 assert( 0 != sheet_delete );
2354
2355 // after deleting cells, check for empty chords & sheets, and delete those too
2356 for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit )
2357 {
2358 Range tmp_ents;
2359 result = mbImpl->get_entities_by_handle( *rit, tmp_ents );
2360 if( MB_SUCCESS != result ) return result;
2361 if( tmp_ents.empty() )
2362 {
2363 result = mbImpl->delete_entities( &( *rit ), 1 );
2364 if( MB_SUCCESS != result ) return result;
2365 }
2366 else if( *rit == sheet_delete )
2367 {
2368 // delete the sheet
2369 result = mbImpl->delete_entities( &( *rit ), 1 );
2370 if( MB_SUCCESS != result ) return result;
2371 }
2372 }
2373
2374 // now just to be safe, add the hexes bridge-adjacent across vertices
2375 // to the hexes we already have
2376 Range tmp_hexes;
2377 MeshTopoUtil mtu( mbImpl );
2378 for( Range::iterator rit = hexes.begin(); rit != hexes.end(); ++rit )
2379 {
2380 result = mtu.get_bridge_adjacencies( *rit, 0, 3, tmp_hexes );
2381 if( MB_SUCCESS != result ) return result;
2382 }
2383 hexes.merge( tmp_hexes );
2384
2385 return MB_SUCCESS;
2386 }
2387
2388 //! returns true if all vertices are dual to hexes (not faces)
2389 bool DualTool::is_blind( const EntityHandle chord_or_sheet )
2390 {
2391 // must be an entity set
2392 if( TYPE_FROM_HANDLE( chord_or_sheet ) != MBENTITYSET ) return false;
2393
2394 // get the vertices
2395 Range verts, ents;
2396 ErrorCode result = mbImpl->get_entities_by_handle( chord_or_sheet, ents );
2397 if( MB_SUCCESS != result || ents.empty() ) return false;
2398
2399 result = mbImpl->get_adjacencies( ents, 0, false, verts, Interface::UNION );
2400 if( MB_SUCCESS != result || verts.empty() ) return false;
2401
2402 for( Range::iterator rit = verts.begin(); rit != verts.end(); ++rit )
2403 {
2404 // get dual entity for this vertex
2405 EntityHandle dual_ent = get_dual_entity( *rit );
2406 if( 0 == dual_ent ) continue;
2407 if( TYPE_FROM_HANDLE( dual_ent ) == MBQUAD ) return false;
2408 }
2409
2410 // if none of the vertices' duals were quads, chord_or_sheet must be blind
2411 return true;
2412 }
2413
2414 //! given a 1-cell and a chord, return the neighboring vertices on the
2415 //! chord, in the same order as the 1-cell's vertices
2416 ErrorCode DualTool::get_opposite_verts( const EntityHandle middle_edge, const EntityHandle chord, EntityHandle* verts )
2417 {
2418 // get the edges on the chord, in order, and move to middle_edge
2419 std::vector< EntityHandle > chord_edges;
2420 const EntityHandle* connect;
2421 int num_connect;
2422
2423 ErrorCode result = mbImpl->get_entities_by_handle( chord, chord_edges );RR;
2424 std::vector< EntityHandle >::iterator vit = std::find( chord_edges.begin(), chord_edges.end(), middle_edge );
2425 result = mbImpl->get_connectivity( middle_edge, connect, num_connect );RR;
2426
2427 if(
2428 // middle_edge isn't on this chord
2429 vit == chord_edges.end() ||
2430 // chord only has 1 edge
2431 chord_edges.size() == 1 ||
2432 // middle_edge is at beginning or end and chord isn't blind
2433 ( ( vit == chord_edges.begin() || vit == chord_edges.end() - 1 ) && !is_blind( chord ) ) )
2434 return MB_FAILURE;
2435
2436 else if( chord_edges.size() == 2 )
2437 {
2438 // easier if it's a 2-edge blind chord, just get vertices in opposite order
2439 verts[0] = connect[1];
2440 verts[1] = connect[0];
2441 return MB_SUCCESS;
2442 }
2443
2444 // get vertices with the prev edge & subtract vertices of 1-cell
2445 if( vit == chord_edges.begin() )
2446 vit = chord_edges.end() - 1;
2447 else
2448 --vit;
2449 Range dum_connect, middle_connect;
2450 result = mbImpl->get_connectivity( &middle_edge, 1, middle_connect );RR;
2451 result = mbImpl->get_connectivity( &( *vit ), 1, dum_connect );RR;
2452 dum_connect = subtract( dum_connect, middle_connect );
2453 if( dum_connect.size() != 1 )
2454 {
2455 std::cerr << "Trouble traversing chord." << std::endl;
2456 return MB_FAILURE;
2457 }
2458
2459 // put in verts[0]
2460 verts[0] = *dum_connect.begin();
2461
2462 // same with prev edge
2463 ++vit;
2464 if( vit == chord_edges.end() ) vit = chord_edges.begin();
2465 ++vit;
2466 dum_connect.clear();
2467 result = mbImpl->get_connectivity( &( *vit ), 1, dum_connect );RR;
2468 dum_connect = subtract( dum_connect, middle_connect );
2469 if( dum_connect.size() != 1 )
2470 {
2471 std::cerr << "Trouble traversing chord." << std::endl;
2472 return MB_FAILURE;
2473 }
2474
2475 // put in verts[1]
2476 verts[1] = *dum_connect.begin();
2477
2478 // if verts[0] and 1st vertex of 1cell don't have common edge, switch verts
2479 MeshTopoUtil mtu( mbImpl );
2480 if( 0 == mtu.common_entity( verts[0], connect[0], 1 ) )
2481 {
2482 EntityHandle dum_h = verts[0];
2483 verts[0] = verts[1];
2484 verts[1] = dum_h;
2485 }
2486
2487 if( 0 == mtu.common_entity( verts[0], connect[0], 1 ) )
2488 {
2489 std::cerr << "Trouble traversing chord." << std::endl;
2490 return MB_FAILURE;
2491 }
2492
2493 return MB_SUCCESS;
2494 }
2495
2496 ErrorCode DualTool::get_dual_entities( const EntityHandle dual_ent,
2497 Range* dcells,
2498 Range* dedges,
2499 Range* dverts,
2500 Range* dverts_loop,
2501 Range* dedges_loop )
2502 {
2503 ErrorCode result = MB_SUCCESS;
2504
2505 if( NULL != dcells )
2506 {
2507 result = mbImpl->get_entities_by_type( dual_ent, MBPOLYGON, *dcells );
2508 if( MB_SUCCESS != result ) return result;
2509 }
2510
2511 if( NULL != dedges )
2512 {
2513 if( NULL != dcells )
2514 result = mbImpl->get_adjacencies( *dcells, 1, false, *dedges, Interface::UNION );
2515 else
2516 result = mbImpl->get_entities_by_type( dual_ent, MBEDGE, *dedges );
2517
2518 if( MB_SUCCESS != result ) return result;
2519 }
2520
2521 if( NULL != dverts )
2522 {
2523 if( NULL != dcells )
2524 result = mbImpl->get_adjacencies( *dcells, 0, false, *dverts, Interface::UNION );
2525 else if( NULL != dedges )
2526 result = mbImpl->get_adjacencies( *dedges, 0, false, *dverts, Interface::UNION );
2527 else
2528 {
2529 Range all_ents;
2530 result = mbImpl->get_entities_by_handle( dual_ent, all_ents );RR;
2531 result = mbImpl->get_adjacencies( all_ents, 0, false, *dverts, Interface::UNION );
2532 }
2533
2534 if( MB_SUCCESS != result ) return result;
2535 }
2536
2537 if( NULL != dverts_loop && NULL != dverts )
2538 {
2539 static std::vector< EntityHandle > dual_ents;
2540 dual_ents.resize( dverts->size() );
2541 result = mbImpl->tag_get_data( dualEntity_tag(), *dverts, &dual_ents[0] );
2542 if( MB_SUCCESS != result ) return result;
2543 Range::iterator rit;
2544 unsigned int i;
2545 for( rit = dverts->begin(), i = 0; rit != dverts->end(); ++rit, i++ )
2546 if( 0 != dual_ents[i] && mbImpl->type_from_handle( dual_ents[i] ) == MBQUAD ) dverts_loop->insert( *rit );
2547 }
2548
2549 if( NULL != dedges_loop && NULL != dedges )
2550 {
2551 static std::vector< EntityHandle > dual_ents;
2552 dual_ents.resize( dedges->size() );
2553 result = mbImpl->tag_get_data( dualEntity_tag(), *dedges, &dual_ents[0] );
2554 if( MB_SUCCESS != result ) return result;
2555 Range::iterator rit;
2556 unsigned int i;
2557 for( rit = dedges->begin(), i = 0; rit != dedges->end(); ++rit, i++ )
2558 if( 0 != dual_ents[i] && mbImpl->type_from_handle( dual_ents[i] ) == MBEDGE ) dedges_loop->insert( *rit );
2559 }
2560
2561 return result;
2562 }
2563
2564 ErrorCode DualTool::list_entities( const EntityHandle* entities, const int num_entities ) const
2565 {
2566 Range temp_range;
2567 ErrorCode result;
2568 if( NULL == entities && 0 == num_entities )
2569 return mbImpl->list_entities( entities, num_entities );
2570
2571 else if( NULL == entities && 0 < num_entities )
2572 {
2573
2574 // list all entities of all types
2575 std::cout << std::endl;
2576 for( EntityType this_type = MBVERTEX; this_type < MBMAXTYPE; this_type++ )
2577 {
2578 result = mbImpl->get_entities_by_type( 0, this_type, temp_range );
2579 if( MB_SUCCESS != result ) return result;
2580 }
2581 }
2582
2583 else
2584 {
2585 std::copy( entities, entities + num_entities, range_inserter( temp_range ) );
2586 }
2587
2588 return list_entities( temp_range );
2589 }
2590
2591 ErrorCode DualTool::list_entities( const Range& entities ) const
2592 {
2593 // now print each entity, listing the dual information first then calling Interface to do
2594 // the rest
2595 ErrorCode result = MB_SUCCESS, tmp_result;
2596 for( Range::const_iterator iter = entities.begin(); iter != entities.end(); ++iter )
2597 {
2598 EntityType this_type = TYPE_FROM_HANDLE( *iter );
2599 std::cout << CN::EntityTypeName( this_type ) << " " << ID_FROM_HANDLE( *iter ) << ":" << std::endl;
2600
2601 EntityHandle dual_ent = get_dual_entity( *iter );
2602 if( 0 != dual_ent )
2603 {
2604 std::cout << "Dual to " << CN::EntityTypeName( mbImpl->type_from_handle( dual_ent ) ) << " "
2605 << mbImpl->id_from_handle( dual_ent ) << std::endl;
2606 }
2607
2608 if( TYPE_FROM_HANDLE( *iter ) == MBENTITYSET )
2609 {
2610 EntityHandle chord = 0, sheet = 0;
2611 int id;
2612 result = mbImpl->tag_get_data( dualCurve_tag(), &( *iter ), 1, &chord );
2613 if( MB_SUCCESS != result ) return result;
2614 result = mbImpl->tag_get_data( dualSurface_tag(), &( *iter ), 1, &sheet );
2615 if( MB_SUCCESS != result ) return result;
2616 result = mbImpl->tag_get_data( globalId_tag(), &( *iter ), 1, &id );
2617 if( MB_SUCCESS != result ) return result;
2618
2619 if( 0 != chord ) std::cout << "(Dual chord " << id << ")" << std::endl;
2620 if( 0 != sheet ) std::cout << "(Dual sheet " << id << ")" << std::endl;
2621 }
2622
2623 tmp_result = mbImpl->list_entity( *iter );
2624 if( MB_SUCCESS != tmp_result ) result = tmp_result;
2625 }
2626
2627 return result;
2628 }
2629
2630 ErrorCode DualTool::face_shrink( EntityHandle odedge )
2631 {
2632 // some preliminary checking
2633 if( mbImpl->type_from_handle( odedge ) != MBEDGE ) return MB_TYPE_OUT_OF_RANGE;
2634
2635 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
2636
2637 std::cout << "FS(";
2638 print_cell( odedge );
2639 std::cout << ")" << std::endl;
2640
2641 EntityHandle quads[4], hexes[2];
2642 std::vector< EntityHandle > connects[4], side_quads[2];
2643
2644 // get the quads along the chord through the 2 hexes, and the vertices
2645 // for those quads
2646 ErrorCode result = fs_get_quads( odedge, quads, hexes, connects );
2647 if( MB_SUCCESS != result ) return result;
2648
2649 // flip/rotate connect arrays so they align & are same sense
2650 result = fs_check_quad_sense( hexes[0], quads[0], connects );
2651 if( MB_SUCCESS != result ) return result;
2652
2653 // get the quad loops along the "side" surfaces
2654 result = fs_get_quad_loops( hexes, connects, side_quads );
2655 if( MB_SUCCESS != result ) return result;
2656
2657 // ok, done with setup; now delete dual entities affected by this operation,
2658 // which is all the entities adjacent to vertices of dual edge
2659 Range adj_verts, adj_edges, dual_ents, cells1or2;
2660 MeshTopoUtil mtu( mbImpl );
2661 result = mtu.get_bridge_adjacencies( odedge, 0, 1, adj_edges );
2662 if( MB_SUCCESS != result ) return result;
2663 result = mbImpl->get_adjacencies( adj_edges, 0, false, adj_verts, Interface::UNION );
2664 if( MB_SUCCESS != result ) return result;
2665 for( int i = 1; i <= 3; i++ )
2666 {
2667 result = mbImpl->get_adjacencies( adj_verts, i, false, dual_ents, Interface::UNION );
2668 if( MB_SUCCESS != result ) return result;
2669 }
2670
2671 // before deleting dual, grab the 1- and 2-cells
2672 for( Range::iterator rit = dual_ents.begin(); rit != dual_ents.end(); ++rit )
2673 {
2674 int dim = mbImpl->dimension_from_handle( *rit );
2675 if( 1 == dim || 2 == dim ) cells1or2.insert( *rit );
2676 }
2677 Range dual_hps;
2678 for( Range::iterator rit = cells1or2.begin(); rit != cells1or2.end(); ++rit )
2679 dual_hps.insert( get_dual_hyperplane( *rit ) );
2680
2681 dual_ents.insert( odedge );
2682 result = delete_dual_entities( dual_ents );
2683 if( MB_SUCCESS != result ) return result;
2684
2685 // after deleting cells, check for empty chords & sheets, and delete those too
2686 for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit )
2687 {
2688 Range tmp_ents;
2689 result = mbImpl->get_entities_by_handle( *rit, tmp_ents );
2690 if( MB_SUCCESS != result ) return result;
2691 if( tmp_ents.empty() )
2692 {
2693 result = mbImpl->delete_entities( &( *rit ), 1 );
2694 if( MB_SUCCESS != result ) return result;
2695 }
2696 }
2697
2698 // remove any explicit adjacencies between side quads and hexes; don't check
2699 // for error, since there might not be adjacencies
2700 for( int i = 0; i < 4; i++ )
2701 {
2702 for( int j = 0; j < 2; j++ )
2703 {
2704 result = mbImpl->remove_adjacencies( side_quads[j][i], &hexes[j], 1 );
2705 }
2706 }
2707
2708 // make inner ring of vertices
2709 // get centroid of quad2
2710 double q2coords[12], avg[3] = { 0.0, 0.0, 0.0 };
2711 result = mbImpl->get_coords( &connects[1][0], 4, q2coords );
2712 if( MB_SUCCESS != result ) return result;
2713 for( int i = 0; i < 4; i++ )
2714 {
2715 avg[0] += q2coords[3 * i];
2716 avg[1] += q2coords[3 * i + 1];
2717 avg[2] += q2coords[3 * i + 2];
2718 }
2719 avg[0] *= .25;
2720 avg[1] *= .25;
2721 avg[2] *= .25;
2722 // position new vertices
2723 connects[3].resize( 4 );
2724 for( int i = 0; i < 4; i++ )
2725 {
2726 q2coords[3 * i] = avg[0] + .25 * ( q2coords[3 * i] - avg[0] );
2727 q2coords[3 * i + 1] = avg[1] + .25 * ( q2coords[3 * i + 1] - avg[1] );
2728 q2coords[3 * i + 2] = avg[2] + .25 * ( q2coords[3 * i + 2] - avg[2] );
2729 result = mbImpl->create_vertex( &q2coords[3 * i], connects[3][i] );
2730 if( MB_SUCCESS != result ) return result;
2731 }
2732
2733 // ok, now have the 4 connectivity arrays for 4 quads; construct hexes
2734 EntityHandle hconnect[8], new_hexes[4];
2735 int new_hex_ids[4];
2736
2737 for( int i = 0; i < 4; i++ )
2738 {
2739 int i1 = i, i2 = ( i + 1 ) % 4;
2740 hconnect[0] = connects[0][i1];
2741 hconnect[1] = connects[0][i2];
2742 hconnect[2] = connects[3][i2];
2743 hconnect[3] = connects[3][i1];
2744
2745 hconnect[4] = connects[1][i1];
2746 hconnect[5] = connects[1][i2];
2747 hconnect[6] = connects[2][i2];
2748 hconnect[7] = connects[2][i1];
2749
2750 result = mbImpl->create_element( MBHEX, hconnect, 8, new_hexes[i] );
2751 if( MB_SUCCESS != result ) return result;
2752
2753 // test for equiv entities from the side quads, and make explicit adjacencies
2754 // if there are any
2755 for( int j = 0; j < 2; j++ )
2756 {
2757 if( mtu.equivalent_entities( side_quads[j][i] ) )
2758 {
2759 result = mbImpl->add_adjacencies( side_quads[j][i], &new_hexes[i], 1, false );
2760 if( MB_SUCCESS != result ) return result;
2761 }
2762 }
2763
2764 new_hex_ids[i] = ++maxHexId;
2765 }
2766
2767 // set the global id tag on the new hexes
2768 result = mbImpl->tag_set_data( globalId_tag(), new_hexes, 4, new_hex_ids );
2769 if( MB_SUCCESS != result ) return result;
2770
2771 // now fixup other two hexes; start by getting hex through quads 0, 1
2772 // make this first hex switch to the other side, to make the dual look like
2773 // a hex push
2774 int tmp_ids[2];
2775 result = mbImpl->tag_get_data( globalId_tag(), hexes, 2, tmp_ids );
2776 if( MB_SUCCESS != result ) return result;
2777
2778 result = mbImpl->delete_entities( hexes, 2 );
2779 if( MB_SUCCESS != result ) return result;
2780 result = mbImpl->delete_entities( &quads[1], 1 );
2781 if( MB_SUCCESS != result ) return result;
2782 for( int i = 0; i < 4; i++ )
2783 {
2784 hconnect[i] = connects[3][i];
2785 hconnect[4 + i] = connects[2][i];
2786 }
2787 result = mbImpl->create_element( MBHEX, hconnect, 8, hexes[0] );
2788 if( MB_SUCCESS != result ) return result;
2789
2790 for( int i = 0; i < 4; i++ )
2791 {
2792 hconnect[i] = connects[0][i];
2793 hconnect[4 + i] = connects[3][i];
2794 }
2795 result = mbImpl->create_element( MBHEX, hconnect, 8, hexes[1] );
2796 if( MB_SUCCESS != result ) return result;
2797
2798 // check for and fix any explicit adjacencies on either end quad
2799 if( mtu.equivalent_entities( quads[0] ) ) mbImpl->add_adjacencies( quads[0], &hexes[1], 1, false );
2800 if( mtu.equivalent_entities( quads[2] ) ) mbImpl->add_adjacencies( quads[2], &hexes[0], 1, false );
2801
2802 // re-set the global ids for the hexes to what they were
2803 result = mbImpl->tag_set_data( globalId_tag(), hexes, 2, tmp_ids );
2804 if( MB_SUCCESS != result ) return result;
2805
2806 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
2807
2808 // now update the dual
2809 Range tmph;
2810 result = mtu.get_bridge_adjacencies( hexes[0], 0, 3, tmph );
2811 if( MB_SUCCESS != result ) return result;
2812 result = mtu.get_bridge_adjacencies( hexes[1], 0, 3, tmph );
2813 if( MB_SUCCESS != result ) return result;
2814 tmph.insert( hexes[1] );
2815 result = construct_hex_dual( tmph );
2816 if( MB_SUCCESS != result ) return result;
2817
2818 return result;
2819 }
2820
2821 ErrorCode DualTool::fs_get_quad_loops( EntityHandle* hexes,
2822 std::vector< EntityHandle >* connects,
2823 std::vector< EntityHandle >* side_quads )
2824 {
2825 for( int i = 0; i < 4; i++ )
2826 {
2827 for( int j = 0; j < 2; j++ )
2828 {
2829 Range adj_ents, dum_quads;
2830 adj_ents.insert( hexes[j] );
2831 adj_ents.insert( connects[j][i] );
2832 adj_ents.insert( connects[j][( i + 1 ) % 4] );
2833 adj_ents.insert( connects[j + 1][i] );
2834 adj_ents.insert( connects[j + 1][( i + 1 ) % 4] );
2835
2836 ErrorCode result = mbImpl->get_adjacencies( adj_ents, 2, false, dum_quads );
2837 if( MB_SUCCESS != result ) return result;
2838 assert( 1 == dum_quads.size() );
2839 side_quads[j].push_back( *dum_quads.begin() );
2840 }
2841 }
2842
2843 return MB_SUCCESS;
2844 }
2845
2846 ErrorCode DualTool::fs_check_quad_sense( EntityHandle hex0, EntityHandle quad0, std::vector< EntityHandle >* connects )
2847 {
2848 // check sense of 0th quad wrt hex; since sense is out of element,
2849 // switch if quad is NOT reversed wrt hex
2850 int dum1, dum2, sense = 0;
2851 ErrorCode result = mbImpl->side_number( hex0, quad0, dum1, sense, dum2 );
2852 if( MB_SUCCESS != result ) return result;
2853 assert( 0 != sense );
2854 if( 1 == sense )
2855 {
2856 // just switch sense of this one; others will get switched next
2857 EntityHandle dum = connects[0][0];
2858 connects[0][0] = connects[0][2];
2859 connects[0][2] = dum;
2860 }
2861
2862 // check sense of 1st, 2nd quads, rotate if necessary to align connect arrays
2863 int index0 = -1, index2 = -1, sense0 = 0, sense2 = 0;
2864 MeshTopoUtil mtu( mbImpl );
2865 for( int i = 0; i < 4; i++ )
2866 {
2867 if( 0 != mtu.common_entity( connects[0][0], connects[1][i], 1 ) )
2868 {
2869 index0 = i;
2870 if( 0 != mtu.common_entity( connects[0][1], connects[1][( i + 1 ) % 4], 1 ) )
2871 sense0 = 1;
2872 else if( 0 != mtu.common_entity( connects[0][1], connects[1][( i + 4 - 1 ) % 4], 1 ) )
2873 sense0 = -1;
2874 break;
2875 }
2876 }
2877
2878 assert( index0 != -1 && sense0 != 0 );
2879
2880 if( sense0 == -1 )
2881 {
2882 EntityHandle dumh = connects[1][0];
2883 connects[1][0] = connects[1][2];
2884 connects[1][2] = dumh;
2885 if( index0 % 2 == 0 ) index0 = ( index0 + 2 ) % 4;
2886 }
2887
2888 if( index0 != 0 )
2889 {
2890 std::vector< EntityHandle > tmpc;
2891 for( int i = 0; i < 4; i++ )
2892 tmpc.push_back( connects[1][( index0 + i ) % 4] );
2893 connects[1].swap( tmpc );
2894 }
2895
2896 for( int i = 0; i < 4; i++ )
2897 {
2898 if( 0 != mtu.common_entity( connects[1][0], connects[2][i], 1 ) )
2899 {
2900 index2 = i;
2901 if( 0 != mtu.common_entity( connects[1][1], connects[2][( i + 1 ) % 4], 1 ) )
2902 sense2 = 1;
2903 else if( 0 != mtu.common_entity( connects[1][1], connects[2][( i + 4 - 1 ) % 4], 1 ) )
2904 sense2 = -1;
2905 break;
2906 }
2907 }
2908
2909 assert( index2 != -1 && sense2 != 0 );
2910
2911 if( sense2 == -1 )
2912 {
2913 EntityHandle dumh = connects[2][0];
2914 connects[2][0] = connects[2][2];
2915 connects[2][2] = dumh;
2916 if( index2 % 2 == 0 ) index2 = ( index2 + 2 ) % 4;
2917 }
2918
2919 if( index2 != 0 )
2920 {
2921 std::vector< EntityHandle > tmpc;
2922 for( int i = 0; i < 4; i++ )
2923 tmpc.push_back( connects[2][( index2 + i ) % 4] );
2924 connects[2].swap( tmpc );
2925 }
2926
2927 return MB_SUCCESS;
2928 }
2929
2930 //! effect reverse face shrink operation
2931 ErrorCode DualTool::rev_face_shrink( EntityHandle odedge )
2932 {
2933 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
2934
2935 // some preliminary checking
2936 if( mbImpl->type_from_handle( odedge ) != MBEDGE ) return MB_TYPE_OUT_OF_RANGE;
2937
2938 std::cout << "-FS(";
2939 print_cell( odedge );
2940 std::cout << ")" << std::endl;
2941
2942 EntityHandle quads[4], hexes[2];
2943 std::vector< EntityHandle > connects[4], side_quads[2];
2944
2945 // get three quads (shared quad & 2 end quads), hexes, and quad
2946 // connects
2947 ErrorCode result = fs_get_quads( odedge, quads, hexes, connects );
2948 if( MB_SUCCESS != result ) return result;
2949
2950 // adjust sense & rotation so they're aligned, together & wrt first
2951 // hex
2952 result = fs_check_quad_sense( hexes[0], quads[0], connects );
2953 if( MB_SUCCESS != result ) return result;
2954
2955 result = fsr_get_fourth_quad( connects, side_quads );
2956 if( MB_SUCCESS != result )
2957 {
2958 std::cout << "Can't do -FS here, two hexes must be adjacent to ring of 4 hexes." << std::endl;
2959 return result;
2960 }
2961
2962 Range adj_ents, outer_hexes, all_adjs;
2963
2964 // first get the entities connected to interior 4 verts
2965 for( int i = 1; i <= 3; i++ )
2966 {
2967 result = mbImpl->get_adjacencies( &connects[1][0], 4, i, false, adj_ents, Interface::UNION );
2968 if( MB_SUCCESS != result ) return result;
2969 }
2970
2971 // next get all entities adjacent to those; these will have their dual
2972 // entities deleted
2973 for( int i = 0; i < 3; i++ )
2974 {
2975 result = mbImpl->get_adjacencies( adj_ents, i, false, all_adjs, Interface::UNION );
2976 if( MB_SUCCESS != result ) return result;
2977 }
2978
2979 // get the dual entities and delete them
2980 Range dual_ents, dual_hps;
2981 for( Range::iterator rit = all_adjs.begin(); rit != all_adjs.end(); ++rit )
2982 {
2983 EntityHandle this_ent = get_dual_entity( *rit );
2984 dual_ents.insert( this_ent );
2985 }
2986
2987 // before deleting dual, grab the 1- and 2-cells
2988 for( Range::iterator rit = dual_ents.begin(); rit != dual_ents.end(); ++rit )
2989 {
2990 int dim = mbImpl->dimension_from_handle( *rit );
2991 if( 1 == dim || 2 == dim ) dual_hps.insert( get_dual_hyperplane( *rit ) );
2992 }
2993
2994 result = delete_dual_entities( dual_ents );
2995 if( MB_SUCCESS != result ) return result;
2996
2997 // after deleting cells, check for empty chords & sheets, and delete those too
2998 for( Range::iterator rit = dual_hps.begin(); rit != dual_hps.end(); ++rit )
2999 {
3000 Range tmp_ents;
3001 result = mbImpl->get_entities_by_handle( *rit, tmp_ents );
3002 if( MB_SUCCESS != result ) return result;
3003 if( tmp_ents.empty() )
3004 {
3005 result = mbImpl->delete_entities( &( *rit ), 1 );
3006 if( MB_SUCCESS != result ) return result;
3007 }
3008 }
3009
3010 // before re-connecting two hexes, check for existing quad on 4th quad vertices;
3011 // if there is a quad there, need to add explicit adjs to any adj hexes, since
3012 // by definition there'll be another quad on those vertices
3013 bool need_explicit = false;
3014 Range adj_quads;
3015 result = mbImpl->get_adjacencies( &connects[3][0], 4, 2, false, adj_quads );
3016 if( MB_MULTIPLE_ENTITIES_FOUND == result || !adj_quads.empty() )
3017 {
3018 // there's already a quad for these 4 vertices; by definition,
3019 // we'll be creating equivalent entities, so that original quad
3020 // needs explicit adj's to its bounding elements
3021 need_explicit = true;
3022 for( Range::iterator rit = adj_quads.begin(); rit != adj_quads.end(); ++rit )
3023 {
3024 Range adj_hexes;
3025 result = mbImpl->get_adjacencies( &( *rit ), 1, 3, false, adj_hexes );RR;
3026 result = mbImpl->add_adjacencies( *rit, adj_hexes, false );RR;
3027 }
3028 }
3029
3030 // re-connect the two hexes
3031 std::vector< EntityHandle > new_connect;
3032 std::copy( connects[3].begin(), connects[3].end(), std::back_inserter( new_connect ) );
3033 std::copy( connects[2].begin(), connects[2].end(), std::back_inserter( new_connect ) );
3034 result = mbImpl->set_connectivity( hexes[0], &new_connect[0], 8 );
3035 if( MB_SUCCESS != result ) return result;
3036
3037 new_connect.clear();
3038 std::copy( connects[0].begin(), connects[0].end(), std::back_inserter( new_connect ) );
3039 std::copy( connects[3].begin(), connects[3].end(), std::back_inserter( new_connect ) );
3040 result = mbImpl->set_connectivity( hexes[1], &new_connect[0], 8 );
3041 if( MB_SUCCESS != result ) return result;
3042
3043 // test for equiv entities from the side quads, and make explicit
3044 // adjacencies if there are any
3045 MeshTopoUtil mtu( mbImpl );
3046 for( int j = 0; j < 2; j++ )
3047 {
3048 for( int i = 0; i < 4; i++ )
3049 {
3050 if( mtu.equivalent_entities( side_quads[j][i] ) )
3051 {
3052 result = mbImpl->add_adjacencies( side_quads[j][i], &hexes[j], 1, false );
3053 if( MB_SUCCESS != result ) return result;
3054 }
3055 }
3056 }
3057
3058 // remove hexes we want to keep
3059 adj_ents.erase( hexes[0] );
3060 adj_ents.erase( hexes[1] );
3061
3062 // delete the other interior entities
3063 result = mbImpl->delete_entities( adj_ents );
3064 if( MB_SUCCESS != result ) return result;
3065
3066 EntityHandle new_quad;
3067 result = mbImpl->create_element( MBQUAD, &connects[3][0], 4, new_quad );RR;
3068 if( need_explicit )
3069 {
3070 result = mbImpl->add_adjacencies( new_quad, hexes, 2, false );RR;
3071 }
3072
3073 if( debug_ap ) ( (Core*)mbImpl )->check_adjacencies();
3074
3075 // now update the dual
3076 result = construct_hex_dual( hexes, 2 );
3077 if( MB_SUCCESS != result ) return result;
3078
3079 return MB_SUCCESS;
3080 }
3081
3082 ErrorCode DualTool::fsr_get_fourth_quad( std::vector< EntityHandle >* connects,
3083 std::vector< EntityHandle >* side_quads )
3084 {
3085 // given the first three quad connectivities in ordered vectors, get the fourth,
3086 // where the fourth is really the 4 vertices originally shared by the 2 hexes
3087 // before the face shrink on them
3088
3089 // vertex on 4th quad is in quad adj to other 3 verts
3090 for( int i = 0; i < 4; i++ )
3091 {
3092 Range start_verts, tmp_verts, quads;
3093 for( int j = 0; j < 3; j++ )
3094 start_verts.insert( connects[j][i] );
3095 ErrorCode result = mbImpl->get_adjacencies( start_verts, 2, false, quads );
3096 if( MB_SUCCESS != result ) return result;
3097 assert( quads.size() == 1 );
3098 result = mbImpl->get_adjacencies( &( *quads.begin() ), 1, 0, false, tmp_verts );RR;
3099 tmp_verts = subtract( tmp_verts, start_verts );
3100 assert( 1 == tmp_verts.size() );
3101 connects[3].push_back( *tmp_verts.begin() );
3102 }
3103
3104 // now get the side quads
3105 for( int i = 0; i < 4; i++ )
3106 {
3107 Range dum_ents, hexes;
3108 dum_ents.insert( connects[1][i] );
3109 dum_ents.insert( connects[1][( i + 1 ) % 4] );
3110 dum_ents.insert( connects[3][i] );
3111 ErrorCode result = mbImpl->get_adjacencies( dum_ents, 3, false, hexes );
3112 if( MB_SUCCESS != result ) return result;
3113 assert( 1 == hexes.size() );
3114
3115 hexes.insert( connects[0][i] );
3116 hexes.insert( connects[0][( i + 1 ) % 4] );
3117 hexes.insert( connects[3][i] );
3118 hexes.insert( connects[3][( i + 1 ) % 4] );
3119 dum_ents.clear();
3120 result = mbImpl->get_adjacencies( hexes, 2, false, dum_ents );
3121 if( MB_SUCCESS != result ) return result;
3122 assert( dum_ents.size() == 1 );
3123 side_quads[0].push_back( *dum_ents.begin() );
3124
3125 hexes.erase( connects[0][i] );
3126 hexes.erase( connects[0][( i + 1 ) % 4] );
3127 hexes.insert( connects[2][i] );
3128 hexes.insert( connects[2][( i + 1 ) % 4] );
3129 dum_ents.clear();
3130 result = mbImpl->get_adjacencies( hexes, 2, false, dum_ents );
3131 if( MB_SUCCESS != result ) return result;
3132 side_quads[1].push_back( *dum_ents.begin() );
3133 }
3134
3135 return MB_SUCCESS;
3136 }
3137
3138 ErrorCode DualTool::fs_get_quads( EntityHandle odedge,
3139 EntityHandle* quads,
3140 EntityHandle* hexes,
3141 std::vector< EntityHandle >* connects )
3142 {
3143 // need to get the three quads along the chord
3144 EntityHandle chord = get_dual_hyperplane( odedge );
3145 if( 0 == chord ) return MB_FAILURE;
3146
3147 std::vector< EntityHandle > edges;
3148 ErrorCode result = mbImpl->get_entities_by_handle( chord, edges );
3149 if( MB_FAILURE == result ) return result;
3150
3151 std::vector< EntityHandle >::iterator vit = std::find( edges.begin(), edges.end(), odedge );
3152 // shouldn't be first or last edge on chord
3153 if( vit == edges.end() || *edges.begin() == *vit || *edges.rbegin() == *vit ) return MB_FAILURE;
3154
3155 // get quads/connectivity for first 3 quads
3156 quads[0] = get_dual_entity( *( vit - 1 ) );
3157 quads[1] = get_dual_entity( *vit );
3158 quads[2] = get_dual_entity( *( vit + 1 ) );
3159 for( int i = 0; i < 3; i++ )
3160 {
3161 result = mbImpl->get_connectivity( &quads[i], 1, connects[i], true );
3162 if( MB_SUCCESS != result ) return result;
3163 }
3164
3165 Range tmph;
3166 result = mbImpl->get_adjacencies( quads, 2, 3, false, tmph );
3167 if( MB_SUCCESS != result ) return result;
3168 assert( tmph.size() == 1 );
3169 hexes[0] = *tmph.begin();
3170
3171 tmph.clear();
3172 result = mbImpl->get_adjacencies( &quads[1], 2, 3, false, tmph );
3173 if( MB_SUCCESS != result ) return result;
3174 assert( tmph.size() == 1 );
3175 hexes[1] = *tmph.begin();
3176
3177 return MB_SUCCESS;
3178 }
3179
3180 ErrorCode DualTool::delete_whole_dual()
3181 {
3182 // delete dual hyperplanes
3183 Range dual_surfs, dual_curves;
3184 ErrorCode result = this->get_dual_hyperplanes( mbImpl, 2, dual_surfs );RR;
3185 result = mbImpl->delete_entities( dual_surfs );RR;
3186 result = this->get_dual_hyperplanes( mbImpl, 1, dual_curves );RR;
3187 result = mbImpl->delete_entities( dual_curves );RR;
3188
3189 // gather up all dual entities
3190 Range dual_ents;
3191 result = mbImpl->get_entities_by_type_and_tag( 0, MBVERTEX, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
3192 result = mbImpl->get_entities_by_type_and_tag( 0, MBEDGE, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
3193 result = mbImpl->get_entities_by_type_and_tag( 0, MBPOLYGON, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
3194 result =
3195 mbImpl->get_entities_by_type_and_tag( 0, MBPOLYHEDRON, &isDualCellTag, NULL, 1, dual_ents, Interface::UNION );RR;
3196
3197 // delete them, in reverse order of dimension
3198 ErrorCode tmp_result;
3199 for( Range::reverse_iterator rit = dual_ents.rbegin(); rit != dual_ents.rend(); ++rit )
3200 {
3201 tmp_result = mbImpl->delete_entities( &( *rit ), 1 );
3202 if( MB_SUCCESS != tmp_result ) result = tmp_result;
3203 }
3204 RR;
3205
3206 // delete dual-related tags
3207 if( 0 != dualSurfaceTag )
3208 {
3209 tmp_result = mbImpl->tag_delete( dualSurfaceTag );
3210 if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
3211 }
3212 if( 0 != dualCurveTag )
3213 {
3214 tmp_result = mbImpl->tag_delete( dualCurveTag );
3215 if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
3216 }
3217 if( 0 != dualEntityTag )
3218 {
3219 tmp_result = mbImpl->tag_delete( dualEntityTag );
3220 if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
3221 }
3222 if( 0 != extraDualEntityTag )
3223 {
3224 tmp_result = mbImpl->tag_delete( extraDualEntityTag );
3225 if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
3226 }
3227 if( 0 != dualGraphicsPointTag )
3228 {
3229 tmp_result = mbImpl->tag_delete( dualGraphicsPointTag );
3230 if( MB_SUCCESS != tmp_result && MB_TAG_NOT_FOUND != tmp_result ) result = tmp_result;
3231 }
3232
3233 return MB_SUCCESS;
3234 }
3235
3236 ErrorCode DualTool::check_dual_adjs()
3237 {
3238 // check primal-dual correspondence
3239
3240 // get the primal entities
3241 Range pents[4];
3242 ErrorCode result = mbImpl->get_entities_by_type( 0, MBHEX, pents[3] );RR;
3243 for( int i = 2; i >= 0; i-- )
3244 {
3245 result = mbImpl->get_adjacencies( pents[3], 2, false, pents[2], Interface::UNION );RR;
3246 }
3247
3248 // for each primal entity of dimension pd
3249 #define PRENT( ent ) CN::EntityTypeName( TYPE_FROM_HANDLE( ent ) ) << " " << ID_FROM_HANDLE( ent )
3250 ErrorCode overall_result = MB_SUCCESS;
3251 for( int pd = 1; pd <= 3; pd++ )
3252 {
3253 for( Range::iterator prit = pents[pd].begin(); prit != pents[pd].end(); ++prit )
3254 {
3255 // get corresponding dual entity of dimension dd = 3-pd
3256 EntityHandle dual_ent = get_dual_entity( *prit );
3257 if( 0 == dual_ent ) std::cerr << "Problem getting dual entity for " << PRENT( *prit ) << std::endl;
3258
3259 // for each sub dimension sd = 0..pd-1
3260 for( int sd = 0; sd < pd; sd++ )
3261 {
3262 Range R1, R2, R3;
3263 // R1 = entities bounding primal entity of dim sd
3264 result = mbImpl->get_adjacencies( &( *prit ), 1, sd, false, R1 );RR;
3265
3266 // R2 = entities bounded by dual entity, of dim 3-sd
3267 result = mbImpl->get_adjacencies( &dual_ent, 1, 3 - sd, false, R2 );RR;
3268
3269 if( R1.size() != R2.size() )
3270 {
3271 std::cerr << PRENT( *prit ) << ": number of adj ents in "
3272 << "primal/dual don't agree for dimension " << sd << "." << std::endl;
3273 overall_result = MB_FAILURE;
3274 }
3275
3276 // for each entity in R1, get its dual and look for it in R2
3277 for( Range::iterator r1it = R1.begin(); r1it != R1.end(); ++r1it )
3278 {
3279 EntityHandle tmp_dual = get_dual_entity( *r1it );
3280 if( R2.find( tmp_dual ) == R2.end() )
3281 {
3282 std::cerr << PRENT( *prit ) << ": adj entity " << PRENT( *r1it ) << " isn't adjacent in dual."
3283 << std::endl;
3284 overall_result = MB_FAILURE;
3285 }
3286 }
3287 // ditto for R2
3288 for( Range::iterator r2it = R2.begin(); r2it != R2.end(); ++r2it )
3289 {
3290 EntityHandle tmp_prim = get_dual_entity( *r2it );
3291 if( R1.find( tmp_prim ) == R1.end() )
3292 {
3293 std::cerr << PRENT( *prit ) << ": adj entity " << PRENT( *r2it ) << " isn't adjacent in primal."
3294 << std::endl;
3295 overall_result = MB_FAILURE;
3296 }
3297 }
3298 }
3299 }
3300 }
3301
3302 return overall_result;
3303 }
3304
3305 } // namespace moab
1.9.1.