 |
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 #include "moab/MergeMesh.hpp"
2
3 #include "moab/Skinner.hpp"
4 #include "moab/AdaptiveKDTree.hpp"
5 #include "moab/Range.hpp"
6 #include "moab/CartVect.hpp"
7
8 #include "Internals.hpp"
9 #include <vector>
10 #include <algorithm>
11 #include <string>
12 #include <vector>
13 #include <cassert>
14 #include <iostream>
15 #include <iomanip>
16
17 #include <cstdlib>
18
19 namespace moab
20 {
21
22 MergeMesh::MergeMesh( Interface* impl, bool printErrorIn )
23 : mbImpl( impl ), mbMergeTag( 0 ), mergeTol( 0.001 ), mergeTolSq( 0.000001 ), printError( printErrorIn )
24 {
25 }
26
27 MergeMesh::~MergeMesh()
28 {
29 if( mbMergeTag ) mbImpl->tag_delete( mbMergeTag );
30 mbMergeTag = NULL;
31 }
32
33 ErrorCode MergeMesh::merge_entities( EntityHandle* elems,
34 int elems_size,
35 const double merge_tol,
36 const int do_merge,
37 const int update_sets,
38 Tag merge_tag,
39 bool do_higher_dim )
40 {
41 mergeTol = merge_tol;
42 mergeTolSq = merge_tol * merge_tol;
43 Range tmp_elems;
44 tmp_elems.insert( elems, elems + elems_size );
45 ErrorCode result = merge_entities( tmp_elems, merge_tol, do_merge, update_sets, (Tag)merge_tag, do_higher_dim );
46
47 return result;
48 }
49
50 /* This function appears to be not necessary after MOAB conversion
51
52 void MergeMesh::perform_merge(iBase_TagHandle merge_tag)
53 {
54 // put into a range
55 ErrorCode result = perform_merge((Tag) merge_tag);
56 if (result != MB_SUCCESS)
57 throw MKException(iBase_FAILURE, "");
58 }*/
59
60 ErrorCode MergeMesh::merge_entities( Range& elems,
61 const double merge_tol,
62 const int do_merge,
63 const int,
64 Tag merge_tag,
65 bool merge_higher_dim )
66 {
67 // If merge_higher_dim is true, do_merge must also be true
68 if( merge_higher_dim && !do_merge )
69 {
70 return MB_FAILURE;
71 }
72
73 mergeTol = merge_tol;
74 mergeTolSq = merge_tol * merge_tol;
75
76 // get the skin of the entities
77 Skinner skinner( mbImpl );
78 Range skin_range;
79 ErrorCode result = skinner.find_skin( 0, elems, 0, skin_range, false, false );
80 if( MB_SUCCESS != result ) return result;
81
82 // create a tag to mark merged-to entity; reuse tree_root
83 EntityHandle tree_root = 0;
84 if( 0 == merge_tag )
85 {
86 result = mbImpl->tag_get_handle( "__merge_tag", 1, MB_TYPE_HANDLE, mbMergeTag, MB_TAG_DENSE | MB_TAG_EXCL,
87 &tree_root );
88 if( MB_SUCCESS != result ) return result;
89 }
90 else
91 mbMergeTag = merge_tag;
92
93 // build a kd tree with the vertices
94 AdaptiveKDTree kd( mbImpl );
95 result = kd.build_tree( skin_range, &tree_root );
96 if( MB_SUCCESS != result ) return result;
97
98 // find matching vertices, mark them
99 result = find_merged_to( tree_root, kd, mbMergeTag );
100 if( MB_SUCCESS != result ) return result;
101
102 // merge them if requested
103 if( do_merge )
104 {
105 result = perform_merge( mbMergeTag );
106 if( MB_SUCCESS != result ) return result;
107 }
108
109 if( merge_higher_dim && deadEnts.size() != 0 )
110 {
111 result = merge_higher_dimensions( elems );
112 if( MB_SUCCESS != result ) return result;
113 }
114
115 return MB_SUCCESS;
116 }
117
118 ErrorCode MergeMesh::merge_all( EntityHandle meshset, const double merge_tol )
119 {
120 ErrorCode rval;
121 if( 0 == mbMergeTag )
122 {
123 EntityHandle def_val = 0;
124 rval = mbImpl->tag_get_handle( "__merge_tag", 1, MB_TYPE_HANDLE, mbMergeTag, MB_TAG_DENSE | MB_TAG_EXCL,
125 &def_val );MB_CHK_ERR( rval );
126 }
127 // get all entities;
128 // get all vertices connected
129 // build a kdtree
130 // find merged to
131 mergeTol = merge_tol;
132 mergeTolSq = merge_tol * merge_tol;
133
134 // get all vertices
135 Range entities;
136 rval = mbImpl->get_entities_by_handle( meshset, entities, /*recursive*/ true );MB_CHK_ERR( rval );
137 Range sets = entities.subset_by_type( MBENTITYSET );
138 entities = subtract( entities, sets );
139 Range verts;
140 rval = mbImpl->get_connectivity( entities, verts );MB_CHK_ERR( rval );
141
142 // build a kd tree with the vertices
143 AdaptiveKDTree kd( mbImpl );
144 EntityHandle tree_root;
145 rval = kd.build_tree( verts, &tree_root );MB_CHK_ERR( rval );
146 // find matching vertices, mark them
147 rval = find_merged_to( tree_root, kd, mbMergeTag );MB_CHK_ERR( rval );
148
149 rval = perform_merge( mbMergeTag );MB_CHK_ERR( rval );
150
151 if( deadEnts.size() != 0 )
152 {
153 rval = merge_higher_dimensions( entities );MB_CHK_ERR( rval );
154 }
155 return MB_SUCCESS;
156 }
157 ErrorCode MergeMesh::perform_merge( Tag merge_tag )
158 {
159 // we start with an empty range of vertices that are "merged to"
160 // they are used (eventually) for higher dim entities
161 mergedToVertices.clear();
162 ErrorCode result;
163 if( deadEnts.size() == 0 )
164 {
165 if( printError ) std::cout << "\nWarning: Geometries don't have a common face; Nothing to merge" << std::endl;
166 return MB_SUCCESS; // nothing to merge carry on with the program
167 }
168 if( mbImpl->type_from_handle( *deadEnts.begin() ) != MBVERTEX ) return MB_FAILURE;
169 std::vector< EntityHandle > merge_tag_val( deadEnts.size() );
170 Range deadEntsRange;
171 std::copy( deadEnts.rbegin(), deadEnts.rend(), range_inserter( deadEntsRange ) );
172 result = mbImpl->tag_get_data( merge_tag, deadEntsRange, &merge_tag_val[0] );
173 if( MB_SUCCESS != result ) return result;
174
175 std::set< EntityHandle >::iterator rit;
176 unsigned int i;
177 for( rit = deadEnts.begin(), i = 0; rit != deadEnts.end(); ++rit, i++ )
178 {
179 assert( merge_tag_val[i] );
180 if( MBVERTEX == TYPE_FROM_HANDLE( merge_tag_val[i] ) ) mergedToVertices.insert( merge_tag_val[i] );
181 result = mbImpl->merge_entities( merge_tag_val[i], *rit, false, false );
182 if( MB_SUCCESS != result )
183 {
184 return result;
185 }
186 }
187 result = mbImpl->delete_entities( deadEntsRange );
188 return result;
189 }
190 // merge vertices according to an input tag
191 // merge them if the tags are equal
192 struct handle_id
193 {
194 EntityHandle eh;
195 int val;
196 };
197
198 // handle structure comparison function for qsort
199 // if the id is the same , compare the handle.
200 int compare_handle_id( const void* a, const void* b )
201 {
202
203 handle_id* ia = (handle_id*)a;
204 handle_id* ib = (handle_id*)b;
205 if( ia->val == ib->val )
206 {
207 return ( ia->eh < ib->eh ) ? -1 : 1;
208 }
209 else
210 {
211 return ( ia->val - ib->val );
212 }
213 }
214
215 ErrorCode MergeMesh::merge_using_integer_tag( Range& verts, Tag user_tag, Tag merge_tag )
216 {
217 ErrorCode rval;
218 DataType tag_type;
219 rval = mbImpl->tag_get_data_type( user_tag, tag_type );
220 if( rval != MB_SUCCESS || tag_type != MB_TYPE_INTEGER ) return MB_FAILURE;
221
222 std::vector< int > vals( verts.size() );
223 rval = mbImpl->tag_get_data( user_tag, verts, &vals[0] );
224 if( rval != MB_SUCCESS ) return rval;
225
226 if( 0 == merge_tag )
227 {
228 EntityHandle def_val = 0;
229 rval = mbImpl->tag_get_handle( "__merge_tag", 1, MB_TYPE_HANDLE, mbMergeTag, MB_TAG_DENSE | MB_TAG_EXCL,
230 &def_val );
231 if( MB_SUCCESS != rval ) return rval;
232 }
233 else
234 mbMergeTag = merge_tag;
235
236 std::vector< handle_id > handles( verts.size() );
237 int i = 0;
238 for( Range::iterator vit = verts.begin(); vit != verts.end(); ++vit )
239 {
240 handles[i].eh = *vit;
241 handles[i].val = vals[i];
242 i++;
243 }
244 // std::sort(handles.begin(), handles.end(), compare_handle_id);
245 qsort( &handles[0], handles.size(), sizeof( handle_id ), compare_handle_id );
246 i = 0;
247 while( i < (int)verts.size() - 1 )
248 {
249 handle_id first = handles[i];
250 int j = i + 1;
251 while( j < (int)verts.size() && handles[j].val == first.val )
252 {
253 rval = mbImpl->tag_set_data( mbMergeTag, &( handles[j].eh ), 1, &( first.eh ) );
254 if( rval != MB_SUCCESS ) return rval;
255 deadEnts.insert( handles[j].eh );
256 j++;
257 }
258 i = j;
259 }
260
261 rval = perform_merge( mbMergeTag );
262
263 return rval;
264 }
265 ErrorCode MergeMesh::find_merged_to( EntityHandle& tree_root, AdaptiveKDTree& tree, Tag merge_tag )
266 {
267 AdaptiveKDTreeIter iter;
268
269 // evaluate vertices in this leaf
270 Range leaf_range, leaf_range2;
271 std::vector< EntityHandle > sorted_leaves;
272 std::vector< double > coords;
273 std::vector< EntityHandle > merge_tag_val, leaves_out;
274
275 ErrorCode result = tree.get_tree_iterator( tree_root, iter );
276 if( MB_SUCCESS != result ) return result;
277 while( result == MB_SUCCESS )
278 {
279 sorted_leaves.push_back( iter.handle() );
280 result = iter.step();
281 }
282 if( result != MB_ENTITY_NOT_FOUND ) return result;
283 std::sort( sorted_leaves.begin(), sorted_leaves.end() );
284
285 std::vector< EntityHandle >::iterator it;
286 for( it = sorted_leaves.begin(); it != sorted_leaves.end(); ++it )
287 {
288
289 leaf_range.clear();
290 result = mbImpl->get_entities_by_handle( *it, leaf_range );
291 if( MB_SUCCESS != result ) return result;
292 coords.resize( 3 * leaf_range.size() );
293 merge_tag_val.resize( leaf_range.size() );
294 result = mbImpl->get_coords( leaf_range, &coords[0] );
295 if( MB_SUCCESS != result ) return result;
296 result = mbImpl->tag_get_data( merge_tag, leaf_range, &merge_tag_val[0] );
297 if( MB_SUCCESS != result ) return result;
298 Range::iterator rit;
299 unsigned int i;
300 bool inleaf_merged, outleaf_merged = false;
301 unsigned int lr_size = leaf_range.size();
302
303 for( i = 0, rit = leaf_range.begin(); i != lr_size; ++rit, i++ )
304 {
305 if( 0 != merge_tag_val[i] ) continue;
306 CartVect from( &coords[3 * i] );
307 inleaf_merged = false;
308
309 // check close-by leaves too
310 leaves_out.clear();
311 result =
312 tree.distance_search( from.array(), mergeTol, leaves_out, mergeTol, 1.0e-6, NULL, NULL, &tree_root );
313 leaf_range2.clear();
314 for( std::vector< EntityHandle >::iterator vit = leaves_out.begin(); vit != leaves_out.end(); ++vit )
315 {
316 if( *vit > *it )
317 { // if we haven't visited this leaf yet in the outer loop
318 result = mbImpl->get_entities_by_handle( *vit, leaf_range2, Interface::UNION );
319 if( MB_SUCCESS != result ) return result;
320 }
321 }
322 if( !leaf_range2.empty() )
323 {
324 coords.resize( 3 * ( lr_size + leaf_range2.size() ) );
325 merge_tag_val.resize( lr_size + leaf_range2.size() );
326 result = mbImpl->get_coords( leaf_range2, &coords[3 * lr_size] );
327 if( MB_SUCCESS != result ) return result;
328 result = mbImpl->tag_get_data( merge_tag, leaf_range2, &merge_tag_val[lr_size] );
329 if( MB_SUCCESS != result ) return result;
330 outleaf_merged = false;
331 }
332
333 // check other verts in this leaf
334 for( unsigned int j = i + 1; j < merge_tag_val.size(); j++ )
335 {
336 EntityHandle to_ent = j >= lr_size ? leaf_range2[j - lr_size] : leaf_range[j];
337
338 if( *rit == to_ent ) continue;
339
340 if( ( from - CartVect( &coords[3 * j] ) ).length_squared() < mergeTolSq )
341 {
342 merge_tag_val[j] = *rit;
343 if( j < lr_size )
344 {
345 inleaf_merged = true;
346 }
347 else
348 {
349 outleaf_merged = true;
350 }
351 deadEnts.insert( to_ent );
352 }
353 }
354 if( outleaf_merged )
355 {
356 result = mbImpl->tag_set_data( merge_tag, leaf_range2, &merge_tag_val[leaf_range.size()] );
357 if( MB_SUCCESS != result ) return result;
358 outleaf_merged = false;
359 }
360 if( inleaf_merged )
361 {
362 result = mbImpl->tag_set_data( merge_tag, leaf_range, &merge_tag_val[0] );
363 if( MB_SUCCESS != result ) return result;
364 }
365 }
366 }
367 return MB_SUCCESS;
368 }
369
370 // Determine which higher dimensional entities should be merged
371 ErrorCode MergeMesh::merge_higher_dimensions( Range& elems )
372 {
373 // apply a different strategy
374 // look at the vertices that were merged to, earlier, and find all entities adjacent to them
375 // elems (input) are used just for initial connectivity
376 ErrorCode result;
377 Range verts;
378 result = mbImpl->get_connectivity( elems, verts );
379 if( MB_SUCCESS != result ) return result;
380
381 // all higher dim entities that will be merged will be connected to the vertices that were
382 // merged earlier; we will look at these vertices only
383 Range mergedToVertsRange;
384 std::copy( mergedToVertices.rbegin(), mergedToVertices.rend(), range_inserter( mergedToVertsRange ) );
385 Range vertsOfInterest = intersect( mergedToVertsRange, verts );
386 // Go through each dimension
387 Range possibleEntsToMerge, conn, matches, moreDeadEnts;
388
389 for( int dim = 1; dim < 3; dim++ )
390 {
391 moreDeadEnts.clear();
392 possibleEntsToMerge.clear();
393 result = mbImpl->get_adjacencies( vertsOfInterest, dim, false, possibleEntsToMerge, Interface::UNION );
394 if( MB_SUCCESS != result ) return result;
395 // Go through each possible entity and see if it shares vertices with another entity of same
396 // dimension
397 for( Range::iterator pit = possibleEntsToMerge.begin(); pit != possibleEntsToMerge.end(); ++pit )
398 {
399 EntityHandle eh = *pit; // possible entity to be matched
400 conn.clear();
401 // Get the vertices connected to it in a range
402
403 result = mbImpl->get_connectivity( &eh, 1, conn );
404 if( MB_SUCCESS != result ) return result;
405 matches.clear();
406 // now retrieve all entities connected to all conn vertices
407 result = mbImpl->get_adjacencies( conn, dim, false, matches, Interface::INTERSECT );
408 if( MB_SUCCESS != result ) return result;
409 if( matches.size() > 1 )
410 {
411 for( Range::iterator matchIt = matches.begin(); matchIt != matches.end(); ++matchIt )
412 {
413 EntityHandle to_remove = *matchIt;
414 if( to_remove != eh )
415 {
416 moreDeadEnts.insert( to_remove );
417 result = mbImpl->merge_entities( eh, to_remove, false, false );
418 if( result != MB_SUCCESS ) return result;
419 possibleEntsToMerge.erase( to_remove );
420 }
421 }
422 }
423 }
424 // Delete the entities of dimension dim
425 result = mbImpl->delete_entities( moreDeadEnts );
426 if( result != MB_SUCCESS ) return result;
427 }
428 return MB_SUCCESS;
429 #if 0
430 Range skinEnts, adj, matches, moreDeadEnts;
431 ErrorCode result;
432 Skinner skinner(mbImpl);
433 //Go through each dimension
434 for (int dim = 1; dim < 3; dim++)
435 {
436 skinEnts.clear();
437 moreDeadEnts.clear();
438 result = skinner.find_skin(0, elems, dim, skinEnts, false, false);
439 //Go through each skin entity and see if it shares adjacancies with another entity
440 for (Range::iterator skinIt = skinEnts.begin();
441 skinIt != skinEnts.end(); ++skinIt)
442 {
443 adj.clear();
444 //Get the adjacencies 1 dimension lower
445 result = mbImpl->get_adjacencies(&(*skinIt), 1, dim - 1, false, adj);
446 if (result != MB_SUCCESS)
447 return result;
448 //See what other entities share these adjacencies
449 matches.clear();
450 result = mbImpl->get_adjacencies(adj, dim, false, matches,
451 Interface::INTERSECT);
452 if (result != MB_SUCCESS)
453 return result;
454 //If there is more than one entity, then we have some to merge and erase
455 if (matches.size() > 1)
456 {
457 for (Range::iterator matchIt = matches.begin();
458 matchIt != matches.end(); ++matchIt)
459 {
460 if (*matchIt != *skinIt)
461 {
462 moreDeadEnts.insert(*matchIt);
463 result = mbImpl->merge_entities(*skinIt, *matchIt, false, false);
464 if (result != MB_SUCCESS)
465 return result;
466 skinEnts.erase(*matchIt);
467 }
468 }
469 }
470 }
471 //Delete the entities
472 result = mbImpl->delete_entities(moreDeadEnts);
473 if (result != MB_SUCCESS)
474 return result;
475 }
476 return MB_SUCCESS;
477 #endif
478 }
479
480 } // End namespace moab
1.9.1.