MeshKit  1.0
OneDefectRemeshing.cpp
Go to the documentation of this file.
00001 #include "meshkit/QuadCleanUp.hpp"
00002 
00003 using namespace Jaal;
00004 
00005 size_t OneDefectPatch :: MAX_FACES_ALLOWED = 500;
00006 
00007 size_t OneDefectPatch :: num_boundaries  = 0;
00008 size_t OneDefectPatch :: num_3_patches   = 0;
00009 size_t OneDefectPatch :: num_4_patches   = 0;
00010 size_t OneDefectPatch :: num_5_patches   = 0;
00011 size_t OneDefectPatch :: disk_remeshable = 0;
00012 double OneDefectPatch :: exec_time       = 0.0;
00013 
00015 int TriRemeshTemplate::remesh(Mesh *msh,
00016                               NodeSequence &anodes,
00017                               NodeSequence &bnodes,
00018                               NodeSequence &cnodes, int *partition)
00019 {
00020      mesh = msh;
00021 
00022      // First thing to do is to clear the existing record.
00023      newnodes.clear();
00024      newfaces.clear();
00025 
00026      // We need atleast three nodes on each side ...
00027      if (anodes.size() < 3) return 1;
00028      if (bnodes.size() < 3) return 1;
00029      if (cnodes.size() < 3) return 1;
00030 
00031      if (partition == NULL) {
00032           segments[0] = anodes.size() - 1;
00033           segments[1] = bnodes.size() - 1;
00034           segments[2] = cnodes.size() - 1;
00035 
00036           if (!Face::is_3_sided_convex_loop_quad_meshable(segments, partSegments))
00037                return 1;
00038      } else {
00039           for (int i = 0; i < 6; i++)
00040                partSegments[i] = partition[i];
00041      }
00042 
00043 
00044      int err;
00045      if (anodes.back() == bnodes.back()) {
00046           reverse(bnodes.begin(), bnodes.end());
00047           swap(partSegments[2], partSegments[3]);
00048      }
00049 
00050      if (anodes.front() == cnodes.front()) {
00051           reverse(cnodes.begin(), cnodes.end());
00052           swap(partSegments[4], partSegments[5]);
00053      }
00054 
00055      // Ensure that input is closed loop of three sides ...
00056      assert(anodes.front() == cnodes.back());
00057      assert(bnodes.front() == anodes.back());
00058      assert(cnodes.front() == bnodes.back());
00059 
00060      // Ensure that segments are valid ...
00061      assert((int)anodes.size() == partSegments[0] + partSegments[1] + 1);
00062      assert((int)bnodes.size() == partSegments[2] + partSegments[3] + 1);
00063      assert((int)cnodes.size() == partSegments[4] + partSegments[5] + 1);
00064 
00065      // Split Side "A" nodes into a1nodes and b1nodes.
00066      int a1 = partSegments[0];
00067      err = split_stl_vector(anodes, a1 + 1, a1nodes, b1nodes);
00068      if (err) return 1;
00069 
00070      // Split Side "B" nodes into a2nodes and b2nodes.
00071      int a2 = partSegments[2];
00072      err = split_stl_vector(bnodes, a2 + 1, a2nodes, b2nodes);
00073      if (err) return 1;
00074 
00075      // Split Side "C" nodes into a3nodes and b3nodes.
00076      int a3 = partSegments[4];
00077      err = split_stl_vector(cnodes, a3 + 1, a3nodes, b3nodes);
00078      if (err) return 1;
00079 
00080      // Splitting nodes on each side
00081      Vertex *ca = a1nodes.back();
00082      Vertex *cb = a2nodes.back();
00083      Vertex *cc = a3nodes.back();
00084 
00085      // We are not if the centroid will be in the patch, take the risk.
00086      Vertex *co = Face::centroid(ca, cb, cc);
00087 
00088      mesh->addNode(co);
00089      newnodes.push_back(co);
00090 
00091      linear_interpolation(mesh, co, ca, a2 + 1, oa_nodes);
00092      size_t nSize = oa_nodes.size();
00093      for (size_t i = 1; i < nSize - 1; i++)
00094           newnodes.push_back(oa_nodes[i]);
00095 
00096      linear_interpolation(mesh, co, cb, a3 + 1, ob_nodes);
00097      nSize = ob_nodes.size();
00098      for (size_t i = 1; i < nSize - 1; i++)
00099           newnodes.push_back(ob_nodes[i]);
00100 
00101      linear_interpolation(mesh, co, cc, a1 + 1, oc_nodes);
00102      nSize = oc_nodes.size();
00103      for (size_t i = 1; i < nSize - 1; i++)
00104           newnodes.push_back(oc_nodes[i]);
00105 
00106      err = quadtemplate.remesh(mesh, b2nodes, a3nodes, oc_nodes, ob_nodes);
00107      if (!err) addNewElements( quadtemplate.newnodes, quadtemplate.newfaces);
00108 
00109      if (!err) {
00110           err = quadtemplate.remesh(mesh, a1nodes, oa_nodes, oc_nodes, b3nodes);
00111           if (!err) addNewElements( quadtemplate.newnodes, quadtemplate.newfaces);
00112      }
00113 
00114      if (!err) {
00115           err = quadtemplate.remesh(mesh, a2nodes, ob_nodes, oa_nodes, b1nodes);
00116           if (!err) addNewElements( quadtemplate.newnodes, quadtemplate.newfaces);
00117      }
00118 
00119      if (err) {
00120           discard();
00121           return 1;
00122      }
00123 
00124      // Check for any inversion of the element, if there is inversion,
00125      // undo everthing (i.e. remove new nodes and faces).
00126      //
00127      nSize = newfaces.size();
00128      for (size_t i = 0; i < nSize; i++) {
00129           if (newfaces[i]->concaveAt() >= 0) {
00130                discard();
00131                return 2;
00132           }
00133      }
00134 
00135      return 0;
00136 }
00137 
00139 
00140 int
00141 QuadRemeshTemplate :: remesh()
00142 {
00143      size_t nSize = boundnodes.size();
00144      assert( nSize == size_t(2 * nx + 2 * ny - 4));
00145 
00146      qnodes.resize(nx * ny);
00147 
00148      //
00149      // Put the boundary nodes on the structured mesh: The orientation is
00150      // Counter clockwise ( south->east->north->west );
00151      //
00152      int offset, index = 0;
00153 
00154      // South Side ...
00155      index = 0;
00156      for (int i = 0; i < nx; i++) {
00157           offset = i;
00158           qnodes[i] = boundnodes[index++];
00159      }
00160 
00161      // East Side
00162      for (int j = 1; j < ny; j++) {
00163           offset = j * nx + (nx - 1);
00164           qnodes[offset] = boundnodes[index++];
00165      }
00166 
00167      // North Side
00168      for (int i = nx - 2; i >= 0; i--) {
00169           offset = (ny - 1) * nx + i;
00170           qnodes[offset] = boundnodes[index++];
00171      }
00172 
00173      // West Side
00174      for (int j = ny - 2; j >= 1; j--) {
00175           offset = j*nx;
00176           qnodes[j * nx] = boundnodes[index++];
00177      }
00178 
00179      // New internal nodes must be generated.(either new or from  the trash bin)
00180      mesh->objects_from_pool( (nx-2)*(ny-2), newnodes );
00181 
00182      // Assign the coordinates using Transfinite-interpolation method. For 3D
00183      // Coordinates must be projected on the surface.
00184      //
00185      index = 0;
00186      for (int j = 1; j < ny - 1; j++) {
00187           for (int i = 1; i < nx - 1; i++) {
00188                offset = j * nx + i;
00189                qnodes[offset] = newnodes[index++];
00190                set_tfi_coords(i, j, nx, ny, qnodes); // Coordinates values
00191           }
00192      }
00193 
00194      // All new faces must be generated. ( Objects are either new or from the trash bin)
00195      mesh->objects_from_pool( (nx-1)*(ny-1), newfaces );
00196 
00197      NodeSequence qc(4);
00198 
00199      // Create new faces ...
00200      index = 0;
00201      for (int j = 0; j < ny - 1; j++) {
00202           for (int i = 0; i < nx - 1; i++) {
00203                int offset = j * nx + i;
00204                qc[0] = qnodes[offset];
00205                qc[1] = qnodes[offset + 1];
00206                qc[2] = qnodes[offset + 1 + nx];
00207                qc[3] = qnodes[offset + nx];
00208                Face *face = newfaces[index++];
00209                face->setNodes(qc);
00210                face->setVisitMark(1);
00211           }
00212      }
00213 
00214      // Update the mesh ...
00215      nSize = newnodes.size();
00216      for (size_t i = 0; i < nSize; i++)
00217           mesh->reactivate( newnodes[i] );
00218 
00219      nSize = newfaces.size();
00220      for (size_t i = 0; i < nSize; i++)
00221           mesh->reactivate( newfaces[i] );
00222 
00223      // Check for any inversion of the element, if there is inversion,
00224      // undo everthing (i.e. remove new nodes and faces).
00225      //
00226      nSize = newfaces.size();
00227      for (size_t i = 0; i < nSize; i++) {
00228           if (newfaces[i]->concaveAt() >= 0) {
00229                discard();
00230                return 1;
00231           }
00232      }
00233 
00234      return 0;
00235 }
00236 
00238 
00239 int
00240 QuadRemeshTemplate ::remesh(Mesh *msh,
00241                             NodeSequence &anodes,
00242                             NodeSequence &bnodes,
00243                             NodeSequence &cnodes,
00244                             NodeSequence &dnodes )
00245 {
00246      mesh = msh;
00247 
00248      newnodes.clear();
00249      newfaces.clear();
00250 
00251      assert(anodes.size() == cnodes.size());
00252      assert(bnodes.size() == dnodes.size());
00253 
00254      nx = anodes.size();
00255      ny = bnodes.size();
00256 
00257      boundnodes.resize(2 * nx + 2 * ny - 4);
00258 
00259      int index = 0;
00260 
00261      // First edge is always the reference edge which decides the orientation.
00262      // Append anodes..
00263      size_t nSize = anodes.size();
00264      for (size_t i = 0; i < nSize; i++)
00265           boundnodes[index++] = anodes[i];
00266 
00267      // Append bnodes ...
00268      if (bnodes.front() != anodes.back())
00269           reverse(bnodes.begin(), bnodes.end());
00270 
00271      assert(anodes.back() == bnodes.front());
00272 
00273      nSize = bnodes.size();
00274      for (size_t i = 1; i < nSize; i++)
00275           boundnodes[index++] = bnodes[i];
00276 
00277      // Append cnodes ...
00278      if (cnodes.front() != bnodes.back())
00279           reverse(cnodes.begin(), cnodes.end());
00280 
00281      assert(bnodes.back() == cnodes.front());
00282 
00283      nSize = cnodes.size();
00284      for (size_t i = 1; i < nSize; i++)
00285           boundnodes[index++] = cnodes[i];
00286 
00287      // Append dnodes ...
00288      if (dnodes.front() != cnodes.back())
00289           reverse(dnodes.begin(), dnodes.end());
00290 
00291      assert(cnodes.back() == dnodes.front());
00292      nSize = dnodes.size();
00293      for (size_t i = 1; i < nSize; i++)
00294           boundnodes[index++] = dnodes[i];
00295 
00296      // Ensure that loop is closed ...
00297      assert(anodes.front() == dnodes.back());
00298 
00299      return remesh();
00300 }
00301 
00303 
00304 int
00305 PentaRemeshTemplate :: remesh(Mesh *msh,
00306                               NodeSequence &anodes,
00307                               NodeSequence &bnodes,
00308                               NodeSequence &cnodes,
00309                               NodeSequence &dnodes,
00310                               NodeSequence &enodes,
00311                               int *partition)
00312 {
00313      //
00314      // Implicit in this implementation is the orientation of the newly formed
00315      // faces which is decided by the orientation of the boundary. User need
00316      // not check for the consistency of the mesh. It will be, if the initial
00317      // mesh is consistent.
00318      //
00319      mesh = msh;
00320 
00321      newnodes.clear();
00322      newfaces.clear();
00323 
00324      if (partition == NULL) {
00325           segments[0] = anodes.size() - 1;
00326           segments[1] = bnodes.size() - 1;
00327           segments[2] = cnodes.size() - 1;
00328           segments[3] = dnodes.size() - 1;
00329           segments[4] = enodes.size() - 1;
00330           if (!Face::is_5_sided_convex_loop_quad_meshable(segments, partSegments)) {
00331                cout << "Warning: Penta patch not quad meshable " << endl;
00332                return 1;
00333           }
00334      } else {
00335           for (int i = 0; i < 10; i++)
00336                partSegments[i] = partition[i];
00337      }
00338 
00339      int err;
00340      // Ensure that input is closed loop of three sides ...
00341      assert(anodes.front() == enodes.back());
00342      assert(bnodes.front() == anodes.back());
00343      assert(cnodes.front() == bnodes.back());
00344      assert(dnodes.front() == cnodes.back());
00345      assert(enodes.front() == dnodes.back());
00346 
00347      // Ensure that segments are valid ...
00348      assert((int)anodes.size() == partSegments[0] + partSegments[1] + 1);
00349      assert((int)bnodes.size() == partSegments[2] + partSegments[3] + 1);
00350      assert((int)cnodes.size() == partSegments[4] + partSegments[5] + 1);
00351      assert((int)dnodes.size() == partSegments[6] + partSegments[7] + 1);
00352      assert((int)enodes.size() == partSegments[8] + partSegments[9] + 1);
00353 
00354      // Split Side "A" nodes into a1nodes and b1nodes.
00355      int a1 = partSegments[0];
00356      err = split_stl_vector(anodes, a1 + 1, a1nodes, b1nodes);
00357      if (err) return 1;
00358 
00359      // Split Side "B" nodes into a2nodes and b2nodes.
00360      int a2 = partSegments[2];
00361      err = split_stl_vector(bnodes, a2 + 1, a2nodes, b2nodes);
00362      if (err) return 1;
00363 
00364      // Split Side "C" nodes into a3nodes and b3nodes.
00365      int a3 = partSegments[4];
00366      err = split_stl_vector(cnodes, a3 + 1, a3nodes, b3nodes);
00367      if (err) return 1;
00368 
00369      // Split Side "C" nodes into a3nodes and b3nodes.
00370      int a4 = partSegments[6];
00371      split_stl_vector(dnodes, a4 + 1, a4nodes, b4nodes);
00372      if (err) return 1;
00373 
00374      // Split Side "C" nodes into a3nodes and b3nodes.
00375      int a5 = partSegments[8];
00376      err = split_stl_vector(enodes, a5 + 1, a5nodes, b5nodes);
00377      if (err) return 1;
00378 
00379      // Splitting nodes on each side
00380      Vertex *ca = a1nodes.back();
00381      Vertex *cb = a2nodes.back();
00382      Vertex *cc = a3nodes.back();
00383      Vertex *cd = a4nodes.back();
00384      Vertex *ce = a5nodes.back();
00385 
00386      // Centoid may not be within the patch. But take the risk.
00387      Vertex *co = Face::centroid(ca, cb, cc, cd, ce);
00388 
00389      mesh->addNode(co);
00390      newnodes.push_back(co);
00391 
00392      size_t nSize;
00393      linear_interpolation(mesh, co, ca, a2 + 1, oa_nodes);
00394      nSize = oa_nodes.size();
00395      for (size_t i = 1; i < nSize - 1; i++) {
00396           newnodes.push_back(oa_nodes[i]);
00397      }
00398 
00399      linear_interpolation(mesh, co, cb, a3 + 1, ob_nodes);
00400      nSize = ob_nodes.size();
00401      for (size_t i = 1; i < nSize - 1; i++) {
00402           newnodes.push_back(ob_nodes[i]);
00403      }
00404 
00405      linear_interpolation(mesh, co, cc, a4 + 1, oc_nodes);
00406      nSize = oc_nodes.size();
00407      for (size_t i = 1; i < nSize - 1; i++) {
00408           newnodes.push_back(oc_nodes[i]);
00409      }
00410 
00411      linear_interpolation(mesh, co, cd, a5 + 1, od_nodes);
00412      nSize = od_nodes.size();
00413      for (size_t i = 1; i < nSize - 1; i++) {
00414           newnodes.push_back(od_nodes[i]);
00415      }
00416 
00417      linear_interpolation(mesh, co, ce, a1 + 1, oe_nodes);
00418      nSize = oe_nodes.size();
00419      for (size_t i = 1; i < nSize - 1; i++) {
00420           newnodes.push_back(oe_nodes[i]);
00421      }
00422 
00423      err = quadtemplate.remesh(mesh, a1nodes, oa_nodes, oe_nodes, b5nodes);
00424      if (!err) addNewElements( quadtemplate.newnodes, quadtemplate.newfaces);
00425 
00426      if (!err) {
00427           err = quadtemplate.remesh(mesh, a2nodes, ob_nodes, oa_nodes, b1nodes);
00428           if (!err) addNewElements( quadtemplate.newnodes, quadtemplate.newfaces);
00429      }
00430 
00431      if (!err) {
00432           err = quadtemplate.remesh(mesh, a3nodes, oc_nodes, ob_nodes, b2nodes);
00433           if (!err) addNewElements( quadtemplate.newnodes, quadtemplate.newfaces);
00434      }
00435 
00436      if (!err) {
00437           err = quadtemplate.remesh(mesh, a4nodes, od_nodes, oc_nodes, b3nodes );
00438           if (!err) addNewElements( quadtemplate.newnodes, quadtemplate.newfaces);
00439      }
00440 
00441      if (!err) {
00442           err = quadtemplate.remesh(mesh, a5nodes, oe_nodes, od_nodes, b4nodes);
00443           if (!err) addNewElements( quadtemplate.newnodes, quadtemplate.newfaces);
00444      }
00445 
00446      if (err) {
00447           discard();
00448           return 1;
00449      }
00450 
00451      // Check for any inversion of the element, if there is inversion,
00452      // undo everthing (i.e. remove new nodes and faces).
00453      //
00454      nSize = newfaces.size();
00455      for (size_t i = 0; i < nSize; i++) {
00456           if (newfaces[i]->concaveAt() >= 0) {
00457                discard();
00458                return 2;
00459           }
00460      }
00461 
00462      return 0;
00463 }
00464 
00466 
00467 void OneDefectPatch::setAttributes()
00468 {
00469      size_t nSize;
00470 
00471      nSize = mesh->getSize(0);
00472      for (size_t i = 0; i < nSize; i++) {
00473           Vertex *v = mesh->getNodeAt(i);
00474           v->setAttribute("DefectPatch", 0);
00475      }
00476 
00477      nSize = mesh->getSize(2);
00478      for (size_t i = 0; i < nSize; i++) {
00479           Face *f = mesh->getFaceAt(i);
00480           f->setAttribute("DefectPatch", 0);
00481      }
00482 
00483      nSize = nodepath.size();
00484      for (size_t i = 0; i < nSize; i++)
00485           nodepath[i]->setAttribute("DefectPatch", 1);
00486 
00487      nSize = bound_nodes.size();
00488      for (size_t i = 0; i <  nSize; i++)
00489           bound_nodes[i]->setAttribute("DefectPatch", 3);
00490 
00491      NodeSet::const_iterator vit;
00492      for (vit = corners.begin(); vit != corners.end(); ++vit)
00493           (*vit)->setAttribute("DefectPatch", 2);
00494 
00495      FaceSet::const_iterator it;
00496      for (it = faces.begin(); it != faces.end(); ++it)
00497           (*it)->setAttribute("DefectPatch", 1);
00498 
00499 }
00500 
00502 
00503 bool OneDefectPatch::isSafe()
00504 {
00505      FaceSet::const_iterator it;
00506      for (it = faces.begin(); it != faces.end(); ++it)
00507           if((*it)->isRemoved()) return 0;
00508 
00509      // Check if the previous updates modified the splitting node degree.
00510      if (quad_splitting_node) {
00511           int nSize = quad_splitting_node->getNumRelations(2);
00512           if ( nSize != quad_splitting_node_degree) return 0;
00513      }
00514 
00515      return 1;
00516 }
00518 
00519 void OneDefectPatch::start_boundary_loop_from(Vertex *vmid)
00520 {
00521      assert(corners.find(vmid) != corners.end());
00522 
00523      NodeSequence::iterator middle;
00524      middle = find(bound_nodes.begin(), bound_nodes.end(), vmid);
00525      assert(middle != bound_nodes.end());
00526 
00527      std::rotate(bound_nodes.begin(), middle, bound_nodes.end());
00528      assert(bound_nodes[0] == vmid);
00529 
00530      set_boundary_segments();
00531 }
00533 
00534 void OneDefectPatch::get_bound_nodes(const Vertex *src, const Vertex *dst, NodeSequence &seq)
00535 {
00536      int start_pos = getPosOf(src);
00537      int end_pos = getPosOf(dst);
00538      int nsize = bound_nodes.size();
00539 
00540      assert(nsize > 1);
00541 
00542      if (end_pos == 0) end_pos = nsize;
00543      assert(end_pos > start_pos);
00544 
00545      seq.resize(end_pos - start_pos + 1);
00546      int index = 0;
00547      for (int i = start_pos; i <= end_pos; i++)
00548           seq[index++] = bound_nodes[i % nsize];
00549 }
00550 
00552 
00553 bool OneDefectPatch::has_irregular_node_on_first_segment() const
00554 {
00555      int start_pos = cornerPos[0];
00556      int end_pos = start_pos + segSize[0] - 1;
00557 
00558      for (int i = 1; i < end_pos - 1; i++)
00559           if( bound_nodes[i]->getNumRelations(2) < 4 ) return 1;
00560 
00561      return 0;
00562 }
00563 
00565 
00566 bool OneDefectPatch::is_simply_connected()
00567 {
00568      size_t V = inner_nodes.size() + bound_nodes.size();
00569      size_t F = faces.size();
00570 
00571      vector<Edge> edges;
00572      edges.reserve(4 * F);
00573 
00574      FaceSet::const_iterator it;
00575      for (it = faces.begin(); it != faces.end(); ++it) {
00576           Face *face = *it;
00577           for (int i = 0; i < 4; i++) {
00578                Vertex *v0 = face->getNodeAt(i);
00579                Vertex *v1 = face->getNodeAt(i + 1);
00580                Edge edge(v0, v1);
00581                int found = 0;
00582                size_t nSize =  edges.size();
00583                for (size_t j = 0; j < nSize; j++) {
00584                     if (edges[j].isSameAs(edge)) {
00585                          found = 1;
00586                          break;
00587                     }
00588                }
00589                if (!found) edges.push_back(edge);
00590           }
00591      }
00592 
00593      size_t E = edges.size();
00594 
00595      if (V - E + F == 1) return 1;
00596      return 0;
00597 }
00599 
00600 int OneDefectPatch::reorient_4_sided_loop()
00601 {
00602      //Always remeshable, nothing has to be done...
00603      if ((segSize[0] == segSize[2]) && (segSize[1] == segSize[3])) return 0;
00604 
00606      // Defination:  A four sided convex loop has four segments.
00607      // Objectives:  A four sided convex loop must be orietned such that
00608      //   1.  First segment must be smaller than the third one, because
00609      //       we need to create a triangle patch based at segment#3.
00610      //
00611      //   2.  If there are two choices, then the side having irregular
00612      //       node must be given higher priority. ( Here irregulaty means
00613      //       that vertex valency < 4 ).
00614      //
00615      //   3.  A side having less number of nodes on the first segment than
00616      //       the third is given preference.
00617      //
00618      // Pre-Conditions  :  A loop must be oriented ( CW or CCW ).
00619      //
00620      // Date: 17th Nov. 2010.
00622 
00623      Vertex *start_corner = NULL;
00624 
00625      if (segSize[0] == segSize[2]) {
00626           if (min(segSize[1], segSize[3]) == 2) return 1;
00627           //  Either Segment 2 or 3 must be starting node
00628           if (segSize[1] < segSize[3])
00629                start_corner = bound_nodes[ cornerPos[1] ];
00630           else
00631                start_corner = bound_nodes[ cornerPos[3] ];
00632           start_boundary_loop_from(start_corner);
00633      }
00634 
00635      if (min(segSize[0], segSize[2]) == 2) return 1;
00636 
00637      if (segSize[2] < segSize[0]) {
00638           start_corner = bound_nodes[ cornerPos[2] ];
00639           start_boundary_loop_from(start_corner);
00640      }
00641 
00642      // By this stage, the loop must be reoriented correctly.
00643      assert(segSize[0] < segSize[2]);
00644 
00645      // Great, we found one irregular node on the first boundary...
00646      if (has_irregular_node_on_first_segment()) return 1;
00647 
00648      // If the segment 2 and 3 have same size, Alas, nothing can be done.
00649      if (segSize[1] == segSize[3]) return 1;
00650 
00651      if (min(segSize[1], segSize[3]) == 2) return 1;
00652 
00653      if (segSize[3] < segSize[1]) {
00654           start_corner = bound_nodes[ cornerPos[3] ];
00655           start_boundary_loop_from(start_corner);
00656      } else {
00657           start_corner = bound_nodes[ cornerPos[1] ];
00658           start_boundary_loop_from(start_corner);
00659      }
00660 
00661      //
00662      // Note that we didn't check for irregular node here. So if this segment
00663      // has at least one irregular node, then we are lucky. Otherwise decision
00664      // to remesh it done based wthere remeshing will result in the reduction
00665      // of irregular nodes in patch.
00666      //
00667      assert(segSize[0] < segSize[2]);
00668      return 0;
00669 }
00670 
00672 
00673 int OneDefectPatch::remesh_3_sided_patch()
00674 {
00675      Vertex *c0, *c1, *c2;
00676 
00677      c0 = bound_nodes[ cornerPos[0] ];
00678      c1 = bound_nodes[ cornerPos[1] ];
00679      c2 = bound_nodes[ cornerPos[2] ];
00680 
00681      get_bound_nodes(c0, c1, anodes);
00682      get_bound_nodes(c1, c2, bnodes);
00683      get_bound_nodes(c2, c0, cnodes);
00684 
00685      int err = template3.remesh(mesh, anodes, bnodes, cnodes, partSegments);
00686      if( !err ) {
00687           num_3_patches++;
00688           for( size_t i = 0; i < template3.newnodes.size(); i++)
00689                newnodes.push_back( template3.newnodes[i] );
00690           for( size_t i = 0; i < template3.newfaces.size(); i++)
00691                newfaces.push_back( template3.newfaces[i] );
00692      }
00693      return err;
00694 }
00696 
00697 int OneDefectPatch::remesh_4_sided_patch()
00698 {
00699      int err;
00700      size_t nSize;
00701      Vertex *c0, *c1, *c2, *c3, *ta, *tb, *tc;
00702 
00703      newnodes.clear();
00704      newfaces.clear();
00705 
00706      // Collect landmark nodes on the patch..
00707      c0 = bound_nodes[ cornerPos[0] ];
00708      c1 = bound_nodes[ cornerPos[1] ];
00709      c2 = bound_nodes[ cornerPos[2] ];
00710      c3 = bound_nodes[ cornerPos[3] ];
00711      ta = quad_splitting_node;
00712 
00713      if (ta == NULL) {
00714           if ((segSize[0] == segSize[2]) && (segSize[1] == segSize[3])) {
00715                get_bound_nodes(c0, c1, anodes);
00716                get_bound_nodes(c1, c2, bnodes);
00717                get_bound_nodes(c2, c3, cnodes);
00718                get_bound_nodes(c3, c0, dnodes);
00719                err = template4.remesh(mesh, anodes, bnodes, cnodes, dnodes);
00720                if (!err) {
00721                    for( size_t i = 0; i < template4.newnodes.size(); i++)
00722                     newnodes.push_back( template4.newnodes[i] );
00723                    for( size_t i = 0; i < template4.newfaces.size(); i++)
00724                     newfaces.push_back( template4.newfaces[i] );
00725                 }
00726                return err;
00727           }
00728      }
00729 
00730      assert(ta != NULL);
00731 
00732      get_bound_nodes(c0, ta, a1nodes);
00733      get_bound_nodes(ta, c1, a2nodes);
00734 
00735      tb = bound_nodes[cornerPos[2] + a2nodes.size() - 1 ];
00736      tc = bound_nodes[cornerPos[3] - a1nodes.size() + 1];
00737      get_bound_nodes(tb, tc, c0nodes);
00738 
00739      linear_interpolation(mesh, ta, tb, segSize[1], abnodes);
00740 
00741      nSize = abnodes.size();
00742      for (size_t  i = 1; i < nSize - 1; i++)
00743           newnodes.push_back(abnodes[i]);
00744 
00745      linear_interpolation(mesh, tc, ta, segSize[3], canodes);
00746 
00747      nSize = canodes.size();
00748      for (size_t i = 1; i < nSize - 1; i++)
00749           newnodes.push_back(canodes[i]);
00750 
00751      get_bound_nodes(tb, tc, bcnodes);
00752      get_bound_nodes(c1, c2, b1nodes);
00753      get_bound_nodes(c2, tb, c2nodes);
00754      get_bound_nodes(tc, c3, c1nodes);
00755      get_bound_nodes(c3, c0, d1nodes);
00756 
00757      err = template3.remesh(mesh, abnodes, bcnodes, canodes, partSegments);
00758      if (!err) {
00759           for( size_t i = 0; i < template3.newnodes.size(); i++)
00760                newnodes.push_back( template3.newnodes[i] );
00761           for( size_t i = 0; i < template3.newfaces.size(); i++)
00762                newfaces.push_back( template3.newfaces[i] );
00763      }
00764 
00765      if (!err) {
00766           err = template4.remesh(mesh, a1nodes, canodes, c1nodes, d1nodes);
00767           if (!err) {
00768                for( size_t i = 0; i < template4.newnodes.size(); i++)
00769                     newnodes.push_back( template4.newnodes[i] );
00770                for( size_t i = 0; i < template4.newfaces.size(); i++)
00771                     newfaces.push_back( template4.newfaces[i] );
00772           }
00773      }
00774 
00775      if (!err) {
00776           err = template4.remesh(mesh, a2nodes, b1nodes, c2nodes, abnodes);
00777           if (!err) {
00778                for( size_t i = 0; i < template4.newnodes.size(); i++)
00779                     newnodes.push_back( template4.newnodes[i] );
00780                for( size_t i = 0; i < template4.newfaces.size(); i++)
00781                     newfaces.push_back( template4.newfaces[i] );
00782           }
00783      }
00784 
00785      if (err) {
00786           nSize = newfaces.size();
00787           for (size_t i = 0; i < nSize; i++)
00788                mesh->remove(newfaces[i]);
00789 
00790           nSize = newnodes.size();
00791           for (size_t i = 0; i < nSize; i++)
00792                mesh->remove(newnodes[i]);
00793 
00794           newnodes.clear();
00795           newfaces.clear();
00796           return 2;
00797      }
00798 
00799      if( !err ) num_4_patches++;
00800 
00801      return err;
00802 }
00803 
00805 
00806 int OneDefectPatch::remesh_5_sided_patch()
00807 {
00808      Vertex *c0, *c1, *c2, *c3, *c4;
00809      c0 = bound_nodes[ cornerPos[0] ];
00810      c1 = bound_nodes[ cornerPos[1] ];
00811      c2 = bound_nodes[ cornerPos[2] ];
00812      c3 = bound_nodes[ cornerPos[3] ];
00813      c4 = bound_nodes[ cornerPos[4] ];
00814 
00815      get_bound_nodes(c0, c1, anodes);
00816      get_bound_nodes(c1, c2, bnodes);
00817      get_bound_nodes(c2, c3, cnodes);
00818      get_bound_nodes(c3, c4, dnodes);
00819      get_bound_nodes(c4, c0, enodes);
00820 
00821      int err = template5.remesh(mesh, anodes, bnodes, cnodes, dnodes, enodes, partSegments);
00822 
00823      if( !err)  {
00824           num_5_patches++;
00825           for( size_t i = 0; i < template5.newnodes.size(); i++)
00826                newnodes.push_back( template5.newnodes[i] );
00827           for( size_t i = 0; i < template5.newfaces.size(); i++)
00828                newfaces.push_back( template5.newfaces[i] );
00829      }
00830 
00831      return err;
00832 }
00833 
00835 
00836 bool OneDefectPatch::is_quad_breakable_at(const Vertex *ta)
00837 {
00838      // Landmark nodes on the boundary. There are seven such nodes.
00839      Vertex *c0, *c1, *tb, *tc;
00840      //     Vertex  *c2, *c3;
00841      c0 = bound_nodes[ cornerPos[0] ];
00842      c1 = bound_nodes[ cornerPos[1] ];
00843      //     c2 = bound_nodes[ cornerPos[2] ];
00844      //     c3 = bound_nodes[ cornerPos[3] ];
00845 
00846      get_bound_nodes(c0, ta, a1nodes);
00847      get_bound_nodes(ta, c1, a2nodes);
00848 
00849      tb = bound_nodes[cornerPos[2] + a2nodes.size() - 1 ];
00850      tc = bound_nodes[cornerPos[3] - a1nodes.size() + 1];
00851      get_bound_nodes(tb, tc, c0nodes);
00852 
00853      int segments[3];
00854      segments[0] = c0nodes.size() - 1;
00855      segments[1] = segSize[3] - 1;
00856      segments[2] = segSize[2] - 1;
00857 
00858      if (Face::is_3_sided_convex_loop_quad_meshable(segments, partSegments))
00859           return 0;
00860 
00861      return 1;
00862 }
00864 
00865 bool OneDefectPatch::is_4_sided_convex_loop_quad_meshable()
00866 {
00867      assert(corners.size() == 4);
00868 
00869      // It is always quadrangulable.
00870      if ((segSize[0] == segSize[2]) && (segSize[1] == segSize[3])) return 1;
00871 
00872      int nr0 = count_irregular_nodes(0);
00873      int nr1 = count_irregular_nodes(1);
00874 
00875      //
00876      // If there are only two irregular nodes in the patch ( inner + boundary)
00877      // then it may be quadrangable, but it will not reduce number of irregular
00878      // nodes ( one irregular node and one irregular node on the boundary
00879      // will be created. Since there is no reduction in #of irregular nodes,
00880      // it is not profitable, and we just return it.
00881      //
00882      // If there are more than two irrgular nodes, then we can ensure that
00883      // that there will at most two irregular nodes in the patch ( one inside
00884      // and probably one on the boundary.
00885      //
00886 
00887      if (nr0 + nr1 == 2) return 0;
00888 
00889      // The loop is reoriented so that the first edge is smaller than the
00890      // opposite edge ....
00891      if (reorient_4_sided_loop()) return 0;
00892 
00893      // We need at least 3 points on the first seqment ...
00894      if (segSize[0] < 3) return 0;
00895 
00896      // Landmark nodes on the boundary. There are seven such nodes.
00897      Vertex *c0, *c1, *ta, *tb, *tc;
00898      //     Vertex *c2, *c3;
00899      c0 = bound_nodes[ cornerPos[0] ];
00900      c1 = bound_nodes[ cornerPos[1] ];
00901      //     c2 = bound_nodes[ cornerPos[2] ];
00902      //     c3 = bound_nodes[ cornerPos[3] ];
00903 
00904      // First let us prioritieze vertices on the first segment that will
00905      // be tried for splitting the region. First of all, if #of irregular
00906      // nodes inside the region is 2, then we must have one irregular
00907      // node on the boundary. But even if the #irregular nodes inside
00908      // the region > 2, we still prefer irregular node in the boundary
00909      // as it will reduce irregularity by one.
00910 
00911      // Both end nodes on the segment are not included in splitting.
00912 
00913      // Method: Start node(1,...n-1) and check it could lead to acceptable
00914      // splitting. A proper splitting means that the triangle sandwitched
00915      // between two quads regions is "quad-meshable" with one defect.
00916      //
00917      // Future Work: It is possible that the triangle region is not "one
00918      // defect" meshable, but with more than one defect. We will extend
00919      // this case later.
00920      //
00921 
00922      NodeSequence trynodes;
00923      for (int i = 1; i < segSize[0] - 1; i++) {
00924           if (bound_nodes[i]->getNumRelations(2) == 3)
00925                trynodes.push_back(bound_nodes[i]);
00926      }
00927 
00928      if (trynodes.empty() && nr0 == 2) {
00929           // Patch not remeshable because there are no irregular node on the boundary
00930           return 0;
00931      }
00932 
00933      if (nr0 > 2) {
00934           for (int i = 1; i < segSize[0] - 1; i++) {
00935                if (bound_nodes[i]->getNumRelations(2) == 4)
00936                     trynodes.push_back(bound_nodes[i]);
00937           }
00938      }
00939 
00940      quad_splitting_node = NULL;
00941 
00942      int segments[3];
00943      int nSize = trynodes.size();
00944      for (int i = 0; i < nSize; i++) {
00945           ta = trynodes[i];
00946           get_bound_nodes(c0, ta, a1nodes);
00947           get_bound_nodes(ta, c1, a2nodes);
00948 
00949           tb = bound_nodes[cornerPos[2] + a2nodes.size() - 1 ];
00950           tc = bound_nodes[cornerPos[3] - a1nodes.size() + 1];
00951 
00952           get_bound_nodes(tb, tc, bcnodes);
00953           segments[0] = segSize[1] - 1;
00954           segments[1] = bcnodes.size() - 1;
00955           segments[2] = segSize[3] - 1;
00956 
00957           if (Face::is_3_sided_convex_loop_quad_meshable(segments, partSegments)) {
00958                quad_splitting_node = ta;
00959                quad_splitting_node_degree = ta->getNumRelations(2);
00960                return 1;
00961           }
00962      }
00963 
00964      return 0;
00965 }
00966 
00968 
00969 int OneDefectPatch::init_blob()
00970 {
00971      nodes.clear();
00972 
00973      inner_nodes.clear();
00974      bound_nodes.clear();
00975 
00976      faces.clear();
00977      inner_faces.clear();
00978 
00979      relations02.clear();
00980 
00981      int err;
00982      size_t nSize;
00983      if (apex) {
00984           err = expand_blob( apex );
00985           if( err ) return 1;
00986      } else {
00987           nSize = nodepath.size();
00988           assert( nSize > 0);
00989           for (size_t i = 0; i < nSize; i++) {
00990                Vertex *v = nodepath[i];
00991                err = expand_blob(v );
00992                if( err ) return 1;
00993           }
00994      }
00995      return 0;
00996 }
00997 
00999 
01000 size_t OneDefectPatch::count_irregular_nodes(int where)
01001 {
01002      NodeSequence::const_iterator it;
01003      assert(mesh->getAdjTable(0, 2));
01004 
01005      size_t ncount = 0;
01006 
01007      if (where == 0) {
01008           for (it = inner_nodes.begin(); it != inner_nodes.end(); ++it) {
01009                Vertex *v = *it;
01010                if( v->isActive() && (v->getNumRelations(2) != 4)) ncount++;
01011           }
01012      } else {
01013           for (size_t i = 0; i < bound_nodes.size(); i++) {
01014                Vertex *v = bound_nodes[i];
01015                if( v->isActive() )
01016                     if (!v->isBoundary() && v->getNumRelations(2) == 3) ncount++;
01017           }
01018      }
01019 
01020      return ncount;
01021 }
01023 int
01024 OneDefectPatch::update_boundary()
01025 {
01026      StopWatch wc;
01027 
01028      boundary.clear();
01029      bound_nodes.clear();
01030      inner_nodes.clear();
01031 
01032      wc.start();
01033 
01034      // We need to rebuild relations locally to identfy corners and boundary.
01035      FaceSet::const_iterator fiter;
01036 
01037      assert( faces.size() > 0);
01038 
01039      // A boundary edge must have exactly one face neighbor...
01040      FaceSequence faceneighs;
01041      for (fiter = faces.begin(); fiter != faces.end(); ++fiter) {
01042           Face *face = *fiter;
01043           if( inner_faces.find(face) == inner_faces.end() ) {
01044                int ncount = 0;
01045                for (int j = 0; j < 4; j++) {
01046                     Vertex *v0 = face->getNodeAt(j + 0);
01047                     Vertex *v1 = face->getNodeAt(j + 1);
01048                     faceneighs.clear();
01049                     assert(relations02[v0].size() > 0);
01050                     assert(relations02[v1].size() > 0);
01051                     set_intersection(relations02[v0].begin(), relations02[v0].end(),
01052                                      relations02[v1].begin(), relations02[v1].end(),
01053                                      back_inserter(faceneighs));
01054                     if (faceneighs.size() == 1) {
01055                          ncount = 1;
01056                          Edge newedge(v0, v1);
01057                          bound_nodes.push_back(v0);
01058                          bound_nodes.push_back(v1);
01059                          boundary.push_back(newedge);
01060                     }
01061                }
01062                if( ncount == 0) inner_faces.insert(face);
01063           }
01064      }
01065      int nSize = bound_nodes.size()/2;
01066      assert( nSize > 0);
01067 
01068      sort( bound_nodes.begin(), bound_nodes.end() );
01069 
01070      int simple= 1;
01071      for( int i = 0; i < nSize; i++)  {
01072           if( bound_nodes[2*i] != bound_nodes[2*i+1] )  {
01073                simple= 0;
01074                break;
01075           }
01076      }
01077 
01078      if( simple ) {
01079           sort( bound_nodes.begin(), bound_nodes.end() );
01080           NodeSequence::iterator it;
01081           it  = unique( bound_nodes.begin(), bound_nodes.end() );
01082           bound_nodes.erase( it, bound_nodes.end() );
01083           set_difference(nodes.begin(), nodes.end(),
01084                          bound_nodes.begin(), bound_nodes.end(),
01085                          back_inserter(inner_nodes));
01086      }
01087 
01088      wc.stop();
01089      exec_time += wc.getSeconds();
01090 
01091      if( !simple) return 1;
01092 
01093      return 0;
01094 }
01095 
01097 
01098 int
01099 OneDefectPatch::finalize_boundary()
01100 {
01101      corners.clear();
01102 
01103      // Sequence the chain and start from one of the corner...
01104      int err = Mesh::make_chain(boundary);
01105      if (err) return 2;
01106 
01107      size_t nSize = boundary.size();
01108 
01109      Vertex *vertex;
01110      for (size_t k = 0; k < nSize; k++) {
01111           vertex = boundary[k].getNodeAt(0);
01112           if (relations02[vertex].size() == 1) corners.insert(vertex);
01113 
01114           vertex = boundary[k].getNodeAt(1);
01115           if (relations02[vertex].size() == 1) corners.insert(vertex);
01116      }
01117 
01118      if( corners.size() < 3) return 3;
01119 
01120      // Start the chain from one of the corners.
01121      err = Mesh::rotate_chain(boundary, *corners.begin());
01122      if (err) return 4;
01123 
01124      bound_nodes.resize( nSize );
01125      for (size_t k = 0; k < nSize; k++)
01126           bound_nodes[k] = boundary[k].getNodeAt(0);
01127 
01128      // Split the boundary loop into segments.
01129      // (i.e. End of the segments are the corners identified earlier )
01130      set_boundary_segments();
01131 
01132      return 0;
01133 }
01134 
01135 
01137 
01138 void OneDefectPatch::set_boundary_segments()
01139 {
01140      size_t nSize = corners.size();
01141      // Although this stage will not come in this algorithm...
01142      if ( nSize == 0) return;
01143 
01144      cornerPos.resize( nSize + 1);
01145 
01146      NodeSet::const_iterator it;
01147      int index = 0;
01148      for (it = corners.begin(); it != corners.end(); ++it) {
01149           cornerPos[index++] = getPosOf(*it);
01150      }
01151 
01152      cornerPos[corners.size()] = bound_nodes.size();
01153 
01154      sort(cornerPos.begin(), cornerPos.end());
01155 
01156      segSize.resize( nSize );
01157 
01158      for (size_t i = 0; i < nSize; i++)
01159           segSize[i] = cornerPos[(i + 1)] - cornerPos[i] + 1;
01160 }
01161 
01163 
01164 int OneDefectPatch::get_topological_outer_angle(Vertex *vertex)
01165 {
01166      FaceSequence vfaces;
01167      vertex->getRelations( vfaces );
01168 
01169      // How many faces are outside the regions.
01170      int ncount = 0;
01171      int nSize  = vfaces.size();
01172      for (int i = 0; i < nSize; i++)
01173           if (faces.find(vfaces[i]) == faces.end()) ncount++;
01174 
01175      return ncount - 2;
01176 }
01178 
01179 int OneDefectPatch::expand_blob(Vertex *vertex)
01180 {
01181      FaceSequence vfaces;
01182      vertex->getRelations( vfaces );
01183 
01184      int nSize = vfaces.size();
01185      assert( nSize > 0);
01186      for (int i = 0; i < nSize; i++) {
01187           Face *face = vfaces[i];
01188           assert(face->isActive() );
01189           if( face->isVisited() ) return 1;
01190           faces.insert(face);
01191           for (int j = 0; j < 4; j++) {
01192                Vertex *vf = face->getNodeAt(j);
01193                relations02[vf].insert(face);
01194                nodes.insert( vf );
01195           }
01196      }
01197      return 0;
01198 }
01199 
01201 
01202 int OneDefectPatch::expand_blob()
01203 {
01204      int err;
01205      size_t before_expansion = faces.size();
01206 
01207      size_t nSize = bound_nodes.size();
01208      assert( nSize > 0);
01209 
01210      for (size_t i = 0; i < nSize; i++) {
01211           err = expand_blob(bound_nodes[i]);
01212           if( err ) return 1;
01213      }
01214 
01215      size_t after_expansion = faces.size();
01216 
01217      // bad luck; no expansion took place....
01218      if( after_expansion == before_expansion ) return 1;
01219 
01220      return 0;
01221 
01222 }
01223 
01224 int
01225 OneDefectPatch::build_remeshable_boundary()
01226 {
01227      assert(mesh->getAdjTable(0, 2));
01228 
01230      // A valid boundary patch has the following propeties:
01231      //
01232      // 1) The loop is closed.
01233      // 2) The region is simply connected. The Euler characteristic must be
01234      //    one, which mean there are no holes in the region.
01235      // 3) The loop has 3, 4 or 5 corner nodes.
01236      // 4) The loop form almost convex boundary.
01237      //
01238      // If these conditions are met, we may try to remesh the region with
01239      // quadrilateral elements. Sometimes it may not be possible to remesh
01240      // a reason, or the resulting elements have unacceptable quality that
01241      // we would like to avoid.
01242      //
01243      // A resulting mesh MUST decrease the irregular nodes count.
01245      int err, topo_convex_region;
01246 
01247      err = init_blob();
01248      if( err ) return 1;
01249 
01250      size_t  nSize;
01251      while (1) {
01252           err = update_boundary();
01253           if (err) return 2;
01254 
01255           if (faces.size() > MAX_FACES_ALLOWED ) return 3;
01256 
01257           // Expand the blob till topological set is found ..
01258           topo_convex_region = 1;
01259           nSize = bound_nodes.size();
01260           assert( nSize > 0);
01261           for (size_t i = 0; i < nSize; i++) {
01262                if (!bound_nodes[i]->isBoundary()) {
01263                     int topo_angle = get_topological_outer_angle(bound_nodes[i]);
01264                     if (topo_angle < 0) {
01265                          err = expand_blob(bound_nodes[i]);
01266                          if( err ) return 1;
01267                          topo_convex_region = 0;
01268                     }
01269                }
01270           }
01271 
01272           // A valid patch must contain atleast two irregular nodes ...
01273           if( topo_convex_region ) {
01274                if (count_irregular_nodes(0) < 2) {
01275                     topo_convex_region = 0;
01276                     err = expand_blob();
01277                     if( err ) return 1;
01278                }
01279           }
01280 
01281 
01282           // Now finalize the patch for remeshing...
01283           if (topo_convex_region) {
01284                err = finalize_boundary();
01285                if( err ) {
01286                     topo_convex_region = 0;
01287                     err = expand_blob();
01288                     if( err ) return 1;
01289                }
01290           }
01291 
01292           // check the patch for remeshing..
01293           if (topo_convex_region) {
01294                int nsides = corners.size();
01295 
01296                if (nsides > 5) return 6;
01297 
01298                int segments[5], meshable = 0;
01299                switch (nsides) {
01300                case 3:
01301                     segments[0] = segSize[0] - 1;
01302                     segments[1] = segSize[1] - 1;
01303                     segments[2] = segSize[2] - 1;
01304                     meshable = Face::is_3_sided_convex_loop_quad_meshable(segments, partSegments);
01305                     break;
01306 
01307                case 4:
01308                     meshable = is_4_sided_convex_loop_quad_meshable();
01309                     break;
01310 
01311                case 5:
01312                     segments[0] = segSize[0] - 1;
01313                     segments[1] = segSize[1] - 1;
01314                     segments[2] = segSize[2] - 1;
01315                     segments[3] = segSize[3] - 1;
01316                     segments[4] = segSize[4] - 1;
01317                     meshable = Face::is_5_sided_convex_loop_quad_meshable(segments, partSegments);
01318                     break;
01319                }
01320                // If the region is meshable, it is good, it will be remeshed in the next step...
01321                if (meshable) {
01322                     if (!is_simply_connected()) return 4;
01323                     return 0;
01324                }
01325 
01326                /*
01327                               if ( bound_nodes.size()%2 == 0 && count_irregular_nodes(0) > 4) {
01328                                    cout << "BoundNodes " << bound_nodes.size() << " " << count_irregular_nodes(0) << endl;
01329                                    disk_remeshable++;
01330                               }
01331                */
01332 
01333                err = expand_blob();
01334                if( err ) return 1;
01335           }
01336      }
01337 
01338      // You should never come at this stage...
01339      return 1;
01340 }
01341 
01343 
01344 void OneDefectPatch::rollback()
01345 {
01346      NodeSequence::const_iterator it;
01347      for (it = inner_nodes.begin(); it != inner_nodes.end(); ++it)
01348           mesh->reactivate(*it);
01349 
01350      FaceSet::const_iterator fiter;
01351      for (fiter = faces.begin(); fiter != faces.end(); ++fiter)
01352           mesh->reactivate(*fiter);
01353 
01354      size_t nSize = newfaces.size();
01355      for (size_t i = 0; i < nSize; i++)
01356           mesh->remove(newfaces[i]);
01357      newfaces.clear();
01358 
01359      nSize = newnodes.size();
01360      for (size_t i = 0; i < nSize; i++)
01361           mesh->remove(newnodes[i]);
01362      newnodes.clear();
01363 
01364 }
01365 
01367 
01368 void OneDefectPatch::pre_remesh()
01369 {
01370      //
01371      // Before remeshing deactivate all the internal faces and nodes so that all
01372      // the relations are cleared. They are not removed because if something goes
01373      // wrong, the data structures are restored using rollback operations.
01374      //
01375      FaceSet::const_iterator fiter;
01376      for (fiter = faces.begin(); fiter != faces.end(); ++fiter)
01377           mesh->deactivate(*fiter);
01378 
01379      irregular_nodes_removed.clear();
01380      NodeSequence::const_iterator niter;
01381      for (niter = inner_nodes.begin(); niter != inner_nodes.end(); ++niter) {
01382           if (!QuadCleanUp::isRegular(*niter)) irregular_nodes_removed.push_back(*niter);
01383           mesh->deactivate(*niter);
01384      }
01385 
01386      // At this stage, we will have an empty patch ...
01387 }
01388 
01390 
01391 void OneDefectPatch::post_remesh()
01392 {
01393      size_t nSize;
01394 
01395      assert( faces.size() );
01396 
01397      // Remove the old faces ( will go to trash bin );
01398      FaceSet::const_iterator fiter;
01399      for (fiter = faces.begin(); fiter != faces.end(); ++fiter)
01400           mesh->remove(*fiter);
01401 
01402      assert( inner_nodes.size() );
01403      // Remove the old nodes.. (will go to trash bin );
01404      NodeSequence::const_iterator niter;
01405      for (niter = inner_nodes.begin(); niter != inner_nodes.end(); ++niter)
01406           mesh->remove(*niter);
01407 
01408      // Patch contains new nodes new faces now ...
01409      faces.clear();
01410      nSize = newfaces.size();
01411      assert( nSize );
01412      for (size_t i = 0; i < nSize; i++) {
01413           faces.insert(newfaces[i]);
01414      }
01415      newfaces.clear();
01416 
01417      inner_nodes.clear();
01418 
01419      nSize = newnodes.size();
01420      //assert( nSize );
01421      for (size_t i = 0; i < nSize; i++)
01422           inner_nodes.push_back(newnodes[i]);
01423      newnodes.clear();
01424 
01425      int ncount = 0;
01426      new_defective_node = NULL;
01427      for (niter = inner_nodes.begin(); niter != inner_nodes.end(); ++niter) {
01428           Vertex *v = *niter;
01429           if (!QuadCleanUp::isRegular(v)) {
01430                ncount++;
01431                new_defective_node = v;
01432           }
01433      }
01434 
01435      if( ncount > 2) {
01436           cout << "Fatal Error:There must be at the most one defective node in the patch" << endl;
01437           exit(0);
01438      }
01439 
01440      nSize = bound_nodes.size();
01441      assert( nSize );
01442      for( size_t i = 0; i < nSize; i++) {
01443           if( QuadCleanUp::isDoublet( bound_nodes[i] ) )  {
01444                cout << "Warning: Doublet created at the patch boundary: " << endl;
01445           }
01446      }
01447 
01448 }
01449 
01451 
01452 int OneDefectPatch::remesh()
01453 {
01454      // Check if any node or element is modifed by previous operations.
01455      if (!isSafe()) return 1;
01456 
01457      newnodes.clear();
01458      newfaces.clear();
01459 
01460      if (corners.size() > 5) cout << "Warning: Corners more than 5 " << endl;
01461 
01462 #ifdef DEBUG
01463      int nirregular0 = count_irregular_nodes(0) + count_irregular_nodes(1);
01464 #endif
01465 
01466      // Do some pre-processing before remeshing the patch
01467      pre_remesh();
01468 
01469      int err = 1;
01470      switch (corners.size()) {
01471      case 3:
01472           err = remesh_3_sided_patch();
01473           break;
01474      case 4:
01475           err = remesh_4_sided_patch();
01476           break;
01477      case 5:
01478           err = remesh_5_sided_patch();
01479           break;
01480      }
01481 
01482      if (err) {
01483           rollback();
01484           return 1;
01485      }
01486 
01487      // Do some post-processing of the patch
01488      post_remesh();
01489 
01490 #ifdef DEBUG
01491      int nirregular1 = count_irregular_nodes(0) + count_irregular_nodes(1);
01492 
01493      // Our method must stricly decrease the number of irregular nodes,
01494      // otherwise it will be go into infinite loop.
01495 
01496      assert(nirregular1 < nirregular0);
01497 #endif
01498 
01499      return err;
01500 }
01501 
01503 
01504 OneDefectPatch *QuadCleanUp::build_one_defect_patch(Vertex *vertex)
01505 {
01506 
01507      if( defective_patch == NULL )
01508           defective_patch = new OneDefectPatch(mesh);
01509 
01510      defective_patch->clear();
01511 
01512      if( !vertex->isActive() || isRegular(vertex) ) return NULL;
01513 
01514      if( djkpath == NULL ) {
01515           djkpath = new DijkstraShortestPath(mesh, 1);
01516      }
01517 
01518      MeshFilter *filter = new FirstIrregularNode();
01519      const NodeSequence &nodepath = djkpath->getPath(vertex, filter);
01520      delete filter;
01521 
01522      if (isRegular(nodepath.front())) return NULL;
01523      if (isRegular(nodepath.back()))  return NULL;
01524 
01525      defective_patch->set_initial_path( nodepath );
01526 
01527      int err = defective_patch->build_remeshable_boundary();
01528 
01529      if (err) return NULL;
01530 
01531      return defective_patch;
01532 }
01533 
01535 
01536 int QuadCleanUp::remesh_defective_patches()
01537 {
01538      int relexist0 = mesh->build_relations(0, 0);
01539      int relexist2 = mesh->build_relations(0, 2);
01540 
01541      mesh->search_boundary();
01542 
01543      if( !mesh->is_consistently_oriented() )
01544           mesh->make_consistently_oriented();
01545 
01546      djkpath = new DijkstraShortestPath(mesh, 1);
01547 
01548      assert( djkpath);
01549 
01550      assert(mesh->getAdjTable(0, 2));
01551      OneDefectPatch *patch;
01552 
01553      Jaal::MeshOptimization mopt;
01554 
01555      NodeSequence nextSeq;
01556      mesh->get_irregular_nodes(irregular_nodes, 4);
01557 
01558      // Record execution time ...
01559      StopWatch wc[4];
01560 
01561      vector<double> patchCircularity;
01562      vector<int>    patchIrregularity;
01563      vector<int>    meshIrregularity;
01564 
01565      meshIrregularity.push_back( irregular_nodes.size() );
01566 
01567      int ncount1 = 0;
01568      int nfailures = 0;
01569      size_t  num_curr_irregular_nodes = irregular_nodes.size();
01570 
01571      for( int i = 4; i < 12; i++) {
01572           OneDefectPatch :: MAX_FACES_ALLOWED = ( int)pow(2,i);
01573           cout << "Remeshing with maximum patch size :  " << (int)pow(2,i) << endl;
01574           cout << "#of irregular nodes in the mesh " << irregular_nodes.size() << endl;
01575 
01576           wc[3].start();
01577 
01578           while (1) {
01579 
01580                // for creating disjoint faces ..
01581                size_t numFaces = mesh->getSize(2);
01582                for( size_t i = 0; i < numFaces; i++) {
01583                     Face *f = mesh->getFaceAt(i);
01584                     f->setVisitMark(0);
01585                }
01586 
01587                nextSeq.clear();
01588                size_t ncount2 = 0;
01589                size_t nSize = irregular_nodes.size();
01590 
01591                for( size_t i = 0; i < nSize; i++) {
01592                     Vertex *vertex = irregular_nodes[i];
01593                     if( vertex->isActive() ) {
01594                          wc[0].start();
01595                          patch = build_one_defect_patch(vertex);
01596                          wc[0].stop();
01597                          if (patch) {
01598                               wc[1].start();
01599                               int err = patch->remesh();
01600                               wc[1].stop();
01601                               if (!err) {
01602                                    if( ncount1% 100 == 0) cout << "patch count : " << ncount1 << endl;
01603                                    double pq = patch->get_isoperimetic_quotient();
01604                                    patchCircularity.push_back(pq);
01605 
01606                                    int ir = patch->get_irregular_nodes_removed().size();
01607                                    patchIrregularity.push_back(ir);
01608 
01609                                    num_curr_irregular_nodes -= ir;
01610                                    Vertex *vtx = patch->get_new_defective_node();
01611                                    if (vtx) {
01612                                         nextSeq.push_back(vtx);
01613                                         num_curr_irregular_nodes += 1;
01614                                    }
01615 
01616 /*
01617                                    num_curr_irregular_nodes = mesh->count_irregular_nodes(4);
01618                                    meshIrregularity.push_back( num_curr_irregular_nodes);
01619 */
01620 
01621                                    ncount2++;
01622                                    ncount1++;
01623                               } else {
01624                                    nextSeq.push_back( vertex );  // Because patch was not remeshed
01625                                    nfailures++;
01626                               }
01627                          } else
01628                               nextSeq.push_back( vertex ); // Because invalid patch occured.
01629                     }
01630                }
01631                irregular_nodes.swap(nextSeq);
01632 
01633                if (ncount2) {
01634 //                  wc[2].start();
01635 //                  mopt.shape_optimize(mesh, MeshOptimization::QUASI_NEWTON, 10);
01636 //                  wc[2].stop();
01637                }
01638                if( ncount2 == 0) break;
01639           }
01640 
01641           wc[2].start();
01642 //        mopt.shape_optimize(mesh, MeshOptimization::QUASI_NEWTON, 10);
01643           wc[2].stop();
01644 
01645           wc[3].stop();
01646           cout << "Execution Summary: " << endl;
01647           cout << "    Time for searching patches  : " << wc[0].getSeconds()  << endl;
01648           cout << "    Time for remeshing patches  : " << wc[1].getSeconds()  << endl;
01649           cout << "    Time for shape optimization : " << wc[2].getSeconds()  << endl;
01650           cout << "    Total Execution time        : " << wc[3].getSeconds() << endl;
01651      }
01652 
01653      // Perform some aggresive good optimization in the last...
01654      mopt.shape_optimize(mesh, MeshOptimization::QUASI_NEWTON, 100);
01655 
01656      cout << "#of times boundaries created : " << OneDefectPatch::num_boundaries << endl;
01657      cout << "Exec time boundaries created : " << OneDefectPatch::exec_time      << endl;
01658      cout << "3 Sided Patches : " << OneDefectPatch::num_3_patches << endl;
01659      cout << "4 Sided Patches : " << OneDefectPatch::num_4_patches << endl;
01660      cout << "5 Sided Patches : " << OneDefectPatch::num_5_patches << endl;
01661      cout << "n Sides Patches : " << OneDefectPatch::disk_remeshable << endl;
01662      cout << "#Remesh failure : " << nfailures << endl;
01663 
01664      cout << "Execution Summary: " << endl;
01665      cout << "    Time for searching patches  : " << wc[0].getSeconds()  << endl;
01666      cout << "    Time for remeshing patches  : " << wc[1].getSeconds()  << endl;
01667      cout << "    Time for shape optimization : " << wc[2].getSeconds()  << endl;
01668 
01669      ofstream ofile1( "PatchCircularity.dat", ios::out);
01670      for( size_t i = 0; i < patchCircularity.size(); i++)
01671           ofile1 << i+1 << " " << patchCircularity[i] << endl;
01672      ofile1.close();
01673 
01674      ofstream ofile2( "PatchIrregulaity.dat", ios::out);
01675      for( size_t i = 0; i < patchIrregularity.size(); i++)
01676           ofile2 << i+1 << " " << patchIrregularity[i] << endl;
01677      ofile2.close();
01678 
01679      ofstream ofile3( "MeshIrregulaity.dat", ios::out);
01680      for( size_t i = 0; i < meshIrregularity.size(); i++)
01681           ofile3 << i+1 << " " << meshIrregularity[i] << endl;
01682      ofile3.close();
01683 
01684      if (!relexist0)
01685           mesh->clear_relations(0, 0);
01686 
01687      if (!relexist2)
01688           mesh->clear_relations(0, 2);
01689 
01690      return 0;
01691 }
01692 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines