MeshKit  1.0
MBVolOp.cpp
Go to the documentation of this file.
00001 /*
00002  * MBVolOp.cpp
00003  *
00004  *  Created on: Jan 13, 2012
00005  */
00006 
00007 #include "meshkit/MBVolOp.hpp"
00008 #include "meshkit/MBSplitOp.hpp"
00009 #include "meshkit/ModelEnt.hpp"
00010 #include "moab/CartVect.hpp"
00011 
00012 //#include "moab/Core.hpp"
00013 //#include "meshkit/FBiGeom.hpp"
00014 
00015 
00016 namespace MeshKit {
00017 
00018 
00019 //Entity Type initialization for splitting; no mesh output
00020 moab::EntityType MBVolOp_tps[] = { moab::MBMAXTYPE }; // no mesh, really
00021 const moab::EntityType* MBVolOp::output_types()
00022 {
00023   return MBVolOp_tps;
00024 }
00025 
00026 MBVolOp::MBVolOp(MKCore *mk_core, const MEntVector &me_vec) :
00027   MeshScheme(mk_core, me_vec)
00028 {
00029   _direction[0]=_direction[1]=0.;
00030   _direction[2]=1.0;
00031   _pGTT = NULL;
00032   _rootSet = 0; // it means no set yet, although 0 mean root set in moab
00033   _fbe = NULL;
00034 }
00035 
00036 MBVolOp::~MBVolOp() {
00037   // TODO Auto-generated destructor stub
00038 }
00039 void MBVolOp::setup_this()
00040 {
00041   // construct an FBEngine object, used more like a container, and to trigger the weaving
00042   // it is not really required; this object is based on gtt object build from top and bottom faces
00043   // ( which are extracted from model entities of dimension 2)
00044   // so, involve FBEngine just because we have something we need there
00045   // collect the top and bottom faces from ment vector
00046   int nTotSurf = this->mentSelection.size();
00047   std::cout << " total number of faces:" << nTotSurf << "\n";
00048 
00049   // grab all surfaces, and duplicate model using gtt, to get new gsets, that will be
00050   // continued with volume sets in the end
00051   // establish the loops from faces
00052   MEntSelection::iterator mit;
00053 
00054   moab::Interface * MBI = mk_core()->moab_instance();
00055   moab::GeomTopoTool gtt(MBI, true);// to retrieve the gsets
00056   moab::EntityHandle mset;
00057   // these should all be of dimension 2, faces
00058   std::vector<moab::EntityHandle> vSurfaces;
00059 
00060   moab::ErrorCode rval;
00061 
00062   for (mit = mentSelection.begin(); mit != mentSelection.end(); mit++) {
00063     mset = (mit->first)->mesh_handle();
00064     // get the globalId tag
00065     vSurfaces.push_back(mset);
00066   }
00067 
00068   rval = gtt.duplicate_model(_pGTT, &vSurfaces);
00069   MBERRCHK(rval, MBI);
00070   // now, _pGTT will have the new sets, which will form the basis for the volume creation
00071   // new gsets will be added to this _pGTT, and this will (unfortunately) create an OBB tree too,
00072   // although we do not really need it (actually, 2 obb trees that are not needed!) bummer!
00073   // the result will be the model root set of the new GTT, stored in result range
00074   // that corresponds to the first model ent
00075   ModelEnt * firstMe = (*(mentSelection.begin()) ).first;
00076   moab::Range & resultRange = mentSelection[firstMe];
00077   _rootSet = _pGTT->get_root_model_set();
00078   resultRange.insert(_rootSet);
00079   _fbe = new moab::FBEngine(MBI, _pGTT, false);// smooth=false, not needed here (are you sure?)
00080   //we will use FBEngine for querying; although ModelEnt would have helped
00081   // still, the model ent works with geometry adjacency; in this case, there is no geometry
00082   // we need to work directly with MOAB;
00083   // not pretty :(
00084   // we build all this infrastructure and I am not using it
00085 
00086 
00087   establish_mapping();
00088   return;
00089 }
00090 void MBVolOp::establish_mapping()
00091 {
00092   // the first half of the surfaces are bottom, next are top
00093   // establish the connection between top and bottom entities
00094   // use the direction for vertices, and tangent direction for edges
00095 
00096   // these are all vertices from current gtt
00097   moab::Interface * MBI = mk_core()->moab_instance();
00098   moab::Range verticeSets= _pGTT->geoRanges()[0];
00099   int num_vertices = verticeSets.size();
00100   std::vector<moab::CartVect> coordVert;
00101   coordVert.resize(num_vertices);
00102   std::vector<int> corrV;
00103   corrV.resize(num_vertices);
00104   moab::ErrorCode rval = moab::MB_SUCCESS;
00105   std::map<moab::EntityHandle, int> indexVertex;
00106   for (int i=0; i<num_vertices; i++)
00107   {
00108     rval = _fbe->getVtxCoord(verticeSets[i], &(coordVert[i][0]),
00109         &(coordVert[i][1]), &(coordVert[i][2]));
00110     MBERRCHK(rval, MBI);
00111     corrV[i]=-1;// no correspondence yet
00112     indexVertex[verticeSets[i]] = i;
00113   }
00114 
00115   // decide with an expensive linear search, what vertices correspond to what
00116   moab::CartVect dirr(_direction);
00117   dirr.normalize();
00118   int j=-1;
00119   for ( j=0; j<num_vertices; j++)
00120   {
00121     if (corrV[j]!=-1)
00122       continue; // we already know about this one
00123     moab::CartVect & v1 = coordVert[j];
00124     int minIndex = -1;
00125     double minVal = 1.e38; // HUGE
00126     for (int k=0; k<num_vertices; k++)
00127     {
00128       if (j==k || corrV[k]>-1)
00129         continue;
00130       moab::CartVect product = (v1-coordVert[k])*dirr;
00131       double valProd = product.length_squared();
00132       if (valProd<minVal)
00133       {
00134         minVal = valProd;
00135         minIndex=k;
00136       }
00137     }
00138     std::cout<<"j: "<< j <<" val min:" << minVal << " min index:" << minIndex << "\n";
00139     if (minIndex==-1)
00140     {
00141       std::cout<< "Error Stop it. j: " << j << "\n";
00142       continue;
00143     }
00144     // make sure that the lower index is on bottom
00145     double dotProd = (coordVert[minIndex]-coordVert[j])%dirr;
00146     if (dotProd < 0)
00147     {
00148       std::cout<<"wrong orientation, bottom vertices should be of lower index\n";
00149       MBERRCHK(moab::MB_FAILURE, MBI);// bail out from here, stop it!
00150     }
00151 
00152     corrV[j] = minIndex;
00153     corrV[minIndex] = j;
00154   }
00155   // check that we do not have any left vertices
00156   for (j=0; j<num_vertices; j++)
00157   {
00158     if (corrV[j]==-1)
00159     {
00160       std::cout<<"Error in finding a corresponding vertex to j:"<<j <<"\n";
00161       MBERRCHK(moab::MB_FAILURE, MBI);// bail out from here, stop it!
00162     }
00163     vertexMap[verticeSets[j]]= verticeSets[corrV[j]];
00164   }
00165 
00166   // so now we have mappings between vertices
00167   // we need to build mappings between edges and faces.
00168   // we usually start with bottom and continue to the top
00169   moab::Range edgeSets= _pGTT->geoRanges()[1];
00170   int numEdges = edgeSets.size();
00171   std::vector<moab::CartVect> edgeTangents;
00172   edgeTangents.resize(numEdges*2);// we need 2 tangents per edge
00173   // for each vertex, compute the tangents at ends, projected on a plan perpendicular to the direction
00174   // (usually, xy plan...)
00175   // tangents at both ends!// we know that the u range is from 0 to 1
00176   std::vector<int> corrE;
00177   corrE.resize(numEdges);
00178   std::map<moab::EntityHandle, int> indexEdge;
00179   for (j=0; j<numEdges; j++)
00180   {
00181     corrE[j] = -1;// no correspondence yet
00182     indexEdge[edgeSets[j]] = j;
00183     // careful about the FBEngine, it is not smooth
00184     std::vector<moab::EntityHandle> meshEdges;
00185     rval = MBI->get_entities_by_type(edgeSets[j], moab::MBEDGE,  meshEdges);
00186     MBERRCHK(rval, MBI);
00187     // get first and last edges, connectivity and compute tangent
00188     // it will be used to map edges
00189     // first edge, last edge
00190     moab::EntityHandle firstMeshEdge = meshEdges[0];
00191     moab::EntityHandle lastMeshEdge = meshEdges[meshEdges.size()-1];
00192     const moab::EntityHandle  * conn=NULL;
00193     int nnodes;
00194     rval = MBI->get_connectivity(firstMeshEdge, conn, nnodes);
00195     MBERRCHK(rval, MBI);
00196     moab::CartVect posi[2];
00197     rval = MBI->get_coords(conn, 2, &(posi[0][0]));
00198     MBERRCHK(rval, MBI);
00199     posi[0] = posi[1]-posi[0]; //
00200     posi[0] = posi[0]*dirr;
00201     edgeTangents[2*j] = posi[0]*dirr;// this is in the plane
00202     edgeTangents[2*j].normalize();// it should be non null, but accidents happen :)
00203     rval = MBI->get_connectivity(lastMeshEdge, conn, nnodes);
00204     MBERRCHK(rval, MBI);
00205     rval = MBI->get_coords(conn, 2, &(posi[0][0]));
00206     MBERRCHK(rval, MBI);
00207     posi[0] = posi[1]-posi[0]; //
00208     posi[0] = posi[0]*dirr;
00209     edgeTangents[2*j+1] = posi[0]*dirr;// this is in the plane
00210     edgeTangents[2*j+1].normalize();
00211   }
00212   // now try to match edges based on their vertices matching, and start and end tangents matching
00213   for (j = 0; j<numEdges; j++)
00214   {
00215     if (corrE[j]>=0)
00216       continue; // we already have a correspondent edge for it
00217     moab::Range adjVertices;
00218     rval = _fbe-> getEntAdj(edgeSets[j], /*vertex type*/0, adjVertices);
00219     MBERRCHK(rval, MBI);
00220     if (adjVertices.size()==2)
00221     {
00222       int sense=0;
00223       rval = _fbe->getEgVtxSense(edgeSets[j], adjVertices[0], adjVertices[1], sense);
00224       MBERRCHK(rval, MBI);
00225       // find the edges that are adjacent to map(v0) and map(v1)
00226       moab::EntityHandle mapV0 = vertexMap [adjVertices[0]];
00227       moab::EntityHandle mapV1 = vertexMap [adjVertices[1]];
00228       // get all edges adjacent to both vertices
00229       moab::Range adjEdges0, adjEdges1;
00230       rval = _fbe-> getEntAdj(mapV0, /*edge type*/1, adjEdges0);
00231       MBERRCHK(rval, MBI);
00232       rval = _fbe-> getEntAdj(mapV1, /*edge type*/1, adjEdges1);
00233       adjEdges0 = intersect(adjEdges0, adjEdges1);
00234       // the mapped edge should be in this range
00235       moab::CartVect & t0=edgeTangents[2*j];
00236       moab::CartVect & t1=edgeTangents[2*j+1];
00237       for (moab::Range::iterator rit= adjEdges0.begin(); rit!=adjEdges0.end(); rit++)
00238       {
00239         moab::EntityHandle candidateMapEdge = *rit;
00240         int sense1=0;
00241         int indexCandEdge = indexEdge[candidateMapEdge];
00242         if (indexCandEdge == j)
00243           // error
00244           MBERRCHK(moab::MB_FAILURE, MBI);
00245         moab::CartVect & tm0=edgeTangents[2*indexCandEdge];
00246         moab::CartVect & tm1=edgeTangents[2*indexCandEdge+1];
00247         rval = _fbe->getEgVtxSense(candidateMapEdge, mapV0, mapV1, sense1);
00248         MBERRCHK(rval, MBI);
00249 
00250         if (sense*sense1>0)// same sense
00251         {
00252           if ( (t0%tm0 >= 0.99 ) && (t1%tm1>=0.99))// same sense, almost..
00253           {
00254             corrE[j]=indexCandEdge;
00255             corrE[indexCandEdge] = j;
00256           }
00257         }
00258         else
00259         {
00260           // different senses, check the opposite tangents
00261           if ( (t0%tm0 <=-0.99) && (t1%tm1<=-0.99) )
00262           {
00263             // how do we store the opposite senses?
00264             corrE[j]=indexCandEdge;
00265             corrE[indexCandEdge] = j;
00266           }
00267         }
00268       }
00269       if (corrE[j]==-1)
00270         MBERRCHK(moab::MB_FAILURE, MBI);// can't find a corresponding edge
00271 
00272     }
00273   }
00274   // so we have edge matching for every edge
00275   for (j=0; j<numEdges; j++)
00276   {
00277     if (corrE[j]==-1)
00278     {
00279       std::cout<<"Error in finding a corresponding edge: "<<j <<"\n";
00280       MBERRCHK(moab::MB_FAILURE, MBI);// bail out from here, stop it!
00281     }
00282     std::cout<< "edge j=" << j << "  mapped to edge " << corrE[j] << "\n";
00283     edgeMap[edgeSets[j]]= edgeSets[corrE[j]];
00284   }
00285   // now, we still need to separate top and bottom faces, somehow
00286   // we assume faces were stenciled correctly with the same polylines
00287   // start from bottom vertices (first half of the vertices is on bottom)
00288   moab::Range bottomFaces;
00289   for (j=0; j<num_vertices/2; j++)// we know that half are on bottom, we verified that already
00290   {
00291     moab::Range faces;
00292     rval = _fbe-> getEntAdj(verticeSets[j], /*face type*/2, faces);
00293     MBERRCHK(rval, MBI);
00294     bottomFaces.merge(faces);
00295   }
00296   // now find the mapped face for each of these, based on mapped edges
00297   for (j=0; j<(int)bottomFaces.size(); j++)
00298   {
00299     moab::Range edges;
00300     rval = _fbe-> getEntAdj(bottomFaces[j], /*edge type*/1, edges);
00301     MBERRCHK(rval, MBI);
00302     moab::Range mappedEdges;
00303     // get now the face connected to the mapped edges
00304     int i=0;
00305     // find the face with all those edges among the other faces
00306 
00307     moab::Range mapFaces;
00308     rval = _fbe-> getEntAdj(edgeMap[edges[0]], /*edge type*/2, mapFaces);
00309     MBERRCHK(rval, MBI);
00310     for (i=1; i<(int)edges.size(); i++)
00311     {
00312       moab::Range faces;
00313       rval = _fbe-> getEntAdj(edgeMap[edges[i]], /*edge type*/2, faces);
00314       mapFaces = intersect(mapFaces, faces);
00315     }
00316     if (mapFaces.size()!=1)
00317     {
00318       std::cout<<"Can't find unique mapped face to face index " << j << "\n";
00319       MBERRCHK(moab::MB_FAILURE, MBI);// bail out from here, stop it!
00320     }
00321     faceMap[bottomFaces[j]] = mapFaces[0];
00322     faceMap[mapFaces[0]]=bottomFaces[j];
00323     std::cout << "face j:" << j << " set:" << MBI->id_from_handle(bottomFaces[j]) << " mapped to "
00324         << MBI->id_from_handle(mapFaces[0]) << "\n";
00325   }
00326 
00327 }
00328 
00329 void MBVolOp::execute_this()
00330 {
00331   // now, we have the maps between top and bottom faces established, also for edges and vertices
00332   // first build edges between corresponding vertices, then faces between edges, then
00333   // finally, volumes
00334   // we start from bottom towards top
00335   moab::ErrorCode rval = moab::MB_SUCCESS;
00336   moab::Interface * MBI = mk_core()->moab_instance();
00337 
00338   moab::Range verticeSets=_pGTT->geoRanges()[0];
00339   int num_vertices = verticeSets.size();
00340   moab::Range bottomEdges;
00341   int j=0;
00342   for (j=0; j<num_vertices/2; j++)// we know that half are on bottom, we verified that already
00343   {
00344     moab::Range edges;
00345     rval = _fbe-> getEntAdj(verticeSets[j], /*edge type*/1, edges);
00346     MBERRCHK(rval, MBI);
00347     bottomEdges.merge(edges);
00348   }
00349   // we need to decide bottom faces before we create new faces, edges, etc
00350   moab::Range bottomFaces;
00351   for (j=0; j<num_vertices/2; j++)// we know that half are on bottom, we verified that already
00352   {
00353     moab::Range faces;
00354     rval = _fbe-> getEntAdj(verticeSets[j], /*face type*/2, faces);
00355     MBERRCHK(rval, MBI);
00356     bottomFaces.merge(faces);
00357   }
00358 
00359   // we have verified that the first half are from bottom surfaces
00360   // for each bottom face we will create a volume;
00361   // start with the edges, create weaving faces between edges (no new vertices!)
00362   // now, for each bottom edge, create a weaving face
00363   std::vector<moab::EntityHandle> newFaces;
00364   newFaces.resize(bottomEdges.size());
00365   for (j=0; j<(int)bottomEdges.size(); j++)
00366   {
00367     rval = _fbe->weave_lateral_face_from_edges(bottomEdges[j], edgeMap[bottomEdges[j]],  _direction, newFaces[j]);
00368     MBERRCHK(rval, MBI);
00369   }
00370   // now, to create volumes, we need to loop over bottom faces and add volumes one by one, and the
00371   // orientation of faces in them
00372 
00373   int volumeMatId = 0;// just give a mat id, for debugging, mostly
00374   moab::Tag matTag;
00375   rval = MBI->tag_get_handle("MATERIAL_SET", 1, moab::MB_TYPE_INTEGER, matTag);
00376   MBERRCHK(rval, MBI);
00377 
00378   for (j = 0; j<(int)bottomFaces.size(); j++)
00379   {
00380     // create volume here , finally!!!
00381     moab::EntityHandle volume;
00382     rval = MBI->create_meshset(moab::MESHSET_SET, volume);
00383     MBERRCHK(rval, MBI);
00384     volumeMatId++;
00385     rval= MBI->tag_set_data(matTag, &volume, 1, &volumeMatId);
00386     MBERRCHK(rval, MBI);
00387     moab::EntityHandle botFace = bottomFaces[j];
00388     moab::EntityHandle topFace = faceMap[botFace];
00389     // start copy
00390     // get the edges of bot face
00391     rval= MBI->add_parent_child(volume, botFace);
00392     MBERRCHK(rval, MBI);
00393 
00394     rval= MBI->add_parent_child(volume, topFace);
00395     MBERRCHK(rval, MBI);
00396 
00397     rval = _pGTT->add_geo_set(volume, 3);
00398     MBERRCHK(rval, MBI);
00399 
00400     // set senses
00401     // bottom face is negatively oriented, its normal is toward interior of the volume
00402     rval = _pGTT->set_sense(botFace, volume, -1);
00403     MBERRCHK(rval, MBI);
00404 
00405     // the top face is positively oriented
00406     rval = _pGTT->set_sense(topFace, volume, 1);
00407     MBERRCHK(rval, MBI);
00408 
00409     // the children should be in the same direction
00410     //   get the side edges of each face, and form lateral faces, along direction
00411     std::vector<moab::EntityHandle> edges1;
00412 
00413     rval = MBI->get_child_meshsets(botFace, edges1); // no hops
00414     MBERRCHK(rval, MBI);
00415 
00416     for (unsigned int i = 0; i < edges1.size(); ++i)
00417     {
00418       // the orientation of edge in face  will give the sense of face in volume
00419       int indexB = bottomEdges.index(edges1[i]);
00420       if (indexB<=-1)
00421         MBERRCHK(moab::MB_FAILURE, MBI);
00422       int sense = 1;
00423       rval = _pGTT->get_sense(edges1[i], botFace, sense);
00424       MBERRCHK(rval, MBI);
00425       rval=MBI->add_parent_child(volume, newFaces[indexB]);
00426       MBERRCHK(rval, MBI);
00427 
00428       // set sense as the one decided from the bottom edge sense within the bottom face
00429       rval = _pGTT->set_sense(newFaces[indexB], volume, sense);
00430       MBERRCHK(rval, MBI);
00431     }
00432     // end copy
00433   }
00434 
00435 
00436   delete _fbe;
00437   delete _pGTT; // when we are done, remove the _pGTT;
00438   // at this point, the result will be the model root set of the first ment of the model
00439   // ( (*(mentSelection.begin()) ).second )
00440 }
00441 
00442 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines