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