guismoothsurface.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 "guismoothsurface.h"
  24 #include "swaptriangles.h"
  25 #include "createspecialmapping.h"
  26 #include "vertexdelegate.h"
  27 #include "settingssheet.h"
  28 #include "laplacesmoother.h"
  29
  30 #include <vtkSmoothPolyDataFilter.h>
  31 #include <vtksmoothpolydatafilter2.h>
  32 #include <vtkWindowedSincPolyDataFilter.h>
  33
  34 #include <vtkLongArray.h>
  35 #include <QtGui>
  36
  37 #include <QTextStream>
  38 #include <stdio.h>
  39
  40 #include "vtkEgNormalExtrusion.h"
  41 #include "containertricks.h"
  42
  43 #include <iostream>
  44 #include <fstream>
  45 #include <vtkIdList.h>
  46 #include <vtkCell.h>
  47 #include <cmath>
  48 #include <vtkCellLocator.h>
  49 #include <vtkFloatArray.h>
  50 #include <vtkCellArray.h>
  51
  52 #include "vtkeggridsmoothpolydatafilter.h"
  53 #include "vtkeggridwindowedsincpolydatafilter.h"
  54
  55 ///////////////////////////////////////////
  56 /* Here is how we we get QTextStreams that look like iostreams */
  57 QTextStream Qcin(stdin, QIODevice::ReadOnly);
  58 QTextStream Qcout(stdout, QIODevice::WriteOnly);
  59 QTextStream Qcerr(stderr, QIODevice::WriteOnly);
  60 ///////////////////////////////////////////
  61
  62 //////////////////////////////////////////////
  63
  64 int cout_vtkWindowedSincPolyDataFilter(vtkWindowedSincPolyDataFilter* smooth)
  65 {
  66   cout<<"NumberOfIterations="<<smooth->GetNumberOfIterations()<<endl;
  67   cout<<"PassBand="<<smooth->GetPassBand()<<endl;
  68   cout<<"FeatureEdgeSmoothing="<<smooth->GetFeatureEdgeSmoothing()<<endl;
  69   cout<<"FeatureAngle="<<smooth->GetFeatureAngle()<<endl;
  70   cout<<"EdgeAngle="<<smooth->GetEdgeAngle()<<endl;
  71   cout<<"BoundarySmoothing="<<smooth->GetBoundarySmoothing()<<endl;
  72   cout<<"GenerateErrorScalars="<<smooth->GetGenerateErrorScalars()<<endl;
  73   cout<<"GenerateErrorVectors="<<smooth->GetGenerateErrorVectors()<<endl;
  74   return(0);
  75 }
  76
  77 int cout_vtkSmoothPolyDataFilter(vtkSmoothPolyDataFilter* smooth)
  78 {
  79   cout<<"GetConvergence="<<smooth->GetConvergence()<<endl;
  80   cout<<"GetNumberOfIterations="<<smooth->GetNumberOfIterations()<<endl;
  81   cout<<"GetRelaxationFactor="<<smooth->GetRelaxationFactor()<<endl;
  82   cout<<"GetFeatureEdgeSmoothing="<<smooth->GetFeatureEdgeSmoothing()<<endl;
  83   cout<<"GetFeatureAngle="<<smooth->GetFeatureAngle()<<endl;
  84   cout<<"GetEdgeAngle="<<smooth->GetEdgeAngle()<<endl;
  85   cout<<"GetBoundarySmoothing="<<smooth->GetBoundarySmoothing()<<endl;
  86   cout<<"GetGenerateErrorScalars="<<smooth->GetGenerateErrorScalars()<<endl;
  87   cout<<"GetGenerateErrorVectors="<<smooth->GetGenerateErrorVectors()<<endl;
  88   return(0);
  89 }
  90
  91 ///////////////////////////////////////////
  92
  93 int GuiSmoothSurface::readSettings()
  94 {
  95   local_qset=new QSettings("enGits","enGrid_smoothsurface");
  96   current_filename=local_qset->value("Filename", "").toString();
  97   ui.SmoothMethod->setCurrentIndex(local_qset->value("Method", 0).toInt());
  98   ui.spinBox_NumberOfSmoothIterations->setValue(local_qset->value("NumberOfSmoothIterations", 20).toInt());
  99   ui.spinBox_maxiter_density->setValue(local_qset->value("maxiter_density", 1000).toInt());
 100   ui.spinBox_DebugLevel->setValue(local_qset->value("DebugLevel", 0).toInt());
 101
 102   ui.doubleSpinBox_Convergence_meshdensity->setValue(local_qset->value("Convergence_meshdensity", 0.000001).toDouble());
 103   ui.checkBox_insert_FP->setCheckState(int2CheckState(local_qset->value("insert_FP", 2).toInt()));
 104   ui.checkBox_insert_EP->setCheckState(int2CheckState(local_qset->value("insert_EP", 2).toInt()));
 105   ui.checkBox_remove_FP->setCheckState(int2CheckState(local_qset->value("remove_FP", 2).toInt()));
 106   ui.checkBox_remove_EP->setCheckState(int2CheckState(local_qset->value("remove_EP", 2).toInt()));
 107   ui.checkBox_Swap->setCheckState(int2CheckState(local_qset->value("DoSwap", 2).toInt()));
 108   ui.checkBox_LaplaceSmoothing->setCheckState(int2CheckState(local_qset->value("DoLaplaceSmoothing", 2).toInt()));
 109
 110   if(local_qset->value("DensityUnit_is_length", false).toBool()){
 111     ui.radioButton_length->toggle();
 112   }
 113   else{
 114     ui.radioButton_density->toggle();
 115   }
 116
 117   int size;
 118   size = local_qset->beginReadArray("list_BC");
 119   if(ui.listWidget->count()==size)
 120   {
 121     for (int i = 0; i < size; ++i) {
 122       local_qset->setArrayIndex(i);
 123       Qt::CheckState x=int2CheckState(local_qset->value("state").toInt());
 124       ui.listWidget->item(i)->setCheckState(x);
 125     }
 126     local_qset->endArray();
 127   }
 128
 129   size = local_qset->beginReadArray("list_BC_Source");
 130   if(ui.listWidget_Source->count()==size)
 131   {
 132     for (int i = 0; i < size; ++i) {
 133       local_qset->setArrayIndex(i);
 134       Qt::CheckState x=int2CheckState(local_qset->value("state").toInt());
 135       ui.listWidget_Source->item(i)->setCheckState(x);
 136     }
 137     local_qset->endArray();
 138   }
 139   return(0);
 140 }
 141
 142 int GuiSmoothSurface::writeSettings()
 143 {
 144   local_qset=new QSettings("enGits","enGrid_smoothsurface");
 145   local_qset->setValue("Filename", current_filename);
 146   local_qset->setValue("Method", ui.SmoothMethod->currentIndex());
 147   local_qset->setValue("NumberOfSmoothIterations", ui.spinBox_NumberOfSmoothIterations->value());
 148   local_qset->setValue("maxiter_density", ui.spinBox_maxiter_density->value());
 149   local_qset->setValue("DebugLevel", ui.spinBox_DebugLevel->value());
 150   local_qset->setValue("Convergence_meshdensity", ui.doubleSpinBox_Convergence_meshdensity->value());
 151   local_qset->setValue("insert_FP", ui.checkBox_insert_FP->checkState());
 152   local_qset->setValue("insert_EP", ui.checkBox_insert_EP->checkState());
 153   local_qset->setValue("remove_FP", ui.checkBox_remove_FP->checkState());
 154   local_qset->setValue("remove_EP", ui.checkBox_remove_EP->checkState());
 155   local_qset->setValue("DoSwap", ui.checkBox_Swap->checkState());
 156   local_qset->setValue("DoLaplaceSmoothing", ui.checkBox_LaplaceSmoothing->checkState());
 157   local_qset->setValue("DensityUnit_is_length",ui.radioButton_length->isChecked());
 158
 159   QList<Qt::CheckState> list;
 160
 161   for (int i = 0; i < ui.listWidget->count(); ++i) {
 162     list << ui.listWidget->item(i)->checkState();
 163   };
 164   local_qset->beginWriteArray("list_BC");
 165   for (int i = 0; i < list.size(); ++i) {
 166     local_qset->setArrayIndex(i);
 167     local_qset->setValue("state", list.at(i));
 168   }
 169   local_qset->endArray();
 170
 171   list.clear();
 172   for (int i = 0; i < ui.listWidget_Source->count(); ++i) {
 173     list << ui.listWidget_Source->item(i)->checkState();
 174   };
 175   local_qset->beginWriteArray("list_BC_Source");
 176   for (int i = 0; i < list.size(); ++i) {
 177     local_qset->setArrayIndex(i);
 178     local_qset->setValue("state", list.at(i));
 179   }
 180   local_qset->endArray();
 181
 182   return(0);
 183 }
 184
 185 ///////////////////////////////////////////
 186
 187 void GuiSmoothSurface::before()
 188 {
 189
 190   tableWidget=new SettingsSheet();
 191   ui.verticalLayout_SettingsSheet->addWidget(tableWidget);
 192
 193   populateBoundaryCodes(ui.listWidget);
 194   populateBoundaryCodes(ui.listWidget_Source);
 195
 196   ui.SmoothMethod->addItem("Method 0: vtkSmoothPolyDataFilter smoothing");
 197   ui.SmoothMethod->addItem("Method 1: vtkWindowedSincPolyDataFilter smoothing");
 198   ui.SmoothMethod->addItem("Method 2: edge subdivision");
 199   ui.SmoothMethod->addItem("Method 3: swap triangles");
 200   ui.SmoothMethod->addItem("Method 4: center subdivision");
 201   ui.SmoothMethod->addItem("Method 5: boundary refinement");
 202   ui.SmoothMethod->addItem("Method 6: Laplacian smoothing");
 203   ui.SmoothMethod->addItem("Method 7: VertexAvgDist test");
 204   ui.SmoothMethod->addItem("Method 8: Create mesh density map");
 205   ui.SmoothMethod->addItem("Method 9: vtkWindowedSincPolyDataFilter smoothing");
 206   ui.SmoothMethod->addItem("Method 10: Super smoothing :)");
 207   ui.SmoothMethod->addItem("Method 11: Update current mesh density + node types");
 208   ui.SmoothMethod->addItem("Method 12: Delete all possible points :)");
 209   ui.SmoothMethod->addItem("Method 13");
 210   ui.SmoothMethod->addItem("Method 14");
 211   ui.SmoothMethod->addItem("Method 15");
 212   ui.SmoothMethod->addItem("Method 16");
 213   ui.SmoothMethod->addItem("Method 17");
 214   ui.SmoothMethod->addItem("Method 18");
 215
 216   //Load settings
 217   readSettings();
 218
 219   if(ui.listWidget->count()>0) ui.lineEdit_BoundaryCode-> setText(ui.listWidget->item(0)->text());
 220   else ui.lineEdit_BoundaryCode-> setText("42");
 221   ui.spinBox_NumberOfSubdivisions->setValue(1);
 222
 223   vtkSmoothPolyDataFilter* smooth=vtkSmoothPolyDataFilter::New();
 224   vtkWindowedSincPolyDataFilter* smooth2=vtkWindowedSincPolyDataFilter::New();
 225
 226   cout_vtkSmoothPolyDataFilter(smooth);
 227   cout_vtkWindowedSincPolyDataFilter(smooth2);
 228
 229   ui.doubleSpinBox_Convergence->setValue(smooth->GetConvergence());
 230 //   ui.spinBox_NumberOfIterations->setValue(smooth->GetNumberOfIterations());
 231 //   ui.spinBox_NumberOfIterations->setValue(1000);
 232   ui.spinBox_NumberOfIterations->setValue(1);
 233   QString tmp;
 234   ui.lineEdit_RelaxationFactor->setText(tmp.setNum(smooth->GetRelaxationFactor()));
 235   ui.doubleSpinBox_PassBand->setValue(smooth2->GetPassBand());
 236 //   ui.checkBox_FeatureEdgeSmoothing->setCheckState(int2CheckState(smooth->GetFeatureEdgeSmoothing()));
 237   ui.checkBox_FeatureEdgeSmoothing->setCheckState(Qt::Checked);
 238 //   ui.checkBox_FeatureEdgeSmoothing->setFlags(Qt::ItemIsTristate | ui.checkBox_FeatureEdgeSmoothing->flags);
 239   ui.doubleSpinBox_FeatureAngle->setValue(smooth->GetFeatureAngle());
 240   ui.doubleSpinBox_EdgeAngle->setValue(smooth->GetEdgeAngle());
 241   ui.checkBox_BoundarySmoothing->setCheckState(int2CheckState(smooth->GetBoundarySmoothing()));
 242 /*  ui.checkBox_GenerateErrorScalars->setCheckState(int2CheckState(smooth->GetGenerateErrorScalars()));
 243   ui.checkBox_GenerateErrorVectors->setCheckState(int2CheckState(smooth->GetGenerateErrorVectors()));*/
 244   ui.checkBox_GenerateErrorScalars->setCheckState(int2CheckState(2));
 245   ui.checkBox_GenerateErrorVectors->setCheckState(int2CheckState(2));
 246
 247 /*  ui.doubleSpinBox_VTK_SIMPLE_VERTEX->setValue(-1);
 248   ui.doubleSpinBox_VTK_FIXED_VERTEX->setValue(-1);
 249   ui.doubleSpinBox_VTK_FEATURE_EDGE_VERTEX->setValue(-1);
 250   ui.doubleSpinBox_VTK_BOUNDARY_EDGE_VERTEX->setValue(-1);*/
 251
 252   int row=0;
 253   int column=0;
 254   QTableWidgetItem *newItem = new QTableWidgetItem(tr("%1").arg((row+1)*(column+1)));
 255   tableWidget->setItem(row, column, newItem);
 256
 257   cout<<"tableWidget->rowCount()="<<tableWidget->rowCount()<<endl;
 258   cout<<"tableWidget->columnCount()="<<tableWidget->columnCount()<<endl;
 259   int Nrow,Ncol;
 260   Nrow=tableWidget->rowCount();
 261   Ncol=tableWidget->columnCount();
 262
 263   QList<QString> list;
 264   list << "VTK_SIMPLE_VERTEX"
 265     <<"VTK_FIXED_VERTEX"
 266     <<"VTK_FEATURE_EDGE_VERTEX"
 267     <<"VTK_BOUNDARY_EDGE_VERTEX"
 268     <<"any";
 269
 270   QList<QString> list2;
 271   list2 << "yes"
 272     <<"no"
 273     <<"any";
 274
 275   Nbc=ui.listWidget-> count ();
 276   tableWidget->setColumnCount(Nbc+3);
 277   tableWidget->setItemDelegate(new VertexDelegate(Nbc, list));
 278
 279   QStringList L;
 280   for(int i=0;i<Nbc;i++)
 281   {
 282 //     Qcout<<"BASE!!!="<<ui.listWidget->item(i)->text()<<endl;
 283     L<<ui.listWidget->item(i)->text();
 284   }
 285   L<<"Vertex Type";
 286   L<<"Nodelist";
 287   L<<"Mesh Density";
 288   tableWidget->setHorizontalHeaderLabels(L);
 289   tableWidget->resizeColumnsToContents();
 290
 291   if (!current_filename.isEmpty() && !tableWidget->readFile(current_filename,0)) {
 292     cout<<"Loading failed"<<endl;
 293   }
 294
 295   connect(ui.pushButton_AddSet, SIGNAL(clicked()), this, SLOT(AddSet()));
 296   connect(ui.pushButton_RemoveSet, SIGNAL(clicked()), this, SLOT(RemoveSet()));
 297   connect(ui.pushButton_TestSet, SIGNAL(clicked()), this, SLOT(TestSet()));
 298   connect(ui.pushButton_Load, SIGNAL(clicked()), this, SLOT(Load()));
 299   connect(ui.pushButton_Save, SIGNAL(clicked()), this, SLOT(Save()));
 300   connect(ui.pushButton_SelectAll_BC, SIGNAL(clicked()), this, SLOT(SelectAll_BC()));
 301   connect(ui.pushButton_ClearAll_BC, SIGNAL(clicked()), this, SLOT(ClearAll_BC()));
 302   connect(ui.pushButton_SelectAll_Source, SIGNAL(clicked()), this, SLOT(SelectAll_Source()));
 303   connect(ui.pushButton_ClearAll_Source, SIGNAL(clicked()), this, SLOT(ClearAll_Source()));
 304 };
 305
 306 void GuiSmoothSurface::Load()
 307 {
 308   current_filename = QFileDialog::getOpenFileName(this,tr("Open SettingsSheet"), ".",tr("SettingsSheet files (*.sp)"));
 309   if (!current_filename.isEmpty() && !tableWidget->readFile(current_filename)) {
 310     cout<<"Loading failed"<<endl;
 311   }
 312 }
 313 void GuiSmoothSurface::Save()
 314 {
 315   current_filename = QFileDialog::getSaveFileName(this,tr("Save SettingsSheet as..."), ".",tr("SettingsSheet files (*.sp)"));
 316   if (!current_filename.isEmpty() && !tableWidget->writeFile(current_filename)) {
 317     cout<<"Saving failed"<<endl;
 318   }
 319 }
 320
 321 void GuiSmoothSurface::SelectAll_BC()
 322 {
 323   for (int i = 0; i < ui.listWidget->count(); ++i) {
 324     ui.listWidget->item(i)->setCheckState(Qt::Checked);
 325   };
 326 }
 327 void GuiSmoothSurface::ClearAll_BC()
 328 {
 329   for (int i = 0; i < ui.listWidget->count(); ++i) {
 330     ui.listWidget->item(i)->setCheckState(Qt::Unchecked);
 331   };
 332 }
 333 void GuiSmoothSurface::SelectAll_Source()
 334 {
 335   for (int i = 0; i < ui.listWidget_Source->count(); ++i) {
 336     ui.listWidget_Source->item(i)->setCheckState(Qt::Checked);
 337   };
 338 }
 339 void GuiSmoothSurface::ClearAll_Source()
 340 {
 341   for (int i = 0; i < ui.listWidget_Source->count(); ++i) {
 342     ui.listWidget_Source->item(i)->setCheckState(Qt::Unchecked);
 343   };
 344 }
 345
 346 void GuiSmoothSurface::TestSet()
 347 {
 348   cout<<"Testing set"<<endl;
 349   GetSet();
 350 }
 351
 352 //This is where we get the user defined mesh densities
 353 QVector <VertexMeshDensity> GuiSmoothSurface::GetSet()
 354 {
 355   cout<<"Getting set"<<endl;
 356   QVector <VertexMeshDensity> VMDvector;
 357   cout<<"VMDvector:"<<VMDvector<<endl;
 358
 359   int N_VMD=tableWidget->rowCount();
 360   VMDvector.resize(N_VMD);
 361   cout<<"VMDvector.size()="<<VMDvector.size()<<endl;
 362   for(int i=0;i<N_VMD;i++)
 363   {
 364     for(int j=0;j<Nbc;j++)
 365     {
 366       if(tableWidget->item(i,j)->checkState()) VMDvector[i].BClist.push_back(tableWidget->horizontalHeaderItem(j)->text().toInt());
 367     }
 368     VMDvector[i].type=Str2VertexType(tableWidget->item(i,Nbc)->text());
 369     VMDvector[i].SetNodes(tableWidget->item(i,Nbc+1)->text());
 370     if(ui.radioButton_density->isChecked()){
 371       VMDvector[i].density=tableWidget->item(i,Nbc+2)->text().toDouble();
 372     }
 373     else{
 374       cout<<"ze_density="<<1.0/(tableWidget->item(i,Nbc+2)->text().toDouble())<<endl;
 375       VMDvector[i].density=1.0/(tableWidget->item(i,Nbc+2)->text().toDouble());
 376     }
 377   }
 378   cout<<"VMDvector:"<<VMDvector<<endl;
 379
 380   return(VMDvector);
 381 }
 382
 383 void GuiSmoothSurface::AddSet()
 384 {
 385   cout<<"Adding set"<<endl;
 386   int row=tableWidget->rowCount();
 387   tableWidget->insertRow(row);
 388
 389   int Nbc=ui.listWidget->count();
 390   for(int i=0;i<Nbc;i++)
 391   {
 392     TriStateTableWidgetItem *newBC = new TriStateTableWidgetItem();
 393     newBC->setFlags(Qt::ItemIsTristate | Qt::ItemIsEnabled | Qt::ItemIsUserCheckable);
 394     tableWidget->setItem(row, i, newBC);
 395   }
 396   QTableWidgetItem *item;
 397   item = new QTableWidgetItem("any");
 398   tableWidget->setItem(row, Nbc, item);
 399   item = new QTableWidgetItem("");
 400   tableWidget->setItem(row, Nbc+1, item);
 401   item = new QTableWidgetItem("-1");
 402   tableWidget->setItem(row, Nbc+2, item);
 403   tableWidget->resizeColumnsToContents();
 404 }
 405
 406 void GuiSmoothSurface::RemoveSet()
 407 {
 408   cout<<"Removing set"<<endl;
 409   tableWidget->removeRow(tableWidget->currentRow());
 410   tableWidget->resizeColumnsToContents();
 411 }
 412
 413 int GuiSmoothSurface::DisplayErrorScalars(vtkPolyDataAlgorithm* algo)
 414 {
 415
 416   cout<<"==============="<<endl;
 417   cout<<"ErrorScalars:"<<endl;
 418   int N1,N2;
 419   double x1[3], x2[3], x3[3], l1[3], l2[3];
 420   double dist;
 421   int numPts=0;
 422   N1=algo->GetOutput()->GetPointData()->GetNumberOfArrays();
 423 //   cout<<"nb of arrays="<<N1<<endl;
 424   algo->GetOutput();//vtkPolyData
 425 //   cout<<algo->GetOutput()->GetPointData()<<endl;//vtkPointData*
 426   vtkFloatArray *newScalars = vtkFloatArray::New();
 427   newScalars=(vtkFloatArray *)algo->GetOutput()->GetPointData()->GetArray(1);
 428   N1=newScalars->GetNumberOfComponents();
 429   N2=newScalars->GetNumberOfTuples();
 430   cout<<"Number of components=N1="<<N1<<endl;
 431   cout<<"Number of tuples=N2="<<N2<<endl;
 432   for (int i=0; i<N2; i++)
 433   {
 434     dist=newScalars->GetComponent(i-1,1);//strange, but works. O.o
 435     cout<<"dist["<<i<<"]="<<dist<<endl;
 436   }
 437
 438   cout<<"==============="<<endl;
 439   return(0);
 440 }
 441
 442 int GuiSmoothSurface::DisplayErrorVectors(vtkPolyDataAlgorithm* algo)
 443 {
 444   cout<<"==============="<<endl;
 445   cout<<"ErrorVectors:"<<endl;
 446 //   int N1,N2;
 447 /*  vec3_t xx;
 448   algo->GetOutput()->GetPoint(0, xx.data());*/
 449   int N1=algo->GetOutput()->GetPointData()->GetVectors()->GetNumberOfComponents();
 450   int N2=algo->GetOutput()->GetPointData()->GetVectors()->GetNumberOfTuples();
 451   cout<<"Number of components=N1="<<N1<<endl;
 452   cout<<"Number of tuples=N2="<<N2<<endl;
 453 /*  vtkPointData* newPointData = vtkPointData::New();
 454   newPointData=algo->GetOutput()->GetPointData();
 455   cout<<"Are you aware?:"<<newPointData->IsA("vtkDataArray")<<endl;*/
 456
 457 /*  void vtkFieldData::GetTuple         (       const vtkIdType          i,
 458                                         double *        tuple
 459                                 )
 460     Copy the ith tuple value into a user provided tuple array. Make sure that you've allocated enough space for the copy.
 461       Deprecated:
 462       as of VTK 5.2. Using this method for FieldData having arrays that are not subclasses of vtkDataArray may yield unexpected results. */
 463   //Yes, indeed, very unexpected... And what should we do instead?
 464
 465   for(vtkIdType i=0;i<N2;i++)
 466   {
 467     cout<<"WTF!"<<endl;
 468     double tuple[4];
 469     algo->GetOutput()->GetPointData()->GetTuple(i,tuple);
 470     cout<<"tuple["<<tuple[0]<<"]=("<<tuple[1]<<","<<tuple[2]<<","<<tuple[3]<<")"<<endl;//TODO: This works, but seems incorrect
 471   }
 472   cout<<"==============="<<endl;
 473   return(0);
 474 }
 475
 476 void GuiSmoothSurface::operate()
 477 {
 478   //Save settings
 479   writeSettings();
 480
 481   if(!current_filename.isEmpty()) tableWidget->writeFile(current_filename);
 482
 483   cout<<"METHOD "<<ui.SmoothMethod->currentIndex()<<endl;
 484   //can't use switch case because dynamic variables seem to be forbidden inside case statements
 485   //////////////////////////////////////////////////////////////////////////////////////////////
 486   if(ui.SmoothMethod->currentIndex()==0)//vtkSmoothPolyDataFilter smoothing
 487   {
 488     //preparations
 489     QSet<int> bcs;
 490     getSelectedItems(ui.listWidget, bcs);
 491     QVector<vtkIdType> cells;
 492     getSurfaceCells(bcs, cells, grid);
 493     EG_VTKSP(vtkPolyData, pdata);
 494     addToPolyData(cells, pdata, grid);
 495     EG_VTKSP(vtkSmoothPolyDataFilter, smooth);
 496
 497 //     EG_VTKSP(vtkEgGridSmoothPolyDataFilter, smooth2);
 498 //     vtkSmartPointer<vtkEgGridSmoothPolyDataFilter> smooth2 = vtkSmartPointer<vtkEgGridSmoothPolyDataFilter>::New();
 499
 500     cout_vtkSmoothPolyDataFilter(smooth);
 501
 502     //configure vtkSmoothPolyDataFilter
 503     smooth->SetInput(pdata);
 504
 505     smooth->SetConvergence (ui.doubleSpinBox_Convergence->value());
 506     smooth->SetNumberOfIterations (ui.spinBox_NumberOfIterations->value());
 507     smooth->SetRelaxationFactor (ui.lineEdit_RelaxationFactor->text().toDouble());
 508     smooth->SetFeatureEdgeSmoothing (ui.checkBox_FeatureEdgeSmoothing->checkState());
 509     smooth->SetFeatureAngle (ui.doubleSpinBox_FeatureAngle->value());
 510     smooth->SetEdgeAngle (ui.doubleSpinBox_EdgeAngle->value());
 511     smooth->SetBoundarySmoothing (ui.checkBox_BoundarySmoothing->checkState());
 512     smooth->SetGenerateErrorScalars (ui.checkBox_GenerateErrorScalars->checkState());
 513     smooth->SetGenerateErrorVectors (ui.checkBox_GenerateErrorVectors->checkState());
 514
 515     QSet<int> bcs_Source;
 516     getSelectedItems(ui.listWidget,bcs_Source);
 517     QVector<vtkIdType> cells_Source;
 518     getSurfaceCells(bcs_Source, cells_Source, grid);
 519     EG_VTKSP(vtkPolyData, pdata_Source);
 520     addToPolyData(cells_Source, pdata_Source, grid);
 521     smooth->SetSource (pdata_Source);
 522
 523     cout_vtkSmoothPolyDataFilter(smooth);
 524
 525     //smooth
 526     smooth->Update();
 527
 528     if(ui.checkBox_GenerateErrorScalars->checkState()) DisplayErrorScalars(smooth);
 529     if(ui.checkBox_GenerateErrorVectors->checkState()) DisplayErrorVectors(smooth);
 530
 531     //copy smoothed grid to main grid
 532     EG_VTKDCN(vtkLongArray_t, node_index, pdata, "node_index");
 533     for (vtkIdType i = 0; i < smooth->GetOutput()->GetNumberOfPoints(); ++i) {
 534       vec3_t x;
 535       smooth->GetOutput()->GetPoints()->GetPoint(i, x.data());
 536       vtkIdType nodeId = node_index->GetValue(i);
 537       grid->GetPoints()->SetPoint(nodeId, x.data());
 538     };
 539     updateActors();
 540   }
 541   //////////////////////////////////////////////////////////////////////////////////////////////
 542   else if(ui.SmoothMethod->currentIndex()==1)//vtkWindowedSincPolyDataFilter smoothing
 543   {
 544     QSet<int> bcs;
 545     getSelectedItems(ui.listWidget, bcs);
 546     QVector<vtkIdType> cells;
 547     getSurfaceCells(bcs, cells, grid);
 548     EG_VTKSP(vtkPolyData, pdata);
 549     addToPolyData(cells, pdata, grid);
 550     EG_VTKSP(vtkWindowedSincPolyDataFilter, smooth);
 551
 552     cout_vtkWindowedSincPolyDataFilter(smooth);
 553
 554     smooth->SetInput(pdata);
 555
 556     smooth->SetNumberOfIterations(ui.spinBox_NumberOfIterations->value());
 557     smooth->SetPassBand(ui.doubleSpinBox_PassBand->value());
 558     smooth->SetFeatureEdgeSmoothing(ui.checkBox_FeatureEdgeSmoothing->checkState());
 559     smooth->SetFeatureAngle(ui.doubleSpinBox_FeatureAngle->value());
 560     smooth->SetEdgeAngle(ui.doubleSpinBox_EdgeAngle->value());
 561     smooth->SetBoundarySmoothing(ui.checkBox_BoundarySmoothing->checkState());
 562     smooth->SetGenerateErrorScalars(ui.checkBox_GenerateErrorScalars->checkState());
 563     smooth->SetGenerateErrorVectors(ui.checkBox_GenerateErrorVectors->checkState());
 564
 565     cout_vtkWindowedSincPolyDataFilter(smooth);
 566
 567     smooth->Update();
 568
 569     if(ui.checkBox_GenerateErrorScalars->checkState()) DisplayErrorScalars(smooth);
 570     if(ui.checkBox_GenerateErrorVectors->checkState()) DisplayErrorVectors(smooth);
 571
 572     EG_VTKDCN(vtkLongArray_t, node_index, pdata, "node_index");
 573     for (vtkIdType i = 0; i < smooth->GetOutput()->GetNumberOfPoints(); ++i) {
 574       vec3_t x;
 575       smooth->GetOutput()->GetPoints()->GetPoint(i, x.data());
 576       vtkIdType nodeId = node_index->GetValue(i);
 577       grid->GetPoints()->SetPoint(nodeId, x.data());
 578     };
 579     updateActors();
 580   }
 581   //////////////////////////////////////////////////////////////////////////////////////////////
 582   else if(ui.SmoothMethod->currentIndex()==2)//edge subdivision
 583   {
 584     cout_grid(cout,grid);
 585
 586     QSet<int> bcs;
 587     getSelectedItems(ui.listWidget, bcs);
 588
 589     int N_iter=ui.spinBox_NumberOfSubdivisions->value();
 590     for(int i_iter=0;i_iter<N_iter;i_iter++)
 591     {
 592       cout<<"i_iter="<<i_iter<<endl;
 593
 594       QVector<vtkIdType> SelectedCells;
 595       getSurfaceCells(bcs, SelectedCells, grid);
 596       QVector<vtkIdType> AllCells;
 597       getAllSurfaceCells(AllCells,grid);
 598
 599       createCellToCell(AllCells, c2c, grid);
 600
 601       int N_points=grid->GetNumberOfPoints();
 602       int N_cells=grid->GetNumberOfCells();
 603
 604       QMap< pair<vtkIdType,vtkIdType>, vtkIdType> midpoint_map;
 605 //       QMap<double, int> midpoint_map;
 606       int N_extmidpoints=0;
 607       int N_intmidpoints=0;
 608
 609       int N_newpoints=0;
 610       int N_newcells=0;
 611       vtkIdType nodeId = N_points;
 612       foreach(vtkIdType id_cell, SelectedCells)
 613       {
 614         vtkIdType type_cell = grid->GetCellType(id_cell);
 615         int N_neighbours=c2c[id_cell].size();
 616         for(int i=0;i<N_neighbours;i++)
 617         {
 618           vtkIdType id_neighbour=c2c[id_cell][i];
 619           if(id_neighbour<0)
 620           {
 621             N_extmidpoints++;
 622           }
 623           else
 624           {
 625             midpoint_map[OrderedPair(id_cell,id_neighbour)]=nodeId; nodeId++;
 626           }
 627         }
 628         if (type_cell == VTK_TRIANGLE)
 629         {
 630           N_newcells+=3;
 631         }
 632         if (type_cell == VTK_QUAD)
 633         {
 634           N_newcells+=3;
 635           N_intmidpoints++;
 636         }
 637       }
 638
 639       int N_c2cmidpoints=midpoint_map.size();
 640       N_newpoints=N_c2cmidpoints+N_extmidpoints+N_intmidpoints;
 641
 642       cout<<"N_c2cmidpoints="<<N_c2cmidpoints<<endl;
 643       cout<<"N_extmidpoints="<<N_extmidpoints<<endl;
 644       cout<<"N_intmidpoints="<<N_intmidpoints<<endl;
 645       cout<<"N_newpoints="<<N_newpoints<<endl;
 646       cout<<"N_newcells="<<N_newcells<<endl;
 647       cout<<"N_cells="<<N_cells<<endl;
 648       cout<<"N_points="<<N_points<<endl;
 649
 650       EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 651       allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 652       makeCopyNoAlloc(grid, grid_tmp);
 653
 654       EG_VTKDCC(vtkIntArray, cell_code, grid_tmp, "cell_code");
 655
 656       midpoint_map.clear();//clear midpoint_map
 657       nodeId=N_points;//reset nodeId
 658
 659       foreach(vtkIdType id_cell, SelectedCells)
 660       {
 661         vtkIdType N_pts, *pts;
 662         grid->GetCellPoints(id_cell, N_pts, pts);
 663         vtkIdType intmidpoint;
 664         vtkIdType edgemidpoints[4];
 665         vec3_t M[4];
 666
 667         vtkIdType type_cell = grid->GetCellType(id_cell);
 668         int N_neighbours=c2c[id_cell].size();
 669         vec3_t corner[4];
 670         for(int i=0;i<N_neighbours;i++)
 671         {
 672           grid->GetPoints()->GetPoint(pts[i], corner[i].data());
 673         }
 674
 675         for(int i=0;i<N_neighbours;i++)
 676         {
 677           vtkIdType id_neighbour=c2c[id_cell][i];
 678           if(id_neighbour<0)
 679           {
 680             M[i]=0.5*(corner[i]+corner[(i+1)%N_neighbours]);
 681             addPoint(grid_tmp,nodeId,M[i].data());
 682             edgemidpoints[i]=nodeId;
 683             nodeId++;
 684           }
 685           else
 686           {
 687             if(midpoint_map.contains(OrderedPair(id_cell,id_neighbour)))
 688             {
 689               //pt already exists!
 690 /*              cout<<"pt already exists!: i="<<i<<" midpoint_map[OrderedPair(id_cell,id_neighbour)]="<<midpoint_map[OrderedPair(id_cell,id_neighbour)]<<endl;*/
 691               edgemidpoints[i]=midpoint_map[OrderedPair(id_cell,id_neighbour)];
 692             }
 693             else
 694             {
 695               M[i]=0.5*(corner[i]+corner[(i+1)%N_neighbours]);
 696               addPoint(grid_tmp,nodeId,M[i].data());
 697               midpoint_map[OrderedPair(id_cell,id_neighbour)]=nodeId;
 698               edgemidpoints[i]=nodeId;
 699               nodeId++;
 700             }
 701           }
 702         }
 703         if (type_cell == VTK_TRIANGLE)
 704         {
 705           vtkIdType pts_triangle[4][3];
 706           pts_triangle[0][0]=pts[0];//A;
 707           pts_triangle[0][1]=edgemidpoints[0];
 708           pts_triangle[0][2]=edgemidpoints[2];
 709           pts_triangle[1][0]=pts[1];//B;
 710           pts_triangle[1][1]=edgemidpoints[1];
 711           pts_triangle[1][2]=edgemidpoints[0];
 712           pts_triangle[2][0]=pts[2];//C;
 713           pts_triangle[2][1]=edgemidpoints[2];
 714           pts_triangle[2][2]=edgemidpoints[1];
 715           pts_triangle[3][0]=edgemidpoints[0];
 716           pts_triangle[3][1]=edgemidpoints[1];
 717           pts_triangle[3][2]=edgemidpoints[2];
 718
 719           vtkIdType newCellId;
 720           newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[0]);
 721           cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 722           newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[1]);
 723           cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 724           newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[2]);
 725           cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 726           grid_tmp->ReplaceCell(id_cell , 3, pts_triangle[3]);
 727           cell_code->SetValue(id_cell, ui.lineEdit_BoundaryCode->text().toInt());
 728         }
 729         if (type_cell == VTK_QUAD)
 730         {
 731           vec3_t C=0.25*(corner[0]+corner[1]+corner[2]+corner[3]);
 732           addPoint(grid_tmp,nodeId,C.data());
 733           intmidpoint=nodeId;
 734           nodeId++;
 735
 736           vtkIdType pts_quad[4][4];
 737           pts_quad[0][0]=pts[0];
 738           pts_quad[0][1]=edgemidpoints[0];
 739           pts_quad[0][2]=intmidpoint;
 740           pts_quad[0][3]=edgemidpoints[3];
 741           pts_quad[1][0]=pts[1];
 742           pts_quad[1][1]=edgemidpoints[1];
 743           pts_quad[1][2]=intmidpoint;
 744           pts_quad[1][3]=edgemidpoints[0];
 745           pts_quad[2][0]=pts[2];
 746           pts_quad[2][1]=edgemidpoints[2];
 747           pts_quad[2][2]=intmidpoint;
 748           pts_quad[2][3]=edgemidpoints[1];
 749           pts_quad[3][0]=pts[3];
 750           pts_quad[3][1]=edgemidpoints[3];
 751           pts_quad[3][2]=intmidpoint;
 752           pts_quad[3][3]=edgemidpoints[2];
 753
 754           vtkIdType newCellId;
 755           newCellId = grid_tmp->InsertNextCell(VTK_QUAD,4,pts_quad[0]);
 756           cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 757           newCellId = grid_tmp->InsertNextCell(VTK_QUAD,4,pts_quad[1]);
 758           cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 759           newCellId = grid_tmp->InsertNextCell(VTK_QUAD,4,pts_quad[2]);
 760           cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 761           grid_tmp->ReplaceCell(id_cell , 4, pts_quad[3]);
 762           cell_code->SetValue(id_cell, ui.lineEdit_BoundaryCode->text().toInt());
 763         }
 764       }
 765
 766 //       cout_grid(cout,grid_tmp,true,true,true,true);
 767       cout<<"Copying..."<<endl;
 768       makeCopy(grid_tmp,grid);
 769       cout<<"Copy successful"<<endl;
 770     }
 771     cout_grid(cout,grid);
 772     updateActors();
 773   }
 774   //////////////////////////////////////////////////////////////////////////////////////////////
 775   else if(ui.SmoothMethod->currentIndex()==3)//swap triangles
 776   {
 777     cout_grid(cout,grid);
 778
 779     QSet<int> bcs;
 780     getSelectedItems(ui.listWidget, bcs);
 781
 782     QSet<int> bcs_complement=complementary_bcs(bcs,grid,cells);
 783
 784     cout<<"bcs="<<bcs<<endl;
 785     cout<<"bcs_complement="<<bcs_complement<<endl;
 786
 787     SwapTriangles swap;
 788     swap.setGrid(grid);
 789     swap.setBoundaryCodes(bcs_complement);
 790     swap();
 791
 792     cout_grid(cout,grid);
 793     updateActors();
 794   }
 795   //////////////////////////////////////////////////////////////////////////////////////////////
 796   else if(ui.SmoothMethod->currentIndex()==4)//center subdivision
 797   {
 798     cout_grid(cout,grid);
 799
 800     QSet<int> bcs;
 801     getSelectedItems(ui.listWidget, bcs);
 802     QSet<int> bcs_complement=complementary_bcs(bcs,grid,cells);
 803     cout<<"bcs="<<bcs<<endl;
 804     cout<<"bcs_complement="<<bcs_complement<<endl;
 805
 806     int N_iter=ui.spinBox_NumberOfSubdivisions->value();
 807     for(int i_iter=0;i_iter<N_iter;i_iter++)
 808     {
 809       int N_points=grid->GetNumberOfPoints();
 810       int N_cells=grid->GetNumberOfCells();
 811
 812       QVector<vtkIdType> cells;
 813       getSurfaceCells(bcs, cells, grid);
 814
 815       int N_newcells=0;
 816       foreach(int id_cell, cells)
 817       {
 818         vtkIdType type_cell = grid->GetCellType(id_cell);
 819         if (type_cell == VTK_TRIANGLE) N_newcells+=2;
 820         if (type_cell == VTK_QUAD) N_newcells+=3;
 821       }
 822
 823       cout<<"N_newcells="<<N_newcells<<endl;
 824
 825       EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 826       allocateGrid(grid_tmp,N_cells+N_newcells,N_points+1*N_cells);
 827       makeCopyNoAlloc(grid, grid_tmp);
 828       EG_VTKDCC(vtkIntArray, cell_code, grid_tmp, "cell_code");
 829
 830       vtkIdType nodeId=N_points;
 831       foreach(int id_cell, cells)
 832       {
 833         vtkIdType N_pts, *pts;
 834         grid->GetCellPoints(id_cell, N_pts, pts);
 835         vec3_t C(0,0,0);
 836
 837         vtkIdType type_cell = grid->GetCellType(id_cell);
 838         int N_neighbours=N_pts;
 839         cout<<"N_neighbours="<<N_neighbours<<endl;
 840         vec3_t corner[4];
 841         vtkIdType pts_triangle[4][3];
 842         for(int i=0;i<N_neighbours;i++)
 843         {
 844           grid->GetPoints()->GetPoint(pts[i], corner[i].data());
 845           C+=(1/(double)N_neighbours)*corner[i];
 846         }
 847         addPoint(grid_tmp,nodeId,C.data());
 848         vtkIdType intmidpoint=nodeId;
 849         nodeId++;
 850
 851         for(int i=0;i<N_neighbours;i++)
 852         {
 853           pts_triangle[i][0]=pts[i];
 854           pts_triangle[i][1]=pts[(i+1)%N_neighbours];
 855           pts_triangle[i][2]=intmidpoint;
 856           if(i==0)
 857           {
 858             grid_tmp->ReplaceCell(id_cell , 3, pts_triangle[0]);
 859             cell_code->SetValue(id_cell, ui.lineEdit_BoundaryCode->text().toInt());
 860           }
 861           else
 862           {
 863             vtkIdType newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[i]);
 864             cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 865           }
 866         }
 867       }
 868
 869       makeCopy(grid_tmp,grid);
 870     }
 871     cout_grid(cout,grid,true,true,true,true);
 872     updateActors();
 873   }
 874   //////////////////////////////////////////////////////////////////////////////////////////////
 875   else if(ui.SmoothMethod->currentIndex()==5)//boundary refinement
 876   {
 877     cout_grid(cout,grid);
 878
 879     QSet<int> bcs;
 880     getSelectedItems(ui.listWidget, bcs);
 881     QSet<int> bcs_complement=complementary_bcs(bcs,grid,cells);
 882     cout<<"bcs="<<bcs<<endl;
 883     cout<<"bcs_complement="<<bcs_complement<<endl;
 884
 885     int N_iter=ui.spinBox_NumberOfIterations->value();
 886     for(int i_iter=0;i_iter<N_iter;i_iter++){
 887       int N_points=grid->GetNumberOfPoints();
 888       int N_cells=grid->GetNumberOfCells();
 889
 890       QVector<vtkIdType> SelectedCells;
 891       getSurfaceCells(bcs, SelectedCells, grid);
 892       QVector<vtkIdType> AllCells;
 893       getAllSurfaceCells(AllCells,grid);
 894
 895       EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 896
 897 //       QSet <int> cells_to_split;
 898       QVector <stencil_t> StencilVector;
 899
 900       QMap <vtkIdType,bool> marked;
 901
 902       createCellMapping(AllCells, _cells, grid);
 903       createCellToCell(AllCells, c2c, grid);
 904       setCells(AllCells);
 905
 906       cout<<"AllCells.size()="<<AllCells.size()<<endl;
 907       cout<<"SelectedCells.size()="<<SelectedCells.size()<<endl;
 908
 909       int N_newcells=0;
 910       int N_newpoints=0;
 911       foreach(vtkIdType id_cell, SelectedCells)
 912       {
 913         cout<<"==>id_cell="<<id_cell<<endl;
 914         int bc0=cell_code->GetValue(id_cell);
 915         if(!marked[id_cell])
 916         {
 917           vtkIdType N_pts, *pts;
 918           grid->GetCellPoints(id_cell, N_pts, pts);
 919           int count=0;
 920           for(int i=0;i<N_pts;i++)
 921           {
 922             int bc1=cell_code->GetValue(c2c[id_cell][i]);
 923             if(bc0!=bc1) count++;
 924           }
 925           if(count>0)//cell is near at least one neighbour with different cell code
 926           {
 927             int SideToSplit = getLongestSide(id_cell,grid);
 928             cout<<"SideToSplit="<<SideToSplit<<endl;
 929             cout<<"c2c[id_cell][SideToSplit]="<<c2c[id_cell][SideToSplit]<<endl;
 930             for(int i=0;i<3;i++) cout<<"c2c[id_cell]["<<i<<"]="<<c2c[id_cell][i]<<endl;
 931             stencil_t S=getStencil(id_cell,SideToSplit);
 932             if(S.valid){//there is a neighbour cell
 933               if(!marked[S.id_cell2])
 934               {
 935 /*                cells_to_split.insert(S.id_cell1);
 936                 cells_to_split.insert(S.id_cell2);*/
 937                 cout<<"marked["<<S.id_cell1<<"]=true;"<<endl;
 938                 cout<<"marked["<<S.id_cell2<<"]=true;"<<endl;
 939                 marked[S.id_cell1]=true;
 940                 marked[S.id_cell2]=true;
 941                 StencilVector.push_back(S);
 942                 N_newpoints++;
 943                 N_newcells+=2;
 944               }
 945             }
 946             else{//there is no neighbour cell
 947 //               cells_to_split.insert(S.id_cell1);
 948               cout<<"marked["<<S.id_cell1<<"]=true;"<<endl;
 949               marked[S.id_cell1]=true;
 950               StencilVector.push_back(S);
 951               N_newpoints++;
 952               N_newcells+=1;
 953             }
 954           }
 955         }
 956       }
 957
 958 /*      cout<<"cells_to_split.size()="<<cells_to_split.size()<<endl;
 959       cout<<cells_to_split<<endl;*/
 960       cout<<"StencilVector.size()="<<StencilVector.size()<<endl;
 961       cout<<"N_newpoints="<<N_newpoints<<endl;
 962       cout<<"N_newcells="<<N_newcells<<endl;
 963
 964       EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 965       allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 966       makeCopyNoAlloc(grid, grid_tmp);
 967       EG_VTKDCC(vtkIntArray, cell_code_tmp, grid_tmp, "cell_code");
 968
 969       vtkIdType nodeId=N_points;
 970       foreach(stencil_t S,StencilVector)
 971       {
 972         cout<<S<<endl;
 973         vtkIdType N_pts, *pts;
 974         vec3_t A,B;
 975         grid_tmp->GetPoint(S.p[1],A.data());
 976         grid_tmp->GetPoint(S.p[3],B.data());
 977         vec3_t M=0.5*(A+B);
 978         addPoint(grid_tmp,nodeId,M.data());
 979
 980         vtkIdType pts_triangle[4][3];
 981
 982         if(S.valid){//there is a neighbour cell
 983           for(int i=0;i<4;i++)
 984           {
 985             pts_triangle[i][0]=S.p[i];
 986             pts_triangle[i][1]=S.p[(i+1)%4];
 987             pts_triangle[i][2]=nodeId;
 988           }
 989
 990           int bc1=cell_code_tmp->GetValue(S.id_cell1);
 991           int bc2=cell_code_tmp->GetValue(S.id_cell2);
 992
 993           grid_tmp->ReplaceCell(S.id_cell1 , 3, pts_triangle[0]);
 994           cell_code_tmp->SetValue(S.id_cell1, bc1);
 995
 996           grid_tmp->ReplaceCell(S.id_cell2 , 3, pts_triangle[1]);
 997           cell_code_tmp->SetValue(S.id_cell2, bc2);
 998
 999           vtkIdType newCellId;
1000           newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[2]);
1001           cell_code_tmp->SetValue(newCellId, bc2);
1002           newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[3]);
1003           cell_code_tmp->SetValue(newCellId, bc1);
1004         }
1005         else{//there is no neighbour cell
1006           pts_triangle[0][0]=S.p[0];
1007           pts_triangle[0][1]=S.p[1];
1008           pts_triangle[0][2]=nodeId;
1009           pts_triangle[3][0]=S.p[3];
1010           pts_triangle[3][1]=S.p[0];
1011           pts_triangle[3][2]=nodeId;
1012
1013           int bc1=cell_code_tmp->GetValue(S.id_cell1);
1014
1015           grid_tmp->ReplaceCell(S.id_cell1 , 3, pts_triangle[0]);
1016           cell_code_tmp->SetValue(S.id_cell1, bc1);
1017
1018           vtkIdType newCellId;
1019           newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[3]);
1020           cell_code_tmp->SetValue(newCellId, bc1);
1021         }
1022
1023         nodeId++;
1024       }
1025
1026       makeCopy(grid_tmp,grid);
1027     }//end of i_iter loop
1028 //     cout_grid(cout,grid,true,true,true,true);
1029     updateActors();
1030   }
1031   //////////////////////////////////////////////////////////////////////////////////////////////
1032   else if(ui.SmoothMethod->currentIndex()==6)//Laplacian smoothing
1033   {
1034     QSet<int> bcs;
1035     getSelectedItems(ui.listWidget, bcs);
1036
1037     LaplaceSmoother toto;
1038     toto.SetInput(bcs,grid);
1039     toto.SetNumberOfIterations(ui.spinBox_NumberOfIterations->value());
1040
1041     toto();
1042
1043     updateActors();
1044   }
1045   //////////////////////////////////////////////////////////////////////////////////////////////
1046   else if(ui.SmoothMethod->currentIndex()==7)//VertexAvgDist test
1047   {
1048     cout_grid(cout,grid);
1049
1050     QSet<int> bcs;
1051     getSelectedItems(ui.listWidget, bcs);
1052     QSet<int> bcs_complement=complementary_bcs(bcs,grid,cells);
1053     cout<<"bcs="<<bcs<<endl;
1054     cout<<"bcs_complement="<<bcs_complement<<endl;
1055
1056     int N_points=grid->GetNumberOfPoints();
1057     int N_cells=grid->GetNumberOfCells();
1058
1059     QVector<vtkIdType> SelectedCells;
1060     getSurfaceCells(bcs, SelectedCells, grid);
1061     QVector<vtkIdType> AllCells;
1062     getAllSurfaceCells(AllCells,grid);
1063
1064     QSet <vtkIdType> SelectedNodes;
1065     getSurfaceNodes(bcs,SelectedNodes,grid);
1066     createNodeToNode(cells, nodes, _nodes, n2n, grid);
1067
1068     foreach(vtkIdType node,SelectedNodes)
1069     {
1070       cout<<"node="<<node<<" VertexAvgDist="<<CurrentVertexAvgDist(node,n2n,grid)<<endl;
1071     }
1072
1073     updateActors();
1074   }
1075   //////////////////////////////////////////////////////////////////////////////////////////////
1076   else if(ui.SmoothMethod->currentIndex()==8)//Create mesh density map
1077   {
1078     cout_grid(cout,grid);
1079
1080     QSet<int> bcs;
1081     getSelectedItems(ui.listWidget, bcs);
1082     QSet<int> bcs_complement=complementary_bcs(bcs,grid,cells);
1083     cout<<"bcs="<<bcs<<endl;
1084     cout<<"bcs_complement="<<bcs_complement<<endl;
1085
1086     int N_points=grid->GetNumberOfPoints();
1087     int N_cells=grid->GetNumberOfCells();
1088
1089     QVector<vtkIdType> SelectedCells;
1090     getSurfaceCells(bcs, SelectedCells, grid);
1091     QVector<vtkIdType> AllCells;
1092     getAllSurfaceCells(AllCells,grid);
1093
1094     EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
1095     EG_VTKDCN(vtkDoubleArray, node_meshdensity, grid, "node_meshdensity");
1096
1097     QSet <vtkIdType> SelectedNodes;
1098     getSurfaceNodes(bcs,SelectedNodes,grid);
1099     createNodeToNode(cells, nodes, _nodes, n2n, grid);
1100
1101     foreach(vtkIdType node,SelectedNodes)
1102     {
1103       double L=CurrentVertexAvgDist(node,n2n,grid);
1104       double D=1./L;
1105       cout<<"node="<<node<<" VertexAvgDist="<<L<<" Net density="<<D<<endl;
1106       node_meshdensity->SetValue(node, D);
1107     }
1108
1109     int N_iter=ui.spinBox_NumberOfIterations->value();
1110     for(int i_iter=0;i_iter<N_iter;i_iter++)
1111     {
1112       foreach(vtkIdType node,SelectedNodes)
1113       {
1114         double D=DesiredMeshDensity(node,n2n,grid);
1115         double L=1./D;
1116         cout<<"node="<<node<<" VertexAvgDist="<<L<<" Net density="<<D<<endl;
1117         node_meshdensity->SetValue(node, D);
1118       }
1119     }
1120
1121     updateActors();
1122   }
1123   //////////////////////////////////////////////////////////////////////////////////////////////
1124   else if(ui.SmoothMethod->currentIndex()==9)//vtkWindowedSincPolyDataFilter smoothing
1125   {
1126     QSet<int> bcs;
1127     getSelectedItems(ui.listWidget, bcs);
1128     QVector<vtkIdType> cells;
1129     getSurfaceCells(bcs, cells, grid);
1130     EG_VTKSP(vtkPolyData, input);
1131     addToPolyData(cells, input, grid);
1132
1133     EG_VTKSP(vtkSmoothPolyDataFilter, smooth);
1134     smooth->SetInput(input);
1135
1136
1137     int j, k;
1138     vtkIdType npts = 0;
1139     vtkIdType *pts = 0;
1140
1141     vtkCellArray *inVerts, *inLines, *inPolys, *inStrips;
1142
1143     cout<<"input->GetVerts()="<<input->GetVerts()<<endl;
1144     cout<<"input->GetVerts()->GetSize()="<<input->GetVerts()->GetSize()<<endl;
1145     inVerts=input->GetVerts();
1146     inVerts->InitTraversal();
1147     cout<<"inVerts->GetSize()="<<inVerts->GetSize()<<endl;
1148     cout<<"inVerts->GetNextCell(npts,pts)="<<inVerts->GetNextCell(npts,pts)<<endl;
1149     cout<<"inVerts->GetNumberOfCells()="<<inVerts->GetNumberOfCells()<<endl;
1150
1151     cout<<"input->GetVerts()->GetSize()="<<input->GetVerts()->GetSize()<<endl;
1152     cout<<"input->GetVerts()->GetNextCell(npts,pts)="<<input->GetVerts()->GetNextCell(npts,pts)<<endl;
1153     cout<<"input->GetVerts()->GetNumberOfCells()="<<input->GetVerts()->GetNumberOfCells()<<endl;
1154
1155   // check vertices first. Vertices are never smoothed_--------------
1156   for (inVerts=input->GetVerts(), inVerts->InitTraversal();
1157   inVerts->GetNextCell(npts,pts); )
1158     {
1159         cout<<"npts="<<npts<<endl;
1160     for (j=0; j<npts; j++)
1161       {
1162         cout<<"pts["<<j<<"]="<<pts[j]<<endl;
1163       }
1164     }
1165
1166     updateActors();
1167   }
1168   //////////////////////////////////////////////////////////////////////////////////////////////
1169   else if(ui.SmoothMethod->currentIndex()==10)// super smoothing
1170   {
1171     QSet<int> bcs;
1172     getSelectedItems(ui.listWidget, bcs);
1173
1174     QVector <VertexMeshDensity> VMDvector=GetSet();
1175
1176     CreateSpecialMapping toto;
1177
1178     toto.SetInput(bcs,grid);
1179     toto.SetVertexMeshDensityVector(VMDvector);
1180     toto.SetConvergence (ui.doubleSpinBox_Convergence->value());
1181     toto.SetNumberOfIterations (ui.spinBox_NumberOfIterations->value());
1182     toto.SetRelaxationFactor (ui.lineEdit_RelaxationFactor->text().toDouble());
1183     toto.SetFeatureEdgeSmoothing (ui.checkBox_FeatureEdgeSmoothing->checkState());
1184     toto.SetFeatureAngle (ui.doubleSpinBox_FeatureAngle->value());
1185     toto.SetEdgeAngle (ui.doubleSpinBox_EdgeAngle->value());
1186     toto.SetBoundarySmoothing (ui.checkBox_BoundarySmoothing->checkState());
1187     toto.SetGenerateErrorScalars (ui.checkBox_GenerateErrorScalars->checkState());
1188     toto.SetGenerateErrorVectors (ui.checkBox_GenerateErrorVectors->checkState());
1189
1190 /*    toto.Set_SV_value(ui.doubleSpinBox_VTK_SIMPLE_VERTEX->value());
1191     toto.Set_FV_value(ui.doubleSpinBox_VTK_FIXED_VERTEX->value());
1192     toto.Set_FEV_value(ui.doubleSpinBox_VTK_FEATURE_EDGE_VERTEX->value());
1193     toto.Set_BEV_value(ui.doubleSpinBox_VTK_BOUNDARY_EDGE_VERTEX->value());*/
1194
1195     toto.SetConvergence_meshdensity(ui.doubleSpinBox_Convergence_meshdensity->value());
1196
1197     toto.Set_insert_FP(ui.checkBox_insert_FP->checkState());
1198     toto.Set_insert_EP(ui.checkBox_insert_EP->checkState());
1199     toto.Set_remove_FP(ui.checkBox_remove_FP->checkState());
1200     toto.Set_remove_EP(ui.checkBox_remove_EP->checkState());
1201     toto.DoSwap=ui.checkBox_Swap->checkState();
1202     toto.DoLaplaceSmoothing=ui.checkBox_LaplaceSmoothing->checkState();
1203
1204     toto.N_SmoothIterations=ui.spinBox_NumberOfSmoothIterations->value();
1205     toto.maxiter_density=ui.spinBox_maxiter_density->value();
1206     DebugLevel=ui.spinBox_DebugLevel->value();
1207
1208     toto.Process();
1209
1210     updateActors();
1211   }
1212   //////////////////////////////////////////////////////////////////////////////////////////////
1213   else if(ui.SmoothMethod->currentIndex()==11)// Update current mesh density + node types
1214   {
1215     QSet<int> bcs;
1216     getSelectedItems(ui.listWidget, bcs);
1217     CreateSpecialMapping toto;
1218     toto.SetInput(bcs,grid);
1219
1220     SetConvergence(ui.doubleSpinBox_Convergence->value());
1221     SetFeatureEdgeSmoothing(ui.checkBox_FeatureEdgeSmoothing->checkState());
1222     SetFeatureAngle(ui.doubleSpinBox_FeatureAngle->value());
1223     SetEdgeAngle(ui.doubleSpinBox_EdgeAngle->value());
1224     SetBoundarySmoothing(ui.checkBox_BoundarySmoothing->checkState());
1225
1226     UpdateMeshDensity();
1227     UpdateNodeType();
1228     updateActors();
1229   }
1230   //////////////////////////////////////////////////////////////////////////////////////////////
1231   else if(ui.SmoothMethod->currentIndex()==12)// Delete all possible points
1232   {
1233     QSet<int> bcs;
1234     getSelectedItems(ui.listWidget, bcs);
1235     CreateSpecialMapping toto;
1236
1237     SetConvergence(ui.doubleSpinBox_Convergence->value());
1238     SetFeatureEdgeSmoothing(ui.checkBox_FeatureEdgeSmoothing->checkState());
1239     SetFeatureAngle(ui.doubleSpinBox_FeatureAngle->value());
1240     SetEdgeAngle(ui.doubleSpinBox_EdgeAngle->value());
1241     SetBoundarySmoothing(ui.checkBox_BoundarySmoothing->checkState());
1242
1243     int N_newpoints;
1244     int N_newcells;
1245
1246     vtkIdType N_points=grid->GetNumberOfPoints();
1247     vtkIdType N_cells=grid->GetNumberOfCells();
1248
1249     bool Global_DelResult=true;
1250     while(Global_DelResult)
1251     {
1252       Global_DelResult=false;
1253       vtkIdType DeadNode=0;
1254       while(DeadNode<grid->GetNumberOfPoints())
1255       {
1256         bool Local_DelResult=true;
1257         while(Local_DelResult)
1258         {
1259           Local_DelResult=DeletePoint_2(grid,DeadNode,N_newpoints,N_newcells);
1260           if(Local_DelResult) Global_DelResult=true;
1261         }
1262         DeadNode++;
1263       }
1264     }
1265
1266   }
1267   //////////////////////////////////////////////////////////////////////////////////////////////
1268   else
1269   {
1270     cout<<"UNKNOWN METHOD"<<endl;
1271     QSet<int> bcs;
1272     getSelectedItems(ui.listWidget, bcs);
1273     QVector<vtkIdType> cells;
1274     getSurfaceCells(bcs, cells, grid);
1275     setCells(cells);
1276     cout<<"cells="<<cells<<endl;
1277     cout<<"_cells="<<_cells<<endl;
1278     cout<<"nodes="<<nodes<<endl;
1279     cout<<"_nodes="<<_nodes<<endl;
1280 /*    QVector< QSet< int > >    n2c;
1281     QVector< QSet< int > >      n2n;
1282     QVector< QVector< int > >   c2c;*/
1283     cout<<"n2c="<<n2c<<endl;
1284     cout<<"n2n="<<n2n<<endl;
1285     cout<<"c2c="<<c2c<<endl;
1286   }
1287 };