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