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_specified_density", vtkIntArray);
 743   EGVTKOBJECT_COPYNODEDATA("node_meshdensity",  vtkDoubleArray);
 744   EGVTKOBJECT_COPYNODEDATA("node_meshdensity_current",  vtkDoubleArray);
 745   EGVTKOBJECT_COPYNODEDATA("node_type",  vtkCharArray);
 746 };
 747
 748 #define EGVTKOBJECT_CREATECELLFIELD(FIELD,TYPE,OW) \
 749 if (!grid->GetCellData()->GetArray(FIELD)) { \
 750 EG_VTKSP(TYPE, var); \
 751 var->SetName(FIELD); \
 752 var->SetNumberOfValues(Ncells); \
 753 grid->GetCellData()->AddArray(var); \
 754 } else if (OW) { \
 755 EG_VTKDCC(TYPE, var, grid, FIELD); \
 756 var->SetNumberOfValues(Ncells); \
 757 };
 758
 759 #define EGVTKOBJECT_CREATENODEFIELD(FIELD,TYPE,OW) \
 760 if (!grid->GetPointData()->GetArray(FIELD)) { \
 761 EG_VTKSP(TYPE, var); \
 762 var->SetName(FIELD); \
 763 var->SetNumberOfValues(Nnodes); \
 764 grid->GetPointData()->AddArray(var); \
 765 for (int i = 0; i < grid->GetNumberOfPoints(); ++i) var->SetValue(i,0); \
 766 } else if (OW) { \
 767 EG_VTKDCN(TYPE, var, grid, FIELD); \
 768 var->SetNumberOfValues(Nnodes); \
 769 for (int i = 0; i < grid->GetNumberOfPoints(); ++i) var->SetValue(i,0); \
 770 };
 771
 772 void EgVtkObject::createBasicFields
 773 (
 774   vtkUnstructuredGrid *grid,
 775   vtkIdType            Ncells,
 776   vtkIdType            Nnodes,
 777   bool                 overwrite
 778 )
 779 {
 780   createBasicNodeFields(grid, Nnodes, overwrite);
 781   createBasicCellFields(grid, Ncells, overwrite);
 782 };
 783
 784 void EgVtkObject::createBasicCellFields
 785 (
 786   vtkUnstructuredGrid *grid,
 787   vtkIdType            Ncells,
 788   bool                 overwrite
 789 )
 790 {
 791   EGVTKOBJECT_CREATECELLFIELD("vtk_type" ,     vtkIntArray, overwrite);
 792   EGVTKOBJECT_CREATECELLFIELD("cell_code",     vtkIntArray, overwrite);
 793   EGVTKOBJECT_CREATECELLFIELD("cell_index",    vtkLongArray_t, overwrite);
 794   EGVTKOBJECT_CREATECELLFIELD("cell_err_tet",  vtkDoubleArray, overwrite);
 795   EGVTKOBJECT_CREATECELLFIELD("cell_err_pria", vtkDoubleArray, overwrite);
 796   EGVTKOBJECT_CREATECELLFIELD("cell_err_prib", vtkDoubleArray, overwrite);
 797   EGVTKOBJECT_CREATECELLFIELD("cell_err_pric", vtkDoubleArray, overwrite);
 798   EGVTKOBJECT_CREATECELLFIELD("cell_err_prid", vtkDoubleArray, overwrite);
 799   EGVTKOBJECT_CREATECELLFIELD("cell_err_prie", vtkDoubleArray, overwrite);
 800   EGVTKOBJECT_CREATECELLFIELD("cell_err_prif", vtkDoubleArray, overwrite);
 801   EGVTKOBJECT_CREATECELLFIELD("cell_VA",       vtkDoubleArray, overwrite);
 802 };
 803
 804 void EgVtkObject::createBasicNodeFields
 805 (
 806   vtkUnstructuredGrid *grid,
 807   vtkIdType            Nnodes,
 808   bool                 overwrite
 809 )
 810 {
 811   EGVTKOBJECT_CREATENODEFIELD("node_status", vtkIntArray,    overwrite);
 812   EGVTKOBJECT_CREATENODEFIELD("node_layer",  vtkIntArray,    overwrite);
 813   EGVTKOBJECT_CREATENODEFIELD("node_index",  vtkLongArray_t, overwrite);
 814   EGVTKOBJECT_CREATENODEFIELD("node_specified_density", vtkIntArray,    overwrite);
 815   EGVTKOBJECT_CREATENODEFIELD("node_meshdensity",  vtkDoubleArray, overwrite);
 816   EGVTKOBJECT_CREATENODEFIELD("node_meshdensity_current",  vtkDoubleArray, overwrite);
 817   EGVTKOBJECT_CREATENODEFIELD("node_type",  vtkCharArray, overwrite);
 818 };
 819
 820 void EgVtkObject::allocateGrid
 821 (
 822   vtkUnstructuredGrid *grid,
 823   vtkIdType            Ncells,
 824   vtkIdType            Nnodes
 825 )
 826 {
 827   EG_VTKSP(vtkPoints,points);
 828   points->SetNumberOfPoints(Nnodes);
 829   grid->SetPoints(points);
 830   grid->Allocate(Ncells,max(vtkIdType(1),Ncells/10));
 831   createBasicFields(grid, Ncells, Nnodes);
 832 };
 833
 834 vec3_t EgVtkObject::cellCentre(vtkUnstructuredGrid *grid, vtkIdType id_cell)
 835 {
 836   vec3_t x,xc(0,0,0);
 837   vtkIdType *pts, N_pts;
 838   grid->GetCellPoints(id_cell, N_pts, pts);
 839   double f = 1.0/N_pts;
 840   for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
 841     grid->GetPoint(pts[i_pts], x.data());
 842     xc += f*x;
 843   };
 844   return xc;
 845 };
 846
 847 void EgVtkObject::getRestCells(vtkUnstructuredGrid      *grid,
 848                                const QVector<vtkIdType> &cells,
 849                                QVector<vtkIdType>       &rest_cells)
 850 {
 851   QVector<bool> is_in_cells(grid->GetNumberOfCells(), false);
 852   foreach (vtkIdType id_cell, cells) {
 853     is_in_cells[id_cell] = true;
 854   };
 855   rest_cells.resize(grid->GetNumberOfCells() - cells.size());
 856   int i_rest_cells = 0;
 857   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 858     if (!is_in_cells[id_cell]) {
 859       rest_cells[i_rest_cells] = id_cell;
 860       ++i_rest_cells;
 861     };
 862   };
 863 };
 864
 865 void EgVtkObject::makeCopy(vtkUnstructuredGrid *src, vtkUnstructuredGrid *dst)
 866 {
 867   allocateGrid(dst, src->GetNumberOfCells(), src->GetNumberOfPoints());
 868   for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
 869     vec3_t x;
 870     src->GetPoints()->GetPoint(id_node, x.data());
 871     dst->GetPoints()->SetPoint(id_node, x.data());
 872     copyNodeData(src, id_node, dst, id_node);
 873   };
 874   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 875     vtkIdType N_pts, *pts;
 876     vtkIdType type_cell = src->GetCellType(id_cell);
 877     src->GetCellPoints(id_cell, N_pts, pts);
 878     vtkIdType id_new_cell = dst->InsertNextCell(type_cell, N_pts, pts);
 879     copyCellData(src, id_cell, dst, id_new_cell);
 880   };
 881 };
 882
 883 void EgVtkObject::makeCopyNoAlloc(vtkUnstructuredGrid *src, vtkUnstructuredGrid *dst)
 884 {
 885   for (vtkIdType id_node = 0; id_node < src->GetNumberOfPoints(); ++id_node) {
 886     vec3_t x;
 887     src->GetPoints()->GetPoint(id_node, x.data());
 888     dst->GetPoints()->SetPoint(id_node, x.data());
 889     copyNodeData(src, id_node, dst, id_node);
 890   };
 891   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 892     vtkIdType N_pts, *pts;
 893     vtkIdType type_cell = src->GetCellType(id_cell);
 894     src->GetCellPoints(id_cell, N_pts, pts);
 895     vtkIdType id_new_cell = dst->InsertNextCell(type_cell, N_pts, pts);
 896     copyCellData(src, id_cell, dst, id_new_cell);
 897   };
 898 };
 899
 900 // void EgVtkObject::makeCopyNoAllocFiltered(vtkUnstructuredGrid *src, vtkUnstructuredGrid *dst, vector <bool> DeadNode, QVector <QSet <vtkIdType>> newCells)
 901 void EgVtkObject::makeCopyNoAllocFiltered(vtkUnstructuredGrid *src, vtkUnstructuredGrid *dst, vector <bool> DeadNode)
 902 {
 903   vtkIdType src_N_points=src->GetNumberOfPoints();
 904   vector <vtkIdType> OffSet(src_N_points);
 905   vtkIdType dst_id_node=0;
 906   for (vtkIdType src_id_node = 0; src_id_node < src_N_points; ++src_id_node) {
 907     if(!DeadNode[src_id_node])
 908     {
 909       vec3_t x;
 910       src->GetPoints()->GetPoint(src_id_node, x.data());
 911       dst->GetPoints()->SetPoint(dst_id_node, x.data());
 912       copyNodeData(src, src_id_node, dst, dst_id_node);
 913       OffSet[src_id_node]=src_id_node-dst_id_node;
 914       dst_id_node++;
 915     }
 916     else
 917     {
 918       //search closest node
 919     }
 920   };
 921   for (vtkIdType id_cell = 0; id_cell < src->GetNumberOfCells(); ++id_cell) {
 922     vtkIdType N_pts, *src_pts, *dst_pts;
 923     vtkIdType type_cell = src->GetCellType(id_cell);
 924     src->GetCellPoints(id_cell, N_pts, src_pts);
 925     src->GetCellPoints(id_cell, N_pts, dst_pts);
 926     bool DeadCell=false;
 927     for(int i=0;i<N_pts;i++)
 928     {
 929       if(DeadNode[src_pts[i]]) {DeadCell=true;}
 930       dst_pts[i]=src_pts[i]-OffSet[src_pts[i]];
 931     }
 932     if(!DeadCell)
 933     {
 934       vtkIdType id_new_cell = dst->InsertNextCell(type_cell, N_pts, dst_pts);
 935       copyCellData(src, id_cell, dst, id_new_cell);
 936     }
 937   };
 938 }
 939
 940 int EgVtkObject::findVolumeCell
 941 (
 942   vtkUnstructuredGrid      *grid,
 943   vtkIdType                id_surf,
 944   const QVector<int>       _nodes,
 945   const QVector<vtkIdType> cells,
 946   const QVector<int>       _cells,
 947   QVector<QSet<int> >       &n2c
 948 )
 949 {
 950   vtkIdType N_pts, *pts;
 951   if (_cells.size()) N_pts = N_pts; // dummy statement to get rid of compiler warning ...
 952   grid->GetCellPoints(id_surf, N_pts, pts);
 953   QVector<QSet<int> > inters(N_pts-1);
 954   setIntersection(n2c[_nodes[pts[0]]], n2c[_nodes[pts[1]]], inters[0]);
 955   int i_pts = 2;
 956   while (i_pts < N_pts) {
 957     setIntersection(inters[i_pts-2], n2c[_nodes[pts[i_pts]]], inters[i_pts-1]);
 958     ++i_pts;
 959   };
 960   if (inters[N_pts-2].size() == 0) {
 961     return -1;
 962   } else if (inters[N_pts-2].size() > 2) {
 963     EG_BUG;
 964   };
 965   vtkIdType id_vol = -1;
 966   foreach (int i_cells, inters[N_pts-2]) {
 967     if (cells[i_cells] != id_surf) {
 968       id_vol = cells[i_cells];
 969     };
 970   };
 971   return id_vol;
 972 };
 973
 974 void EgVtkObject::setBoundaryCodes(const QSet<int> &bcs)
 975 {
 976   boundary_codes.clear();
 977   int bc;
 978   foreach(bc, bcs) {
 979     boundary_codes.insert(bc);
 980   };
 981 };
 982
 983 void EgVtkObject::createIndices(vtkUnstructuredGrid *grid)
 984 {
 985   if (!grid->GetCellData()->GetArray("cell_index")) {
 986     EG_VTKSP(vtkLongArray_t, var);
 987     var->SetName("cell_index");
 988     var->SetNumberOfValues(grid->GetNumberOfCells());
 989     grid->GetCellData()->AddArray(var);
 990   } else {
 991     EG_VTKDCC(vtkLongArray_t, var, grid, "cell_index");
 992     var->SetNumberOfValues(grid->GetNumberOfCells());
 993   };
 994   EG_VTKDCC(vtkLongArray_t, cell_index, grid, "cell_index");
 995   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 996     cell_index->SetValue(id_cell, id_cell);
 997   };
 998
 999   if (!grid->GetCellData()->GetArray("vtk_type")) {
1000     EG_VTKSP(vtkIntArray, var);
1001     var->SetName("vtk_type");
1002     var->SetNumberOfValues(grid->GetNumberOfCells());
1003     grid->GetCellData()->AddArray(var);
1004   } else {
1005     EG_VTKDCC(vtkIntArray, var, grid, "vtk_type");
1006     var->SetNumberOfValues(grid->GetNumberOfCells());
1007   };
1008   EG_VTKDCC(vtkIntArray, vtk_type, grid, "vtk_type");
1009   for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
1010     vtk_type->SetValue(id_cell, grid->GetCellType(id_cell));
1011   };
1012
1013   if (!grid->GetCellData()->GetArray("node_index")) {
1014     EG_VTKSP(vtkLongArray_t, var);
1015     var->SetName("node_index");
1016     var->SetNumberOfValues(grid->GetNumberOfPoints());
1017     grid->GetPointData()->AddArray(var);
1018   } else {
1019     EG_VTKDCC(vtkLongArray_t, var, grid, "node_index");
1020     var->SetNumberOfValues(grid->GetNumberOfPoints());
1021   };
1022   EG_VTKDCN(vtkLongArray_t, node_index, grid, "node_index");
1023   for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
1024     node_index->SetValue(id_node, id_node);
1025   };
1026 };
1027
1028 BoundaryCondition EgVtkObject::getBC(int bc)
1029 {
1030   return GuiMainWindow::pointer()->getBC(bc);
1031 };
1032
1033 int EgVtkObject::getSet(QString group, QString key, int value, int& variable)
1034 {
1035   QSettings *qset = GuiMainWindow::settings();
1036   QString typed_key = "int/" + key;
1037   qset->beginGroup(group);
1038   //if key=value pair not found in settings file, write it
1039   if (!qset->contains(typed_key)) qset->setValue(typed_key,value);
1040   //read key value from settings file and assign it to variable
1041   variable = (qset->value(typed_key,variable)).toInt();
1042   qset->endGroup();
1043   return(variable);
1044 }
1045
1046 double EgVtkObject::getSet(QString group, QString key, double value, double& variable)
1047 {
1048   QSettings *qset = GuiMainWindow::settings();
1049   QString typed_key = "double/" + key;
1050   qset->beginGroup(group);
1051   //if key=value pair not found in settings file, write it
1052   if (!qset->contains(typed_key)) qset->setValue(typed_key,value);
1053   //read key value from settings file and assign it to variable
1054   variable = (qset->value(typed_key,variable)).toDouble();
1055   qset->endGroup();
1056   return(variable);
1057 }
1058
1059 bool EgVtkObject::getSet(QString group, QString key, bool value, bool& variable)
1060 {
1061   QSettings *qset = GuiMainWindow::settings();
1062   QString typed_key = "bool/" + key;
1063   qset->beginGroup(group);
1064   Qt::CheckState state = (Qt::CheckState) ( value ? 2 : 0 );
1065   //if key=value pair not found in settings file, write it
1066   if (!qset->contains(typed_key)) qset->setValue(typed_key,state);
1067   //read key value from settings file and assign it to variable
1068   variable = (qset->value(typed_key,variable)).toBool();
1069   qset->endGroup();
1070   return(variable);
1071 }
1072
1073 ///////////////////////////////////////////
1074 int cout_grid(ostream &stream, vtkUnstructuredGrid *grid, bool npoints, bool ncells, bool points, bool cells)
1075 {
1076   stream<<"============="<<endl;
1077   if(npoints) stream << "grid->GetNumberOfPoints()=" << grid->GetNumberOfPoints() << endl;
1078   if(ncells) stream << "grid->GetNumberOfCells()=" << grid->GetNumberOfCells() << endl;
1079   if(points) {
1080     for (vtkIdType i = 0; i < grid->GetNumberOfPoints(); ++i) {
1081       vec3_t x;
1082       grid->GetPoint(i, x.data());
1083       stream << "Vertex " << i << " = " << x << endl;
1084     };
1085   }
1086   if(cells) {
1087     EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
1088     for (vtkIdType i = 0; i < grid->GetNumberOfCells(); ++i) {
1089       vtkCell *C = (vtkCell *) vtkCell::New();
1090       C=grid->GetCell(i);
1091       vtkIdType npts=C->GetNumberOfPoints();
1092       vtkIdType* pts;
1093       grid->GetCellPoints(i, npts, pts);
1094       stream << "Cell " << i << " = ";
1095       for(int j=0;j<npts;j++) stream << pts[j] << " ";
1096       stream << "boundary_code=" << cell_code->GetValue(i);
1097       stream << endl;
1098     };
1099   }
1100   stream<<"============="<<endl;
1101   return 0;
1102 }
1103
1104 ///////////////////////////////////////////
1105 int addPoint(vtkUnstructuredGrid* a_grid,vtkIdType index,vec3_t x)
1106 {
1107   a_grid->GetPoints()->SetPoint(index, x.data());
1108   return(0);
1109 }
1110
1111 int addCell(vtkUnstructuredGrid* a_grid, vtkIdType A, vtkIdType B, vtkIdType C, int bc)
1112 {
1113   vtkIdType npts=3;
1114   vtkIdType pts[3];
1115   pts[0]=A;
1116   pts[1]=B;
1117   pts[2]=C;
1118   vtkIdType newCellId = a_grid->InsertNextCell(VTK_TRIANGLE,npts,pts);
1119   EG_VTKDCC(vtkIntArray, cell_code, a_grid, "cell_code");
1120   cell_code->SetValue(newCellId, bc);
1121   return(0);
1122 }
1123
1124
1125 ///////////////////////////////////////////
1126
1127 int getShortestSide(vtkIdType a_id_cell,vtkUnstructuredGrid* a_grid)
1128 {
1129   vtkIdType N_pts, *pts;
1130   a_grid->GetCellPoints(a_id_cell, N_pts, pts);
1131   vec3_t* x=new vec3_t[N_pts];
1132   for(int i=0;i<N_pts;i++) a_grid->GetPoints()->GetPoint(pts[i], x[i].data());
1133   int id_minlen=0;
1134   double minlen=(x[1]-x[0]).abs();
1135   for(int i=1;i<N_pts;i++)
1136   {
1137     double len=(x[(i+1)%N_pts]-x[i]).abs();
1138     if(len<minlen){
1139       minlen=len;
1140       id_minlen=i;
1141     }
1142   }
1143   delete x;
1144   return(id_minlen);
1145 }
1146
1147 int getLongestSide(vtkIdType a_id_cell,vtkUnstructuredGrid* a_grid)
1148 {
1149   vtkIdType N_pts, *pts;
1150   a_grid->GetCellPoints(a_id_cell, N_pts, pts);
1151   vec3_t* x=new vec3_t[N_pts];
1152   for(int i=0;i<N_pts;i++) a_grid->GetPoints()->GetPoint(pts[i], x[i].data());
1153   int id_maxlen=0;
1154   double maxlen=(x[1]-x[0]).abs();
1155   cout<<"maxlen="<<maxlen<<endl;
1156   for(int i=1;i<N_pts;i++)
1157   {
1158     double len=(x[(i+1)%N_pts]-x[i]).abs();
1159     cout<<"len["<<i<<"]="<<len<<endl;
1160     if(len>maxlen){
1161       maxlen=len;
1162       id_maxlen=i;
1163     }
1164   }
1165   delete x;
1166   return(id_maxlen);
1167 }
1168
1169 //get number of the edge corresponding to node1-node2
1170 int getSide(vtkIdType a_id_cell,vtkUnstructuredGrid* a_grid,vtkIdType a_id_node1,vtkIdType a_id_node2)
1171 {
1172   vtkIdType N_pts, *pts;
1173   a_grid->GetCellPoints(a_id_cell, N_pts, pts);
1174   QVector <vtkIdType> edge(2);
1175
1176   int n=0;
1177   for(int i=0;i<N_pts;i++)
1178   {
1179     if(pts[i]==a_id_node1) { edge[0]=i;n++;}
1180     if(pts[i]==a_id_node2) { edge[1]=i;n++;}
1181   }
1182   if(n!=2){
1183     EG_BUG;
1184     return(-1);
1185   }
1186   qSort(edge.begin(),edge.end());
1187   if(edge[0]==0 && edge[1]==N_pts-1) return(N_pts-1);
1188   else return(edge[0]);
1189 }
1190 ///////////////////////////////////////////
1191
1192 QSet <int> complementary_bcs(QSet <int> &bcs, vtkUnstructuredGrid *a_grid, QVector <vtkIdType> &a_cells)
1193 {
1194   QSet <int> bcs_complement;
1195   EG_VTKDCC(vtkIntArray, bc, a_grid, "cell_code");
1196   foreach (vtkIdType id_cell, a_cells) {
1197     int code=bc->GetValue(id_cell);
1198     if (!bcs.contains(code)) {
1199       bcs_complement.insert(code);
1200     };
1201   }
1202   return(bcs_complement);
1203 }
1204
1205 ///////////////////////////////////////////
1206
1207 QString cell2str(vtkIdType id_cell,vtkUnstructuredGrid* grid)
1208 {
1209   QString tmp;
1210   tmp.setNum(id_cell);
1211   QString txt = "id_cell=" + tmp;
1212
1213   vtkIdType N_pts, *pts;
1214   grid->GetCellPoints(id_cell, N_pts, pts);
1215
1216   txt += " [";
1217   for (int i_pts = 0; i_pts < N_pts; ++i_pts) {
1218     tmp.setNum(pts[i_pts]);
1219     txt += tmp;
1220     if (i_pts < N_pts-1) {
1221       txt += ",";
1222     };
1223   };
1224   txt += "]";
1225   return(txt);
1226 }
1227
1228
1229
1230 ///////////////////////////////////////////
1231
1232 Qt::CheckState int2CheckState(int a)
1233 {
1234   if(a==0) return(Qt::Unchecked);
1235   if(a==1) return(Qt::PartiallyChecked);
1236   else return(Qt::Checked);
1237 }
1238
1239 int CheckState2int(Qt::CheckState a)
1240 {
1241   if(a==Qt::Unchecked) return(0);
1242   if(a==Qt::PartiallyChecked) return(1);
1243   else return(2);
1244 }
1245
1246 // ///////////////////////////////////////////
1247 // /* Here is how we we get QTextStreams that look like iostreams */
1248 // Qcin=QTextStream(stdin, QIODevice::ReadOnly);
1249 // Qcout=QTextStream(stdout, QIODevice::WriteOnly);
1250 // Qcerr=QTextStream(stderr, QIODevice::WriteOnly);
1251 // ///////////////////////////////////////////
1252
1253 pair<vtkIdType,vtkIdType> OrderedPair(vtkIdType a, vtkIdType b)
1254 {
1255   vtkIdType x=min(a,b);
1256   vtkIdType y=max(a,b);
1257   return(pair<vtkIdType,vtkIdType>(x,y));
1258 }
1259
1260 vtkIdType nextcell(vtkIdType a_cell, vtkIdType a_node, QVector< QVector< int > > &a_c2c, vtkUnstructuredGrid *a_grid)
1261 {
1262   vtkIdType N_pts, *pts;
1263   a_grid->GetCellPoints(a_cell, N_pts, pts);
1264
1265   int i;
1266   for(i=0;i<N_pts;i++)
1267   {
1268     if(pts[i]==a_node) break;
1269   }
1270   return a_c2c[a_cell][i];
1271 }
1272
1273 const char* VertexType2Str(char T)
1274 {
1275   if(T==VTK_SIMPLE_VERTEX) return("VTK_SIMPLE_VERTEX");
1276   if(T==VTK_FIXED_VERTEX) return("VTK_FIXED_VERTEX");
1277   if(T==VTK_FEATURE_EDGE_VERTEX) return("VTK_FEATURE_EDGE_VERTEX");
1278   if(T==VTK_BOUNDARY_EDGE_VERTEX) return("VTK_BOUNDARY_EDGE_VERTEX");
1279   else return("Unknown vertex type");
1280 }
1281
1282 char Str2VertexType(QString S)
1283 {
1284   if(S=="VTK_SIMPLE_VERTEX") return((char)0);
1285   if(S=="VTK_FIXED_VERTEX") return((char)1);
1286   if(S=="VTK_FEATURE_EDGE_VERTEX") return((char)2);
1287   if(S=="VTK_BOUNDARY_EDGE_VERTEX") return((char)3);
1288   else return((char)-1);
1289 }
1290
1291 const char* vertex_type(char T)
1292 {
1293   if(T==VTK_SIMPLE_VERTEX) return("VTK_SIMPLE_VERTEX");
1294   if(T==VTK_FIXED_VERTEX) return("VTK_FIXED_VERTEX");
1295   if(T==VTK_FEATURE_EDGE_VERTEX) return("VTK_FEATURE_EDGE_VERTEX");
1296   if(T==VTK_BOUNDARY_EDGE_VERTEX) return("VTK_BOUNDARY_EDGE_VERTEX");
1297   else return("Unknown vertex type");
1298 }