guismoothsurface.cpp

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