MeshKit  1.0
iMeshRefine2D.cpp
Go to the documentation of this file.
00001 #include "meshkit/MeshRefine2D.hpp"
00002 
00004 
00005 int MeshRefine2D :: initialize()
00006 {
00007   int err, result, namelen;
00008   numfacesRefined = 0;
00009 
00010   insertedNodes.clear();
00011   insertedFaces.clear();
00012 
00013   const char *tag1 = "VERTEX_ON_EDGE";
00014   namelen = strlen( tag1 );
00015   iMesh_createTag(mesh, tag1, 1, iBase_ENTITY_HANDLE, &vertex_on_edge_tag, &result, namelen); 
00016 
00017   const char *tag2 = "RemoveFace";
00018   namelen = strlen( tag2 );
00019   iMesh_createTag(mesh, tag2, 1, iBase_INTEGER,  &remove_tag, &result, namelen); 
00020 
00021   const char *tag3 = "Boundary";
00022   namelen = strlen( tag3 );
00023   iMesh_createTag(mesh, tag3, 1, iBase_INTEGER,  &boundary_tag, &result, namelen); 
00024 
00025   iMesh_getTagHandle(mesh, "GLOBAL_ID", &globalID_tag, &err, strlen("GLOBAL_ID") );
00026 
00027   iMesh_getRootSet(mesh, &rootSet, &err);
00028 
00029   return 0;
00030 }
00031 
00033 
00034 int MeshRefine2D :: finalize() 
00035 {
00036   int err;
00037 
00038   SimpleArray<iBase_EntityHandle> faces;
00039   iMesh_getEntities( mesh, rootSet, iBase_FACE, iMesh_ALL_TOPOLOGIES, 
00040                        ARRAY_INOUT(faces), &err); assert( !err );
00041 
00042   int val;
00043   for( size_t i = 0; i < faces.size(); i++) {
00044        iMesh_getIntData( mesh, faces[i], remove_tag, &val, &err);
00045        if( err == 0 && val == 1 ) {
00046            iMesh_rmvTag(mesh, faces[i], remove_tag, &err);
00047            iMesh_deleteEnt( mesh, faces[i], &err);
00048        }
00049   }
00050   iMesh_destroyTag(mesh, remove_tag, 1, &err);
00051 
00052   SimpleArray<iBase_EntityHandle> nodes;
00053   iMesh_getEntities( mesh, rootSet, iBase_VERTEX, iMesh_ALL_TOPOLOGIES, 
00054                        ARRAY_INOUT(nodes), &err);
00055 
00056   for( size_t i = 0; i < nodes.size(); i++) 
00057        iMesh_setIntData(mesh, nodes[i], globalID_tag, i, &err);
00058 
00059   return 0;
00060 }
00061 
00063 
00064 bool
00065 MeshRefine2D::searchEdge( iBase_EntityHandle v1, iBase_EntityHandle v2,
00066                           iBase_EntityHandle &edgehandle ) const
00067 {
00068    int err;
00069    iBase_EntityHandle vmin = min(v1,v2);
00070    iBase_EntityHandle vmax = max(v1,v2);
00071 
00072    SimpleArray<iBase_EntityHandle> edges;
00073    iMesh_getEntAdj(mesh, vmin, iBase_EDGE, ARRAY_INOUT(edges), &err);
00074 
00075    SimpleArray<iBase_EntityHandle> edgenodes;
00076    for( size_t i = 0; i < edges.size(); i++) {
00077         iMesh_getEntAdj(mesh, edges[i], iBase_VERTEX, ARRAY_INOUT(edgenodes), &err);
00078         if( edgenodes[0] == vmin && edgenodes[1] == vmax )  {
00079             edgehandle = edges[i];
00080             return 1;
00081         }
00082    }
00083    return 0;
00084 }
00085 
00087 
00088 bool 
00089 MeshRefine2D:: allow_edge_refinement( iBase_EntityHandle &edgehandle ) const
00090 {
00091    int err, bound_val;
00092    iMesh_getIntData( mesh, edgehandle, boundary_tag, &bound_val, &err);
00093 
00094    // If err == 0, means edge is boundary 
00095    if( !err && !boundary_split_flag ) return 0;
00096 
00097    return 1;
00098 }
00099 
00101 
00102 iBase_EntityHandle 
00103 MeshRefine2D::create_new_edge( iBase_EntityHandle v1, iBase_EntityHandle v2)
00104 {
00105   int status, err;
00106   iBase_EntityHandle edgehandle;
00107   static vector<iBase_EntityHandle> pnodes(2);
00108 
00109   pnodes[0] = std::min(v1,v2);
00110   pnodes[1] = std::max(v1,v2);
00111   iMesh_createEnt(mesh, iMesh_LINE_SEGMENT, &pnodes[0], 2, &edgehandle, &status, &err);
00112 
00113   return edgehandle;
00114 }
00115 
00117 
00118 Point3D 
00119 MeshRefine2D:: edge_centroid( iBase_EntityHandle v1, iBase_EntityHandle v2) const
00120 {
00121   int  err;
00122   double x, y, z;
00123   Point3D p3d;
00124 
00125   p3d[0] = 0.0;
00126   p3d[1] = 0.0;
00127   p3d[2] = 0.0;
00128 
00129   iMesh_getVtxCoord(mesh, v1, &x, &y, &z, &err);
00130   p3d[0] += x; p3d[1] += y; p3d[2] += z;
00131 
00132   iMesh_getVtxCoord(mesh, v2, &x, &y, &z, &err);
00133   p3d[0] += x; p3d[1] += y; p3d[2] += z;
00134 
00135   p3d[0] *= 0.5;
00136   p3d[1] *= 0.5;
00137   p3d[2] *= 0.5;
00138 
00139   return p3d;
00140 }
00141 
00143 
00144 double  
00145 MeshRefine2D:: edge_length( iBase_EntityHandle v1, iBase_EntityHandle v2) const
00146 {
00147   int  err;
00148 
00149   double x0, y0, z0;
00150   iMesh_getVtxCoord(mesh, v1, &x0, &y0, &z0, &err);
00151 
00152   double x1, y1, z1;
00153   iMesh_getVtxCoord(mesh, v2, &x1, &y1, &z1, &err);
00154 
00155   double dx = x1-x0;
00156   double dy = y1-y0;
00157   double dz = z1-z0;
00158 
00159   double len = sqrt(dx*dx + dy*dy + dz*dz );
00160 
00161   return len;
00162 }
00163 
00165 
00166 Point3D 
00167 MeshRefine2D:: face_centroid( iBase_EntityHandle facehandle) const
00168 {
00169   int  err;
00170   double  x, y, z;
00171 
00172   SimpleArray<iBase_EntityHandle> facenodes;
00173   iMesh_getEntAdj(mesh, facehandle, iBase_VERTEX, ARRAY_INOUT(facenodes), &err);
00174 
00175   size_t numNodes = facenodes.size(); assert( numNodes > 2 );
00176 
00177   Point3D p3d;
00178   p3d[0] = 0.0;
00179   p3d[1] = 0.0;
00180   p3d[2] = 0.0;
00181   for(size_t i =  0; i < numNodes; i++) {
00182       iMesh_getVtxCoord(mesh, facenodes[i], &x, &y, &z, &err);
00183       p3d[0] += x;
00184       p3d[1] += y;
00185       p3d[2] += z;
00186   }
00187 
00188   p3d[0] /= (double)numNodes;
00189   p3d[1] /= (double)numNodes;
00190   p3d[2] /= (double)numNodes;
00191 
00192   return p3d;
00193 }
00195 
00196 double  
00197 MeshRefine2D:: face_aspect_ratio( iBase_EntityHandle facehandle) const
00198 {
00199   int  err;
00200   double  x, y, z;
00201 
00202   SimpleArray<iBase_EntityHandle> facenodes;
00203   iMesh_getEntAdj(mesh, facehandle, iBase_VERTEX, ARRAY_INOUT(facenodes), &err);
00204 
00205   size_t numNodes = facenodes.size(); assert( numNodes > 2 );
00206 
00207   vector<double> elen;
00208   elen.resize( numNodes );
00209 
00210   Point3D p3d;
00211   for(size_t i =  0; i < numNodes; i++) {
00212       iBase_EntityHandle  v0 = facenodes[i];
00213       iBase_EntityHandle  v1 = facenodes[(i+1)%numNodes];
00214       elen[i] = edge_length(v0,v1);
00215   }
00216 
00217   double minlen = *min_element( elen.begin(), elen.end() );
00218   double maxlen = *max_element( elen.begin(), elen.end() );
00219 
00220   return minlen/maxlen;
00221 }
00223 
00224 int  
00225 MeshRefine2D::setVertexOnEdge( iBase_EntityHandle v1, iBase_EntityHandle v2,
00226                                iBase_EntityHandle &vmid ) 
00227 {
00228   int err;
00229   iBase_EntityHandle edgehandle;
00230 
00231   // Case I: Edge doesn't exist then create and set the tag value
00232   if( searchEdge(v1,v2,edgehandle) == 0) 
00233   {
00234       edgehandle = create_new_edge(v1,v2);
00235       if( allow_edge_refinement(edgehandle) ) 
00236       {
00237           Point3D p3d = edge_centroid(v1,v2);
00238           iMesh_createVtx(mesh, p3d[0], p3d[1], p3d[2], &vmid, &err);
00239           iMesh_setEHData(mesh, edgehandle, vertex_on_edge_tag, vmid, &err);
00240           return 0;
00241       }
00242       return 1;
00243   } 
00244 
00245   //Edge exist but not the tag value.
00246   if( allow_edge_refinement( edgehandle) )  
00247   {
00248       iMesh_getEHData(mesh, edgehandle, vertex_on_edge_tag, &vmid, &err);
00249       if( err ) {
00250           Point3D p3d = edge_centroid(v1,v2);
00251           iMesh_createVtx(mesh, p3d[0], p3d[1], p3d[2], &vmid, &err);
00252           iMesh_setEHData(mesh, edgehandle, vertex_on_edge_tag, vmid, &err);
00253       }
00254       return 0;
00255   }
00256 
00257   return 1;
00258 }
00259 
00261 
00262 int MeshRefine2D::getVertexOnEdge( iBase_EntityHandle v1, iBase_EntityHandle v2,
00263                                    iBase_EntityHandle &hangVertex ) const
00264 {
00265    int err;
00266    iBase_EntityHandle edgehandle;
00267 
00268    searchEdge( v1, v2, edgehandle) ;
00269    iMesh_getEHData(mesh, edgehandle, vertex_on_edge_tag, &hangVertex, &err);
00270 
00271    if( err ) hangVertex = 0;
00272 
00273    return err;
00274 }
00275 
00277 
00278 int Sqrt3Refine2D :: execute()
00279 {
00280 /*
00281    CentroidRefine2D refine(mesh);
00282    EdgeFlip eflip(mesh);
00283 
00284    for ( int itime = 0; itime < numIterations; itime++) {
00285        refine.execute();
00286        eflip.execute();
00287    }
00288 */
00289 
00290    return 0;
00291 }
00292 
00294 
00295 int LongestEdgeRefine2D :: initialize() 
00296 {
00297   MeshRefine2D::initialize();
00298   return 0;
00299 }
00300 
00302 
00303 int LongestEdgeRefine2D :: atomicOp( const iBase_EntityHandle &oldface)
00304 {
00305   int err;
00306   SimpleArray<iBase_EntityHandle> eConnect;
00307   iMesh_getEntAdj(mesh, oldface, iBase_VERTEX, ARRAY_INOUT(eConnect), &err);
00308 
00309   iBase_EntityHandle  v1, v2, vmid;
00310   vector<double> elen(3);
00311 
00312   double maxlen = 0.0;
00313   for( int i = 0; i < 3; i++) {
00314       v1      = eConnect[(i+1)%3];
00315       v2      = eConnect[(i+2)%3];
00316       elen[i] = edge_length( v1, v2);
00317       maxlen  = std::max(elen[i], maxlen);
00318   }
00319 
00320   for( int i = 0; i < 3; i++) {
00321       if( elen[i]/maxlen > 0.90 ) 
00322        {
00323            v1 = eConnect[(i+1)%3];
00324            v2 = eConnect[(i+2)%3];
00325            setVertexOnEdge(v1,v2, vmid);
00326         }
00327    }
00328    return 0;
00329 }
00330 
00332 
00333 int LongestEdgeRefine2D::execute() 
00334 {
00335   int err;
00336   initialize();
00337 
00338   SimpleArray<iBase_EntityHandle> faces;
00339   iMesh_getEntities( mesh, rootSet, iBase_FACE, iMesh_ALL_TOPOLOGIES, 
00340                        ARRAY_INOUT(faces), &err);
00341 
00342   size_t ncount = 0;
00343   for( size_t i = 0; i < faces.size(); i++) {
00344         double ratio  = face_aspect_ratio( faces[i] );
00345         if( ratio < cutOffAspectRatio) {
00346              err = atomicOp( faces[i] );
00347              if( !err ) ncount++;
00348         }
00349   }
00350 
00351   if( ncount ) {
00352   ConsistencyRefine2D consistency;
00353   consistency.setMesh(mesh);
00354   consistency.execute();
00355   finalize();
00356   } else
00357     cout << "Warning: No Edge was refined " << endl;
00358 
00359   return 0;
00360 }
00361 
00363 
00364 
00365 int ConsistencyRefine2D :: initialize() 
00366 {
00367   hangingVertex.resize(3);
00368   edge0.set(0);
00369   edge1.set(1);
00370   edge2.set(2);
00371 
00372   insertedNodes.clear();
00373   insertedFaces.clear();
00374 
00375   MeshRefine2D::initialize();
00376 
00377   return 0;
00378 }
00379 
00380 //#############################################################################
00381 
00382 int  ConsistencyRefine2D :: execute()
00383 {
00384     initialize();
00385 
00386     makeConsistent();
00387 
00388     finalize();
00389 
00390     return 0;
00391 }
00392 
00393 //#############################################################################
00394 
00395 void ConsistencyRefine2D :: subDivideQuad2Tri( const vector<iBase_EntityHandle> &connect)
00396 {
00397   int err, status;
00398   assert( connect.size() == 4 );
00399   //********************************************************************
00400   // Subdivide a quadrilateral cell into two triangles. We can choose
00401   // either quadrilateral diagonal for spliiting, but we choose to
00402   // select quadrilateral which gives, maximum of minimum aspect
00403   // ratio of the resulting two triangles ....
00404   //*******************************************************************
00405   double diagonal[2];
00406 
00407   iBase_EntityHandle v0 = connect[0]; 
00408   iBase_EntityHandle v1 = connect[1]; 
00409   iBase_EntityHandle v2 = connect[2]; 
00410   iBase_EntityHandle v3 = connect[3]; 
00411 
00412   diagonal[0]  = edge_length( v0, v3 );
00413   diagonal[1]  = edge_length( v1, v2 );
00414 
00415   iBase_EntityHandle facehandle;
00416   vector<iBase_EntityHandle> pnodes(3);
00417 
00418   if( diagonal[0] < diagonal[1] ) {
00419     pnodes[0] = v0;
00420     pnodes[1] = v1;
00421     pnodes[2] = v3;
00422     iMesh_createEnt(mesh, iMesh_TRIANGLE, &pnodes[0], 3, &facehandle, &status, &err);
00423     insertedFaces.push_back( facehandle );
00424 
00425     pnodes[0] = v0;
00426     pnodes[1] = v3;
00427     pnodes[2] = v2;
00428     iMesh_createEnt(mesh, iMesh_TRIANGLE, &pnodes[0], 3, &facehandle, &status, &err);
00429     insertedFaces.push_back( facehandle );
00430   } else {
00431     pnodes[0] = v0;
00432     pnodes[1] = v1;
00433     pnodes[2] = v2;
00434     iMesh_createEnt(mesh, iMesh_TRIANGLE, &pnodes[0], 3, &facehandle, &status, &err);
00435     insertedFaces.push_back( facehandle );
00436 
00437     pnodes[0] = v1;
00438     pnodes[1] = v3;
00439     pnodes[2] = v2;
00440 
00441     iMesh_createEnt(mesh, iMesh_TRIANGLE, &pnodes[0], 3, &facehandle, &status, &err);
00442     insertedFaces.push_back( facehandle );
00443   }
00444 }
00445 
00446 //#############################################################################
00447 
00448 void ConsistencyRefine2D :: makeConsistent1( const iBase_EntityHandle &oldface)
00449 {
00450   int err, status;
00451   //------------------------------------------------------------------
00452   // When only one edge is inconsistent, we can direcly join the
00453   // hanging node to the opposite node of the triangle. Therefore
00454   // one additional triangle is generated with this case.
00455   //------------------------------------------------------------------
00456   SimpleArray<iBase_EntityHandle> pnodes;
00457   iMesh_getEntAdj(mesh, oldface, iBase_VERTEX, ARRAY_INOUT(pnodes), &err);
00458 
00459   if( pnodes.size() != 3 ) return;
00460 
00461   iBase_EntityHandle n1 = pnodes[0];
00462   iBase_EntityHandle n2 = pnodes[1];
00463   iBase_EntityHandle n3 = pnodes[2];
00464 
00465   int val = 1;
00466   iMesh_setIntData(mesh, oldface, remove_tag, val, &err);
00467 
00468   iBase_EntityHandle facehandle;
00469 
00470   // If the only hanging vertex lies of the edge0. i.e. opposite to the
00471   // vertex 0 of the existing triangle.
00472   if( bitvec == edge0) {
00473     pnodes[0] = n1;
00474     pnodes[1] = n2;
00475     pnodes[2] = hangingVertex[0];
00476 
00477     iMesh_createEnt(mesh, iMesh_TRIANGLE, &pnodes[0], 3, &facehandle, &status, &err);
00478     insertedFaces.push_back(facehandle);
00479 
00480     // Create a new triangle ....
00481     pnodes[0] = n3;
00482     pnodes[1] = n1;
00483     pnodes[2] = hangingVertex[0];
00484     iMesh_createEnt(mesh, iMesh_TRIANGLE, &pnodes[0], 3, &facehandle, &status, &err);
00485     insertedFaces.push_back(facehandle);
00486     return;
00487   }
00488 
00489   // If the only hanging vertex lies of the edge1. i.e. opposite to the
00490   // vertex 1 of the existing triangle.
00491   if( bitvec == edge1) 
00492   {
00493     // Replace the existing triangle ....
00494     pnodes[0] = n2;
00495     pnodes[1] = n3;
00496     pnodes[2] = hangingVertex[1];
00497     iMesh_createEnt(mesh, iMesh_TRIANGLE, &pnodes[0], 3, &facehandle, &status, &err);
00498     insertedFaces.push_back(facehandle);
00499 
00500     // Create a new triangle ....
00501     pnodes[0] = n2;
00502     pnodes[1] = hangingVertex[1];
00503     pnodes[2] = n1;
00504 
00505     iMesh_createEnt(mesh, iMesh_TRIANGLE, &pnodes[0], 3, &facehandle, &status, &err);
00506     insertedFaces.push_back(facehandle);
00507 
00508     return;
00509   }
00510 
00511   // If the only hanging vertex lies of the edge2. i.e. opposite to the
00512   // vertex 2 of the existing triangle.
00513   if( bitvec == edge2) {
00514     // Replace the existing triangle ....
00515     pnodes[0] = n3;
00516     pnodes[1] = n1;
00517     pnodes[2] = hangingVertex[2];
00518 
00519     iMesh_createEnt(mesh, iMesh_TRIANGLE, &pnodes[0], 3, &facehandle, &status, &err);
00520     insertedFaces.push_back(facehandle);
00521 
00522     // Create a new triangle ....
00523     pnodes[0] = n3;
00524     pnodes[1] = hangingVertex[2];
00525     pnodes[2] = n2;
00526 
00527     iMesh_createEnt(mesh, iMesh_TRIANGLE, &pnodes[0], 3, &facehandle, &status, &err);
00528     insertedFaces.push_back(facehandle);
00529 
00530     return;
00531   }
00532 
00533 }
00534 
00535 //#############################################################################
00536 
00537 void ConsistencyRefine2D :: refineEdge0(const iBase_EntityHandle &oldface)
00538 {
00539   int err, status;
00540   iBase_EntityHandle facehandle;
00541 
00542   SimpleArray<iBase_EntityHandle> eConnect;
00543   iMesh_getEntAdj(mesh, oldface, iBase_VERTEX, ARRAY_INOUT(eConnect), &err);
00544   assert( eConnect.size() == 3 );
00545 
00546   iBase_EntityHandle n1 = eConnect[0];
00547   iBase_EntityHandle n2 = eConnect[1];
00548   iBase_EntityHandle n3 = eConnect[2];
00549 
00550   static vector<iBase_EntityHandle> tnodes(3);
00551   tnodes[0] = n1;
00552   tnodes[1] = hangingVertex[2];
00553   tnodes[2] = hangingVertex[1];
00554   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tnodes[0], 3, &facehandle, &status, &err);
00555   insertedFaces.push_back(facehandle);
00556 
00557   // One Quadrilateral is created, divide into 2 triangles.
00558   static vector<iBase_EntityHandle> qnodes(4);
00559   qnodes[0] =  hangingVertex[2];
00560   qnodes[1] =  n2;
00561   qnodes[2] =  hangingVertex[1];
00562   qnodes[3] =  n3;
00563 
00564   subDivideQuad2Tri(qnodes );
00565 }
00566 
00567 //#############################################################################
00568 
00569 void ConsistencyRefine2D :: refineEdge1(const iBase_EntityHandle &oldface)
00570 {
00571   int err, status;
00572   iBase_EntityHandle facehandle;
00573 
00574   SimpleArray<iBase_EntityHandle> eConnect;
00575   iMesh_getEntAdj(mesh, oldface, iBase_VERTEX, ARRAY_INOUT(eConnect), &err);
00576   assert( eConnect.size() == 3 );
00577 
00578   iBase_EntityHandle n1 = eConnect[0];
00579   iBase_EntityHandle n2 = eConnect[1];
00580   iBase_EntityHandle n3 = eConnect[2];
00581 
00582   static vector<iBase_EntityHandle> tnodes(3);
00583   tnodes[0] = n2;
00584   tnodes[1] = hangingVertex[0];
00585   tnodes[2] = hangingVertex[2];
00586 
00587   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tnodes[0], 3, &facehandle, &status, &err);
00588   insertedFaces.push_back( facehandle );
00589 
00590   // One Quadrilateral is created, divide into 2 triangles.
00591   static vector<iBase_EntityHandle> qnodes(4);
00592   qnodes[0] =  n1;
00593   qnodes[1] =  hangingVertex[2];
00594   qnodes[2] =  n3;
00595   qnodes[3] =  hangingVertex[0];
00596 
00597   subDivideQuad2Tri( qnodes );
00598 }
00599 
00600 
00601 //#############################################################################
00602 
00603 void ConsistencyRefine2D :: refineEdge2(const iBase_EntityHandle &oldface)
00604 {
00605   int err, status;
00606   iBase_EntityHandle facehandle;
00607 
00608   SimpleArray<iBase_EntityHandle> eConnect;
00609   iMesh_getEntAdj(mesh, oldface, iBase_VERTEX, ARRAY_INOUT(eConnect), &err);
00610   assert( eConnect.size() == 3 );
00611 
00612   iBase_EntityHandle n1 = eConnect[0];
00613   iBase_EntityHandle n2 = eConnect[1];
00614   iBase_EntityHandle n3 = eConnect[2];
00615 
00616   static vector<iBase_EntityHandle> tnodes(3);
00617   tnodes[0] = n3;
00618   tnodes[1] = hangingVertex[1];
00619   tnodes[2] = hangingVertex[0];
00620   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tnodes[0], 3, &facehandle, &status, &err);
00621   insertedFaces.push_back( facehandle );
00622 
00623   static vector<iBase_EntityHandle> qnodes(4);
00624   qnodes[0] =  n1;
00625   qnodes[1] =  n2;
00626   qnodes[2] =  hangingVertex[1];
00627   qnodes[3] =  hangingVertex[0];
00628 
00629   subDivideQuad2Tri( qnodes );
00630 }
00631 
00632 
00633 //#############################################################################
00634 
00635 void ConsistencyRefine2D :: makeConsistent2( const iBase_EntityHandle &oldface)
00636 {
00637   //--------------------------------------------------------------------
00638   // When there are two edges which are inconsistent, then we create
00639   // one triangle and one quadrilateral. This quadrilateral is further
00640   // divided into 2 triangle, which produces better aspect ratio.
00641   // Therefore, three triangles are generated in this procedure.
00642   //--------------------------------------------------------------------
00643   // Find out which edge is consistent ...
00644   bitvec.flip();
00645 
00646   int err, val = 1;
00647   iMesh_setIntData(mesh, oldface, remove_tag, val, &err);
00648 
00649   if( bitvec == edge0) {
00650     refineEdge0(oldface); 
00651     return;
00652   }
00653 
00654   if( bitvec == edge1) {
00655     refineEdge1(oldface); 
00656     return;
00657   }
00658 
00659   if( bitvec == edge2) {
00660     refineEdge2(oldface); 
00661     return;
00662   }
00663 
00664 }
00665 
00666 //#############################################################################
00667 
00668 void ConsistencyRefine2D :: makeConsistent3( const iBase_EntityHandle &oldface)
00669 {
00670   int err, status;
00671 
00672   int val = 1;
00673   iMesh_setIntData(mesh, oldface, remove_tag, val, &err);
00674 
00675   SimpleArray<iBase_EntityHandle> eConnect;
00676   iMesh_getEntAdj(mesh, oldface, iBase_VERTEX, ARRAY_INOUT(eConnect), &err);
00677   assert( eConnect.size() == 3 );
00678 
00679   iBase_EntityHandle n1 = eConnect[0];
00680   iBase_EntityHandle n2 = eConnect[1];
00681   iBase_EntityHandle n3 = eConnect[2];
00682 
00683   iBase_EntityHandle facehandle;
00684 
00685   static vector<iBase_EntityHandle> tnodes(3);
00686   // First Triangle  Using 0(old), 1,2(new): Replace existing triangle
00687   tnodes[0] = n1;
00688   tnodes[1] = hangingVertex[2];
00689   tnodes[2] = hangingVertex[1];
00690   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tnodes[0], 3, &facehandle, &status, &err);
00691   insertedFaces.push_back(facehandle);
00692 
00693   // Second Triangle  Using 1(old), 2,0(new): Create new triangle
00694   tnodes[0] = n2;
00695   tnodes[1] = hangingVertex[0];
00696   tnodes[2] = hangingVertex[2];
00697   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tnodes[0], 3, &facehandle, &status, &err);
00698   insertedFaces.push_back(facehandle);
00699 
00700   // Second Triangle  Using 2(old), 1,0(new): Create new triangle
00701   tnodes[0] = n3;
00702   tnodes[1] = hangingVertex[1];
00703   tnodes[2] = hangingVertex[0];
00704   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tnodes[0], 3, &facehandle, &status, &err);
00705   insertedFaces.push_back(facehandle);
00706 
00707   //  All new only : Create new triangle
00708   tnodes[0] = hangingVertex[0];
00709   tnodes[1] = hangingVertex[1];
00710   tnodes[2] = hangingVertex[2];
00711   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tnodes[0], 3, &facehandle, &status, &err);
00712   insertedFaces.push_back(facehandle);
00713 }
00714 
00715 //#############################################################################
00716 
00717 void ConsistencyRefine2D :: checkFaceConsistency( const iBase_EntityHandle &oldface ) 
00718 {
00719   int err;
00720   SimpleArray<iBase_EntityHandle> eConnect;
00721   iMesh_getEntAdj(mesh, oldface, iBase_VERTEX, ARRAY_INOUT(eConnect), &err);
00722 
00723   iBase_EntityHandle v1,v2;
00724 
00725   bitvec.reset();
00726   for( int i = 0; i < 3; i++) {
00727     v1  = eConnect[ (i+1)%3 ];
00728     v2  = eConnect[ (i+2)%3 ];
00729     getVertexOnEdge( v1, v2, hangingVertex[i] );
00730     if( hangingVertex[i] ) {
00731         bitvec.set(i);
00732     }
00733   }
00734 
00735 }
00736 
00737 //#############################################################################
00738 
00739 int ConsistencyRefine2D :: atomicOp(const iBase_EntityHandle &oldface)
00740 {
00741     int err, val;
00742     iMesh_getIntData(mesh, oldface, remove_tag, &val, &err);
00743 
00744     if( err == 0 && val == 1 ) return 1;
00745     
00746     checkFaceConsistency( oldface );
00747 
00748     switch( bitvec.count() )
00749     {
00750         case 1:
00751           numfacesRefined++;
00752           makeConsistent1( oldface );
00753           break;
00754         case 2:
00755           numfacesRefined++;
00756           makeConsistent2( oldface );
00757           break;
00758         case 3:
00759           numfacesRefined++;
00760           makeConsistent3( oldface );
00761           break;
00762      }
00763      return 0;
00764  }
00765 
00766 //#############################################################################
00767 
00768 void ConsistencyRefine2D :: makeConsistent()
00769 {
00770   //**********************************************************************
00771   // The previous step, will leave some hanging nodes. So adjacent cells 
00772   // will be forced to refined. All the hanging nodes are attached 
00773   // directly to the opposite node of the triangle, thus creating two 
00774   // triangle. This may produce some bad triangle, which could be 
00775   // improved using edge swapping algorithm. In this process, only new 
00776   // edges are created and  no new vertices are introduced.
00777   //**********************************************************************
00778 
00779   int err;
00780   SimpleArray<iBase_EntityHandle> faces;
00781   iMesh_getEntities( mesh, rootSet, iBase_FACE, iMesh_ALL_TOPOLOGIES, 
00782                        ARRAY_INOUT(faces), &err);
00783 
00784   for( size_t i = 0; i < faces.size(); i++) 
00785         atomicOp( faces[i] );
00786 }
00787 
00788 
00789 //#############################################################################
00790 
00791 int CentroidRefine2D::refine_tri(const iBase_EntityHandle &oldface)
00792 {
00793   int err, status;
00794   iBase_EntityHandle facehandle, vcenter;
00795 
00796   Point3D p3d = face_centroid(oldface);
00797   iMesh_createVtx(mesh, p3d[0], p3d[1], p3d[2], &vcenter, &err);
00798 
00799   SimpleArray<iBase_EntityHandle> eConnect;
00800   iMesh_getEntAdj(mesh, oldface, iBase_VERTEX, ARRAY_INOUT(eConnect), &err);
00801 
00802   assert( eConnect.size() == 3 );
00803 
00804   iBase_EntityHandle v1 = eConnect[0];
00805   iBase_EntityHandle v2 = eConnect[1];
00806   iBase_EntityHandle v3 = eConnect[2];
00807 
00808   static vector<iBase_EntityHandle> tconn(3);
00809 
00810   tconn[0] = vcenter;
00811   tconn[1] = v1;
00812   tconn[2] = v2;
00813   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tconn[0], 3, &facehandle, &status, &err);
00814   insertedFaces.push_back( facehandle );
00815 
00816   tconn[0] = vcenter;
00817   tconn[1] = v2;
00818   tconn[2] = v3;
00819   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tconn[0], 3, &facehandle, &status, &err);
00820   insertedFaces.push_back( facehandle );
00821 
00822   tconn[0] = vcenter;
00823   tconn[1] = v3;
00824   tconn[2] = v1;
00825   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tconn[0], 3, &facehandle, &status, &err);
00826   insertedFaces.push_back( facehandle );
00827 
00828   int val = 1;
00829   iMesh_setIntData(mesh, oldface, remove_tag, val, &err);
00830 
00831   return 0;
00832 }
00833 
00834 
00836 
00837 int  CentroidRefine2D::refine_quad(const iBase_EntityHandle &oldface)
00838 {
00839   int status, err;
00840   iBase_EntityHandle  facehandle, vcenter;
00841 
00842   Point3D p3d = face_centroid(oldface);
00843   iMesh_createVtx(mesh, p3d[0], p3d[1], p3d[2], &vcenter, &err);
00844 
00845   SimpleArray<iBase_EntityHandle> eConnect;
00846   iMesh_getEntAdj(mesh, oldface, iBase_VERTEX, ARRAY_INOUT(eConnect), &err);
00847 
00848   iBase_EntityHandle v0 = eConnect[0];
00849   iBase_EntityHandle v1 = eConnect[1];
00850   iBase_EntityHandle v2 = eConnect[2];
00851   iBase_EntityHandle v3 = eConnect[3];
00852 
00853   vector<iBase_EntityHandle> tconn(3);
00854 
00855   tconn[0] = vcenter;
00856   tconn[1] = v0;
00857   tconn[2] = v1;
00858   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tconn[0], 3, &facehandle, &status, &err);
00859   insertedFaces.push_back( facehandle );
00860 
00861   tconn[0] = vcenter;
00862   tconn[1] = v1;
00863   tconn[2] = v3;
00864   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tconn[0], 3, &facehandle, &status, &err);
00865   insertedFaces.push_back( facehandle );
00866 
00867   tconn[0] = vcenter;
00868   tconn[1] = v3;
00869   tconn[2] = v2;
00870   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tconn[0], 3, &facehandle, &status, &err);
00871   insertedFaces.push_back( facehandle );
00872 
00873   tconn[0] = vcenter;
00874   tconn[1] = v2;
00875   tconn[2] = v0;
00876   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tconn[0], 3, &facehandle, &status, &err);
00877   insertedFaces.push_back( facehandle );
00878 
00879   int val = 1;
00880   iMesh_setIntData(mesh, oldface, remove_tag, val, &err);
00881 
00882   return 0;
00883 }
00884 
00885 
00887 
00888 int CentroidRefine2D::atomicOp(const iBase_EntityHandle &oldface)
00889 {
00890    int err, face_topo;
00891    iMesh_getEntTopo(mesh, oldface, &face_topo, &err);
00892    
00893    switch( face_topo ) 
00894    {
00895      case iMesh_TRIANGLE:
00896          refine_tri(oldface);
00897          break;
00898      case iMesh_QUADRILATERAL:
00899          refine_quad(oldface);
00900          break;
00901      default:
00902          cout << "Warning: Element not supported for refinement " << endl;
00903          break;
00904   }
00905 
00906   return 0;
00907 }
00908 
00910 
00911 int CentroidRefine2D::initialize()
00912 {
00913   int err, result;
00914 
00915   const char *tag = "RemoveFace";
00916   int namelen = strlen( tag );
00917   iMesh_createTag(mesh, tag, 1, iBase_INTEGER,  &remove_tag, &result, namelen); 
00918 
00919   iMesh_getTagHandle(mesh, "GLOBAL_ID", &globalID_tag, &err, strlen("GLOBAL_ID") );
00920 
00921   iMesh_getRootSet(mesh, &rootSet, &err);
00922 
00923   insertedFaces.clear();
00924 }
00925 
00927 
00928 
00929 int CentroidRefine2D::execute()
00930 {
00931   int err;
00932 
00933   initialize();
00934 
00935   SimpleArray<iBase_EntityHandle> faces;
00936   iMesh_getEntities( mesh, rootSet, iBase_FACE, iMesh_ALL_TOPOLOGIES, 
00937                        ARRAY_INOUT(faces), &err);
00938 
00939   for( int i = 0; i < faces.size(); i++)
00940       atomicOp( faces[i]  );
00941 
00942   finalize();
00943 
00944   return 0;
00945 }
00947 
00948 int ObtuseRefine2D :: initialize()
00949 {
00950    return 0;
00951 };
00952 
00954 
00955 int ObtuseRefine2D :: atomicOp( const iBase_EntityHandle &facehandle)
00956 {
00957     int err;
00958     SimpleArray<iBase_EntityHandle> facenodes;
00959     iMesh_getEntAdj(mesh, facehandle, iBase_VERTEX, ARRAY_INOUT(facenodes), &err);
00960 
00961     Point3D pv0, pv1, pv2;
00962     Point3D vec1, vec2;
00963 
00964     iBase_EntityHandle vmid;
00965     double angle;
00966     for( int i = 0; i < 3; i++) {
00967          getVertexCoords( facenodes[(i+0)%3], pv0);
00968          getVertexCoords( facenodes[(i+1)%3], pv1);
00969          getVertexCoords( facenodes[(i+2)%3], pv2);
00970          vec1 = Math::create_vector( pv2, pv0);
00971          vec2 = Math::create_vector( pv1, pv0);
00972          angle = Math::getVectorAngle(vec1, vec2);
00973          if( angle > cutoffAngle) {
00974             setVertexOnEdge(facenodes[(i+1)%2], facenodes[(i+2)%2], vmid );
00975             return 0;
00976          }
00977     }
00978 
00979     return 1;
00980 }
00981 
00983 int ObtuseRefine2D :: execute()
00984 {
00985   int err;
00986 
00987   initialize();
00988 
00989   SimpleArray<iBase_EntityHandle> faces;
00990   iMesh_getEntities( mesh, rootSet, iBase_FACE, iMesh_ALL_TOPOLOGIES, 
00991                        ARRAY_INOUT(faces), &err);
00992 
00993   size_t ncount = 0;
00994   for( size_t i = 0; i < faces.size(); i++) {
00995        err = atomicOp( faces[i]  );
00996        if( !err ) ncount++;
00997   }
00998 
00999   ConsistencyRefine2D refine;
01000   if( ncount ) {
01001       refine.setMesh(mesh);
01002       refine.execute();
01003       finalize();
01004   } else 
01005       cout << "Warning: No triangle was refined " << endl;
01006 
01007   return 0;
01008 }
01009 
01011 
01012 int Refine2D14 :: initialize() 
01013 {
01014   MeshRefine2D::initialize();
01015   return 0;
01016 }
01018 
01019 int Refine2D14::refine_quad(const iBase_EntityHandle &oldface)
01020 {
01021   int status, err;
01022   iBase_EntityHandle  facehandle, vcenter;
01023 
01024   Point3D p3d = face_centroid(oldface);
01025   iMesh_createVtx(mesh, p3d[0], p3d[1], p3d[2], &vcenter, &err);
01026 
01027   SimpleArray<iBase_EntityHandle> eConnect;
01028   iMesh_getEntAdj(mesh, oldface, iBase_VERTEX, ARRAY_INOUT(eConnect), &err);
01029   assert( eConnect.size() == 4 );
01030   
01031   iBase_EntityHandle v0 = eConnect[0];
01032   iBase_EntityHandle v1 = eConnect[1];
01033   iBase_EntityHandle v2 = eConnect[2];
01034   iBase_EntityHandle v3 = eConnect[3];
01035 
01036   iBase_EntityHandle v01, v23, v02, v13;
01037   err = setVertexOnEdge(v0,v1, v01);
01038   err = setVertexOnEdge(v2,v3, v23);
01039   err = setVertexOnEdge(v0,v2, v02);
01040   err = setVertexOnEdge(v1,v3, v13);
01041   
01042   static vector<iBase_EntityHandle> qconn(4);
01043 
01044   qconn[0] = v0;
01045   qconn[1] = v01;
01046   qconn[2] = v02;
01047   qconn[3] = vcenter;
01048   iMesh_createEnt(mesh, iMesh_QUADRILATERAL, &qconn[0], 4, &facehandle, &status, &err);
01049   insertedFaces.push_back( facehandle );
01050 
01051   qconn[0] = v01;
01052   qconn[1] = v1;
01053   qconn[2] = vcenter;
01054   qconn[3] = v13;
01055   iMesh_createEnt(mesh, iMesh_QUADRILATERAL, &qconn[0], 4, &facehandle, &status, &err);
01056   insertedFaces.push_back( facehandle );
01057 
01058   qconn[0] = v02;
01059   qconn[1] = vcenter;
01060   qconn[2] = v2;
01061   qconn[3] = v23;
01062   iMesh_createEnt(mesh, iMesh_QUADRILATERAL, &qconn[0], 4, &facehandle, &status, &err);
01063   insertedFaces.push_back( facehandle );
01064 
01065   qconn[0] = vcenter;
01066   qconn[1] = v13;
01067   qconn[2] = v23;
01068   qconn[3] = v3;
01069   iMesh_createEnt(mesh, iMesh_QUADRILATERAL, &qconn[0], 4, &facehandle, &status, &err);
01070   insertedFaces.push_back( facehandle );
01071 
01072   int val = 1;
01073   iMesh_setIntData(mesh, oldface, remove_tag, val, &err);
01074 
01075   return 0;
01076 }
01077 
01079 
01080 int Refine2D14:: refine_tri( const iBase_EntityHandle &oldface)
01081 {
01082   int err, status;
01083   iBase_EntityHandle  facehandle;
01084 
01085   SimpleArray<iBase_EntityHandle> eConnect;
01086   iMesh_getEntAdj(mesh, oldface, iBase_VERTEX, ARRAY_INOUT(eConnect), &err);
01087   assert( eConnect.size() == 3 );
01088 
01089   iBase_EntityHandle v0 = eConnect[0];
01090   iBase_EntityHandle v1 = eConnect[1];
01091   iBase_EntityHandle v2 = eConnect[2];
01092 
01093   iBase_EntityHandle v01, v12, v20;
01094  
01095   err = setVertexOnEdge(v0,v1, v01);    
01096   err = setVertexOnEdge(v1,v2, v12);    
01097   err = setVertexOnEdge(v2,v0, v20);  
01098 
01099   vector<iBase_EntityHandle> tconn(3);
01100 
01101   tconn[0] = v0;
01102   tconn[1] = v01;
01103   tconn[2] = v20;
01104   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tconn[0], 3, &facehandle, &status, &err);
01105   insertedFaces.push_back( facehandle );
01106 
01107   tconn[0] = v01;
01108   tconn[1] = v1;
01109   tconn[2] = v12;
01110   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tconn[0], 3, &facehandle, &status, &err);
01111   insertedFaces.push_back( facehandle );
01112 
01113   tconn[0] = v12;
01114   tconn[1] = v2;
01115   tconn[2] = v20;
01116   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tconn[0], 3, &facehandle, &status, &err);
01117   insertedFaces.push_back( facehandle );
01118 
01119   tconn[0] = v01;
01120   tconn[1] = v12;
01121   tconn[2] = v20;
01122   iMesh_createEnt(mesh, iMesh_TRIANGLE, &tconn[0], 3, &facehandle, &status, &err);
01123   insertedFaces.push_back( facehandle );
01124 
01125   int val = 1;
01126   iMesh_setIntData(mesh, oldface, remove_tag, val, &err);
01127 
01128   return 0;
01129 }
01130 
01132 
01133 int Refine2D14::atomicOp(const iBase_EntityHandle &oldface)
01134 {
01135    int err, face_topo;
01136    iMesh_getEntTopo(mesh, oldface, &face_topo, &err);
01137    
01138    switch( face_topo ) 
01139    {
01140      case iMesh_TRIANGLE:
01141          refine_tri(oldface);
01142          break;
01143      case iMesh_QUADRILATERAL:
01144          refine_quad(oldface);
01145          break;
01146      default:
01147          cout << "Warning: Element not supported for refinement " << endl;
01148          break;
01149   }
01150   return 0;
01151 }
01152 
01154 
01155 int Refine2D14 ::execute() 
01156 {
01157   int err;
01158   initialize();
01159 
01160   SimpleArray<iBase_EntityHandle> faces;
01161   iMesh_getEntities( mesh, rootSet, iBase_FACE, iMesh_ALL_TOPOLOGIES, 
01162                        ARRAY_INOUT(faces), &err);
01163 
01164   for( size_t i = 0; i < faces.size(); i++)  
01165        atomicOp( faces[i]  );
01166 
01167   ConsistencyRefine2D refine;
01168 
01169   refine.setMesh(mesh);
01170   refine.execute();
01171 
01172   MeshRefine2D::finalize();
01173 
01174   return 0;
01175 }
01176 
01178 
01179 int GradeRefine2D :: initialize()
01180 {
01181    return 0;
01182 }
01183 
01185 
01186 int GradeRefine2D :: atomicOp( const iBase_EntityHandle &apexVertex)
01187 {
01188 /*
01189     SimpleArray<iBase_EntityHandle> vneighs;
01190     iMesh_getEntAdj(mesh, apexVertex, iBase_VERTEX, ARRAY_INOUT(vneighs), &err);
01191 
01192     size_t numNeighs = vneighs.size();
01193 
01194     if( numNeighs == 0 ) return 0;
01195 
01196     vector<double> elen;
01197     elen.resize( numNeighs);
01198 
01199     for( int i = 0; i < numNeighs; i++) 
01200          elen[i] = length( vertex, vneighs[i] );
01201 
01202     sort( elen.begin(), elen.end() );
01203 
01204     double median_value = elen[numNeighs/2];
01205 
01206     for( int i = 0; i < numNeighs; i++) {
01207         if( elen[i] > 0.90*median_value) 
01208             setVertexOnEdge( vertex, vneighs[i]);
01209     }
01210 */
01211     return 1;
01212 }
01214 
01215 int GradeRefine2D :: finalize()
01216 {  
01217   return 0; 
01218 }
01219 
01221 
01222 int GradeRefine2D :: execute()
01223 {
01224 /*
01225     SimpleArray<iBase_EntityHandle> nodes;
01226     iMesh_getEntities( mesh, rootSet, iBase_VERTEX, iMesh_ALL_TOPOLOGIES, 
01227                        ARRAY_INOUT(nodes), &err);
01228 
01229     initialize();
01230 
01231     for( int i = 0; i < nodes.size(); i++) 
01232          atomicOp( nodes[i] );
01233 
01234     ConsistencyRefine2D refine(mesh, edgemap);
01235     refine.execute();
01236 */
01237 
01238     finalize();
01239 }
01240 
01242 
01243 #ifdef TEST_MESHREFINE
01244 int main(int argc, char **argv)
01245 {
01246     iMesh_Instance mesh = read_off_file( "model.off");
01247 
01248 //  LongestEdgeRefine2D  meshrefine;
01249     Refine2D14  meshrefine;
01250     meshrefine.setBoundarySplitFlag(0);
01251     meshrefine.setMesh( mesh );
01252     meshrefine.execute();
01253 
01254     write_off_file( mesh, "refine.off");
01255     
01256     return 0;
01257 }
01258 #endif
01259 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines