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