src/libengrid/meshpartition.h

   1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   2 // +                                                                      +
   3 // + This file is part of enGrid.                                         +
   4 // +                                                                      +
   5 // + Copyright 2008-2014 enGits GmbH                                      +
   6 // +                                                                      +
   7 // + enGrid is free software: you can redistribute it and/or modify       +
   8 // + it under the terms of the GNU General Public License as published by +
   9 // + the Free Software Foundation, either version 3 of the License, or    +
  10 // + (at your option) any later version.                                  +
  11 // +                                                                      +
  12 // + enGrid is distributed in the hope that it will be useful,            +
  13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  15 // + GNU General Public License for more details.                         +
  16 // +                                                                      +
  17 // + You should have received a copy of the GNU General Public License    +
  18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  19 // +                                                                      +
  20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  21
  22 #ifndef MESHPARTITION_H
  23 #define MESHPARTITION_H
  24
  25 class MeshPartition;
  26
  27 #include "egvtkobject.h"
  28
  29 class MeshPartition : public EgVtkObject
  30 {
  31
  32 private: // attributes
  33
  34
  35   vtkUnstructuredGrid*   m_Grid;   ///< the grid underlying this mesh partition
  36   QVector<vtkIdType>     m_Cells;  ///< all cells of the mesh partition
  37   QVector<int>           m_LCells; ///< inverse indexing for the cells
  38   QVector<vtkIdType>     m_Nodes;  ///< all nodes of the mesh partition
  39   QVector<int>           m_LNodes; ///< inverse indexing for the nodes
  40   QVector<QVector<int> > m_N2C;    ///< node to cell information
  41   QVector<QVector<int> > m_N2BC;   ///< node to boundary code information
  42   QVector<QVector<int> > m_N2N;    ///< node to node information
  43   QVector<QVector<int> > m_C2C;    ///< cell to cell information
  44
  45   int m_CellsStamp;  ///< "time"-stamp
  46   int m_LCellsStamp; ///< "time"-stamp
  47   int m_NodesStamp;  ///< "time"-stamp
  48   int m_LNodesStamp; ///< "time"-stamp
  49   int m_N2NStamp;    ///< "time"-stamp
  50   int m_N2CStamp;    ///< "time"-stamp
  51   int m_N2BCStamp;   ///< "time"-stamp
  52   int m_C2CStamp;    ///< "time"-stamp
  53
  54   bool m_TrackGrid;              ///< flag to determine if grid should be tracked (make sure that all cells are always included)
  55   unsigned long int m_GridMTime; ///< VTK's modification time of the underlying grid
  56
  57 private: // methods
  58
  59   void createNodeToBC();
  60
  61   void resetTimeStamps();
  62   void checkCells();
  63   void checkNodes();
  64   void checkLCells();
  65   void checkLNodes();
  66   void checkN2N();
  67   void checkN2C();
  68   void checkN2BC();
  69   void checkC2C();
  70
  71 public: // methods
  72
  73   /// Create an empty (undefined) mesh partition
  74   MeshPartition();
  75
  76   /**
  77    * Create a mesh partition with the grid set. Optionally all cells can be selected.
  78    * @param grid the grid to use
  79    * @param use_all_cells if set to true all cells will be selected;
  80    */
  81   MeshPartition(vtkUnstructuredGrid *grid, bool use_all_cells = false);
  82
  83   /**
  84    * Create a mesh partition from a global volume definition
  85    * @param volume_name the name of the volume
  86    */
  87   MeshPartition(QString volume_name);
  88
  89   /**
  90    * Set the grid.
  91    * @param a pointer to the grid
  92    */
  93   void setGrid(vtkUnstructuredGrid *grid) { m_Grid = grid; }
  94
  95   /**
  96    * Set the grid and make sure all cells are always included (automatic tracking).
  97    * @param a pointer to the grid
  98    */
  99   void trackGrid(vtkUnstructuredGrid *grid);
 100
 101   /**
 102    * Access to the grid.
 103    * @return a pointer to the grid
 104    */
 105   vtkUnstructuredGrid* getGrid() const { return m_Grid; }
 106
 107   /**
 108    * Define the mesh partition by defining its cells.
 109    * @param cls the cells of the subset
 110    */
 111   template <class C>
 112   void setCells(const C& cls);
 113
 114   /**
 115    * Define the mesh partition by defining its nodes.
 116    * @param nds the nodes of the subset
 117    */
 118   template <class C>
 119   void setNodes(const C& nds);
 120
 121   /**
 122    * Define the mesh partition by defining boundary codes.
 123    * @param bcs the boundary codes of the subset
 124    */
 125   template <class C>
 126   void setBCs(const C& bcs);
 127
 128   /**
 129    * Define the mesh partition by defining all its cells.
 130    */
 131   void setAllCells();
 132
 133   /**
 134    * Define the mesh partition by giving a symbolic volume name.
 135    * The grid will be changed to the default (main) grid that is currently loaded into ENGRID.
 136    * @param volume_name the symbolic volume name
 137    */
 138   void setVolume(QString volume_name);
 139
 140   /**
 141    * Define the mesh partition as the remainder of an existing partition.
 142    * @param part the existing partition
 143    */
 144   void setRemainder(const MeshPartition& part);
 145
 146   const QVector<vtkIdType>&     getCells() const; ///< Access to the cell indices
 147   const QVector<int>&           getLocalCells();  ///< Access to the local cell indices
 148   const QVector<vtkIdType>&     getNodes();       ///< Access to the node indices
 149   const QVector<int>&           getLocalNodes();  ///< Access to the local node indices
 150   const QVector<QVector<int> >& getN2N();         ///< Access to the local node to node structure
 151   const QVector<QVector<int> >& getN2C();         ///< Access to the local node to cell structure
 152   const QVector<QVector<int> >& getC2C();         ///< Access to the local cell to cell structure
 153
 154   void setVolumeOrientation();   ///< change the face orientation to match the volume definition
 155   void setOriginalOrientation(); ///< change the orientation to match the original orientation
 156
 157   /**
 158    * Copy the partition to a VTK grid.
 159    * @param new_grid the grid to copy the partition to (will be resized accordingly).
 160    */
 161   void extractToVtkGrid(vtkUnstructuredGrid *new_grid);
 162
 163   /**
 164    * Add another partition to this one.
 165    * At the moment overlapping partitions on two different grids will not be handled well.
 166    * If both partitions do not have the same underlying grid the grid will be extended in order
 167    * to add the other partition.
 168    * @param part the partition to add
 169    * @param tol the tolerance to identify duplicate nodes
 170    *        (negative values denote a relative tolerance -- relative to the smallest edge length)
 171    */
 172   void addPartition(const MeshPartition& part, double tol = -1e-3);
 173
 174   /**
 175    * compute the smallest edge length of the partition
 176    * @return the smallest edge length
 177    */
 178   double getSmallestEdgeLength() const;
 179
 180   /**
 181    * Get the number of nodes in the partition.
 182    * @return the number of nodes
 183    */
 184   int getNumberOfNodes();
 185
 186   /**
 187    * Get the number of cells in the partition.
 188    * @return the number of cells
 189    */
 190   int getNumberOfCells();
 191
 192   /**
 193    * Get the average length of all surface edges connected to this node.
 194    * @param id_node the node ID of the node in question
 195    * @return the average length of all connected surface edges
 196    */
 197   double getAverageSurfaceEdgeLength(vtkIdType id_node);
 198
 199   /**
 200    * @brief compute the minimal and maximal edge length of a surface stencil
 201    * A surface stencil consists of all surface elements which have a single node in common.
 202    * @param id_node the node in common
 203    * @param l_min on return this will hold the minimal edge length
 204    * @param l_max on return this will hold the maximal edge length
 205    */
 206   void computeMinAndMaxSurfaceStencilEdgeLengths(vtkIdType id_node, double &l_min, double &l_max);
 207
 208   /**
 209    * @brief get the minimal edge length of a surface stencil
 210    * A surface stencil consists of all surface elements which have a single node in common.
 211    * @param id_node the node in common
 212    * @return the minimal edge length
 213    */
 214   double getMinSurfaceStencilEdgeLength(vtkIdType id_node);
 215
 216   /**
 217    * @brief get the maximal edge length of a surface stencil
 218    * A surface stencil consists of all surface elements which have a single node in common.
 219    * @param id_node the node in common
 220    * @return the maximal edge length
 221    */
 222   double getMaxSurfaceStencilEdgeLength(vtkIdType id_node);
 223
 224   vtkIdType getVolumeCell(vtkIdType id_face);
 225
 226   int       localNode(vtkIdType id_node);
 227   vtkIdType globalNode(int i);
 228   int       localCell(vtkIdType id_cell);
 229   vtkIdType globalCell(int i);
 230
 231
 232   int       n2nLSize(int i_nodes);
 233   int       n2nLL(int i_nodes, int j);
 234   vtkIdType n2nLG(int i_nodes, int j);
 235   int       n2nGSize(vtkIdType id_node);
 236   int       n2nGL(vtkIdType id_node, int j);
 237   vtkIdType n2nGG(vtkIdType id_node, int j);
 238   int       n2cLSize(int i_nodes);
 239   int       n2cLL(int i_nodes, int j);
 240   vtkIdType n2cLG(int i_nodes, int j);
 241   int       n2cGSize(vtkIdType id_node);
 242   int       n2cGL(vtkIdType id_node, int j);
 243   vtkIdType n2cGG(vtkIdType id_node, int j);
 244   int       c2cLSize(int i_cells);
 245   int       c2cLL(int i_cells, int j);
 246   vtkIdType c2cLG(int i_cells, int j);
 247   int       c2cGSize(vtkIdType id_cell);
 248   int       c2cGL(vtkIdType id_cell, int j);
 249   vtkIdType c2cGG(vtkIdType id_cell, int j);
 250   int       n2bcLSize(int i_nodes);
 251   int       n2bcL(int i_nodes, int j);
 252   int       n2bcGSize(vtkIdType id_node);
 253   int       n2bcG(vtkIdType id_node, int j);
 254
 255   bool hasNeighNode(vtkIdType id_node, vtkIdType id_neigh);
 256   bool hasBC(vtkIdType id_node, int bc);
 257
 258   /**
 259    * Compute the normal vector of a node.
 260    * @param id_node the global ID of the node
 261    * @return the normalised normal vector
 262    */
 263   vec3_t globalNormal(vtkIdType id_node);
 264
 265   template <typename C>
 266   void getGlobalN2N(vtkIdType id_node, C& cont);
 267
 268   int getNumberOfFeatureNeighbours(vtkIdType id_node);
 269
 270   template <typename C>
 271   void getEdgeFaces(vtkIdType id_node1, vtkIdType id_node2, C &edge_faces);
 272
 273   int getEdgeType(vtkIdType id_node1, vtkIdType id_node2);
 274
 275   /**
 276    * @brief compute topological distance between two nodes
 277    * @param id_node1 index of the first node
 278    * @param id_node2 index of the second node
 279    * @param max_dist maximal search distance
 280    * @param restriction_type (0: no restriction, 1: only surface nodes, 2: only edge nodes)
 281    * @return the number of edges for the shortest connection between the two nodes
 282    */
 283   int computeTopoDistance(vtkIdType id_node1, vtkIdType id_node2, int max_dist, int restriction_type);
 284
 285   template <typename C>
 286   void getCommonNodes(vtkIdType id_cell1, vtkIdType id_cell2, C &common_nodes);
 287
 288 };
 289
 290
 291 template <class C>
 292 inline void MeshPartition::setCells(const C& cls)
 293 {
 294   m_Cells.resize(cls.size());
 295   qCopy(cls.begin(), cls.end(), m_Cells.begin());
 296   ++m_CellsStamp;
 297 }
 298
 299 template <class C>
 300 inline void MeshPartition::setNodes(const C& nds)
 301 {
 302   QList<vtkIdType> cls;
 303   QVector<bool> node_inside(m_Grid->GetNumberOfPoints(), false);
 304   foreach (vtkIdType id_node, nds) {
 305     node_inside[id_node] = true;
 306   }
 307   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 308     vtkIdType N_pts, *pts;
 309     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 310     bool append_cell = true;
 311     for (int i = 0; i < N_pts; ++i) {
 312       if (!node_inside[pts[i]]) {
 313         append_cell = false;
 314         break;
 315       }
 316     }
 317     if (append_cell) {
 318       cls.append(id_cell);
 319     }
 320   }
 321   setCells(cls);
 322 }
 323
 324 template <class C>
 325 inline void MeshPartition::setBCs(const C& bcs)
 326 {
 327   QList<vtkIdType> cls;
 328   EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
 329   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 330     foreach (int bc, bcs) {
 331       if (cell_code->GetValue(id_cell) == bc) {
 332         cls.append(id_cell);
 333         break;
 334       }
 335     }
 336   }
 337   setCells(cls);
 338 }
 339
 340 inline void MeshPartition::setAllCells()
 341 {
 342   QVector <vtkIdType> all_cells;
 343   getAllCells(all_cells, m_Grid);
 344   this->setCells(all_cells);
 345 }
 346
 347 inline void MeshPartition::checkCells()
 348 {
 349   if (m_Grid->GetMTime() > m_GridMTime) {
 350     setAllCells();
 351     m_GridMTime = m_Grid->GetMTime();
 352   }
 353 }
 354
 355 inline void MeshPartition::checkNodes()
 356 {
 357   if (m_TrackGrid) {
 358     checkCells();
 359   }
 360   if (m_CellsStamp > m_NodesStamp) {
 361     getNodesFromCells(m_Cells, m_Nodes, m_Grid);
 362     m_NodesStamp = m_CellsStamp;
 363   }
 364 }
 365
 366 inline void MeshPartition::checkLCells()
 367 {
 368   if (m_TrackGrid) {
 369     checkCells();
 370   }
 371   if (m_CellsStamp > m_LCellsStamp) {
 372     createCellMapping(m_Cells, m_LCells, m_Grid);
 373     m_LCellsStamp = m_CellsStamp;
 374   }
 375 }
 376
 377 inline void MeshPartition::checkLNodes()
 378 {
 379   checkNodes();
 380   if (m_NodesStamp > m_LNodesStamp) {
 381     createNodeMapping(m_Nodes, m_LNodes, m_Grid);
 382     m_LNodesStamp = m_NodesStamp;
 383   }
 384 }
 385
 386 inline void MeshPartition::checkN2N()
 387 {
 388   checkLNodes();
 389   if (m_LNodesStamp > m_N2NStamp) {
 390     createNodeToNode(m_Cells, m_Nodes, m_LNodes, m_N2N, m_Grid);
 391     m_N2NStamp = m_LNodesStamp;
 392   }
 393 }
 394
 395 inline void MeshPartition::checkN2C()
 396 {
 397   checkLNodes();
 398   if (m_LNodesStamp > m_N2CStamp) {
 399     createNodeToCell(m_Cells, m_Nodes, m_LNodes, m_N2C, m_Grid);
 400     m_N2CStamp = m_LNodesStamp;
 401   }
 402 }
 403
 404 inline void MeshPartition::checkN2BC()
 405 {
 406   checkN2C();
 407   if (m_N2CStamp > m_N2BCStamp) {
 408     createNodeToBC();
 409     m_N2BCStamp = m_N2CStamp;
 410   }
 411 }
 412
 413 inline void MeshPartition::checkC2C()
 414 {
 415   checkLCells();
 416   if (m_CellsStamp > m_C2CStamp) {
 417     createCellToCell(m_Cells, m_C2C, m_Grid);
 418     m_C2CStamp = m_CellsStamp;
 419   }
 420 }
 421
 422 inline const QVector<vtkIdType>& MeshPartition::getCells() const
 423 {
 424   return m_Cells;
 425 }
 426
 427 inline const QVector<int>& MeshPartition::getLocalCells()
 428 {
 429   checkLCells();
 430   return m_LCells;
 431 }
 432
 433 inline const QVector<vtkIdType>& MeshPartition::getNodes()
 434 {
 435   checkNodes();
 436   return m_Nodes;
 437 }
 438
 439 inline const QVector<int>& MeshPartition::getLocalNodes()
 440 {
 441   checkLNodes();
 442   return m_LNodes;
 443 }
 444
 445 inline const QVector<QVector<int> >& MeshPartition::getN2N()
 446 {
 447   checkN2N();
 448   return m_N2N;
 449 }
 450
 451 inline const QVector<QVector<int> >& MeshPartition::getN2C()
 452 {
 453   checkN2C();
 454   return m_N2C;
 455 }
 456
 457 inline const QVector<QVector<int> >& MeshPartition::getC2C()
 458 {
 459   checkC2C();
 460   return m_C2C;
 461 }
 462
 463 inline int MeshPartition::n2nLSize(int i_nodes)
 464 {
 465   checkN2N();
 466   return m_N2N[i_nodes].size();
 467 }
 468
 469 inline int MeshPartition::n2nLL(int i_nodes, int j)
 470 {
 471   checkN2N();
 472   return m_N2N[i_nodes][j];
 473 }
 474
 475 inline vtkIdType MeshPartition::n2nLG(int i_nodes, int j)
 476 {
 477   checkN2N();
 478   return m_Nodes[m_N2N[i_nodes][j]];
 479 }
 480
 481 inline int MeshPartition::n2nGSize(vtkIdType id_node)
 482 {
 483   checkN2N();
 484   return m_N2N[m_LNodes[id_node]].size();
 485 }
 486
 487 inline int MeshPartition::n2nGL(vtkIdType id_node, int j)
 488 {
 489   checkN2N();
 490   return m_N2N[m_LNodes[id_node]][j];
 491 }
 492
 493 inline vtkIdType MeshPartition::n2nGG(vtkIdType id_node, int j)
 494 {
 495   checkN2N();
 496   return m_Nodes[m_N2N[m_LNodes[id_node]][j]];
 497 }
 498
 499 inline int MeshPartition::n2cLSize(int i_nodes)
 500 {
 501   checkN2C();
 502   return m_N2C[i_nodes].size();
 503 }
 504
 505 inline int MeshPartition::n2cLL(int i_nodes, int j)
 506 {
 507   checkN2C();
 508   return m_N2C[i_nodes][j];
 509 }
 510
 511 inline vtkIdType MeshPartition::n2cLG(int i_nodes, int j)
 512 {
 513   checkN2C();
 514   int i_cell = m_N2C[i_nodes][j];
 515   if(i_cell<0) return(-1);
 516   else return m_Cells[i_cell];
 517 }
 518
 519 inline int MeshPartition::n2cGSize(vtkIdType id_node)
 520 {
 521   checkN2C();
 522   return m_N2C[m_LNodes[id_node]].size();
 523 }
 524
 525 inline int MeshPartition::n2cGL(vtkIdType id_node, int j)
 526 {
 527   checkN2C();
 528   return m_N2C[m_LNodes[id_node]][j];
 529 }
 530
 531 inline vtkIdType MeshPartition::n2cGG(vtkIdType id_node, int j)
 532 {
 533   checkN2C();
 534   int i_cell = m_N2C[m_LNodes[id_node]][j];
 535   if(i_cell<0) return(-1);
 536   else return m_Cells[i_cell];
 537 }
 538
 539 inline int MeshPartition::c2cLSize(int i_cells)
 540 {
 541   checkC2C();
 542   return m_C2C[i_cells].size();
 543 }
 544
 545 inline int MeshPartition::c2cLL(int i_cells, int j)
 546 {
 547   checkC2C();
 548   return m_C2C[i_cells][j];
 549 }
 550
 551 inline vtkIdType MeshPartition::c2cLG(int i_cells, int j)
 552 {
 553   checkC2C();
 554   int i_cell = m_C2C[i_cells][j];
 555   if(i_cell<0) return(-1);
 556   else return m_Cells[i_cell];
 557 }
 558
 559 inline int MeshPartition::c2cGSize(vtkIdType id_cell)
 560 {
 561   checkC2C();
 562   checkLCells();
 563   return m_C2C[m_LCells[id_cell]].size();
 564 }
 565
 566 inline int MeshPartition::c2cGL(vtkIdType id_cell, int j)
 567 {
 568   checkC2C();
 569   checkLCells();
 570   return m_C2C[m_LCells[id_cell]][j];
 571 }
 572
 573 inline vtkIdType MeshPartition::c2cGG(vtkIdType id_cell, int j)
 574 {
 575   checkC2C();
 576   checkLCells();
 577   int i_cell = m_C2C[m_LCells[id_cell]][j];
 578   if(i_cell<0) return(-1);
 579   else return m_Cells[i_cell];
 580 }
 581
 582 inline int MeshPartition::getNumberOfCells()
 583 {
 584   return m_Cells.size();
 585 }
 586
 587 inline int MeshPartition::getNumberOfNodes()
 588 {
 589   checkNodes();
 590   return m_Nodes.size();
 591 }
 592
 593 inline int MeshPartition::localNode(vtkIdType id_node)
 594 {
 595   checkLNodes();
 596   return m_LNodes[id_node];
 597 }
 598
 599 inline vtkIdType MeshPartition::globalNode(int i)
 600 {
 601   checkNodes();
 602   return m_Nodes[i];
 603 }
 604
 605 inline int MeshPartition::localCell(vtkIdType id_cell)
 606 {
 607   checkLCells();
 608   return m_LCells[id_cell];
 609 }
 610
 611 inline vtkIdType MeshPartition::globalCell(int i)
 612 {
 613   if(i<0) return(-1);
 614   else return m_Cells[i];
 615 }
 616
 617 inline int MeshPartition::n2bcLSize(int i_nodes)
 618 {
 619   checkN2BC();
 620   return m_N2BC[i_nodes].size();
 621 }
 622
 623 inline int MeshPartition::n2bcL(int i_nodes, int j)
 624 {
 625   checkN2BC();
 626   return m_N2BC[i_nodes][j];
 627 }
 628
 629 inline int MeshPartition::n2bcGSize(vtkIdType id_node)
 630 {
 631   checkN2BC();
 632   return m_N2BC[m_LNodes[id_node]].size();
 633 }
 634
 635 inline int MeshPartition::n2bcG(vtkIdType id_node, int j)
 636 {
 637   checkN2BC();
 638   return m_N2BC[m_LNodes[id_node]][j];
 639 }
 640
 641 template <typename C>
 642 void MeshPartition::getGlobalN2N(vtkIdType id_node, C& cont)
 643 {
 644   cont.clear();
 645   for (int i = 0; i < n2nGSize(id_node); ++i) {
 646     cont << n2nGG(id_node, i);
 647   }
 648 }
 649
 650 template <typename C>
 651 void MeshPartition::getEdgeFaces(vtkIdType id_node1, vtkIdType id_node2, C &edge_faces)
 652 {
 653   edge_faces.clear();
 654   for (int i = 0; i < n2cGSize(id_node1); ++i) {
 655     vtkIdType id_cell = n2cGG(id_node1, i);
 656     if (isSurface(id_cell, m_Grid)) {
 657       vtkIdType num_pts, *pts;
 658       m_Grid->GetCellPoints(id_cell, num_pts, pts);
 659       for (int j = 0; j < num_pts; ++j) {
 660         if (pts[j] == id_node2) {
 661           edge_faces << id_cell;
 662           break;
 663         }
 664       }
 665     }
 666   }
 667 }
 668
 669 template <typename C>
 670 void MeshPartition::getCommonNodes(vtkIdType id_cell1, vtkIdType id_cell2, C &common_nodes)
 671 {
 672   common_nodes.clear();
 673   QSet<vtkIdType> nodes1, nodes2;
 674   vtkIdType num_pts, *pts;
 675   m_Grid->GetCellPoints(id_cell1, num_pts, pts);
 676   for (int i = 0; i < num_pts; ++i) {
 677     nodes1.insert(pts[i]);
 678   }
 679   m_Grid->GetCellPoints(id_cell2, num_pts, pts);
 680   for (int i = 0; i < num_pts; ++i) {
 681     nodes2.insert(pts[i]);
 682   }
 683   nodes1.intersect(nodes2);
 684   common_nodes.resize(nodes1.size());
 685   qCopy(nodes1.begin(), nodes1.end(), common_nodes.begin());
 686 }
 687
 688 #endif // MESHPARTITION_H