MeshKit  1.0
Doublet.cpp
Go to the documentation of this file.
00001 #include "meshkit/QuadCleanUp.hpp"
00002 
00003 using namespace Jaal;
00004 
00006 
00007 void
00008 Jaal::set_doublet_tag(Mesh *mesh, const string &name)
00009 {
00010      int relexist2 = mesh->build_relations(0, 2);
00011 
00012      mesh->search_boundary();
00013 
00014      size_t numnodes = mesh->getSize(0);
00015      for (size_t i = 0; i < numnodes; i++) {
00016           Vertex *vertex = mesh->getNodeAt(i);
00017           if( vertex->isActive() ) {
00018                if (QuadCleanUp::isDoublet(vertex)) {
00019                     vertex->setAttribute(name, 1);
00020                } else
00021                     vertex->setAttribute(name, 0);
00022           }
00023      }
00024 
00025      if (!relexist2)
00026           mesh->clear_relations(0, 2);
00027 }
00028 
00030 
00031 bool
00032 Doublet::isSafe() const
00033 {
00034      assert(vertex);
00035 
00036      if (vertex->isRemoved()) return 0;
00037 
00038      FaceSequence apexfaces;
00039      vertex->getRelations( apexfaces );
00040 
00041      int nSize = apexfaces.size();
00042      for (int i = 0; i < nSize; i++) {
00043           Face *face = apexfaces[i];
00044           if (face->isRemoved()) return 0;
00045           if (face->isVisited()) return 0;
00046      }
00047      return 1;
00048 }
00049 
00051 
00052 void
00053 Doublet::makeShield()
00054 {
00055      FaceSequence apexfaces;
00056      vertex->getRelations( apexfaces );
00057 
00058      size_t nSize = apexfaces.size();
00059      for (size_t i = 0; i < nSize; i++) {
00060           Face *face = apexfaces[i];
00061           face->setVisitMark(1);
00062      }
00063 }
00064 
00066 
00067 vector<Doublet>
00068 QuadCleanUp::search_interior_doublets()
00069 {
00070      //
00072      // An interior doublet is a vertex, which is shared by two face neighbours.
00073      // They are undesirables in the quadmesh as it would mean the angle is 180
00074      // between some adjacent edges...
00075      //
00077 
00078      size_t numnodes = mesh->getSize(0);
00079 
00080      mesh->search_boundary();
00081 
00082      assert(mesh->getAdjTable(0, 2));
00083 
00084      size_t numfaces = mesh->getSize(2);
00085      for (size_t i = 0; i < numfaces; i++) {
00086           Face *face = mesh->getFaceAt(i);
00087           face->setVisitMark(0);
00088      }
00089 
00090      vDoublets.clear();
00091 
00092      for (size_t i = 0; i < numnodes; i++) {
00093           Vertex *vertex = mesh->getNodeAt(i);
00094           if( !vertex->isRemoved() ) {
00095                if (isDoublet(vertex)) {
00096                     Doublet newdoublet(mesh, vertex);
00097                     if (newdoublet.isSafe()) {
00098                          newdoublet.makeShield();
00099                          vDoublets.push_back(newdoublet);
00100                     }
00101                }
00102           }
00103      }
00104 
00105      return vDoublets;
00106 }
00107 
00109 
00110 Vertex*
00111 QuadCleanUp::insert_doublet(Face *face, Vertex *v0, Vertex *v2)
00112 {
00113      //Create new vertex at the center of (v0,v2)
00114      Point3D p3d;
00115      Vertex::mid_point(v0, v2, p3d);
00116 
00117      Vertex *doublet = Vertex::newObject();
00118      doublet->setXYZCoords(p3d);
00119 
00120      int pos = face->getPosOf( v0 );
00121      if( pos < 0) return NULL;
00122 
00123      assert( v2 == face->getNodeAt( (pos+2) ) );
00124 
00125      Vertex *o1 = face->getNodeAt( pos+1 );
00126      Vertex *o2 = face->getNodeAt( pos+3 );
00127 
00128      //  Creating a doublet in the mesh changes:
00129      //  (1)  insert new node
00130      //  (2)  one old face is removed
00131      //  (3)  two new faces inserted.
00132      //
00133 
00134      NodeSequence connect(4);
00135 
00136      mesh->addNode(doublet);
00137      connect[0] = doublet;
00138      connect[1] = v0;
00139      connect[2] = o1;
00140      connect[3] = v2;
00141 
00142      Face *newquad1 = Face::newObject();
00143      newquad1->setNodes(connect);
00144      mesh->addFace(newquad1);
00145 
00146      connect[0] = doublet;
00147      connect[1] = v2;
00148      connect[2] = o2;
00149      connect[3] = v0;
00150 
00151      Face *newquad2 = Face::newObject();
00152      newquad2->setNodes(connect);
00153      mesh->addFace(newquad2);
00154 
00155      mesh->remove( face );
00156 
00157      return doublet;
00158 }
00159 
00161 
00162 Vertex*
00163 QuadCleanUp::insert_doublet(Face *face)
00164 {
00165      if (face->getSize(0) != 4) return NULL;
00166      Vertex *v0 = face->getNodeAt(0);
00167      Vertex *v2 = face->getNodeAt(2);
00168      return insert_doublet(face, v0, v2);
00169 }
00170 
00172 
00173 Vertex*
00174 QuadCleanUp::insert_boundary_doublet(Face *face)
00175 {
00177      //
00178      //                                  X  ( Internal Vertex)
00179      //                               .      .
00180      //                            .           .
00181      //                           .               .
00182      //                         .                   .
00183      //               ********X...........X...........X************************
00184      //        Boundary                 Singlet                Boundary
00185      //
00186      //   There is one quad on the boundary with three nodes on the boundary and
00187      //   one internal nodes. In order to remove the singlet node, we artifically
00188      //   create one doublet between the internal node and the singlet node.
00190 
00191      if (!face->isBoundary())
00192           return NULL;
00193 
00194      if (face->getSize(0) != 4)
00195           return NULL;
00196 
00197      int ncount = 0;
00198      for (int i = 0; i < 4; i++) {
00199           Vertex *v = face->getNodeAt(i);
00200           ncount += v->isBoundary();
00201      }
00202 
00203      if (ncount != 3)
00204           return NULL;
00205 
00206      Vertex *v0 = NULL, *v2 = NULL;
00207      for (int i = 0; i < 4; i++) {
00208           Vertex *v = face->getNodeAt(i);
00209           if (!v->isBoundary()) {
00210                v0 = face->getNodeAt((i + 0) % 4);
00211                v2 = face->getNodeAt((i + 2) % 4);
00212           }
00213      }
00214 
00215      return insert_doublet(face, v0, v2);
00216 }
00217 
00219 
00220 int
00221 Doublet::remove()
00222 {
00223      if (vertex->isRemoved() || vertex->isBoundary()) return 1;
00224 
00225      FaceSequence neighs;
00226      vertex->getRelations( neighs );
00227      assert(neighs.size() > 0);
00228 
00229      if (neighs.size() != 2) return 1;
00230 
00231      Face *face0 = neighs[0];
00232      Face *face1 = neighs[1];
00233 
00234      assert(face0 != face1 );
00235 
00236      if (face0->isRemoved() || face1->isRemoved() )
00237           return 1;
00238 
00239      Vertex *d1 = NULL, *d2 = NULL, *o1 = NULL, *o2 = NULL;
00240 
00241      NodeSequence connect = face0->getNodes();
00242      if (connect.size() != 4)
00243           return 2;
00244 
00245      int nC = connect.size();
00246      for (int i = 0; i < nC; i++) {
00247           if (connect[i] == vertex) {
00248                d1 = connect[(i + 1) % 4];
00249                o1 = connect[(i + 2) % 4];
00250                d2 = connect[(i + 3) % 4];
00251                break;
00252           }
00253      }
00254 
00255      connect = face1->getNodes();
00256      if (connect.size() != 4)
00257           return 2;
00258 
00259      nC = connect.size();
00260      for (int i = 0; i < nC; i++) {
00261           if (connect[i] == vertex) {
00262                o2 = connect[(i + 2) % 4];
00263                break;
00264           }
00265      }
00266 
00267      assert(d1);
00268      assert(d2);
00269      assert(o1);
00270      assert(o2);
00271      assert(o1 != o2);
00272 
00273      // Change the connectivity of face0 and face1 is removed.
00274      mesh->deactivate( face0 );
00275      mesh->deactivate( face1 );
00276 
00277      connect[0] = d1;
00278      connect[1] = o1;
00279      connect[2] = d2;
00280      connect[3] = o2;
00281 
00282      // Reuse one of the face
00283      face0->setNodes( connect );
00284      mesh->reactivate( face0 );
00285 
00286      // Other face and doublet die.
00287      mesh->remove( face1  );
00288      mesh->remove( vertex );
00289 
00290      return 0;
00291 }
00292 
00294 
00295 int
00296 QuadCleanUp::remove_doublets_once()
00297 {
00298      search_interior_doublets();
00299 
00300      int ncount = 0;
00301      size_t nSize = vDoublets.size();
00302 
00303      for (size_t i = 0; i < nSize; i++) {
00304           int err = vDoublets[i].remove();
00305           if (!err) ncount++;
00306      }
00307 
00308 
00309      return ncount;
00310 }
00311 
00313 
00314 int
00315 QuadCleanUp::remove_interior_doublets()
00316 {
00317      StopWatch swatch;
00318      swatch.start();
00319 
00320      int relexist2 = mesh->build_relations(0, 2);
00321 
00322      mesh->search_boundary();
00323 
00324      // It is possible that removal of doublet may create singlets on the boundary
00325      // so, it is better to call doublets first and then call to singlet removal next.
00326      int ncount1 = 0;
00327      while (1) {
00328           int ncount2 = remove_doublets_once();
00329           if (ncount2 == 0) break;
00330           ncount1 += ncount2;
00331      }
00332 
00333      swatch.stop();
00334 
00335      if (!relexist2)
00336           mesh->clear_relations(0, 2);
00337 
00338      cout << "Info: #Doublets removed : " << ncount1 << endl;
00339      cout << "Info:  Doublet removal execution time : " << swatch.getSeconds() << endl;
00340 
00341      return ncount1;
00342 }
00343 
00345 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines