operation.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008,2009 Oliver Gloth                                     +
   7 // +                                                                      +
   8 // + enGrid is free software: you can redistribute it and/or modify       +
   9 // + it under the terms of the GNU General Public License as published by +
  10 // + the Free Software Foundation, either version 3 of the License, or    +
  11 // + (at your option) any later version.                                  +
  12 // +                                                                      +
  13 // + enGrid is distributed in the hope that it will be useful,            +
  14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  16 // + GNU General Public License for more details.                         +
  17 // +                                                                      +
  18 // + You should have received a copy of the GNU General Public License    +
  19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  20 // +                                                                      +
  21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  22 //
  23 #include "operation.h"
  24 #include "guimainwindow.h"
  25 #include "vtkTriangleFilter.h"
  26 #include "vtkInformation.h"
  27 #include "vtkInformationVector.h"
  28 #include "vtkObjectFactory.h"
  29 #include "vtkPointData.h"
  30 #include "vtkPolyData.h"
  31 #include "vtkPolygon.h"
  32 #include "vtkStreamingDemandDrivenPipeline.h"
  33 #include "vtkCellArray.h"
  34 #include "vtkCellData.h"
  35 #include "vtkCellLocator.h"
  36 #include "vtkFloatArray.h"
  37 #include "vtkMath.h"
  38 #include <vtkCharArray.h>
  39
  40 #include "geometrytools.h"
  41 using namespace GeometryTools;
  42
  43 #include <QApplication>
  44
  45 QSet<Operation*> Operation::garbage_operations;
  46
  47 void Operation::collectGarbage()
  48 {
  49   QSet<Operation*> delete_operations;
  50   foreach (Operation *op, garbage_operations) {
  51     if (!op->getThread().isRunning()) {
  52       delete_operations.insert(op);
  53       cout << "deleting Operation " << op << endl;
  54       delete op;
  55     };
  56   };
  57   foreach (Operation *op, delete_operations) {
  58     garbage_operations.remove(op);
  59   };
  60 };
  61
  62 Operation::Operation()
  63 {
  64   grid = NULL;
  65   gui = false;
  66   err = NULL;
  67   autoset = true;
  68 };
  69
  70 Operation::~Operation()
  71 {
  72   if (err) {
  73     err->display();
  74     delete err;
  75   };
  76 };
  77
  78 void Operation::del()
  79 {
  80   garbage_operations.insert(this);
  81 };
  82
  83 void OperationThread::run()
  84 {
  85   try {
  86     GuiMainWindow::lock();
  87     GuiMainWindow::pointer()->setBusy();
  88     op->operate();
  89   } catch (Error err) {
  90     op->err = new Error();
  91     *(op->err) = err;
  92   };
  93   GuiMainWindow::unlock();
  94   GuiMainWindow::pointer()->setIdle();
  95 };
  96
  97 void Operation::operator()()
  98 {
  99   if (gui) {
 100     if (GuiMainWindow::tryLock()) {
 101       checkGrid();
 102       thread.setOperation(this);
 103       GuiMainWindow::unlock();
 104       thread.start(QThread::LowPriority);
 105     } else {
 106       QMessageBox::warning(NULL, "not permitted", "Operation is not permitted while background process is running!");
 107     };
 108   } else {
 109     checkGrid();
 110     operate();
 111   };
 112   /*
 113   checkGrid();
 114   operate();
 115   */
 116 };
 117
 118 void Operation::setAllCells()
 119 {
 120   QVector<vtkIdType> all_cells;
 121   getAllCells(all_cells, grid);
 122   setCells(all_cells);
 123 };
 124
 125 void Operation::setAllVolumeCells()
 126 {
 127   QVector<vtkIdType> cells;
 128   getAllVolumeCells(cells, grid);
 129   setCells(cells);
 130 };
 131
 132 void Operation::setAllSurfaceCells()
 133 {
 134   QVector<vtkIdType> cells;
 135   getAllSurfaceCells(cells, grid);
 136   setCells(cells);
 137 };
 138
 139 void Operation::initMapping()
 140 {
 141   nodes_map.resize(nodes.size());
 142   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 143     nodes_map[i_nodes] = nodes[i_nodes];
 144   };
 145   cells_map.resize(cells.size());
 146   for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
 147     cells_map[i_cells] = cells[i_cells];
 148   };
 149 };
 150
 151 void Operation::checkGrid()
 152 {
 153   if (grid == NULL) {
 154     grid = GuiMainWindow::pointer()->getGrid();
 155   };
 156   if ((cells.size() == 0) && autoset) {
 157     setAllCells();
 158   };
 159 };
 160
 161 void Operation::updateActors()
 162 {
 163   mainWindow()->updateActors();
 164 };
 165
 166 GuiMainWindow* Operation::mainWindow()
 167 {
 168   return GuiMainWindow::pointer();
 169 };
 170
 171 void Operation::populateBoundaryCodes(QListWidget *lw)
 172 {
 173   QSet<int> bcs;
 174   mainWindow()->getAllBoundaryCodes(bcs);
 175   foreach(int bc, bcs) {
 176     QListWidgetItem *lwi = new QListWidgetItem(lw);
 177     lwi->setCheckState(Qt::Unchecked);
 178     QString text = "";
 179     QTextStream ts(&text);
 180     ts << bc;
 181     lwi->setText(text);
 182     lwi->setFlags(Qt::ItemIsUserCheckable | Qt::ItemIsEnabled);
 183   };
 184 };
 185
 186 stencil_t Operation::getStencil(vtkIdType id_cell1, int j1)
 187 {
 188   stencil_t S;
 189   S.valid = true;
 190   S.id_cell1 = id_cell1;
 191   if (c2c[_cells[id_cell1]][j1] != -1) {
 192     S.id_cell2 = cells[c2c[_cells[id_cell1]][j1]];
 193     if (grid->GetCellType(S.id_cell2) != VTK_TRIANGLE) {
 194       EG_BUG;
 195     };
 196     vtkIdType N1, N2, *pts1, *pts2;
 197     grid->GetCellPoints(S.id_cell1, N1, pts1);
 198     grid->GetCellPoints(S.id_cell2, N2, pts2);
 199     if      (j1 == 0) { S.p[0] = pts1[2]; S.p[1] = pts1[0]; S.p[3] = pts1[1]; }
 200     else if (j1 == 1) { S.p[0] = pts1[0]; S.p[1] = pts1[1]; S.p[3] = pts1[2]; }
 201     else if (j1 == 2) { S.p[0] = pts1[1]; S.p[1] = pts1[2]; S.p[3] = pts1[0]; };
 202     bool p2 = false;
 203     if (c2c[_cells[S.id_cell2]][0] != -1) {
 204       if (cells[c2c[_cells[S.id_cell2]][0]] == S.id_cell1) {
 205         S.p[2] = pts2[2];
 206         p2 = true;
 207       };
 208     };
 209     if (c2c[_cells[S.id_cell2]][1] != -1) {
 210       if (cells[c2c[_cells[S.id_cell2]][1]] == S.id_cell1) {
 211         S.p[2] = pts2[0];
 212         p2 = true;
 213       };
 214     };
 215     if (c2c[_cells[S.id_cell2]][2] != -1) {
 216       if (cells[c2c[_cells[S.id_cell2]][2]] == S.id_cell1) {
 217         S.p[2] = pts2[1];
 218         p2 = true;
 219       };
 220     };
 221     if (!p2) {
 222       EG_BUG;
 223     };
 224   } else {
 225     S.valid = false;
 226     S.id_cell2 = -1;
 227     vtkIdType N1, *pts1;
 228     grid->GetCellPoints(S.id_cell1, N1, pts1);
 229     if      (j1 == 0) { S.p[0] = pts1[2]; S.p[1] = pts1[0]; S.p[3] = pts1[1]; }
 230     else if (j1 == 1) { S.p[0] = pts1[0]; S.p[1] = pts1[1]; S.p[3] = pts1[2]; }
 231     else if (j1 == 2) { S.p[0] = pts1[1]; S.p[1] = pts1[2]; S.p[3] = pts1[0]; };
 232   };
 233   return S;
 234 };
 235
 236 ostream& operator<<(ostream &out, stencil_t S)
 237 {
 238   out<<"S.id_cell1="<<S.id_cell1<<" ";
 239   out<<"S.id_cell2="<<S.id_cell2<<" ";
 240   out<<"S.valid="<<S.valid<<" ";
 241   out<<"[";
 242   for(int i=0;i<4;i++){
 243     out<<S.p[i];
 244     if(i!=3) out<<",";
 245   }
 246   out<<"]";
 247   return(out);
 248 }
 249
 250 //////////////////////////////////////////////
 251 double CurrentVertexAvgDist(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 252 {
 253   double total_dist=0;
 254   double avg_dist=0;
 255   int N=n2n[a_vertex].size();
 256   vec3_t C;
 257   a_grid->GetPoint(a_vertex, C.data());
 258   foreach(int i,n2n[a_vertex])
 259   {
 260     vec3_t M;
 261     a_grid->GetPoint(i, M.data());
 262     total_dist+=(M-C).abs();
 263   }
 264   avg_dist=total_dist/(double)N;
 265   return(avg_dist);
 266 }
 267
 268 double CurrentMeshDensity(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 269 {
 270   double total_dist=0;
 271   double avg_dist=0;
 272   int N=n2n[a_vertex].size();
 273   vec3_t C;
 274   a_grid->GetPoint(a_vertex, C.data());
 275   foreach(int i,n2n[a_vertex])
 276   {
 277     vec3_t M;
 278     a_grid->GetPoint(i, M.data());
 279     total_dist+=(M-C).abs();
 280   }
 281   avg_dist=total_dist/(double)N;
 282   double avg_density=1./avg_dist;
 283   return(avg_density);
 284 }
 285
 286 double DesiredVertexAvgDist(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 287 {
 288   double total_dist=0;
 289   double avg_dist=0;
 290   EG_VTKDCN(vtkDoubleArray, node_meshdensity, a_grid, "node_meshdensity");
 291   int N=n2n[a_vertex].size();
 292   foreach(int i,n2n[a_vertex])
 293   {
 294     total_dist+=1./node_meshdensity->GetValue(i);
 295   }
 296   avg_dist=total_dist/(double)N;
 297   return(avg_dist);
 298 }
 299
 300 double DesiredMeshDensity(vtkIdType a_vertex,QVector< QSet< int > > &n2n,vtkUnstructuredGrid *a_grid)
 301 {
 302   double total_density=0;
 303   double avg_density=0;
 304   EG_VTKDCN(vtkDoubleArray, node_meshdensity, a_grid, "node_meshdensity");
 305   int N=n2n[a_vertex].size();
 306   foreach(int i,n2n[a_vertex])
 307   {
 308     total_density+=node_meshdensity->GetValue(i);
 309   }
 310   avg_density=total_density/(double)N;
 311   return(avg_density);
 312 }
 313
 314 ///////////////////////////////////////////
 315 vtkIdType Operation::getClosestNode(vtkIdType a_id_node,vtkUnstructuredGrid* a_grid)
 316 {
 317   vec3_t C;
 318   a_grid->GetPoint(a_id_node,C.data());
 319   vtkIdType id_minlen=-1;
 320   double minlen=-1;
 321   foreach(vtkIdType neighbour,n2n[a_id_node])
 322   {
 323     vec3_t M;
 324     a_grid->GetPoint(neighbour,M.data());
 325     double len=(M-C).abs();
 326     if(minlen<0 or len<minlen)
 327     {
 328       minlen=len;
 329       id_minlen=neighbour;
 330     }
 331   }
 332   return(id_minlen);
 333 }
 334
 335 vtkIdType Operation::getFarthestNode(vtkIdType a_id_node,vtkUnstructuredGrid* a_grid)
 336 {
 337   vec3_t C;
 338   a_grid->GetPoint(a_id_node,C.data());
 339   vtkIdType id_maxlen=-1;
 340   double maxlen=-1;
 341   foreach(vtkIdType neighbour,n2n[a_id_node])
 342   {
 343     vec3_t M;
 344     a_grid->GetPoint(neighbour,M.data());
 345     double len=(M-C).abs();
 346     if(maxlen<0 or len>maxlen)
 347     {
 348       maxlen=len;
 349       id_maxlen=neighbour;
 350     }
 351   }
 352   return(id_maxlen);
 353 }
 354
 355 bool Operation::SwapCells(vtkUnstructuredGrid* a_grid, stencil_t S)
 356 {
 357   bool swap = false;
 358   if (S.valid) {
 359     vec3_t x3[4], x3_0(0,0,0);
 360     vec2_t x[4];
 361     for (int k = 0; k < 4; ++k) {
 362       a_grid->GetPoints()->GetPoint(S.p[k], x3[k].data());
 363       x3_0 += x3[k];
 364     };
 365     vec3_t n1 = triNormal(x3[0], x3[1], x3[3]);
 366     vec3_t n2 = triNormal(x3[1], x3[2], x3[3]);
 367     n1.normalise();
 368     n2.normalise();
 369     if ( (n1*n2) > 0.8) {
 370       vec3_t n = n1 + n2;
 371       n.normalise();
 372       vec3_t ex = orthogonalVector(n);
 373       vec3_t ey = ex.cross(n);
 374       for (int k = 0; k < 4; ++k) {
 375         x[k] = vec2_t(x3[k]*ex, x3[k]*ey);
 376       };
 377       vec2_t r1, r2, r3, u1, u2, u3;
 378       r1 = 0.5*(x[0] + x[1]); u1 = turnLeft(x[1] - x[0]);
 379       r2 = 0.5*(x[1] + x[2]); u2 = turnLeft(x[2] - x[1]);
 380       r3 = 0.5*(x[1] + x[3]); u3 = turnLeft(x[3] - x[1]);
 381       double k, l;
 382       vec2_t xm1, xm2;
 383       bool ok = true;
 384       if (intersection(k, l, r1, u1, r3, u3)) {
 385         xm1 = r1 + k*u1;
 386         if (intersection(k, l, r2, u2, r3, u3)) {
 387           xm2 = r2 + k*u2;
 388         } else {
 389           ok = false;
 390         };
 391       } else {
 392         ok = false;
 393         swap = true;
 394       };
 395       if (ok) {
 396         if ((xm1 - x[2]).abs() < (xm1 - x[0]).abs()) {
 397           swap = true;
 398         };
 399         if ((xm2 - x[0]).abs() < (xm2 - x[2]).abs()) {
 400           swap = true;
 401         };
 402       };
 403     };
 404   };
 405   if (swap) {
 406     vtkIdType new_pts1[3], new_pts2[3];
 407     new_pts1[0] = S.p[1];
 408     new_pts1[1] = S.p[2];
 409     new_pts1[2] = S.p[0];
 410     new_pts2[0] = S.p[2];
 411     new_pts2[1] = S.p[3];
 412     new_pts2[2] = S.p[0];
 413     a_grid->ReplaceCell(S.id_cell1, 3, new_pts1);
 414     a_grid->ReplaceCell(S.id_cell2, 3, new_pts2);
 415   };
 416   return(swap);
 417 }
 418
 419 void Operation::quad2triangle(vtkUnstructuredGrid* src,vtkIdType quadcell)
 420 {
 421   vtkIdType type_cell = grid->GetCellType(quadcell);
 422   cout<<"It's a "<<type_cell<<endl;
 423   if(type_cell==VTK_QUAD)
 424   {
 425     cout_grid(cout,src,true,true,true,true);
 426     EG_VTKSP(vtkUnstructuredGrid, dst);
 427     //src grid info
 428     int N_points=src->GetNumberOfPoints();
 429     int N_cells=src->GetNumberOfCells();
 430     allocateGrid(dst,N_cells+1,N_points);
 431
 432     for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
 433       vec3_t x;
 434       src->GetPoints()->GetPoint(id_node, x.data());
 435       dst->GetPoints()->SetPoint(id_node, x.data());
 436       copyNodeData(src, id_node, dst, id_node);
 437     };
 438     for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 439       vtkIdType N_pts, *pts;
 440       vtkIdType type_cell = src->GetCellType(id_cell);
 441       src->GetCellPoints(id_cell, N_pts, pts);
 442       vtkIdType id_new_cell;
 443       vtkIdType id_new_cell1;
 444       vtkIdType id_new_cell2;
 445       if(id_cell!=quadcell)
 446       {
 447         id_new_cell = dst->InsertNextCell(type_cell, N_pts, pts);
 448         copyCellData(src, id_cell, dst, id_new_cell);
 449       }
 450       else
 451       {
 452         vtkIdType triangle1[3];
 453         vtkIdType triangle2[3];
 454         triangle1[0]=pts[1];
 455         triangle1[1]=pts[3];
 456         triangle1[2]=pts[0];
 457         triangle2[0]=pts[3];
 458         triangle2[1]=pts[1];
 459         triangle2[2]=pts[2];
 460         id_new_cell1 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle1);
 461         copyCellData(src, id_cell, dst, id_new_cell1);
 462         id_new_cell2 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle2);
 463         copyCellData(src, id_cell, dst, id_new_cell2);
 464         stencil_t S;
 465         S.id_cell1=id_new_cell1;
 466         S.id_cell2=id_new_cell2;
 467         S.p[0]=pts[0];
 468         S.p[1]=pts[1];
 469         S.p[2]=pts[2];
 470         S.p[3]=pts[3];
 471         S.valid=true;
 472         SwapCells(dst,S);
 473       }
 474     };
 475     cout_grid(cout,dst,true,true,true,true);
 476     makeCopy(dst, src);
 477   }//end of if quad
 478 }
 479
 480 void Operation::quad2triangle(vtkUnstructuredGrid* src,vtkIdType quadcell,vtkIdType MovingPoint)
 481 {
 482   vtkIdType type_cell = grid->GetCellType(quadcell);
 483   cout<<"It's a "<<type_cell<<endl;
 484   if(type_cell==VTK_QUAD)
 485   {
 486     cout_grid(cout,src,true,true,true,true);
 487     EG_VTKSP(vtkUnstructuredGrid, dst);
 488     //src grid info
 489     int N_points=src->GetNumberOfPoints();
 490     int N_cells=src->GetNumberOfCells();
 491     allocateGrid(dst,N_cells+1,N_points);
 492
 493     for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
 494       vec3_t x;
 495       src->GetPoints()->GetPoint(id_node, x.data());
 496       dst->GetPoints()->SetPoint(id_node, x.data());
 497       copyNodeData(src, id_node, dst, id_node);
 498     };
 499     for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 500       vtkIdType N_pts, *pts;
 501       src->GetCellPoints(id_cell, N_pts, pts);
 502       vtkIdType type_cell = src->GetCellType(id_cell);
 503       vtkIdType id_new_cell;
 504       vtkIdType id_new_cell1;
 505       vtkIdType id_new_cell2;
 506       if(id_cell!=quadcell)
 507       {
 508         id_new_cell = dst->InsertNextCell(type_cell, N_pts, pts);
 509         copyCellData(src, id_cell, dst, id_new_cell);
 510       }
 511       else
 512       {
 513         vtkIdType triangle1[3];
 514         vtkIdType triangle2[3];
 515         if(MovingPoint==pts[0] || MovingPoint==pts[2])
 516         {
 517           triangle1[0]=pts[1];
 518           triangle1[1]=pts[3];
 519           triangle1[2]=pts[0];
 520           triangle2[0]=pts[3];
 521           triangle2[1]=pts[1];
 522           triangle2[2]=pts[2];
 523         }
 524         else
 525         {
 526           triangle1[0]=pts[2];
 527           triangle1[1]=pts[0];
 528           triangle1[2]=pts[1];
 529           triangle2[0]=pts[0];
 530           triangle2[1]=pts[2];
 531           triangle2[2]=pts[3];
 532         }
 533         id_new_cell1 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle1);
 534         copyCellData(src, id_cell, dst, id_new_cell1);
 535         id_new_cell2 = dst->InsertNextCell(VTK_TRIANGLE, 3, triangle2);
 536         copyCellData(src, id_cell, dst, id_new_cell2);
 537       }
 538     };
 539     cout_grid(cout,dst,true,true,true,true);
 540     makeCopy(dst, src);
 541   }//end of if quad
 542 }
 543
 544 int Operation::NumberOfCommonPoints(vtkIdType node1, vtkIdType node2, bool& IsTetra)
 545 {
 546 //   QVector< QSet< int > >     n2n
 547   QSet <int> node1_neighbours=n2n[node1];
 548   QSet <int> node2_neighbours=n2n[node2];
 549   QSet <int> intersection=node1_neighbours.intersect(node2_neighbours);
 550   int N=intersection.size();
 551   IsTetra=false;
 552   if(N==2)
 553   {
 554     QSet<int>::const_iterator p1=intersection.begin();
 555     QSet<int>::const_iterator p2=p1+1;
 556     vtkIdType intersection1=_nodes[*p1];
 557     vtkIdType intersection2=_nodes[*p2];
 558     if(n2n[intersection1].contains(intersection2))//if there's an edge between intersection1 and intersection2
 559     {
 560       //check if (node1,intersection1,intersection2) and (node2,intersection1,intersection2) are defined as cells!
 561   //     QVector< QSet< int > >         n2c
 562       QSet< int > S1=n2c[intersection1];
 563       QSet< int > S2=n2c[intersection2];
 564       QSet< int > Si=S1.intersect(S2);
 565       int counter=0;
 566       foreach(vtkIdType C,Si){
 567         vtkIdType N_pts, *pts;
 568         grid->GetCellPoints(C, N_pts, pts);
 569         for(int i=0;i<N_pts;i++)
 570         {
 571           if(pts[i]==node1 || pts[i]==node2) counter++;
 572         }
 573       }
 574       if(counter>=2) IsTetra=true;
 575     }
 576   }
 577   return(N);
 578 }
 579
 580 // vtkIdType Operation::FindSnapPoint(vtkUnstructuredGrid *src, vtkIdType DeadNode)
 581 // {
 582 //   return(0);
 583 // }
 584
 585 bool Operation::DeletePoint(vtkUnstructuredGrid *src, vtkIdType DeadNode)
 586 {
 587   EG_VTKSP(vtkUnstructuredGrid, dst);
 588
 589     //src grid info
 590   int N_points=src->GetNumberOfPoints();
 591   int N_cells=src->GetNumberOfCells();
 592   int N_newpoints=-1;
 593   int N_newcells=0;
 594
 595   if(DeadNode<0 || DeadNode>=N_points)
 596   {
 597     cout<<"Warning: Point out of range: DeadNode="<<DeadNode<<" N_points="<<N_points<<endl;
 598     return(false);
 599   }
 600
 601   QSet <vtkIdType> DeadCells;
 602   QSet <vtkIdType> MutatedCells;
 603   QSet <vtkIdType> MutilatedCells;
 604
 605   vtkIdType SnapPoint=-1;
 606   //Find closest point to DeadNode
 607 //   vtkIdType SnapPoint = getClosestNode(DeadNode,src);//DeadNode moves to SnapPoint
 608
 609   foreach(vtkIdType PSP, n2n[DeadNode])
 610   {
 611     bool IsValidSnapPoint=true;
 612
 613     cout<<"====>PSP="<<PSP<<endl;
 614     bool IsTetra=true;
 615     if(NumberOfCommonPoints(DeadNode,PSP,IsTetra)>2)//common point check
 616     {
 617       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 618       IsValidSnapPoint=false;
 619     }
 620     if(IsTetra)//tetra check
 621     {
 622       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 623       IsValidSnapPoint=false;
 624     }
 625
 626     //count number of points and cells to remove + analyse cell transformations
 627     N_newpoints=-1;
 628     N_newcells=0;
 629     DeadCells.clear();
 630     MutatedCells.clear();
 631     MutilatedCells.clear();
 632     foreach(vtkIdType C, n2c[DeadNode])//loop through potentially dead cells
 633     {
 634       cout<<"C="<<C<<endl;
 635       //get points around cell
 636       vtkIdType N_pts, *pts;
 637       src->GetCellPoints(C, N_pts, pts);
 638
 639       bool ContainsSnapPoint=false;
 640       bool invincible=false;
 641       for(int i=0;i<N_pts;i++)
 642       {
 643         cout<<"pts["<<i<<"]="<<pts[i]<<" and PSP="<<PSP<<endl;
 644         if(pts[i]==PSP) {ContainsSnapPoint=true;}
 645         if(pts[i]!=DeadNode && pts[i]!=PSP &&  n2c[pts[i]].size()<=1) invincible=true;
 646       }
 647       if(ContainsSnapPoint)
 648       {
 649         if(N_pts==3)//dead cell
 650         {
 651           if(invincible)//Check that empty lines aren't left behind when a cell is killed
 652           {
 653             cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 654             IsValidSnapPoint=false;
 655           }
 656           else
 657           {
 658             DeadCells.insert(C);
 659             N_newcells-=1;
 660             cout<<"cell "<<C<<" has been pwned!"<<endl;
 661           }
 662         }
 663   /*      else if(N_pts==4)//mutilated cell
 664         {
 665           MutilatedCells.insert(C);
 666           cout<<"cell "<<C<<" has lost a limb!"<<endl;
 667         }*/
 668         else
 669         {
 670           cout<<"RED ALERT: Xenomorph detected!"<<endl;
 671           EG_BUG;
 672         }
 673       }
 674       else
 675       {
 676         vtkIdType src_N_pts, *src_pts;
 677         src->GetCellPoints(C, src_N_pts, src_pts);
 678
 679         if(src_N_pts!=3)
 680         {
 681           cout<<"RED ALERT: Xenomorph detected!"<<endl;
 682           EG_BUG;
 683         }
 684
 685         vtkIdType OldTriangle[3];
 686         vtkIdType NewTriangle[3];
 687
 688         for(int i=0;i<src_N_pts;i++)
 689         {
 690           OldTriangle[i]=src_pts[i];
 691           NewTriangle[i]=( (src_pts[i]==DeadNode) ? PSP : src_pts[i] );
 692         }
 693         vec3_t Old_N= triNormal(src, OldTriangle[0], OldTriangle[1], OldTriangle[2]);
 694         vec3_t New_N= triNormal(src, NewTriangle[0], NewTriangle[1], NewTriangle[2]);
 695         double OldArea=Old_N.abs();
 696         double NewArea=New_N.abs();
 697         double scal=Old_N*New_N;
 698         double cross=(Old_N.cross(New_N)).abs();//double-cross on Nar Shadaa B-)
 699
 700         cout<<"OldArea="<<OldArea<<endl;
 701         cout<<"NewArea="<<NewArea<<endl;
 702         cout<<"scal="<<scal<<endl;
 703         cout<<"cross="<<cross<<endl;
 704
 705         if(Old_N*New_N<Old_N*Old_N*1./100.)//area + inversion check
 706         {
 707           cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 708           IsValidSnapPoint=false;
 709         }
 710
 711   /*      if(NewArea<OldArea*1./100.)
 712         {
 713           cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 714           IsValidSnapPoint=false;
 715         }
 716
 717         if(abs(cross)>10e-4)
 718         {
 719           cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 720           IsValidSnapPoint=false;
 721         }*/
 722
 723         //mutated cell
 724         MutatedCells.insert(C);
 725         cout<<"cell "<<C<<" has been infected!"<<endl;
 726       }
 727     }
 728
 729     if(N_cells+N_newcells<=0)//survivor check
 730     {
 731       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 732       IsValidSnapPoint=false;
 733     }
 734 /*    if(EmptyVolume(DeadNode,PSP))//simplified volume check
 735     {
 736       cout<<"Sorry, but you are not allowed to move point "<<DeadNode<<" to point "<<PSP<<"."<<endl;
 737       IsValidSnapPoint=false;
 738     }*/
 739     if(IsValidSnapPoint) {SnapPoint=PSP; break;}
 740   }//end of loop through potential SnapPoints
 741
 742   cout<<"===>SNAPPOINT="<<SnapPoint<<endl;
 743   if(SnapPoint<0) {cout<<"Sorry no possible SnapPoint found."<<endl; return(false);}
 744
 745   //allocate
 746   cout<<"N_points="<<N_points<<endl;
 747   cout<<"N_cells="<<N_cells<<endl;
 748   cout<<"N_newpoints="<<N_newpoints<<endl;
 749   cout<<"N_newcells="<<N_newcells<<endl;
 750   allocateGrid(dst,N_cells+N_newcells,N_points+N_newpoints);
 751
 752   //vector used to redefine the new point IDs
 753   QVector <vtkIdType> OffSet(N_points);
 754
 755   //copy undead points
 756   vtkIdType dst_id_node=0;
 757   for (vtkIdType src_id_node = 0; src_id_node < N_points; src_id_node++) {//loop through src points
 758     if(src_id_node!=DeadNode)//if the node isn't dead, copy it
 759     {
 760       vec3_t x;
 761       src->GetPoints()->GetPoint(src_id_node, x.data());
 762       dst->GetPoints()->SetPoint(dst_id_node, x.data());
 763       copyNodeData(src, src_id_node, dst, dst_id_node);
 764       OffSet[src_id_node]=src_id_node-dst_id_node;
 765       dst_id_node++;
 766     }
 767   };
 768
 769   cout<<"DeadCells="<<DeadCells<<endl;
 770   cout<<"MutatedCells="<<MutatedCells<<endl;
 771   cout<<"MutilatedCells="<<MutilatedCells<<endl;
 772
 773   //Copy undead cells
 774   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {//loop through src cells
 775     if(!DeadCells.contains(id_cell))//if the cell isn't dead
 776     {
 777       vtkIdType src_N_pts, *src_pts;
 778       vtkIdType dst_N_pts, *dst_pts;
 779       src->GetCellPoints(id_cell, src_N_pts, src_pts);
 780
 781       vtkIdType type_cell = src->GetCellType(id_cell);
 782       cout<<"-->id_cell="<<id_cell<<endl;
 783       for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
 784 //       src->GetCellPoints(id_cell, dst_N_pts, dst_pts);
 785       dst_N_pts=src_N_pts;
 786       dst_pts=new vtkIdType[dst_N_pts];
 787       if(MutatedCells.contains(id_cell))//mutated cell
 788       {
 789         cout<<"processing mutated cell "<<id_cell<<endl;
 790         for(int i=0;i<src_N_pts;i++)
 791         {
 792           if(src_pts[i]==DeadNode) {
 793             cout<<"SnapPoint="<<SnapPoint<<endl;
 794             cout<<"OffSet[SnapPoint]="<<OffSet[SnapPoint]<<endl;
 795             cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
 796             dst_pts[i]=SnapPoint-OffSet[SnapPoint];
 797           }
 798           else dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
 799         }
 800         cout<<"--->dst_pts:"<<endl;
 801         for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
 802
 803       }
 804       else if(MutilatedCells.contains(id_cell))//mutilated cell
 805       {
 806         cout<<"processing mutilated cell "<<id_cell<<endl;
 807
 808         if(type_cell==VTK_QUAD) {
 809           type_cell=VTK_TRIANGLE;
 810           dst_N_pts-=1;
 811         }
 812         else {cout<<"FATAL ERROR: Unknown mutilated cell detected! It is not a quad! Potential xenomorph infestation!"<<endl;EG_BUG;}
 813         //merge points
 814         int j=0;
 815         for(int i=0;i<src_N_pts;i++)
 816         {
 817           if(src_pts[i]==SnapPoint) { dst_pts[j]=SnapPoint-OffSet[SnapPoint];j++; }//SnapPoint
 818           else if(src_pts[i]!=DeadNode) { dst_pts[j]=src_pts[i]-OffSet[src_pts[i]];j++; }//pre-snap/dead + post-snap/dead
 819           //do nothing in case of DeadNode
 820         }
 821       }
 822       else//normal cell
 823       {
 824         cout<<"processing normal cell "<<id_cell<<endl;
 825         for(int i=0;i<src_N_pts;i++)
 826         {
 827           dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
 828         }
 829       }
 830       //copy the cell
 831       vtkIdType id_new_cell = dst->InsertNextCell(type_cell, dst_N_pts, dst_pts);
 832       copyCellData(src, id_cell, dst, id_new_cell);
 833       cout<<"===Copying cell "<<id_cell<<" to "<<id_new_cell<<"==="<<endl;
 834       cout<<"src_pts:"<<endl;
 835       for(int i=0;i<src_N_pts;i++) cout<<"src_pts["<<i<<"]="<<src_pts[i]<<endl;
 836       cout<<"dst_pts:"<<endl;
 837       for(int i=0;i<dst_N_pts;i++) cout<<"dst_pts["<<i<<"]="<<dst_pts[i]<<endl;
 838       cout<<"OffSet="<<OffSet<<endl;
 839       cout<<"===Copying cell end==="<<endl;
 840       delete dst_pts;
 841     }
 842   };
 843 //   cout_grid(cout,dst,true,true,true,true);
 844   makeCopy(dst, src);
 845   return(true);
 846 }
 847
 848 bool Operation::EmptyVolume(vtkIdType DeadNode, vtkIdType PSP)
 849 {
 850   c2c[DeadNode];
 851   c2c[PSP];
 852   return(true);
 853 }
 854
 855 vec3_t Operation::GetCenter(vtkIdType cellId, double& R)
 856 {
 857   vtkIdType *pts, Npts;
 858   grid->GetCellPoints(cellId, Npts, pts);
 859
 860   vec3_t x(0,0,0);
 861   for (vtkIdType i = 0; i < Npts; ++i) {
 862     vec3_t xp;
 863     grid->GetPoints()->GetPoint(pts[i], xp.data());
 864     x += double(1)/Npts * xp;
 865   };
 866
 867   R = 1e99;
 868   for (vtkIdType i = 0; i < Npts; ++i) {
 869     vec3_t xp;
 870     grid->GetPoints()->GetPoint(pts[i], xp.data());
 871     R = min(R, 0.25*(xp-x).abs());
 872   };
 873
 874   return(x);
 875 }
 876
 877 bool Operation::getNeighbours(vtkIdType Boss, vtkIdType& Peon1, vtkIdType& Peon2, int BC)
 878 {
 879   QVector <vtkIdType> Peons;
 880
 881   QSet <int> S1=n2c[Boss];
 882   cout<<"S1="<<S1<<endl;
 883   foreach(vtkIdType PN,n2n[Boss])
 884   {
 885     cout<<"PN="<<PN<<endl;
 886     QSet <int> S2=n2c[PN];
 887     cout<<"S2="<<S2<<endl;
 888     QSet <int> Si=S2.intersect(S1);
 889     cout<<"PN="<<PN<<" Si="<<Si<<endl;
 890     if(Si.size()<2)//only one common cell
 891     {
 892       Peons.push_back(PN);
 893     }
 894     else
 895     {
 896       QSet <int> bc_set;
 897       foreach(vtkIdType C,Si)
 898       {
 899         EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 900         int bc=cell_code->GetValue(C);
 901         cout<<"C="<<C<<" bc="<<bc<<endl;
 902         bc_set.insert(bc);
 903       }
 904       if(bc_set.size()>1)//2 different boundary codes
 905       {
 906         Peons.push_back(PN);
 907       }
 908     }
 909   }
 910   if(Peons.size()==2)
 911   {
 912     Peon1=Peons[0];
 913     Peon2=Peons[1];
 914     return(true);
 915   }
 916   else
 917   {
 918     cout<<"FATAL ERROR: number of neighbours != 2"<<endl;
 919     EG_BUG;
 920   }
 921   return(false);
 922 }
 923
 924 int Operation::UpdateMeshDensity()
 925 {
 926   cout<<"===UpdateMeshDensity START==="<<endl;
 927
 928   getAllSurfaceCells(cells,grid);
 929   EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 930   EG_VTKDCN(vtkDoubleArray, node_meshdensity, grid, "node_meshdensity");
 931   getNodesFromCells(cells, nodes, grid);
 932   setGrid(grid);
 933   setCells(cells);
 934
 935   if(DebugLevel>5) cout<<"cells.size()="<<cells.size()<<endl;
 936
 937   EG_VTKDCN(vtkDoubleArray, node_meshdensity_current, grid, "node_meshdensity_current");
 938   foreach(vtkIdType node,nodes)
 939   {
 940     double L=CurrentVertexAvgDist(node,n2n,grid);
 941     double D=1./L;
 942     node_meshdensity_current->SetValue(node, D);
 943   }
 944   cout<<"===UpdateMeshDensity END==="<<endl;
 945   return(0);
 946 }
 947
 948 // Special structure for marking vertices
 949 typedef struct _vtkMeshVertex
 950 {
 951   char      type;
 952   vtkIdList *edges; // connected edges (list of connected point ids)
 953 } vtkMeshVertex, *vtkMeshVertexPtr;
 954
 955 int Operation::UpdateNodeType()
 956 {
 957 //   cout<<"===UpdateNodeType START==="<<endl;
 958
 959   getAllSurfaceCells(cells,grid);
 960   if(DebugLevel>5) cout<<"cells.size()="<<cells.size()<<endl;
 961
 962   EG_VTKSP(vtkPolyData, pdata);
 963   //   addToPolyData(m_SelectedCells, pdata, grid);
 964   addToPolyData(cells, pdata, grid);
 965
 966   vtkPolyData* input=pdata;
 967
 968   vtkPolyData *source = 0;
 969
 970   vtkIdType numPts, numCells, i, numPolys, numStrips;
 971   int j, k;
 972   vtkIdType npts = 0;
 973   vtkIdType *pts = 0;
 974   vtkIdType p1, p2;
 975   double x[3], y[3], deltaX[3], xNew[3], conv, maxDist, dist, factor;
 976   double x1[3], x2[3], x3[3], l1[3], l2[3];
 977   double CosFeatureAngle; //Cosine of angle between adjacent polys
 978   double CosEdgeAngle; // Cosine of angle between adjacent edges
 979   double closestPt[3], dist2, *w = NULL;
 980   int iterationNumber, abortExecute;
 981   vtkIdType numSimple=0, numBEdges=0, numFixed=0, numFEdges=0;
 982   vtkPolyData *inMesh, *Mesh;
 983   vtkPoints *inPts;
 984   vtkTriangleFilter *toTris=NULL;
 985   vtkCellArray *inVerts, *inLines, *inPolys, *inStrips;
 986   vtkPoints *newPts;
 987   vtkMeshVertexPtr Verts;
 988   vtkCellLocator *cellLocator=NULL;
 989
 990     // Check input
 991     //
 992   numPts=input->GetNumberOfPoints();
 993   numCells=input->GetNumberOfCells();
 994   if (numPts < 1 || numCells < 1)
 995   {
 996     cout<<"No data to smooth!"<<endl;
 997     return 1;
 998   }
 999
1000   CosFeatureAngle =
1001     cos((double) vtkMath::DegreesToRadians() * this->FeatureAngle);
1002   CosEdgeAngle = cos((double) vtkMath::DegreesToRadians() * this->EdgeAngle);
1003
1004   if(DebugLevel>5) {
1005     cout<<"Smoothing " << numPts << " vertices, " << numCells
1006       << " cells with:\n"
1007       << "\tConvergence= " << this->Convergence << "\n"
1008       << "\tIterations= " << this->NumberOfIterations << "\n"
1009       << "\tRelaxation Factor= " << this->RelaxationFactor << "\n"
1010       << "\tEdge Angle= " << this->EdgeAngle << "\n"
1011       << "\tBoundary Smoothing " << (this->BoundarySmoothing ? "On\n" : "Off\n")
1012       << "\tFeature Edge Smoothing " << (this->FeatureEdgeSmoothing ? "On\n" : "Off\n")
1013       << "\tError Scalars " << (this->GenerateErrorScalars ? "On\n" : "Off\n")
1014       << "\tError Vectors " << (this->GenerateErrorVectors ? "On\n" : "Off\n")<<endl;
1015   }
1016     // Peform topological analysis. What we're gonna do is build a connectivity
1017     // array of connected vertices. The outcome will be one of three
1018     // classifications for a vertex: VTK_SIMPLE_VERTEX, VTK_FIXED_VERTEX. or
1019     // VTK_EDGE_VERTEX. Simple vertices are smoothed using all connected
1020     // vertices. FIXED vertices are never smoothed. Edge vertices are smoothed
1021     // using a subset of the attached vertices.
1022     //
1023   if(DebugLevel>5) cout<<"===>Analyze topology==="<<endl;
1024   if(DebugLevel>5) cout<<"Analyzing topology..."<<endl;
1025   if(DebugLevel>5) cout<<"0:numPts="<<numPts<<endl;
1026   Verts = new vtkMeshVertex[numPts];
1027   for (i=0; i<numPts; i++)
1028   {
1029     if(DebugLevel>5) cout<<"0:VTK_SIMPLE_VERTEX"<<endl;
1030     Verts[i].type = VTK_SIMPLE_VERTEX; //can smooth
1031     Verts[i].edges = NULL;
1032   }
1033
1034   inPts = input->GetPoints();
1035   conv = this->Convergence * input->GetLength();
1036
1037     // check vertices first. Vertices are never smoothed_--------------
1038   for (inVerts=input->GetVerts(), inVerts->InitTraversal();
1039        inVerts->GetNextCell(npts,pts); )
1040   {
1041     for (j=0; j<npts; j++)
1042     {
1043       if(DebugLevel>5) cout<<"pts[j]="<<pts[j]<<"->vertices:VTK_FIXED_VERTEX 0"<<endl;
1044       Verts[pts[j]].type = VTK_FIXED_VERTEX;
1045     }
1046   }
1047
1048   if(DebugLevel>5) cout<<"==check lines=="<<endl;
1049     // now check lines. Only manifold lines can be smoothed------------
1050   for (inLines=input->GetLines(), inLines->InitTraversal();
1051        inLines->GetNextCell(npts,pts); )
1052   {
1053     for (j=0; j<npts; j++)
1054     {
1055       if(DebugLevel>5) cout<<"pts[j]="<<pts[j]<<"->lines"<<endl;
1056       if ( Verts[pts[j]].type == VTK_SIMPLE_VERTEX )
1057       {
1058         if ( j == (npts-1) ) //end-of-line marked FIXED
1059         {
1060           if(DebugLevel>5) cout<<"pts[j]="<<pts[j]<<"2:VTK_FIXED_VERTEX 1"<<endl;
1061           Verts[pts[j]].type = VTK_FIXED_VERTEX;
1062         }
1063         else if ( j == 0 ) //beginning-of-line marked FIXED
1064         {
1065           if(DebugLevel>5) cout<<"pts[j]="<<pts[j]<<"3:VTK_FIXED_VERTEX 2"<<endl;
1066           Verts[pts[0]].type = VTK_FIXED_VERTEX;
1067           inPts->GetPoint(pts[0],x2);
1068           inPts->GetPoint(pts[1],x3);
1069         }
1070         else //is edge vertex (unless already edge vertex!)
1071         {
1072           if(DebugLevel>5) cout<<"pts[j]="<<pts[j]<<"4:VTK_FEATURE_EDGE_VERTEX"<<endl;
1073           Verts[pts[j]].type = VTK_FEATURE_EDGE_VERTEX;
1074           Verts[pts[j]].edges = vtkIdList::New();
1075           Verts[pts[j]].edges->SetNumberOfIds(2);
1076           Verts[pts[j]].edges->SetId(0,pts[j-1]);
1077           Verts[pts[j]].edges->SetId(1,pts[j+1]);
1078         }
1079       } //if simple vertex
1080
1081       else if ( Verts[pts[j]].type == VTK_FEATURE_EDGE_VERTEX )
1082       { //multiply connected, becomes fixed!
1083         if(DebugLevel>5) cout<<"pts[j]="<<pts[j]<<"5:VTK_FIXED_VERTEX 3"<<endl;
1084         Verts[pts[j]].type = VTK_FIXED_VERTEX;
1085         Verts[pts[j]].edges->Delete();
1086         Verts[pts[j]].edges = NULL;
1087       }
1088
1089     } //for all points in this line
1090   } //for all lines
1091
1092   if(DebugLevel>5) cout<<"==polygons and triangle strips=="<<endl;
1093     // now polygons and triangle strips-------------------------------
1094   inPolys=input->GetPolys();
1095   numPolys = inPolys->GetNumberOfCells();
1096   inStrips=input->GetStrips();
1097   numStrips = inStrips->GetNumberOfCells();
1098
1099   if(DebugLevel>5) cout<<"numPolys="<<numPolys<<endl;
1100   if(DebugLevel>5) cout<<"numStrips="<<numStrips<<endl;
1101
1102
1103   if ( numPolys > 0 || numStrips > 0 )
1104   { //build cell structure
1105     vtkCellArray *polys;
1106     vtkIdType cellId;
1107     int numNei, nei, edge;
1108     vtkIdType numNeiPts;
1109     vtkIdType *neiPts;
1110     double normal[3], neiNormal[3];
1111     vtkIdList *neighbors;
1112
1113     neighbors = vtkIdList::New();
1114     neighbors->Allocate(VTK_CELL_SIZE);
1115
1116     inMesh = vtkPolyData::New();
1117     inMesh->SetPoints(inPts);
1118     inMesh->SetPolys(inPolys);
1119     Mesh = inMesh;
1120
1121     if ( (numStrips = inStrips->GetNumberOfCells()) > 0 )
1122     { // convert data to triangles
1123       inMesh->SetStrips(inStrips);
1124       toTris = vtkTriangleFilter::New();
1125       toTris->SetInput(inMesh);
1126       toTris->Update();
1127       Mesh = toTris->GetOutput();
1128     }
1129
1130     Mesh->BuildLinks(); //to do neighborhood searching
1131     polys = Mesh->GetPolys();
1132
1133     for (cellId=0, polys->InitTraversal(); polys->GetNextCell(npts,pts);
1134          cellId++)
1135     {
1136       if(DebugLevel>5) cout<<"->cellId="<<cellId<<endl;
1137       for (i=0; i < npts; i++)
1138       {
1139         if(DebugLevel>5) cout<<"-->i="<<i<<endl;
1140         p1 = pts[i];
1141         p2 = pts[(i+1)%npts];
1142
1143         if ( Verts[p1].edges == NULL )
1144         {
1145           Verts[p1].edges = vtkIdList::New();
1146           Verts[p1].edges->Allocate(16,6);
1147         }
1148         if ( Verts[p2].edges == NULL )
1149         {
1150           Verts[p2].edges = vtkIdList::New();
1151           Verts[p2].edges->Allocate(16,6);
1152         }
1153
1154         Mesh->GetCellEdgeNeighbors(cellId,p1,p2,neighbors);
1155         numNei = neighbors->GetNumberOfIds();
1156         if(DebugLevel>5) cout<<"-->numNei="<<numNei<<endl;
1157
1158         edge = VTK_SIMPLE_VERTEX;
1159         if ( numNei == 0 )
1160         {
1161           edge = VTK_BOUNDARY_EDGE_VERTEX;
1162         }
1163
1164         else if ( numNei >= 2 )
1165         {
1166             // check to make sure that this edge hasn't been marked already
1167           for (j=0; j < numNei; j++)
1168           {
1169             if ( neighbors->GetId(j) < cellId )
1170             {
1171               break;
1172             }
1173           }
1174           if ( j >= numNei )
1175           {
1176             edge = VTK_FEATURE_EDGE_VERTEX;
1177           }
1178         }
1179
1180         else if ( numNei == 1 && (nei=neighbors->GetId(0)) > cellId )
1181         {
1182           vtkPolygon::ComputeNormal(inPts,npts,pts,normal);
1183           Mesh->GetCellPoints(nei,numNeiPts,neiPts);
1184           vtkPolygon::ComputeNormal(inPts,numNeiPts,neiPts,neiNormal);
1185
1186           if ( this->FeatureEdgeSmoothing &&
1187                vtkMath::Dot(normal,neiNormal) <= CosFeatureAngle )
1188           {
1189             edge = VTK_FEATURE_EDGE_VERTEX;
1190           }
1191         }
1192         else // a visited edge; skip rest of analysis
1193         {
1194           continue;
1195         }
1196
1197         if ( edge && Verts[p1].type == VTK_SIMPLE_VERTEX )
1198         {
1199           Verts[p1].edges->Reset();
1200           Verts[p1].edges->InsertNextId(p2);
1201           Verts[p1].type = edge;
1202         }
1203         else if ( (edge && Verts[p1].type == VTK_BOUNDARY_EDGE_VERTEX) ||
1204                   (edge && Verts[p1].type == VTK_FEATURE_EDGE_VERTEX) ||
1205                   (!edge && Verts[p1].type == VTK_SIMPLE_VERTEX ) )
1206         {
1207           Verts[p1].edges->InsertNextId(p2);
1208           if ( Verts[p1].type && edge == VTK_BOUNDARY_EDGE_VERTEX )
1209           {
1210             Verts[p1].type = VTK_BOUNDARY_EDGE_VERTEX;
1211           }
1212         }
1213
1214         if ( edge && Verts[p2].type == VTK_SIMPLE_VERTEX )
1215         {
1216           Verts[p2].edges->Reset();
1217           Verts[p2].edges->InsertNextId(p1);
1218           Verts[p2].type = edge;
1219         }
1220         else if ( (edge && Verts[p2].type == VTK_BOUNDARY_EDGE_VERTEX ) ||
1221                   (edge && Verts[p2].type == VTK_FEATURE_EDGE_VERTEX) ||
1222                   (!edge && Verts[p2].type == VTK_SIMPLE_VERTEX ) )
1223         {
1224           Verts[p2].edges->InsertNextId(p1);
1225           if ( Verts[p2].type && edge == VTK_BOUNDARY_EDGE_VERTEX )
1226           {
1227             Verts[p2].type = VTK_BOUNDARY_EDGE_VERTEX;
1228           }
1229         }
1230       }
1231     }
1232
1233     inMesh->Delete();
1234     if (toTris) {toTris->Delete();}
1235
1236     neighbors->Delete();
1237   }//if strips or polys
1238
1239     //post-process edge vertices to make sure we can smooth them
1240   for (i=0; i<numPts; i++)
1241   {
1242     if ( Verts[i].type == VTK_SIMPLE_VERTEX )
1243     {
1244       numSimple++;
1245     }
1246
1247     else if ( Verts[i].type == VTK_FIXED_VERTEX )
1248     {
1249       numFixed++;
1250     }
1251
1252     else if ( Verts[i].type == VTK_FEATURE_EDGE_VERTEX ||
1253               Verts[i].type == VTK_BOUNDARY_EDGE_VERTEX )
1254     { //see how many edges; if two, what the angle is
1255
1256       if ( !this->BoundarySmoothing &&
1257            Verts[i].type == VTK_BOUNDARY_EDGE_VERTEX )
1258       {
1259         if(DebugLevel>5) cout<<"Verts[i].type = VTK_FIXED_VERTEX; 4"<<endl;
1260         Verts[i].type = VTK_FIXED_VERTEX;
1261         numBEdges++;
1262       }
1263
1264       else if ( (npts = Verts[i].edges->GetNumberOfIds()) != 2 )
1265       {
1266         if(DebugLevel>5) cout<<"Verts["<<i<<"].type = VTK_FIXED_VERTEX; 5"<<endl;
1267         Verts[i].type = VTK_FIXED_VERTEX;
1268         numFixed++;
1269       }
1270
1271       else //check angle between edges
1272       {
1273         inPts->GetPoint(Verts[i].edges->GetId(0),x1);
1274         inPts->GetPoint(i,x2);
1275         inPts->GetPoint(Verts[i].edges->GetId(1),x3);
1276
1277         for (k=0; k<3; k++)
1278         {
1279           l1[k] = x2[k] - x1[k];
1280           l2[k] = x3[k] - x2[k];
1281         }
1282         if ( vtkMath::Normalize(l1) >= 0.0 &&
1283              vtkMath::Normalize(l2) >= 0.0 &&
1284              vtkMath::Dot(l1,l2) < CosEdgeAngle)
1285         {
1286           if(DebugLevel>5) cout<<"Verts["<<i<<"].type = VTK_FIXED_VERTEX; 6"<<endl;
1287           Verts[i].type = VTK_FIXED_VERTEX;
1288           numFixed++;
1289         }
1290         else
1291         {
1292           if ( Verts[i].type == VTK_FEATURE_EDGE_VERTEX )
1293           {
1294             numFEdges++;
1295           }
1296           else
1297           {
1298             numBEdges++;
1299           }
1300         }
1301       }//if along edge
1302     }//if edge vertex
1303   }//for all points
1304
1305   if(DebugLevel>5) {
1306     cout<<"Found\n\t" << numSimple << " simple vertices\n\t"
1307       << numFEdges << " feature edge vertices\n\t"
1308       << numBEdges << " boundary edge vertices\n\t"
1309       << numFixed << " fixed vertices\n\t"<<endl;
1310     cout<<"1:numPts="<<numPts<<endl;
1311   }
1312
1313   for (i=0; i<numPts; i++)
1314   {
1315     if(DebugLevel>5) cout<<"Verts["<<i<<"].type="<<VertexType2Str(Verts[i].type)<<endl;
1316     if(Verts[i].edges != NULL && (npts = Verts[i].edges->GetNumberOfIds()) > 0)
1317     {
1318       for (j=0; j<npts; j++)
1319       {
1320         if(DebugLevel>5) cout<<"Verts["<<i<<"].edges->GetId("<<j<<")="<<Verts[i].edges->GetId(j)<<endl;
1321       }//for all connected points
1322     }
1323   }
1324
1325   //Copy node type info from Verts
1326   EG_VTKDCN(vtkCharArray, node_type, grid, "node_type");
1327   if(DebugLevel>5) cout<<"nodes.size()="<<nodes.size()<<endl;
1328   foreach(vtkIdType node,nodes)
1329   {
1330     if(DebugLevel>5) cout<<"Verts["<<node<<"].type="<<VertexType2Str(Verts[node].type)<<endl;
1331     node_type->SetValue(node,Verts[node].type);
1332   }
1333
1334   //free up connectivity storage
1335   for (int i=0; i<numPts; i++)
1336   {
1337     if ( Verts[i].edges != NULL )
1338     {
1339       Verts[i].edges->Delete();
1340       Verts[i].edges = NULL;
1341     }
1342   }
1343   delete [] Verts;
1344
1345   return(0);
1346 }
1347 //End of UpdateNodeType