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");
 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     }
 380     VMDvector[i].type=Str2VertexType(tableWidget->item(i,Nbc)->text());
 381     VMDvector[i].SetNodes(tableWidget->item(i,Nbc+1)->text());
 382     if(ui.radioButton_density->isChecked()){
 383       VMDvector[i].density=tableWidget->item(i,Nbc+2)->text().toDouble();
 384     }
 385     else{
 386       cout<<"ze_density="<<1.0/(tableWidget->item(i,Nbc+2)->text().toDouble())<<endl;
 387       VMDvector[i].density=1.0/(tableWidget->item(i,Nbc+2)->text().toDouble());
 388     }
 389   }
 390   cout<<"VMDvector:"<<VMDvector<<endl;
 391
 392   return(VMDvector);
 393 }
 394
 395 void GuiSmoothSurface::AddSet()
 396 {
 397   cout<<"Adding set"<<endl;
 398   int row=tableWidget->rowCount();
 399   tableWidget->insertRow(row);
 400
 401   int Nbc=ui.listWidget->count();
 402   for(int i=0;i<Nbc;i++)
 403   {
 404     TriStateTableWidgetItem *newBC = new TriStateTableWidgetItem();
 405     newBC->setFlags(Qt::ItemIsTristate | Qt::ItemIsEnabled | Qt::ItemIsUserCheckable);
 406     tableWidget->setItem(row, i, newBC);
 407   }
 408   QTableWidgetItem *item;
 409   item = new QTableWidgetItem("any");
 410   tableWidget->setItem(row, Nbc, item);
 411   item = new QTableWidgetItem("");
 412   tableWidget->setItem(row, Nbc+1, item);
 413   item = new QTableWidgetItem("-1");
 414   tableWidget->setItem(row, Nbc+2, item);
 415   tableWidget->resizeColumnsToContents();
 416 }
 417
 418 void GuiSmoothSurface::RemoveSet()
 419 {
 420   cout<<"Removing set"<<endl;
 421   tableWidget->removeRow(tableWidget->currentRow());
 422   tableWidget->resizeColumnsToContents();
 423 }
 424
 425 int GuiSmoothSurface::DisplayErrorScalars(vtkPolyDataAlgorithm* algo)
 426 {
 427
 428   cout<<"==============="<<endl;
 429   cout<<"ErrorScalars:"<<endl;
 430   int N1,N2;
 431   double x1[3], x2[3], x3[3], l1[3], l2[3];
 432   double dist;
 433   int numPts=0;
 434   N1=algo->GetOutput()->GetPointData()->GetNumberOfArrays();
 435 //   cout<<"nb of arrays="<<N1<<endl;
 436   algo->GetOutput();//vtkPolyData
 437 //   cout<<algo->GetOutput()->GetPointData()<<endl;//vtkPointData*
 438   vtkFloatArray *newScalars = vtkFloatArray::New();
 439   newScalars=(vtkFloatArray *)algo->GetOutput()->GetPointData()->GetArray(1);
 440   N1=newScalars->GetNumberOfComponents();
 441   N2=newScalars->GetNumberOfTuples();
 442   cout<<"Number of components=N1="<<N1<<endl;
 443   cout<<"Number of tuples=N2="<<N2<<endl;
 444   for (int i=0; i<N2; i++)
 445   {
 446     dist=newScalars->GetComponent(i-1,1);//strange, but works. O.o
 447     cout<<"dist["<<i<<"]="<<dist<<endl;
 448   }
 449
 450   cout<<"==============="<<endl;
 451   return(0);
 452 }
 453
 454 int GuiSmoothSurface::DisplayErrorVectors(vtkPolyDataAlgorithm* algo)
 455 {
 456   cout<<"==============="<<endl;
 457   cout<<"ErrorVectors:"<<endl;
 458 //   int N1,N2;
 459 /*  vec3_t xx;
 460   algo->GetOutput()->GetPoint(0, xx.data());*/
 461   int N1=algo->GetOutput()->GetPointData()->GetVectors()->GetNumberOfComponents();
 462   int N2=algo->GetOutput()->GetPointData()->GetVectors()->GetNumberOfTuples();
 463   cout<<"Number of components=N1="<<N1<<endl;
 464   cout<<"Number of tuples=N2="<<N2<<endl;
 465 /*  vtkPointData* newPointData = vtkPointData::New();
 466   newPointData=algo->GetOutput()->GetPointData();
 467   cout<<"Are you aware?:"<<newPointData->IsA("vtkDataArray")<<endl;*/
 468
 469 /*  void vtkFieldData::GetTuple         (       const vtkIdType          i,
 470                                         double *        tuple
 471                                 )
 472     Copy the ith tuple value into a user provided tuple array. Make sure that you've allocated enough space for the copy.
 473       Deprecated:
 474       as of VTK 5.2. Using this method for FieldData having arrays that are not subclasses of vtkDataArray may yield unexpected results. */
 475   //Yes, indeed, very unexpected... And what should we do instead?
 476
 477   for(vtkIdType i=0;i<N2;i++)
 478   {
 479     cout<<"WTF!"<<endl;
 480     double tuple[4];
 481     algo->GetOutput()->GetPointData()->GetTuple(i,tuple);
 482     cout<<"tuple["<<tuple[0]<<"]=("<<tuple[1]<<","<<tuple[2]<<","<<tuple[3]<<")"<<endl;//TODO: This works, but seems incorrect
 483   }
 484   cout<<"==============="<<endl;
 485   return(0);
 486 }
 487
 488 void GuiSmoothSurface::operate()
 489 {
 490   cout<<"Saving settings..."<<endl;
 491   //Save settings
 492   writeSettings();
 493
 494   Qcout<<"current_filename="<<current_filename<<endl;
 495   Qcout<<"Saving settings as "+current_filename+".sp..."<<endl;
 496 //   if(!current_settingssheet_name.isEmpty()) tableWidget->writeFile(current_settingssheet_name);
 497   if(!tableWidget->writeFile(current_filename+".sp")) cout<<"Saving settingssheet failed."<<endl;
 498
 499   cout<<"METHOD "<<ui.SmoothMethod->currentIndex()<<endl;
 500   //can't use switch case because dynamic variables seem to be forbidden inside case statements
 501   //////////////////////////////////////////////////////////////////////////////////////////////
 502   if(ui.SmoothMethod->currentIndex()==0)//vtkSmoothPolyDataFilter smoothing
 503   {
 504     //preparations
 505     QSet<int> bcs;
 506     getSelectedItems(ui.listWidget, bcs);
 507     QVector<vtkIdType> cells;
 508     getSurfaceCells(bcs, cells, grid);
 509     EG_VTKSP(vtkPolyData, pdata);
 510     addToPolyData(cells, pdata, grid);
 511     EG_VTKSP(vtkSmoothPolyDataFilter, smooth);
 512
 513 //     EG_VTKSP(vtkEgGridSmoothPolyDataFilter, smooth2);
 514 //     vtkSmartPointer<vtkEgGridSmoothPolyDataFilter> smooth2 = vtkSmartPointer<vtkEgGridSmoothPolyDataFilter>::New();
 515
 516     cout_vtkSmoothPolyDataFilter(smooth);
 517
 518     //configure vtkSmoothPolyDataFilter
 519     smooth->SetInput(pdata);
 520
 521     smooth->SetConvergence (ui.doubleSpinBox_Convergence->value());
 522     smooth->SetNumberOfIterations (ui.spinBox_NumberOfIterations->value());
 523     smooth->SetRelaxationFactor (ui.lineEdit_RelaxationFactor->text().toDouble());
 524     smooth->SetFeatureEdgeSmoothing (ui.checkBox_FeatureEdgeSmoothing->checkState());
 525     smooth->SetFeatureAngle (ui.doubleSpinBox_FeatureAngle->value());
 526     smooth->SetEdgeAngle (ui.doubleSpinBox_EdgeAngle->value());
 527     smooth->SetBoundarySmoothing (ui.checkBox_BoundarySmoothing->checkState());
 528     smooth->SetGenerateErrorScalars (ui.checkBox_GenerateErrorScalars->checkState());
 529     smooth->SetGenerateErrorVectors (ui.checkBox_GenerateErrorVectors->checkState());
 530
 531     QSet<int> bcs_Source;
 532     getSelectedItems(ui.listWidget,bcs_Source);
 533     QVector<vtkIdType> cells_Source;
 534     getSurfaceCells(bcs_Source, cells_Source, grid);
 535     EG_VTKSP(vtkPolyData, pdata_Source);
 536     addToPolyData(cells_Source, pdata_Source, grid);
 537     smooth->SetSource (pdata_Source);
 538
 539     cout_vtkSmoothPolyDataFilter(smooth);
 540
 541     //smooth
 542     smooth->Update();
 543
 544     if(ui.checkBox_GenerateErrorScalars->checkState()) DisplayErrorScalars(smooth);
 545     if(ui.checkBox_GenerateErrorVectors->checkState()) DisplayErrorVectors(smooth);
 546
 547     //copy smoothed grid to main grid
 548     EG_VTKDCN(vtkLongArray_t, node_index, pdata, "node_index");
 549     for (vtkIdType i = 0; i < smooth->GetOutput()->GetNumberOfPoints(); ++i) {
 550       vec3_t x;
 551       smooth->GetOutput()->GetPoints()->GetPoint(i, x.data());
 552       vtkIdType nodeId = node_index->GetValue(i);
 553       grid->GetPoints()->SetPoint(nodeId, x.data());
 554     };
 555     updateActors();
 556   }
 557   //////////////////////////////////////////////////////////////////////////////////////////////
 558   else if(ui.SmoothMethod->currentIndex()==1)//vtkWindowedSincPolyDataFilter smoothing
 559   {
 560     QSet<int> bcs;
 561     getSelectedItems(ui.listWidget, bcs);
 562     QVector<vtkIdType> cells;
 563     getSurfaceCells(bcs, cells, grid);
 564     EG_VTKSP(vtkPolyData, pdata);
 565     addToPolyData(cells, pdata, grid);
 566     EG_VTKSP(vtkWindowedSincPolyDataFilter, smooth);
 567
 568     cout_vtkWindowedSincPolyDataFilter(smooth);
 569
 570     smooth->SetInput(pdata);
 571
 572     smooth->SetNumberOfIterations(ui.spinBox_NumberOfIterations->value());
 573     smooth->SetPassBand(ui.doubleSpinBox_PassBand->value());
 574     smooth->SetFeatureEdgeSmoothing(ui.checkBox_FeatureEdgeSmoothing->checkState());
 575     smooth->SetFeatureAngle(ui.doubleSpinBox_FeatureAngle->value());
 576     smooth->SetEdgeAngle(ui.doubleSpinBox_EdgeAngle->value());
 577     smooth->SetBoundarySmoothing(ui.checkBox_BoundarySmoothing->checkState());
 578     smooth->SetGenerateErrorScalars(ui.checkBox_GenerateErrorScalars->checkState());
 579     smooth->SetGenerateErrorVectors(ui.checkBox_GenerateErrorVectors->checkState());
 580
 581     cout_vtkWindowedSincPolyDataFilter(smooth);
 582
 583     smooth->Update();
 584
 585     if(ui.checkBox_GenerateErrorScalars->checkState()) DisplayErrorScalars(smooth);
 586     if(ui.checkBox_GenerateErrorVectors->checkState()) DisplayErrorVectors(smooth);
 587
 588     EG_VTKDCN(vtkLongArray_t, node_index, pdata, "node_index");
 589     for (vtkIdType i = 0; i < smooth->GetOutput()->GetNumberOfPoints(); ++i) {
 590       vec3_t x;
 591       smooth->GetOutput()->GetPoints()->GetPoint(i, x.data());
 592       vtkIdType nodeId = node_index->GetValue(i);
 593       grid->GetPoints()->SetPoint(nodeId, x.data());
 594     };
 595     updateActors();
 596   }
 597   //////////////////////////////////////////////////////////////////////////////////////////////
 598   else if(ui.SmoothMethod->currentIndex()==2)//edge subdivision
 599   {
 600     cout_grid(cout,grid);
 601
 602     QSet<int> bcs;
 603     getSelectedItems(ui.listWidget, bcs);
 604
 605     int N_iter=ui.spinBox_NumberOfSubdivisions->value();
 606     for(int i_iter=0;i_iter<N_iter;i_iter++)
 607     {
 608       cout<<"i_iter="<<i_iter<<endl;
 609
 610       QVector<vtkIdType> SelectedCells;
 611       getSurfaceCells(bcs, SelectedCells, grid);
 612       QVector<vtkIdType> AllCells;
 613       getAllSurfaceCells(AllCells,grid);
 614
 615       createCellToCell(AllCells, c2c, grid);
 616
 617       int N_points=grid->GetNumberOfPoints();
 618       int N_cells=grid->GetNumberOfCells();
 619
 620       QMap< pair<vtkIdType,vtkIdType>, vtkIdType> midpoint_map;
 621 //       QMap<double, int> midpoint_map;
 622       int N_extmidpoints=0;
 623       int N_intmidpoints=0;
 624
 625       int N_newpoints=0;
 626       int N_newcells=0;
 627       vtkIdType nodeId = N_points;
 628       foreach(vtkIdType id_cell, SelectedCells)
 629       {
 630         vtkIdType type_cell = grid->GetCellType(id_cell);
 631         int N_neighbours=c2c[id_cell].size();
 632         for(int i=0;i<N_neighbours;i++)
 633         {
 634           vtkIdType id_neighbour=c2c[id_cell][i];
 635           if(id_neighbour<0)
 636           {
 637             N_extmidpoints++;
 638           }
 639           else
 640           {
 641             midpoint_map[OrderedPair(id_cell,id_neighbour)]=nodeId; nodeId++;
 642           }
 643         }
 644         if (type_cell == VTK_TRIANGLE)
 645         {
 646           N_newcells+=3;
 647         }
 648         if (type_cell == VTK_QUAD)
 649         {
 650           N_newcells+=3;
 651           N_intmidpoints++;
 652         }
 653       }
 654
 655       int N_c2cmidpoints=midpoint_map.size();
 656       N_newpoints=N_c2cmidpoints+N_extmidpoints+N_intmidpoints;
 657
 658       cout<<"N_c2cmidpoints="<<N_c2cmidpoints<<endl;
 659       cout<<"N_extmidpoints="<<N_extmidpoints<<endl;
 660       cout<<"N_intmidpoints="<<N_intmidpoints<<endl;
 661       cout<<"N_newpoints="<<N_newpoints<<endl;
 662       cout<<"N_newcells="<<N_newcells<<endl;
 663       cout<<"N_cells="<<N_cells<<endl;
 664       cout<<"N_points="<<N_points<<endl;
 665
 666       EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 667       allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 668       makeCopyNoAlloc(grid, grid_tmp);
 669
 670       EG_VTKDCC(vtkIntArray, cell_code, grid_tmp, "cell_code");
 671
 672       midpoint_map.clear();//clear midpoint_map
 673       nodeId=N_points;//reset nodeId
 674
 675       foreach(vtkIdType id_cell, SelectedCells)
 676       {
 677         vtkIdType N_pts, *pts;
 678         grid->GetCellPoints(id_cell, N_pts, pts);
 679         vtkIdType intmidpoint;
 680         vtkIdType edgemidpoints[4];
 681         vec3_t M[4];
 682
 683         vtkIdType type_cell = grid->GetCellType(id_cell);
 684         int N_neighbours=c2c[id_cell].size();
 685         vec3_t corner[4];
 686         for(int i=0;i<N_neighbours;i++)
 687         {
 688           grid->GetPoints()->GetPoint(pts[i], corner[i].data());
 689         }
 690
 691         for(int i=0;i<N_neighbours;i++)
 692         {
 693           vtkIdType id_neighbour=c2c[id_cell][i];
 694           if(id_neighbour<0)
 695           {
 696             M[i]=0.5*(corner[i]+corner[(i+1)%N_neighbours]);
 697             addPoint(grid_tmp,nodeId,M[i].data());
 698             edgemidpoints[i]=nodeId;
 699             nodeId++;
 700           }
 701           else
 702           {
 703             if(midpoint_map.contains(OrderedPair(id_cell,id_neighbour)))
 704             {
 705               //pt already exists!
 706 /*              cout<<"pt already exists!: i="<<i<<" midpoint_map[OrderedPair(id_cell,id_neighbour)]="<<midpoint_map[OrderedPair(id_cell,id_neighbour)]<<endl;*/
 707               edgemidpoints[i]=midpoint_map[OrderedPair(id_cell,id_neighbour)];
 708             }
 709             else
 710             {
 711               M[i]=0.5*(corner[i]+corner[(i+1)%N_neighbours]);
 712               addPoint(grid_tmp,nodeId,M[i].data());
 713               midpoint_map[OrderedPair(id_cell,id_neighbour)]=nodeId;
 714               edgemidpoints[i]=nodeId;
 715               nodeId++;
 716             }
 717           }
 718         }
 719         if (type_cell == VTK_TRIANGLE)
 720         {
 721           vtkIdType pts_triangle[4][3];
 722           pts_triangle[0][0]=pts[0];//A;
 723           pts_triangle[0][1]=edgemidpoints[0];
 724           pts_triangle[0][2]=edgemidpoints[2];
 725           pts_triangle[1][0]=pts[1];//B;
 726           pts_triangle[1][1]=edgemidpoints[1];
 727           pts_triangle[1][2]=edgemidpoints[0];
 728           pts_triangle[2][0]=pts[2];//C;
 729           pts_triangle[2][1]=edgemidpoints[2];
 730           pts_triangle[2][2]=edgemidpoints[1];
 731           pts_triangle[3][0]=edgemidpoints[0];
 732           pts_triangle[3][1]=edgemidpoints[1];
 733           pts_triangle[3][2]=edgemidpoints[2];
 734
 735           vtkIdType newCellId;
 736           newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[0]);
 737           cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 738           newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[1]);
 739           cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 740           newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[2]);
 741           cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 742           grid_tmp->ReplaceCell(id_cell , 3, pts_triangle[3]);
 743           cell_code->SetValue(id_cell, ui.lineEdit_BoundaryCode->text().toInt());
 744         }
 745         if (type_cell == VTK_QUAD)
 746         {
 747           vec3_t C=0.25*(corner[0]+corner[1]+corner[2]+corner[3]);
 748           addPoint(grid_tmp,nodeId,C.data());
 749           intmidpoint=nodeId;
 750           nodeId++;
 751
 752           vtkIdType pts_quad[4][4];
 753           pts_quad[0][0]=pts[0];
 754           pts_quad[0][1]=edgemidpoints[0];
 755           pts_quad[0][2]=intmidpoint;
 756           pts_quad[0][3]=edgemidpoints[3];
 757           pts_quad[1][0]=pts[1];
 758           pts_quad[1][1]=edgemidpoints[1];
 759           pts_quad[1][2]=intmidpoint;
 760           pts_quad[1][3]=edgemidpoints[0];
 761           pts_quad[2][0]=pts[2];
 762           pts_quad[2][1]=edgemidpoints[2];
 763           pts_quad[2][2]=intmidpoint;
 764           pts_quad[2][3]=edgemidpoints[1];
 765           pts_quad[3][0]=pts[3];
 766           pts_quad[3][1]=edgemidpoints[3];
 767           pts_quad[3][2]=intmidpoint;
 768           pts_quad[3][3]=edgemidpoints[2];
 769
 770           vtkIdType newCellId;
 771           newCellId = grid_tmp->InsertNextCell(VTK_QUAD,4,pts_quad[0]);
 772           cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 773           newCellId = grid_tmp->InsertNextCell(VTK_QUAD,4,pts_quad[1]);
 774           cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 775           newCellId = grid_tmp->InsertNextCell(VTK_QUAD,4,pts_quad[2]);
 776           cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 777           grid_tmp->ReplaceCell(id_cell , 4, pts_quad[3]);
 778           cell_code->SetValue(id_cell, ui.lineEdit_BoundaryCode->text().toInt());
 779         }
 780       }
 781
 782 //       cout_grid(cout,grid_tmp,true,true,true,true);
 783       cout<<"Copying..."<<endl;
 784       makeCopy(grid_tmp,grid);
 785       cout<<"Copy successful"<<endl;
 786     }
 787     cout_grid(cout,grid);
 788     updateActors();
 789   }
 790   //////////////////////////////////////////////////////////////////////////////////////////////
 791   else if(ui.SmoothMethod->currentIndex()==3)//swap triangles
 792   {
 793     cout_grid(cout,grid);
 794
 795     QSet<int> bcs;
 796     getSelectedItems(ui.listWidget, bcs);
 797
 798     QSet<int> bcs_complement=complementary_bcs(bcs,grid,cells);
 799
 800     cout<<"bcs="<<bcs<<endl;
 801     cout<<"bcs_complement="<<bcs_complement<<endl;
 802
 803     SwapTriangles swap;
 804     swap.setGrid(grid);
 805     swap.setBoundaryCodes(bcs_complement);
 806     swap();
 807
 808     cout_grid(cout,grid);
 809     updateActors();
 810   }
 811   //////////////////////////////////////////////////////////////////////////////////////////////
 812   else if(ui.SmoothMethod->currentIndex()==4)//center subdivision
 813   {
 814     cout_grid(cout,grid);
 815
 816     QSet<int> bcs;
 817     getSelectedItems(ui.listWidget, bcs);
 818     QSet<int> bcs_complement=complementary_bcs(bcs,grid,cells);
 819     cout<<"bcs="<<bcs<<endl;
 820     cout<<"bcs_complement="<<bcs_complement<<endl;
 821
 822     int N_iter=ui.spinBox_NumberOfSubdivisions->value();
 823     for(int i_iter=0;i_iter<N_iter;i_iter++)
 824     {
 825       int N_points=grid->GetNumberOfPoints();
 826       int N_cells=grid->GetNumberOfCells();
 827
 828       QVector<vtkIdType> cells;
 829       getSurfaceCells(bcs, cells, grid);
 830
 831       int N_newcells=0;
 832       foreach(int id_cell, cells)
 833       {
 834         vtkIdType type_cell = grid->GetCellType(id_cell);
 835         if (type_cell == VTK_TRIANGLE) N_newcells+=2;
 836         if (type_cell == VTK_QUAD) N_newcells+=3;
 837       }
 838
 839       cout<<"N_newcells="<<N_newcells<<endl;
 840
 841       EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 842       allocateGrid(grid_tmp,N_cells+N_newcells,N_points+1*N_cells);
 843       makeCopyNoAlloc(grid, grid_tmp);
 844       EG_VTKDCC(vtkIntArray, cell_code, grid_tmp, "cell_code");
 845
 846       vtkIdType nodeId=N_points;
 847       foreach(int id_cell, cells)
 848       {
 849         vtkIdType N_pts, *pts;
 850         grid->GetCellPoints(id_cell, N_pts, pts);
 851         vec3_t C(0,0,0);
 852
 853         vtkIdType type_cell = grid->GetCellType(id_cell);
 854         int N_neighbours=N_pts;
 855         cout<<"N_neighbours="<<N_neighbours<<endl;
 856         vec3_t corner[4];
 857         vtkIdType pts_triangle[4][3];
 858         for(int i=0;i<N_neighbours;i++)
 859         {
 860           grid->GetPoints()->GetPoint(pts[i], corner[i].data());
 861           C+=(1/(double)N_neighbours)*corner[i];
 862         }
 863         addPoint(grid_tmp,nodeId,C.data());
 864         vtkIdType intmidpoint=nodeId;
 865         nodeId++;
 866
 867         for(int i=0;i<N_neighbours;i++)
 868         {
 869           pts_triangle[i][0]=pts[i];
 870           pts_triangle[i][1]=pts[(i+1)%N_neighbours];
 871           pts_triangle[i][2]=intmidpoint;
 872           if(i==0)
 873           {
 874             grid_tmp->ReplaceCell(id_cell , 3, pts_triangle[0]);
 875             cell_code->SetValue(id_cell, ui.lineEdit_BoundaryCode->text().toInt());
 876           }
 877           else
 878           {
 879             vtkIdType newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[i]);
 880             cell_code->SetValue(newCellId, ui.lineEdit_BoundaryCode->text().toInt());
 881           }
 882         }
 883       }
 884
 885       makeCopy(grid_tmp,grid);
 886     }
 887     cout_grid(cout,grid,true,true,true,true);
 888     updateActors();
 889   }
 890   //////////////////////////////////////////////////////////////////////////////////////////////
 891   else if(ui.SmoothMethod->currentIndex()==5)//boundary refinement
 892   {
 893     cout_grid(cout,grid);
 894
 895     QSet<int> bcs;
 896     getSelectedItems(ui.listWidget, bcs);
 897     QSet<int> bcs_complement=complementary_bcs(bcs,grid,cells);
 898     cout<<"bcs="<<bcs<<endl;
 899     cout<<"bcs_complement="<<bcs_complement<<endl;
 900
 901     int N_iter=ui.spinBox_NumberOfIterations->value();
 902     for(int i_iter=0;i_iter<N_iter;i_iter++){
 903       int N_points=grid->GetNumberOfPoints();
 904       int N_cells=grid->GetNumberOfCells();
 905
 906       QVector<vtkIdType> SelectedCells;
 907       getSurfaceCells(bcs, SelectedCells, grid);
 908       QVector<vtkIdType> AllCells;
 909       getAllSurfaceCells(AllCells,grid);
 910
 911       EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 912
 913 //       QSet <int> cells_to_split;
 914       QVector <stencil_t> StencilVector;
 915
 916       QMap <vtkIdType,bool> marked;
 917
 918       createCellMapping(AllCells, _cells, grid);
 919       createCellToCell(AllCells, c2c, grid);
 920       setCells(AllCells);
 921
 922       cout<<"AllCells.size()="<<AllCells.size()<<endl;
 923       cout<<"SelectedCells.size()="<<SelectedCells.size()<<endl;
 924
 925       int N_newcells=0;
 926       int N_newpoints=0;
 927       foreach(vtkIdType id_cell, SelectedCells)
 928       {
 929         cout<<"==>id_cell="<<id_cell<<endl;
 930         int bc0=cell_code->GetValue(id_cell);
 931         if(!marked[id_cell])
 932         {
 933           vtkIdType N_pts, *pts;
 934           grid->GetCellPoints(id_cell, N_pts, pts);
 935           int count=0;
 936           for(int i=0;i<N_pts;i++)
 937           {
 938             int bc1=cell_code->GetValue(c2c[id_cell][i]);
 939             if(bc0!=bc1) count++;
 940           }
 941           if(count>0)//cell is near at least one neighbour with different cell code
 942           {
 943             int SideToSplit = getLongestSide(id_cell,grid);
 944             cout<<"SideToSplit="<<SideToSplit<<endl;
 945             cout<<"c2c[id_cell][SideToSplit]="<<c2c[id_cell][SideToSplit]<<endl;
 946             for(int i=0;i<3;i++) cout<<"c2c[id_cell]["<<i<<"]="<<c2c[id_cell][i]<<endl;
 947             stencil_t S=getStencil(id_cell,SideToSplit);
 948             if(S.valid){//there is a neighbour cell
 949               if(!marked[S.id_cell2])
 950               {
 951 /*                cells_to_split.insert(S.id_cell1);
 952                 cells_to_split.insert(S.id_cell2);*/
 953                 cout<<"marked["<<S.id_cell1<<"]=true;"<<endl;
 954                 cout<<"marked["<<S.id_cell2<<"]=true;"<<endl;
 955                 marked[S.id_cell1]=true;
 956                 marked[S.id_cell2]=true;
 957                 StencilVector.push_back(S);
 958                 N_newpoints++;
 959                 N_newcells+=2;
 960               }
 961             }
 962             else{//there is no neighbour cell
 963 //               cells_to_split.insert(S.id_cell1);
 964               cout<<"marked["<<S.id_cell1<<"]=true;"<<endl;
 965               marked[S.id_cell1]=true;
 966               StencilVector.push_back(S);
 967               N_newpoints++;
 968               N_newcells+=1;
 969             }
 970           }
 971         }
 972       }
 973
 974 /*      cout<<"cells_to_split.size()="<<cells_to_split.size()<<endl;
 975       cout<<cells_to_split<<endl;*/
 976       cout<<"StencilVector.size()="<<StencilVector.size()<<endl;
 977       cout<<"N_newpoints="<<N_newpoints<<endl;
 978       cout<<"N_newcells="<<N_newcells<<endl;
 979
 980       EG_VTKSP(vtkUnstructuredGrid,grid_tmp);
 981       allocateGrid(grid_tmp,N_cells+N_newcells,N_points+N_newpoints);
 982       makeCopyNoAlloc(grid, grid_tmp);
 983       EG_VTKDCC(vtkIntArray, cell_code_tmp, grid_tmp, "cell_code");
 984
 985       vtkIdType nodeId=N_points;
 986       foreach(stencil_t S,StencilVector)
 987       {
 988         cout<<S<<endl;
 989         vtkIdType N_pts, *pts;
 990         vec3_t A,B;
 991         grid_tmp->GetPoint(S.p[1],A.data());
 992         grid_tmp->GetPoint(S.p[3],B.data());
 993         vec3_t M=0.5*(A+B);
 994         addPoint(grid_tmp,nodeId,M.data());
 995
 996         vtkIdType pts_triangle[4][3];
 997
 998         if(S.valid){//there is a neighbour cell
 999           for(int i=0;i<4;i++)
1000           {
1001             pts_triangle[i][0]=S.p[i];
1002             pts_triangle[i][1]=S.p[(i+1)%4];
1003             pts_triangle[i][2]=nodeId;
1004           }
1005
1006           int bc1=cell_code_tmp->GetValue(S.id_cell1);
1007           int bc2=cell_code_tmp->GetValue(S.id_cell2);
1008
1009           grid_tmp->ReplaceCell(S.id_cell1 , 3, pts_triangle[0]);
1010           cell_code_tmp->SetValue(S.id_cell1, bc1);
1011
1012           grid_tmp->ReplaceCell(S.id_cell2 , 3, pts_triangle[1]);
1013           cell_code_tmp->SetValue(S.id_cell2, bc2);
1014
1015           vtkIdType newCellId;
1016           newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[2]);
1017           cell_code_tmp->SetValue(newCellId, bc2);
1018           newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[3]);
1019           cell_code_tmp->SetValue(newCellId, bc1);
1020         }
1021         else{//there is no neighbour cell
1022           pts_triangle[0][0]=S.p[0];
1023           pts_triangle[0][1]=S.p[1];
1024           pts_triangle[0][2]=nodeId;
1025           pts_triangle[3][0]=S.p[3];
1026           pts_triangle[3][1]=S.p[0];
1027           pts_triangle[3][2]=nodeId;
1028
1029           int bc1=cell_code_tmp->GetValue(S.id_cell1);
1030
1031           grid_tmp->ReplaceCell(S.id_cell1 , 3, pts_triangle[0]);
1032           cell_code_tmp->SetValue(S.id_cell1, bc1);
1033
1034           vtkIdType newCellId;
1035           newCellId = grid_tmp->InsertNextCell(VTK_TRIANGLE,3,pts_triangle[3]);
1036           cell_code_tmp->SetValue(newCellId, bc1);
1037         }
1038
1039         nodeId++;
1040       }
1041
1042       makeCopy(grid_tmp,grid);
1043     }//end of i_iter loop
1044 //     cout_grid(cout,grid,true,true,true,true);
1045     updateActors();
1046   }
1047   //////////////////////////////////////////////////////////////////////////////////////////////
1048   else if(ui.SmoothMethod->currentIndex()==6)//Laplacian smoothing
1049   {
1050     QSet<int> bcs;
1051     getSelectedItems(ui.listWidget, bcs);
1052
1053     LaplaceSmoother toto;
1054     toto.SetInput(bcs,grid);
1055     toto.SetNumberOfIterations(ui.spinBox_NumberOfIterations->value());
1056
1057     toto();
1058
1059     updateActors();
1060   }
1061   //////////////////////////////////////////////////////////////////////////////////////////////
1062   else if(ui.SmoothMethod->currentIndex()==7)//VertexAvgDist test
1063   {
1064     cout_grid(cout,grid);
1065
1066     QSet<int> bcs;
1067     getSelectedItems(ui.listWidget, bcs);
1068     QSet<int> bcs_complement=complementary_bcs(bcs,grid,cells);
1069     cout<<"bcs="<<bcs<<endl;
1070     cout<<"bcs_complement="<<bcs_complement<<endl;
1071
1072     int N_points=grid->GetNumberOfPoints();
1073     int N_cells=grid->GetNumberOfCells();
1074
1075     QVector<vtkIdType> SelectedCells;
1076     getSurfaceCells(bcs, SelectedCells, grid);
1077     QVector<vtkIdType> AllCells;
1078     getAllSurfaceCells(AllCells,grid);
1079
1080     QSet <vtkIdType> SelectedNodes;
1081     getSurfaceNodes(bcs,SelectedNodes,grid);
1082     createNodeToNode(cells, nodes, _nodes, n2n, grid);
1083
1084     foreach(vtkIdType node,SelectedNodes)
1085     {
1086       cout<<"node="<<node<<" VertexAvgDist="<<CurrentVertexAvgDist(node,n2n,grid)<<endl;
1087     }
1088
1089     updateActors();
1090   }
1091   //////////////////////////////////////////////////////////////////////////////////////////////
1092   else if(ui.SmoothMethod->currentIndex()==8)//Create mesh density map
1093   {
1094     cout_grid(cout,grid);
1095
1096     QSet<int> bcs;
1097     getSelectedItems(ui.listWidget, bcs);
1098     QSet<int> bcs_complement=complementary_bcs(bcs,grid,cells);
1099     cout<<"bcs="<<bcs<<endl;
1100     cout<<"bcs_complement="<<bcs_complement<<endl;
1101
1102     int N_points=grid->GetNumberOfPoints();
1103     int N_cells=grid->GetNumberOfCells();
1104
1105     QVector<vtkIdType> SelectedCells;
1106     getSurfaceCells(bcs, SelectedCells, grid);
1107     QVector<vtkIdType> AllCells;
1108     getAllSurfaceCells(AllCells,grid);
1109
1110     EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
1111     EG_VTKDCN(vtkDoubleArray, node_meshdensity, grid, "node_meshdensity");
1112
1113     QSet <vtkIdType> SelectedNodes;
1114     getSurfaceNodes(bcs,SelectedNodes,grid);
1115     createNodeToNode(cells, nodes, _nodes, n2n, grid);
1116
1117     foreach(vtkIdType node,SelectedNodes)
1118     {
1119       double L=CurrentVertexAvgDist(node,n2n,grid);
1120       double D=1./L;
1121       cout<<"node="<<node<<" VertexAvgDist="<<L<<" Net density="<<D<<endl;
1122       node_meshdensity->SetValue(node, D);
1123     }
1124
1125     int N_iter=ui.spinBox_maxiter_density->value();
1126     for(int i_iter=0;i_iter<N_iter;i_iter++)
1127     {
1128       foreach(vtkIdType node,SelectedNodes)
1129       {
1130         double D=DesiredMeshDensity(node,n2n,grid);
1131         double L=1./D;
1132         cout<<"node="<<node<<" VertexAvgDist="<<L<<" Net density="<<D<<endl;
1133         node_meshdensity->SetValue(node, D);
1134       }
1135     }
1136
1137     updateActors();
1138   }
1139   //////////////////////////////////////////////////////////////////////////////////////////////
1140   else if(ui.SmoothMethod->currentIndex()==9)//vtkWindowedSincPolyDataFilter smoothing
1141   {
1142     QSet<int> bcs;
1143     getSelectedItems(ui.listWidget, bcs);
1144     QVector<vtkIdType> cells;
1145     getSurfaceCells(bcs, cells, grid);
1146     EG_VTKSP(vtkPolyData, input);
1147     addToPolyData(cells, input, grid);
1148
1149     EG_VTKSP(vtkSmoothPolyDataFilter, smooth);
1150     smooth->SetInput(input);
1151
1152
1153     int j, k;
1154     vtkIdType npts = 0;
1155     vtkIdType *pts = 0;
1156
1157     vtkCellArray *inVerts, *inLines, *inPolys, *inStrips;
1158
1159     cout<<"input->GetVerts()="<<input->GetVerts()<<endl;
1160     cout<<"input->GetVerts()->GetSize()="<<input->GetVerts()->GetSize()<<endl;
1161     inVerts=input->GetVerts();
1162     inVerts->InitTraversal();
1163     cout<<"inVerts->GetSize()="<<inVerts->GetSize()<<endl;
1164     cout<<"inVerts->GetNextCell(npts,pts)="<<inVerts->GetNextCell(npts,pts)<<endl;
1165     cout<<"inVerts->GetNumberOfCells()="<<inVerts->GetNumberOfCells()<<endl;
1166
1167     cout<<"input->GetVerts()->GetSize()="<<input->GetVerts()->GetSize()<<endl;
1168     cout<<"input->GetVerts()->GetNextCell(npts,pts)="<<input->GetVerts()->GetNextCell(npts,pts)<<endl;
1169     cout<<"input->GetVerts()->GetNumberOfCells()="<<input->GetVerts()->GetNumberOfCells()<<endl;
1170
1171   // check vertices first. Vertices are never smoothed_--------------
1172   for (inVerts=input->GetVerts(), inVerts->InitTraversal();
1173   inVerts->GetNextCell(npts,pts); )
1174     {
1175         cout<<"npts="<<npts<<endl;
1176     for (j=0; j<npts; j++)
1177       {
1178         cout<<"pts["<<j<<"]="<<pts[j]<<endl;
1179       }
1180     }
1181
1182     updateActors();
1183   }
1184   //////////////////////////////////////////////////////////////////////////////////////////////
1185   else if(ui.SmoothMethod->currentIndex()==10)// super smoothing
1186   {
1187     QSet<int> bcs;
1188     getSelectedItems(ui.listWidget, bcs);
1189
1190     QVector <VertexMeshDensity> VMDvector=GetSet();
1191
1192     CreateSpecialMapping toto;
1193
1194     toto.SetInput(bcs,grid);
1195     toto.SetVertexMeshDensityVector(VMDvector);
1196     toto.SetConvergence (ui.doubleSpinBox_Convergence->value());
1197     toto.SetNumberOfIterations (ui.spinBox_NumberOfIterations->value());
1198     toto.SetRelaxationFactor (ui.lineEdit_RelaxationFactor->text().toDouble());
1199     toto.SetFeatureEdgeSmoothing (ui.checkBox_FeatureEdgeSmoothing->checkState());
1200     toto.SetFeatureAngle (ui.doubleSpinBox_FeatureAngle->value());
1201     toto.SetEdgeAngle (ui.doubleSpinBox_EdgeAngle->value());
1202     toto.SetBoundarySmoothing (ui.checkBox_BoundarySmoothing->checkState());
1203     toto.SetGenerateErrorScalars (ui.checkBox_GenerateErrorScalars->checkState());
1204     toto.SetGenerateErrorVectors (ui.checkBox_GenerateErrorVectors->checkState());
1205
1206 /*    toto.Set_SV_value(ui.doubleSpinBox_VTK_SIMPLE_VERTEX->value());
1207     toto.Set_FV_value(ui.doubleSpinBox_VTK_FIXED_VERTEX->value());
1208     toto.Set_FEV_value(ui.doubleSpinBox_VTK_FEATURE_EDGE_VERTEX->value());
1209     toto.Set_BEV_value(ui.doubleSpinBox_VTK_BOUNDARY_EDGE_VERTEX->value());*/
1210
1211     toto.SetConvergence_meshdensity(ui.doubleSpinBox_Convergence_meshdensity->value());
1212
1213     toto.Set_insert_FP(ui.checkBox_insert_FP->checkState());
1214     toto.Set_insert_EP(ui.checkBox_insert_EP->checkState());
1215     toto.Set_remove_FP(ui.checkBox_remove_FP->checkState());
1216     toto.Set_remove_EP(ui.checkBox_remove_EP->checkState());
1217     toto.DoSwap=ui.checkBox_Swap->checkState();
1218     toto.DoLaplaceSmoothing=ui.checkBox_LaplaceSmoothing->checkState();
1219
1220     toto.N_SmoothIterations=ui.spinBox_NumberOfSmoothIterations->value();
1221     toto.maxiter_density=ui.spinBox_maxiter_density->value();
1222     DebugLevel=ui.spinBox_DebugLevel->value();
1223
1224     toto.Process();
1225
1226     updateActors();
1227   }
1228   //////////////////////////////////////////////////////////////////////////////////////////////
1229   else if(ui.SmoothMethod->currentIndex()==11)// Update current mesh density + node types
1230   {
1231     QSet<int> bcs;
1232     getSelectedItems(ui.listWidget, bcs);
1233     CreateSpecialMapping toto;
1234     toto.SetInput(bcs,grid);
1235
1236     SetConvergence(ui.doubleSpinBox_Convergence->value());
1237     SetFeatureEdgeSmoothing(ui.checkBox_FeatureEdgeSmoothing->checkState());
1238     SetFeatureAngle(ui.doubleSpinBox_FeatureAngle->value());
1239     SetEdgeAngle(ui.doubleSpinBox_EdgeAngle->value());
1240     SetBoundarySmoothing(ui.checkBox_BoundarySmoothing->checkState());
1241
1242     UpdateMeshDensity();
1243     UpdateNodeType();
1244     updateActors();
1245   }
1246   //////////////////////////////////////////////////////////////////////////////////////////////
1247   else if(ui.SmoothMethod->currentIndex()==12)// Delete all possible points
1248   {
1249     QSet<int> bcs;
1250     getSelectedItems(ui.listWidget, bcs);
1251     CreateSpecialMapping toto;
1252
1253     SetConvergence(ui.doubleSpinBox_Convergence->value());
1254     SetFeatureEdgeSmoothing(ui.checkBox_FeatureEdgeSmoothing->checkState());
1255     SetFeatureAngle(ui.doubleSpinBox_FeatureAngle->value());
1256     SetEdgeAngle(ui.doubleSpinBox_EdgeAngle->value());
1257     SetBoundarySmoothing(ui.checkBox_BoundarySmoothing->checkState());
1258
1259     int N_newpoints;
1260     int N_newcells;
1261
1262     vtkIdType N_points=grid->GetNumberOfPoints();
1263     vtkIdType N_cells=grid->GetNumberOfCells();
1264
1265     bool Global_DelResult=true;
1266     while(Global_DelResult)
1267     {
1268       Global_DelResult=false;
1269       vtkIdType DeadNode=0;
1270       while(DeadNode<grid->GetNumberOfPoints())
1271       {
1272         bool Local_DelResult=true;
1273         while(Local_DelResult)
1274         {
1275           Local_DelResult=DeletePoint_2(grid,DeadNode,N_newpoints,N_newcells);
1276           if(Local_DelResult) Global_DelResult=true;
1277         }
1278         DeadNode++;
1279       }
1280     }
1281
1282   }
1283   //////////////////////////////////////////////////////////////////////////////////////////////
1284   else
1285   {
1286     cout<<"UNKNOWN METHOD"<<endl;
1287     QSet<int> bcs;
1288     getSelectedItems(ui.listWidget, bcs);
1289     QVector<vtkIdType> cells;
1290     getSurfaceCells(bcs, cells, grid);
1291     setCells(cells);
1292     cout<<"cells="<<cells<<endl;
1293     cout<<"_cells="<<_cells<<endl;
1294     cout<<"nodes="<<nodes<<endl;
1295     cout<<"_nodes="<<_nodes<<endl;
1296 /*    QVector< QSet< int > >    n2c;
1297     QVector< QSet< int > >      n2n;
1298     QVector< QVector< int > >   c2c;*/
1299     cout<<"n2c="<<n2c<<endl;
1300     cout<<"n2n="<<n2n<<endl;
1301     cout<<"c2c="<<c2c<<endl;
1302   }
1303 };