MOAB: Mesh Oriented datABase  (version 5.3.1)
MergeMesh.cpp
Go to the documentation of this file.
00001 #include "moab/MergeMesh.hpp"
00002 
00003 #include "moab/Skinner.hpp"
00004 #include "moab/AdaptiveKDTree.hpp"
00005 #include "moab/Range.hpp"
00006 #include "moab/CartVect.hpp"
00007 
00008 #include "Internals.hpp"
00009 #include <vector>
00010 #include <algorithm>
00011 #include <string>
00012 #include <vector>
00013 #include <cassert>
00014 #include <iostream>
00015 #include <iomanip>
00016 
00017 #include <cstdlib>
00018 
00019 namespace moab
00020 {
00021 
00022 MergeMesh::MergeMesh( Interface* impl, bool printErrorIn )
00023     : mbImpl( impl ), mbMergeTag( 0 ), mergeTol( 0.001 ), mergeTolSq( 0.000001 ), printError( printErrorIn )
00024 {
00025 }
00026 
00027 MergeMesh::~MergeMesh()
00028 {
00029     if( mbMergeTag ) mbImpl->tag_delete( mbMergeTag );
00030     mbMergeTag = NULL;
00031 }
00032 
00033 ErrorCode MergeMesh::merge_entities( EntityHandle* elems, int elems_size, const double merge_tol, const int do_merge,
00034                                      const int update_sets, Tag merge_tag, bool do_higher_dim )
00035 {
00036     mergeTol   = merge_tol;
00037     mergeTolSq = merge_tol * merge_tol;
00038     Range tmp_elems;
00039     tmp_elems.insert( elems, elems + elems_size );
00040     ErrorCode result = merge_entities( tmp_elems, merge_tol, do_merge, update_sets, (Tag)merge_tag, do_higher_dim );
00041 
00042     return result;
00043 }
00044 
00045 /*  This function appears to be not necessary after MOAB conversion
00046 
00047  void MergeMesh::perform_merge(iBase_TagHandle merge_tag)
00048  {
00049  // put into a range
00050  ErrorCode result = perform_merge((Tag) merge_tag);
00051  if (result != MB_SUCCESS)
00052  throw MKException(iBase_FAILURE, "");
00053  }*/
00054 
00055 ErrorCode MergeMesh::merge_entities( Range& elems, const double merge_tol, const int do_merge, const int, Tag merge_tag,
00056                                      bool merge_higher_dim )
00057 {
00058     // If merge_higher_dim is true, do_merge must also be true
00059     if( merge_higher_dim && !do_merge ) { return MB_FAILURE; }
00060 
00061     mergeTol   = merge_tol;
00062     mergeTolSq = merge_tol * merge_tol;
00063 
00064     // get the skin of the entities
00065     Skinner skinner( mbImpl );
00066     Range skin_range;
00067     ErrorCode result = skinner.find_skin( 0, elems, 0, skin_range, false, false );
00068     if( MB_SUCCESS != result ) return result;
00069 
00070     // create a tag to mark merged-to entity; reuse tree_root
00071     EntityHandle tree_root = 0;
00072     if( 0 == merge_tag )
00073     {
00074         result = mbImpl->tag_get_handle( "__merge_tag", 1, MB_TYPE_HANDLE, mbMergeTag, MB_TAG_DENSE | MB_TAG_EXCL,
00075                                          &tree_root );
00076         if( MB_SUCCESS != result ) return result;
00077     }
00078     else
00079         mbMergeTag = merge_tag;
00080 
00081     // build a kd tree with the vertices
00082     AdaptiveKDTree kd( mbImpl );
00083     result = kd.build_tree( skin_range, &tree_root );
00084     if( MB_SUCCESS != result ) return result;
00085 
00086     // find matching vertices, mark them
00087     result = find_merged_to( tree_root, kd, mbMergeTag );
00088     if( MB_SUCCESS != result ) return result;
00089 
00090     // merge them if requested
00091     if( do_merge )
00092     {
00093         result = perform_merge( mbMergeTag );
00094         if( MB_SUCCESS != result ) return result;
00095     }
00096 
00097     if( merge_higher_dim && deadEnts.size() != 0 )
00098     {
00099         result = merge_higher_dimensions( elems );
00100         if( MB_SUCCESS != result ) return result;
00101     }
00102 
00103     return MB_SUCCESS;
00104 }
00105 
00106 ErrorCode MergeMesh::merge_all( EntityHandle meshset, const double merge_tol )
00107 {
00108     ErrorCode rval;
00109     if( 0 == mbMergeTag )
00110     {
00111         EntityHandle def_val = 0;
00112         rval = mbImpl->tag_get_handle( "__merge_tag", 1, MB_TYPE_HANDLE, mbMergeTag, MB_TAG_DENSE | MB_TAG_EXCL,
00113                                        &def_val );MB_CHK_ERR( rval );
00114     }
00115     // get all entities;
00116     // get all vertices connected
00117     // build a kdtree
00118     // find merged to
00119     mergeTol   = merge_tol;
00120     mergeTolSq = merge_tol * merge_tol;
00121 
00122     // get all vertices
00123     Range entities;
00124     rval = mbImpl->get_entities_by_handle( meshset, entities, /*recursive*/ true );MB_CHK_ERR( rval );
00125     Range sets = entities.subset_by_type( MBENTITYSET );
00126     entities   = subtract( entities, sets );
00127     Range verts;
00128     rval = mbImpl->get_connectivity( entities, verts );MB_CHK_ERR( rval );
00129 
00130     // build a kd tree with the vertices
00131     AdaptiveKDTree kd( mbImpl );
00132     EntityHandle tree_root;
00133     rval = kd.build_tree( verts, &tree_root );MB_CHK_ERR( rval );
00134     // find matching vertices, mark them
00135     rval = find_merged_to( tree_root, kd, mbMergeTag );MB_CHK_ERR( rval );
00136 
00137     rval = perform_merge( mbMergeTag );MB_CHK_ERR( rval );
00138 
00139     if( deadEnts.size() != 0 )
00140     {
00141         rval = merge_higher_dimensions( entities );MB_CHK_ERR( rval );
00142     }
00143     return MB_SUCCESS;
00144 }
00145 ErrorCode MergeMesh::perform_merge( Tag merge_tag )
00146 {
00147     // we start with an empty range of vertices that are "merged to"
00148     // they are used (eventually) for higher dim entities
00149     mergedToVertices.clear();
00150     ErrorCode result;
00151     if( deadEnts.size() == 0 )
00152     {
00153         if( printError ) std::cout << "\nWarning: Geometries don't have a common face; Nothing to merge" << std::endl;
00154         return MB_SUCCESS;  // nothing to merge carry on with the program
00155     }
00156     if( mbImpl->type_from_handle( *deadEnts.begin() ) != MBVERTEX ) return MB_FAILURE;
00157     std::vector< EntityHandle > merge_tag_val( deadEnts.size() );
00158     Range deadEntsRange;
00159     std::copy( deadEnts.rbegin(), deadEnts.rend(), range_inserter( deadEntsRange ) );
00160     result = mbImpl->tag_get_data( merge_tag, deadEntsRange, &merge_tag_val[0] );
00161     if( MB_SUCCESS != result ) return result;
00162 
00163     std::set< EntityHandle >::iterator rit;
00164     unsigned int i;
00165     for( rit = deadEnts.begin(), i = 0; rit != deadEnts.end(); ++rit, i++ )
00166     {
00167         assert( merge_tag_val[i] );
00168         if( MBVERTEX == TYPE_FROM_HANDLE( merge_tag_val[i] ) ) mergedToVertices.insert( merge_tag_val[i] );
00169         result = mbImpl->merge_entities( merge_tag_val[i], *rit, false, false );
00170         if( MB_SUCCESS != result ) { return result; }
00171     }
00172     result = mbImpl->delete_entities( deadEntsRange );
00173     return result;
00174 }
00175 // merge vertices according to an input tag
00176 // merge them if the tags are equal
00177 struct handle_id
00178 {
00179     EntityHandle eh;
00180     int val;
00181 };
00182 
00183 // handle structure comparison function for qsort
00184 // if the id is the same , compare the handle.
00185 int compare_handle_id( const void* a, const void* b )
00186 {
00187 
00188     handle_id* ia = (handle_id*)a;
00189     handle_id* ib = (handle_id*)b;
00190     if( ia->val == ib->val ) { return ( ia->eh < ib->eh ) ? -1 : 1; }
00191     else
00192     {
00193         return ( ia->val - ib->val );
00194     }
00195 }
00196 
00197 ErrorCode MergeMesh::merge_using_integer_tag( Range& verts, Tag user_tag, Tag merge_tag )
00198 {
00199     ErrorCode rval;
00200     DataType tag_type;
00201     rval = mbImpl->tag_get_data_type( user_tag, tag_type );
00202     if( rval != MB_SUCCESS || tag_type != MB_TYPE_INTEGER ) return MB_FAILURE;
00203 
00204     std::vector< int > vals( verts.size() );
00205     rval = mbImpl->tag_get_data( user_tag, verts, &vals[0] );
00206     if( rval != MB_SUCCESS ) return rval;
00207 
00208     if( 0 == merge_tag )
00209     {
00210         EntityHandle def_val = 0;
00211         rval = mbImpl->tag_get_handle( "__merge_tag", 1, MB_TYPE_HANDLE, mbMergeTag, MB_TAG_DENSE | MB_TAG_EXCL,
00212                                        &def_val );
00213         if( MB_SUCCESS != rval ) return rval;
00214     }
00215     else
00216         mbMergeTag = merge_tag;
00217 
00218     std::vector< handle_id > handles( verts.size() );
00219     int i = 0;
00220     for( Range::iterator vit = verts.begin(); vit != verts.end(); ++vit )
00221     {
00222         handles[i].eh  = *vit;
00223         handles[i].val = vals[i];
00224         i++;
00225     }
00226     // std::sort(handles.begin(), handles.end(), compare_handle_id);
00227     qsort( &handles[0], handles.size(), sizeof( handle_id ), compare_handle_id );
00228     i = 0;
00229     while( i < (int)verts.size() - 1 )
00230     {
00231         handle_id first = handles[i];
00232         int j           = i + 1;
00233         while( j < (int)verts.size() && handles[j].val == first.val  )
00234         {
00235             rval = mbImpl->tag_set_data( mbMergeTag, &( handles[j].eh ), 1, &( first.eh ) );
00236             if( rval != MB_SUCCESS ) return rval;
00237             deadEnts.insert( handles[j].eh );
00238             j++;
00239         }
00240         i = j;
00241     }
00242 
00243     rval = perform_merge( mbMergeTag );
00244 
00245     return rval;
00246 }
00247 ErrorCode MergeMesh::find_merged_to( EntityHandle& tree_root, AdaptiveKDTree& tree, Tag merge_tag )
00248 {
00249     AdaptiveKDTreeIter iter;
00250 
00251     // evaluate vertices in this leaf
00252     Range leaf_range, leaf_range2;
00253     std::vector< EntityHandle > sorted_leaves;
00254     std::vector< double > coords;
00255     std::vector< EntityHandle > merge_tag_val, leaves_out;
00256 
00257     ErrorCode result = tree.get_tree_iterator( tree_root, iter );
00258     if( MB_SUCCESS != result ) return result;
00259     while( result == MB_SUCCESS )
00260     {
00261         sorted_leaves.push_back( iter.handle() );
00262         result = iter.step();
00263     }
00264     if( result != MB_ENTITY_NOT_FOUND ) return result;
00265     std::sort( sorted_leaves.begin(), sorted_leaves.end() );
00266 
00267     std::vector< EntityHandle >::iterator it;
00268     for( it = sorted_leaves.begin(); it != sorted_leaves.end(); ++it )
00269     {
00270 
00271         leaf_range.clear();
00272         result = mbImpl->get_entities_by_handle( *it, leaf_range );
00273         if( MB_SUCCESS != result ) return result;
00274         coords.resize( 3 * leaf_range.size() );
00275         merge_tag_val.resize( leaf_range.size() );
00276         result = mbImpl->get_coords( leaf_range, &coords[0] );
00277         if( MB_SUCCESS != result ) return result;
00278         result = mbImpl->tag_get_data( merge_tag, leaf_range, &merge_tag_val[0] );
00279         if( MB_SUCCESS != result ) return result;
00280         Range::iterator rit;
00281         unsigned int i;
00282         bool inleaf_merged, outleaf_merged = false;
00283         unsigned int lr_size = leaf_range.size();
00284 
00285         for( i = 0, rit = leaf_range.begin(); i != lr_size; ++rit, i++ )
00286         {
00287             if( 0 != merge_tag_val[i] ) continue;
00288             CartVect from( &coords[3 * i] );
00289             inleaf_merged = false;
00290 
00291             // check close-by leaves too
00292             leaves_out.clear();
00293             result =
00294                 tree.distance_search( from.array(), mergeTol, leaves_out, mergeTol, 1.0e-6, NULL, NULL, &tree_root );
00295             leaf_range2.clear();
00296             for( std::vector< EntityHandle >::iterator vit = leaves_out.begin(); vit != leaves_out.end(); ++vit )
00297             {
00298                 if( *vit > *it )
00299                 {  // if we haven't visited this leaf yet in the outer loop
00300                     result = mbImpl->get_entities_by_handle( *vit, leaf_range2, Interface::UNION );
00301                     if( MB_SUCCESS != result ) return result;
00302                 }
00303             }
00304             if( !leaf_range2.empty() )
00305             {
00306                 coords.resize( 3 * ( lr_size + leaf_range2.size() ) );
00307                 merge_tag_val.resize( lr_size + leaf_range2.size() );
00308                 result = mbImpl->get_coords( leaf_range2, &coords[3 * lr_size] );
00309                 if( MB_SUCCESS != result ) return result;
00310                 result = mbImpl->tag_get_data( merge_tag, leaf_range2, &merge_tag_val[lr_size] );
00311                 if( MB_SUCCESS != result ) return result;
00312                 outleaf_merged = false;
00313             }
00314 
00315             // check other verts in this leaf
00316             for( unsigned int j = i + 1; j < merge_tag_val.size(); j++ )
00317             {
00318                 EntityHandle to_ent = j >= lr_size ? leaf_range2[j - lr_size] : leaf_range[j];
00319 
00320                 if( *rit == to_ent ) continue;
00321 
00322                 if( ( from - CartVect( &coords[3 * j] ) ).length_squared() < mergeTolSq )
00323                 {
00324                     merge_tag_val[j] = *rit;
00325                     if( j < lr_size ) { inleaf_merged = true; }
00326                     else
00327                     {
00328                         outleaf_merged = true;
00329                     }
00330                     deadEnts.insert( to_ent );
00331                 }
00332             }
00333             if( outleaf_merged )
00334             {
00335                 result = mbImpl->tag_set_data( merge_tag, leaf_range2, &merge_tag_val[leaf_range.size()] );
00336                 if( MB_SUCCESS != result ) return result;
00337                 outleaf_merged = false;
00338             }
00339             if( inleaf_merged )
00340             {
00341                 result = mbImpl->tag_set_data( merge_tag, leaf_range, &merge_tag_val[0] );
00342                 if( MB_SUCCESS != result ) return result;
00343             }
00344         }
00345     }
00346     return MB_SUCCESS;
00347 }
00348 
00349 // Determine which higher dimensional entities should be merged
00350 ErrorCode MergeMesh::merge_higher_dimensions( Range& elems )
00351 {
00352     // apply a different strategy
00353     // look at the vertices that were merged to, earlier, and find all entities adjacent to them
00354     // elems (input) are used just for initial connectivity
00355     ErrorCode result;
00356     Range verts;
00357     result = mbImpl->get_connectivity( elems, verts );
00358     if( MB_SUCCESS != result ) return result;
00359 
00360     // all higher dim entities that will be merged will be connected to the vertices that were
00361     // merged earlier; we will look at these vertices only
00362     Range mergedToVertsRange;
00363     std::copy( mergedToVertices.rbegin(), mergedToVertices.rend(), range_inserter( mergedToVertsRange ) );
00364     Range vertsOfInterest = intersect( mergedToVertsRange, verts );
00365     // Go through each dimension
00366     Range possibleEntsToMerge, conn, matches, moreDeadEnts;
00367 
00368     for( int dim = 1; dim < 3; dim++ )
00369     {
00370         moreDeadEnts.clear();
00371         possibleEntsToMerge.clear();
00372         result = mbImpl->get_adjacencies( vertsOfInterest, dim, false, possibleEntsToMerge, Interface::UNION );
00373         if( MB_SUCCESS != result ) return result;
00374         // Go through each possible entity and see if it shares vertices with another entity of same
00375         // dimension
00376         for( Range::iterator pit = possibleEntsToMerge.begin(); pit != possibleEntsToMerge.end(); ++pit )
00377         {
00378             EntityHandle eh = *pit;  // possible entity to be matched
00379             conn.clear();
00380             // Get the vertices connected to it in a range
00381 
00382             result = mbImpl->get_connectivity( &eh, 1, conn );
00383             if( MB_SUCCESS != result ) return result;
00384             matches.clear();
00385             // now retrieve all entities connected to all conn vertices
00386             result = mbImpl->get_adjacencies( conn, dim, false, matches, Interface::INTERSECT );
00387             if( MB_SUCCESS != result ) return result;
00388             if( matches.size() > 1 )
00389             {
00390                 for( Range::iterator matchIt = matches.begin(); matchIt != matches.end(); ++matchIt )
00391                 {
00392                     EntityHandle to_remove = *matchIt;
00393                     if( to_remove != eh )
00394                     {
00395                         moreDeadEnts.insert( to_remove );
00396                         result = mbImpl->merge_entities( eh, to_remove, false, false );
00397                         if( result != MB_SUCCESS ) return result;
00398                         possibleEntsToMerge.erase( to_remove );
00399                     }
00400                 }
00401             }
00402         }
00403         // Delete the entities of dimension dim
00404         result = mbImpl->delete_entities( moreDeadEnts );
00405         if( result != MB_SUCCESS ) return result;
00406     }
00407     return MB_SUCCESS;
00408 #if 0
00409   Range skinEnts, adj, matches, moreDeadEnts;
00410   ErrorCode result;
00411   Skinner skinner(mbImpl);
00412   //Go through each dimension
00413   for (int dim = 1; dim < 3; dim++)
00414   {
00415     skinEnts.clear();
00416     moreDeadEnts.clear();
00417     result = skinner.find_skin(0, elems, dim, skinEnts, false, false);
00418     //Go through each skin entity and see if it shares adjacancies with another entity
00419     for (Range::iterator skinIt = skinEnts.begin();
00420         skinIt != skinEnts.end(); ++skinIt)
00421     {
00422       adj.clear();
00423       //Get the adjacencies 1 dimension lower
00424       result = mbImpl->get_adjacencies(&(*skinIt), 1, dim - 1, false, adj);
00425       if (result != MB_SUCCESS)
00426         return result;
00427       //See what other entities share these adjacencies
00428       matches.clear();
00429       result = mbImpl->get_adjacencies(adj, dim, false, matches,
00430           Interface::INTERSECT);
00431       if (result != MB_SUCCESS)
00432         return result;
00433       //If there is more than one entity, then we have some to merge and erase
00434       if (matches.size() > 1)
00435       {
00436         for (Range::iterator matchIt = matches.begin();
00437             matchIt != matches.end(); ++matchIt)
00438         {
00439           if (*matchIt != *skinIt)
00440           {
00441             moreDeadEnts.insert(*matchIt);
00442             result = mbImpl->merge_entities(*skinIt, *matchIt, false, false);
00443             if (result != MB_SUCCESS)
00444               return result;
00445             skinEnts.erase(*matchIt);
00446           }
00447         }
00448       }
00449     }
00450     //Delete the entities
00451     result = mbImpl->delete_entities(moreDeadEnts);
00452     if (result != MB_SUCCESS)
00453       return result;
00454   }
00455   return MB_SUCCESS;
00456 #endif
00457 }
00458 
00459 }  // End namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines