MeshKit  1.0
Mesh.cpp
Go to the documentation of this file.
00001 #include <iomanip>
00002 
00003 #include "meshkit/Mesh.hpp"
00004 #include "meshkit/basic_math.hpp"
00005 #include <limits>
00006 using namespace std;
00007 using namespace Jaal;
00008 
00010 
00011 size_t Vertex::global_id = 0;
00012 
00014 
00015 void
00016 Vertex::mid_point(const Vertex *v0, const Vertex *v1, Point3D &pmid, double alpha)
00017 {
00018      const Point3D &p0 = v0->getXYZCoords();
00019      const Point3D &p1 = v1->getXYZCoords();
00020 
00021      pmid[0] = (1 - alpha) * p0[0] + alpha * p1[0];
00022      pmid[1] = (1 - alpha) * p0[1] + alpha * p1[1];
00023      pmid[2] = (1 - alpha) * p0[2] + alpha * p1[2];
00024 }
00025 
00027 
00028 Vertex*
00029 Vertex::mid_node(const Vertex *v0, const Vertex *v1, double alpha)
00030 {
00031      Point3D xyz;
00032      Vertex::mid_point(v0, v1, xyz, alpha);
00033 
00034      Vertex *vmid = Vertex::newObject();
00035      vmid->setXYZCoords(xyz);
00036      return vmid;
00037 }
00038 
00040 
00041 double
00042 Vertex::length(const Vertex *v0, const Vertex *v1)
00043 {
00044      const Point3D &p0 = v0->getXYZCoords();
00045      const Point3D &p1 = v1->getXYZCoords();
00046 
00047      double dx = p0[0] - p1[0];
00048      double dy = p0[1] - p1[1];
00049      double dz = p0[2] - p1[2];
00050 
00051      return sqrt(dx * dx + dy * dy + dz * dz);
00052 }
00053 
00055 
00056 double
00057 Vertex::length2(const Vertex *v0, const Vertex *v1)
00058 {
00059      const Point3D &p0 = v0->getXYZCoords();
00060      const Point3D &p1 = v1->getXYZCoords();
00061 
00062      double dx = p0[0] - p1[0];
00063      double dy = p0[1] - p1[1];
00064      double dz = p0[2] - p1[2];
00065 
00066      return dx * dx + dy * dy + dz * dz;
00067 }
00068 
00070 
00071 int
00072 Face::quad_tessalate(const NodeSequence &orgNodes, NodeSequence &rotatedNodes)
00073 {
00074      double A013 = tri_area(orgNodes[0]->getXYZCoords(),
00075                             orgNodes[1]->getXYZCoords(),
00076                             orgNodes[3]->getXYZCoords());
00077 
00078      double A123 = tri_area(orgNodes[1]->getXYZCoords(),
00079                             orgNodes[2]->getXYZCoords(),
00080                             orgNodes[3]->getXYZCoords());
00081 
00082      double A012 = tri_area(orgNodes[0]->getXYZCoords(),
00083                             orgNodes[1]->getXYZCoords(),
00084                             orgNodes[2]->getXYZCoords());
00085 
00086      double A023 = tri_area(orgNodes[0]->getXYZCoords(),
00087                             orgNodes[2]->getXYZCoords(),
00088                             orgNodes[3]->getXYZCoords());
00089 
00090      rotatedNodes.resize(4);
00091      if (fabs(A013) + fabs(A123) < fabs(A012) + fabs(A023)) {
00092           rotatedNodes[0] = orgNodes[1];
00093           rotatedNodes[1] = orgNodes[2];
00094           rotatedNodes[2] = orgNodes[3];
00095           rotatedNodes[3] = orgNodes[0];
00096           return 1;
00097      }
00098      rotatedNodes[0] = orgNodes[0];
00099      rotatedNodes[1] = orgNodes[1];
00100      rotatedNodes[2] = orgNodes[2];
00101      rotatedNodes[3] = orgNodes[3];
00102      return 0;
00103 }
00105 
00106 int
00107 Face::triangulate( vector<Face> &trifaces ) const
00108 {
00109      trifaces.clear();
00110      NodeSequence rotatedNodes(4), tconnect(3);
00111      if (this->getSize(0) == 4) {
00112           quad_tessalate(connect, rotatedNodes);
00113           trifaces.resize(2);
00114           tconnect[0] = rotatedNodes[0];
00115           tconnect[1] = rotatedNodes[1];
00116           tconnect[2] = rotatedNodes[2];
00117           trifaces[0].setNodes(tconnect);
00118 
00119           tconnect[0] = rotatedNodes[0];
00120           tconnect[1] = rotatedNodes[2];
00121           tconnect[2] = rotatedNodes[3];
00122           trifaces[1].setNodes(tconnect);
00123      }
00124      return 0;
00125 }
00127 
00128 int
00129 Face::get_interior_angles( vector<double> &vals ) const
00130 {
00131      vals.clear();
00132 
00133      NodeSequence rotatedNodes;
00134      map<Vertex*, double> mapangles;
00135 
00136      int nsize = connect.size();
00137      for (int i = 0; i < nsize; i++)
00138           mapangles[ connect[i] ] = 0.0;
00139      quad_tessalate(connect, rotatedNodes);
00140 
00141      Array3D tangles;
00142 
00143      const Point3D &p0 = rotatedNodes[0]->getXYZCoords();
00144      const Point3D &p1 = rotatedNodes[1]->getXYZCoords();
00145      const Point3D &p2 = rotatedNodes[2]->getXYZCoords();
00146      const Point3D &p3 = rotatedNodes[3]->getXYZCoords();
00147 
00148      Math::getTriAngles(p0, p1, p2, tangles);
00149      mapangles[ rotatedNodes[0]] += tangles[0];
00150      mapangles[ rotatedNodes[1]] += tangles[1];
00151      mapangles[ rotatedNodes[2]] += tangles[2];
00152 
00153      Math::getTriAngles(p0, p2, p3, tangles);
00154      mapangles[ rotatedNodes[0]] += tangles[0];
00155      mapangles[ rotatedNodes[2]] += tangles[1];
00156      mapangles[ rotatedNodes[3]] += tangles[2];
00157 
00158      vals.resize(nsize);
00159      for (int i = 0; i < nsize; i++)
00160           vals[i] = mapangles[ connect[i] ];
00161      return 0;
00162 }
00164 
00165 PNode
00166 Face::opposite_node(const PFace tri, PNode n1, PNode n2)
00167 {
00168      PNode tn0 = tri->getNodeAt(0);
00169      PNode tn1 = tri->getNodeAt(1);
00170      PNode tn2 = tri->getNodeAt(2);
00171 
00172      if (tn0 == n1 && tn1 == n2)
00173           return tn2;
00174      if (tn0 == n2 && tn1 == n1)
00175           return tn2;
00176 
00177      if (tn1 == n1 && tn2 == n2)
00178           return tn0;
00179      if (tn1 == n2 && tn2 == n1)
00180           return tn0;
00181 
00182      if (tn2 == n1 && tn0 == n2)
00183           return tn1;
00184      if (tn2 == n2 && tn0 == n1)
00185           return tn1;
00186 
00187      cout << " Warning: You should not come here " << endl;
00188      cout << " Face " << tn0 << " " << tn1 << " " << tn2 << endl;
00189      cout << " search for " << n1 << "  " << n2 << endl;
00190      exit(0);
00191 }
00192 
00194 
00195 PNode
00196 Face::opposite_node(const PFace face, PNode n1)
00197 {
00198      int pos = face->getPosOf(n1);
00199      if (pos >= 0) {
00200           int size = face->getSize(0);
00201           if (size == 4)
00202                return face->getNodeAt(pos + 2);
00203      }
00204      return NULL;
00205 }
00206 
00208 
00209 int Face::is_cyclic_quad(const Point3D &p0, const Point3D &p1, const Point3D &p2,
00210                          const Point3D &p3)
00211 {
00212      //
00213      // http://en.wikipedia.org/wiki/Ptolemy%27s_theorem
00214      //
00215 
00216      double d02 = Math::length(p0, p2);
00217      double d13 = Math::length(p1, p3);
00218      double d01 = Math::length(p0, p1);
00219      double d12 = Math::length(p1, p2);
00220      double d23 = Math::length(p2, p3);
00221      double d30 = Math::length(p3, p0);
00222 
00223      double lval = d02*d13;
00224      double rval = d01 * d23 + d12*d30;
00225 
00226      if (fabs(lval - rval) < 1.0E-15) return 1;
00227 
00228      return 0;
00229 }
00230 
00232 
00233 void
00234 Face::opposite_nodes(const PFace quad, PNode n1, PNode n2,
00235                      PNode &n3, PNode &n4)
00236 {
00237      PNode qn0 = quad->getNodeAt(0);
00238      PNode qn1 = quad->getNodeAt(1);
00239      PNode qn2 = quad->getNodeAt(2);
00240      PNode qn3 = quad->getNodeAt(3);
00241 
00242      if ((qn0 == n1 && qn1 == n2) || (qn0 == n2 && qn1 == n1)) {
00243           n3 = qn2;
00244           n4 = qn3;
00245           return;
00246      }
00247 
00248      if ((qn1 == n1 && qn2 == n2) || (qn1 == n2 && qn2 == n1)) {
00249           n3 = qn0;
00250           n4 = qn3;
00251           return;
00252      }
00253 
00254      if ((qn2 == n1 && qn3 == n2) || (qn2 == n2 && qn3 == n1)) {
00255           n3 = qn0;
00256           n4 = qn1;
00257           return;
00258      }
00259 
00260      if ((qn3 == n1 && qn0 == n2) || (qn3 == n2 && qn0 == n1)) {
00261           n3 = qn1;
00262           n4 = qn2;
00263           return;
00264      }
00265 
00266      cout << " Warning: You should not come here " << endl;
00267      cout << " search for " << n1 << "  " << n2 << endl;
00268      exit(0);
00269 }
00271 
00272 PFace
00273 Face::create_quad(const PFace t1, const PFace t2, int replace)
00274 {
00275      NodeSequence connect;
00276      PNode commonnodes[3];
00277 
00278      connect = t1->getNodes();
00279 
00280      int index = 0;
00281      for (int i = 0; i < 3; i++) {
00282           if (t2->hasNode(connect[i]))
00283                commonnodes[index++] = connect[i];
00284      }
00285 
00286      if(index != 2) {
00287           cout << "Fatal Error:  Faces don't match for comoon edge " << endl;
00288           exit(0);
00289      }
00290 
00291      PNode ot1 = Face::opposite_node(t1, commonnodes[0], commonnodes[1]);
00292      PNode ot2 = Face::opposite_node(t2, commonnodes[0], commonnodes[1]);
00293 
00294      connect.resize(4);
00295      connect[0] = ot1;
00296      connect[1] = commonnodes[0];
00297      connect[2] = ot2;
00298      connect[3] = commonnodes[1];
00299 
00300      if (!replace) {
00301           Face *qface = new Face;
00302           qface->setNodes(connect);
00303           return qface;
00304      }
00305 
00306      t1->setNodes(connect);
00307      t2->setStatus(MeshEntity::REMOVE);
00308      return t1;
00309 }
00310 
00312 
00313 int
00314 Face::hexagon_2_quads(const NodeSequence &hexnodes, FaceSequence &quads, int offset)
00315 {
00316 #ifdef DEBUG
00317      if( hexnodes.size() != 6) {
00318           cout << "Fatal Error: Given polygon is not hex " << endl;
00319           exit(0);
00320      }
00321 #endif
00322 
00323      PFace face1 = Face::newObject();
00324      PFace face2 = Face::newObject();
00325 
00326      NodeSequence connect(4);
00327 
00328      for (int i = 0; i < 3; i++) {
00329           connect[0] = hexnodes[ (i + offset + 0) % 6];
00330           connect[1] = hexnodes[ (i + offset + 1) % 6];
00331           connect[2] = hexnodes[ (i + offset + 2) % 6];
00332           connect[3] = hexnodes[ (i + offset + 3) % 6];
00333           face1->setNodes(connect);
00334 
00335           connect[0] = hexnodes[ (i + offset + 3) % 6];
00336           connect[1] = hexnodes[ (i + offset + 4) % 6];
00337           connect[2] = hexnodes[ (i + offset + 5) % 6];
00338           connect[3] = hexnodes[ (i + offset + 6) % 6];
00339           face2->setNodes(connect);
00340           if (face1->isConvex() && face2->isConvex()) {
00341                quads.resize(2);
00342                quads[0] = face1;
00343                quads[1] = face2;
00344                return 0;
00345           }
00346      }
00347      quads.clear();
00348      delete face1;
00349      delete face2;
00350      return 1;
00351 }
00352 
00354 
00355 void
00356 Face::getAvgPos( Point3D &pc) const
00357 {
00358      pc[0] = 0.0;
00359      pc[1] = 0.0;
00360      pc[2] = 0.0;
00361 
00362      int nsize = connect.size();
00363      for (int inode = 0; inode < nsize; inode++) {
00364           Vertex *v = connect[inode];
00365           const Point3D &p3d = v->getXYZCoords();
00366           pc[0] += p3d[0];
00367           pc[1] += p3d[1];
00368           pc[2] += p3d[2];
00369      }
00370 
00371      double multby = 1.0/(double)nsize;
00372      pc[0] *=  multby;
00373      pc[1] *=  multby;
00374      pc[2] *=  multby;
00375 }
00376 
00378 
00379 Vec3D
00380 Face::normal(const Point3D &p0, const Point3D &p1, const Point3D &p2)
00381 {
00382      Vec3D normal, p1p0, p2p0;
00383      Math::create_vector(p2, p0, p2p0);
00384      Math::create_vector(p1, p0, p1p0);
00385      Math::cross_product(p2p0, p1p0, normal);
00386 
00387      double mag = Math::magnitude(normal);
00388      normal[0] /= mag;
00389      normal[1] /= mag;
00390      normal[2] /= mag;
00391 
00392      return normal;
00393 }
00394 
00396 
00397 Vec3D
00398 Face::normal(const Vertex *v0, const Vertex *v1, const Vertex *v2)
00399 {
00400      Vec3D n;
00401      Math::normal(v0->getXYZCoords(), v1->getXYZCoords(), v2->getXYZCoords(), n);
00402      return n;
00403 }
00404 
00406 
00407 double
00408 Face::tri_area(const Point3D &p0, const Point3D &p1, const Point3D &p2)
00409 {
00411      // Ref: http://mathworld.wolfram.com/HeronsFormula.html
00412      //
00413      // Heron's formular  A = sqrt(s(s-a)(s-b)(s-c)) is expensive because if require
00414      // three square roots for each a,b, and c.
00415      // Instead we will use alternate formula of "Heron" which avoids three
00416      // expensive square roots.
00417      //
00418      // Sorting is done to reduce the truncation errors. Very similar to Kahan's
00419      // Original idea.
00421 
00422      double d[3];
00423      d[0] = Math::length2(p1, p2);
00424      d[1] = Math::length2(p2, p0);
00425      d[2] = Math::length2(p0, p1);
00426 
00427      std::sort(d, d + 3); // May be we should have optimized version than STL one
00428 
00429      double a2 = d[0];
00430      double b2 = d[1];
00431      double c2 = d[2];
00432 
00433      double area = 0.25 * sqrt(4 * a2 * b2 - (a2 + b2 - c2)*(a2 + b2 - c2));
00434      return area;
00435 }
00436 
00438 
00439 double
00440 Face::quad_area(const Point3D &p0, const Point3D &p1,
00441                 const Point3D &p2, const Point3D &p3)
00442 {
00444      // For explanation of some amazing proofs and theorems, please refer to the
00445      // following site. This implementation is based on this article.
00446      //
00447      // http://softsurfer.com/Archive/algorithm_0101/algorithm_0101.htm#Quadrilaterals
00448      //
00449      // For Bretschneider's Formula: refer to
00450      // http://mathworld.wolfram.com/BretschneidersFormula.html
00451      // Given a general quadrilateral with sides of lengths a, b, c, and d, the area
00452      // is given by
00453      //              K = 1/4sqrt(4p^2q^2-(b^2+d^2-a^2-c^2)^2)
00454      //
00455      // It seems that the first method is better because it doesn't require
00456      // expensive 6 lenghts. ( 4 sides + 2 diagonal ). But a proper analysis needs
00457      // be done. But probably, the most amazing thing about the formula is that
00458      // it is valid for 3D quadrilateral and for both convex and concave.
00459      // But why it handles concavity, I am not quite sure.
00460      //
00461      // Chaman Singh Verma
00462      // 16th Feb 2010.
00464 
00465      Vec3D d0d1, v2v0, v3v1;
00466      Math::create_vector(p2, p0, v2v0);
00467      Math::create_vector(p3, p1, v3v1);
00468      Math::cross_product(v2v0, v3v1, d0d1);
00469 
00470      double area = 0.5 * Math::magnitude(d0d1);
00471      return area;
00472 }
00474 
00475 bool
00476 Face::is_convex_quad(const Point3D &p0, const Point3D &p1,
00477                      const Point3D &p2, const Point3D &p3)
00478 {
00479      double qarea = quad_area(p0, p1, p2, p3);
00480 
00481      double tarea1, tarea2;
00482      tarea1 = tri_area(p0, p1, p2);
00483      tarea2 = tri_area(p0, p2, p3);
00484      if (fabs(tarea1 + tarea2 - qarea) > 1.0E-10) return 0;
00485 
00486      tarea1 = tri_area(p0, p1, p3);
00487      tarea2 = tri_area(p1, p2, p3);
00488      if (fabs(tarea1 + tarea2 - qarea) > 1.0E-10) return 0;
00489 
00490      return 1;
00491 }
00492 
00494 
00495 int
00496 Face::getRelations( NodeSequence &vneighs)
00497 {
00498      vneighs.clear();
00499      set<Vertex*> vset;
00500 
00501      int nnodes = connect.size();
00502 
00503      for (int i = 0; i < nnodes; i++) {
00504           Vertex *vertex = connect[i];
00505           vertex->getRelations( vneighs );
00506           int numneighs = vneighs.size();
00507           for (int j = 0; j < numneighs; j++)
00508                vset.insert(vneighs[j]);
00509      }
00510 
00511      for (int i = 0; i < nnodes; i++)
00512           vset.erase(connect[i]);
00513 
00514      NodeSequence vresult;
00515      if (!vset.empty()) {
00516           set<Vertex*>::const_iterator it;
00517           size_t index = 0;
00518           for (it = vset.begin(); it != vset.end(); ++it)
00519                vresult[index++] = *it;
00520      }
00521      return 0;
00522 }
00523 
00525 
00526 double
00527 Face::getAspectRatio()
00528 {
00529      int nSize = connect.size();
00530 
00531      double minlen = std::numeric_limits< double >::max();
00532      double maxlen = 0.0;
00533 
00534      for (int i = 0; i < nSize; i++) {
00535           Vertex *v0 = connect[(i + 0) % nSize];
00536           Vertex *v1 = connect[(i + 1) % nSize];
00537           double len2 = Vertex::length2(v0, v1);
00538           if (len2 > maxlen) maxlen = len2;
00539           if (len2 < minlen) minlen = len2;
00540      }
00541 
00542      return sqrt(minlen / maxlen);
00543 }
00544 
00546 
00547 int
00548 Face::getRelations02( FaceSequence &faceneighs)
00549 {
00550      faceneighs.clear();
00551 
00552      FaceSequence vneighs;
00553 
00554      int nSize = connect.size();
00555      for (int i = 0; i < nSize; i++) {
00556           Vertex *v0 = connect[i];
00557           v0->getRelations( vneighs );
00558           int numneighs = vneighs.size();
00559           for (int  j = 0; j < numneighs; j++) {
00560                if (vneighs[j] != this) {
00561                     if (find(faceneighs.begin(), faceneighs.end(), vneighs[j])
00562                               == faceneighs.end())
00563                          faceneighs.push_back(vneighs[j]);
00564                }
00565           }
00566      }
00567      return 0;
00568 }
00570 
00571 int
00572 Face::getRelations12( FaceSequence &faceneighs)
00573 {
00574      faceneighs.clear();
00575      FaceSequence edgeneighs;
00576 
00577      int nSize = connect.size();
00578      for (int i = 0; i < nSize; i++) {
00579           Vertex *v0 = getNodeAt(i + 0);
00580           Vertex *v1 = getNodeAt(i + 1);
00581           Mesh::getRelations112(v0, v1, edgeneighs);
00582           int numneighs = edgeneighs.size();
00583           for (int j = 0; j <  numneighs; j++) {
00584                if (edgeneighs[j] != this) {
00585                     if (find(faceneighs.begin(), faceneighs.end(), edgeneighs[j])
00586                               == faceneighs.end())
00587                          faceneighs.push_back(edgeneighs[j]);
00588                }
00589           }
00590      }
00591      return 0;
00592 }
00593 
00595 
00596 bool
00597 Face::isValid() const
00598 {
00599      if (!isActive()) return 0;
00600 
00601      int nSize = connect.size();
00602      for (int i = 0; i < nSize; i++)
00603           if (!connect[i]->isActive()) return 0;
00604 
00605      for (int i = 0; i < nSize; i++) {
00606           int ncount = 0;
00607           for (int j = 0; j < nSize; j++)
00608                if (connect[i] == connect[j]) ncount++;
00609           if (ncount != 1) return 0;
00610      }
00611      return 1;
00612 }
00613 
00615 
00616 bool
00617 Face::has_boundary_edge() const
00618 {
00619      int nSize = connect.size();
00620      FaceSequence neighs;
00621      for (int i = 0; i < nSize; i++) {
00622           Vertex *v0 = getNodeAt(i + 0);
00623           Vertex *v1 = getNodeAt(i + 1);
00624           if (v0->isBoundary() && v1->isBoundary()) {
00625                Mesh::getRelations112(v0, v1, neighs);
00626                if (neighs.size() == 1) return 1;
00627           }
00628      }
00629      return 0;
00630 }
00632 
00633 void
00634 Face::bilinear_weights(double xi, double eta, vector<double> &weight)
00635 {
00636      weight.resize(4);
00637 
00638 #ifdef DEBUG
00639      assert(xi >= -1.0 && xi <= 1.0);
00640      assert(eta >= -1.0 && eta <= 1.0);
00641 #endif
00642 
00643      double coeff = 1.0 / 4.0;
00644 
00645      weight[0] = coeff * (1.0 - xi)*(1.0 - eta);
00646      weight[1] = coeff * (1.0 + xi)*(1.0 - eta);
00647      weight[2] = coeff * (1.0 + xi)*(1.0 + eta);
00648      weight[3] = coeff * (1.0 - xi)*(1.0 + eta);
00649 }
00650 
00652 
00653 double
00654 Face::linear_interpolation(const vector<double> &x, const vector<double> &w)
00655 {
00656      size_t numNodes = x.size();
00657      assert(w.size() == numNodes);
00658 
00659      double sum = 0.0;
00660      for (size_t i = 0; i < numNodes; i++) sum += x[i] * w[i];
00661 
00662      return sum;
00663 }
00664 
00666 
00667 int
00668 Face::concaveAt() const
00669 {
00670      int nsize = getSize(0);
00671      Vertex *vertex;
00672      double x[5], y[5], z[5], triarea;
00673 
00674      Point3D xyz;
00675 
00676      for (int i = 0; i < nsize; i++) {
00677           vertex = getNodeAt(i);
00678           xyz = vertex->getXYZCoords();
00679           x[0] = xyz[0];
00680           y[0] = xyz[1];
00681           z[0] = xyz[2];
00682 
00683           vertex = getNodeAt(i + 1 );
00684           xyz = vertex->getXYZCoords();
00685           x[1] = xyz[0];
00686           y[1] = xyz[1];
00687           z[1] = xyz[2];
00688 
00689           vertex = getNodeAt(i + nsize - 1 );
00690           xyz = vertex->getXYZCoords();
00691           x[2] = xyz[0];
00692           y[2] = xyz[1];
00693           z[2] = xyz[2];
00694 
00695           triarea = PolygonArea3D(3, x, y, z);
00696           if (triarea < 0.0) return i;
00697      }
00698      return -1;
00699 }
00700 
00702 double Vertex:: point_orient( const Point3D &p0, const Point3D &p1, const Point3D &qpoint)
00703 {
00704      double x[5], y[5], z[5];
00705 
00706      x[0] = p0[0];
00707      x[1] = p1[0];
00708      x[2] = qpoint[0];
00709 
00710      y[0] = p0[1];
00711      y[1] = p1[1];
00712      y[2] = qpoint[1];
00713 
00714      z[0] = p0[2];
00715      z[1] = p1[2];
00716      z[2] = qpoint[2];
00717 
00718      return PolygonArea3D(3, x, y, z);
00719 }
00721 
00722 bool
00723 Face::isSimple() const
00724 {
00725      int nnodes =  getSize(0);
00726      if( nnodes == 3 ) return 1;
00727 
00728      double d0, d1, d2, d3;
00729      if( nnodes == 4 ) {
00730           const Point3D &v0 = getNodeAt(0)->getXYZCoords();
00731           const Point3D &v1 = getNodeAt(1)->getXYZCoords();
00732           const Point3D &v2 = getNodeAt(2)->getXYZCoords();
00733           const Point3D &v3 = getNodeAt(3)->getXYZCoords();
00734 
00735           d0  = Vertex::point_orient( v0, v1, v2);
00736           d1  = Vertex::point_orient( v0, v1, v3);
00737 
00738           d2  = Vertex::point_orient( v2, v3, v0);
00739           d3  = Vertex::point_orient( v2, v3, v1);
00740 
00741           if( (d0*d1 < 0.0) && (d2*d3 < 0.0 ) ) return 0;
00742 
00743           d0  = Vertex::point_orient( v1, v2, v0);
00744           d1  = Vertex::point_orient( v1, v2, v3);
00745 
00746           d2  = Vertex::point_orient( v0, v3, v1);
00747           d3  = Vertex::point_orient( v0, v3, v2);
00748 
00749           if( (d0*d1 < 0.0) && (d2*d3 < 0.0 ) ) return 0;
00750 
00751           return 1;
00752      }
00753 
00754      cout << "Warning: General Polygons not supported yet " << endl;
00755 
00756      return 0;
00757 
00758 }
00759 
00761 
00762 int
00763 Face::refine_quad15(NodeSequence &newnodes, FaceSequence &newfaces)
00764 {
00765      if (!isConvex())
00766           return refine_concave_quad15(newnodes, newfaces);
00767 
00768      return refine_convex_quad15(newnodes, newfaces);
00769 }
00770 
00772 
00773 int
00774 Face::refine_concave_quad15(NodeSequence &newnodes, FaceSequence &newfaces)
00775 {
00776      newnodes.resize(4);
00777      newfaces.resize(5);
00778 
00779      NodeSequence rotatedNodes(4), tconnect(3);
00780 
00781      Point3D xyz;
00782 
00783      quad_tessalate(connect, rotatedNodes);
00784      Vertex::mid_point(rotatedNodes[0], rotatedNodes[2], xyz, 1.0 / 3.0);
00785      newnodes[0] = Vertex::newObject();
00786      newnodes[0]->setXYZCoords(xyz);
00787 
00788      Vertex::mid_point(rotatedNodes[0], rotatedNodes[2], xyz, 2.0 / 3.0);
00789      newnodes[2] = Vertex::newObject();
00790      newnodes[2]->setXYZCoords(xyz);
00791 
00792      Face triface;
00793      tconnect[0] = rotatedNodes[0];
00794      tconnect[1] = rotatedNodes[1];
00795      tconnect[2] = rotatedNodes[2];
00796      triface.setNodes(tconnect);
00797      triface.getAvgPos( xyz );
00798      newnodes[1] = Vertex::newObject();
00799      newnodes[1]->setXYZCoords(xyz);
00800 
00801      tconnect[0] = rotatedNodes[0];
00802      tconnect[1] = rotatedNodes[2];
00803      tconnect[2] = rotatedNodes[3];
00804      triface.setNodes(tconnect);
00805      triface.getAvgPos( xyz );
00806      newnodes[3] = Vertex::newObject();
00807      newnodes[3]->setXYZCoords(xyz);
00808 
00809      NodeSequence connect(4);
00810 
00811      connect[0] = rotatedNodes[0];
00812      connect[1] = rotatedNodes[1];
00813      connect[2] = newnodes[1];
00814      connect[3] = newnodes[0];
00815      newfaces[0] = Face::newObject();
00816      newfaces[0]->setNodes(connect);
00817 
00818      connect[0] = rotatedNodes[1];
00819      connect[1] = rotatedNodes[2];
00820      connect[2] = newnodes[2];
00821      connect[3] = newnodes[1];
00822      newfaces[1] = Face::newObject();
00823      newfaces[1]->setNodes(connect);
00824 
00825      connect[0] = rotatedNodes[2];
00826      connect[1] = rotatedNodes[3];
00827      connect[2] = newnodes[3];
00828      connect[3] = newnodes[2];
00829      newfaces[2] = Face::newObject();
00830      newfaces[2]->setNodes(connect);
00831 
00832      connect[0] = rotatedNodes[3];
00833      connect[1] = rotatedNodes[0];
00834      connect[2] = newnodes[0];
00835      connect[3] = newnodes[3];
00836      newfaces[3] = Face::newObject();
00837      newfaces[3]->setNodes(connect);
00838 
00839      connect[0] = newnodes[0];
00840      connect[1] = newnodes[1];
00841      connect[2] = newnodes[2];
00842      connect[3] = newnodes[3];
00843      newfaces[4] = Face::newObject();
00844      newfaces[4]->setNodes(connect);
00845 
00846      return 0;
00847 }
00848 
00850 
00851 int
00852 Face::refine_convex_quad15(NodeSequence &newnodes, FaceSequence &newfaces)
00853 {
00854      newnodes.resize(4);
00855      newfaces.resize(5);
00856 
00857      Point3D xyz;
00858      vector<double> weight;
00859      vector<double> xc(4), yc(4), zc(4);
00860      for (int i = 0; i < 4; i++) {
00861           Vertex *v = this->getNodeAt(i);
00862           xyz = v->getXYZCoords();
00863           xc[i] = xyz[0];
00864           yc[i] = xyz[1];
00865           zc[i] = xyz[2];
00866      }
00867 
00868      double dl = 2.0 / 3.0;
00869 
00870      bilinear_weights(-1.0 + dl, -1.0 + dl, weight);
00871      xyz[0] = linear_interpolation(xc, weight);
00872      xyz[1] = linear_interpolation(yc, weight);
00873      xyz[2] = linear_interpolation(zc, weight);
00874      newnodes[0] = Vertex::newObject();
00875      newnodes[0]->setXYZCoords(xyz);
00876 
00877      bilinear_weights(-1.0 + 2 * dl, -1.0 + dl, weight);
00878      xyz[0] = linear_interpolation(xc, weight);
00879      xyz[1] = linear_interpolation(yc, weight);
00880      xyz[2] = linear_interpolation(zc, weight);
00881      newnodes[1] = Vertex::newObject();
00882      newnodes[1]->setXYZCoords(xyz);
00883 
00884      bilinear_weights(-1.0 + 2 * dl, -1.0 + 2 * dl, weight);
00885      xyz[0] = linear_interpolation(xc, weight);
00886      xyz[1] = linear_interpolation(yc, weight);
00887      xyz[2] = linear_interpolation(zc, weight);
00888      newnodes[2] = Vertex::newObject();
00889      newnodes[2]->setXYZCoords(xyz);
00890 
00891      bilinear_weights(-1.0 + dl, -1.0 + 2.0 * dl, weight);
00892      xyz[0] = linear_interpolation(xc, weight);
00893      xyz[1] = linear_interpolation(yc, weight);
00894      xyz[2] = linear_interpolation(zc, weight);
00895      newnodes[3] = Vertex::newObject();
00896      newnodes[3]->setXYZCoords(xyz);
00897 
00898      NodeSequence connect(4), oldnodes;
00899      oldnodes = this->getNodes();
00900 
00901      connect[0] = oldnodes[0];
00902      connect[1] = oldnodes[1];
00903      connect[2] = newnodes[1];
00904      connect[3] = newnodes[0];
00905      newfaces[0] = Face::newObject();
00906      newfaces[0]->setNodes(connect);
00907 
00908      connect[0] = oldnodes[1];
00909      connect[1] = oldnodes[2];
00910      connect[2] = newnodes[2];
00911      connect[3] = newnodes[1];
00912      newfaces[1] = Face::newObject();
00913      newfaces[1]->setNodes(connect);
00914 
00915      connect[0] = oldnodes[2];
00916      connect[1] = oldnodes[3];
00917      connect[2] = newnodes[3];
00918      connect[3] = newnodes[2];
00919      newfaces[2] = Face::newObject();
00920      newfaces[2]->setNodes(connect);
00921 
00922      connect[0] = oldnodes[3];
00923      connect[1] = oldnodes[0];
00924      connect[2] = newnodes[0];
00925      connect[3] = newnodes[3];
00926      newfaces[3] = Face::newObject();
00927      newfaces[3]->setNodes(connect);
00928 
00929      connect[0] = newnodes[0];
00930      connect[1] = newnodes[1];
00931      connect[2] = newnodes[2];
00932      connect[3] = newnodes[3];
00933      newfaces[4] = Face::newObject();
00934      newfaces[4]->setNodes(connect);
00935 
00936      return 0;
00937 }
00938 
00940 
00941 Vertex *
00942 Face::centroid(const Vertex *v0, const Vertex *v1, const Vertex *v2)
00943 {
00944      Vertex *vc = Vertex::newObject();
00945 
00946      Point3D pc, xyz;
00947      pc[0] = 0.0;
00948      pc[1] = 0.0;
00949      pc[2] = 0.0;
00950 
00951      xyz = v0->getXYZCoords();
00952      pc[0] += xyz[0];
00953      pc[1] += xyz[1];
00954      pc[2] += xyz[2];
00955 
00956      xyz = v1->getXYZCoords();
00957      pc[0] += xyz[0];
00958      pc[1] += xyz[1];
00959      pc[2] += xyz[2];
00960 
00961      xyz = v2->getXYZCoords();
00962      pc[0] += xyz[0];
00963      pc[1] += xyz[1];
00964      pc[2] += xyz[2];
00965 
00966      pc[0] /= 3.0;
00967      pc[1] /= 3.0;
00968      pc[2] /= 3.0;
00969 
00970      vc->setXYZCoords(pc);
00971      return vc;
00972 }
00973 
00975 
00976 Vertex *
00977 Face::centroid(const Vertex *v0, const Vertex *v1, const Vertex *v2,
00978                const Vertex *v3, const Vertex *v4)
00979 {
00980      Vertex *vc = Vertex::newObject();
00981 
00982      Point3D pc, xyz;
00983      pc[0] = 0.0;
00984      pc[1] = 0.0;
00985      pc[2] = 0.0;
00986 
00987      xyz = v0->getXYZCoords();
00988      pc[0] += xyz[0];
00989      pc[1] += xyz[1];
00990      pc[2] += xyz[2];
00991 
00992      xyz = v1->getXYZCoords();
00993      pc[0] += xyz[0];
00994      pc[1] += xyz[1];
00995      pc[2] += xyz[2];
00996 
00997      xyz = v2->getXYZCoords();
00998      pc[0] += xyz[0];
00999      pc[1] += xyz[1];
01000      pc[2] += xyz[2];
01001 
01002      xyz = v3->getXYZCoords();
01003      pc[0] += xyz[0];
01004      pc[1] += xyz[1];
01005      pc[2] += xyz[2];
01006 
01007      xyz = v4->getXYZCoords();
01008      pc[0] += xyz[0];
01009      pc[1] += xyz[1];
01010      pc[2] += xyz[2];
01011 
01012      pc[0] /= 5.0;
01013      pc[1] /= 5.0;
01014      pc[2] /= 5.0;
01015 
01016      vc->setXYZCoords(pc);
01017      return vc;
01018 }
01019 
01021 
01022 int
01023 Face::is_3_sided_convex_loop_quad_meshable(const int *segments, int *partsegments)
01024 {
01025      double M[6][6], rhs[6];
01026 
01027      // octave:1> A = [ 0  0 -1  1  0  0;
01028      //                -1  0  0  0  1  0;
01029      //                 0 -1  0  0  0  1;
01030      //                1  0  0  1  0  0;
01031      //                 0  1  0  0  1  0;
01032      //                 0  0  1  0  0  1]
01033      //A =
01034      //
01035      //   0   0  -1   1   0   0
01036      //  -1   0   0   0   1   0
01037      //   0  -1   0   0   0   1
01038      //   1   0   0   1   0   0
01039      //   0   1   0   0   1   0
01040      //   0   0   1   0   0   1
01041      //  octave:2> inv(A)
01042      //  ans =
01043 
01044      //  -0.5  -0.5   0.5   0.5   0.5  -0.5
01045      //   0.5  -0.5  -0.5  -0.5   0.5   0.5
01046      //  -0.5   0.5  -0.5   0.5  -0.5   0.5
01047      //   0.5   0.5  -0.5   0.5  -0.5   0.5
01048      //  -0.5   0.5   0.5   0.5   0.5  -0.5
01049      //   0.5  -0.5   0.5  -0.5   0.5   0.5
01050 
01051      //   First Row
01052      //  -0.5  -0.5   0.5   0.5   0.5  -0.5
01053      M[0][0] = -0.5;
01054      M[0][1] = -0.5;
01055      M[0][2] = 0.5;
01056      M[0][3] = 0.5;
01057      M[0][4] = 0.5;
01058      M[0][5] = -0.5;
01059 
01060      //   Second Row
01061      //   0.5  -0.5  -0.5  -0.5   0.5   0.5
01062      M[1][0] = 0.5;
01063      M[1][1] = -0.5;
01064      M[1][2] = -0.5;
01065      M[1][3] = -0.5;
01066      M[1][4] = 0.5;
01067      M[1][5] = 0.5;
01068 
01069      //  Third Row
01070      //  -0.5   0.5  -0.5   0.5  -0.5   0.5
01071      M[2][0] = -0.5;
01072      M[2][1] = 0.5;
01073      M[2][2] = -0.5;
01074      M[2][3] = 0.5;
01075      M[2][4] = -0.5;
01076      M[2][5] = 0.5;
01077 
01078      // Forth Row
01079      //   0.5   0.5  -0.5   0.5  -0.5   0.5
01080      M[3][0] = 0.5;
01081      M[3][1] = 0.5;
01082      M[3][2] = -0.5;
01083      M[3][3] = 0.5;
01084      M[3][4] = -0.5;
01085      M[3][5] = 0.5;
01086 
01087      //   Fifth Row
01088      //  -0.5   0.5   0.5   0.5   0.5  -0.5
01089      M[4][0] = -0.5;
01090      M[4][1] = 0.5;
01091      M[4][2] = 0.5;
01092      M[4][3] = 0.5;
01093      M[4][4] = 0.5;
01094      M[4][5] = -0.5;
01095 
01096      //  Sixth Row
01097      //   0.5  -0.5   0.5  -0.5   0.5   0.5
01098      M[5][0] = 0.5;
01099      M[5][1] = -0.5;
01100      M[5][2] = 0.5;
01101      M[5][3] = -0.5;
01102      M[5][4] = 0.5;
01103      M[5][5] = 0.5;
01104 
01105      rhs[0] = 0.0;
01106      rhs[1] = 0.0;
01107      rhs[2] = 0.0;
01108 
01109      rhs[3] = segments[0];
01110      rhs[4] = segments[1];
01111      rhs[5] = segments[2];
01112 
01113      vector<int> x(6);
01114      for (int i = 0; i < 6; i++) {
01115           double sum = 0.0;
01116           for (int j = 0; j < 6; j++)
01117                sum += M[i][j] * rhs[j];
01118           x[i] = (int) sum;
01119      }
01120 
01121      for (int i = 0; i < 6; i++)
01122           if (x[i] < 1) return 0;
01123 
01124      if (x[0] + x[3] != rhs[3]) return 0;
01125      partsegments[0] = x[0];
01126      partsegments[1] = x[3];
01127 
01128      if (x[1] + x[4] != rhs[4]) return 0;
01129      partsegments[2] = x[1];
01130      partsegments[3] = x[4];
01131 
01132      if (x[2] + x[5] != rhs[5]) return 0;
01133      partsegments[4] = x[2];
01134      partsegments[5] = x[5];
01135 
01136      return 1;
01137 }
01138 
01140 
01141 int
01142 Face::is_5_sided_convex_loop_quad_meshable(const int *segments, int *partSegments)
01143 {
01144      double M[10][10], rhs[10];
01145 
01146      //  Equations:
01147      //      b0 -a2   = 0
01148      //      b1 -a3   = 0
01149      //      b2 -a4   = 0
01150      //      b3 -a0   = 0
01151      //      b4 -a1   = 0
01152      //      a0 + b0  = s0
01153      //      a1 + b1  = s1
01154      //      a2 + b2  = s2
01155      //      a3 + b3  = s3
01156      //      a4 + b4  = s4
01157 
01158      // For more details; See Guy Bunin's paper.
01159      // M =
01160      //   a0  a1 a2  a3  a4  b0   b1 b2   b3  b4
01161      //   0   0  -1   0   0   1   0   0   0   0
01162      //   0   0   0  -1   0   0   1   0   0   0
01163      //   0   0   0   0  -1   0   0   1   0   0
01164      //  -1   0   0   0   0   0   0   0   1   0
01165      //   0  -1   0   0   0   0   0   0   0   1
01166      //   1   0   0   0   0   1   0   0   0   0
01167      //   0   1   0   0   0   0   1   0   0   0
01168      //   0   0   1   0   0   0   0   1   0   0
01169      //   0   0   0   1   0   0   0   0   1   0
01170      //   0   0   0   0   1   0   0   0   0   1
01171 
01172      // octave:38> inv(M)
01173      // ans =
01174      //  -0.5   0.5   0.5  -0.5  -0.5   0.5  -0.5  -0.5   0.5   0.5
01175      //  -0.5  -0.5   0.5   0.5  -0.5   0.5   0.5  -0.5  -0.5   0.5
01176      //  -0.5  -0.5  -0.5   0.5   0.5   0.5   0.5   0.5  -0.5  -0.5
01177      //   0.5  -0.5  -0.5  -0.5   0.5  -0.5   0.5   0.5   0.5  -0.5
01178      //   0.5   0.5  -0.5  -0.5  -0.5  -0.5  -0.5   0.5   0.5   0.5
01179      //   0.5  -0.5  -0.5   0.5   0.5   0.5   0.5   0.5  -0.5  -0.5
01180      //   0.5   0.5  -0.5  -0.5   0.5  -0.5   0.5   0.5   0.5  -0.5
01181      //   0.5   0.5   0.5  -0.5  -0.5  -0.5  -0.5   0.5   0.5   0.5
01182      //  -0.5   0.5   0.5   0.5  -0.5   0.5  -0.5  -0.5   0.5   0.5
01183      //  -0.5  -0.5   0.5   0.5   0.5   0.5   0.5  -0.5  -0.5   0.5
01184 
01185 
01186      //   First Row
01187      //  -0.5   0.5   0.5  -0.5  -0.5   0.5  -0.5  -0.5   0.5   0.5
01188      M[0][0] = -0.5;
01189      M[0][1] = 0.5;
01190      M[0][2] = 0.5;
01191      M[0][3] = -0.5;
01192      M[0][4] = -0.5;
01193      M[0][5] = 0.5;
01194      M[0][6] = -0.5;
01195      M[0][7] = -0.5;
01196      M[0][8] = 0.5;
01197      M[0][9] = 0.5;
01198 
01199      //  -0.5  -0.5   0.5   0.5  -0.5   0.5   0.5  -0.5  -0.5   0.5
01200      M[1][0] = -0.5;
01201      M[1][1] = -0.5;
01202      M[1][2] = 0.5;
01203      M[1][3] = 0.5;
01204      M[1][4] = -0.5;
01205      M[1][5] = 0.5;
01206      M[1][6] = 0.5;
01207      M[1][7] = -0.5;
01208      M[1][8] = -0.5;
01209      M[1][9] = 0.5;
01210 
01211 
01212      //  -0.5  -0.5  -0.5   0.5   0.5   0.5   0.5   0.5  -0.5  -0.5
01213      M[2][0] = -0.5;
01214      M[2][1] = -0.5;
01215      M[2][2] = -0.5;
01216      M[2][3] = 0.5;
01217      M[2][4] = 0.5;
01218      M[2][5] = 0.5;
01219      M[2][6] = 0.5;
01220      M[2][7] = 0.5;
01221      M[2][8] = -0.5;
01222      M[2][9] = -0.5;
01223 
01224      //   0.5  -0.5  -0.5  -0.5   0.5  -0.5   0.5   0.5   0.5  -0.5
01225      M[3][0] = 0.5;
01226      M[3][1] = -0.5;
01227      M[3][2] = -0.5;
01228      M[3][3] = -0.5;
01229      M[3][4] = 0.5;
01230      M[3][5] = -0.5;
01231      M[3][6] = 0.5;
01232      M[3][7] = 0.5;
01233      M[3][8] = 0.5;
01234      M[3][9] = -0.5;
01235 
01236      //   0.5   0.5  -0.5  -0.5  -0.5  -0.5  -0.5   0.5   0.5   0.5
01237      M[4][0] = 0.5;
01238      M[4][1] = 0.5;
01239      M[4][2] = -0.5;
01240      M[4][3] = -0.5;
01241      M[4][4] = -0.5;
01242      M[4][5] = -0.5;
01243      M[4][6] = -0.5;
01244      M[4][7] = 0.5;
01245      M[4][8] = 0.5;
01246      M[4][9] = 0.5;
01247 
01248      //   0.5  -0.5  -0.5   0.5   0.5   0.5   0.5   0.5  -0.5  -0.5
01249      M[5][0] = 0.5;
01250      M[5][1] = -0.5;
01251      M[5][2] = -0.5;
01252      M[5][3] = 0.5;
01253      M[5][4] = 0.5;
01254      M[5][5] = 0.5;
01255      M[5][6] = 0.5;
01256      M[5][7] = 0.5;
01257      M[5][8] = -0.5;
01258      M[5][9] = -0.5;
01259 
01260      //   0.5   0.5  -0.5  -0.5   0.5  -0.5   0.5   0.5   0.5  -0.5
01261      M[6][0] = 0.5;
01262      M[6][1] = 0.5;
01263      M[6][2] = -0.5;
01264      M[6][3] = -0.5;
01265      M[6][4] = 0.5;
01266      M[6][5] = -0.5;
01267      M[6][6] = 0.5;
01268      M[6][7] = 0.5;
01269      M[6][8] = 0.5;
01270      M[6][9] = -0.5;
01271 
01272      //   0.5   0.5   0.5  -0.5  -0.5  -0.5  -0.5   0.5   0.5   0.5
01273      M[7][0] = 0.5;
01274      M[7][1] = 0.5;
01275      M[7][2] = 0.5;
01276      M[7][3] = -0.5;
01277      M[7][4] = -0.5;
01278      M[7][5] = -0.5;
01279      M[7][6] = -0.5;
01280      M[7][7] = 0.5;
01281      M[7][8] = 0.5;
01282      M[7][9] = 0.5;
01283 
01284      //  -0.5   0.5   0.5   0.5  -0.5   0.5  -0.5  -0.5   0.5   0.5
01285      M[8][0] = -0.5;
01286      M[8][1] = 0.5;
01287      M[8][2] = 0.5;
01288      M[8][3] = 0.5;
01289      M[8][4] = -0.5;
01290      M[8][5] = 0.5;
01291      M[8][6] = -0.5;
01292      M[8][7] = -0.5;
01293      M[8][8] = 0.5;
01294      M[8][9] = 0.5;
01295 
01296      //  -0.5  -0.5   0.5   0.5   0.5   0.5   0.5  -0.5  -0.5   0.5
01297      M[9][0] = -0.5;
01298      M[9][1] = -0.5;
01299      M[9][2] = 0.5;
01300      M[9][3] = 0.5;
01301      M[9][4] = 0.5;
01302      M[9][5] = 0.5;
01303      M[9][6] = 0.5;
01304      M[9][7] = -0.5;
01305      M[9][8] = -0.5;
01306      M[9][9] = 0.5;
01307 
01308      rhs[0] = 0.0;
01309      rhs[1] = 0.0;
01310      rhs[2] = 0.0;
01311      rhs[3] = 0.0;
01312      rhs[4] = 0.0;
01313 
01314      rhs[5] = segments[0];
01315      rhs[6] = segments[1];
01316      rhs[7] = segments[2];
01317      rhs[8] = segments[3];
01318      rhs[9] = segments[4];
01319 
01320      vector<int> x(10);
01321      for (int i = 0; i < 10; i++) {
01322           double sum = 0.0;
01323           for (int j = 0; j < 10; j++)
01324                sum += M[i][j] * rhs[j];
01325           x[i] = (int) sum;
01326      }
01327 
01328      for (int i = 0; i < 10; i++)
01329           if (x[i] < 1) return 0;
01330 
01331      if (x[0] != x[8]) return 0;
01332      if (x[0] + x[5] != rhs[5]) return 0;
01333      partSegments[0] = x[0];
01334      partSegments[1] = x[5];
01335 
01336      if (x[1] != x[9]) return 0;
01337      if (x[1] + x[6] != rhs[6]) return 0;
01338      partSegments[2] = x[1];
01339      partSegments[3] = x[6];
01340 
01341      if (x[2] != x[5]) return 0;
01342      if (x[2] + x[7] != rhs[7]) return 0;
01343      partSegments[4] = x[2];
01344      partSegments[5] = x[7];
01345 
01346      if (x[3] != x[6]) return 0;
01347      if (x[3] + x[8] != rhs[8]) return 0;
01348      partSegments[6] = x[3];
01349      partSegments[7] = x[8];
01350 
01351      if (x[4] != x[7]) return 0;
01352      if (x[4] + x[9] != rhs[9]) return 0;
01353      partSegments[8] = x[4];
01354      partSegments[9] = x[9];
01355 
01356      return 1;
01357 }
01358 
01360 int Mesh::refine_tri_edge( Vertex *v0, Vertex *v1, int numSegments,
01361                            vector<Vertex*> &newnodes, vector<Face*> &newfaces)
01362 {
01363      assert( getAdjTable(0,2) );
01364 
01365      size_t nSize;
01366      (void) nSize;
01367      newnodes.clear();
01368      newfaces.clear();
01369      assert( numSegments >= 2);
01370 
01371      FaceSequence efaces;
01372      Mesh::getRelations112( v0, v1, efaces);
01373 
01374      int numFaces = efaces.size();
01375      if( numFaces == 0) return 1;
01376 
01377      assert( numFaces <= 2 );
01378 
01379      newnodes.reserve(numSegments-1);
01380      newfaces.reserve(2*numSegments -numFaces ); // minuse because faces will be reused..
01381 
01382      double dt = 1.0/(double)numSegments;
01383 
01384      NodeSequence edgenodes;
01385      edgenodes.reserve(numSegments+1);
01386      edgenodes.push_back(v0);
01387      Point3D xyz;
01388      for( int i = 1; i < numSegments; i++) {
01389           Vertex::mid_point(v0, v1, xyz, i*dt);
01390           Vertex *vmid = Vertex::newObject();
01391           vmid->setXYZCoords( xyz );
01392           newnodes.push_back(vmid);
01393           edgenodes.push_back(vmid);
01394           addNode( vmid );
01395      }
01396      edgenodes.push_back(v1);
01397 
01398      int pos0, pos1, pos2;
01399      Vertex *ov1 =  Face::opposite_node(efaces[0], v0, v1);
01400      assert( ov1 );
01401 
01402      deactivate( efaces[0] );
01403 
01404      nSize = edgenodes.size();
01405      vector<Vertex*> conn(3);
01406 
01407      pos0 = efaces[0]->getPosOf(v0);
01408      pos1 = efaces[0]->getPosOf(v1);
01409      pos2 = efaces[0]->getPosOf(ov1);
01410      for( int i = 0; i < numSegments; i++) {
01411           conn[pos0] = edgenodes[i];
01412           conn[pos1] = edgenodes[i+1];
01413           conn[pos2] = ov1;
01414           if( i == 0) {
01415                efaces[0]->setNodes( conn );
01416                reactivate( efaces[0] );
01417           } else {
01418                Face *f = Face::newObject();
01419                f->setNodes( conn );
01420                addFace(f);
01421                newfaces.push_back(f);
01422           }
01423      }
01424 
01425 
01426      if( numFaces > 1 ) {
01427           Vertex *ov2 =  Face::opposite_node(efaces[1], v0, v1);
01428           assert( ov2 );
01429           deactivate( efaces[1] );
01430 
01431           nSize = edgenodes.size();
01432           pos0 = efaces[1]->getPosOf(v0);
01433           pos1 = efaces[1]->getPosOf(v1);
01434           pos2 = efaces[1]->getPosOf(ov2);
01435           for( int i = 0; i < numSegments; i++) {
01436                conn[pos0] = edgenodes[i];
01437                conn[pos1] = edgenodes[i+1];
01438                conn[pos2] = ov2;
01439                if( i == 0) {
01440                     efaces[1]->setNodes( conn );
01441                     reactivate( efaces[1] );
01442                } else {
01443                     Face *f = Face::newObject();
01444                     f->setNodes( conn );
01445                     addFace(f);
01446                     newfaces.push_back(f);
01447                }
01448           }
01449      }
01450      return 0;
01451 }
01452 
01454 
01455 int Mesh::collapse_tri_edge( Vertex *v0, Vertex *v1)
01456 {
01457      assert( getAdjTable(0,2) );
01458 
01459      if( v0->isBoundary() && v1->isBoundary() ) return 1;
01460 
01461      FaceSequence efaces;
01462      Mesh::getRelations112( v0, v1, efaces);
01463 
01464      assert( efaces.size() == 2 );
01465      if( efaces[0]->isRemoved() || efaces[1]->isRemoved() ) return 1;
01466 
01467      Vertex *keepVertex   = NULL;
01468      Vertex *removeVertex = NULL;
01469 
01470      if( v0->isBoundary() ) {
01471           keepVertex   = v0;
01472           removeVertex = v1;
01473      }
01474 
01475      if( v1->isBoundary() ) {
01476           keepVertex   = v1;
01477           removeVertex = v0;
01478      }
01479 
01480      Point3D p3d;
01481      if( !v0->isBoundary()  &&  !v1->isBoundary() ) {
01482           keepVertex   = v0;
01483           removeVertex = v1;
01484           Vertex::mid_point( v0, v1, p3d);
01485           keepVertex->setXYZCoords( p3d );
01486      }
01487 
01488      remove( efaces[0] );
01489      remove( efaces[1] );
01490 
01491      FaceSequence v1faces;
01492      removeVertex->getRelations(v1faces);
01493      int nSize = v1faces.size();
01494      for( int i = 0; i < nSize; i++) {
01495           Face *f = v1faces[i];
01496           assert( !f->isRemoved() );
01497           f->replaceNode(removeVertex, keepVertex);
01498           reactivate( f );
01499      }
01500 
01501      removeVertex->clearRelations(0);
01502      removeVertex->clearRelations(2);
01503 
01504      remove( removeVertex );
01505 
01506      return 0;
01507 
01508 }
01510 int
01511 Mesh::refine_quad15(Face *face)
01512 {
01513      if( face == NULL ) return 1;
01514      if( !face->isActive() ) return 2;
01515 
01516      NodeSequence newnodes;
01517      FaceSequence newfaces;
01518      face->refine_quad15(newnodes, newfaces);
01519 
01520      int numnodes =  newnodes.size();
01521      for (int i = 0; i < numnodes; i++)
01522           addNode(newnodes[i]);
01523 
01524      int numfaces = newfaces.size();
01525      for (int i = 0; i < numfaces; i++)
01526           addFace(newfaces[i]);
01527 
01528      remove(face);
01529 
01530      return 0;
01531 }
01533 Mesh* Mesh :: getPartMesh( int p)
01534 {
01535      NodeSet vset;
01536      FaceSequence  pfaces;
01537      size_t numfaces = getSize(2);
01538      int ptid;
01539      for( size_t i = 0; i < numfaces; i++) {
01540           Face *f = getFaceAt(i);
01541           if( f->isActive() ) {
01542                f->getAttribute("Partition", ptid );
01543                if( ptid == p) {
01544                     pfaces.push_back(f);
01545                     int nnodes = f->getSize(0);
01546                     for( int j = 0; j < nnodes; j++)
01547                          vset.insert( f->getNodeAt(j));
01548                }
01549           }
01550      }
01551 
01552      Mesh *sbmesh = new Mesh;
01553      NodeSet::const_iterator it;
01554      for( it = vset.begin(); it != vset.end(); ++it)
01555           sbmesh->addNode( *it);
01556 
01557      numfaces = pfaces.size();
01558      for( size_t i = 0; i < numfaces; i++)
01559           sbmesh->addFace( pfaces[i] );
01560      return sbmesh;
01561 }
01562 
01563 
01564 int  Mesh :: getNumPartitions(int e )
01565 {
01566      int pid;
01567      if( e == 0) {
01568           set<int> vg;
01569           size_t nSize = getSize(0);
01570           for( size_t i = 0; i < nSize; i++) {
01571                Vertex *v = getNodeAt(i);
01572                if( v->isActive() ) {
01573                     v->getAttribute("Partition", pid );
01574                     vg.insert( pid );
01575                }
01576           }
01577           return vg.size();
01578      }
01579      if( e == 2) {
01580           set<int> fg;
01581           size_t nSize = getSize(2);
01582           for( size_t i = 0; i < nSize; i++) {
01583                Face *f = getFaceAt(i);
01584                if( f->isActive() ) {
01585                     f->getAttribute("Partition", pid );
01586                     fg.insert( pid );
01587                }
01588           }
01589           return fg.size();
01590      }
01591      return 0;
01592 }
01593 
01594 
01595 Jaal::NodeSequence
01596 Mesh::boundary_chain_nodes(Vertex *v0, Vertex *v1)
01597 {
01598      NodeSequence bndnodes;
01599 
01600      FaceSequence neighs;
01601      Mesh::getRelations112(v0, v1, neighs);
01602 
01603      if (neighs.size() != 2) return bndnodes;
01604 
01605      vector<Edge> bndedges;
01606 
01607      bndedges.reserve(6);
01608      Edge sharededge(v0, v1);
01609 
01610      for (int i = 0; i < 2; i++) {
01611           for (int j = 0; j < 4; j++) {
01612                Vertex *vf0 = neighs[i]->getNodeAt(j + 0);
01613                Vertex *vf1 = neighs[i]->getNodeAt(j + 1);
01614                Edge edge(vf0, vf1);
01615                if (!edge.isSameAs(sharededge))
01616                     bndedges.push_back(edge);
01617           }
01618      }
01619 
01620      int err = Mesh::make_chain(bndedges);
01621      if( err ) return bndnodes;
01622 
01623      bndnodes.reserve(6);
01624      bndnodes.push_back(bndedges[0].getNodeAt(0));
01625      bndnodes.push_back(bndedges[0].getNodeAt(1));
01626 
01627      size_t nSize = bndedges.size();
01628      for (size_t i = 1; i < nSize - 1; i++) {
01629           Vertex *v0 = bndedges[i].getNodeAt(0);
01630           Vertex *v1 = bndedges[i].getNodeAt(1);
01631           if (v0 == bndnodes[i]) {
01632                bndnodes.push_back(v1);
01633           } else if (v1 == bndnodes[i])
01634                bndnodes.push_back(v0);
01635           else {
01636                cout << "Error in bound edges : " << endl;
01637                bndnodes.clear();
01638                return bndnodes;
01639           }
01640      }
01641      return bndnodes;
01642 }
01643 
01645 
01646 int
01647 Mesh::check_convexity()
01648 {
01649      size_t numfaces = getSize(2);
01650 
01651      size_t nconvex = 0, nconcave = 0;
01652      for (size_t i = 0; i < numfaces; i++) {
01653           Face *face = getFaceAt(i);
01654           if( face->isActive() ) {
01655                if (face->isConvex())
01656                     nconvex++;
01657                else
01658                     nconcave++;
01659           }
01660      }
01661      cout << "Info: Total Number of faces " << numfaces << " #Convex : " << nconvex << " #Concave " << nconcave << endl;
01662      return 0;
01663 }
01664 
01666 
01667 void
01668 Mesh::getMinMaxFaceArea(double &minarea, double &maxarea)
01669 {
01670      size_t numfaces = getSize(2);
01671 
01672      vector<double> area(numfaces);
01673      for (size_t i = 0; i < numfaces; i++) {
01674           Face *face = getFaceAt(i);
01675           if( face->isActive() ) area[i] = face->getArea();
01676      }
01677      minarea = *min_element(area.begin(), area.end());
01678      maxarea = *max_element(area.begin(), area.end());
01679 }
01680 
01681 
01683 
01684 int
01685 Mesh::getCoordsArray( vector<double> &vcoords, vector<size_t> &l2g)
01686 {
01687      size_t numnodes = getSize(0);
01688 
01689      vcoords.clear();
01690      l2g.clear();
01691 
01692      vcoords.reserve(3 * numnodes);
01693      l2g.reserve(numnodes);
01694 
01695      for (size_t i = 0; i < numnodes; i++) {
01696           Vertex *v = getNodeAt(i);
01697           if( v->isActive() ) {
01698                l2g.push_back( v->getID() );
01699                const Point3D &xyz = v->getXYZCoords();
01700                vcoords.push_back( xyz[0] );
01701                vcoords.push_back( xyz[1] );
01702                vcoords.push_back( xyz[2] );
01703           }
01704      }
01705      return 0;
01706 }
01707 
01709 
01710 int
01711 Mesh::getNodesArray( vector<size_t> &nodearray, vector<int> &elmtopo)
01712 {
01713      elmtopo.clear();
01714      nodearray.clear();
01715 
01716      size_t numfaces = getSize(2);
01717      size_t numnodes = getSize(0);
01718 
01719      int topo = isHomogeneous();
01720      if (topo) nodearray.reserve(topo * numfaces);
01721 
01722      map<size_t,size_t> nodemap;
01723 
01724      size_t index = 0;
01725      for (size_t i = 0; i < numnodes; i++) {
01726           Vertex *v = getNodeAt(i);
01727           if( v->isActive() ) nodemap[v->getID()] = index++;
01728      }
01729 
01730      elmtopo.reserve( numfaces );
01731 
01732      for (size_t i = 0; i < numfaces; i++) {
01733           Face *face = getFaceAt(i);
01734           if( face->isActive() ) {
01735                elmtopo.push_back( face->getSize(0) );
01736                for (int j = 0; j < face->getSize(0); j++) {
01737                     Vertex *v = face->getNodeAt(j);
01738                     int lid   = nodemap[v->getID()];
01739                     nodearray.push_back(lid);
01740                }
01741           }
01742      }
01743      return 0;
01744 }
01745 
01747 
01748 Mesh*
01749 Mesh::deep_copy()
01750 {
01751      std::map<Vertex*, Vertex*> vmap;
01752 
01753      Mesh *newmesh = new Mesh;
01754      size_t numnodes = getSize(0);
01755      newmesh->reserve(numnodes, 0);
01756      for (size_t i = 0; i < numnodes; i++) {
01757           Vertex *vold = getNodeAt(i);
01758           Vertex *vnew = Vertex::newObject();
01759           vnew->setXYZCoords(vold->getXYZCoords());
01760           vmap[vold] = vnew;
01761           newmesh->addNode(vnew);
01762      }
01763 
01764      size_t numfaces = getSize(2);
01765      newmesh->reserve(numnodes, 2);
01766      NodeSequence connect;
01767      for (size_t i = 0; i < numfaces; i++) {
01768           Face *fold = getFaceAt(i);
01769           connect  = fold->getNodes();
01770           int nv   = connect.size();
01771           for (int j = 0; j < nv; j++)
01772                connect[j] = vmap[connect[j]];
01773 
01774           Face *fnew = Face::newObject();
01775           fnew->setNodes(connect);
01776           newmesh->addFace(fnew);
01777      }
01778 
01779      return newmesh;
01780 }
01781 
01783 
01784 int
01785 Mesh::setCoordsArray(const vector<double> &vcoords, const vector<size_t> &l2g)
01786 {
01787      size_t numnodes = getSize(0);
01788 
01789      Point3D xyz;
01790      int index = 0;
01791      for (size_t i = 0; i < numnodes; i++) {
01792           Vertex *v = getNodeAt(i);
01793           if( v->isActive() ) {
01794                assert( v->getID() == l2g[index] );
01795                xyz[0] = vcoords[3 * index + 0];
01796                xyz[1] = vcoords[3 * index + 1];
01797                xyz[2] = vcoords[3 * index + 2];
01798                v->setXYZCoords(xyz);
01799                index++;
01800           }
01801      }
01802      return 0;
01803 }
01804 
01806 
01807 int
01808 Mesh::make_chain(vector<Edge> &boundedges)
01809 {
01810      size_t nSize = boundedges.size();
01811 
01812      Vertex *curr_vertex = boundedges[0].getNodeAt(1);
01813 
01814      for( size_t i = 1; i < nSize; i++) {
01815           for( size_t  j = i; j < nSize; j++) {
01816                Vertex *v0 = boundedges[j].getNodeAt(0);
01817                Vertex *v1 = boundedges[j].getNodeAt(1);
01818                assert( v0->isActive() && v1->isActive() );
01819                if (v0 == curr_vertex) {
01820                     curr_vertex = v1;
01821                     if( j > i) swap( boundedges[i], boundedges[j] );
01822                     break;
01823                }
01824                if (v1 == curr_vertex) {
01825                     curr_vertex = v0;
01826                     boundedges[j].setNodes(v1, v0);
01827                     if( j > i ) swap( boundedges[i], boundedges[j] );
01828                     break;
01829                }
01830           }
01831      }
01832 
01833      if( !is_closed_chain( boundedges )  ) return 1;
01834 
01835      return 0;
01836 }
01837 
01839 
01840 int
01841 Mesh::is_closeable_chain(const vector<Edge> &boundedges)
01842 {
01843      std::map<Vertex*, set<Vertex*> > relations00;
01844 
01845      for (size_t i = 0; i < boundedges.size(); i++) {
01846           Vertex *v0 = boundedges[i].getNodeAt(0);
01847           Vertex *v1 = boundedges[i].getNodeAt(1);
01848           relations00[v0].insert(v1);
01849           relations00[v1].insert(v0);
01850      }
01851 
01852      std::map<Vertex*, set<Vertex*> > ::const_iterator it;
01853      for (it = relations00.begin(); it != relations00.end(); ++it) {
01854           Vertex *v = it->first;
01855           if (relations00[v].size() != 2) return 0;
01856      }
01857 
01858      return 1;
01859 }
01860 
01862 
01863 int
01864 Mesh::is_closed_chain(const vector<Edge> &boundedges)
01865 {
01866      Vertex *first_vertex = boundedges[0].getNodeAt(0);
01867      Vertex *curr_vertex  = boundedges[0].getNodeAt(1);
01868 
01869      for (size_t i = 1; i < boundedges.size(); i++) {
01870           if (boundedges[i].getNodeAt(0) != curr_vertex) return 0;
01871           curr_vertex = boundedges[i].getNodeAt(1);
01872      }
01873      if (curr_vertex != first_vertex) return 0;
01874 
01875      return 1;
01876 }
01878 
01879 int
01880 Mesh::rotate_chain(vector<Edge> &boundedges, Vertex *first_vertex)
01881 {
01882      if (!is_closeable_chain(boundedges)) return 1;
01883 
01884      vector<Edge>::iterator middle, it;
01885      for( it = boundedges.begin(); it != boundedges.end(); ++it) {
01886           if (it->getNodeAt(0) == first_vertex) {
01887                middle = it;
01888 
01889                break;
01890           }
01891      }
01892 
01893      if( it == boundedges.end() ) {
01894           assert( first_vertex );
01895           cout << "Fatal Error in rotate Chain : First vertex " << first_vertex->getID() << endl;
01896           exit(0);
01897      }
01898 
01899      std::rotate( boundedges.begin(), middle, boundedges.end() );
01900 
01901      assert( boundedges[0].getNodeAt(0) == first_vertex );
01902 
01903      return 0;
01904 }
01906 
01907 int
01908 Mesh::getRelations102(const PNode vtx0, const PNode vtx1, FaceSequence &faceneighs)
01909 {
01910      faceneighs.clear();
01911 
01912      assert( vtx0->isActive() && vtx1->isActive() );
01913 
01914      FaceSequence v0faces, v1faces;
01915      vtx0->getRelations( v0faces );
01916      vtx1->getRelations( v1faces );
01917 
01918      if (v0faces.empty() || v1faces.empty()) {
01919           cout << "Warning: Vertex-Faces relations are empty " << endl;
01920           return 1;
01921      }
01922 
01923      FaceSet vset;
01924      for (size_t i = 0; i < v0faces.size(); i++)
01925           vset.insert(v0faces[i]);
01926 
01927      for (size_t i = 0; i < v1faces.size(); i++)
01928           vset.insert(v1faces[i]);
01929 
01930      FaceSet::const_iterator it;
01931 
01932      if (vset.size()) {
01933           faceneighs.resize(vset.size());
01934           int index = 0;
01935           for (it = vset.begin(); it != vset.end(); ++it)
01936                faceneighs[index++] = *it;
01937      }
01938 
01939      return 0;
01940 }
01941 
01943 
01944 int
01945 Mesh::getRelations112(const PNode vtx0, const PNode vtx1, FaceSequence &faceneighs)
01946 {
01947      faceneighs.clear();
01948 
01949      assert( vtx0->isActive() && vtx1->isActive() );
01950 
01951      FaceSequence v0faces, v1faces;
01952      vtx0->getRelations( v0faces );
01953      vtx1->getRelations( v1faces );
01954 
01955      if (v0faces.empty() || v1faces.empty()) {
01956           cout << "Warning: Vertex-Faces relations are empty " << endl;
01957           return 1;
01958      }
01959 
01960      for (size_t i = 0; i < v0faces.size(); i++)
01961           assert(v0faces[i]->isActive());
01962 
01963      for (size_t i = 0; i < v0faces.size() - 1; i++)
01964           assert(v0faces[i] < v0faces[i + 1]);
01965 
01966      for (size_t i = 0; i < v1faces.size(); i++)
01967           assert(v1faces[i]->isActive());
01968 
01969      for (size_t i = 0; i < v1faces.size() - 1; i++)
01970           assert(v1faces[i] < v1faces[i + 1]);
01971 
01972      set_intersection(v0faces.begin(), v0faces.end(), v1faces.begin(),
01973                       v1faces.end(), back_inserter(faceneighs));
01974 
01975      return 0;
01976 }
01977 
01979 
01980 size_t
01981 Mesh::count_edges()
01982 {
01983      int relexist = build_relations(0, 0);
01984 
01985      size_t numnodes = getSize(0);
01986      NodeSequence neighs;
01987 
01988      size_t ncount = 0;
01989      for (size_t i = 0; i < numnodes; i++) {
01990           Vertex *vertex = nodes[i];
01991           if( vertex->isActive() )  {
01992                vertex->getRelations( neighs );
01993                size_t nSize = neighs.size();
01994                for (size_t j = 0; j < nSize; j++)
01995                     if (vertex > neighs[j]) ncount++;
01996           }
01997      }
01998 
01999      if (!relexist) clear_relations(0, 0);
02000 
02001      return ncount;
02002 }
02003 
02005 
02006 void Mesh::build_edges()
02007 {
02008      int relexist0 = build_relations(0, 0);
02009 
02010      size_t numnodes = getSize(0);
02011 
02012      edges.clear();
02013 
02014      NodeSequence neighs;
02015      for (size_t i = 0; i < numnodes; i++) {
02016           Vertex *vertex = nodes[i];
02017           if( vertex->isActive() ) {
02018                vertex->getRelations( neighs );
02019                for (size_t j = 0; j < neighs.size(); j++)
02020                     if (vertex < neighs[j]) {
02021                          Edge *newedge = new Edge( vertex, neighs[j]);
02022                          assert(newedge);
02023                          edges.push_back(newedge);
02024                     }
02025           }
02026      }
02027 
02028      if (!relexist0) clear_relations(0, 0);
02029 
02030      adjTable[1][0] = 1;
02031 }
02032 
02034 
02035 void
02036 Mesh::prune()
02037 {
02038      size_t nSize;
02039 
02040      NodeSequence vneighs;
02041 
02042      nSize = faces.size();
02043      for (size_t i = 0; i < nSize; i++) {
02044           Face *f = faces[i];
02045           if (!f->isActive()) {
02046                int nnodes = f->getSize(0);
02047                for( int j = 0; j < nnodes; j++) {
02048                     Vertex *v = f->getNodeAt(j);
02049                     v->setStatus( MeshEntity::REMOVE );
02050                }
02051           }
02052      }
02053 
02054      nSize = faces.size();
02055      for (size_t i = 0; i < nSize; i++) {
02056           Face *f = faces[i];
02057           if (f->isActive()) {
02058                int nnodes = f->getSize(0);
02059                for( int j = 0; j < nnodes; j++) {
02060                     Vertex *v = f->getNodeAt(j);
02061                     v->setStatus( MeshEntity::ACTIVE );
02062                }
02063           }
02064      }
02065 
02066      if (adjTable[0][0]) {
02067           nSize = nodes.size();
02068           for (size_t i = 0; i < nSize; i++) {
02069                Vertex *vi = nodes[i];
02070                if( !vi->isActive() ) {
02071                     vi->getRelations( vneighs );
02072                     for (size_t j = 0; j < vneighs.size(); j++)
02073                          vneighs[j]->removeRelation(vi);
02074                     vi->clearRelations(0);
02075                }
02076           }
02077      }
02078 
02079      if (adjTable[0][2]) {
02080           nSize = faces.size();
02081           for (size_t i = 0; i < nSize; i++) {
02082                Face *f = faces[i];
02083                if( !f->isActive() ) {
02084                     for (int j = 0; j < f->getSize(0); j++) {
02085                          Vertex *v = f->getNodeAt(j);
02086                          v->removeRelation(f);
02087                     }
02088                }
02089           }
02090      }
02091 
02092      NodeSequence::iterator vend;
02093      vend = remove_if(nodes.begin(), nodes.end(), EntityRemovedPred());
02094      nodes.erase(vend, nodes.end());
02095 
02096      nSize = nodes.size();
02097      for (size_t i = 0; i < nSize; i++)
02098           assert(nodes[i]->isActive());
02099 
02100      EdgeSequence::iterator eend;
02101      eend = remove_if(edges.begin(), edges.end(), EntityRemovedPred());
02102      edges.erase(eend, edges.end());
02103 
02104      nSize = edges.size();
02105      for(size_t i = 0; i < nSize; i++)
02106           assert(edges[i]->isActive());
02107 
02108      FaceSequence::iterator fend;
02109      fend = remove_if(faces.begin(), faces.end(), EntityRemovedPred());
02110      faces.erase(fend, faces.end());
02111 
02112      nSize = faces.size();
02113      for (size_t i = 0; i < nSize; i++)
02114           assert(faces[i]->isActive());
02115 
02116      enumerate(0);
02117      enumerate(2);
02118 }
02119 
02121 
02122 bool
02123 Mesh::isPruned() const
02124 {
02125      for (size_t i = 0; i < nodes.size(); i++)
02126           if (!nodes[i]->isActive()) return 0;
02127 
02128      for (size_t i = 0; i < faces.size(); i++)
02129           if (!faces[i]->isActive()) return 0;
02130 
02131      return 1;
02132 }
02133 
02135 
02136 void
02137 Mesh::collect_garbage()
02138 {
02139      if( !isPruned() ) prune();
02140 
02141      list<Face*>::const_iterator fiter;
02142      for (fiter = garbageFaces.begin(); fiter != garbageFaces.end(); ++fiter) {
02143           Face *face = *fiter;
02144           assert(face);
02145           if (face->isRemoved()) delete face;
02146      }
02147      garbageFaces.clear();
02148 
02149      if( nodes.empty() ) {
02150           list<Vertex*>::const_iterator viter;
02151           for (viter = garbageNodes.begin(); viter != garbageNodes.end(); ++viter) {
02152                Vertex *vertex = *viter;
02153                assert(vertex);
02154                if (vertex->isRemoved()) delete vertex;
02155           }
02156           garbageNodes.clear();
02157      }
02158 }
02159 
02161 
02162 void
02163 Mesh::enumerate(int etype)
02164 {
02165      size_t index = 0;
02166 
02167      NodeSequence::const_iterator viter;
02168      if (etype == 0) {
02169           for (viter = nodes.begin(); viter != nodes.end(); ++viter) {
02170                Vertex *vertex = *viter;
02171                if( vertex->isActive() ) vertex->setID(index++);
02172           }
02173      }
02174 
02175      FaceSequence::const_iterator fiter;
02176      if (etype == 2) {
02177           for (fiter = faces.begin(); fiter != faces.end(); ++fiter) {
02178                Face *face = *fiter;
02179                if( face->isActive() ) {
02180                     face->setID(index++);
02181                     for( int i  = 0; i < face->getSize(0); i++) {
02182                          if( !face->getNodeAt(i)->isActive() )
02183                               cout << "Error: Face is active, but one of the nodes is not" << endl;
02184                     }
02185                }
02186           }
02187      }
02188 }
02189 
02191 
02192 size_t
02193 Mesh::getBoundarySize(int d)
02194 {
02195      if (boundary_known == 0) search_boundary();
02196 
02197      size_t nSize, ncount = 0;
02198 
02199      if (d == 0) {
02200           nSize = nodes.size();
02201           for (size_t i = 0; i < nSize; i++)
02202                if (nodes[i]->isActive() && nodes[i]->isBoundary())
02203                     ncount++;
02204      }
02205 
02206      if (d == 2) {
02207           nSize = faces.size();
02208           for (size_t i = 0; i < nSize; i++)
02209                if (faces[i]->isActive() && faces[i]->isBoundary())
02210                     ncount++;
02211      }
02212 
02213      return ncount;
02214 }
02215 
02217 
02218 int Mesh :: getEulerCharacteristic()
02219 {
02220      size_t nSize;
02221 
02222      size_t  V = 0;
02223      size_t  E = 0;
02224      size_t  F = 0;
02225 
02226      nSize = getSize(0);
02227      for( size_t i = 0; i < nSize; i++)
02228           if( nodes[i]->isActive() ) V++;
02229 
02230      nSize = getSize(2);
02231      for( size_t i = 0; i < nSize; i++)
02232           if( faces[i]->isActive() ) F++;
02233 
02234      E = count_edges();
02235 
02236      return F - E + V;
02237 }
02239 
02240 
02241 int
02242 Mesh::isHomogeneous() const
02243 {
02244      int maxnodes = 0;
02245      size_t nSize = faces.size();
02246      for (size_t i = 0; i < nSize; i++) {
02247           if (faces[i]->isActive())
02248                maxnodes = max(maxnodes, faces[i]->getSize(0));
02249      }
02250 
02251      return maxnodes;
02252 }
02253 
02255 
02256 int
02257 Mesh::saveAs(const string &s)
02258 {
02259      MeshExporter mexp;
02260      return mexp.saveAs(this, s);
02261 }
02263 
02264 Vertex*
02265 Mesh::nearest_neighbour(const Vertex *myself, double &mindist)
02266 {
02267      Vertex* nearest;
02268      assert(getAdjTable(0, 0));
02269 
02270      mindist = std::numeric_limits< double >::max();
02271      nearest = NULL;
02272 
02273      NodeSequence neighs;
02274      myself->getRelations( neighs );
02275 
02276      for (size_t i = 0; i < neighs.size(); i++) {
02277           double d = Vertex::length2(myself, neighs[i]);
02278           if (d < mindist) {
02279                mindist = d;
02280                nearest = neighs[i];
02281           }
02282      }
02283 
02284      mindist = sqrt(mindist);
02285 
02286      return nearest;
02287 }
02288 
02290 int
02291 Mesh::build_relations01(bool rebuild)
02292 {
02293      if (rebuild) clear_relations(0, 1);
02294 
02295      int status = adjTable[0][1];
02296 
02297      if( status == 1 ) return 1;
02298 
02299      size_t numedges = getSize(1);
02300      for (size_t i = 0; i < numedges; i++) {
02301           Edge *edge = getEdgeAt(i);
02302           if( edge->isActive() ) {
02303                Vertex *v0 = edge->getNodeAt(0);
02304                Vertex *v1 = edge->getNodeAt(1);
02305                v0->addRelation(edge);
02306                v1->addRelation(edge);
02307           }
02308      }
02309 
02310      adjTable[0][1] = 1;
02311 
02312      return status;
02313 }
02314 
02315 int
02316 Mesh::build_relations02(bool rebuild)
02317 {
02318      if (rebuild) clear_relations(0, 2);
02319 
02320      int status = adjTable[0][2];
02321 
02322      if( status == 1 ) return 1;
02323 
02324      size_t numfaces = getSize(2);
02325      for (size_t iface = 0; iface < numfaces; iface++) {
02326           Face *face = getFaceAt(iface);
02327           if( face ) {
02328                if( face->isActive() ) {
02329                     int nf = face->getSize(0);
02330                     for (int j = 0; j < nf; j++) {
02331                          Vertex *vtx = face->getNodeAt(j);
02332                          vtx->addRelation(face);
02333                     }
02334                }
02335           }
02336      }
02337 
02338      adjTable[0][2] = 1;
02339 
02340      return status;
02341 }
02342 
02344 
02345 int
02346 Mesh::build_relations00(bool rebuild)
02347 {
02348      if (rebuild) clear_relations(0, 0);
02349 
02350      int status = adjTable[0][0];
02351      if( status == 1 ) return 1;
02352 
02353      size_t numfaces = getSize(2);
02354      for (size_t iface = 0; iface < numfaces; iface++) {
02355           Face *face = getFaceAt(iface);
02356           if( face ) {
02357                if( face->isActive() ) {
02358                     size_t nnodes = face->getSize(0);
02359                     for (size_t j = 0; j < nnodes; j++) {
02360                          Vertex *v0 = face->getNodeAt(j);
02361                          Vertex *v1 = face->getNodeAt(j + 1);
02362                          v0->addRelation(v1);
02363                          v1->addRelation(v0);
02364                     }
02365                }
02366           }
02367      }
02368      adjTable[0][0] = 1;
02369 
02370      return status;
02371 }
02372 
02374 
02375 int
02376 Mesh::build_relations12(bool rebuild)
02377 {
02378      if (rebuild) clear_relations(1, 2);
02379      int relexist2 = build_relations(0,2);
02380 
02381      int status = adjTable[1][2];
02382      if( status == 1 ) return 1;
02383 
02384      FaceSequence neighs;
02385      size_t numfaces = getSize(1);
02386      for (size_t iedge = 0; iedge < numfaces; iedge++) {
02387           Edge *edge = getEdgeAt(iedge);
02388           if( edge->isActive() ) {
02389                Vertex *v0 = edge->getNodeAt(0);
02390                Vertex *v1 = edge->getNodeAt(1);
02391                Mesh::getRelations112(v0, v1, neighs);
02392                for( size_t j = 0; j < neighs.size(); j++)
02393                     edge->addRelation( neighs[j] );
02394           }
02395      }
02396      adjTable[1][2] = 1;
02397 
02398      if( !relexist2 )  clear_relations(0,2);
02399 
02400      return status;
02401 }
02402 
02403 
02405 
02406 void
02407 Mesh::clear_relations(int src, int dst)
02408 {
02409      size_t numnodes = getSize(0);
02410 
02411      if (src == 0 && dst == 0) {
02412           for (size_t i = 0; i < numnodes; i++) {
02413                Vertex *vtx = getNodeAt(i);
02414                vtx->clearRelations(0);
02415           }
02416           adjTable[0][0] = 0;
02417      }
02418 
02419      if (src == 0 && dst == 2) {
02420           for (size_t i = 0; i < numnodes; i++) {
02421                Vertex *vtx = getNodeAt(i);
02422                vtx->clearRelations(2);
02423           }
02424           adjTable[0][2] = 0;
02425      }
02426 }
02427 
02429 
02430 int
02431 Mesh::search_boundary()
02432 {
02433      if (boundary_known == 1) return 1;
02434 
02435      if (!isPruned()) prune();
02436 
02437      int relexist = build_relations(0, 2);
02438      size_t numfaces = getSize(2);
02439 
02440      int bmark;
02441      FaceSequence neighs;
02442      for (size_t iface = 0; iface < numfaces; iface++) {
02443           Face *face = getFaceAt(iface);
02444           if( face->isActive() ) {
02445                size_t nnodes = face->getSize(0);
02446                for (size_t j = 0; j < nnodes; j++) {
02447                     Vertex *v0 = face->getNodeAt(j);
02448                     Vertex *v1 = face->getNodeAt(j + 1);
02449                     Mesh::getRelations112(v0, v1, neighs);
02450 
02451                     if (neighs.size() == 1) {
02452                          bmark = face->getBoundaryMark();
02453                          bmark = max(1, bmark);
02454                          face->setBoundaryMark(bmark);
02455 
02456                          bmark = v0->getBoundaryMark();
02457                          bmark = max(1, bmark);
02458                          v0->setBoundaryMark(bmark);
02459 
02460                          bmark = v1->getBoundaryMark();
02461                          bmark = max(1, bmark);
02462                          v1->setBoundaryMark(bmark);
02463                     }
02464                }
02465           }
02466      }
02467 
02468      if (!relexist)
02469           clear_relations(0, 2);
02470 
02471      boundary_known = 1;
02472 
02473      return 0;
02474 
02475 }
02477 
02478 Jaal::FaceSequence
02479 Mesh::filter(int facetype) const
02480 {
02481      FaceSequence::const_iterator it;
02482      size_t ncount = 0;
02483      for (it = faces.begin(); it != faces.end(); ++it) {
02484           Face *face = *it;
02485           if (face->getType() == facetype)
02486                ncount++;
02487      }
02488 
02489      FaceSequence tmpfaces;
02490      if (ncount) {
02491           tmpfaces.resize(ncount);
02492           size_t index = 0;
02493           for (it = faces.begin(); it != faces.end(); ++it) {
02494                Face *face = *it;
02495                if (face->getType() == facetype)
02496                     tmpfaces[index++] = face;
02497           }
02498      }
02499 
02500      return tmpfaces;
02501 }
02502 
02504 
02505 bool
02506 Mesh::isSimple()
02507 {
02508      int simple = 1;
02509      int relexist = build_relations(0, 2);
02510 
02511      size_t numfaces = getSize(2);
02512 
02513      FaceSequence neighs;
02514      for (size_t iface = 0; iface < numfaces; iface++) {
02515           Face *face = getFaceAt(iface);
02516           if( face->isActive() ) {
02517                assert(face);
02518                int nnodes = face->getSize(0);
02519                for (int j = 0; j < nnodes; j++) {
02520                     Vertex *v0 = face->getNodeAt(j);
02521                     Vertex *v1 = face->getNodeAt(j + 1);
02522                     Mesh::getRelations112(v0, v1, neighs);
02523                     if (neighs.size() > 2) {
02524                          simple = 0;
02525                          break;
02526                     }
02527                }
02528           }
02529      }
02530 
02531      if (!relexist)
02532           clear_relations(0, 2);
02533 
02534      return simple;
02535 
02536 }
02537 
02539 
02540 size_t Mesh::count_irregular_nodes(int degree_of_regular_node)
02541 {
02542      int relexist = build_relations(0, 2);
02543 
02544      search_boundary();
02545 
02546      size_t numnodes = getSize(0);
02547      size_t ncount = 0;
02548      for (size_t i = 0; i < numnodes; i++) {
02549           Vertex *vertex = getNodeAt(i);
02550           if( vertex->isActive() ) {
02551                int nSize = vertex->getNumRelations(2);
02552                if (!vertex->isBoundary()) {
02553                     if ( nSize != degree_of_regular_node)
02554                          ncount++;
02555                }
02556           }
02557      }
02558 
02559      if (!relexist)
02560           clear_relations(0, 2);
02561 
02562      return ncount;
02563 }
02565 
02566 bool
02567 Mesh::is_consistently_oriented()
02568 {
02569      int consistent = 1;
02570      int relexist = build_relations(0, 2);
02571 
02572      size_t numfaces = getSize(2);
02573 
02574      FaceSequence neighs;
02575      for (size_t iface = 0; iface < numfaces; iface++) {
02576           Face *face = getFaceAt(iface);
02577           if( face->isActive() ) {
02578                int nnodes = face->getSize(0);
02579                for (int j = 0; j < nnodes; j++) {
02580                     Vertex *v0 = face->getNodeAt(j);
02581                     Vertex *v1 = face->getNodeAt(j + 1);
02582                     Mesh::getRelations112(v0, v1, neighs);
02583                     if (neighs.size() == 2) {
02584                          assert(neighs[0]->isActive());
02585                          assert(neighs[1]->isActive());
02586                          int dir1 = neighs[0]->getOrientation(v0, v1);
02587                          int dir2 = neighs[1]->getOrientation(v0, v1);
02588                          if (dir1 * dir2 == 1) {
02589                               cout << "Warning: Mesh is not consistently oriented " << endl;
02590                               cout << neighs[0]->getID() << " Face 1: ";
02591                               for (int k = 0; k < neighs[0]->getSize(0); k++)
02592                                    cout << neighs[0]->getNodeAt(k)->getID() << " ";
02593                               cout << endl;
02594                               cout << neighs[1]->getID() << " Face 2: ";
02595                               for (int k = 0; k < neighs[1]->getSize(0); k++)
02596                                    cout << neighs[1]->getNodeAt(k)->getID() << " ";
02597                               cout << endl;
02598                               consistent = 0;
02599                               break;
02600                          }
02601                     }
02602                }
02603                if (!consistent)
02604                     break;
02605           }
02606      }
02607 
02608      if (!relexist)
02609           clear_relations(0, 2);
02610 
02611      return consistent;
02612 }
02613 
02615 
02616 void
02617 Mesh::make_consistently_oriented()
02618 {
02619      int relexist2 = build_relations(0, 2);
02620 
02621      Face *face = NULL;
02622      deque<Face*> faceQ;
02623      FaceSequence neighs;
02624 
02625      size_t numfaces = getSize(2);
02626 
02627      for (size_t iface = 0; iface < numfaces; iface++) {
02628           face = getFaceAt(iface);
02629           face->setID(iface);
02630           face->setVisitMark(0);
02631      }
02632 
02633      for (size_t iface = 0; iface < numfaces; iface++) {
02634           face = getFaceAt( iface );
02635           if( face->isActive() ) break;
02636      }
02637 
02638      faceQ.push_back(face);
02639 
02640      while (!faceQ.empty()) {
02641           Face *face = faceQ.front();
02642           faceQ.pop_front();
02643           if (!face->isVisited()) {
02644                face->setVisitMark(1);
02645                int nnodes = face->getSize(0);
02646                for (int j = 0; j < nnodes; j++) {
02647                     Vertex *v0 = face->getNodeAt(j);
02648                     Vertex *v1 = face->getNodeAt(j + 1);
02649                     Mesh::getRelations112(v0, v1, neighs);
02650                     if (neighs.size() == 2) {
02651                          int dir1 = neighs[0]->getOrientation(v0, v1);
02652                          int dir2 = neighs[1]->getOrientation(v0, v1);
02653                          if (dir1 * dir2 == 1) {
02654                               if (!neighs[0]->isVisited() && neighs[1]->isVisited())
02655                                    neighs[0]->reverse();
02656 
02657                               if (!neighs[1]->isVisited() && neighs[0]->isVisited())
02658                                    neighs[1]->reverse();
02659                          }
02660                          faceQ.push_back(neighs[0]);
02661                          faceQ.push_back(neighs[1]);
02662                     }
02663                }
02664           }
02665      }
02666 
02667      for (size_t iface = 0; iface < numfaces; iface++) {
02668           face = getFaceAt(iface);
02669           if( face->isActive() && !face->isVisited()) {
02670                cout << "Error: not visited : " << face->getID() << " "
02671                     << face->isVisited() << endl;
02672           }
02673      }
02674 
02675      if( !relexist2) clear_relations(0, 2);
02676 }
02677 
02679 
02680 int
02681 Mesh::getNumComponents(bool stop_at_interface)
02682 {
02683      build_relations(0, 2);
02684 
02685      Face *face = NULL;
02686      deque<Face*> faceQ;
02687      FaceSequence neighs;
02688 
02689      size_t numfaces = getSize(2);
02690 
02691      for (size_t iface = 0; iface < numfaces; iface++) {
02692           face = getFaceAt(iface);
02693           face->setID(iface);
02694           face->setVisitMark(0);
02695      }
02696 
02697      int numComponents = 0;
02698 
02699      while (1) {
02700           face = NULL;
02701           faceQ.clear();
02702           for (size_t iface = 0; iface < numfaces; iface++) {
02703                face = getFaceAt(iface);
02704                if (!face->isVisited()) {
02705                     face->setAttribute("Component", numComponents);
02706                     faceQ.push_back(face);
02707                     break;
02708                }
02709           }
02710 
02711           if (faceQ.empty())
02712                break;
02713 
02714           while (!faceQ.empty()) {
02715                Face *face = faceQ.front();
02716                faceQ.pop_front();
02717                if (!face->isVisited()) {
02718                     face->setAttribute( "Component", numComponents);
02719                     face->setVisitMark(1);
02720                     int nnodes = face->getSize(0);
02721                     for (int j = 0; j < nnodes; j++) {
02722                          Vertex *v0 = face->getNodeAt(j);
02723                          Vertex *v1 = face->getNodeAt(j + 1);
02724 
02725                          int proceed = 1;
02726                          if (stop_at_interface) {
02727                               if (v0->isConstrained() && v1->isConstrained()) {
02728                                    Edge edge(v0, v1);
02729                                    if (hasFeatureEdge(edge)) proceed = 0;
02730                               }
02731                          }
02732 
02733                          if (proceed) {
02734                               Mesh::getRelations112(v0, v1, neighs);
02735                               if (neighs.size() == 2) {
02736                                    faceQ.push_back(neighs[0]);
02737                                    faceQ.push_back(neighs[1]);
02738                               }
02739                          }
02740 
02741                     }
02742                }
02743           } // Complete one Component
02744           numComponents++;
02745      }
02746 
02747      for (size_t iface = 0; iface < numfaces; iface++) {
02748           face = getFaceAt(iface);
02749           if (!face->isVisited())
02750                cout << "Error: not visited : " << face->getID() << " "
02751                     << face->isVisited() << endl;
02752      }
02753 
02754      clear_relations(0, 2);
02755      return numComponents;
02756 }
02757 
02759 
02760 int Mesh::getNodes( NodeList &nodelist) const
02761 {
02762      nodelist.clear();
02763      size_t nSize = nodes.size();
02764      for (size_t i = 0; i < nSize; i++)  {
02765           Vertex *v = getNodeAt(i);
02766           if( v->isActive() ) nodelist.push_back( v );
02767      }
02768      return 0;
02769 }
02770 
02771 int Mesh::getFaces( FaceList &facelist) const
02772 {
02773      facelist.clear();
02774      size_t nSize = faces.size();
02775      for (size_t i = 0; i < nSize; i++) {
02776           Face *f = getFaceAt(i);
02777           if( f->isActive() ) facelist.push_back(f);
02778      }
02779      return 0;
02780 }
02781 
02782 Mesh* Mesh::getComponent(int id)
02783 {
02784      Mesh *submesh = new Mesh;
02785 
02786      size_t numnodes = getSize(0);
02787      for (size_t i = 0; i < numnodes; i++) {
02788           Vertex *v = getNodeAt(i);
02789           v->setVisitMark(0);
02790      }
02791 
02792      size_t numfaces = getSize(2);
02793 
02794      int pid;
02795      for (size_t iface = 0; iface < numfaces; iface++) {
02796           Face *face = getFaceAt(iface);
02797           if( face->isActive() ) {
02798                face->getAttribute("Partition", pid);
02799                if ( pid  == id) {
02800                     for (int j = 0; j < face->getSize(0); j++) {
02801                          Vertex *v = face->getNodeAt(j);
02802                          if (!v->isVisited()) {
02803                               submesh->addNode(v);
02804                               v->setVisitMark(1);
02805                          }
02806                     }
02807                     submesh->addFace(face);
02808                }
02809           }
02810      }
02811      return submesh;
02812 }
02814 
02815 Mesh *
02816 Jaal::struct_tri_grid(int nx, int ny)
02817 {
02818      Mesh *trimesh = new Mesh;
02819 
02820      double dx = 2.0 / (nx - 1);
02821      double dy = 2.0 / (ny - 1);
02822 
02823      Point3D xyz;
02824 
02825      int index = 0;
02826      for (int j = 0; j < ny; j++) {
02827           for (int i = 0; i < nx; i++) {
02828                xyz[0] = -1.0 + i * dx;
02829                xyz[1] = -1.0 + j * dy;
02830                xyz[2] = 0.0;
02831                Vertex *vnew = Vertex::newObject();
02832                vnew->setID(index++);
02833                vnew->setXYZCoords(xyz);
02834                trimesh->addNode(vnew);
02835           }
02836      }
02837 
02838      NodeSequence connect(3);
02839      index = 0;
02840      Face *newtri;
02841      for (int j = 0; j < ny - 1; j++) {
02842           for (int i = 0; i < nx - 1; i++) {
02843                int n0 = j * nx + i;
02844                int n1 = n0 + 1;
02845                int n2 = n1 + nx;
02846                int n3 = n0 + nx;
02847                connect[0] = trimesh->getNodeAt(n0);
02848                connect[1] = trimesh->getNodeAt(n1);
02849                connect[2] = trimesh->getNodeAt(n2);
02850                newtri = Face::newObject();
02851                newtri->setNodes(connect);
02852                trimesh->addFace(newtri);
02853 
02854                connect[0] = trimesh->getNodeAt(n0);
02855                connect[1] = trimesh->getNodeAt(n2);
02856                connect[2] = trimesh->getNodeAt(n3);
02857                newtri = Face::newObject();
02858                newtri->setNodes(connect);
02859                trimesh->addFace(newtri);
02860           }
02861      }
02862      return trimesh;
02863 }
02864 
02866 
02867 Mesh *
02868 Jaal::struct_quad_grid(int nx, int ny)
02869 {
02870      Mesh *quadmesh = new Mesh;
02871 
02872      double dx = 2.0 / (nx - 1);
02873      double dy = 2.0 / (ny - 1);
02874 
02875      Point3D xyz;
02876 
02877      int index = 0;
02878      for (int j = 0; j < ny; j++) {
02879           for (int i = 0; i < nx; i++) {
02880                xyz[0] = -1.0 + i * dx;
02881                xyz[1] = -1.0 + j * dy;
02882                xyz[2] = 0.0;
02883                Vertex *vnew = Vertex::newObject();
02884                vnew->setID(index++);
02885                vnew->setXYZCoords(xyz);
02886                quadmesh->addNode(vnew);
02887           }
02888      }
02889 
02890      NodeSequence connect(4);
02891      index = 0;
02892      Face *newquad;
02893      for (int j = 0; j < ny - 1; j++) {
02894           for (int i = 0; i < nx - 1; i++) {
02895                int n0 = j * nx + i;
02896                int n1 = n0 + 1;
02897                int n2 = n1 + nx;
02898                int n3 = n0 + nx;
02899                connect[0] = quadmesh->getNodeAt(n0);
02900                connect[1] = quadmesh->getNodeAt(n1);
02901                connect[2] = quadmesh->getNodeAt(n2);
02902                connect[3] = quadmesh->getNodeAt(n3);
02903                newquad = Face::newObject();
02904                newquad->setNodes(connect);
02905                quadmesh->addFace(newquad);
02906           }
02907      }
02908      return quadmesh;
02909 }
02910 
02912 
02913 void
02914 expand_strip(Face *prevface, Vertex *v0, Vertex *v1, list<Face*> &strip)
02915 {
02916      FaceSequence neighs;
02917      Mesh::getRelations112(v0, v1, neighs);
02918 
02919      Vertex *vn0, *vn1; // Next-Edge Nodes;
02920 
02921      Face *nextface;
02922      if (neighs.size() == 2) {
02923           if (!neighs[0]->isVisited() && neighs[1]->isVisited()) {
02924                nextface = neighs[0];
02925                if (nextface->getSize(0) == 4) {
02926                     Face::opposite_nodes(nextface, v0, v1, vn0, vn1);
02927                     nextface->setVisitMark(1);
02928                     strip.push_back(nextface);
02929                     expand_strip(nextface, vn0, vn1, strip);
02930                }
02931           }
02932           if (neighs[0]->isVisited() && !neighs[1]->isVisited()) {
02933                nextface = neighs[1];
02934                if (nextface->getSize(0) == 4) {
02935                     Face::opposite_nodes(nextface, v0, v1, vn0, vn1);
02936                     nextface->setVisitMark(1);
02937                     strip.push_back(nextface);
02938                     expand_strip(nextface, vn0, vn1, strip);
02939                }
02940           }
02941      }
02942 }
02944 
02945 #ifdef CSV
02946 void
02947 Mesh::get_quad_strips(Face *rootface, FaceSequence &strip1,
02948                       FaceSequence &strip2)
02949 {
02950 
02951      Vertex *v0, *v1;
02952 
02953      list<Face*> strip01, strip12, strip23, strip03;
02954      list<Face*>::const_iterator it;
02955      size_t numfaces = getSize(2);
02956 
02957      // Strip Starting from edge 0-1
02958      for (size_t i = 0; i < numfaces; i++) {
02959           Face *face = getFaceAt(i);
02960           face->setVisitMark(0);
02961      }
02962 
02963      v0 = rootface->getNodeAt(0);
02964      v1 = rootface->getNodeAt(1);
02965 
02966      rootface->setVisitMark(1);
02967      strip01.push_back(rootface);
02968      expand_strip(rootface, v0, v1, strip01);
02969      for (it = strip01.begin(); it != strip01.end(); ++it)
02970           strip1.push_back(*it);
02971 
02972      // Strip Starting from edge 2-3
02973      for (size_t i = 0; i < numfaces; i++) {
02974           Face *face = getFaceAt(i);
02975           face->setVisitMark(0);
02976      }
02977 
02978      v0 = rootface->getNodeAt(2);
02979      v1 = rootface->getNodeAt(3);
02980 
02981      rootface->setVisitMark(1);
02982      strip23.push_back(rootface);
02983      expand_strip(rootface, v0, v1, strip23);
02984      for (it = strip23.begin(); it != strip23.end(); ++it)
02985           strip1.push_back(*it);
02986 
02987      // Strip Starting from edge 1-2
02988      for (size_t i = 0; i < numfaces; i++) {
02989           Face *face = getFaceAt(i);
02990           face->setVisitMark(0);
02991      }
02992 
02993      v0 = rootface->getNodeAt(1);
02994      v1 = rootface->getNodeAt(2);
02995 
02996      rootface->setVisitMark(1);
02997      strip12.push_back(rootface);
02998      expand_strip(rootface, v0, v1, strip12);
02999      for (it = strip12.begin(); it != strip12.end(); ++it)
03000           strip2.push_back(*it);
03001 
03002      // Strip Starting from edge 0-3
03003      for (size_t i = 0; i < numfaces; i++) {
03004           Face *face = getFaceAt(i);
03005           face->setVisitMark(0);
03006      }
03007 
03008      v0 = rootface->getNodeAt(0);
03009      v1 = rootface->getNodeAt(3);
03010 
03011      rootface->setVisitMark(1);
03012      strip12.push_back(rootface);
03013      expand_strip(rootface, v0, v1, strip03);
03014 
03015      for (it = strip03.begin(); it != strip03.end(); ++it)
03016           strip2.push_back(*it);
03017 }
03018 #endif
03019 
03020 
03021 int
03022 Mesh::get_bound_faces( FaceSequence &result, int bound_what)
03023 {
03024      result.clear();
03025 
03026      int relexist2 = build_relations(0, 2);
03027 
03028      assert(getAdjTable(0, 2));
03029 
03030      search_boundary();
03031 
03032      size_t numfaces = getSize(2);
03033 
03034      set<Face*> bfaces;
03035 
03036      if (bound_what == 0) {
03037           for (size_t i = 0; i < numfaces; i++) {
03038                Face *face = getFaceAt(i);
03039                if (face->hasBoundaryNode() && face->isActive() )
03040                     bfaces.insert(face);
03041           }
03042      }
03043 
03044      if (bound_what == 1) {
03045           for (size_t i = 0; i < numfaces; i++) {
03046                Face *face = getFaceAt(i);
03047                if (face->has_boundary_edge() && face->isActive() )
03048                     bfaces.insert(face);
03049           }
03050      }
03051 
03052      size_t nSize = bfaces.size();
03053 
03054      if (nSize) {
03055           result.resize(nSize);
03056           set<Face*>::const_iterator it;
03057 
03058           size_t index = 0;
03059           for (it = bfaces.begin(); it != bfaces.end(); ++it)
03060                result[index++] = *it;
03061      }
03062 
03063      if (!relexist2)
03064           clear_relations(0, 2);
03065 
03066      return 0;
03067 }
03068 
03070 
03071 int Mesh::get_bound_nodes( NodeSequence &seq)
03072 {
03073      search_boundary();
03074 
03075      size_t numnodes = getSize(0);
03076 
03077      seq.clear() ;
03078      for (size_t i = 0; i < numnodes; i++) {
03079           Vertex *v = getNodeAt(i);
03080           if (v->isBoundary() && v->isActive() ) seq.push_back(v);
03081      }
03082      return 0;
03083 }
03085 
03086 int Mesh::get_irregular_nodes(NodeSequence &seq, int regular_count, int from_where)
03087 {
03088      seq.clear();
03089 
03090      int relexist2 = build_relations(0, 2);
03091      search_boundary();
03092 
03093      size_t numnodes = getSize(0);
03094      int    nSize;
03095 
03096 
03097      // Query from the boundary nodes ...
03098      if (from_where == 1) {
03099           for (size_t i = 0; i < numnodes; i++) {
03100                Vertex *v = getNodeAt(i);
03101                if( v->isActive() ) {
03102                     nSize  = v->getNumRelations(2);
03103                     if (v->isBoundary() && nSize != regular_count) seq.push_back(v);
03104                }
03105           }
03106      }
03107 
03108      // Query from the internal nodes ...
03109      if (from_where == 0) {
03110           for (size_t i = 0; i < numnodes; i++) {
03111                Vertex *v = getNodeAt(i);
03112                if( v->isActive() ) {
03113                     nSize  = v->getNumRelations(2);
03114                     if (!v->isBoundary() &&  nSize != regular_count) seq.push_back(v);
03115                }
03116           }
03117      }
03118 
03119      if (!relexist2)
03120           clear_relations(0, 2);
03121 
03122      return 0;
03123 }
03124 
03126 
03127 #ifdef CSV
03128 void
03129 Mesh::set_strip_markers()
03130 {
03131      FaceSequence bound_faces;
03132      get_bound_faces(bound_faces, 1);
03133 
03134      size_t numfaces = getSize(2);
03135      for (size_t i = 0; i < numfaces; i++) {
03136           Face *face = getFaceAt(i);
03137           face->setAttribute("Strip", 0);
03138      }
03139 
03140      FaceSequence strip1, strip2;
03141      int id = 0;
03142      for (size_t i = 0; i < bound_faces.size(); i++) {
03143           strip1.clear();
03144           strip2.clear();
03145           get_quad_strips(bound_faces[i], strip1, strip2);
03146           id++;
03147           for (size_t j = 0; j < strip1.size(); j++)
03148                strip1[j]->setAttribute( "Strip", id);
03149           /*
03150            id++;
03151            for( int j = 0; j < strip2.size(); j++)
03152            strip2[j]->setTag( id );
03153            */
03154      }
03155 }
03156 #endif
03157 
03159 
03160 vector<int>
03161 Mesh::get_topological_statistics(int entity, bool sorted)
03162 {
03163      int relexist = build_relations(0, 2);
03164 
03165      assert(getAdjTable(0, 2));
03166 
03167      int numnodes = getSize(0);
03168 
03169      vector<int> degree;
03170      degree.reserve(numnodes);
03171      for (int i = 0; i < numnodes; i++) {
03172           Vertex *v = getNodeAt(i);
03173           if( !v->isRemoved() )
03174                degree.push_back(v->getNumRelations(2) );
03175      }
03176 
03177      int mindegree = *min_element(degree.begin(), degree.end());
03178      int maxdegree = *max_element(degree.begin(), degree.end());
03179 
03180      cout << " Mesh Topological Quality : " << endl;
03181 
03182      cout << " ************************ " << endl;
03183      cout << " Degree         FaceCount " << endl;
03184      cout << " ************************ " << endl;
03185 
03186      for (int i = mindegree; i <= maxdegree; i++) {
03187           int ncount = 0;
03188           for (size_t j = 0; j < degree.size(); j++)
03189                if (degree[j] == i)
03190                     ncount++;
03191           cout << setw(5) << i << setw(15) << ncount << endl;
03192      }
03193 
03194      if (!relexist)
03195           clear_relations(0, 2);
03196 
03197      return degree;
03198 }
03199 
03201 
03202 size_t
03203 Mesh::setNodeWavefront(int layerid)
03204 {
03205      NodeSequence vnodes;
03206      assert(layerid >= 0);
03207 
03208      int relexist = build_relations(0, 0);
03209 
03210      size_t numNodes = getSize(0);
03211 
03212      size_t ncount = 0;
03213      int lid = 0;
03214 
03215      if (layerid == 0) {
03216           search_boundary();
03217           for (size_t i = 0; i < numNodes; i++) {
03218                Vertex *vertex = getNodeAt(i);
03219                if (vertex->isBoundary()) {
03220                     vertex->setAttribute("Layer", 0);
03221                     ncount++;
03222                }
03223           }
03224      } else {
03225           for (size_t i = 0; i < numNodes; i++) {
03226                Vertex *vertex = getNodeAt(i);
03227                vertex->getAttribute("Layer", lid );
03228                if ( lid == layerid - 1) {
03229                     vertex->getRelations( vnodes );
03230                     for (size_t j = 0; j < vnodes.size(); j++) {
03231                          vnodes[j]->getAttribute("Layer", lid);
03232                          if (lid >= 0 && lid <= layerid - 1) continue;
03233                          vnodes[j]->setAttribute("LayerID", layerid );
03234                          ncount++;
03235                     }
03236                }
03237           }
03238      }
03239 
03240      if (!relexist)
03241           clear_relations(0, 0);
03242 
03243      return ncount;
03244 }
03245 
03247 
03248 int
03249 Mesh::setNodeWavefront()
03250 {
03251      int relexist = build_relations(0, 0);
03252      search_boundary();
03253 
03254      size_t numNodes = getSize(0);
03255      for (size_t i = 0; i < numNodes; i++) {
03256           Vertex *v = getNodeAt(i);
03257           v->setAttribute("Layer", -1);
03258           v->setVisitMark(0);
03259      }
03260 
03261      NodeSequence vertexQ;
03262      NodeSet  nextQ;
03263      for (size_t i = 0; i < numNodes; i++) {
03264           Vertex *v = getNodeAt(i);
03265           if (v->isBoundary() && !v->isRemoved() ) {
03266                v->setAttribute("Layer", 0);
03267                v->setVisitMark(1);
03268                vertexQ.push_back(v);
03269           }
03270      }
03271 
03272      if (vertexQ.empty())
03273           cout << "Warning: No boundary detected " << endl;
03274 
03275      NodeSequence vneighs;
03276      int layerid = 1;
03277      size_t nSize;
03278      while (!vertexQ.empty()) {
03279           nextQ.clear();
03280           nSize = vertexQ.size();
03281           for (size_t j = 0; j < nSize; j++) {
03282                Vertex *currVertex = vertexQ[j];
03283                currVertex->getRelations( vneighs );
03284                for (size_t i = 0; i < vneighs.size(); i++) {
03285                     Vertex *vn = vneighs[i];
03286                     if (!vn->isVisited()) nextQ.insert(vn);
03287                }
03288           }
03289 
03290           nSize = nextQ.size();
03291           if( nSize == 0) break;
03292 
03293           NodeSet::const_iterator it;
03294           for( it = nextQ.begin(); it != nextQ.end(); ++it) {
03295                Vertex *currVertex = *it;
03296                currVertex->setAttribute("Layer", layerid);
03297                currVertex->setVisitMark(1);
03298           }
03299 
03300           vertexQ.clear();
03301           vertexQ.reserve( nSize );
03302           for( it = nextQ.begin(); it != nextQ.end(); ++it)
03303                vertexQ.push_back( *it );
03304           layerid++;
03305      }
03306 
03307      if (!relexist)
03308           clear_relations(0, 0);
03309 
03310      cout << "Info: Number of layers in node front : " << layerid - 1 << endl;
03311 
03312      return layerid - 1;
03313 }
03314 
03316 
03317 size_t
03318 Mesh::setFaceWavefront(int layerid)
03319 {
03320      assert(layerid >= 0);
03321 
03322      int relexist = build_relations(0, 2);
03323 
03324      size_t numfaces = getSize(2);
03325 
03326      size_t ncount = 0;
03327      int lid = 0;
03328 
03329      FaceSequence vfaces;
03330      if (layerid == 0) {
03331           for (size_t i = 0; i < numfaces; i++) {
03332                Face *face = getFaceAt(i);
03333                int nsize = face->getSize(0);
03334                for (int j = 0; j < nsize; j++) {
03335                     Vertex *v0 = face->getNodeAt(j + 0);
03336                     Vertex *v1 = face->getNodeAt(j + 1);
03337                     Mesh::getRelations112(v0, v1, vfaces);
03338                     if (vfaces.size() == 1) {
03339                          face->setAttribute("Layer", 0);
03340                          ncount++;
03341                     }
03342                }
03343           }
03344      } else {
03345           for (size_t i = 0; i < numfaces; i++) {
03346                Face *face = getFaceAt(i);
03347                face->getAttribute("Layer", lid);
03348                if(lid == layerid - 1) {
03349                     face->getRelations12( vfaces );
03350                     for (size_t j = 0; j < vfaces.size(); j++) {
03351                          vfaces[j]->getAttribute("Layer", lid );
03352                          if (lid >= 0 && lid < layerid) continue;
03353                          vfaces[j]->setAttribute("Layer", layerid);
03354                          ncount++;
03355                     }
03356                }
03357           }
03358      }
03359 
03360      if (!relexist)
03361           clear_relations(0, 2);
03362 
03363      return ncount;
03364 }
03365 
03367 
03368 int
03369 Mesh::setFaceWavefront()
03370 {
03371      int relexist = build_relations(0, 2);
03372 
03373      search_boundary();
03374 
03375      size_t numFaces = getSize(2);
03376 
03377      for (size_t i = 0; i < numFaces; i++) {
03378           Face *f = getFaceAt(i);
03379           f->setAttribute("Layer", 0);
03380           f->setVisitMark(0);
03381      }
03382 
03383      FaceSequence faceQ, nextQ, neighs;
03384      for (size_t i = 0; i < numFaces; i++) {
03385           Face *f = getFaceAt(i);
03386           if( !f->isRemoved() ) {
03387                f->setAttribute("Layer", 1);
03388                if (f->has_boundary_edge()) {
03389                     f->setAttribute("Layer", 0);
03390                     f->setVisitMark(1);
03391                     faceQ.push_back(f);
03392                }
03393           }
03394      }
03395 
03396      int layerid = 1;
03397      size_t nSize;
03398 
03399      while (!faceQ.empty()) {
03400           nextQ.clear();
03401           nextQ.reserve(faceQ.size());
03402           nSize = faceQ.size();
03403           for (size_t j = 0; j < nSize; j++) {
03404                Face *currFace = faceQ[j];
03405                currFace->getRelations12( neighs );
03406                for (size_t i = 0; i < neighs.size(); i++) {
03407                     Face *nf = neighs[i];
03408                     if (!nf->isVisited()) nextQ.push_back(nf);
03409                }
03410           }
03411 
03412           nSize = nextQ.size();
03413           for (size_t i = 0; i < nSize; i++) {
03414                Face *f = nextQ[i];
03415                f->setAttribute("Layer", layerid);
03416                f->setVisitMark(1);
03417           }
03418 
03419           layerid++;
03420           faceQ = nextQ;
03421      }
03422 
03423      if (!relexist)
03424           clear_relations(0, 2);
03425 
03426      cout << "Info: Number of layers in face front : " << layerid - 1 << endl;
03427 
03428      return layerid - 1;
03429 }
03431 
03432 int
03433 Mesh::setWavefront(int forwhat)
03434 {
03435      if (forwhat == 0)
03436           return setNodeWavefront();
03437 
03438      if (forwhat == 2)
03439           return setFaceWavefront();
03440 
03441      return 0;
03442 }
03443 
03445 
03446 int
03447 Mesh::verify_front_ordering(int mentity)
03448 {
03449      int l1, l2;
03450      size_t numfaces = getSize(2);
03451      if (mentity == 0) {
03452           for (size_t i = 0; i < numfaces; i++) {
03453                Face *face = getFaceAt(i);
03454                int nsize = face->getSize(0);
03455                for (int j = 0; j < nsize; j++) {
03456                     Vertex *v0 = face->getNodeAt( j + 0 );
03457                     Vertex *v1 = face->getNodeAt( j + 1 );
03458                     v0->getAttribute("Layer", l1 );
03459                     v1->getAttribute("Layer", l2 );
03460                     assert(l1 >= 0 && l2 >= 0);
03461                     if (abs(l2 - l1) > 1) return 1;
03462                }
03463           }
03464      }
03465 
03466      FaceSequence neighs;
03467 
03468      if (mentity == 2) {
03469           int relexist2 = build_relations(0, 2);
03470           for (size_t i = 0; i < numfaces; i++) {
03471                Face *face = getFaceAt(i);
03472                face->getRelations12( neighs );
03473                face->getAttribute("Layer", l1);
03474                int nf  = neighs.size();
03475                assert(l1 >= 0);
03476                for (int j = 0; j < nf; j++) {
03477                     neighs[j]->getAttribute("Layer", l2);
03478                     assert(l2 >= 0);
03479                     if (abs(l2 - l1) > 1) return 1;
03480                }
03481           }
03482           if (!relexist2) clear_relations(0, 2);
03483      }
03484 
03485      return 0;
03486 }
03487 
03489 
03490 int Mesh::remove_unreferenced_nodes()
03491 {
03492      size_t numnodes = getSize(0);
03493      for (size_t i = 0; i < numnodes; i++) {
03494           Vertex *v = getNodeAt(i);
03495           v->setVisitMark(0);
03496      }
03497 
03498      size_t numfaces = getSize(2);
03499      for (size_t i = 0; i < numfaces; i++) {
03500           Face *face = getFaceAt(i);
03501           for (int j = 0; j < face->getSize(0); j++) {
03502                Vertex *v = face->getNodeAt(j);
03503                v->setVisitMark(1);
03504           }
03505      }
03506 
03507      for (size_t i = 0; i < numnodes; i++) {
03508           Vertex *v = getNodeAt(i);
03509           if (!v->isVisited()) cout << " Has Unreferenced Node " << endl;
03510      }
03511 
03512      /*
03513         for( size_t i = 0; i < numnodes; i++) {
03514              Vertex *v = getNodeAt(i);
03515              if( !v->isVisited()) v->setRemoveMark(1);
03516         }
03517       */
03518 
03519      return 0;
03520 }
03522 
03523 double
03524 Mesh::getSurfaceArea()
03525 {
03526      double facearea, sumArea = 0.0;
03527 
03528      size_t numfaces = getSize(2);
03529      double minarea = std::numeric_limits< double >::max();
03530      double maxarea = 0.0;
03531      for (size_t i = 0; i < numfaces; i++) {
03532           Face *face = getFaceAt(i);
03533           facearea = face->getArea();
03534           sumArea += facearea;
03535           if (facearea < minarea) minarea = facearea;
03536           if (facearea > maxarea) maxarea = facearea;
03537      }
03538 
03539      /*
03540         cout << "Info:   Min face Area : " << minarea << endl;
03541         cout << "Info:   Max face Area : " << maxarea << endl;
03542       */
03543      return sumArea;
03544 }
03546 
03547 NodeSequence
03548 Mesh::chain_nodes(const vector<Edge> &bndedges)
03549 {
03550      NodeSequence bndnodes;
03551      bndnodes.push_back(bndedges[0].getNodeAt(0));
03552      bndnodes.push_back(bndedges[0].getNodeAt(1));
03553 
03554      size_t nSize = bndedges.size();
03555      for (size_t i = 1; i < nSize - 1; i++) {
03556           Vertex *v0 = bndedges[i].getNodeAt(0);
03557           Vertex *v1 = bndedges[i].getNodeAt(1);
03558           if (v0 == bndnodes[i]) {
03559                bndnodes.push_back(v1);
03560           } else if (v1 == bndnodes[i])
03561                bndnodes.push_back(v0);
03562           else {
03563                cout << "Error in bound edges : " << endl;
03564                bndnodes.clear();
03565                return bndnodes;
03566           }
03567      }
03568      return bndnodes;
03569 }
03570 
03571 vector<double>
03572 Mesh::getAspectRatio(bool sorted)
03573 {
03574      size_t numfaces = getSize(2);
03575 
03576      vector<double> quality;
03577      quality.resize(numfaces);
03578 
03579      for (size_t i = 0; i < numfaces; i++) {
03580           Face *face = getFaceAt(i);
03581           quality[i] = face->getAspectRatio();
03582      }
03583 
03584      if (sorted)
03585           std::sort(quality.begin(), quality.end());
03586 
03587      double minval = *min_element(quality.begin(), quality.end());
03588      double maxval = *max_element(quality.begin(), quality.end());
03589 
03590      cout << "Info: Minimum Aspect Ratio  " << minval << endl;
03591      cout << "Info: Maximum Aspect Ratio  " << maxval << endl;
03592 
03593      return quality;
03594 }
03596 
03597 void
03598 Mesh::emptyAll()
03599 {
03600      nodes.clear();
03601      edges.clear();
03602      faces.clear();
03603 }
03604 
03606 
03607 void Mesh::deleteNodes()
03608 {
03609      size_t nSize = nodes.size();
03610      for (size_t i = 0; i < nSize; i++) {
03611           if( !nodes[i]->isRemoved() ) delete nodes[i];
03612      }
03613      nodes.clear();
03614 }
03615 
03616 void Mesh::deleteEdges()
03617 {
03618      size_t nSize = edges.size();
03619      for (size_t i = 0; i < nSize; i++) {
03620           if( !edges[i]->isRemoved() ) delete edges[i];
03621      }
03622      edges.clear();
03623 }
03624 
03625 void Mesh::deleteFaces()
03626 {
03627      size_t nSize = faces.size();
03628      for (size_t i = 0; i < nSize; i++) {
03629           if( !faces[i]->isRemoved() ) delete faces[i];
03630      }
03631      faces.clear();
03632 }
03633 
03634 void
03635 Mesh::deleteAll()
03636 {
03637      for (int i = 0; i < 4; i++)
03638           for (int j = 0; j < 4; j++)
03639                clear_relations(i, j);
03640 
03641      deleteFaces();
03642      deleteEdges();
03643      deleteNodes();
03644 
03645      boundary_known = 0;
03646 }
03647 
03649 
03650 int
03651 Face ::refine_quad14(NodeSequence &newnodes, FaceSequence &newfaces)
03652 {
03653      assert(newnodes.size() == 5);
03654 
03655      Point3D xyz;
03656      this->getAvgPos( xyz );
03657      newnodes[4] = Vertex::newObject();
03658      newnodes[4]->setXYZCoords(xyz);
03659 
03660      newfaces.resize(4);
03661      NodeSequence qC(4);
03662      for (int i = 0; i < 4; i++)
03663           newfaces[i] = Face::newObject();
03664 
03665      qC[0] = this->getNodeAt(0);
03666      qC[1] = newnodes[0];
03667      qC[2] = newnodes[4];
03668      qC[3] = newnodes[3];
03669      newfaces[0]->setNodes(qC);
03670 
03671      qC[0] = newnodes[0];
03672      qC[1] = this->getNodeAt(1);
03673      qC[2] = newnodes[1];
03674      qC[3] = newnodes[4];
03675      newfaces[1]->setNodes(qC);
03676 
03677      qC[0] = newnodes[1];
03678      qC[1] = this->getNodeAt(2);
03679      qC[2] = newnodes[2];
03680      qC[3] = newnodes[4];
03681      newfaces[2]->setNodes(qC);
03682 
03683      qC[0] = newnodes[3];
03684      qC[1] = newnodes[4];
03685      qC[2] = newnodes[2];
03686      qC[3] = this->getNodeAt(3);
03687      newfaces[3]->setNodes(qC);
03688      return 0;
03689 }
03691 
03692 int
03693 Mesh::refine_quads14()
03694 {
03695      build_edges();
03696 
03697      std::multimap<Vertex*, Edge*> edgemap;
03698      Point3D xyz;
03699      for (size_t i = 0; i < edges.size(); i++) {
03700           Vertex *v0 = edges[i]->getNodeAt(0);
03701           Vertex *v1 = edges[i]->getNodeAt(1);
03702 
03703           Vertex::mid_point(v0, v1, xyz);
03704           Vertex *vnew = Vertex::newObject();
03705           vnew->setXYZCoords(xyz);
03706 
03707           edges[i]->addRelation(vnew);
03708           Vertex *vhash = edges[i]->getHashNode();
03709           edgemap.insert(std::make_pair(vhash, edges[i]));
03710           this->addNode(vnew);
03711      }
03712 
03713      NodeSequence newnodes(5);
03714      FaceSequence newfaces(4);
03715      NodeSequence vneighs;
03716 
03717      multimap<Vertex*, Edge*>::const_iterator ilower, iupper, iter;
03718 
03719      size_t numfaces = this->getSize(2);
03720      for (size_t iface = 0; iface < numfaces; iface++) {
03721           Face *face = this->getFaceAt(iface);
03722           if( !face->isRemoved() ) {
03723                for (int j = 0; j < 4; j++) {
03724                     Vertex *v0 = face->getNodeAt((j + 0));
03725                     Vertex *v1 = face->getNodeAt((j + 1));
03726                     Edge edge(v0, v1);
03727                     Vertex *vh = edge.getHashNode();
03728                     ilower = edgemap.lower_bound(vh);
03729                     iupper = edgemap.upper_bound(vh);
03730                     Vertex *vmid = NULL;
03731                     for (iter = ilower; iter != iupper; ++iter) {
03732                          Edge *existing_edge = iter->second;
03733                          if (existing_edge->isSameAs(edge)) {
03734                               existing_edge->getRelations( vneighs );
03735                               assert(vneighs.size() == 1);
03736                               vmid = vneighs[0];
03737                               break;
03738                          }
03739                     }
03740                     assert(vmid);
03741                     newnodes[j] = vmid;
03742                }
03743                face->refine_quad14(newnodes, newfaces);
03744                this->addNode(newnodes[4]); // Only the center node, others already registered
03745                for (int j = 0; j < 4; j++)
03746                     this->addFace(newfaces[j]);
03747                face->setStatus(MeshEntity::REMOVE);
03748           }
03749      }
03750 
03751      prune();
03752      enumerate(0);
03753      enumerate(2);
03754      collect_garbage();
03755 
03756      return 0;
03757 }
03759 
03760 int
03761 Mesh::refine_quads15()
03762 {
03763      FaceSequence newfaces;
03764      NodeSequence newnodes;
03765 
03766      int nSize;
03767      size_t numfaces = getSize(2);
03768      for (size_t iface = 0; iface < numfaces; iface++) {
03769           Face *face = getFaceAt(iface);
03770           if( !face->isRemoved() ) {
03771                face->refine_quad15(newnodes, newfaces);
03772 
03773                nSize = newnodes.size();
03774                for (int  j = 0; j < nSize; j++)
03775                     addNode(newnodes[j]);
03776 
03777                nSize = newfaces.size();
03778                for (int j = 0; j < nSize; j++)
03779                     addFace(newfaces[j]);
03780                face->setStatus(MeshEntity::REMOVE);
03781           }
03782      }
03783 
03784      prune();
03785      enumerate(0);
03786      enumerate(2);
03787      collect_garbage();
03788 
03789      return 0;
03790 }
03792 
03793 /*
03794 int
03795 Mesh::get_quality_statistics(const string &fname)
03796 {
03797     ofstream ofile(fname.c_str(), ios::out);
03798     if (ofile.fail())
03799     {
03800         cout << "Warning: Cann't open file " << fname << endl;
03801         return 1;
03802     }
03803 
03804     // Collect Element Area informtion
03805     size_t numfaces = getSize(2);
03806     vector<double> quality(numfaces);
03807 
03808     double minval, maxval, stddev, avgval, medianval;
03809 
03810     for (size_t i = 0; i < numfaces; i++)
03811     {
03812         Face *face = getFaceAt(i);
03813         quality[i] = face->getArea();
03814     }
03815 
03816     sort(quality.begin(), quality.end());
03817     minval = quality.front();
03818     maxval = quality.back();
03819 
03820     ofile << "# ********************************************* " << endl;
03821     ofile << "# Measure            Element Area        " << endl;
03822     ofile << "# Num of Elements  " << numfaces << endl;
03823     ofile << "# Min              " << minval << endl;
03824     ofile << "# Max              " << maxval << endl;
03825     ofile << "# Average          " << avgval << endl;
03826     ofile << "# MeanVal          " << medianval << endl;
03827     //ofile <<  "# StdDev           " <<   stddev     << endl;
03828     for (size_t i = 0; i < numfaces; i++)
03829         ofile << quality[i] << endl;
03830     ofile << " ********************************************* " << endl;
03831 
03832 
03833     for (size_t i = 0; i < numfaces; i++)
03834     {
03835         Face *face = getFaceAt(i);
03836         quality[i] = face->getAspectRatio();
03837     }
03838 
03839     sort(quality.begin(), quality.end());
03840     minval = quality.front();
03841     maxval = quality.back();
03842 
03843     ofile << "# Measure            Aspect Ratio      " << endl;
03844     ofile << "# Num of Elements  " << numfaces << endl;
03845     ofile << "# Min              " << minval << endl;
03846     ofile << "# Max              " << maxval << endl;
03847     ofile << "# Average          " << avgval << endl;
03848     ofile << "# MeanVal          " << medianval << endl;
03849     ofile << "# StdDev           " << stddev << endl;
03850 
03851     for (size_t i = 0; i < numfaces; i++)
03852         ofile << quality[i] << endl;
03853     ofile << " ********************************************* " << endl;
03854 }
03855  */
03856 
03858 double Vertex::getFeatureLength() const
03859 {
03860      if (!isBoundary()) return std::numeric_limits< double >::max();
03861 
03862      NodeSequence vneighs;
03863      getRelations( vneighs );
03864 
03865      assert(!vneighs.empty());
03866 
03867      double minlen = std::numeric_limits< double >::max();
03868 
03869      int  nSize = vneighs.size();
03870      for (int j = 0; j < nSize; j++) {
03871           if (vneighs[j]->isBoundary()) {
03872                const Point3D &p0 = getXYZCoords();
03873                const Point3D &p1 = vneighs[j]->getXYZCoords();
03874                minlen = min(minlen, Math::length(p0, p1));
03875           }
03876      }
03877      return minlen;
03878 }
03879 
03881 
03882 /*
03883 void
03884 Mesh::setFeatureLength()
03885 {
03886     int relexist0 = build_relations(0, 0);
03887     int relexist2 = build_relations(0, 2);
03888 
03889     search_boundary();
03890 
03891     size_t numnodes = getSize(0);
03892     vector<Vertex*> vneighs, bndnodes, sidenodes;
03893 
03894     for (size_t i = 0; i < numnodes; i++)
03895     {
03896         Vertex *vertex = getNodeAt(i);
03897         if (vertex->isBoundary())
03898         {
03899             vneighs = vertex->getRelations0();
03900             double minlen = std::numeric_limits< double >::max();
03901             for (int j = 0; j < vneighs.size(); j++)
03902             {
03903                 if (vneighs[j]->isBoundary())
03904                 {
03905                     Point3D p0 = vertex->getXYZCoords();
03906                     Point3D p1 = vneighs[j]->getXYZCoords();
03907                     minlen = min(minlen, Math::length(p0, p1));
03908                 }
03909             }
03910             vertex->setFeatureLength(minlen);
03911         }
03912     }
03913 
03914     if (!relexist0)
03915         clear_relations(0, 0);
03916 
03917     if (!relexist2)
03918         clear_relations(0, 2);
03919 }
03920  */
03921 
03923 
03924 int
03925 Mesh::hasDuplicates(int what)
03926 {
03927      map<Vertex*, FaceSequence> mapfaces;
03928 
03929      size_t numfaces = getSize(2);
03930      assert(numfaces);
03931 
03932      for (size_t i = 0; i < numfaces; i++) {
03933           Face *face = getFaceAt(i);
03934           Vertex *vertex = face->getHashNode();
03935           mapfaces[vertex].push_back(face);
03936      }
03937 
03938      for (size_t i = 0; i < numfaces; i++) {
03939           Face *face = getFaceAt(i);
03940           const FaceSequence &hashfaces = mapfaces[face->getHashNode() ];
03941           size_t ncount = 0;
03942           for (size_t j = 0; j < hashfaces.size(); j++)
03943                if (hashfaces[j]->isSameAs(face)) ncount++;
03944           if (ncount != 1) {
03945                cout << "Warning: Mesh has some duplicate faces " << endl;
03946                return 1;
03947           }
03948      }
03949      return 0;
03950 }
03951 
03953 
03954 size_t
03955 Mesh::count_concave_faces()
03956 {
03957      size_t numfaces = getSize(2);
03958      size_t ncount = 0;
03959      for (size_t i = 0; i < numfaces; i++) {
03960           Face *f = getFaceAt(i);
03961           if(f->isActive()  &&  !f->isConvex() ) ncount++;
03962      }
03963      return ncount;
03964 }
03965 
03967 
03968 size_t
03969 Mesh::count_inverted_faces()
03970 {
03971      size_t numfaces = getSize(2);
03972      size_t ncount = 0;
03973      for (size_t i = 0; i < numfaces; i++) {
03974           Face *f = getFaceAt(i);
03975           if( !f->isRemoved()  && (f->concaveAt() >= 0)) ncount++;
03976      }
03977      return ncount;
03978 }
03980 
03981 void
03982 Mesh::normalize()
03983 {
03984      size_t numnodes = nodes.size();
03985 
03986      double xmin, xmax, ymin, ymax, zmin, zmax;
03987      Point3D xyz;
03988      xyz = nodes[0]->getXYZCoords();
03989 
03990      xmin = xyz[0];
03991      xmax = xyz[0];
03992      ymin = xyz[1];
03993      ymax = xyz[1];
03994      zmin = xyz[2];
03995      zmax = xyz[2];
03996 
03997      for (size_t i = 0; i < numnodes; i++) {
03998           Vertex *v = nodes[i];
03999           if( v->isActive() ) {
04000                xyz = v->getXYZCoords();
04001                xmin = min(xmin, xyz[0]);
04002                xmax = max(xmax, xyz[0]);
04003                ymin = min(ymin, xyz[1]);
04004                ymax = max(ymax, xyz[1]);
04005                zmin = min(zmin, xyz[2]);
04006                zmax = max(zmax, xyz[2]);
04007           }
04008      }
04009 
04010      double xlen = fabs(xmax - xmin);
04011      double ylen = fabs(ymax - ymin);
04012      double zlen = fabs(zmax - zmin);
04013      double scale = max(max(xlen, ylen), zlen);
04014 
04015      for (size_t i = 0; i < numnodes; i++) {
04016           Vertex *v = nodes[i];
04017           if( v->isActive() ) {
04018                xyz = v->getXYZCoords();
04019                xyz[0] /= scale;
04020                xyz[1] /= scale;
04021                xyz[2] /= scale;
04022                v->setXYZCoords(xyz);
04023           }
04024      }
04025 }
04026 
04028 
04029 BoundingBox
04030 Mesh::getBoundingBox() const
04031 {
04032      BoundingBox box;
04033 
04034      size_t numnodes = nodes.size();
04035 
04036      double xmin, xmax, ymin, ymax, zmin, zmax;
04037      Point3D xyz;
04038      xyz = nodes[0]->getXYZCoords();
04039 
04040      xmin = xyz[0];
04041      xmax = xyz[0];
04042      ymin = xyz[1];
04043      ymax = xyz[1];
04044      zmin = xyz[2];
04045      zmax = xyz[2];
04046 
04047      for (size_t i = 0; i < numnodes; i++) {
04048           Vertex *v = nodes[i];
04049           if( v->isActive() ) {
04050                xyz = nodes[i]->getXYZCoords();
04051                xmin = min(xmin, xyz[0]);
04052                xmax = max(xmax, xyz[0]);
04053                ymin = min(ymin, xyz[1]);
04054                ymax = max(ymax, xyz[1]);
04055                zmin = min(zmin, xyz[2]);
04056                zmax = max(zmax, xyz[2]);
04057           }
04058      }
04059 
04060      xyz[0] = xmin;
04061      xyz[1] = ymin;
04062      xyz[2] = zmin;
04063      box.setLowerLeftCorner(xyz);
04064 
04065      xyz[0] = xmax;
04066      xyz[1] = ymax;
04067      xyz[2] = zmax;
04068      box.setUpperRightCorner(xyz);
04069 
04070      return box;
04071 }
04072 
04073 double
04074 Mesh::getLength(int dir) const
04075 {
04076      size_t numnodes = nodes.size();
04077 
04078      double xmin, xmax, ymin, ymax, zmin, zmax;
04079      Point3D xyz;
04080      xyz = nodes[0]->getXYZCoords();
04081 
04082      xmin = xyz[0];
04083      xmax = xyz[0];
04084      ymin = xyz[1];
04085      ymax = xyz[1];
04086      zmin = xyz[2];
04087      zmax = xyz[2];
04088 
04089      for (size_t i = 0; i < numnodes; i++) {
04090           xyz = nodes[i]->getXYZCoords();
04091           xmin = min(xmin, xyz[0]);
04092           xmax = max(xmax, xyz[0]);
04093           ymin = min(ymin, xyz[1]);
04094           ymax = max(ymax, xyz[1]);
04095           zmin = min(zmin, xyz[2]);
04096           zmax = max(zmax, xyz[2]);
04097      }
04098 
04099      double xlen = fabs(xmax - xmin);
04100      double ylen = fabs(ymax - ymin);
04101      double zlen = fabs(zmax - zmin);
04102 
04103      if (dir == 0) return xlen;
04104      if (dir == 1) return ylen;
04105      if (dir == 2) return zlen;
04106 
04107      cout << "Warning: Invalid direction : " << dir << endl;
04108      return 0.0;
04109 }
04110 
04112 
04113 void
04114 Jaal::linear_interpolation(Mesh *mesh, Vertex *v0, Vertex *v1, int n, NodeSequence &newnodes)
04115 {
04116      assert(n >= 2);
04117 
04118      const Point3D &xyz0 = v0->getXYZCoords();
04119      const Point3D &xyz1 = v1->getXYZCoords();
04120 
04121      NodeSequence poolnodes;
04122      mesh->objects_from_pool(n-2, poolnodes);
04123 
04124      newnodes.resize(n);
04125 
04126      newnodes[0] = v0;
04127      newnodes[n - 1] = v1;
04128 
04129      double dt = 2.0 / (double) (n - 1);
04130 
04131      Point3D xyzt;
04132      int index = 0;
04133      for (int i = 1; i < n - 1; i++) {
04134           double t = -1.0 + i*dt;
04135           xyzt[0] = TFI::linear_interpolation(t, xyz0[0], xyz1[0]);
04136           xyzt[1] = TFI::linear_interpolation(t, xyz0[1], xyz1[1]);
04137           xyzt[2] = TFI::linear_interpolation(t, xyz0[2], xyz1[2]);
04138           newnodes[i] = poolnodes[index++];
04139           newnodes[i]->setStatus( MeshEntity::ACTIVE);
04140           newnodes[i]->setXYZCoords(xyzt);
04141      }
04142 }
04143 
04145 
04146 void
04147 Jaal::set_tfi_coords(int i, int j, int nx, int ny, vector<Vertex*> &qnodes)
04148 {
04149      int offset;
04150 
04151      offset = 0;
04152      const Point3D &v00 = qnodes[offset]->getXYZCoords();
04153 
04154      offset = i;
04155      const Point3D &vr0 = qnodes[offset]->getXYZCoords();
04156 
04157      offset = (nx - 1);
04158      const Point3D &v10 = qnodes[offset]->getXYZCoords();
04159 
04160      offset = j*nx;
04161      const Point3D &v0s = qnodes[offset]->getXYZCoords();
04162 
04163      offset = j * nx + (nx - 1);
04164      const Point3D &v1s = qnodes[offset]->getXYZCoords();
04165 
04166      offset = (ny - 1) * nx;
04167      const Point3D &v01 = qnodes[offset]->getXYZCoords();
04168 
04169      offset = (ny - 1) * nx + i;
04170      const Point3D &vr1 = qnodes[offset]->getXYZCoords();
04171 
04172      offset = (ny - 1) * nx + (nx - 1);
04173      const Point3D &v11 = qnodes[offset]->getXYZCoords();
04174 
04175      Point3D vrs;
04176 
04177      double dr = 2.0 / (double) (nx - 1);
04178      double ds = 2.0 / (double) (ny - 1);
04179 
04180      double r = -1.0 + i*dr;
04181      double s = -1.0 + j*ds;
04182      for (int k = 0; k < 3; k++) {
04183           vrs[k] = TFI::transfinite_blend(r, s,
04184                                           v00[k], v10[k], v11[k], v01[k],
04185                                           vr0[k], v1s[k], vr1[k], v0s[k]);
04186      }
04187      offset = j * nx + i;
04188      qnodes[offset]->setXYZCoords(vrs);
04189 
04190 }
04192 
04193 int
04194 QTrack::advance_single_step(int endat)
04195 {
04197      //                   **********************
04198      //                   *         *          *
04199      //                   *         * Next     *
04200      //                   *         *          *
04201      //           Avoid   **********************  Avoid
04202      //                   *         *          *
04203      //                   *         * Current  *
04204      //                   *         *          *
04205      //                   *         *          *
04206      //                   **********************
04207      //                            Source
04208      // A Source vertex and Current edge is chosen.
04209      // We want to avoid two edges and want to select "Next" edge.
04211      Vertex *v0, *v1, *v2, *v3, *v4;
04212 
04213      NodeSet vset;
04214 
04215      size_t index = sequence.size();
04216      v0 = sequence[index - 2];
04217      v1 = sequence[index - 1];
04218      v0->setVisitMark(1);
04219 
04220      if (endat == END_AT_CROSSINGS && v1->isVisited()) return 0;
04221      if (v1->isBoundary()) return 0;
04222 
04223      v1->setVisitMark(1);
04224      NodeSequence vneighs;
04225      v1->getRelations( vneighs );
04226      if (vneighs.size() != 4) return 0;
04227 
04228      FaceSequence adjFaces;
04229      Mesh::getRelations112(v0, v1, adjFaces);
04230      assert(adjFaces.size() == 2);
04231      v2 = Face::opposite_node(adjFaces[0], v0);
04232      v3 = Face::opposite_node(adjFaces[1], v0);
04233 
04234      vset.clear();
04235      vset.insert(vneighs[0]);
04236      vset.insert(vneighs[1]);
04237      vset.insert(vneighs[2]);
04238      vset.insert(vneighs[3]);
04239      vset.erase(v0);
04240      vset.erase(v2);
04241      vset.erase(v3);
04242      assert(vset.size() == 1);
04243      v4 = *vset.begin();
04244      sequence.push_back(v4);
04245      return 1;
04246 }
04247 
04249 
04250 void
04251 QTrack::advance(int endat)
04252 {
04253      assert(sequence.size() == 2);
04254 
04255      // Starting node is always irregular ...
04256      assert( sequence[0]->getNumRelations(2) != 4);
04257 
04258      while (1) {
04259           int progress = advance_single_step(endat);
04260           if (!progress) break;
04261      }
04262 
04263 #ifdef DEBUG
04264      Vertex *endvertex;
04265      // Sanity Checking ....
04266      if (endat == END_AT_TERMINALS) {
04267           endvertex = sequence.front();
04268           if (!endvertex->isBoundary()) {
04269                assert( endvertex->getNumRelations(2) != 4);
04270           }
04271 
04272           endvertex = sequence.back();
04273           if (!endvertex->isBoundary()) {
04274                assert( endvertex->getNumRelations(2) != 4);
04275           }
04276 
04277           for (size_t i = 1; i < sequence.size() - 1; i++) {
04278                assert(!sequence[i]->isBoundary());
04279                assert( endvertex->getNumRelations(2) == 4);
04280           }
04281      }
04282 #endif
04283 
04284 }
04285 
04287 
04288 vector<QTrack> Jaal::generate_quad_partitioning(Mesh *mesh)
04289 {
04290      vector<QTrack> qpath;
04291 
04292 #ifdef CHANGE
04293      int nTopo = mesh->isHomogeneous();
04294      if (nTopo != 4) {
04295           cout << "Error: The mesh must be all Quads " << endl;
04296           return qpath;
04297      }
04298 
04299      int relexist2 = mesh->build_relations(0, 2);
04300      int relexist0 = mesh->build_relations(0, 0);
04301 
04302      mesh->search_boundary();
04303 
04304      size_t numnodes = mesh->getSize(0);
04305      for (size_t i = 0; i < numnodes; i++) {
04306           Vertex *vertex = mesh->getNodeAt(i);
04307           vertex->setAttribute("Partition", 0);
04308           vertex->setVisitMark(0);
04309      }
04310 
04311      QTrack qp;
04312      qp.mesh = mesh;
04313 
04314      NodeSequence vnodes;
04315 
04316      int found;
04317      for (size_t i = 0; i < numnodes; i++) {
04318           Vertex *vertex = mesh->getNodeAt(i);
04319           if (!vertex->isBoundary() && !vertex->isRemoved() ) {
04320                vertex->getRelations( vnodes );
04321                int numneighs = vnodes.size();
04322                if(numneighs != 4) {
04323                     for (int j = 0; j < numneighs; j++) {
04324                          qp.sequence.resize(2); // As we know the starting edge
04325                          qp.sequence[0] = vertex;
04326                          qp.sequence[1] = vnodes[j];
04327                          qp.advance(0); // Terminate at irregular nodes only..
04328                          found = 0;
04329                          for (size_t k = 0; k < qpath.size(); k++) {
04330                               if (qpath[k] == qp) {
04331                                    found = 1;
04332                                    break;
04333                               }
04334                          }
04335                          if (!found) qpath.push_back(qp);
04336                     }
04337                }
04338           }
04339      }
04340      sort(qpath.begin(), qpath.end());
04341 
04342      int partid = 1;
04343      for( size_t i = 0; i < qpath.size(); i++) {
04344           for( size_t j = 0; j < qpath[i].sequence.size(); j++)
04345                qpath[i].sequence[j]->setAttribute("Partition", partid );
04346           partid++;
04347      }
04348 
04349      size_t nSize = mesh->getSize(2);
04350      for( size_t i = 0; i < nSize; i++) {
04351           Face *f = mesh->getFaceAt(i);
04352           f->setAttribute( "Partition", 0);
04353      }
04354 
04355      deque<Face*> faceQ;
04356      map<int, int> faceCount;
04357      FaceSequence fneighs;
04358      partid = 1;
04359      int pid;
04360      while(1) {
04361           Face *seedface = NULL;
04362           for( size_t i = 0; i < nSize; i++) {
04363                Face *f = mesh->getFaceAt(i);
04364                if( f->isActive()  ) {
04365                     f->getAttribute("Partition", pid);
04366                     if( pid == 0) {
04367                          seedface = f;
04368                          break;
04369                     }
04370                }
04371           }
04372           if( seedface == NULL ) break;
04373           faceCount[partid] = 0;
04374 
04375           faceQ.clear();
04376           faceQ.push_back(seedface);
04377 
04378           while(!faceQ.empty() ) {
04379 
04380                Face *currface = faceQ.front();
04381                faceQ.pop_front();
04382 
04383                if( currface->isActive() && currface->getGroupID() == 0 ) {
04384                     currface->setAttribute( "Partition", partid) ;
04385                     faceCount[partid]++;
04386                     for( int i = 0; i < 4; i++) {
04387                          Vertex *v0 = currface->getNodeAt(i);
04388                          Vertex *v1 = currface->getNodeAt(i+1);
04389                          if( v0->getGroupID() == 0 || v1->getGroupID() == 0) {
04390                               Mesh::getRelations112( v0, v1, fneighs);
04391                               for( size_t j = 0; j < fneighs.size(); j++)
04392                                    if( fneighs[j]->getGroupID() == 0)
04393                                         faceQ.push_back( fneighs[j] );
04394                          }
04395                     }
04396                }
04397           }
04398           partid++;
04399      }
04400 
04401      cout << "#of Partitions : " << faceCount.size() << endl;
04402      map<int,int>::const_iterator it;
04403      for( it = faceCount.begin(); it != faceCount.end(); ++it)
04404           cout << it->first << " " << it->second << endl;
04405 
04406 
04407      if (!relexist2) mesh->clear_relations(0, 2);
04408      if (!relexist0) mesh->clear_relations(0, 0);
04409 
04410 #endif
04411      return qpath;
04412 }
04413 
04415 
04416 void
04417 Jaal::set_layer_tag(Mesh *mesh, const string &name )
04418 {
04419      assert(1);
04420      //  mesh->setWavefront (0);
04421      int lid;
04422 
04423      size_t numnodes = mesh->getSize(0);
04424      int max_tag_val = -1;
04425      for (size_t i = 0; i < numnodes; i++) {
04426           Vertex *vertex = mesh->getNodeAt(i);
04427           vertex->getAttribute(name, lid );
04428           max_tag_val = max(max_tag_val, lid );
04429      }
04430 
04431      for (size_t i = 0; i < numnodes; i++) {
04432           Vertex *vertex = mesh->getNodeAt(i);
04433           vertex->getAttribute(name, lid );
04434           if (lid < 0)
04435                vertex->setAttribute(name, max_tag_val + 1);
04436           else
04437                vertex->setAttribute(name, lid);
04438      }
04439 
04440      //  mesh->setWavefront (2);
04441      size_t numfaces = mesh->getSize(2);
04442 
04443      max_tag_val = -1;
04444      for (size_t i = 0; i < numfaces; i++) {
04445           Face *face = mesh->getFaceAt(i);
04446           face->getAttribute(name, lid );
04447           max_tag_val = max(max_tag_val, lid);
04448      }
04449 
04450      for (size_t i = 0; i < numfaces; i++) {
04451           Face *face = mesh->getFaceAt(i);
04452           face->getAttribute(name, lid);
04453           if (lid < 0)
04454                face->setAttribute(name, max_tag_val);
04455           else
04456                face->setAttribute(name, lid);
04457      }
04458 }
04459 
04461 
04462 void
04463 Jaal::set_convexity_tag(Mesh *mesh, const string &name)
04464 {
04465      size_t numnodes = mesh->getSize(0);
04466      for (size_t i = 0; i < numnodes; i++) {
04467           Vertex *vertex = mesh->getNodeAt(i);
04468           vertex->setAttribute(name, 1);
04469      }
04470 
04471      size_t numfaces = mesh->getSize(2);
04472      for (size_t i = 0; i < numfaces; i++) {
04473           Face *face = mesh->getFaceAt(i);
04474           face->setAttribute(name, 0);
04475           if (face->isConvex())
04476                face->setAttribute(name, 1);
04477           else {
04478                for (int j = 0; j < face->getSize(0); j++) {
04479                     Vertex *v = face->getNodeAt(j);
04480                     v->setAttribute(name, 0);
04481                }
04482           }
04483 
04484      }
04485 }
04486 
04488 
04489 void
04490 Jaal::set_boundary_tag(Mesh *mesh, const string &name)
04491 {
04492      int relexist = mesh->build_relations(0, 2);
04493      cout << "Set Boundary Tag : " << endl;
04494 
04495      if (!mesh->isBoundaryKnown())
04496           mesh->search_boundary();
04497 
04498      size_t numnodes = mesh->getSize(0);
04499 
04500      for (size_t i = 0; i < numnodes; i++) {
04501           Vertex *vertex = mesh->getNodeAt(i);
04502           if( vertex->isActive() ) {
04503                if (vertex->isBoundary()) {
04504                     vertex->setAttribute(name, 1);
04505                } else
04506                     vertex->setAttribute(name, 0);
04507           }
04508      }
04509 
04510      size_t numfaces = mesh->getSize(2);
04511 
04512      for (size_t i = 0; i < numfaces; i++) {
04513           Face *face = mesh->getFaceAt(i);
04514           if( face->isActive() )
04515                face->setAttribute(name, 0);
04516      }
04517 
04518      for (size_t i = 0; i < numfaces; i++) {
04519           Face *face = mesh->getFaceAt(i);
04520           if( face->isActive() )
04521                if (face->has_boundary_edge())
04522                     face->setAttribute(name, 1);
04523      }
04524 
04525      if (!relexist) mesh->clear_relations(0, 2);
04526 
04527 }
04528 
04530 
04531 void
04532 Jaal::set_constrained_tag(Mesh *mesh)
04533 {
04534      size_t numnodes = mesh->getSize(0);
04535 
04536      for (size_t i = 0; i < numnodes; i++) {
04537           Vertex *vertex = mesh->getNodeAt(i);
04538           if (vertex->isConstrained())
04539                vertex->setAttribute("Constraint", 1);
04540           else
04541                vertex->setAttribute("Constraint", 0);
04542      }
04543 }
04544 
04546 
04547 bool
04548 Face::has_all_bound_nodes() const
04549 {
04550      for (int i = 0; i < getSize(0); i++) {
04551           Vertex *v = getNodeAt(i);
04552           if (!v->isBoundary()) return 0;
04553      }
04554      return 1;
04555 }
04557 
04558 void
04559 Mesh::setNormals()
04560 {
04561      double x[100], y[100], z[100], nx, ny, nz;
04562      Point3D xyz;
04563      Vec3D normal;
04564 
04565      size_t nSize = faces.size();
04566      for (size_t i = 0; i < nSize; i++) {
04567           Face *face = faces[i];
04568           if( face->isActive() ) {
04569                int nsize = face->getSize(0);
04570                for (int j = 0; j < nsize; j++) {
04571                     Vertex *vtx = face->getNodeAt(j);
04572                     xyz = vtx->getXYZCoords();
04573                     x[j] = xyz[0];
04574                     y[j] = xyz[1];
04575                     z[j] = xyz[2];
04576                }
04577                PolygonNormal3D(nsize, x, y, z, &nx, &ny, &nz);
04578                normal[0] = nx;
04579                normal[1] = ny;
04580                normal[2] = nz;
04581                face->setAttribute("Normal", normal);
04582           }
04583      }
04584      int relexist2 = build_relations(0,2);
04585 
04586      FaceSequence vfaces;
04587      nSize = nodes.size();
04588      for (size_t i = 0; i < nSize; i++) {
04589           Vertex *vertex = nodes[i];
04590           if( vertex->isActive() ) {
04591                vertex->getRelations(vfaces);
04592                double nx = 0.0;
04593                double ny = 0.0;
04594                double nz = 0.0;
04595                for( size_t j = 0; j < vfaces.size(); j++) {
04596                     vfaces[j]->getAttribute("Normal", normal);
04597                     nx += normal[0];
04598                     ny += normal[1];
04599                     nz += normal[2];
04600                }
04601                double multby = 1.0/sqrt( nx*nx + ny*ny + nz*nz );
04602                normal[0] = nx*multby;
04603                normal[1] = ny*multby;
04604                normal[2] = nz*multby;
04605                vertex->setAttribute("Normal", normal);
04606           }
04607      }
04608      if( !relexist2 )  clear_relations(0,2);
04609 }
04610 
04612 void Mesh :: objects_from_pool( size_t n, vector<Vertex*> &objects)
04613 {
04614      objects.clear();
04615 
04616      if( n == 0) return;
04617 
04618      objects.reserve( n );
04619 
04620      size_t ncount = 0;
04621      while( !garbageNodes.empty() ) {
04622           Vertex *v = garbageNodes.front();
04623           garbageNodes.pop_front();
04624           if( v->isRemoved() ) {
04625                v->setStatus( MeshEntity::INACTIVE);
04626                objects.push_back(v);
04627                ncount++;
04628                if( ncount == n) break;
04629           }
04630      }
04631 
04632      for( size_t i = ncount; i < n; i++) {
04633           Vertex *v = Vertex::newObject();
04634           v->setStatus( MeshEntity::INACTIVE);
04635           objects.push_back(v);
04636           addNode( v );
04637      }
04638 
04639      assert( objects.size() == n );
04640 }
04641 
04643 
04644 void Mesh :: objects_from_pool( size_t n, vector<Face*> &objects)
04645 {
04646      objects.clear();
04647 
04648      if( n == 0) return;
04649 
04650      objects.reserve( n );
04651 
04652      size_t ncount = 0;
04653      while( !garbageFaces.empty() ) {
04654           Face *f = garbageFaces.front();
04655           garbageFaces.pop_front();
04656           if( f->isRemoved() ) {
04657                f->setStatus( MeshEntity::INACTIVE);
04658                objects.push_back(f);
04659                ncount++;
04660                if( ncount == n) break;
04661           }
04662      }
04663 
04664      for( size_t i = ncount; i < n; i++) {
04665           Face *f = Face::newObject();
04666           f->setStatus( MeshEntity::INACTIVE);
04667           objects.push_back(f);
04668           addFace(f);
04669      }
04670 
04671      assert( objects.size() == n );
04672 }
04674 
04675 int Mesh::get_breadth_first_ordered_nodes(NodeSequence &seq, Vertex *vstart, MeshFilter *filter)
04676 {
04677      assert(vstart != NULL);
04678 
04679      seq.clear();
04680 
04681      int relexist0 = build_relations(0, 0);
04682 
04683      size_t numnodes = getSize(0);
04684 
04685      if (numnodes == 0) return 1;
04686 
04687      for (size_t i = 0; i < numnodes; i++) {
04688           Vertex *v = getNodeAt(i);
04689           v->setVisitMark(0);
04690           v->setAttribute("Layer", 0);
04691      }
04692 
04693      if (vstart == 0) vstart = getNodeAt(0);
04694 
04695      seq.reserve(numnodes);
04696 
04697      list<Vertex*> vertexQ;
04698      vertexQ.push_back(vstart);
04699      NodeSequence vneighs;
04700 
04701      int proceed = 1;
04702      int currlevel;
04703      while (!vertexQ.empty()) {
04704           Vertex *curr_vertex = vertexQ.front();
04705           vertexQ.pop_front();
04706           curr_vertex->getAttribute("Layer", currlevel);
04707           if (filter) {
04708                if (curr_vertex != vstart) proceed = filter->pass(curr_vertex);
04709           }
04710           if (!curr_vertex->isVisited()) {
04711                seq.push_back(curr_vertex);
04712                if (!proceed) break;
04713                curr_vertex->setVisitMark(1);
04714                curr_vertex->getRelations( vneighs );
04715                for (size_t i = 0; i < vneighs.size(); i++) {
04716                     if (!vneighs[i]->isVisited()) {
04717                          vertexQ.push_back(vneighs[i]);
04718                          vneighs[i]->setAttribute("Layer", currlevel + 1);
04719                     }
04720                }
04721           }
04722      }
04723 
04724      if (!relexist0)
04725           clear_relations(0, 0);
04726 
04727      // Free unused memory in sequence...
04728      if (!seq.empty()) NodeSequence(seq).swap(seq);
04729 
04730      return 0;
04731 }
04733 
04734 int Mesh::get_depth_first_ordered_nodes(NodeSequence &seq, Vertex *vstart, MeshFilter *filter)
04735 {
04736      int relexist0 = build_relations(0, 0);
04737 
04738      size_t numnodes = getSize(0);
04739      list<Vertex*> vertexQ;
04740      for (size_t i = 0; i < numnodes; i++) {
04741           Vertex *v = getNodeAt(i);
04742           v->setVisitMark(0);
04743      }
04744 
04745      seq.clear();
04746 
04747      if (vstart == 0) vstart = getNodeAt(0);
04748      vertexQ.push_back(vstart);
04749      NodeSequence vneighs;
04750 
04751      while (!vertexQ.empty()) {
04752           Vertex *curr_vertex = vertexQ.front();
04753           vertexQ.pop_front();
04754           if (!curr_vertex->isVisited()) {
04755                seq.push_back(curr_vertex);
04756                curr_vertex->setVisitMark(1);
04757                curr_vertex->getRelations( vneighs );
04758                for (size_t i = 0; i < vneighs.size(); i++) {
04759                     if (!vneighs[i]->isVisited())
04760                          vertexQ.push_front(vneighs[i]);
04761                }
04762           }
04763      }
04764 
04765      if (!relexist0)
04766           clear_relations(0, 0);
04767 
04768      return 0;
04769 }
04770 
04772 
04773 int Mesh::get_sharp_edges(EdgeSequence &sharp_edges, double creaseAngle)
04774 {
04775      sharp_edges.clear();
04776      /*
04777         if( edges.empty() ) build_edges();
04778 
04779         FaceSequence efaces;
04780         for( size_t i = 0; i < edges.size(): i++) {
04781              PEdge edge = getEdge(i);
04782              Vertex *v1 = edge->getNodeAt(0);
04783              Vertex *v1 = edge->getNodeAt(1);
04784              efaces = Mesh::getRelations112(v0, v1);
04785              if( efaces.size() == 2 ) {
04786                  Vec3D f1normal = efaces[0]->getNormal();
04787                  Vec3D f2normal = efaces[1]->getNormal();
04788                  double angle = Math::getAngle(fn1, fn2, ANGLE_IN_DEGREES);
04789                  if (angle <= 90 && angle >= creaseAngle)
04790                      sharp_edges.push_back(edge->getClone() );
04791                  else if (angle >= 90 && fabs(180 - angle) >= creaseAngle)
04792                      sharp_edges.push_back(edge->getClone() );
04793              }
04794          }
04795       */
04796      return 0;
04797 }
04798 
04800 
04801 void
04802 Jaal::set_inverted_tag(Mesh *mesh, const string &name)
04803 {
04804      size_t numnodes = mesh->getSize(0);
04805      for (size_t i = 0; i < numnodes; i++) {
04806           Vertex *v = mesh->getNodeAt(i);
04807           v->setAttribute(name, 1);
04808      }
04809 
04810      size_t ncount = 0;
04811      size_t numfaces = mesh->getSize(2);
04812      for (size_t i = 0; i < numfaces; i++) {
04813           Face *face = mesh->getFaceAt(i);
04814           face->setAttribute(name, 1);
04815           int pos = face->concaveAt();
04816           if (pos >= 0) {
04817                face->setAttribute(name, 0);
04818                Vertex *v = face->getNodeAt(pos);
04819                v->setAttribute(name, 0);
04820                ncount++;
04821           }
04822      }
04823 }
04824 
04826 
04827 void
04828 Jaal::set_irregular_path_tag(Mesh *mesh, vector<QTrack> &qpath)
04829 {
04830      size_t numnodes = mesh->getSize(0);
04831      Vertex *v;
04832      for (size_t i = 0; i < numnodes; i++) {
04833           v = mesh->getNodeAt(i);
04834           v->setAttribute( "QuadPatch", 0);
04835      }
04836 
04837      for (size_t i = 0; i < qpath.size(); i++) {
04838           v = qpath[i].sequence.front();
04839           v->setAttribute("QuadPatch",1);
04840           v = qpath[i].sequence.back();
04841           v->setAttribute("QuadPatch",1);
04842 
04843           for (size_t j = 1; j < qpath[i].sequence.size() - 1; j++) {
04844                v = qpath[i].sequence[j];
04845                v->setAttribute("QuadPatch",2);
04846           }
04847      }
04848 }
04849 
04851 
04852 int Jaal::SurfPatch::getPosOf(const Vertex *v)
04853 {
04854      for (size_t i = 0; i < bound_nodes.size(); i++)
04855           if (bound_nodes[i] == v) return i;
04856 
04857      cout << "Error: Vertex not found on the boundary " << endl;
04858      exit(0);
04859 
04860      return -1;
04861 }
04862 
04864 
04865 NodeSequence SurfPatch::get_bound_nodes(const Vertex *src, const Vertex *dst)
04866 {
04867      int start_pos = getPosOf(src);
04868      int end_pos = getPosOf(dst);
04869      int nsize = bound_nodes.size();
04870 
04871      if (end_pos == 0) end_pos = nsize;
04872      assert(end_pos > start_pos);
04873 
04874      NodeSequence seq(end_pos - start_pos + 1);
04875      int index = 0;
04876      for (int i = start_pos; i <= end_pos; i++)
04877           seq[index++] = bound_nodes[i % nsize];
04878 
04879      return seq;
04880 }
04881 
04883 
04884 int Jaal::SurfPatch::search_boundary()
04885 {
04886      corners.clear();
04887      boundary.clear();
04888      Vertex *vertex;
04889 
04890      assert(!faces.empty());
04891 
04892      // We need to rebuild relations locally to identfy corners and boundary.
04893      FaceSet::const_iterator fiter;
04894      std::map<Vertex*, FaceSet> relations02;
04895 
04896      for (fiter = faces.begin(); fiter != faces.end(); ++fiter) {
04897           Face *face = *fiter;
04898           for (int j = 0; j < face->getSize(0); j++) {
04899                vertex = face->getNodeAt(j);
04900                relations02[vertex].insert(face);
04901           }
04902      }
04903 
04904      // A boundary edge must have exactly one face neighbor...
04905      FaceSequence faceneighs;
04906      for (fiter = faces.begin(); fiter != faces.end(); ++fiter) {
04907           Face *face = *fiter;
04908           int nnodes = face->getSize(0);
04909           for (int j = 0; j < nnodes; j++) {
04910                Vertex *v0 = face->getNodeAt( j + 0 );
04911                Vertex *v1 = face->getNodeAt( j + 1 );
04912                faceneighs.clear();
04913                assert(relations02[v0].size() > 0);
04914                assert(relations02[v1].size() > 0);
04915                set_intersection(relations02[v0].begin(), relations02[v0].end(),
04916                                 relations02[v1].begin(), relations02[v1].end(),
04917                                 back_inserter(faceneighs));
04918                if (faceneighs.size() == 1) {
04919                     Edge newedge(v0, v1);
04920                     boundary.push_back(newedge);
04921                }
04922           }
04923      }
04924 
04925      // Sequence the chain and start from one of the corner...
04926      int err = Mesh::make_chain(boundary);
04927      if (err) return 2;
04928 
04929      //
04930      // Identify corners in the mesh.
04931      // Should we check only one vertex per edge ?
04932      //
04933      FaceSet neighs;
04934      int boundSize = boundary.size();
04935      for (int k = 0; k < boundSize; k++) {
04936           vertex = boundary[k].getNodeAt(0);
04937           neighs = relations02[vertex];
04938           if (neighs.size() == 1) corners.insert(vertex);
04939 
04940           vertex = boundary[k].getNodeAt(1);
04941           neighs = relations02[vertex];
04942           if (neighs.size() == 1) corners.insert(vertex);
04943      }
04944 
04945      // Start the chain from one of the corners.
04946      err = Mesh::rotate_chain(boundary, *corners.begin());
04947      if (err) return 3;
04948 
04949      // Collect the sequence of boundary nodes...,
04950      bound_nodes.resize( boundSize );
04951      for (int k = 0; k < boundSize; k++) {
04952           bound_nodes[k] = boundary[k].getNodeAt(0); // Only the first node.
04953      }
04954 
04955      //
04956      // Collect the inner nodes of the patch. These nodes will be deleted, if
04957      // the remesh operation is successful...
04958      //
04959      inner_nodes.clear();
04960      for (fiter = faces.begin(); fiter != faces.end(); ++fiter) {
04961           Face *face = *fiter;
04962           int nnodes = face->getSize(0);
04963           for (int j = 0; j < nnodes; j++) {
04964                Vertex *v = face->getNodeAt(j);
04965                if (find(bound_nodes.begin(), bound_nodes.end(), v) == bound_nodes.end())
04966                     inner_nodes.insert(v);
04967           }
04968      }
04969 
04970      // Split the boundary loop into segments.
04971      // (i.e. End of the segments are the corners identified earlier )
04972      set_boundary_segments();
04973 
04974      return 0;
04975 }
04976 
04978 
04979 void Jaal::SurfPatch::set_boundary_segments()
04980 {
04981      // Although this stage will not come in this algorithm...
04982      if (corners.size() == 0) return;
04983 
04984      cornerPos.resize(corners.size() + 1);
04985 
04986      NodeSet::const_iterator it;
04987      int index = 0;
04988      for (it = corners.begin(); it != corners.end(); ++it) {
04989           cornerPos[index++] = getPosOf(*it);
04990      }
04991 
04992      cornerPos[corners.size()] = bound_nodes.size();
04993 
04994      sort(cornerPos.begin(), cornerPos.end());
04995 
04996      segSize.resize(corners.size());
04997 
04998      for (size_t i = 0; i < corners.size(); i++)
04999           segSize[i] = cornerPos[(i + 1)] - cornerPos[i] + 1;
05000 }
05001 
05003 
05004 int Jaal::SurfPatch::reorient_4_sided_loop()
05005 {
05006      //Always remeshable, nothing has to be done...
05007      if ((segSize[0] == segSize[2]) && (segSize[1] == segSize[3])) return 0;
05008 
05010      // Defination:  A four sided convex loop has four segments.
05011      // Objectives:  A four sided convex loop must be orietned such that
05012      //   1.  First segment must be smaller than the third one, because
05013      //       we need to create a triangle patch based at segment#3.
05014      //
05015      //   2.  If there are two choices, then the side having irregular
05016      //       node must be given higher priority. ( Here irregulaty means
05017      //       that vertex valency < 4 ).
05018      //
05019      //   3.  A side having less number of nodes on the first segment than
05020      //       the third is given preference.
05021      //
05022      // Pre-Conditions  :  A loop must be oriented ( CW or CCW ).
05023      //
05024      // Date: 17th Nov. 2010.
05026 
05027      Vertex *start_corner = NULL;
05028 
05029      if (segSize[0] == segSize[2]) {
05030           if (min(segSize[1], segSize[3]) == 2) return 1;
05031           //  Either Segment 2 or 3 must be starting node
05032           if (segSize[1] < segSize[3])
05033                start_corner = bound_nodes[ cornerPos[1] ];
05034           else
05035                start_corner = bound_nodes[ cornerPos[3] ];
05036           start_boundary_loop_from(start_corner);
05037      }
05038 
05039      if (min(segSize[0], segSize[2]) == 2) return 1;
05040 
05041      if (segSize[2] < segSize[0]) {
05042           start_corner = bound_nodes[ cornerPos[2] ];
05043           start_boundary_loop_from(start_corner);
05044      }
05045 
05046      // By this stage, the loop must be reoriented correctly.
05047      assert(segSize[0] < segSize[2]);
05048 
05049      cout << " Careful to change this code " << endl;
05050      // Great, we found one irregular node on the first boundary...
05051      //   if( has_irregular_node_on_first_segment() ) return 1;
05052 
05053      // If the segment 2 and 3 have same size, Alas, nothing can be done.
05054      if (segSize[1] == segSize[3]) return 1;
05055 
05056      if (min(segSize[1], segSize[3]) == 2) return 1;
05057 
05058      if (segSize[3] < segSize[1]) {
05059           start_corner = bound_nodes[ cornerPos[3] ];
05060           start_boundary_loop_from(start_corner);
05061      } else {
05062           start_corner = bound_nodes[ cornerPos[1] ];
05063           start_boundary_loop_from(start_corner);
05064      }
05065 
05066      //
05067      // Note that we didn't check for irregular node here. So if this segment
05068      // has at least one irregular node, then we are lucky. Otherwise decision
05069      // to remesh it done based wthere remeshing will result in the reduction
05070      // of irregular nodes in patch.
05071      //
05072      assert(segSize[0] < segSize[2]);
05073      return 0;
05074 }
05076 
05077 void Jaal::SurfPatch::start_boundary_loop_from(Vertex *vmid)
05078 {
05079      assert(corners.find(vmid) != corners.end());
05080 
05081      NodeSequence::iterator middle;
05082      middle = find(bound_nodes.begin(), bound_nodes.end(), vmid);
05083      assert(middle != bound_nodes.end());
05084 
05085      std::rotate(bound_nodes.begin(), middle, bound_nodes.end());
05086      assert(bound_nodes[0] == vmid);
05087 
05088      set_boundary_segments();
05089 }
05091 
05092 Mesh * Jaal::create_structured_mesh(double *origin, double *length,
05093                                     int *grid_dim, int space_dim)
05094 {
05095      if (space_dim < 2 || space_dim > 3) return NULL;
05096 
05097      Mesh *mesh = new Mesh;
05098 
05099      int nx = grid_dim[0];
05100      int ny = grid_dim[1];
05101      int nz = 1;
05102      if (space_dim == 3) nz = grid_dim[2];
05103 
05104      double xorg, yorg, zorg = 0.0, dx, dy, dz = 0.0;
05105 
05106      xorg = origin[0];
05107      yorg = origin[1];
05108 
05109      if (space_dim == 3) zorg = origin[2];
05110 
05111      dx = length[0] / (double) (nx - 1);
05112      dy = length[1] / (double) (ny - 1);
05113      if (nz > 2) dz = length[2] / (double) (nz - 1);
05114 
05115 
05116      vector<Vertex*> nodes(nx * ny * nz);
05117      int offset;
05118      Point3D xyz;
05119      for (int k = 0; k < nz; k++) {
05120           for (int j = 0; j < ny; j++) {
05121                for (int i = 0; i < nx; i++) {
05122                     offset = k * nx * ny + j * nx + i;
05123                     Vertex *v = Vertex::newObject();
05124                     xyz[0] = xorg + i*dx;
05125                     xyz[1] = yorg + j*dy;
05126                     xyz[2] = zorg + k*dz;
05127                     v->setXYZCoords(xyz);
05128                     v->setID(offset);
05129                     nodes[ offset ] = v;
05130                     mesh->addNode(v);
05131                }
05132           }
05133      }
05134 
05135      NodeSequence connect;
05136 
05137      if (space_dim == 2) {
05138           connect.resize(4);
05139           for (int j = 0; j < ny - 1; j++) {
05140                for (int i = 0; i < nx - 1; i++) {
05141                     offset = j * nx + i;
05142                     connect[0] = nodes[offset];
05143                     connect[1] = nodes[offset + 1];
05144                     connect[2] = nodes[offset + 1 + nx];
05145                     connect[3] = nodes[offset + nx];
05146                     Face *face = Face::newObject();
05147                     face->setNodes(connect);
05148                     mesh->addFace(face);
05149                }
05150           }
05151      }
05152      return mesh;
05153 }
05154 
05156 
05157 int Jaal::quad_concave_tests()
05158 {
05159      NodeSequence nodes = Mesh::generate_nodes(4);
05160      Point3D xyz;
05161 
05162      xyz[0] = -1.0;
05163      xyz[1] = 0.0;
05164      xyz[2] = 0.0;
05165      nodes[0]->setXYZCoords(xyz);
05166 
05167      xyz[0] = 0.0;
05168      xyz[1] = -0.00001;
05169      xyz[2] = 0.0;
05170      nodes[1]->setXYZCoords(xyz);
05171 
05172      xyz[0] = 1.0;
05173      xyz[1] = 0.0;
05174      xyz[2] = 0.0;
05175      nodes[2]->setXYZCoords(xyz);
05176 
05177      xyz[0] = 0.0;
05178      xyz[1] = 1.0;
05179      xyz[2] = 0.0;
05180      nodes[3]->setXYZCoords(xyz);
05181 
05182      Face *face = Face::newObject();
05183      face->setNodes(nodes);
05184 
05185      assert(face->concaveAt() == -1);
05186 
05187      xyz[0] = -0.0;
05188      xyz[1] = 0.000001;
05189      xyz[2] = 0.0;
05190      nodes[1]->setXYZCoords(xyz);
05191      assert(face->concaveAt() == 1);
05192 
05193      return 0;
05194 }
05196 
05197 Mesh* Jaal::quad_to_tri4( Mesh *quadmesh, vector<Vertex*> &steiner)
05198 {
05199      if( quadmesh == NULL ) return NULL;
05200 
05201      Mesh *trimesh = new Mesh;
05202 
05203      size_t numnodes = quadmesh->getSize(0);
05204      size_t numfaces = quadmesh->getSize(2);
05205 
05206      trimesh->reserve( numnodes + numfaces, 0);
05207      trimesh->reserve( 4*numfaces, 2);
05208 
05209      for( size_t i = 0; i < numnodes; i++) {
05210           Vertex *v = quadmesh->getNodeAt(i);
05211           if( v->isActive() ) trimesh->addNode(v);
05212      }
05213      steiner.clear();
05214      steiner.reserve( numfaces );
05215 
05216      Point3D p3d;
05217      Face  *tface;
05218      vector<Vertex*> connect(3);
05219      for( size_t i = 0; i < numfaces; i++) {
05220           Face *f = quadmesh->getFaceAt(i);
05221           if( f->isActive() ) {
05222                int pos = f->concaveAt();
05223                if( pos >= 0) {
05224                     Vertex::mid_point( f->getNodeAt(pos), f->getNodeAt(pos+2), p3d );
05225                } else {
05226                     f->getAvgPos( p3d );
05227                }
05228                Vertex *v0 = Vertex::newObject();
05229                v0->setXYZCoords( p3d );
05230                trimesh->addNode(v0);
05231                steiner.push_back(v0);
05232                connect[0] = v0;
05233                for( int j = 0; j < 4; j++) {
05234                     connect[1] = f->getNodeAt(j);
05235                     connect[2] = f->getNodeAt(j+1);
05236                     tface  = Face::newObject();
05237                     tface->setNodes(connect);
05238                     trimesh->addFace( tface );
05239 //                  assert( tface->concaveAt() < 0);
05240                }
05241           }
05242      }
05243 
05244      return trimesh;
05245 }
05246 
05248 
05249 
05250 Mesh* Jaal::quad_to_tri2( Mesh *quadmesh )
05251 {
05252      if( quadmesh == NULL ) return NULL;
05253 
05254      Mesh *trimesh = new Mesh;
05255 
05256      size_t numnodes = quadmesh->getSize(0);
05257      size_t numfaces = quadmesh->getSize(2);
05258 
05259      trimesh->reserve( numnodes,  0);
05260      trimesh->reserve( 2*numfaces, 2);
05261 
05262      for( size_t i = 0; i < numnodes; i++) {
05263           Vertex *v = quadmesh->getNodeAt(i);
05264           if( v->isActive() ) trimesh->addNode(v);
05265      }
05266 
05267      Face  *tface;
05268      vector<Vertex*> connect(3);
05269      for( size_t i = 0; i < numfaces; i++) {
05270           Face *f = quadmesh->getFaceAt(i);
05271           if( f->isActive() ) {
05272                int pos = f->concaveAt();
05273                if( pos >= 0) {
05274                     connect[0] = f->getNodeAt(pos);
05275                     connect[1] = f->getNodeAt(pos+1);
05276                     connect[2] = f->getNodeAt(pos+2);
05277                     tface  = Face::newObject();
05278                     tface->setNodes(connect);
05279                     trimesh->addFace( tface );
05280 
05281                     connect[0] = f->getNodeAt(pos);
05282                     connect[1] = f->getNodeAt(pos+2);
05283                     connect[2] = f->getNodeAt(pos+3);
05284                     tface  = Face::newObject();
05285                     tface->setNodes(connect);
05286                     trimesh->addFace( tface );
05287                } else {
05288                     connect[0] = f->getNodeAt(0);
05289                     connect[1] = f->getNodeAt(1);
05290                     connect[2] = f->getNodeAt(2);
05291                     tface  = Face::newObject();
05292                     tface->setNodes(connect);
05293                     trimesh->addFace( tface );
05294 
05295                     connect[0] = f->getNodeAt(0);
05296                     connect[1] = f->getNodeAt(2);
05297                     connect[2] = f->getNodeAt(3);
05298                     tface  = Face::newObject();
05299                     tface->setNodes(connect);
05300                     trimesh->addFace( tface );
05301                }
05302           }
05303      }
05304      return trimesh;
05305 }
05306 
05308 
05309 void Jaal::advancing_front_triangle_cleanup( Mesh *mesh)
05310 {
05311      int lid;
05312      int relexist2 = mesh->build_relations(0, 2);
05313      int relexist0 = mesh->build_relations(0, 0);
05314 
05315      mesh->search_boundary();
05316 
05317      size_t numNodes = mesh->getSize(0);
05318      NodeSequence currlayer;
05319      NodeSet   nextlayer;
05320 
05321      for(size_t i = 0; i < numNodes; i++) {
05322           Vertex *v = mesh->getNodeAt(i);
05323           if( v->isBoundary() ) {
05324                v->setAttribute("Layer", 0);
05325                currlayer.push_back(v);
05326           } else
05327                v->setAttribute("Layer",std::numeric_limits< int >::max());
05328      }
05329 
05330      Jaal::MeshOptimization mopt;
05331      size_t nSize;
05332      mesh->make_consistently_oriented();
05333      mopt.shape_optimize(mesh);
05334 
05335      LaplaceNoWeight lw;
05336      LaplaceSmoothing lapsmooth(mesh);
05337      lapsmooth.setWeight(&lw);
05338      lapsmooth.setNumIterations(100);
05339 
05340      int curr_layer_id = 0, l0, l1;
05341      int nirregular0, nirregular1, numNeighs;
05342 
05343      Vertex *v0, *v1, *ov1, *ov2, *ov3;
05344      NodeSequence newnodes;
05345      FaceSequence newfaces, edgeneighs;
05346 
05347      while(1) {
05348           cout << " PROCESS LAYER : " << curr_layer_id << endl;
05349           nSize = currlayer.size();
05350           nirregular0 = 0;
05351           for( size_t i = 0; i < nSize; i++) {
05352                Vertex *v = currlayer[i];
05353                if( !v->isRemoved()  ) {
05354                     int curr_degree  = v->getNumRelations(2);
05355                     int ideal_degree = v->get_ideal_face_degree(3);
05356                     if( curr_degree != ideal_degree ) nirregular0++;
05357                }
05358           }
05359 
05360           FaceSequence v0faces;
05361           nSize = currlayer.size();
05362           for( size_t i = 0; i < nSize; i++) {
05363                Vertex *v = currlayer[i];
05364                if( !v->isRemoved() ) {
05365                     int curr_degree  = v->getNumRelations(2);
05366                     int ideal_degree = v->get_ideal_face_degree( 3 );
05367                     if( curr_degree < ideal_degree ) {
05368                          v->getRelations( v0faces );
05369                          numNeighs = v0faces.size();
05370                          for( int j = 0; j < numNeighs; j++) {
05371                               int pos = v0faces[j]->getPosOf(v);
05372                               v0  = v0faces[j]->getNodeAt( pos + 1 );
05373                               v1  = v0faces[j]->getNodeAt( pos + 2 );
05374                               Mesh::getRelations112( v0, v1, edgeneighs);
05375                               if( edgeneighs.size() == 2 ) {
05376                                    ov1 = Face::opposite_node(edgeneighs[0], v0, v1);
05377                                    ov2 = Face::opposite_node(edgeneighs[1], v0, v1);
05378                                    ov3 = NULL;
05379                                    assert( ov1 != ov2 );
05380                                    if( ov1 == v ) ov3 = ov2;
05381                                    if( ov2 == v ) ov3 = ov1;
05382                                    assert( ov3 );
05383                                    ov3->getAttribute("Layer", l0);
05384                                    if( l0 == std::numeric_limits< int >::max()) {
05385                                         int numSegments = ideal_degree - curr_degree + 1;
05386                                         mesh->refine_tri_edge( v0, v1, numSegments, newnodes, newfaces);
05387                                         for( size_t k = 0; k < newnodes.size(); k++)
05388                                              newnodes[k]->setAttribute("Layer", std::numeric_limits< double >::max());
05389                                         break;
05390                                    }
05391                               }
05392                          }
05393                          curr_degree  = v->getNumRelations(2);
05394                          if( curr_degree != ideal_degree ) {
05395                               cout << "Warning: low ideal vertex degree not achieved : " << v->getID() << endl;
05396                               set_layer_tag(mesh);
05397                               mesh->saveAs("tmp.off");
05398                               exit(0);
05399                          }
05400                     }
05401 
05402                     FaceSequence v1faces;
05403                     if( curr_degree > ideal_degree ) {
05404                          int excess_degree = curr_degree-ideal_degree;
05405                          for( int k = 0; k <  excess_degree; k++) {
05406                               v->getRelations( v1faces );
05407                               numNeighs = v1faces.size();
05408                               for( int j = 0; j < numNeighs; j++) {
05409                                    int pos = v1faces[j]->getPosOf(v);
05410                                    v0  = v1faces[j]->getNodeAt( pos + 1 );
05411                                    v1  = v1faces[j]->getNodeAt( pos + 2 );
05412                                    v0->getAttribute("Layer", l0);
05413                                    v1->getAttribute("Layer", l1);
05414                                    if( l0 == std::numeric_limits< int >::max() && l1 == std::numeric_limits< int >::max()) {
05415                                         Mesh::getRelations112( v0, v1, edgeneighs);
05416                                         assert( edgeneighs.size() == 2 );
05417                                         ov1 = Face::opposite_node(edgeneighs[0], v0, v1);
05418                                         ov2 = Face::opposite_node(edgeneighs[1], v0, v1);
05419                                         ov3 = NULL;
05420                                         assert( ov1 != ov2 );
05421                                         if( ov1 == v ) ov3 = ov2;
05422                                         if( ov2 == v ) ov3 = ov1;
05423                                         assert( ov3 );
05424                                         ov3->getAttribute("Layer", l0);
05425                                         if( l0 == std::numeric_limits< int >::max()) {
05426                                              mesh->collapse_tri_edge(v0, v1);
05427                                              break;
05428                                         }
05429                                    }
05430                               }
05431                          }
05432                          curr_degree  = v->getNumRelations(2);
05433                          if( curr_degree != ideal_degree ) {
05434                               cout << "Warning: high ideal vertex degree not achieved : " << v->getID() << endl;
05435                               set_layer_tag(mesh);
05436                               mesh->saveAs( "tmp.off");
05437                               exit(0);
05438                          }
05439                     }
05440                }
05441           }
05442 
05443           nSize = currlayer.size();
05444           nirregular1 = 0;
05445           for( size_t i = 0; i < nSize; i++) {
05446                Vertex *v = currlayer[i];
05447                if( !v->isRemoved()  ) {
05448                     int curr_degree  = v->getNumRelations(2);
05449                     int ideal_degree = v->get_ideal_face_degree( 3 );
05450                     if( curr_degree != ideal_degree ) nirregular1++;
05451                }
05452 
05453           }
05454           assert( nirregular1 <= nirregular0);
05455 //        lapsmooth.execute();
05456           mopt.shape_optimize(mesh);
05457 
05458           cout << "Layer : " << curr_layer_id << endl;
05459           cout << "# of Irregular nodes before swapping : " << nirregular0 << endl;
05460           cout << "# of Irregular nodes after swapping  : " << nirregular1 << endl;
05461 
05462           nextlayer.clear();
05463           nSize = currlayer.size();
05464           NodeSequence vneighs;
05465           for( size_t i = 0; i < nSize; i++) {
05466                Vertex *v = currlayer[i];
05467                v->getRelations( vneighs );
05468                for( size_t k = 0; k < vneighs.size(); k++) {
05469                     vneighs[k]->getAttribute("Layer", lid );
05470                     if( lid > curr_layer_id ) {
05471                          vneighs[k]->setAttribute("Layer", curr_layer_id+1 );
05472                          nextlayer.insert( vneighs[k] );
05473                     }
05474                }
05475           }
05476           if( nextlayer.empty() ) break;
05477 
05478           NodeSet::const_iterator it;
05479           currlayer.resize(nextlayer.size() );
05480           int index = 0;
05481           for( it = nextlayer.begin(); it != nextlayer.end(); ++it)
05482                currlayer[index++] = *it;
05483           curr_layer_id++;
05484      }
05485 
05486      vector<int>  less_than_ideal, more_than_ideal, total_ideal;
05487 
05488      int numLayers = curr_layer_id;
05489      less_than_ideal.resize( numLayers );
05490      more_than_ideal.resize( numLayers );
05491      total_ideal.resize( numLayers );
05492 
05493      for( int i = 0; i < numLayers; i++) {
05494           less_than_ideal[i] = 0;
05495           more_than_ideal[i] = 0;
05496           total_ideal[i] = 0;
05497      }
05498 
05499      numNodes = mesh->getSize(0);
05500      int final_irregular = 0;
05501      for( size_t i = 0; i < numNodes; i++) {
05502           Vertex *v = mesh->getNodeAt(i);
05503           if( !v->isRemoved()) {
05504                v->getAttribute("Layer", lid);
05505                int curr_degree  = v->getNumRelations(2);
05506                int ideal_degree = v->get_ideal_face_degree(3);
05507                if( curr_degree != ideal_degree ) {
05508                     final_irregular++;
05509                     if( curr_degree < ideal_degree) less_than_ideal[lid]++;
05510                     if( curr_degree > ideal_degree) more_than_ideal[lid]++;
05511                } else
05512                     total_ideal[lid]++;
05513           }
05514      }
05515 
05516      cout << " Layer   Less   More  Ideal " << endl;
05517      for( int i = 0; i < numLayers; i++)
05518           cout << i << setw(10) <<  less_than_ideal[i]
05519                << setw(10) <<  more_than_ideal[i]
05520                << setw(10) <<  total_ideal[i] << endl;
05521      cout << " Final # of irregular nodes : " << final_irregular << endl;
05522 
05523      mopt.shape_optimize(mesh);
05524      if (!relexist2) mesh->clear_relations(0, 2);
05525      if (!relexist0) mesh->clear_relations(0, 0);
05526 }
05527 
05528 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines