MeshKit  1.0
SwapQuadEdges.cpp
Go to the documentation of this file.
00001 #include "meshkit/SwapEdges.hpp"
00002 #include "meshkit/QuadCleanUp.hpp"
00003 
00004 using namespace Jaal;
00006 
00007 int
00008 QuadEdge::commit()
00009 {
00010      if (newFaces.empty()) return 1;
00011 
00012      int numfaces = newFaces.size();
00013      int numnodes = newNodes.size();
00014 
00015      for (int i = 0; i < numfaces; i++)
00016           if (newFaces[i] == NULL) return 2;
00017 
00018      if (adjFaces[0]->isRemoved() || adjFaces[1]->isRemoved()) {
00019           for (int i = 0; i < numnodes; i++)
00020                delete newNodes[i];
00021           newNodes.clear();
00022 
00023           for (int i = 0; i < numfaces; i++)
00024                delete newFaces[i];
00025           newFaces.clear();
00026           return 3;
00027      }
00028 
00029      Vertex *vertex;
00030      for (int i = 0; i < 4; i++) {
00031           vertex = adjFaces[0]->getNodeAt(i);
00032           if (vertex->isVisited()) return 4;
00033 
00034           vertex = adjFaces[1]->getNodeAt(i);
00035           if (vertex->isVisited()) return 4;
00036      }
00037 
00038      mesh->remove( adjFaces[0] );
00039      mesh->remove( adjFaces[1] );
00040 
00041      assert(mesh);
00042 
00043      for (int i = 0; i < numnodes; i++)
00044           mesh->addNode(newNodes[i]);
00045 
00046      for (int i = 0; i < numfaces; i++)
00047           mesh->addFace(newFaces[i]);
00048 
00049      for (int i = 0; i < 4; i++) {
00050           vertex = adjFaces[0]->getNodeAt(i);
00051           vertex->setVisitMark(1);
00052           vertex = adjFaces[1]->getNodeAt(i);
00053           vertex->setVisitMark(1);
00054      }
00055 
00056      return 0;
00057 }
00058 
00060 
00061 bool
00062 SwapQuadEdge::is_topologically_valid_swap(int d1, int d2, int d3, int d4)
00063 {
00064      if (d1 < 4 || d2 < 4) return 0;
00065      if ((d1 > 4) && (d2 > 4) && (d3 == 3) && (d4 == 3)) return 1;
00066      if ((d1 == 5) && (d2 == 5) && (d3 == 3) && (d4 == 4)) return 1;
00067      if ((d1 == 5) && (d2 == 5) && (d3 == 4) && (d4 == 3)) return 1;
00068      if (max(d1, d2) > max(d3 + 1, d4 + 1)) return 1;
00069      return 0;
00070 }
00072 
00073 int
00074 SwapQuadEdge::getPosOf(const Vertex *vertex) const
00075 {
00076      for (int i = 0; i < 6; i++)
00077           if (bound_nodes[i] == vertex) return i;
00078 
00079      return -1;
00080 }
00082 
00083 int
00084 SwapQuadEdge::build_boundary()
00085 {
00086      assert(mesh);
00087 
00088      assert(connect[0]);
00089      assert(connect[1]);
00090 
00091      adjFaces.resize(2);
00092      adjFaces[0] = NULL;
00093      adjFaces[1] = NULL;
00094 
00095      assert(mesh->getAdjTable(0, 2));
00096 
00097      // Since the degree of each node of an existing edge decreases by
00098      // one, by restricting the swapping to degree greater than 3 will
00099      // ensure that no doublet is created.
00100      //
00101      int nsize1 = connect[0]->getNumRelations(2);
00102 
00103      if (connect[0]->isBoundary()) {
00104           if (nsize1 < 3) return 1;
00105      } else {
00106           if (nsize1 < 4) return 1;
00107      }
00108 
00109      int nsize2 = connect[1]->getNumRelations(2);
00110 
00111      if (connect[1]->isBoundary()) {
00112           if (nsize2 < 3) return 1;
00113      } else {
00114           if (nsize2 < 4) return 1;
00115      }
00116 
00117      // Possibility of doublet creation is ruled out now..
00118      FaceSequence nghs;
00119      Mesh::getRelations112(connect[0], connect[1], nghs);
00120 
00121      if (nghs.size() != 2) return 1;
00122 
00123      if (firstFace == NULL) {
00124           adjFaces[0] = nghs[0];
00125           adjFaces[1] = nghs[1];
00126      }
00127 
00128      if (nghs[0] == firstFace) {
00129           adjFaces[0] = nghs[0];
00130           adjFaces[1] = nghs[1];
00131      }
00132 
00133      if (nghs[1] == firstFace) {
00134           adjFaces[0] = nghs[1];
00135           adjFaces[1] = nghs[0];
00136      }
00137 
00138      int err = get_boundary_nodes_chain();
00139 
00140      if (err) return 2;
00141 
00142      return 0;
00143 }
00144 
00146 
00147 int
00148 SwapQuadEdge::get_boundary_nodes_chain()
00149 {
00150      vector<Edge> bndedges;
00151      bndedges.reserve(6);
00152 
00153      Edge sharededge(connect[0], connect[1]);
00154      for (int i = 0; i < 2; i++) {
00155           for (int j = 0; j < 4; j++) {
00156                Vertex *vf0 = adjFaces[i]->getNodeAt(j + 0);
00157                Vertex *vf1 = adjFaces[i]->getNodeAt(j + 1);
00158                Edge edge(vf0, vf1);
00159                if (!edge.isSameAs(sharededge))
00160                     bndedges.push_back(edge);
00161           }
00162      }
00163 
00164      Mesh::make_chain(bndedges);
00165      Mesh::rotate_chain(bndedges, connect[0]);
00166      bound_nodes = Mesh::chain_nodes(bndedges);
00167      assert(bound_nodes.size() == 6);
00168      return 0;
00169 }
00170 
00172 
00173 void SwapQuadEdge::update_front()
00174 {
00175      cout << "Update front : " << endl;
00176      exit(0);
00177 
00178      /*
00179           size_t nSize;
00180           if (check_fronts) {
00181                NodeSequence vneighs;
00182                while (1) {
00183                     int values_updated = 0;
00184                     for (int i = 0; i < 6; i++) {
00185                          bound_nodes[i]->getRelations( vneighs );
00186                          int minid = bound_nodes[i]->getLayerID();
00187                          nSize   = vneighs.size();
00188                          for (size_t j = 0; j < nSize; j++)
00189                               minid = min(minid, vneighs[j]->getLayerID() + 1);
00190 
00191                          if (bound_nodes[i]->getLayerID() != minid) {
00192                               bound_nodes[i]->setLayerID(minid);
00193                               values_updated = 1;
00194                          }
00195                     }
00196                     if (!values_updated) break;
00197                }
00198           }
00199      */
00200 }
00201 
00203 
00204 void
00205 SwapQuadEdge::backup()
00206 {
00207      bkp_data.diagonalConnect[0] = connect[0];
00208      bkp_data.diagonalConnect[1] = connect[1];
00209      bkp_data.face1Connect = adjFaces[0]->getNodes();
00210      bkp_data.face2Connect = adjFaces[1]->getNodes();
00211      for (int i = 0; i < 6; i++) {
00212           Vertex *v = bound_nodes[i];
00213           bkp_data.pCoords[v] = v->getXYZCoords();
00214      }
00215 }
00217 
00218 void
00219 SwapQuadEdge::rollback()
00220 {
00221      // We need to manual delete the relation
00222      connect[0]->removeRelation(connect[1]);
00223      connect[1]->removeRelation(connect[0]);
00224 
00225      // remove all relations by deactivating the elements ...
00226      mesh->deactivate(adjFaces[0]);
00227      mesh->deactivate(adjFaces[1]);
00228 
00229      // Change the connectivity of the two quads and reactivate the element...
00230      adjFaces[0]->setNodes(bkp_data.face1Connect);
00231      adjFaces[1]->setNodes(bkp_data.face2Connect);
00232 
00233      mesh->reactivate(adjFaces[0]);
00234      mesh->reactivate(adjFaces[1]);
00235 
00236      connect[0] = bkp_data.diagonalConnect[0];
00237      connect[1] = bkp_data.diagonalConnect[1];
00238 
00239      for (int i = 0; i < 6; i++) {
00240           Vertex *v = bound_nodes[i];
00241           v->setXYZCoords(bkp_data.pCoords[v]);
00242      }
00243 
00244      assert(!adjFaces[0]->isRemoved());
00245      assert(!adjFaces[1]->isRemoved());
00246 
00247 //     update_front();
00248 }
00250 
00251 int
00252 SwapQuadEdge::make_new_diagonal_at(int pos, bool bound_check)
00253 {
00254      //
00255      // Given closednodes ( exactly six nodes) in an order ( clockwise or
00256      // counter-clockwise) and create two quads having common diagonal between
00257      // (pos, pos+3).
00258      //
00259      assert(pos >= 0 && pos < 6);
00260 
00261      if (bound_check) {
00262           if (bound_nodes[(pos + 0) % 6]->isBoundary()) return 1;
00263           if (bound_nodes[(pos + 3) % 6]->isBoundary()) return 1;
00264      }
00265 
00266      // First perform full backup of the data structures.
00267      backup();
00268 
00269      // Remove the vertex-vertex relationship.
00270      connect[0]->removeRelation(connect[1]);
00271      connect[1]->removeRelation(connect[0]);
00272 
00273      mesh->deactivate(adjFaces[0]);
00274      mesh->deactivate(adjFaces[1]);
00275 
00276      // Change the connectivity of the two quads.
00277      NodeSequence qConnect(4);
00278      qConnect[0] = bound_nodes[(pos + 0) % 6];
00279      qConnect[1] = bound_nodes[(pos + 1) % 6];
00280      qConnect[2] = bound_nodes[(pos + 2) % 6];
00281      qConnect[3] = bound_nodes[(pos + 3) % 6];
00282      connect[0] = qConnect[0];
00283      adjFaces[0]->setNodes(qConnect);
00284      mesh->reactivate(adjFaces[0]);
00285 
00286      qConnect[0] = bound_nodes[(pos + 3) % 6];
00287      qConnect[1] = bound_nodes[(pos + 4) % 6];
00288      qConnect[2] = bound_nodes[(pos + 5) % 6];
00289      qConnect[3] = bound_nodes[(pos + 6) % 6];
00290      connect[1] = qConnect[0];
00291      adjFaces[1]->setNodes(qConnect);
00292      mesh->reactivate(adjFaces[1]);
00293 
00294      assert(!adjFaces[0]->isRemoved());
00295      assert(!adjFaces[1]->isRemoved());
00296 
00297      int con1 = adjFaces[0]->concaveAt();
00298      int con2 = adjFaces[1]->concaveAt();
00299 //  int con1 = adjFaces[0]->isSimple();
00300 //  int con2 = adjFaces[1]->isSimple();
00301 //  if( con1 && con2 ) return  0;
00302 
00303      if( con1 >= 0 || con2 >=0) {
00304           rollback();
00305           return 1;
00306      }
00307      return 0;
00308 
00309      // The change may produce concave or tangled mesh. Do some local
00310      // Laplacian smoothing at the six nodes only. Hopefully, the elements
00311      // will be acceptable. If not, we need to rollback this operation.
00312 
00313      LaplaceLengthWeight lw;
00314      LaplaceSmoothing lapsmooth(mesh);
00315      lapsmooth.setWeight(&lw);
00316      lapsmooth.setNumIterations(10);
00317      lapsmooth.localized_at(bound_nodes);
00318 
00319      // Look at the neighbors, if some elements is inverted. Here we create
00320      // a set of neighbours, to avoid duplicate checking.
00321      FaceSet neighSet;
00322      FaceSequence vfaces;
00323      for (int i = 0; i < 6; i++) {
00324           Vertex *v = bound_nodes[i];
00325           v->getRelations( vfaces );
00326           size_t nSize = vfaces.size();
00327           for (size_t j = 0; j < nSize; j++)
00328                neighSet.insert(vfaces[j]);
00329      }
00330      assert(!neighSet.empty());
00331 
00332      FaceSet::const_iterator it;
00333      for (it = neighSet.begin(); it != neighSet.end(); ++it) {
00334           Face *face = *it;
00335           int pos = face->concaveAt();
00336           if (pos >= 0) {
00337                Vertex *v = face->getNodeAt(pos);
00338                if (!v->isBoundary()) {
00339                     rollback();
00340                     return 1;
00341                }
00342           }
00343      }
00344      update_front();
00345 
00346      return 0;
00347 }
00348 
00350 
00351 int
00352 SwapQuadEdge::apply_reduce_degree_rule()
00353 {
00354      if (build_boundary() != 0) return 1;
00355 
00356      NodeSequence neighs;
00357 
00358      connect[0]->getRelations( neighs );
00359      int d1 = neighs.size();
00360      if (d1 < 4) return 1;
00361 
00362      connect[1]->getRelations( neighs );
00363      int d2 = neighs.size();
00364      if (d2 < 4) return 1;
00365 
00366      int start_pos = getPosOf(connect[0]);
00367      assert(start_pos >= 0);
00368 
00369      for (int i = 0; i < 2; i++) {
00370           int pos = (start_pos + i + 1) % 6;
00371           Vertex *v0 = bound_nodes[ pos ];
00372           if (v0->isBoundary()) continue;
00373           v0->getRelations( neighs );
00374           int d3 = neighs.size();
00375 
00376           Vertex *v3 = bound_nodes[(pos + 3) % 6];
00377           if (v3->isBoundary()) continue;
00378           v3->getRelations( neighs );
00379           int d4 = neighs.size();
00380 
00381           if (SwapQuadEdge::is_topologically_valid_swap(d1, d2, d3, d4)) {
00382                int err = make_new_diagonal_at(pos);
00383                if (!err) return 0;
00384           }
00385      }
00386      return 1;
00387 }
00388 
00390 
00391 int
00392 SwapQuadEdge::apply_concave_rule()
00393 {
00394      if (build_boundary() != 0) return 1;
00395 
00396      int start_pos = getPosOf(connect[0]);
00397      for (int i = 0; i < 2; i++) {
00398           int err = make_new_diagonal_at((start_pos + 1 + i) % 6);
00399           if (!err) return 0;
00400      }
00401      return 1;
00402 }
00403 
00405 
00406 int
00407 SwapQuadEdge::apply_bound_rule()
00408 {
00409      if (!connect[0]->isBoundary()) return 1;
00410      if (connect[1]->isBoundary()) return 1;
00411 
00412      if (build_boundary() != 0) return 2;
00413 
00414      for (int i = 0; i < 3; i++) {
00415           Vertex *v0 = bound_nodes[(i + 0) % 6];
00416           Vertex *v3 = bound_nodes[(i + 3) % 6];
00417           if ((!v0->isBoundary()) && (!v3->isBoundary())) {
00418                int err = make_new_diagonal_at(i, 0);
00419                if (!err) return 0;
00420           }
00421      }
00422 
00423      return 3;
00424 }
00425 
00427 
00428 int
00429 SwapQuadEdge::apply_advance_front_rule()
00430 {
00431      if (build_boundary() != 0) return 1;
00432 
00433      /*
00434           int layerid = connect[0]->getLayerID();
00435           for (int i = 0; i < 3; i++) {
00436                Vertex *v0 = bound_nodes[(i + 0) % 6];
00437                Vertex *v3 = bound_nodes[(i + 3) % 6];
00438                if ((v0->getLayerID() > layerid) && (v3->getLayerID() > layerid)) {
00439                     int err = make_new_diagonal_at(i);
00440                     if (!err) return 0;
00441                }
00442           }
00443      */
00444      return 3;
00445 }
00446 
00448 
00449 int
00450 SwapQuadEdge::apply_singlet_rule(Vertex *singlet)
00451 {
00452      // The first node must be boundary and the other node must be internal..
00453      assert(connect[0]->isBoundary());
00454 
00455      if (connect[1]->isBoundary()) return 1;
00456 
00457      if (build_boundary() != 0) return 1;
00459      // Objective :: Swap quads common diagonal such that the resulting diagonal
00460      //              contains the Singlet node.
00462      assert(QuadCleanUp::isSinglet(singlet));
00463 
00464      if (adjFaces[0] == NULL || adjFaces[1] == NULL) return 2;
00465 
00466      assert(QuadCleanUp::hasSinglet(adjFaces[0]));
00467 
00468      // Make sure that other face doesn't have a singlet
00469      if (QuadCleanUp::hasSinglet(adjFaces[1])) return 3;
00470 
00471      //Find the position of singlet in the contour ...
00472      int pos = this->getPosOf(singlet);
00473      assert(pos >= 0);
00474 
00475      // Create new quads whose common diagonal must contain the singlet.
00476      return make_new_diagonal_at(pos, 0);
00477 }
00478 
00480 
00481 int
00482 SwapQuadEdge::apply_deficient_rule(Vertex *vertex)
00483 {
00484 
00485 #ifdef NEEDS_RECONSIDER
00486      assert(firstFace);
00487      if (vertex->isBoundary()) return 1;
00488 
00489 
00490      int d1 = connect[0]->getNumRelations( 2 );
00491      int d2 = connect[1]->getNumRelations( 2 );
00492 
00493      // Note 1: Just near the deficient node, there must be atleast 5 nodes in
00494      // order to swap the edge,otherwise, things may not converge, on the
00495      // same level, the deficient node will keep changing the location.
00496      //
00497      // Note 2: The node connect[1] is opposite to the deficient node and
00498      // most probably on the next level. If we keep moving the deficient
00499      // node inside the domain, then it is fine to have its degree to 4, so
00500      // that the deficient node will be generated there.
00501      if( d2 < 4 ) return 1;
00502 
00503      FaceSequence faces;
00504      if (d1 == 4) {
00505           vertex->getRelations( faces );
00506           if (faces.size() != 3) return 1;
00507           int layerid = connect[0]->getLayerID();
00508           Face *f0 = firstFace;
00509           Mesh::getRelations112(connect[0], connect[1], faces);
00510           assert(faces.size() == 2);
00511           Face *f1 = NULL;
00512           if (faces[0] == f0) f1 = faces[1];
00513           if (faces[1] == f0) f1 = faces[0];
00514           assert(f1);
00515 
00516 
00517           int pos = f1->getPosOf(connect[0]);
00518           Vertex *vopp = f1->getNodeAt(pos + 2);
00519           if (connect[1]->getLayerID() > layerid && vopp->getLayerID() > layerid) {
00520                set_no_tags(mesh);
00521                connect[1]->setTag(1);
00522                vopp->setTag(1);
00523                mesh->saveAs("dbg.dat");
00524                Break();
00525                SwapQuadEdge edge(mesh, connect[1], vopp, f1);
00526                int err = edge.apply_deficient_rule(connect[0]);
00527                if (err) return 1;
00528           }
00529      }
00530      exit(0);
00531 
00532      d1 = connect[0]->getNumRelations(2);
00533      d2 = connect[1]->getNumRelations(2);
00534 
00535      // Don't create doublet at two nodes...
00536      if( d1 == 3 || d2 == 3 ) return 1;
00537 
00538      // Having these conditions set, now build the structure.
00539      if (build_boundary() != 0) return 3;
00540 
00541      update_front();
00542 
00543      //Find the position of triplet in the contour ...
00544      int pos = this->getPosOf(vertex);
00545      assert(pos == 1 || pos == 5);
00546 
00547      if (check_fronts) {
00548           int layerid = vertex->getLayerID();
00549           if (bound_nodes[ (pos + 3) % 6]->getLayerID() <= layerid) return 1;
00550 
00551           if( connect[0]->getLayerID() <= layerid && d1 < 5 ) return 1;
00552           if( connect[1]->getLayerID() <= layerid && d2 < 5 ) return 1;
00553      }
00554 
00555      // Create new quads whose common diagonal must contain the singlet.
00556      int err = make_new_diagonal_at(pos);
00557 
00558      /*
00559          if( err  ) {
00560              for (int i = 0; i < 6; i++) {
00561                   NodeSequence vneighs = bound_nodes[i]->getRelations0();
00562                   int minid = bound_nodes[i]->getLayerID();
00563                   for( int j = 0; j < vneighs.size(); j++)
00564                        minid = min( minid, vneighs[j]->getLayerID() );
00565                  assert( bound_nodes[i]->getLayerID() == minid+1);
00566               }
00567          }
00568      */
00569 
00570      return err;
00571 #endif
00572     return 1;
00573 
00574 }
00575 
00577 
00578 int
00579 QuadCleanUp::swap_concave_faces()
00580 {
00581      int relexist0 = mesh->build_relations(0, 0);
00582      int relexist2 = mesh->build_relations(0, 2);
00583 
00584      mesh->search_boundary();
00585 
00586      size_t numfaces = mesh->getSize(2);
00587      for (size_t i = 0; i < numfaces; i++) {
00588           Face *face = mesh->getFaceAt(i);
00589           face->setVisitMark(0);
00590      }
00591 
00592      size_t numnodes = mesh->getSize(0);
00593      for (size_t i = 0; i < numnodes; i++) {
00594           Vertex *vertex = mesh->getNodeAt(i);
00595           vertex->setVisitMark(0);
00596      }
00597 
00598      size_t ncount = 0;
00599      for (size_t iface = 0; iface < numfaces; iface++) {
00600           Face *face = mesh->getFaceAt(iface);
00601           int pos = face->concaveAt();
00602           if (pos >= 0) {
00603                for (int i = 0; i < 2; i++) {
00604                     Vertex *v0 = face->getNodeAt(pos + 1 + i);
00605                     Vertex *v1 = face->getNodeAt(pos + 2 + i);
00606                     SwapQuadEdge edge(mesh, v0, v1);
00607                     int err = edge.apply_concave_rule();
00608                     if (!err) {
00609                          ncount++;
00610                          break;
00611                     }
00612                }
00613           }
00614      }
00615 
00616      if (ncount) {
00617           mesh->prune();
00618           mesh->enumerate(0);
00619           mesh->enumerate(2);
00620           cout << "Info: # of swapped edges " << ncount << endl;
00621      }
00622 
00623      if (!relexist0) mesh->clear_relations(0, 0);
00624      if (!relexist2) mesh->clear_relations(0, 2);
00625 
00626      return ncount;
00627 }
00628 
00630 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines