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