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