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