MeshKit  1.0
AdfrontClean.cpp
Go to the documentation of this file.
00001 #include "meshkit/QuadCleanUp.hpp"
00002 
00003 using namespace Jaal;
00004 
00006 
00007 int
00008 QuadCleanUp::apply_advance_front_singlet_rule(Vertex *singlet)
00009 {
00010     int err;
00011 
00012     int lid;
00013     singlet->getAttribute("Layer", lid);
00014     if ( lid  != 0) return 1;
00015     if (!QuadCleanUp::isSinglet(singlet)) return 1;
00016 
00017     FaceSequence vfaces = singlet->getRelations2();
00018     assert(vfaces.size() == 1);
00019 
00020     Face *face = vfaces[0];
00021     int pos = face->getPosOf(singlet);
00022 
00023     Vertex *v1 = face->getNodeAt(pos + 1);
00024     Vertex *v2 = face->getNodeAt(pos + 2);
00025     Vertex *v3 = face->getNodeAt(pos + 3);
00026 
00027     SwapQuadEdge edge1(mesh, v1, v2);
00028     err = edge1.apply_singlet_rule(singlet);
00029     if (!err) return 0;
00030 
00031     SwapQuadEdge edge2(mesh, v3, v2);
00032     err = edge2.apply_singlet_rule(singlet);
00033     if (!err) return 0;
00034 
00035     return 1;
00036 }
00037 
00039 
00040 int
00041 QuadCleanUp::apply_advance_front_triplet_rule(Vertex *vertex)
00042 {
00043     //****************************************************************************
00044     // Implementation ideas:
00045     //
00046     // This is the case of 5-434 transition, Where 5 degree node is opposite
00047     // to 3 degree node and with first swap, it becomes 4-335 or 4-534 and with
00048     // the second swap, it becomes 3-444 and therefore, the 3 degree node moves
00049     // inside the domain.
00050     //
00051     // Having degree 5 is a somewhat strict condition, because degree four is
00052     // sufficient, but it will create doublet, and as per our design goals, we
00053     // do not want to introduce any new doublet in the mesh. In future, we may
00054     // change the conditions, if that is useful.
00055     //
00056     //****************************************************************************
00057     if (vertex->isBoundary()) return 1;
00058 
00059     NodeSequence vnodes = vertex->getRelations0();
00060     if (vnodes.size() != 3) return 1;
00061 
00062     int layerid = vertex->getLayerID();
00063     FaceSequence vfaces = vertex->getRelations2();
00064 
00065     int pos;
00066     // Find out which face which contains the higher level node than the given node.
00067     Face *face = NULL;
00068     for (int j = 0; j < 3; j++) {
00069         Face *f = vfaces[j];
00070         pos = f->getPosOf( vertex );
00071         Vertex *v = f->getNodeAt(  pos+2 );
00072         if( v->getLayerID() > layerid ) {
00073             face = f;
00074            break;
00075         }
00076     }
00077 
00078     if( face == NULL ) return 1;
00079 
00080     pos = face->getPosOf(vertex);
00081     Vertex *v1 = face->getNodeAt((pos + 1) % 4);
00082     Vertex *v2 = face->getNodeAt((pos + 2) % 4);
00083     Vertex *v3 = face->getNodeAt((pos + 3) % 4);
00084 
00085      // If v1 and v3 are at the lower level, they are already processed.
00086      // v2 must always be higher.
00087      if (v1->getLayerID() <  layerid   || 
00088          v2->getLayerID() <= layerid   ||
00089          v3->getLayerID() <  layerid ) return 1;
00090 
00091      int d1 = v1->getRelations2().size();
00092      int d2 = v2->getRelations2().size();
00093      int d3 = v3->getRelations2().size();
00094 
00095      // Avoids creating doublet at v2
00096      if( d2 < 4 ) return 1;
00097 
00098      // Avoids creating doublet at v1 or v3 
00099      if (max(d1, d3) < 4 ) return 1;
00100 
00101      SwapQuadEdge edge1(mesh, v1, v2, face);
00102      edge1.modify_fronts(1);
00103      int err1 = edge1.apply_deficient_rule(vertex);
00104      if (!err1) return 0;
00105 
00106      SwapQuadEdge edge2(mesh, v3, v2, face);
00107      edge2.modify_fronts(1);
00108      int err2 = edge2.apply_deficient_rule(vertex);
00109      if (!err2 ) return 0;
00110 
00111     return 1;
00112 }
00113 
00115 
00116 int
00117 QuadCleanUp::apply_advance_front_bridge_rule(Vertex *v0, Vertex *v1)
00118 {
00119     int layerid = v0->getLayerID();
00120     if (v1->getLayerID() != layerid) return 1;
00121 
00122     // Although it is checked once again, but it is essential.
00123     FaceSequence adjFaces;
00124     adjFaces = v0->getRelations2();
00125     if (adjFaces.size() != 3) return 2;
00126 
00127     adjFaces = v1->getRelations2();
00128     if (adjFaces.size() != 3) return 3;
00129 
00130     // Our assumption is that all the edges are simple.
00131     Mesh::getRelations112(v0, v1, adjFaces);
00132     if (adjFaces.size() != 2) return 4;
00133 
00134     if (adjFaces[0]->getLayerID() == adjFaces[1]->getLayerID()) return 5;
00135 
00136     // Check for the internal face.
00137     Face *internal_face = NULL;
00138     internal_face = adjFaces[0]->getLayerID() > adjFaces[1]->getLayerID() ? adjFaces[0] : adjFaces[1];
00139 
00140     // Swap the opposite edge:
00141     Vertex *v2, *v3;
00142     Face::opposite_nodes(internal_face, v0, v1, v2, v3);
00143 
00144     // Try swapping at the first node.
00145     int err;
00146     SwapQuadEdge edge1(mesh, v2, v3, internal_face);
00147     err = edge1.apply_deficient_rule(v0);
00148     if (!err) return 0;
00149 
00150     // Try swapping at the second node.
00151     SwapQuadEdge edge2(mesh, v2, v3, internal_face);
00152     err = edge2.apply_deficient_rule(v1);
00153     if (!err) return 0;
00154 
00155     return 1;
00156 
00157 }
00158 
00160 
00161 int
00162 QuadCleanUp::apply_advance_front_excess_rule(Vertex *vertex)
00163 {
00164     if (vertex->isBoundary()) return 1;
00165 
00166     int layerid = vertex->getLayerID();
00167     NodeSequence vneighs = vertex->getRelations0();
00168     int degree = vneighs.size();
00169 
00170     if (degree < 5) return 1;
00171 
00172     int ncount = 0;
00173     for (int k = 0; k < degree; k++)
00174         if (vneighs[k]->getLayerID() > layerid) ncount++;
00175 
00176     if (ncount < 2) return 1;
00177 
00178     for (int k = 0; k < degree; k++)
00179     {
00180         if (vneighs[k]->getLayerID() > layerid)
00181         {
00182             SwapQuadEdge edge(mesh, vertex, vneighs[k]);
00183             int err = edge.apply_advance_front_rule();
00184             if (!err) return 0;
00185         }
00186     }
00187 
00188     return 1;
00189 }
00190 
00192 
00193 int
00194 QuadCleanUp::advance_front_edges_swap_once(int layerid)
00195 {
00196     if (mesh->getAdjTable(0, 0))
00197         mesh->clear_relations(0, 0);
00198 
00199     if (mesh->getAdjTable(0, 2))
00200         mesh->clear_relations(0, 2);
00201 
00202     int relexist0 = mesh->build_relations(0, 0);
00203     int relexist2 = mesh->build_relations(0, 2);
00204 
00205     //
00206     // The input mesh should be doublet free and the output mesh will
00207     // be doublet free.
00208     //
00209     vector<Doublet> doublets = search_interior_doublets();
00210     assert(doublets.empty());
00211 
00212     // If the boundary is unknown ...
00213     if (layerid == 0) mesh->search_boundary();
00214 
00215     // We need atleast two layer to work with.
00216     size_t numnodes = mesh->getSize(0);
00217 
00218     // Check how many iiregular nodes in the present layer...
00219     size_t num_irregular_nodes_before = 0;
00220     for (size_t i = 0; i < numnodes; i++)
00221     {
00222         Vertex *vertex = mesh->getNodeAt(i);
00223         if (vertex->getLayerID() == layerid)
00224         {
00225             if (vertex->getRelations2().size() != 4) 
00226                 num_irregular_nodes_before++;
00227         }
00228     }
00229 
00230     size_t ncount = 0;
00231 
00232     //
00233     // If the layer is boundary, then highest priority must be given to remove
00234     // the singlets. Many singlets can be removed by swapping the edges. There are
00235     // some cases, where "Swapping" may not remove singlets. There are two such
00236     // scenarios, which are handled by calling "remove_boundary_singlet" method.
00237     //
00238 
00239 #ifdef REMOVE_LATER
00240     for (size_t i = 0; i < numnodes; i++)
00241     {
00242         Vertex *vertex = mesh->getNodeAt(i);
00243         if (vertex->isBoundary())
00244         {
00245             err = apply_advance_front_singlet_rule(vertex);
00246             if (!err) ncount++;
00247         }
00248     }
00249 
00250     // Second priority is given to bridges in the mesh. Removal of bridges
00251     // create a three degree node, which is the next case.
00252 
00253     if (ncount == 0)
00254     {
00255         /*
00256            vector<Edge33> bridges = search_edges33();
00257            for (size_t i = 0; i < bridges.size(); i++){
00258                 Vertex *v0 = bridges[i].getNodeAt(0);
00259                 Vertex *v1 = bridges[i].getNodeAt(1);
00260                 if( v0->getLayerID() == layerid  && v1->getLayerID() == layerid ) {
00261                     err = apply_advance_front_bridge_rule(v0, v1);
00262                     if( !err ) ncount++;
00263                     v0->setTag(3);
00264                     v1->setTag(3);
00265                     mesh->saveAs("dbg.dat");
00266                     exit(0);
00267                 }
00268            }
00269          */
00270     }
00271 #endif
00272 
00273     // Third Priority is given to degree three nodes. if the adjacent nodes are
00274     // regular, then two swaps are required. In the first step, the degree of the
00275     // adjacent node is increased by swapping from outer layer, and then swapping
00276     // is done to make 3 degree node to regular node.
00277 
00278     if (ncount == 0)
00279     {
00280         for (size_t i = 0; i < numnodes; i++)
00281         {
00282             Vertex *vertex = mesh->getNodeAt(i);
00283 
00284             if (vertex->getLayerID() == layerid)
00285             {
00286                 int err = apply_advance_front_triplet_rule(vertex);
00287                 if (!err) ncount++;
00288             }
00289         }
00290     }
00291 
00292     // Final case, this is the most intuitive and general; swapping is done to
00293     // reduce the vertex degree.
00294     if (ncount == 0)
00295     {
00296         for (size_t i = 0; i < numnodes; i++)
00297         {
00298             Vertex *vertex = mesh->getNodeAt(i);
00299             if (vertex->getLayerID() == layerid) {
00300                 int err = apply_advance_front_excess_rule(vertex);
00301                 if (!err) ncount++;
00302             }
00303         }
00304     }
00305 
00306 #ifdef DEBUG
00307     if (ncount)
00308       cout << "Info: Layer " << layerid << " number of layer edges swapped " << ncount << endl;
00309 
00310     if( ncount ) {
00311         set_no_tags(mesh);
00312         size_t num_irregular_nodes_after = 0;
00313         for (size_t i = 0; i < numnodes; i++)
00314         {
00315              Vertex *vertex = mesh->getNodeAt(i);
00316              vertex->setTag(1);
00317              if (vertex->getLayerID() == layerid)
00318              {
00319                  if (vertex->getRelations2().size() != 4) {
00320                      vertex->setTag(2);
00321                      num_irregular_nodes_after++;
00322                  }
00323                  NodeSequence vnghs = vertex->getRelations0();
00324                  int minid = MAXINT;
00325                  for( int j = 0; j < vnghs.size(); j++)
00326                       minid = min( minid, vnghs[j]->getLayerID() );
00327                  vertex->setTag(2);
00328                  mesh->saveAs( "dbg.dat");
00329                  assert( vertex->getLayerID() == minid+1 );
00330             }
00331       }
00332 
00333       cout << "Edges Swaaped " << num_irregular_nodes_before 
00334                                << " " << num_irregular_nodes_after << endl;
00335       assert( num_irregular_nodes_before > num_irregular_nodes_after);
00336       Break();
00337     }
00338 #endif
00339 
00340     if (!relexist0)
00341         mesh->clear_relations(0, 0);
00342 
00343     if (!relexist2)
00344         mesh->clear_relations(0, 2);
00345 
00346     mesh->collect_garbage();
00347 
00348     doublets = search_interior_doublets();
00349     assert(doublets.empty());
00350 
00351     return ncount;
00352 }
00353 
00355 
00356 void
00357 QuadCleanUp::advancing_front_edges_swap()
00358 {
00359     //****************************************************************************
00360     // Implementation ideas:
00361     // In general, irregular nodes are scattered around in the domain ( even after
00362     // initial cleanup operations. With the advance front swapping, we try to
00363     // (1) clean the mesh starting from the boundary (2) cluster irregular nodes
00364     // deep inside the mesh ( far from the boundary ) with the hope that clustering
00365     // will provide more opportunities for cleanup operations ( for example, enable
00366     // diamonds, bridges etc).
00367     //
00368     // There are two side-effects of this procedure.
00369     // 1)  It may create very high valance nodes in the last layers, which we may
00370     //     not be able to clean. ( this is not a stopper, as we can call vertex
00371     //     reduction modules, which will again scatter the irregular nodes.
00372     // 2)  The position of irregular nodes may depends on the density of mesh.
00373     //
00374     //
00375     //***************************************************************************
00376     int nfronts = mesh->setWavefront(0);
00377 
00378 #ifdef DEBUG
00379     cout << "Info:  Before Advance front swapping " << endl;
00380     map<int,size_t> nodedegree;
00381     map<int,size_t>::const_iterator it;
00382     for (int ilayer = 1; ilayer < nfronts;  ilayer++) {
00383          nodedegree.clear();
00384          for (size_t i = 0; i < numnodes; i++)
00385          {
00386               Vertex *vertex = mesh->getNodeAt(i);
00387               if (vertex->getLayerID() == ilayer) {
00388                   int nd = vertex->getRelations2().size();
00389                   nodedegree[nd]++;
00390               }
00391         }
00392         cout << "In Layer " <<   ilayer << " Degree Distribution" << endl;
00393         for( it = nodedegree.begin(); it != nodedegree.end(); ++it)
00394              cout << it->first <<  "  " << it->second << endl;
00395     }
00396 #endif
00397 
00399 
00400     for (int ilayer = 1; ilayer < nfronts;  ilayer++)
00401         advance_front_edges_swap_once(ilayer);
00402 
00404 
00405 #ifdef DEBUG
00406     nfronts = mesh->setWavefront(0);
00407     cout << "Info:  After Advance front swapping " << endl;
00408     for (int ilayer = 1; ilayer < nfronts;  ilayer++) {
00409          nodedegree.clear();
00410          for (size_t i = 0; i < numnodes; i++)
00411          {
00412               Vertex *vertex = mesh->getNodeAt(i);
00413               if (vertex->getLayerID() == ilayer) {
00414                   int nd = vertex->getRelations2().size();
00415                   nodedegree[nd]++;
00416               }
00417         }
00418         cout << "In Layer " <<   ilayer << " Degree Distribution" << endl;
00419         for( it = nodedegree.begin(); it != nodedegree.end(); ++it)
00420              cout << it->first <<  "  " << it->second << endl;
00421     }
00422 #endif
00423 
00424 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines