MeshKit  1.0
MergeMesh.cpp
Go to the documentation of this file.
00001 #include "meshkit/MergeMesh.hpp"
00002 
00003 #include <algorithm>
00004 #include <cstdlib>
00005 #include <cstring>
00006 #include <functional>
00007 #include <string>
00008 #include <vector>
00009 #include <cassert>
00010 #include <iostream>
00011 
00012 #include "moab/Skinner.hpp"
00013 #include "moab/AdaptiveKDTree.hpp"
00014 #include "moab/Range.hpp"
00015 #include "moab/CartVect.hpp"
00016 
00017 namespace MeshKit
00018 {
00019   // static registration of this  mesh scheme
00020   moab::EntityType MergeMesh_tps[] = { moab::MBVERTEX,
00021                                        moab::MBEDGE,
00022                                        moab::MBTRI,
00023                                        moab::MBHEX,
00024                                        moab::MBMAXTYPE};
00025 
00026   const moab::EntityType* MergeMesh::output_types()
00027   { return MergeMesh_tps; }
00028   
00029   MergeMesh::MergeMesh(MKCore *mkcore, const MEntVector &me_vec)
00030     : MeshScheme(mkcore, me_vec),
00031       imeshImpl(mkcore->imesh_instance()), mbImpl (mkcore->moab_instance())
00032   {
00033     mergeTol = 1e-4;
00034     updateSets = 0;
00035     doMerge = 1;
00036     mergeTag = NULL;
00037   }
00038 
00039   MergeMesh::~MergeMesh()
00040   {
00041     if (mbMergeTag) mbImpl->tag_delete(mbMergeTag);
00042   }
00043 
00044   bool MergeMesh::add_modelent(ModelEnt *model_ent)
00045   {
00046     return MeshOp::add_modelent(model_ent);
00047   }
00048 
00049   void MergeMesh::set_merge_params(double merge_tol, int updatesets,
00050                                    int domerge, iBase_TagHandle merge_tag)
00051   {
00052     mergeTol = merge_tol;
00053     updateSets = updatesets;
00054     doMerge = domerge;
00055     mergeTag = merge_tag;
00056   }
00057 
00058   void MergeMesh::setup_this()
00059   {
00060     std::cout << "IN SETUPTHIS MERGEMESH " << std::endl;
00061 
00062   }
00063 
00064   void MergeMesh::execute_this()
00065   {
00066     std::cout << "IN EXTHIS MERGEMESH " << std::endl;
00067 
00068     std::vector<iBase_EntityHandle> ents(mentSelection.size());
00069     int dim = 3;
00070     for (MEntSelection::iterator mit = mentSelection.begin();
00071          mit != mentSelection.end(); mit++) {
00072       ModelEnt *me = mit->first;
00073       dim  = me->dimension();
00074     }
00075 
00076     if(dim == 3){
00077       imeshImpl->getEntities(0, iBase_REGION, iMesh_ALL_TOPOLOGIES,
00078                              ents);
00079     }
00080     else if (dim == 2){
00081       imeshImpl->getEntities(0, iBase_FACE, iMesh_ALL_TOPOLOGIES,
00082                              ents);
00083     }
00084     merge_entities(&ents[0], ents.size(), mergeTol, doMerge, updateSets,
00085                    mergeTag);
00086   }
00087 
00088   void MergeMesh::merge_entities(iBase_EntityHandle *elems,
00089                                  int elems_size,
00090                                  const double merge_tol,
00091                                  const int do_merge,
00092                                  const int update_sets,
00093                                  iBase_TagHandle merge_tag)  
00094   {
00095     mergeTol = merge_tol;
00096     mergeTolSq = merge_tol*merge_tol;
00097   
00098     moab::Range tmp_elems;
00099     tmp_elems.insert((moab::EntityHandle*)elems, (moab::EntityHandle*)elems + elems_size);
00100     moab::ErrorCode result = merge_entities(tmp_elems, do_merge, update_sets,
00101                                         (moab::Tag)merge_tag);
00102     if (moab::MB_SUCCESS != result)
00103       exit(0);
00104   }
00105 
00106   void MergeMesh::perform_merge(iBase_TagHandle merge_tag) 
00107   {
00108     // put into a range
00109     moab::ErrorCode result = perform_merge((moab::Tag) merge_tag);
00110     if (moab::MB_SUCCESS != result)
00111       exit(0);
00112   }
00113 
00114   moab::ErrorCode MergeMesh::merge_entities(moab::Range &elems,
00115                                         const int do_merge,
00116                                         const int update_sets,
00117                                         moab::Tag merge_tag) 
00118   {
00119     // get the skin of the entities
00120     moab::Skinner skinner(mbImpl);
00121     moab::Range skin_range;
00122     moab::ErrorCode result = skinner.find_skin(0, elems, 0, skin_range);
00123     if (moab::MB_SUCCESS != result) return result;
00124 
00125     // create a tag to mark merged-to entity; reuse tree_root
00126     moab::EntityHandle tree_root = 0;
00127     if (0 == merge_tag) {
00128       //result = mbImpl->tag_create("__merge_tag", sizeof(moab::EntityHandle), 
00129       //                            MB_TAG_DENSE, MB_TYPE_HANDLE,
00130       //                            mbMergeTag, &tree_root);
00131       result = mbImpl->tag_get_handle("__merge_tag", 1, moab::MB_TYPE_HANDLE,
00132                                       mbMergeTag, moab::MB_TAG_DENSE|moab::MB_TAG_CREAT,
00133                                       &tree_root);
00134       if (moab::MB_SUCCESS != result) return result;
00135     }
00136     else mbMergeTag = merge_tag;
00137   
00138     // build a kd tree with the vertices
00139     moab::AdaptiveKDTree kd(mbImpl);
00140     result = kd.build_tree(skin_range, &tree_root);
00141     if (moab::MB_SUCCESS != result) return result;
00142 
00143     // find matching vertices, mark them
00144     result = find_merged_to(kd, tree_root, mbMergeTag);
00145     if (moab::MB_SUCCESS != result) return result;
00146   
00147     // merge them if requested
00148     if (do_merge) {
00149       result = perform_merge(mbMergeTag);
00150       if (moab::MB_SUCCESS != result) return result;
00151     }
00152   
00153     return moab::MB_SUCCESS;
00154   }
00155 
00156   moab::ErrorCode MergeMesh::perform_merge(moab::Tag merge_tag) 
00157   {
00158     if (deadEnts.size()==0){
00159       std::cout << "\nWarning: Geometries don't have a common face; Nothing to merge" << std::endl;
00160       return moab::MB_SUCCESS; //nothing to merge carry on with the program
00161     }
00162     if (mbImpl->type_from_handle(*deadEnts.rbegin()) != moab::MBVERTEX)
00163       return moab::MB_FAILURE;
00164   
00165     std::vector<moab::EntityHandle> merge_tag_val(deadEnts.size());
00166     moab::ErrorCode result = mbImpl->tag_get_data(merge_tag, deadEnts, &merge_tag_val[0]);
00167     if (moab::MB_SUCCESS != result) return result;
00168   
00169     moab::Range::iterator rit;
00170     unsigned int i;
00171     for (rit = deadEnts.begin(), i = 0; rit != deadEnts.end(); rit++, i++) {
00172       assert(merge_tag_val[i]);
00173       result = mbImpl->merge_entities(merge_tag_val[i], *rit, false, false);
00174       if (moab::MB_SUCCESS != result) return result;
00175     }
00176   
00177     return mbImpl->delete_entities(deadEnts);
00178   }
00179 
00180   moab::ErrorCode MergeMesh::find_merged_to(moab::AdaptiveKDTree & tree, moab::EntityHandle &tree_root, moab::Tag merge_tag)
00181   {
00182     moab::AdaptiveKDTreeIter iter;
00183     //moab::AdaptiveKDTree tree(mbImpl);
00184   
00185     // evaluate vertices in this leaf
00186     moab::Range leaf_range, leaf_range2;
00187     std::vector<moab::EntityHandle> sorted_leaves;
00188     std::vector<double> coords;
00189     std::vector<moab::EntityHandle> merge_tag_val, leaves_out;
00190   
00191     moab::ErrorCode result = tree.get_tree_iterator(tree_root, iter);
00192     if (moab::MB_SUCCESS != result) return result;
00193     while (result == moab::MB_SUCCESS) {
00194       sorted_leaves.push_back( iter.handle() );
00195       result = iter.step();
00196     }
00197     if (result != moab::MB_ENTITY_NOT_FOUND)
00198       return result;
00199     std::sort( sorted_leaves.begin(), sorted_leaves.end() );
00200   
00201     std::vector<moab::EntityHandle>::iterator it;
00202     for (it = sorted_leaves.begin(); it != sorted_leaves.end(); ++it) {
00203 
00204       leaf_range.clear();
00205       result = mbImpl->get_entities_by_handle(*it, leaf_range);
00206       if (moab::MB_SUCCESS != result) return result;
00207       coords.resize(3*leaf_range.size());
00208       merge_tag_val.resize(leaf_range.size());
00209       result = mbImpl->get_coords(leaf_range, &coords[0]);
00210       if (moab::MB_SUCCESS != result) return result;
00211       result = mbImpl->tag_get_data(merge_tag, leaf_range, &merge_tag_val[0]);
00212       if (moab::MB_SUCCESS != result) return result;
00213       moab::Range::iterator rit;
00214       unsigned int i;
00215       bool inleaf_merged, outleaf_merged = false;
00216       unsigned int lr_size = leaf_range.size();
00217     
00218       for (i = 0, rit = leaf_range.begin(); i != lr_size; rit++, i++) {
00219         if (0 != merge_tag_val[i]) continue;
00220         moab::CartVect from(&coords[3*i]);
00221         inleaf_merged = false;
00222 
00223         // check close-by leaves too
00224         leaves_out.clear();
00225         result = tree.distance_search( from.array(), mergeTol,
00226                                              leaves_out);
00227         leaf_range2.clear();
00228         for (std::vector<moab::EntityHandle>::iterator vit = leaves_out.begin();
00229              vit != leaves_out.end(); vit++) {
00230           if (*vit > *it) { // if we haven't visited this leaf yet in the outer loop
00231             result = mbImpl->get_entities_by_handle(*vit, leaf_range2, moab::Interface::UNION);
00232             if (moab::MB_SUCCESS != result) return result;
00233           }
00234         }
00235         if (!leaf_range2.empty()) {
00236           coords.resize(3*(lr_size+leaf_range2.size()));
00237           merge_tag_val.resize(lr_size+leaf_range2.size());
00238           result = mbImpl->get_coords(leaf_range2, &coords[3*lr_size]);
00239           if (moab::MB_SUCCESS != result) return result;
00240           result = mbImpl->tag_get_data(merge_tag, leaf_range2, &merge_tag_val[lr_size]);
00241           if (moab::MB_SUCCESS != result) return result;
00242           outleaf_merged = false;
00243         }
00244 
00245         // check other verts in this leaf
00246         for (unsigned int j = i+1; j < merge_tag_val.size(); j++) {
00247           moab::EntityHandle to_ent = j >= lr_size ? leaf_range2[j-lr_size] : 
00248             leaf_range[j];
00249         
00250           if (*rit == to_ent) continue;
00251         
00252           if ((from - moab::CartVect(&coords[3*j])).length_squared() < mergeTolSq) {
00253             merge_tag_val[j] = *rit;
00254             if (j < lr_size){
00255               inleaf_merged = true;}
00256             else{
00257               outleaf_merged = true;}
00258 
00259             deadEnts.insert(to_ent);
00260           }
00261 
00262         }
00263         if (outleaf_merged) {
00264           result = mbImpl->tag_set_data(merge_tag, leaf_range2, &merge_tag_val[leaf_range.size()]);
00265           if (moab::MB_SUCCESS != result) return result;
00266           outleaf_merged = false;
00267         }
00268         if (inleaf_merged) {
00269           result = mbImpl->tag_set_data(merge_tag, leaf_range, &merge_tag_val[0]);
00270           if (moab::MB_SUCCESS != result) return result;
00271         }
00272 
00273       }
00274     }
00275     return moab::MB_SUCCESS;
00276   }
00277 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines