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
  32 void EgVtkObject::computeNormals
  33 (
  34   QVector<vec3_t>     &cell_normals,
  35   QVector<vec3_t>     &node_normals,
  36   QVector<vtkIdType>  &cells,
  37   QVector<vtkIdType>  &nodes,
  38   vtkUnstructuredGrid *grid
  39 )
  40 {
  41   using namespace GeometryTools;
  42
  43   cell_normals.resize(cells.size());
  44   node_normals.fill(vec3_t(0,0,0), nodes.size());
  45   QVector<int> g2s;
  46   createNodeMapping(nodes, g2s, grid);
  47   for (int i_cell = 0; i_cell < cells.count(); ++i_cell) {
  48     vtkIdType id_cell = cells[i_cell];
  49     vtkIdType *pts;
  50     vtkIdType npts;
  51     grid->GetCellPoints(id_cell, npts, pts);
  52     cell_normals[i_cell] = cellNormal(grid, id_cell);
  53     cell_normals[i_cell].normalise();
  54     for (int i_pts = 0; i_pts < npts; ++i_pts) {
  55       if (g2s[pts[i_pts]] != -1) {
  56         node_normals[g2s[pts[i_pts]]] += cell_normals[i_cell];
  57       };
  58     };
  59   };
  60   for (int i_node = 0; i_node < nodes.count(); ++i_node) {
  61     node_normals[i_node].normalise();
  62     //cout << node_normals[i_node] << endl;
  63   };
  64 };
  65
  66 void EgVtkObject::createNodeMapping
  67 (
  68   QVector<vtkIdType>  &nodes,
  69   QVector<int>        &_nodes,
  70   vtkUnstructuredGrid *grid
  71 )
  72 {
  73   _nodes.fill(-1,grid->GetNumberOfPoints());
  74   vtkIdType id;
  75   int i_sub = 0;
  76   foreach(id,nodes) {
  77     _nodes[id] = i_sub;
  78     ++i_sub;
  79   };
  80 };
  81
  82 void EgVtkObject::createCellMapping
  83 (
  84   QVector<vtkIdType>  &cells,
  85   QVector<int>        &_cells,
  86   vtkUnstructuredGrid *grid
  87 )
  88 {
  89   _cells.fill(-1,grid->GetNumberOfCells());
  90   vtkIdType id;
  91   int i_sub = 0;
  92   foreach(id,cells) {
  93     _cells[id] = i_sub;
  94     ++i_sub;
  95   };
  96 };
  97
  98 void EgVtkObject::createNodeToBcMapping
  99 (
 100   QVector<QSet<int> > &bcs,
 101   vtkUnstructuredGrid *grid
 102 )
 103 {
 104   bcs.fill(QSet<int>(), grid->GetNumberOfPoints());
 105   grid->BuildLinks();
 106   EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 107   for (vtkIdType nodeId = 0; nodeId < grid->GetNumberOfPoints(); ++nodeId) {
 108     int Ncells = grid->GetCellLinks()->GetNcells(nodeId);
 109     for (int i = 0; i < Ncells; ++i) {
 110       vtkIdType id_cell = grid->GetCellLinks()->GetCells(nodeId)[i];
 111       vtkIdType ct = grid->GetCellType(id_cell);
 112       if ((ct == VTK_TRIANGLE) || (ct = VTK_QUAD)) {
 113         if (cell_code->GetValue(id_cell) > 0) {
 114           bcs[nodeId].insert(cell_code->GetValue(id_cell));
 115         };
 116       };
 117     };
 118   };
 119 };
 120
 121 void EgVtkObject::createNodeToCell
 122 (
 123   QVector<vtkIdType>  &cells,
 124   QVector<vtkIdType>  &nodes,
 125   QVector<int>        &_nodes,
 126   QVector<QSet<int> > &n2c,
 127   vtkUnstructuredGrid *grid
 128 )
 129 {
 130   n2c.fill(QSet<int>(), nodes.size());
 131   for (vtkIdType i_cells = 0; i_cells < cells.size(); ++i_cells) {
 132     vtkIdType *pts;
 133     vtkIdType  Npts;
 134     grid->GetCellPoints(cells[i_cells], Npts, pts);
 135     for (int i_pts = 0; i_pts < Npts; ++i_pts) {
 136       n2c[_nodes[pts[i_pts]]].insert(i_cells);
 137     };
 138   };
 139 };
 140
 141 void EgVtkObject::addToN2N
 142 (
 143   QVector<QSet<int> > &n2n,
 144   int                  n1,
 145   int                  n2
 146 )
 147 {
 148   n2n[n1].insert(n2);
 149   n2n[n2].insert(n1);
 150 };
 151
 152 void EgVtkObject::createNodeToNode
 153 (
 154   QVector<vtkIdType>  &cells,
 155   QVector<vtkIdType>  &nodes,
 156   QVector<int>        &_nodes,
 157   QVector<QSet<int> > &n2n,
 158   vtkUnstructuredGrid *grid
 159 )
 160 {
 161   n2n.fill(QSet<int>(), nodes.size());
 162   foreach (vtkIdType id_cell, cells) {
 163     vtkIdType *pts;
 164     vtkIdType  Npts;
 165     grid->GetCellPoints(id_cell, Npts, pts);
 166     vector<int> n(Npts);
 167     for (int i = 0; i < Npts; ++i) {
 168       n[i] = _nodes[pts[i]];
 169     };
 170     vtkIdType cellType = grid->GetCellType(id_cell);
 171     if (cellType == VTK_TRIANGLE) {
 172       addToN2N(n2n, n[0], n[1]);
 173       addToN2N(n2n, n[1], n[2]);
 174       addToN2N(n2n, n[2], n[0]);
 175     } else if (cellType == VTK_QUAD) {
 176       addToN2N(n2n, n[0], n[1]);
 177       addToN2N(n2n, n[1], n[2]);
 178       addToN2N(n2n, n[2], n[3]);
 179       addToN2N(n2n, n[3], n[0]);
 180     } else if (cellType == VTK_TETRA) {
 181       addToN2N(n2n, n[0], n[1]);
 182       addToN2N(n2n, n[0], n[2]);
 183       addToN2N(n2n, n[0], n[3]);
 184       addToN2N(n2n, n[1], n[2]);
 185       addToN2N(n2n, n[1], n[3]);
 186       addToN2N(n2n, n[2], n[3]);
 187     } else if (cellType == VTK_PYRAMID) {
 188       addToN2N(n2n, n[0], n[1]);
 189       addToN2N(n2n, n[0], n[3]);
 190       addToN2N(n2n, n[0], n[4]);
 191       addToN2N(n2n, n[1], n[2]);
 192       addToN2N(n2n, n[1], n[4]);
 193       addToN2N(n2n, n[2], n[3]);
 194       addToN2N(n2n, n[2], n[4]);
 195       addToN2N(n2n, n[3], n[4]);
 196     } else if (cellType == VTK_WEDGE) {
 197       addToN2N(n2n, n[0], n[1]);
 198       addToN2N(n2n, n[0], n[2]);
 199       addToN2N(n2n, n[0], n[3]);
 200       addToN2N(n2n, n[1], n[2]);
 201       addToN2N(n2n, n[1], n[4]);
 202       addToN2N(n2n, n[2], n[5]);
 203       addToN2N(n2n, n[3], n[4]);
 204       addToN2N(n2n, n[3], n[5]);
 205       addToN2N(n2n, n[4], n[5]);
 206     } else if (cellType == VTK_HEXAHEDRON) {
 207       addToN2N(n2n, n[0], n[1]);
 208       addToN2N(n2n, n[0], n[3]);
 209       addToN2N(n2n, n[0], n[4]);
 210       addToN2N(n2n, n[1], n[2]);
 211       addToN2N(n2n, n[1], n[5]);
 212       addToN2N(n2n, n[2], n[3]);
 213       addToN2N(n2n, n[2], n[6]);
 214       addToN2N(n2n, n[3], n[7]);
 215       addToN2N(n2n, n[4], n[5]);
 216       addToN2N(n2n, n[4], n[7]);
 217       addToN2N(n2n, n[5], n[6]);
 218       addToN2N(n2n, n[6], n[7]);
 219     };
 220   };
 221 };
 222
 223 void EgVtkObject::getAllCells
 224 (
 225   QVector<vtkIdType>  &cells,
 226   vtkUnstructuredGrid *grid
 227 )
 228 {
 229   int N = 0;
 230   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 231     ++N;
 232   };
 233   cells.resize(N);
 234   N = 0;
 235   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 236     cells[N] = id_cell;
 237     ++N;
 238   };
 239 };
 240
 241 void EgVtkObject::getAllCellsOfType
 242 (
 243   vtkIdType            type,
 244   QVector<vtkIdType>  &cells,
 245   vtkUnstructuredGrid *grid
 246 )
 247 {
 248   int N = 0;
 249   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 250     if (grid->GetCellType(id_cell) == type) {
 251       ++N;
 252     };
 253   };
 254   cells.resize(N);
 255   N = 0;
 256   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 257     if (grid->GetCellType(id_cell) == type) {
 258       cells[N] = id_cell;
 259       ++N;
 260     };
 261   };
 262 };
 263
 264
 265 void EgVtkObject::getAllVolumeCells
 266 (
 267   QVector<vtkIdType>  &cells,
 268   vtkUnstructuredGrid *grid
 269 )
 270 {
 271   int N = 0;
 272   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 273     if (isVolume(id_cell, grid)) {
 274       ++N;
 275     };
 276   };
 277   cells.resize(N);
 278   N = 0;
 279   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 280     if (isVolume(id_cell, grid)) {
 281       cells[N] = id_cell;
 282       ++N;
 283     };
 284   };
 285 };
 286
 287 void EgVtkObject::getAllSurfaceCells
 288 (
 289   QVector<vtkIdType>  &cells,
 290   vtkUnstructuredGrid *grid
 291 )
 292 {
 293   int N = 0;
 294   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 295     if (isSurface(id_cell, grid)) {
 296       ++N;
 297     };
 298   };
 299   cells.resize(N);
 300   N = 0;
 301   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 302     if (isSurface(id_cell, grid)) {
 303       cells[N] = id_cell;
 304       ++N;
 305     };
 306   };
 307 };
 308
 309 void EgVtkObject::getSurfaceCells
 310 (
 311   QSet<int>           &bcs,
 312   QVector<vtkIdType>  &cells,
 313   vtkUnstructuredGrid *grid
 314 )
 315 {
 316   int N = 0;
 317   EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
 318   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 319     if (isSurface(id_cell, grid)) {
 320       if (bcs.contains(cell_code->GetValue(id_cell))) {
 321         ++N;
 322       };
 323     };
 324   };
 325   cells.resize(N);
 326   N = 0;
 327   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 328     if (isSurface(id_cell, grid)) {
 329       if (bcs.contains(cell_code->GetValue(id_cell))) {
 330         cells[N] = id_cell;
 331         ++N;
 332       };
 333     };
 334   };
 335 };
 336
 337 void EgVtkObject::addToC2C
 338 (
 339   vtkIdType               id_cell,
 340   QVector<int>           &_cells,
 341   QVector<QVector<int> > &c2c,
 342   int                     j,
 343   vtkIdList              *nds,
 344   vtkIdList              *cls,
 345   vtkUnstructuredGrid    *grid
 346 )
 347 {
 348   c2c[_cells[id_cell]][j] = -1;
 349   grid->GetCellNeighbors(id_cell, nds, cls);
 350   for (int i = 0; i < cls->GetNumberOfIds(); ++i) {
 351     if (cls->GetId(i) != id_cell) {
 352       if (_cells[cls->GetId(i)] != -1) {
 353         c2c[_cells[id_cell]][j] = _cells[cls->GetId(i)];
 354       };
 355     };
 356   };
 357 };
 358
 359
 360 void EgVtkObject::createCellToCell
 361 (
 362   QVector<vtkIdType>     &cells,
 363   QVector<QVector<int> > &c2c,
 364   vtkUnstructuredGrid    *grid
 365 )
 366 {
 367   // GetCellNeighbors(vtkIdType id_cell, vtkIdList *ptIds, vtkIdList *id_cells)
 368   grid->BuildLinks();
 369   QVector<int> _cells;
 370   createCellMapping(cells, _cells, grid);
 371   c2c.fill(QVector<int>(), cells.size());
 372   EG_VTKSP(vtkIdList, nds);
 373   EG_VTKSP(vtkIdList, cls);
 374   for (int i = 0; i < cells.size(); ++i) {
 375     vtkIdType id_cell = cells[i];
 376     vtkIdType *pts;
 377     vtkIdType  Npts;
 378     grid->GetCellPoints(id_cell, Npts, pts);
 379     if (grid->GetCellType(id_cell) == VTK_TRIANGLE) {
 380       c2c[i].resize(3);
 381       nds->Reset();
 382       nds->InsertNextId(pts[0]);
 383       nds->InsertNextId(pts[1]);
 384       addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
 385       nds->Reset();
 386       nds->InsertNextId(pts[1]);
 387       nds->InsertNextId(pts[2]);
 388       addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
 389       nds->Reset();
 390       nds->InsertNextId(pts[2]);
 391       nds->InsertNextId(pts[0]);
 392       addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
 393     } else if (grid->GetCellType(id_cell) == VTK_QUAD) {
 394       c2c[i].resize(4);
 395       nds->Reset();
 396       nds->InsertNextId(pts[0]);
 397       nds->InsertNextId(pts[1]);
 398       addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
 399       nds->Reset();
 400       nds->InsertNextId(pts[1]);
 401       nds->InsertNextId(pts[2]);
 402       addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
 403       nds->Reset();
 404       nds->InsertNextId(pts[2]);
 405       nds->InsertNextId(pts[3]);
 406       addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
 407       nds->Reset();
 408       nds->InsertNextId(pts[3]);
 409       nds->InsertNextId(pts[0]);
 410       addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
 411     } else if (grid->GetCellType(id_cell) == VTK_TETRA) {
 412       c2c[i].resize(4);
 413       nds->Reset();
 414       nds->InsertNextId(pts[0]);
 415       nds->InsertNextId(pts[1]);
 416       nds->InsertNextId(pts[2]);
 417       addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
 418       nds->Reset();
 419       nds->InsertNextId(pts[0]);
 420       nds->InsertNextId(pts[1]);
 421       nds->InsertNextId(pts[3]);
 422       addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
 423       nds->Reset();
 424       nds->InsertNextId(pts[0]);
 425       nds->InsertNextId(pts[3]);
 426       nds->InsertNextId(pts[2]);
 427       addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
 428       nds->Reset();
 429       nds->InsertNextId(pts[1]);
 430       nds->InsertNextId(pts[2]);
 431       nds->InsertNextId(pts[3]);
 432       addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
 433     } else if (grid->GetCellType(id_cell) == VTK_PYRAMID) {
 434       c2c[i].resize(5);
 435       nds->Reset();
 436       nds->InsertNextId(pts[0]);
 437       nds->InsertNextId(pts[1]);
 438       nds->InsertNextId(pts[2]);
 439       nds->InsertNextId(pts[3]);
 440       addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
 441       nds->Reset();
 442       nds->InsertNextId(pts[0]);
 443       nds->InsertNextId(pts[1]);
 444       nds->InsertNextId(pts[4]);
 445       addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
 446       nds->Reset();
 447       nds->InsertNextId(pts[1]);
 448       nds->InsertNextId(pts[2]);
 449       nds->InsertNextId(pts[4]);
 450       addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
 451       nds->Reset();
 452       nds->InsertNextId(pts[2]);
 453       nds->InsertNextId(pts[3]);
 454       nds->InsertNextId(pts[4]);
 455       addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
 456       nds->Reset();
 457       nds->InsertNextId(pts[3]);
 458       nds->InsertNextId(pts[0]);
 459       nds->InsertNextId(pts[4]);
 460       addToC2C(id_cell, _cells, c2c, 4, nds, cls, grid);
 461     } else if (grid->GetCellType(id_cell) == VTK_WEDGE) {
 462       c2c[i].resize(5);
 463       nds->Reset();
 464       nds->InsertNextId(pts[0]);
 465       nds->InsertNextId(pts[1]);
 466       nds->InsertNextId(pts[2]);
 467       addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
 468       nds->Reset();
 469       nds->InsertNextId(pts[3]);
 470       nds->InsertNextId(pts[4]);
 471       nds->InsertNextId(pts[5]);
 472       addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
 473       nds->Reset();
 474       nds->InsertNextId(pts[0]);
 475       nds->InsertNextId(pts[1]);
 476       nds->InsertNextId(pts[4]);
 477       nds->InsertNextId(pts[3]);
 478       addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
 479       nds->Reset();
 480       nds->InsertNextId(pts[1]);
 481       nds->InsertNextId(pts[4]);
 482       nds->InsertNextId(pts[5]);
 483       nds->InsertNextId(pts[2]);
 484       addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
 485       nds->Reset();
 486       nds->InsertNextId(pts[0]);
 487       nds->InsertNextId(pts[2]);
 488       nds->InsertNextId(pts[5]);
 489       nds->InsertNextId(pts[3]);
 490       addToC2C(id_cell, _cells, c2c, 4, nds, cls, grid);
 491     } else if (grid->GetCellType(id_cell) == VTK_HEXAHEDRON) {
 492       c2c[i].resize(6);
 493       nds->Reset();
 494       nds->InsertNextId(pts[0]);
 495       nds->InsertNextId(pts[3]);
 496       nds->InsertNextId(pts[2]);
 497       nds->InsertNextId(pts[1]);
 498       addToC2C(id_cell, _cells, c2c, 0, nds, cls, grid);
 499       nds->Reset();
 500       nds->InsertNextId(pts[4]);
 501       nds->InsertNextId(pts[5]);
 502       nds->InsertNextId(pts[6]);
 503       nds->InsertNextId(pts[7]);
 504       addToC2C(id_cell, _cells, c2c, 1, nds, cls, grid);
 505       nds->Reset();
 506       nds->InsertNextId(pts[0]);
 507       nds->InsertNextId(pts[1]);
 508       nds->InsertNextId(pts[5]);
 509       nds->InsertNextId(pts[4]);
 510       addToC2C(id_cell, _cells, c2c, 2, nds, cls, grid);
 511       nds->Reset();
 512       nds->InsertNextId(pts[3]);
 513       nds->InsertNextId(pts[7]);
 514       nds->InsertNextId(pts[6]);
 515       nds->InsertNextId(pts[2]);
 516       addToC2C(id_cell, _cells, c2c, 3, nds, cls, grid);
 517       nds->Reset();
 518       nds->InsertNextId(pts[0]);
 519       nds->InsertNextId(pts[4]);
 520       nds->InsertNextId(pts[7]);
 521       nds->InsertNextId(pts[3]);
 522       addToC2C(id_cell, _cells, c2c, 4, nds, cls, grid);
 523       nds->Reset();
 524       nds->InsertNextId(pts[1]);
 525       nds->InsertNextId(pts[2]);
 526       nds->InsertNextId(pts[6]);
 527       nds->InsertNextId(pts[5]);
 528       addToC2C(id_cell, _cells, c2c, 5, nds, cls, grid);
 529     };
 530   };
 531 };
 532
 533 void EgVtkObject::getNodesFromCells
 534 (
 535   QVector<vtkIdType>  &cells,
 536   QVector<vtkIdType>  &nodes,
 537   vtkUnstructuredGrid *grid
 538 )
 539 {
 540   QSet<vtkIdType> ex_nodes;
 541   vtkIdType id_cell;
 542   foreach(id_cell, cells) {
 543     vtkIdType *pts;
 544     vtkIdType  Npts;
 545     grid->GetCellPoints(id_cell, Npts, pts);
 546     for (int i = 0; i < Npts; ++i) {
 547       ex_nodes.insert(pts[i]);
 548     };
 549   };
 550   nodes.resize(ex_nodes.size());
 551   {
 552     int j = 0;
 553     vtkIdType i;
 554     foreach(i,ex_nodes) {
 555       nodes[j] = i;
 556       ++j;
 557     };
 558   };
 559 };
 560
 561 bool EgVtkObject::isVolume(vtkUnstructuredGrid *grid, vtkIdType id_cell)
 562 {
 563   bool isVol = false;
 564   if      (grid->GetCellType(id_cell) == VTK_TETRA)      isVol = true;
 565   else if (grid->GetCellType(id_cell) == VTK_PYRAMID)    isVol = true;
 566   else if (grid->GetCellType(id_cell) == VTK_WEDGE)      isVol = true;
 567   else if (grid->GetCellType(id_cell) == VTK_HEXAHEDRON) isVol = true;
 568   return isVol;
 569 };
 570
 571 bool EgVtkObject::isSurface(vtkUnstructuredGrid *grid, vtkIdType id_cell)
 572 {
 573   bool isSurf = false;
 574   if      (grid->GetCellType(id_cell) == VTK_TRIANGLE) isSurf = true;
 575   else if (grid->GetCellType(id_cell) == VTK_QUAD)     isSurf = true;
 576   return isSurf;
 577 };
 578
 579 void EgVtkObject::UpdateCellIndex(vtkUnstructuredGrid *grid)
 580 {
 581   if (!grid->GetCellData()->GetArray("cell_index")) {
 582     EG_VTKSP(vtkLongArray_t, cell_index);
 583     cell_index->SetName("cell_index");
 584     cell_index->SetNumberOfValues(grid->GetNumberOfCells());
 585     grid->GetCellData()->AddArray(cell_index);
 586   };
 587   EG_VTKDCC(vtkLongArray_t, cell_index, grid, "cell_index");
 588   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 589     cell_index->SetValue(id_cell, id_cell);
 590   };
 591 };
 592
 593 void EgVtkObject::UpdateNodeIndex(vtkUnstructuredGrid *grid)
 594 {
 595   if (!grid->GetPointData()->GetArray("node_index")) {
 596     EG_VTKSP(vtkLongArray_t, node_index);
 597     node_index->SetName("node_index");
 598     node_index->SetNumberOfValues(grid->GetNumberOfPoints());
 599     grid->GetPointData()->AddArray(node_index);
 600   };
 601   EG_VTKDCN(vtkLongArray_t, node_index, grid, "node_index");
 602   for (vtkIdType pointId = 0; pointId < grid->GetNumberOfPoints(); ++pointId) {
 603     node_index->SetValue(pointId, pointId);
 604   };
 605 };
 606
 607 void EgVtkObject::addToPolyData
 608 (
 609   QVector<vtkIdType>  &cells,
 610   vtkPolyData         *pdata,
 611   vtkUnstructuredGrid *grid
 612 )
 613 {
 614   UpdateCellIndex(grid);
 615   UpdateNodeIndex(grid);
 616   QVector<vtkIdType> nodes;
 617   QVector<int>       _nodes;
 618   getNodesFromCells(cells, nodes, grid);
 619   createNodeMapping(nodes, _nodes, grid);
 620   EG_VTKSP(vtkDoubleArray, pcoords);
 621   pcoords->SetNumberOfComponents(3);
 622   pcoords->SetNumberOfTuples(nodes.size());
 623   EG_VTKSP(vtkPoints, points);
 624   points->SetData(pcoords);
 625   pdata->SetPoints(points);
 626   pdata->Allocate(cells.size());
 627   if (!pdata->GetCellData()->GetArray("cell_index")) {
 628     EG_VTKSP(vtkLongArray_t, cell_index);
 629     cell_index->SetName("cell_index");
 630     //cell_index->SetNumberOfValues(cells.size());
 631     pdata->GetCellData()->AddArray(cell_index);
 632   };
 633   if (!pdata->GetPointData()->GetArray("node_index")) {
 634     EG_VTKSP(vtkLongArray_t, node_index);
 635     node_index->SetName("node_index");
 636     //node_index->SetNumberOfValues(nodes.size());
 637     pdata->GetPointData()->AddArray(node_index);
 638   };
 639   EG_VTKDCC(vtkLongArray_t, pd_cell_index, pdata, "cell_index");
 640   EG_VTKDCN(vtkLongArray_t, pd_node_index, pdata, "node_index");
 641   pd_cell_index->SetNumberOfValues(cells.size());
 642   pd_node_index->SetNumberOfValues(nodes.size());
 643   for (int i_cell = 0; i_cell < cells.size(); ++i_cell) {
 644     vtkIdType id_cell = cells[i_cell];
 645     vtkIdType cellType = grid->GetCellType(id_cell);
 646     if ((cellType != VTK_TRIANGLE) && (cellType != VTK_QUAD)) {
 647       EG_ERR_RETURN("unsupported cell type for this operation");
 648     };
 649     vtkIdType Npts, *pts;
 650     grid->GetCellPoints(id_cell, Npts, pts);
 651     vtkIdType *new_pts = new vtkIdType[Npts];
 652     for (int i = 0; i < Npts; ++i) {
 653       new_pts[i] = _nodes[pts[i]];
 654     };
 655     vtkIdType newCellId = pdata->InsertNextCell(cellType, Npts, new_pts);
 656     pd_cell_index->SetValue(newCellId, id_cell);
 657     delete [] new_pts;
 658   };
 659   for (int i_node = 0; i_node < nodes.size(); ++i_node) {
 660     vec3_t x;
 661     grid->GetPoints()->GetPoint(nodes[i_node], x.data());
 662     pdata->GetPoints()->SetPoint(i_node, x.data());
 663     pd_node_index->SetValue(i_node, nodes[i_node]);
 664   };
 665 };
 666
 667 #define EGVTKOBJECT_COPYCELLDATA(FIELD,TYPE) \
 668 { \
 669 if (old_grid->GetCellData()->GetArray(FIELD)) { \
 670 EG_VTKDCC(TYPE, var1, old_grid, FIELD); \
 671 EG_VTKDCC(TYPE, var2, new_grid, FIELD); \
 672 var2->SetValue(newId, var1->GetValue(oldId)); \
 673 }; \
 674 };
 675
 676 void EgVtkObject::copyCellData
 677 (
 678   vtkUnstructuredGrid *old_grid,
 679   vtkIdType            oldId,
 680   vtkUnstructuredGrid *new_grid,
 681   vtkIdType            newId
 682 )
 683 {
 684   EGVTKOBJECT_COPYCELLDATA("vtk_type",      vtkIntArray);
 685   EGVTKOBJECT_COPYCELLDATA("cell_code",     vtkIntArray);
 686   EGVTKOBJECT_COPYCELLDATA("cell_index",    vtkLongArray_t);
 687   EGVTKOBJECT_COPYCELLDATA("cell_err_tet",  vtkDoubleArray);
 688   EGVTKOBJECT_COPYCELLDATA("cell_err_pria", vtkDoubleArray);
 689   EGVTKOBJECT_COPYCELLDATA("cell_err_prib", vtkDoubleArray);
 690   EGVTKOBJECT_COPYCELLDATA("cell_err_pric", vtkDoubleArray);
 691   EGVTKOBJECT_COPYCELLDATA("cell_err_prid", vtkDoubleArray);
 692   EGVTKOBJECT_COPYCELLDATA("cell_err_prie", vtkDoubleArray);
 693   EGVTKOBJECT_COPYCELLDATA("cell_err_prif", vtkDoubleArray);
 694 };
 695
 696 #define EGVTKOBJECT_COPYNODEDATA(FIELD,TYPE) \
 697 { \
 698 if (old_grid->GetPointData()->GetArray(FIELD)) { \
 699 EG_VTKDCN(TYPE, var1, old_grid, FIELD); \
 700 EG_VTKDCN(TYPE, var2, new_grid, FIELD); \
 701 var2->SetValue(newId, var1->GetValue(oldId)); \
 702 }; \
 703 };
 704
 705 void EgVtkObject::copyNodeData
 706 (
 707   vtkUnstructuredGrid *old_grid,
 708   vtkIdType            oldId,
 709   vtkUnstructuredGrid *new_grid,
 710   vtkIdType            newId
 711 )
 712 {
 713   EGVTKOBJECT_COPYNODEDATA("node_status", vtkIntArray);
 714   EGVTKOBJECT_COPYNODEDATA("node_layer",  vtkIntArray);
 715   EGVTKOBJECT_COPYNODEDATA("node_index",  vtkLongArray_t);
 716 };
 717
 718 #define EGVTKOBJECT_CREATECELLFIELD(FIELD,TYPE,OW) \
 719 if (!grid->GetCellData()->GetArray(FIELD)) { \
 720 EG_VTKSP(TYPE, var); \
 721 var->SetName(FIELD); \
 722 var->SetNumberOfValues(Ncells); \
 723 grid->GetCellData()->AddArray(var); \
 724 } else if (OW) { \
 725 EG_VTKDCC(TYPE, var, grid, FIELD); \
 726 var->SetNumberOfValues(Ncells); \
 727 };
 728
 729 #define EGVTKOBJECT_CREATENODEFIELD(FIELD,TYPE,OW) \
 730 if (!grid->GetPointData()->GetArray(FIELD)) { \
 731 EG_VTKSP(TYPE, var); \
 732 var->SetName(FIELD); \
 733 var->SetNumberOfValues(Nnodes); \
 734 grid->GetPointData()->AddArray(var); \
 735 for (int i = 0; i < grid->GetNumberOfPoints(); ++i) var->SetValue(i,0); \
 736 } else if (OW) { \
 737 EG_VTKDCN(TYPE, var, grid, FIELD); \
 738 var->SetNumberOfValues(Nnodes); \
 739 for (int i = 0; i < grid->GetNumberOfPoints(); ++i) var->SetValue(i,0); \
 740 };
 741
 742 void EgVtkObject::createBasicFields
 743 (
 744   vtkUnstructuredGrid *grid,
 745   vtkIdType            Ncells,
 746   vtkIdType            Nnodes,
 747   bool                 overwrite
 748 )
 749 {
 750   createBasicNodeFields(grid, Nnodes, overwrite);
 751   createBasicCellFields(grid, Ncells, overwrite);
 752 };
 753
 754 void EgVtkObject::createBasicCellFields
 755 (
 756   vtkUnstructuredGrid *grid,
 757   vtkIdType            Ncells,
 758   bool                 overwrite
 759 )
 760 {
 761   EGVTKOBJECT_CREATECELLFIELD("vtk_type" ,     vtkIntArray, overwrite);
 762   EGVTKOBJECT_CREATECELLFIELD("cell_code",     vtkIntArray, overwrite);
 763   EGVTKOBJECT_CREATECELLFIELD("cell_index",    vtkLongArray_t, overwrite);
 764   EGVTKOBJECT_CREATECELLFIELD("cell_err_tet",  vtkDoubleArray, overwrite);
 765   EGVTKOBJECT_CREATECELLFIELD("cell_err_pria", vtkDoubleArray, overwrite);
 766   EGVTKOBJECT_CREATECELLFIELD("cell_err_prib", vtkDoubleArray, overwrite);
 767   EGVTKOBJECT_CREATECELLFIELD("cell_err_pric", vtkDoubleArray, overwrite);
 768   EGVTKOBJECT_CREATECELLFIELD("cell_err_prid", vtkDoubleArray, overwrite);
 769   EGVTKOBJECT_CREATECELLFIELD("cell_err_prie", vtkDoubleArray, overwrite);
 770   EGVTKOBJECT_CREATECELLFIELD("cell_err_prif", vtkDoubleArray, overwrite);
 771   EGVTKOBJECT_CREATECELLFIELD("cell_VA",       vtkDoubleArray, overwrite);
 772 };
 773
 774 void EgVtkObject::createBasicNodeFields
 775 (
 776   vtkUnstructuredGrid *grid,
 777   vtkIdType            Nnodes,
 778   bool                 overwrite
 779 )
 780 {
 781   EGVTKOBJECT_CREATENODEFIELD("node_status", vtkIntArray,    overwrite);
 782   EGVTKOBJECT_CREATENODEFIELD("node_layer",  vtkIntArray,    overwrite);
 783   EGVTKOBJECT_CREATENODEFIELD("node_index",  vtkLongArray_t, overwrite);
 784 };
 785
 786 void EgVtkObject::allocateGrid
 787 (
 788   vtkUnstructuredGrid *grid,
 789   vtkIdType            Ncells,
 790   vtkIdType            Nnodes
 791 )
 792 {
 793   EG_VTKSP(vtkPoints,points);
 794   points->SetNumberOfPoints(Nnodes);
 795   grid->SetPoints(points);
 796   grid->Allocate(Ncells,max(vtkIdType(1),Ncells/10));
 797   createBasicFields(grid, Ncells, Nnodes);
 798 };
 799
 800 vec3_t EgVtkObject::cellCentre(vtkUnstructuredGrid *grid, vtkIdType id_cell)
 801 {
 802   vec3_t x,xc(0,0,0);
 803   vtkIdType *pts, N_pts;
 804   grid->GetCellPoints(id_cell, N_pts, pts);
 805   double f = 1.0/N_pts;
 806   for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 807     grid->GetPoint(pts[i_pts], x.data());
 808     xc += f*x;
 809   };
 810   return xc;
 811 };
 812
 813 void EgVtkObject::getRestCells(vtkUnstructuredGrid      *grid,
 814                                const QVector<vtkIdType> &cells,
 815                                QVector<vtkIdType>       &rest_cells)
 816 {
 817   QVector<bool> is_in_cells(grid->GetNumberOfCells(), false);
 818   foreach (vtkIdType id_cell, cells) {
 819     is_in_cells[id_cell] = true;
 820   };
 821   rest_cells.resize(grid->GetNumberOfCells() - cells.size());
 822   int i_rest_cells = 0;
 823   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 824     if (!is_in_cells[id_cell]) {
 825       rest_cells[i_rest_cells] = id_cell;
 826       ++i_rest_cells;
 827     };
 828   };
 829 };
 830
 831 void EgVtkObject::makeCopy(vtkUnstructuredGrid *src, vtkUnstructuredGrid *dst)
 832 {
 833   allocateGrid(dst, src->GetNumberOfCells(), src->GetNumberOfPoints());
 834   for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
 835     vec3_t x;
 836     src->GetPoints()->GetPoint(id_node, x.data());
 837     dst->GetPoints()->SetPoint(id_node, x.data());
 838     copyNodeData(src, id_node, dst, id_node);
 839   };
 840   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 841     vtkIdType N_pts, *pts;
 842     vtkIdType type_cell = src->GetCellType(id_cell);
 843     src->GetCellPoints(id_cell, N_pts, pts);
 844     vtkIdType id_new_cell = dst->InsertNextCell(type_cell, N_pts, pts);
 845     copyCellData(src, id_cell, dst, id_new_cell);
 846   };
 847 };
 848
 849 int EgVtkObject::findVolumeCell
 850 (
 851   vtkUnstructuredGrid      *grid,
 852   vtkIdType                id_surf,
 853   const QVector<int>       _nodes,
 854   const QVector<vtkIdType> cells,
 855   const QVector<int>       _cells,
 856   QVector<QSet<int> >       &n2c
 857 )
 858 {
 859   vtkIdType N_pts, *pts;
 860   grid->GetCellPoints(id_surf, N_pts, pts);
 861   QVector<QSet<int> > inters(N_pts-1);
 862   setIntersection(n2c[_nodes[pts[0]]], n2c[_nodes[pts[1]]], inters[0]);
 863   int i_pts = 2;
 864   while (i_pts < N_pts) {
 865     setIntersection(inters[i_pts-2], n2c[_nodes[pts[i_pts]]], inters[i_pts-1]);
 866     ++i_pts;
 867   };
 868   if (inters[N_pts-2].size() == 0) {
 869     return -1;
 870   } else if (inters[N_pts-2].size() > 2) {
 871     EG_BUG;
 872   };
 873   vtkIdType id_vol = -1;
 874   foreach (int i_cells, inters[N_pts-2]) {
 875     if (cells[i_cells] != id_surf) {
 876       id_vol = cells[i_cells];
 877     };
 878   };
 879   return id_vol;
 880 };
 881
 882 void EgVtkObject::setBoundaryCodes(const QSet<int> &bcs)
 883 {
 884   boundary_codes.clear();
 885   int bc;
 886   foreach(bc, bcs) {
 887     boundary_codes.insert(bc);
 888   };
 889 };
 890
 891 void EgVtkObject::createIndices(vtkUnstructuredGrid *grid)
 892 {
 893   if (!grid->GetCellData()->GetArray("cell_index")) {
 894     EG_VTKSP(vtkLongArray_t, var);
 895     var->SetName("cell_index");
 896     var->SetNumberOfValues(grid->GetNumberOfCells());
 897     grid->GetCellData()->AddArray(var);
 898   } else {
 899     EG_VTKDCC(vtkLongArray_t, var, grid, "cell_index");
 900     var->SetNumberOfValues(grid->GetNumberOfCells());
 901   };
 902   EG_VTKDCC(vtkLongArray_t, cell_index, grid, "cell_index");
 903   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 904     cell_index->SetValue(id_cell, id_cell);
 905   };
 906
 907   if (!grid->GetCellData()->GetArray("vtk_type")) {
 908     EG_VTKSP(vtkIntArray, var);
 909     var->SetName("vtk_type");
 910     var->SetNumberOfValues(grid->GetNumberOfCells());
 911     grid->GetCellData()->AddArray(var);
 912   } else {
 913     EG_VTKDCC(vtkIntArray, var, grid, "vtk_type");
 914     var->SetNumberOfValues(grid->GetNumberOfCells());
 915   };
 916   EG_VTKDCC(vtkIntArray, vtk_type, grid, "vtk_type");
 917   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 918     vtk_type->SetValue(id_cell, grid->GetCellType(id_cell));
 919   };
 920
 921   if (!grid->GetCellData()->GetArray("node_index")) {
 922     EG_VTKSP(vtkLongArray_t, var);
 923     var->SetName("node_index");
 924     var->SetNumberOfValues(grid->GetNumberOfPoints());
 925     grid->GetPointData()->AddArray(var);
 926   } else {
 927     EG_VTKDCC(vtkLongArray_t, var, grid, "node_index");
 928     var->SetNumberOfValues(grid->GetNumberOfPoints());
 929   };
 930   EG_VTKDCN(vtkLongArray_t, node_index, grid, "node_index");
 931   for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
 932     node_index->SetValue(id_node, id_node);
 933   };
 934 };
 935
 936 BoundaryCondition EgVtkObject::getBC(int bc)
 937 {
 938   return GuiMainWindow::pointer()->getBC(bc);
 939 };
 940
 941 int EgVtkObject::getSet(QString group, QString key, int value, int& variable)
 942 {
 943   QSettings *qset = GuiMainWindow::settings();
 944   QString typed_key = "int/" + key;
 945   qset->beginGroup(group);
 946   //if key=value pair not found in settings file, write it
 947   if (!qset->contains(typed_key)) qset->setValue(typed_key,value);
 948   //read key value from settings file and assign it to variable
 949   variable = (qset->value(typed_key,variable)).toInt();
 950   qset->endGroup();
 951   return(variable);
 952 }
 953
 954 double EgVtkObject::getSet(QString group, QString key, double value, double& variable)
 955 {
 956   QSettings *qset = GuiMainWindow::settings();
 957   QString typed_key = "double/" + key;
 958   qset->beginGroup(group);
 959   //if key=value pair not found in settings file, write it
 960   if (!qset->contains(typed_key)) qset->setValue(typed_key,value);
 961   //read key value from settings file and assign it to variable
 962   variable = (qset->value(typed_key,variable)).toDouble();
 963   qset->endGroup();
 964   return(variable);
 965 }
 966
 967 bool EgVtkObject::getSet(QString group, QString key, bool value, bool& variable)
 968 {
 969   QSettings *qset = GuiMainWindow::settings();
 970   QString typed_key = "bool/" + key;
 971   qset->beginGroup(group);
 972   Qt::CheckState state = (Qt::CheckState) ( value ? 2 : 0 );
 973   //if key=value pair not found in settings file, write it
 974   if (!qset->contains(typed_key)) qset->setValue(typed_key,state);
 975   //read key value from settings file and assign it to variable
 976   variable = (qset->value(typed_key,variable)).toBool();
 977   qset->endGroup();
 978   return(variable);
 979 }