egvtkobject.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 "egvtkobject.h"
  24 #include "guimainwindow.h"
  25
  26 #include <vtkCellData.h>
  27 #include <vtkPointData.h>
  28 #include <vtkCellLinks.h>
  29 #include <vtkCellType.h>
  30 #include <vtkIdList.h>
  31 #include <vtkCell.h>
  32 #include <vtkCharArray.h>
  33
  34 void EgVtkObject::computeNormals
  35 (
  36   QVector<vec3_t>     &cell_normals,
  37   QVector<vec3_t>     &node_normals,
  38   QVector<vtkIdType>  &cells,
  39   QVector<vtkIdType>  &nodes,
  40   vtkUnstructuredGrid *grid
  41 )
  42 {
  43   using namespace GeometryTools;
  44
  45   cell_normals.resize(cells.size());
  46   node_normals.fill(vec3_t(0,0,0), nodes.size());
  47   QVector<int> g2s;
  48   createNodeMapping(nodes, g2s, grid);
  49   for (int i_cell = 0; i_cell < cells.count(); ++i_cell) {
  50     vtkIdType id_cell = cells[i_cell];
  51     vtkIdType *pts;
  52     vtkIdType npts;
  53     grid->GetCellPoints(id_cell, npts, pts);
  54     cell_normals[i_cell] = cellNormal(grid, id_cell);
  55     cell_normals[i_cell].normalise();
  56     for (int i_pts = 0; i_pts < npts; ++i_pts) {
  57       if (g2s[pts[i_pts]] != -1) {
  58         node_normals[g2s[pts[i_pts]]] += cell_normals[i_cell];
  59       };
  60     };
  61   };
  62   for (int i_node = 0; i_node < nodes.count(); ++i_node) {
  63     node_normals[i_node].normalise();
  64     //cout << node_normals[i_node] << endl;
  65   };
  66 };
  67
  68 void EgVtkObject::createNodeMapping
  69 (
  70   QVector<vtkIdType>  &nodes,
  71   QVector<int>        &_nodes,
  72   vtkUnstructuredGrid *grid
  73 )
  74 {
  75   _nodes.fill(-1,grid->GetNumberOfPoints());
  76   vtkIdType id;
  77   int i_sub = 0;
  78   foreach(id,nodes) {
  79     _nodes[id] = i_sub;
  80     ++i_sub;
  81   };
  82 };
  83
  84 void EgVtkObject::createCellMapping
  85 (
  86   QVector<vtkIdType>  &cells,
  87   QVector<int>        &_cells,
  88   vtkUnstructuredGrid *grid
  89 )
  90 {
  91   _cells.fill(-1,grid->GetNumberOfCells());
  92   vtkIdType id;
  93   int i_sub = 0;
  94   foreach(id,cells) {
  95     _cells[id] = i_sub;
  96     ++i_sub;
  97   };
  98 };
  99
 100 void EgVtkObject::createNodeToBcMapping
 101 (
 102   QVector<QSet<int> > &bcs,
 103   vtkUnstructuredGrid *grid
 104 )
 105 {
 106   bcs.fill(QSet<int>(), grid->GetNumberOfPoints());
 107   grid->BuildLinks();
 108   EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 109   for (vtkIdType nodeId = 0; nodeId < grid->GetNumberOfPoints(); ++nodeId) {
 110     int Ncells = grid->GetCellLinks()->GetNcells(nodeId);
 111     for (int i = 0; i < Ncells; ++i) {
 112       vtkIdType id_cell = grid->GetCellLinks()->GetCells(nodeId)[i];
 113       vtkIdType ct = grid->GetCellType(id_cell);
 114       if ((ct == VTK_TRIANGLE) || (ct = VTK_QUAD)) {
 115         if (cell_code->GetValue(id_cell) > 0) {
 116           bcs[nodeId].insert(cell_code->GetValue(id_cell));
 117         };
 118       };
 119     };
 120   };
 121 };
 122
 123 void EgVtkObject::createNodeToCell
 124 (
 125   QVector<vtkIdType>  &cells,
 126   QVector<vtkIdType>  &nodes,
 127   QVector<int>        &_nodes,
 128   QVector<QSet<int> > &n2c,
 129   vtkUnstructuredGrid *grid
 130 )
 131 {
 132   n2c.fill(QSet<int>(), nodes.size());
 133   for (vtkIdType i_cells = 0; i_cells < cells.size(); ++i_cells) {
 134     vtkIdType *pts;
 135     vtkIdType  Npts;
 136     grid->GetCellPoints(cells[i_cells], Npts, pts);
 137     for (int i_pts = 0; i_pts < Npts; ++i_pts) {
 138       n2c[_nodes[pts[i_pts]]].insert(i_cells);
 139     };
 140   };
 141 };
 142
 143 void EgVtkObject::addToN2N
 144 (
 145   QVector<QSet<int> > &n2n,
 146   int                  n1,
 147   int                  n2
 148 )
 149 {
 150   n2n[n1].insert(n2);
 151   n2n[n2].insert(n1);
 152 };
 153
 154 void EgVtkObject::createNodeToNode
 155 (
 156   QVector<vtkIdType>  &cells,
 157   QVector<vtkIdType>  &nodes,
 158   QVector<int>        &_nodes,
 159   QVector<QSet<int> > &n2n,
 160   vtkUnstructuredGrid *grid
 161 )
 162 {
 163   n2n.fill(QSet<int>(), nodes.size());
 164   foreach (vtkIdType id_cell, cells) {
 165     vtkIdType *pts;
 166     vtkIdType  Npts;
 167     grid->GetCellPoints(id_cell, Npts, pts);
 168     vector<int> n(Npts);
 169     for (int i = 0; i < Npts; ++i) {
 170       n[i] = _nodes[pts[i]];
 171     };
 172     vtkIdType cellType = grid->GetCellType(id_cell);
 173     if (cellType == VTK_TRIANGLE) {
 174       addToN2N(n2n, n[0], n[1]);
 175       addToN2N(n2n, n[1], n[2]);
 176       addToN2N(n2n, n[2], n[0]);
 177     } else if (cellType == VTK_QUAD) {
 178       addToN2N(n2n, n[0], n[1]);
 179       addToN2N(n2n, n[1], n[2]);
 180       addToN2N(n2n, n[2], n[3]);
 181       addToN2N(n2n, n[3], n[0]);
 182     } else if (cellType == VTK_TETRA) {
 183       addToN2N(n2n, n[0], n[1]);
 184       addToN2N(n2n, n[0], n[2]);
 185       addToN2N(n2n, n[0], n[3]);
 186       addToN2N(n2n, n[1], n[2]);
 187       addToN2N(n2n, n[1], n[3]);
 188       addToN2N(n2n, n[2], n[3]);
 189     } else if (cellType == VTK_PYRAMID) {
 190       addToN2N(n2n, n[0], n[1]);
 191       addToN2N(n2n, n[0], n[3]);
 192       addToN2N(n2n, n[0], n[4]);
 193       addToN2N(n2n, n[1], n[2]);
 194       addToN2N(n2n, n[1], n[4]);
 195       addToN2N(n2n, n[2], n[3]);
 196       addToN2N(n2n, n[2], n[4]);
 197       addToN2N(n2n, n[3], n[4]);
 198     } else if (cellType == VTK_WEDGE) {
 199       addToN2N(n2n, n[0], n[1]);
 200       addToN2N(n2n, n[0], n[2]);
 201       addToN2N(n2n, n[0], n[3]);
 202       addToN2N(n2n, n[1], n[2]);
 203       addToN2N(n2n, n[1], n[4]);
 204       addToN2N(n2n, n[2], n[5]);
 205       addToN2N(n2n, n[3], n[4]);
 206       addToN2N(n2n, n[3], n[5]);
 207       addToN2N(n2n, n[4], n[5]);
 208     } else if (cellType == VTK_HEXAHEDRON) {
 209       addToN2N(n2n, n[0], n[1]);
 210       addToN2N(n2n, n[0], n[3]);
 211       addToN2N(n2n, n[0], n[4]);
 212       addToN2N(n2n, n[1], n[2]);
 213       addToN2N(n2n, n[1], n[5]);
 214       addToN2N(n2n, n[2], n[3]);
 215       addToN2N(n2n, n[2], n[6]);
 216       addToN2N(n2n, n[3], n[7]);
 217       addToN2N(n2n, n[4], n[5]);
 218       addToN2N(n2n, n[4], n[7]);
 219       addToN2N(n2n, n[5], n[6]);
 220       addToN2N(n2n, n[6], n[7]);
 221     };
 222   };
 223 };
 224
 225 void EgVtkObject::getAllCells
 226 (
 227   QVector<vtkIdType>  &cells,
 228   vtkUnstructuredGrid *grid
 229 )
 230 {
 231   int N = 0;
 232   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 233     ++N;
 234   };
 235   cells.resize(N);
 236   N = 0;
 237   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 238     cells[N] = id_cell;
 239     ++N;
 240   };
 241 };
 242
 243 void EgVtkObject::getAllCellsOfType
 244 (
 245   vtkIdType            type,
 246   QVector<vtkIdType>  &cells,
 247   vtkUnstructuredGrid *grid
 248 )
 249 {
 250   int N = 0;
 251   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 252     if (grid->GetCellType(id_cell) == type) {
 253       ++N;
 254     };
 255   };
 256   cells.resize(N);
 257   N = 0;
 258   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 259     if (grid->GetCellType(id_cell) == type) {
 260       cells[N] = id_cell;
 261       ++N;
 262     };
 263   };
 264 };
 265
 266
 267 void EgVtkObject::getAllVolumeCells
 268 (
 269   QVector<vtkIdType>  &cells,
 270   vtkUnstructuredGrid *grid
 271 )
 272 {
 273   int N = 0;
 274   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 275     if (isVolume(id_cell, grid)) {
 276       ++N;
 277     };
 278   };
 279   cells.resize(N);
 280   N = 0;
 281   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 282     if (isVolume(id_cell, grid)) {
 283       cells[N] = id_cell;
 284       ++N;
 285     };
 286   };
 287 };
 288
 289 void EgVtkObject::getAllSurfaceCells
 290 (
 291   QVector<vtkIdType>  &cells,
 292   vtkUnstructuredGrid *grid
 293 )
 294 {
 295   int N = 0;
 296   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 297     if (isSurface(id_cell, grid)) {
 298       ++N;
 299     };
 300   };
 301   cells.resize(N);
 302   N = 0;
 303   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 304     if (isSurface(id_cell, grid)) {
 305       cells[N] = id_cell;
 306       ++N;
 307     };
 308   };
 309 };
 310
 311 void EgVtkObject::getSurfaceCells
 312 (
 313   QSet<int>           &bcs,
 314   QVector<vtkIdType>  &cells,
 315   vtkUnstructuredGrid *grid
 316 )
 317 {
 318   int N = 0;
 319   EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 320   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 321     if (isSurface(id_cell, grid)) {
 322       if (bcs.contains(cell_code->GetValue(id_cell))) {
 323         ++N;
 324       };
 325     };
 326   };
 327   cells.resize(N);
 328   N = 0;
 329   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 330     if (isSurface(id_cell, grid)) {
 331       if (bcs.contains(cell_code->GetValue(id_cell))) {
 332         cells[N] = id_cell;
 333         ++N;
 334       };
 335     };
 336   };
 337 };
 338
 339 void EgVtkObject::getSurfaceNodes
 340 (
 341   QSet<int>           &bcs,
 342   QSet <vtkIdType> &SelectedNodes,
 343   vtkUnstructuredGrid *grid
 344 )
 345 {
 346   QVector<vtkIdType> SelectedCells;
 347   getSurfaceCells(bcs, SelectedCells, grid);
 348
 349   SelectedNodes.clear();
 350   foreach(vtkIdType id_cell, SelectedCells)
 351   {
 352     vtkIdType N_pts, *pts;
 353     grid->GetCellPoints(id_cell, N_pts, pts);
 354     for(int i=0;i<N_pts;i++)
 355     {
 356       SelectedNodes.insert(pts[i]);
 357     }
 358   }
 359 };
 360
 361 void EgVtkObject::addToC2C
 362 (
 363   vtkIdType               id_cell,
 364   QVector<int>           &_cells,
 365   QVector<QVector<int> > &c2c,
 366   int                     j,
 367   vtkIdList              *nds,
 368   vtkIdList              *cls,
 369   vtkUnstructuredGrid    *grid
 370 )
 371 {
 372   c2c[_cells[id_cell]][j] = -1;
 373   grid->GetCellNeighbors(id_cell, nds, cls);
 374   for (int i = 0; i < cls->GetNumberOfIds(); ++i) {
 375     if (cls->GetId(i) != id_cell) {
 376       if (_cells[cls->GetId(i)] != -1) {
 377         c2c[_cells[id_cell]][j] = _cells[cls->GetId(i)];
 378       };
 379     };
 380   };
 381 };
 382
 383
 384 void EgVtkObject::createCellToCell
 385 (
 386   QVector<vtkIdType>     &cells,
 387   QVector<QVector<int> > &c2c,
 388   vtkUnstructuredGrid    *grid
 389 )
 390 {
 391   // GetCellNeighbors(vtkIdType id_cell, vtkIdList *ptIds, vtkIdList *id_cells)
 392   grid->BuildLinks();
 393   QVector<int> _cells;
 394   createCellMapping(cells, _cells, grid);
 395   c2c.fill(QVector<int>(), cells.size());
 396   EG_VTKSP(vtkIdList, nds);
 397   EG_VTKSP(vtkIdList, cls);
 398   for (int i = 0; i < cells.size(); ++i) {
 399     vtkIdType id_cell = cells[i];
 400     vtkIdType *pts;
 401     vtkIdType  Npts;
 402     grid->GetCellPoints(id_cell, Npts, pts);
 403     if (grid->GetCellType(id_cell) == VTK_TRIANGLE) {
 404       c2c[i].resize(3);
 405       nds->Reset();
 406       nds->InsertNextId(pts[0]);
 407       nds->InsertNextId(pts[1]);
 408       addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
 409       nds->Reset();
 410       nds->InsertNextId(pts[1]);
 411       nds->InsertNextId(pts[2]);
 412       addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
 413       nds->Reset();
 414       nds->InsertNextId(pts[2]);
 415       nds->InsertNextId(pts[0]);
 416       addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
 417     } else if (grid->GetCellType(id_cell) == VTK_QUAD) {
 418       c2c[i].resize(4);
 419       nds->Reset();
 420       nds->InsertNextId(pts[0]);
 421       nds->InsertNextId(pts[1]);
 422       addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
 423       nds->Reset();
 424       nds->InsertNextId(pts[1]);
 425       nds->InsertNextId(pts[2]);
 426       addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
 427       nds->Reset();
 428       nds->InsertNextId(pts[2]);
 429       nds->InsertNextId(pts[3]);
 430       addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
 431       nds->Reset();
 432       nds->InsertNextId(pts[3]);
 433       nds->InsertNextId(pts[0]);
 434       addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
 435     } else if (grid->GetCellType(id_cell) == VTK_TETRA) {
 436       c2c[i].resize(4);
 437       nds->Reset();
 438       nds->InsertNextId(pts[0]);
 439       nds->InsertNextId(pts[1]);
 440       nds->InsertNextId(pts[2]);
 441       addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
 442       nds->Reset();
 443       nds->InsertNextId(pts[0]);
 444       nds->InsertNextId(pts[1]);
 445       nds->InsertNextId(pts[3]);
 446       addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
 447       nds->Reset();
 448       nds->InsertNextId(pts[0]);
 449       nds->InsertNextId(pts[3]);
 450       nds->InsertNextId(pts[2]);
 451       addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
 452       nds->Reset();
 453       nds->InsertNextId(pts[1]);
 454       nds->InsertNextId(pts[2]);
 455       nds->InsertNextId(pts[3]);
 456       addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
 457     } else if (grid->GetCellType(id_cell) == VTK_PYRAMID) {
 458       c2c[i].resize(5);
 459       nds->Reset();
 460       nds->InsertNextId(pts[0]);
 461       nds->InsertNextId(pts[1]);
 462       nds->InsertNextId(pts[2]);
 463       nds->InsertNextId(pts[3]);
 464       addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
 465       nds->Reset();
 466       nds->InsertNextId(pts[0]);
 467       nds->InsertNextId(pts[1]);
 468       nds->InsertNextId(pts[4]);
 469       addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
 470       nds->Reset();
 471       nds->InsertNextId(pts[1]);
 472       nds->InsertNextId(pts[2]);
 473       nds->InsertNextId(pts[4]);
 474       addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
 475       nds->Reset();
 476       nds->InsertNextId(pts[2]);
 477       nds->InsertNextId(pts[3]);
 478       nds->InsertNextId(pts[4]);
 479       addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
 480       nds->Reset();
 481       nds->InsertNextId(pts[3]);
 482       nds->InsertNextId(pts[0]);
 483       nds->InsertNextId(pts[4]);
 484       addToC2C(id_cell, _cells, c2c, 4, nds, cls, grid);
 485     } else if (grid->GetCellType(id_cell) == VTK_WEDGE) {
 486       c2c[i].resize(5);
 487       nds->Reset();
 488       nds->InsertNextId(pts[0]);
 489       nds->InsertNextId(pts[1]);
 490       nds->InsertNextId(pts[2]);
 491       addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
 492       nds->Reset();
 493       nds->InsertNextId(pts[3]);
 494       nds->InsertNextId(pts[4]);
 495       nds->InsertNextId(pts[5]);
 496       addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
 497       nds->Reset();
 498       nds->InsertNextId(pts[0]);
 499       nds->InsertNextId(pts[1]);
 500       nds->InsertNextId(pts[4]);
 501       nds->InsertNextId(pts[3]);
 502       addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
 503       nds->Reset();
 504       nds->InsertNextId(pts[1]);
 505       nds->InsertNextId(pts[4]);
 506       nds->InsertNextId(pts[5]);
 507       nds->InsertNextId(pts[2]);
 508       addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
 509       nds->Reset();
 510       nds->InsertNextId(pts[0]);
 511       nds->InsertNextId(pts[2]);
 512       nds->InsertNextId(pts[5]);
 513       nds->InsertNextId(pts[3]);
 514       addToC2C(id_cell, _cells, c2c, 4, nds, cls, grid);
 515     } else if (grid->GetCellType(id_cell) == VTK_HEXAHEDRON) {
 516       c2c[i].resize(6);
 517       nds->Reset();
 518       nds->InsertNextId(pts[0]);
 519       nds->InsertNextId(pts[3]);
 520       nds->InsertNextId(pts[2]);
 521       nds->InsertNextId(pts[1]);
 522       addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
 523       nds->Reset();
 524       nds->InsertNextId(pts[4]);
 525       nds->InsertNextId(pts[5]);
 526       nds->InsertNextId(pts[6]);
 527       nds->InsertNextId(pts[7]);
 528       addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
 529       nds->Reset();
 530       nds->InsertNextId(pts[0]);
 531       nds->InsertNextId(pts[1]);
 532       nds->InsertNextId(pts[5]);
 533       nds->InsertNextId(pts[4]);
 534       addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
 535       nds->Reset();
 536       nds->InsertNextId(pts[3]);
 537       nds->InsertNextId(pts[7]);
 538       nds->InsertNextId(pts[6]);
 539       nds->InsertNextId(pts[2]);
 540       addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
 541       nds->Reset();
 542       nds->InsertNextId(pts[0]);
 543       nds->InsertNextId(pts[4]);
 544       nds->InsertNextId(pts[7]);
 545       nds->InsertNextId(pts[3]);
 546       addToC2C(id_cell, _cells, c2c, 4, nds, cls, grid);
 547       nds->Reset();
 548       nds->InsertNextId(pts[1]);
 549       nds->InsertNextId(pts[2]);
 550       nds->InsertNextId(pts[6]);
 551       nds->InsertNextId(pts[5]);
 552       addToC2C(id_cell, _cells, c2c, 5, nds, cls, grid);
 553     };
 554   };
 555 };
 556
 557 void EgVtkObject::getNodesFromCells
 558 (
 559   QVector<vtkIdType>  &cells,
 560   QVector<vtkIdType>  &nodes,
 561   vtkUnstructuredGrid *grid
 562 )
 563 {
 564   QSet<vtkIdType> ex_nodes;
 565   vtkIdType id_cell;
 566   foreach(id_cell, cells) {
 567     vtkIdType *pts;
 568     vtkIdType  Npts;
 569     grid->GetCellPoints(id_cell, Npts, pts);
 570     for (int i = 0; i < Npts; ++i) {
 571       ex_nodes.insert(pts[i]);
 572     };
 573   };
 574   nodes.resize(ex_nodes.size());
 575   {
 576     int j = 0;
 577     vtkIdType i;
 578     foreach(i,ex_nodes) {
 579       nodes[j] = i;
 580       ++j;
 581     };
 582   };
 583 };
 584
 585 bool EgVtkObject::isVolume(vtkUnstructuredGrid *grid, vtkIdType id_cell)
 586 {
 587   bool isVol = false;
 588   if      (grid->GetCellType(id_cell) == VTK_TETRA)      isVol = true;
 589   else if (grid->GetCellType(id_cell) == VTK_PYRAMID)    isVol = true;
 590   else if (grid->GetCellType(id_cell) == VTK_WEDGE)      isVol = true;
 591   else if (grid->GetCellType(id_cell) == VTK_HEXAHEDRON) isVol = true;
 592   return isVol;
 593 };
 594
 595 bool EgVtkObject::isSurface(vtkUnstructuredGrid *grid, vtkIdType id_cell)
 596 {
 597   bool isSurf = false;
 598   if      (grid->GetCellType(id_cell) == VTK_TRIANGLE) isSurf = true;
 599   else if (grid->GetCellType(id_cell) == VTK_QUAD)     isSurf = true;
 600   return isSurf;
 601 };
 602
 603 void EgVtkObject::UpdateCellIndex(vtkUnstructuredGrid *grid)
 604 {
 605   if (!grid->GetCellData()->GetArray("cell_index")) {
 606     EG_VTKSP(vtkLongArray_t, cell_index);
 607     cell_index->SetName("cell_index");
 608     cell_index->SetNumberOfValues(grid->GetNumberOfCells());
 609     grid->GetCellData()->AddArray(cell_index);
 610   };
 611   EG_VTKDCC(vtkLongArray_t, cell_index, grid, "cell_index");
 612   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 613     cell_index->SetValue(id_cell, id_cell);
 614   };
 615 };
 616
 617 void EgVtkObject::UpdateNodeIndex(vtkUnstructuredGrid *grid)
 618 {
 619   if (!grid->GetPointData()->GetArray("node_index")) {
 620     EG_VTKSP(vtkLongArray_t, node_index);
 621     node_index->SetName("node_index");
 622     node_index->SetNumberOfValues(grid->GetNumberOfPoints());
 623     grid->GetPointData()->AddArray(node_index);
 624   };
 625   EG_VTKDCN(vtkLongArray_t, node_index, grid, "node_index");
 626   for (vtkIdType pointId = 0; pointId < grid->GetNumberOfPoints(); ++pointId) {
 627     node_index->SetValue(pointId, pointId);
 628   };
 629 };
 630
 631 void EgVtkObject::addToPolyData
 632 (
 633   QVector<vtkIdType>  &cells,
 634   vtkPolyData         *pdata,
 635   vtkUnstructuredGrid *grid
 636 )
 637 {
 638   UpdateCellIndex(grid);
 639   UpdateNodeIndex(grid);
 640   QVector<vtkIdType> nodes;
 641   QVector<int>       _nodes;
 642   getNodesFromCells(cells, nodes, grid);
 643   createNodeMapping(nodes, _nodes, grid);
 644   EG_VTKSP(vtkDoubleArray, pcoords);
 645   pcoords->SetNumberOfComponents(3);
 646   pcoords->SetNumberOfTuples(nodes.size());
 647   EG_VTKSP(vtkPoints, points);
 648   points->SetData(pcoords);
 649   pdata->SetPoints(points);
 650   pdata->Allocate(cells.size());
 651   if (!pdata->GetCellData()->GetArray("cell_index")) {
 652     EG_VTKSP(vtkLongArray_t, cell_index);
 653     cell_index->SetName("cell_index");
 654     //cell_index->SetNumberOfValues(cells.size());
 655     pdata->GetCellData()->AddArray(cell_index);
 656   };
 657   if (!pdata->GetPointData()->GetArray("node_index")) {
 658     EG_VTKSP(vtkLongArray_t, node_index);
 659     node_index->SetName("node_index");
 660     //node_index->SetNumberOfValues(nodes.size());
 661     pdata->GetPointData()->AddArray(node_index);
 662   };
 663   EG_VTKDCC(vtkLongArray_t, pd_cell_index, pdata, "cell_index");
 664   EG_VTKDCN(vtkLongArray_t, pd_node_index, pdata, "node_index");
 665   pd_cell_index->SetNumberOfValues(cells.size());
 666   pd_node_index->SetNumberOfValues(nodes.size());
 667   for (int i_cell = 0; i_cell < cells.size(); ++i_cell) {
 668     vtkIdType id_cell = cells[i_cell];
 669     vtkIdType cellType = grid->GetCellType(id_cell);
 670     if ((cellType != VTK_TRIANGLE) && (cellType != VTK_QUAD)) {
 671       EG_ERR_RETURN("unsupported cell type for this operation");
 672     };
 673     vtkIdType Npts, *pts;
 674     grid->GetCellPoints(id_cell, Npts, pts);
 675     vtkIdType *new_pts = new vtkIdType[Npts];
 676     for (int i = 0; i < Npts; ++i) {
 677       new_pts[i] = _nodes[pts[i]];
 678     };
 679     vtkIdType newCellId = pdata->InsertNextCell(cellType, Npts, new_pts);
 680     pd_cell_index->SetValue(newCellId, id_cell);
 681     delete [] new_pts;
 682   };
 683   for (int i_node = 0; i_node < nodes.size(); ++i_node) {
 684     vec3_t x;
 685     grid->GetPoints()->GetPoint(nodes[i_node], x.data());
 686     pdata->GetPoints()->SetPoint(i_node, x.data());
 687     pd_node_index->SetValue(i_node, nodes[i_node]);
 688   };
 689 };
 690
 691 #define EGVTKOBJECT_COPYCELLDATA(FIELD,TYPE) \
 692 { \
 693 if (old_grid->GetCellData()->GetArray(FIELD)) { \
 694 EG_VTKDCC(TYPE, var1, old_grid, FIELD); \
 695 EG_VTKDCC(TYPE, var2, new_grid, FIELD); \
 696 var2->SetValue(newId, var1->GetValue(oldId)); \
 697 }; \
 698 };
 699
 700 void EgVtkObject::copyCellData
 701 (
 702   vtkUnstructuredGrid *old_grid,
 703   vtkIdType            oldId,
 704   vtkUnstructuredGrid *new_grid,
 705   vtkIdType            newId
 706 )
 707 {
 708   EGVTKOBJECT_COPYCELLDATA("vtk_type",      vtkIntArray);
 709   EGVTKOBJECT_COPYCELLDATA("cell_code",     vtkIntArray);
 710   EGVTKOBJECT_COPYCELLDATA("cell_index",    vtkLongArray_t);
 711   EGVTKOBJECT_COPYCELLDATA("cell_err_tet",  vtkDoubleArray);
 712   EGVTKOBJECT_COPYCELLDATA("cell_err_pria", vtkDoubleArray);
 713   EGVTKOBJECT_COPYCELLDATA("cell_err_prib", vtkDoubleArray);
 714   EGVTKOBJECT_COPYCELLDATA("cell_err_pric", vtkDoubleArray);
 715   EGVTKOBJECT_COPYCELLDATA("cell_err_prid", vtkDoubleArray);
 716   EGVTKOBJECT_COPYCELLDATA("cell_err_prie", vtkDoubleArray);
 717   EGVTKOBJECT_COPYCELLDATA("cell_err_prif", vtkDoubleArray);
 718 };
 719
 720 #define EGVTKOBJECT_COPYNODEDATA(FIELD,TYPE) \
 721 { \
 722 if (old_grid->GetPointData()->GetArray(FIELD)) { \
 723 EG_VTKDCN(TYPE, var1, old_grid, FIELD); \
 724 EG_VTKDCN(TYPE, var2, new_grid, FIELD); \
 725 var2->SetValue(newId, var1->GetValue(oldId)); \
 726 }; \
 727 };
 728
 729 void EgVtkObject::copyNodeData
 730 (
 731   vtkUnstructuredGrid *old_grid,
 732   vtkIdType            oldId,
 733   vtkUnstructuredGrid *new_grid,
 734   vtkIdType            newId
 735 )
 736 {
 737   EGVTKOBJECT_COPYNODEDATA("node_status", vtkIntArray);
 738   EGVTKOBJECT_COPYNODEDATA("node_layer",  vtkIntArray);
 739   EGVTKOBJECT_COPYNODEDATA("node_index",  vtkLongArray_t);
 740   EGVTKOBJECT_COPYNODEDATA("node_meshdensity",  vtkDoubleArray);
 741   EGVTKOBJECT_COPYNODEDATA("node_meshdensity_current",  vtkDoubleArray);
 742   EGVTKOBJECT_COPYNODEDATA("node_type",  vtkCharArray);
 743 };
 744
 745 #define EGVTKOBJECT_CREATECELLFIELD(FIELD,TYPE,OW) \
 746 if (!grid->GetCellData()->GetArray(FIELD)) { \
 747 EG_VTKSP(TYPE, var); \
 748 var->SetName(FIELD); \
 749 var->SetNumberOfValues(Ncells); \
 750 grid->GetCellData()->AddArray(var); \
 751 } else if (OW) { \
 752 EG_VTKDCC(TYPE, var, grid, FIELD); \
 753 var->SetNumberOfValues(Ncells); \
 754 };
 755
 756 #define EGVTKOBJECT_CREATENODEFIELD(FIELD,TYPE,OW) \
 757 if (!grid->GetPointData()->GetArray(FIELD)) { \
 758 EG_VTKSP(TYPE, var); \
 759 var->SetName(FIELD); \
 760 var->SetNumberOfValues(Nnodes); \
 761 grid->GetPointData()->AddArray(var); \
 762 for (int i = 0; i < grid->GetNumberOfPoints(); ++i) var->SetValue(i,0); \
 763 } else if (OW) { \
 764 EG_VTKDCN(TYPE, var, grid, FIELD); \
 765 var->SetNumberOfValues(Nnodes); \
 766 for (int i = 0; i < grid->GetNumberOfPoints(); ++i) var->SetValue(i,0); \
 767 };
 768
 769 void EgVtkObject::createBasicFields
 770 (
 771   vtkUnstructuredGrid *grid,
 772   vtkIdType            Ncells,
 773   vtkIdType            Nnodes,
 774   bool                 overwrite
 775 )
 776 {
 777   createBasicNodeFields(grid, Nnodes, overwrite);
 778   createBasicCellFields(grid, Ncells, overwrite);
 779 };
 780
 781 void EgVtkObject::createBasicCellFields
 782 (
 783   vtkUnstructuredGrid *grid,
 784   vtkIdType            Ncells,
 785   bool                 overwrite
 786 )
 787 {
 788   EGVTKOBJECT_CREATECELLFIELD("vtk_type" ,     vtkIntArray, overwrite);
 789   EGVTKOBJECT_CREATECELLFIELD("cell_code",     vtkIntArray, overwrite);
 790   EGVTKOBJECT_CREATECELLFIELD("cell_index",    vtkLongArray_t, overwrite);
 791   EGVTKOBJECT_CREATECELLFIELD("cell_err_tet",  vtkDoubleArray, overwrite);
 792   EGVTKOBJECT_CREATECELLFIELD("cell_err_pria", vtkDoubleArray, overwrite);
 793   EGVTKOBJECT_CREATECELLFIELD("cell_err_prib", vtkDoubleArray, overwrite);
 794   EGVTKOBJECT_CREATECELLFIELD("cell_err_pric", vtkDoubleArray, overwrite);
 795   EGVTKOBJECT_CREATECELLFIELD("cell_err_prid", vtkDoubleArray, overwrite);
 796   EGVTKOBJECT_CREATECELLFIELD("cell_err_prie", vtkDoubleArray, overwrite);
 797   EGVTKOBJECT_CREATECELLFIELD("cell_err_prif", vtkDoubleArray, overwrite);
 798   EGVTKOBJECT_CREATECELLFIELD("cell_VA",       vtkDoubleArray, overwrite);
 799 };
 800
 801 void EgVtkObject::createBasicNodeFields
 802 (
 803   vtkUnstructuredGrid *grid,
 804   vtkIdType            Nnodes,
 805   bool                 overwrite
 806 )
 807 {
 808   EGVTKOBJECT_CREATENODEFIELD("node_status", vtkIntArray,    overwrite);
 809   EGVTKOBJECT_CREATENODEFIELD("node_layer",  vtkIntArray,    overwrite);
 810   EGVTKOBJECT_CREATENODEFIELD("node_index",  vtkLongArray_t, overwrite);
 811   EGVTKOBJECT_CREATENODEFIELD("node_meshdensity",  vtkDoubleArray, overwrite);
 812   EGVTKOBJECT_CREATENODEFIELD("node_meshdensity_current",  vtkDoubleArray, overwrite);
 813   EGVTKOBJECT_CREATENODEFIELD("node_type",  vtkCharArray, overwrite);
 814 };
 815
 816 void EgVtkObject::allocateGrid
 817 (
 818   vtkUnstructuredGrid *grid,
 819   vtkIdType            Ncells,
 820   vtkIdType            Nnodes
 821 )
 822 {
 823   EG_VTKSP(vtkPoints,points);
 824   points->SetNumberOfPoints(Nnodes);
 825   grid->SetPoints(points);
 826   grid->Allocate(Ncells,max(vtkIdType(1),Ncells/10));
 827   createBasicFields(grid, Ncells, Nnodes);
 828 };
 829
 830 vec3_t EgVtkObject::cellCentre(vtkUnstructuredGrid *grid, vtkIdType id_cell)
 831 {
 832   vec3_t x,xc(0,0,0);
 833   vtkIdType *pts, N_pts;
 834   grid->GetCellPoints(id_cell, N_pts, pts);
 835   double f = 1.0/N_pts;
 836   for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 837     grid->GetPoint(pts[i_pts], x.data());
 838     xc += f*x;
 839   };
 840   return xc;
 841 };
 842
 843 void EgVtkObject::getRestCells(vtkUnstructuredGrid      *grid,
 844                                const QVector<vtkIdType> &cells,
 845                                QVector<vtkIdType>       &rest_cells)
 846 {
 847   QVector<bool> is_in_cells(grid->GetNumberOfCells(), false);
 848   foreach (vtkIdType id_cell, cells) {
 849     is_in_cells[id_cell] = true;
 850   };
 851   rest_cells.resize(grid->GetNumberOfCells() - cells.size());
 852   int i_rest_cells = 0;
 853   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 854     if (!is_in_cells[id_cell]) {
 855       rest_cells[i_rest_cells] = id_cell;
 856       ++i_rest_cells;
 857     };
 858   };
 859 };
 860
 861 void EgVtkObject::makeCopy(vtkUnstructuredGrid *src, vtkUnstructuredGrid *dst)
 862 {
 863   allocateGrid(dst, src->GetNumberOfCells(), src->GetNumberOfPoints());
 864   for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
 865     vec3_t x;
 866     src->GetPoints()->GetPoint(id_node, x.data());
 867     dst->GetPoints()->SetPoint(id_node, x.data());
 868     copyNodeData(src, id_node, dst, id_node);
 869   };
 870   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 871     vtkIdType N_pts, *pts;
 872     vtkIdType type_cell = src->GetCellType(id_cell);
 873     src->GetCellPoints(id_cell, N_pts, pts);
 874     vtkIdType id_new_cell = dst->InsertNextCell(type_cell, N_pts, pts);
 875     copyCellData(src, id_cell, dst, id_new_cell);
 876   };
 877 };
 878
 879 void EgVtkObject::makeCopyNoAlloc(vtkUnstructuredGrid *src, vtkUnstructuredGrid *dst)
 880 {
 881   for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
 882     vec3_t x;
 883     src->GetPoints()->GetPoint(id_node, x.data());
 884     dst->GetPoints()->SetPoint(id_node, x.data());
 885     copyNodeData(src, id_node, dst, id_node);
 886   };
 887   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 888     vtkIdType N_pts, *pts;
 889     vtkIdType type_cell = src->GetCellType(id_cell);
 890     src->GetCellPoints(id_cell, N_pts, pts);
 891     vtkIdType id_new_cell = dst->InsertNextCell(type_cell, N_pts, pts);
 892     copyCellData(src, id_cell, dst, id_new_cell);
 893   };
 894 };
 895
 896 // void EgVtkObject::makeCopyNoAllocFiltered(vtkUnstructuredGrid *src, vtkUnstructuredGrid *dst, vector <bool> DeadNode, QVector <QSet <vtkIdType>> newCells)
 897 void EgVtkObject::makeCopyNoAllocFiltered(vtkUnstructuredGrid *src, vtkUnstructuredGrid *dst, vector <bool> DeadNode)
 898 {
 899   vtkIdType src_N_points=src->GetNumberOfPoints();
 900   vector <vtkIdType> OffSet(src_N_points);
 901   vtkIdType dst_id_node=0;
 902   for (vtkIdType src_id_node = 0; src_id_node < src_N_points; ++src_id_node) {
 903     if(!DeadNode[src_id_node])
 904     {
 905       vec3_t x;
 906       src->GetPoints()->GetPoint(src_id_node, x.data());
 907       dst->GetPoints()->SetPoint(dst_id_node, x.data());
 908       copyNodeData(src, src_id_node, dst, dst_id_node);
 909       OffSet[src_id_node]=src_id_node-dst_id_node;
 910       dst_id_node++;
 911     }
 912     else
 913     {
 914       //search closest node
 915     }
 916   };
 917   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 918     vtkIdType N_pts, *src_pts, *dst_pts;
 919     vtkIdType type_cell = src->GetCellType(id_cell);
 920     src->GetCellPoints(id_cell, N_pts, src_pts);
 921     src->GetCellPoints(id_cell, N_pts, dst_pts);
 922     bool DeadCell=false;
 923     for(int i=0;i<N_pts;i++)
 924     {
 925       if(DeadNode[src_pts[i]]) {DeadCell=true;}
 926       dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
 927     }
 928     if(!DeadCell)
 929     {
 930       vtkIdType id_new_cell = dst->InsertNextCell(type_cell, N_pts, dst_pts);
 931       copyCellData(src, id_cell, dst, id_new_cell);
 932     }
 933   };
 934 }
 935
 936 int EgVtkObject::findVolumeCell
 937 (
 938   vtkUnstructuredGrid      *grid,
 939   vtkIdType                id_surf,
 940   const QVector<int>       _nodes,
 941   const QVector<vtkIdType> cells,
 942   const QVector<int>       _cells,
 943   QVector<QSet<int> >       &n2c
 944 )
 945 {
 946   vtkIdType N_pts, *pts;
 947   if (_cells.size()) N_pts = N_pts; // dummy statement to get rid of compiler warning ...
 948   grid->GetCellPoints(id_surf, N_pts, pts);
 949   QVector<QSet<int> > inters(N_pts-1);
 950   setIntersection(n2c[_nodes[pts[0]]], n2c[_nodes[pts[1]]], inters[0]);
 951   int i_pts = 2;
 952   while (i_pts < N_pts) {
 953     setIntersection(inters[i_pts-2], n2c[_nodes[pts[i_pts]]], inters[i_pts-1]);
 954     ++i_pts;
 955   };
 956   if (inters[N_pts-2].size() == 0) {
 957     return -1;
 958   } else if (inters[N_pts-2].size() > 2) {
 959     EG_BUG;
 960   };
 961   vtkIdType id_vol = -1;
 962   foreach (int i_cells, inters[N_pts-2]) {
 963     if (cells[i_cells] != id_surf) {
 964       id_vol = cells[i_cells];
 965     };
 966   };
 967   return id_vol;
 968 };
 969
 970 void EgVtkObject::setBoundaryCodes(const QSet<int> &bcs)
 971 {
 972   boundary_codes.clear();
 973   int bc;
 974   foreach(bc, bcs) {
 975     boundary_codes.insert(bc);
 976   };
 977 };
 978
 979 void EgVtkObject::createIndices(vtkUnstructuredGrid *grid)
 980 {
 981   if (!grid->GetCellData()->GetArray("cell_index")) {
 982     EG_VTKSP(vtkLongArray_t, var);
 983     var->SetName("cell_index");
 984     var->SetNumberOfValues(grid->GetNumberOfCells());
 985     grid->GetCellData()->AddArray(var);
 986   } else {
 987     EG_VTKDCC(vtkLongArray_t, var, grid, "cell_index");
 988     var->SetNumberOfValues(grid->GetNumberOfCells());
 989   };
 990   EG_VTKDCC(vtkLongArray_t, cell_index, grid, "cell_index");
 991   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 992     cell_index->SetValue(id_cell, id_cell);
 993   };
 994
 995   if (!grid->GetCellData()->GetArray("vtk_type")) {
 996     EG_VTKSP(vtkIntArray, var);
 997     var->SetName("vtk_type");
 998     var->SetNumberOfValues(grid->GetNumberOfCells());
 999     grid->GetCellData()->AddArray(var);
1000   } else {
1001     EG_VTKDCC(vtkIntArray, var, grid, "vtk_type");
1002     var->SetNumberOfValues(grid->GetNumberOfCells());
1003   };
1004   EG_VTKDCC(vtkIntArray, vtk_type, grid, "vtk_type");
1005   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
1006     vtk_type->SetValue(id_cell, grid->GetCellType(id_cell));
1007   };
1008
1009   if (!grid->GetCellData()->GetArray("node_index")) {
1010     EG_VTKSP(vtkLongArray_t, var);
1011     var->SetName("node_index");
1012     var->SetNumberOfValues(grid->GetNumberOfPoints());
1013     grid->GetPointData()->AddArray(var);
1014   } else {
1015     EG_VTKDCC(vtkLongArray_t, var, grid, "node_index");
1016     var->SetNumberOfValues(grid->GetNumberOfPoints());
1017   };
1018   EG_VTKDCN(vtkLongArray_t, node_index, grid, "node_index");
1019   for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
1020     node_index->SetValue(id_node, id_node);
1021   };
1022 };
1023
1024 BoundaryCondition EgVtkObject::getBC(int bc)
1025 {
1026   return GuiMainWindow::pointer()->getBC(bc);
1027 };
1028
1029 int EgVtkObject::getSet(QString group, QString key, int value, int& variable)
1030 {
1031   QSettings *qset = GuiMainWindow::settings();
1032   QString typed_key = "int/" + key;
1033   qset->beginGroup(group);
1034   //if key=value pair not found in settings file, write it
1035   if (!qset->contains(typed_key)) qset->setValue(typed_key,value);
1036   //read key value from settings file and assign it to variable
1037   variable = (qset->value(typed_key,variable)).toInt();
1038   qset->endGroup();
1039   return(variable);
1040 }
1041
1042 double EgVtkObject::getSet(QString group, QString key, double value, double& variable)
1043 {
1044   QSettings *qset = GuiMainWindow::settings();
1045   QString typed_key = "double/" + key;
1046   qset->beginGroup(group);
1047   //if key=value pair not found in settings file, write it
1048   if (!qset->contains(typed_key)) qset->setValue(typed_key,value);
1049   //read key value from settings file and assign it to variable
1050   variable = (qset->value(typed_key,variable)).toDouble();
1051   qset->endGroup();
1052   return(variable);
1053 }
1054
1055 bool EgVtkObject::getSet(QString group, QString key, bool value, bool& variable)
1056 {
1057   QSettings *qset = GuiMainWindow::settings();
1058   QString typed_key = "bool/" + key;
1059   qset->beginGroup(group);
1060   Qt::CheckState state = (Qt::CheckState) ( value ? 2 : 0 );
1061   //if key=value pair not found in settings file, write it
1062   if (!qset->contains(typed_key)) qset->setValue(typed_key,state);
1063   //read key value from settings file and assign it to variable
1064   variable = (qset->value(typed_key,variable)).toBool();
1065   qset->endGroup();
1066   return(variable);
1067 }
1068
1069 ///////////////////////////////////////////
1070 int cout_grid(ostream &stream, vtkUnstructuredGrid *grid, bool npoints, bool ncells, bool points, bool cells)
1071 {
1072   stream<<"============="<<endl;
1073   if(npoints) stream << "grid->GetNumberOfPoints()=" << grid->GetNumberOfPoints() << endl;
1074   if(ncells) stream << "grid->GetNumberOfCells()=" << grid->GetNumberOfCells() << endl;
1075   if(points) {
1076     for (vtkIdType i = 0; i < grid->GetNumberOfPoints(); ++i) {
1077       vec3_t x;
1078       grid->GetPoint(i, x.data());
1079       stream << "Vertex " << i << " = " << x << endl;
1080     };
1081   }
1082   if(cells) {
1083     EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
1084     for (vtkIdType i = 0; i < grid->GetNumberOfCells(); ++i) {
1085       vtkCell *C = (vtkCell *) vtkCell::New();
1086       C=grid->GetCell(i);
1087       vtkIdType npts=C->GetNumberOfPoints();
1088       vtkIdType* pts;
1089       grid->GetCellPoints(i, npts, pts);
1090       stream << "Cell " << i << " = ";
1091       for(int j=0;j<npts;j++) stream << pts[j] << " ";
1092       stream << "boundary_code=" << cell_code->GetValue(i);
1093       stream << endl;
1094     };
1095   }
1096   stream<<"============="<<endl;
1097   return 0;
1098 }
1099
1100 ///////////////////////////////////////////
1101 int addPoint(vtkUnstructuredGrid* a_grid,vtkIdType index,vec3_t x)
1102 {
1103   a_grid->GetPoints()->SetPoint(index, x.data());
1104   return(0);
1105 }
1106
1107 int addCell(vtkUnstructuredGrid* a_grid, vtkIdType A, vtkIdType B, vtkIdType C, int bc)
1108 {
1109   vtkIdType npts=3;
1110   vtkIdType pts[3];
1111   pts[0]=A;
1112   pts[1]=B;
1113   pts[2]=C;
1114   vtkIdType newCellId = a_grid->InsertNextCell(VTK_TRIANGLE,npts,pts);
1115   EG_VTKDCC(vtkIntArray, cell_code, a_grid, "cell_code");
1116   cell_code->SetValue(newCellId, bc);
1117   return(0);
1118 }
1119
1120
1121 ///////////////////////////////////////////
1122
1123 int getShortestSide(vtkIdType a_id_cell,vtkUnstructuredGrid* a_grid)
1124 {
1125   vtkIdType N_pts, *pts;
1126   a_grid->GetCellPoints(a_id_cell, N_pts, pts);
1127   vec3_t* x=new vec3_t[N_pts];
1128   for(int i=0;i<N_pts;i++) a_grid->GetPoints()->GetPoint(pts[i], x[i].data());
1129   int id_minlen=0;
1130   double minlen=(x[1]-x[0]).abs();
1131   for(int i=1;i<N_pts;i++)
1132   {
1133     double len=(x[(i+1)%N_pts]-x[i]).abs();
1134     if(len<minlen){
1135       minlen=len;
1136       id_minlen=i;
1137     }
1138   }
1139   delete x;
1140   return(id_minlen);
1141 }
1142
1143 int getLongestSide(vtkIdType a_id_cell,vtkUnstructuredGrid* a_grid)
1144 {
1145   vtkIdType N_pts, *pts;
1146   a_grid->GetCellPoints(a_id_cell, N_pts, pts);
1147   vec3_t* x=new vec3_t[N_pts];
1148   for(int i=0;i<N_pts;i++) a_grid->GetPoints()->GetPoint(pts[i], x[i].data());
1149   int id_maxlen=0;
1150   double maxlen=(x[1]-x[0]).abs();
1151   cout<<"maxlen="<<maxlen<<endl;
1152   for(int i=1;i<N_pts;i++)
1153   {
1154     double len=(x[(i+1)%N_pts]-x[i]).abs();
1155     cout<<"len["<<i<<"]="<<len<<endl;
1156     if(len>maxlen){
1157       maxlen=len;
1158       id_maxlen=i;
1159     }
1160   }
1161   delete x;
1162   return(id_maxlen);
1163 }
1164
1165 //get number of the edge corresponding to node1-node2
1166 int getSide(vtkIdType a_id_cell,vtkUnstructuredGrid* a_grid,vtkIdType a_id_node1,vtkIdType a_id_node2)
1167 {
1168   vtkIdType N_pts, *pts;
1169   a_grid->GetCellPoints(a_id_cell, N_pts, pts);
1170   QVector <vtkIdType> edge(2);
1171
1172   int n=0;
1173   for(int i=0;i<N_pts;i++)
1174   {
1175     if(pts[i]==a_id_node1) { edge[0]=i;n++;}
1176     if(pts[i]==a_id_node2) { edge[1]=i;n++;}
1177   }
1178   if(n!=2){
1179     EG_BUG;
1180     return(-1);
1181   }
1182   qSort(edge.begin(),edge.end());
1183   if(edge[0]==0 && edge[1]==N_pts-1) return(N_pts-1);
1184   else return(edge[0]);
1185 }
1186 ///////////////////////////////////////////
1187
1188 QSet <int> complementary_bcs(QSet <int> &bcs, vtkUnstructuredGrid *a_grid, QVector <vtkIdType> &a_cells)
1189 {
1190   QSet <int> bcs_complement;
1191   EG_VTKDCC(vtkIntArray, bc, a_grid, "cell_code");
1192   foreach (vtkIdType id_cell, a_cells) {
1193     int code=bc->GetValue(id_cell);
1194     if (!bcs.contains(code)) {
1195       bcs_complement.insert(code);
1196     };
1197   }
1198   return(bcs_complement);
1199 }
1200
1201 ///////////////////////////////////////////
1202
1203 QString cell2str(vtkIdType id_cell,vtkUnstructuredGrid* grid)
1204 {
1205   QString tmp;
1206   tmp.setNum(id_cell);
1207   QString txt = "id_cell=" + tmp;
1208
1209   vtkIdType N_pts, *pts;
1210   grid->GetCellPoints(id_cell, N_pts, pts);
1211
1212   txt += " [";
1213   for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
1214     tmp.setNum(pts[i_pts]);
1215     txt += tmp;
1216     if (i_pts < N_pts-1) {
1217       txt += ",";
1218     };
1219   };
1220   txt += "]";
1221   return(txt);
1222 }
1223
1224
1225
1226 ///////////////////////////////////////////
1227
1228 Qt::CheckState int2CheckState(int a)
1229 {
1230   if(a==0) return(Qt::Unchecked);
1231   if(a==1) return(Qt::PartiallyChecked);
1232   else return(Qt::Checked);
1233 }
1234
1235 int CheckState2int(Qt::CheckState a)
1236 {
1237   if(a==Qt::Unchecked) return(0);
1238   if(a==Qt::PartiallyChecked) return(1);
1239   else return(Qt::Checked);
1240 }
1241
1242 // ///////////////////////////////////////////
1243 // /* Here is how we we get QTextStreams that look like iostreams */
1244 // Qcin=QTextStream(stdin, QIODevice::ReadOnly);
1245 // Qcout=QTextStream(stdout, QIODevice::WriteOnly);
1246 // Qcerr=QTextStream(stderr, QIODevice::WriteOnly);
1247 // ///////////////////////////////////////////
1248
1249 pair<vtkIdType,vtkIdType> OrderedPair(vtkIdType a, vtkIdType b)
1250 {
1251   vtkIdType x=min(a,b);
1252   vtkIdType y=max(a,b);
1253   return(pair<vtkIdType,vtkIdType>(x,y));
1254 }
1255
1256 vtkIdType nextcell(vtkIdType a_cell, vtkIdType a_node, QVector< QVector< int > > &a_c2c, vtkUnstructuredGrid *a_grid)
1257 {
1258   vtkIdType N_pts, *pts;
1259   a_grid->GetCellPoints(a_cell, N_pts, pts);
1260
1261   int i;
1262   for(i=0;i<N_pts;i++)
1263   {
1264     if(pts[i]==a_node) break;
1265   }
1266   return a_c2c[a_cell][i];
1267 }
1268
1269 const char* VertexType2Str(char T)
1270 {
1271   if(T==VTK_SIMPLE_VERTEX) return("VTK_SIMPLE_VERTEX");
1272   if(T==VTK_FIXED_VERTEX) return("VTK_FIXED_VERTEX");
1273   if(T==VTK_FEATURE_EDGE_VERTEX) return("VTK_FEATURE_EDGE_VERTEX");
1274   if(T==VTK_BOUNDARY_EDGE_VERTEX) return("VTK_BOUNDARY_EDGE_VERTEX");
1275   else return("Unknown vertex type");
1276 }
1277
1278 char Str2VertexType(QString S)
1279 {
1280   if(S=="VTK_SIMPLE_VERTEX") return((char)0);
1281   if(S=="VTK_FIXED_VERTEX") return((char)1);
1282   if(S=="VTK_FEATURE_EDGE_VERTEX") return((char)2);
1283   if(S=="VTK_BOUNDARY_EDGE_VERTEX") return((char)3);
1284   else return((char)-1);
1285 }
1286
1287 const char* vertex_type(char T)
1288 {
1289   if(T==VTK_SIMPLE_VERTEX) return("VTK_SIMPLE_VERTEX");
1290   if(T==VTK_FIXED_VERTEX) return("VTK_FIXED_VERTEX");
1291   if(T==VTK_FEATURE_EDGE_VERTEX) return("VTK_FEATURE_EDGE_VERTEX");
1292   if(T==VTK_BOUNDARY_EDGE_VERTEX) return("VTK_BOUNDARY_EDGE_VERTEX");
1293   else return("Unknown vertex type");
1294 }