src/libengrid/meshpartition.cpp

   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 #include "meshpartition.h"
  23 #include "guimainwindow.h"
  24
  25 #include <vtkKdTreePointLocator.h>
  26 #include <vtkClipDataSet.h>
  27
  28 MeshPartition::MeshPartition()
  29 {
  30   m_Grid = NULL;
  31   m_TrackGrid = false;
  32   resetTimeStamps();
  33 }
  34
  35 MeshPartition::MeshPartition(vtkUnstructuredGrid *grid, bool use_all_cells)
  36 {
  37   m_TrackGrid = false;
  38   resetTimeStamps();
  39   m_Grid = grid;
  40   if (use_all_cells) {
  41     QVector<vtkIdType> cls(grid->GetNumberOfCells());
  42     for (vtkIdType id_cell = 0; id_cell < cls.size(); ++id_cell) {
  43       cls[id_cell] = id_cell;
  44     }
  45     setCells(cls);
  46   }
  47 }
  48
  49 MeshPartition::MeshPartition(QString volume_name)
  50 {
  51   m_TrackGrid = false;
  52   resetTimeStamps();
  53   setVolume(volume_name);
  54 }
  55
  56
  57 void MeshPartition::resetTimeStamps()
  58 {
  59   m_CellsStamp = 0;
  60   m_LCellsStamp = 0;
  61   m_NodesStamp = 0;
  62   m_LNodesStamp = 0;
  63   m_N2NStamp = 0;
  64   m_N2CStamp = 0;
  65   m_N2BCStamp = 0;
  66   m_C2CStamp = 0;
  67 }
  68
  69 void MeshPartition::trackGrid(vtkUnstructuredGrid *grid)
  70 {
  71   setGrid(grid);
  72   setAllCells();
  73   m_GridMTime = m_Grid->GetMTime();
  74   m_TrackGrid = true;
  75 }
  76
  77 void MeshPartition::setVolume(QString volume_name)
  78 {
  79   m_Grid = GuiMainWindow::pointer()->getGrid();
  80   resetOrientation(m_Grid);
  81   VolumeDefinition V = GuiMainWindow::pointer()->getVol(volume_name);
  82   QList<vtkIdType> cls;
  83   EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
  84   EG_VTKDCC(vtkIntArray, cell_orgdir, m_Grid, "cell_orgdir");
  85   EG_VTKDCC(vtkIntArray, cell_curdir, m_Grid, "cell_curdir");
  86   EG_VTKDCC(vtkIntArray, cell_voldir, m_Grid, "cell_voldir");
  87   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  88     if (isSurface(id_cell, m_Grid)) {
  89       int bc = cell_code->GetValue(id_cell);
  90       cell_voldir->SetValue(id_cell, 0);
  91       if (V.getSign(bc) != 0) {
  92         cls.append(id_cell);
  93         if (V.getSign(bc) == -1) {
  94           cell_voldir->SetValue(id_cell, 1);
  95         }
  96       }
  97     } else {
  98       if (cell_code->GetValue(id_cell) == V.getVC()) {
  99         cls.append(id_cell);
 100       }
 101     }
 102   }
 103   setCells(cls);
 104 }
 105
 106 void MeshPartition::setVolumeOrientation()
 107 {
 108   EG_VTKDCC(vtkIntArray, cell_curdir, m_Grid, "cell_curdir");
 109   EG_VTKDCC(vtkIntArray, cell_voldir, m_Grid, "cell_voldir");
 110   foreach (vtkIdType id_cell, m_Cells) {
 111     if (isSurface(id_cell, m_Grid)) {
 112       if (cell_curdir->GetValue(id_cell) != cell_voldir->GetValue(id_cell)) {
 113         reorientateFace(m_Grid, id_cell);
 114       }
 115     }
 116   }
 117 }
 118
 119 void MeshPartition::setOriginalOrientation()
 120 {
 121   EG_VTKDCC(vtkIntArray, cell_curdir, m_Grid, "cell_curdir");
 122   EG_VTKDCC(vtkIntArray, cell_orgdir, m_Grid, "cell_orgdir");
 123   foreach (vtkIdType id_cell, m_Cells) {
 124     if (isSurface(id_cell, m_Grid)) {
 125       if (cell_curdir->GetValue(id_cell) != cell_orgdir->GetValue(id_cell)) {
 126         reorientateFace(m_Grid, id_cell);
 127       }
 128     }
 129   }
 130 }
 131
 132 void MeshPartition::setRemainder(const MeshPartition& part)
 133 {
 134   setGrid(part.getGrid());
 135   QVector<vtkIdType> rcells;
 136   getRestCells(m_Grid, part.m_Cells, rcells);
 137   setCells(rcells);
 138 }
 139
 140 void MeshPartition::extractToVtkGrid(vtkUnstructuredGrid *new_grid)
 141 {
 142   checkLNodes();
 143   allocateGrid(new_grid, m_Cells.size(), m_Nodes.size());
 144   for (int i_nodes = 0; i_nodes < m_Nodes.size(); ++i_nodes) {
 145     vec3_t x;
 146     m_Grid->GetPoints()->GetPoint(m_Nodes[i_nodes], x.data());
 147     new_grid->GetPoints()->SetPoint(i_nodes, x.data());
 148     copyNodeData(m_Grid, m_Nodes[i_nodes], new_grid, i_nodes);
 149   }
 150   foreach (vtkIdType id_cell, m_Cells) {
 151     /*
 152     vtkIdType N_pts, *pts;
 153     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 154     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 155     QVector<vtkIdType> new_pts(N_pts);
 156     for (int i = 0; i < N_pts; ++i) {
 157       new_pts[i] = m_LNodes[pts[i]];
 158     }
 159     // update for polyhedral cells here
 160     vtkIdType id_new_cell = new_grid->InsertNextCell(type_cell, N_pts, new_pts.data());
 161     */
 162     vtkIdType id_new_cell = copyCell(m_Grid, id_cell, new_grid, m_LNodes);
 163     copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 164   }
 165 }
 166
 167 void MeshPartition::addPartition(const MeshPartition& part, double tol)
 168 {
 169   if (m_Grid == part.m_Grid) {
 170     QVector<bool> cell_marked(m_Grid->GetNumberOfCells(), false);
 171     foreach (vtkIdType id_cell, m_Cells) {
 172       cell_marked[id_cell] = true;
 173     }
 174     foreach (vtkIdType id_cell, part.m_Cells) {
 175       cell_marked[id_cell] = true;
 176     }
 177     QList<vtkIdType> new_cells;
 178     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 179       if (cell_marked[id_cell]) {
 180         new_cells.append(id_cell);
 181       }
 182     }
 183     setCells(new_cells);
 184   } else {
 185     if (tol < 0) {
 186       tol *= -min(getSmallestEdgeLength(), part.getSmallestEdgeLength());
 187     }
 188     tol = max(tol, 1e-30);
 189     cout << "  tol=" << tol << endl;
 190     EG_VTKSP(vtkUnstructuredGrid, new_grid);
 191     EG_VTKSP(vtkKdTreePointLocator,loc);
 192     loc->SetDataSet(m_Grid);
 193     loc->BuildLocator();
 194     QVector<vtkIdType> pnode2node(part.m_Grid->GetNumberOfPoints());
 195     vtkIdType N = m_Grid->GetNumberOfPoints();
 196     Timer T(10);
 197     for (vtkIdType id_pnode = 0; id_pnode < part.m_Grid->GetNumberOfPoints(); ++id_pnode) {
 198       vec3_t xp, x;
 199       part.m_Grid->GetPoint(id_pnode, xp.data());
 200       vtkIdType id_node = loc->FindClosestPoint(xp.data());
 201       m_Grid->GetPoint(id_node, x.data());
 202       if ((x - xp).abs() < tol) {
 203         pnode2node[id_pnode] = id_node;
 204       } else {
 205         pnode2node[id_pnode] = N;
 206         ++N;
 207       }
 208     }
 209     allocateGrid(new_grid, m_Grid->GetNumberOfCells() + part.m_Cells.size(), N);
 210     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 211       vec3_t x;
 212       m_Grid->GetPoint(id_node, x.data());
 213       new_grid->GetPoints()->SetPoint(id_node, x.data());
 214       copyNodeData(m_Grid, id_node, new_grid, id_node);
 215     }
 216     QVector<vtkIdType> part_nodes;
 217     getNodesFromCells(part.m_Cells, part_nodes, part.m_Grid);
 218     foreach (vtkIdType id_pnode, part_nodes) {
 219       vec3_t x;
 220       part.m_Grid->GetPoint(id_pnode, x.data());
 221       new_grid->GetPoints()->SetPoint(pnode2node[id_pnode], x.data());
 222       copyNodeData(part.m_Grid, id_pnode, new_grid, pnode2node[id_pnode]);
 223     }
 224     QList<vtkIdType> new_cells;
 225     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 226       vtkIdType id_new_cell = copyCell(m_Grid, id_cell, new_grid);
 227       copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 228       new_cells.append(id_new_cell);
 229     }
 230     foreach (vtkIdType id_pcell, part.m_Cells) {
 231       vtkIdType id_new_cell = copyCell(part.m_Grid, id_pcell, new_grid, pnode2node);
 232       copyCellData(part.m_Grid, id_pcell, new_grid, id_new_cell);
 233       new_cells.append(id_new_cell);
 234     }
 235     makeCopy(new_grid, m_Grid);
 236   }
 237 }
 238
 239 double MeshPartition::getSmallestEdgeLength() const
 240 {
 241   double L = 1e99;
 242   foreach (vtkIdType id_cell, m_Cells) {
 243     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 244     vtkIdType N_pts, *pts;
 245     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 246     QVector<vec3_t> x(N_pts);
 247     for (int i = 0; i < N_pts; ++i) {
 248       m_Grid->GetPoint(pts[i], x[i].data());
 249     }
 250     if        (type_cell == VTK_TRIANGLE) {
 251       L = min(L, (x[0] - x[1]).abs());
 252       L = min(L, (x[1] - x[2]).abs());
 253       L = min(L, (x[2] - x[0]).abs());
 254     } else if (type_cell == VTK_QUAD) {
 255       L = min(L, (x[0] - x[1]).abs());
 256       L = min(L, (x[1] - x[2]).abs());
 257       L = min(L, (x[2] - x[3]).abs());
 258       L = min(L, (x[3] - x[0]).abs());
 259     } else if (type_cell == VTK_TETRA) {
 260       L = min(L, (x[0] - x[1]).abs());
 261       L = min(L, (x[1] - x[2]).abs());
 262       L = min(L, (x[2] - x[0]).abs());
 263       L = min(L, (x[3] - x[0]).abs());
 264       L = min(L, (x[3] - x[1]).abs());
 265       L = min(L, (x[3] - x[2]).abs());
 266     } else if (type_cell == VTK_PYRAMID) {
 267       L = min(L, (x[0] - x[1]).abs());
 268       L = min(L, (x[1] - x[2]).abs());
 269       L = min(L, (x[2] - x[3]).abs());
 270       L = min(L, (x[3] - x[0]).abs());
 271       L = min(L, (x[4] - x[0]).abs());
 272       L = min(L, (x[4] - x[1]).abs());
 273       L = min(L, (x[4] - x[2]).abs());
 274       L = min(L, (x[4] - x[3]).abs());
 275     } else if (type_cell == VTK_WEDGE) {
 276       L = min(L, (x[0] - x[1]).abs());
 277       L = min(L, (x[1] - x[2]).abs());
 278       L = min(L, (x[2] - x[0]).abs());
 279       L = min(L, (x[3] - x[4]).abs());
 280       L = min(L, (x[4] - x[5]).abs());
 281       L = min(L, (x[5] - x[3]).abs());
 282       L = min(L, (x[0] - x[3]).abs());
 283       L = min(L, (x[1] - x[4]).abs());
 284       L = min(L, (x[2] - x[5]).abs());
 285     } else if (type_cell == VTK_HEXAHEDRON) {
 286       L = min(L, (x[0] - x[1]).abs());
 287       L = min(L, (x[1] - x[2]).abs());
 288       L = min(L, (x[2] - x[3]).abs());
 289       L = min(L, (x[3] - x[0]).abs());
 290       L = min(L, (x[4] - x[5]).abs());
 291       L = min(L, (x[5] - x[6]).abs());
 292       L = min(L, (x[6] - x[7]).abs());
 293       L = min(L, (x[7] - x[4]).abs());
 294       L = min(L, (x[0] - x[4]).abs());
 295       L = min(L, (x[1] - x[5]).abs());
 296       L = min(L, (x[2] - x[6]).abs());
 297       L = min(L, (x[3] - x[7]).abs());
 298     }
 299   }
 300   return L;
 301 }
 302
 303 bool MeshPartition::hasNeighNode(vtkIdType id_node, vtkIdType id_neigh)
 304 {
 305   for (int i = 0; i < n2nGSize(id_node); ++i) {
 306     if (n2nGG(id_node, i) == id_neigh) {
 307       return true;
 308     }
 309   }
 310   return false;
 311 }
 312
 313 void MeshPartition::createNodeToBC()
 314 {
 315   EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
 316   m_N2BC.resize(m_Nodes.size());
 317   for (int i_node = 0; i_node < m_Nodes.size(); ++i_node) {
 318     QSet<int> bcs;
 319     for (int j = 0; j < n2cLSize(i_node); ++j) {
 320       vtkIdType id_cell = n2cLG(i_node, j);
 321       if (isSurface(id_cell, m_Grid)) {
 322         int bc = cell_code->GetValue(n2cLG(i_node, j));
 323         if (bc != 0) {
 324           bcs.insert(cell_code->GetValue(n2cLG(i_node, j)));
 325         }
 326       }
 327     }
 328     m_N2BC[i_node].resize(bcs.size());
 329     int i = 0;
 330     foreach (int bc, bcs) {
 331       m_N2BC[i_node][i++] = bc;
 332     }
 333   }
 334 }
 335
 336 bool MeshPartition::hasBC(vtkIdType id_node, int bc)
 337 {
 338   bool found = false;
 339   for (int j = 0; j < n2bcGSize(id_node); ++j) {
 340     if (n2bcG(id_node, j) == bc) {
 341       found = true;
 342       break;
 343     }
 344   }
 345   return found;
 346 }
 347
 348 vtkIdType MeshPartition::getVolumeCell(vtkIdType id_face)
 349 {
 350   checkLNodes();
 351   checkLCells();
 352   checkN2C();
 353   return findVolumeCell(m_Grid, id_face, m_LNodes, m_Cells, m_LCells, m_N2C);
 354 }
 355
 356 vec3_t MeshPartition::globalNormal(vtkIdType id_node)
 357 {
 358   vec3_t normal(0,0,0);
 359   for (int i = 0; i < n2cGSize(id_node); ++i) {
 360     vtkIdType id_cell = n2cGG(id_node, i);
 361     if (isSurface(id_cell, m_Grid)) {
 362       vtkIdType N_pts, *pts;
 363       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 364       vec3_t a, b, c;
 365       for (int j = 0; j < N_pts; ++j) {
 366         if (pts[j] == id_node) {
 367           m_Grid->GetPoint(pts[j], a.data());
 368           if (j > 0) {
 369             m_Grid->GetPoint(pts[j-1], b.data());
 370           } else {
 371             m_Grid->GetPoint(pts[N_pts-1], b.data());
 372           }
 373           if (j < N_pts - 1) {
 374             m_Grid->GetPoint(pts[j+1], c.data());
 375           } else {
 376             m_Grid->GetPoint(pts[0], c.data());
 377           }
 378         }
 379       }
 380       vec3_t u = b - a;
 381       vec3_t v = c - a;
 382       double alpha = GeometryTools::angle(u, v);
 383       vec3_t n = u.cross(v);
 384       n.normalise();
 385       if (checkVector(n)) {
 386         normal -= alpha*n;
 387       }
 388     }
 389   }
 390   normal.normalise();
 391   return normal;
 392 }
 393
 394 double MeshPartition::getAverageSurfaceEdgeLength(vtkIdType id_node)
 395 {
 396   QSet<vtkIdType> surface_neighbours;
 397   for (int i = 0; i < n2cGSize(id_node); ++i) {
 398     vtkIdType id_cell = n2cGG(id_node, i);
 399     if (isSurface(id_cell, m_Grid)) {
 400       EG_GET_CELL(id_cell, m_Grid);
 401       for (int j = 0; j < num_pts; ++j) {
 402         if (pts[j] != id_node) {
 403           surface_neighbours.insert(pts[j]);
 404         }
 405       }
 406     }
 407   }
 408   double L = 0;
 409   if (surface_neighbours.size() > 0) {
 410     vec3_t x, xn;
 411     m_Grid->GetPoint(id_node, x.data());
 412     foreach (vtkIdType id_neigh, surface_neighbours) {
 413       m_Grid->GetPoint(id_neigh, xn.data());
 414       L += (x - xn).abs();
 415     }
 416     L /= surface_neighbours.size();
 417   }
 418   return L;
 419 }
 420
 421 void MeshPartition::computeMinAndMaxSurfaceStencilEdgeLengths(vtkIdType id_node, double &l_min, double &l_max)
 422 {
 423   l_min = 1e99;
 424   l_max = 0;
 425   for (int i = 0; i < n2cGSize(id_node); ++i) {
 426     vtkIdType id_cell = n2cGG(id_node, i);
 427     if (isSurface(id_cell, m_Grid)) {
 428       vtkIdType *pts, num_pts;
 429       m_Grid->GetCellPoints(id_cell, num_pts, pts);
 430       vec3_t x1, x2;
 431       m_Grid->GetPoint(pts[num_pts - 1], x1.data());
 432       for (int j = 0; j < num_pts; ++j) {
 433         m_Grid->GetPoint(pts[j], x2.data());
 434         double L = (x1 - x2).abs();
 435         l_min = min(L, l_min);
 436         l_max = max(L, l_max);
 437         x1 = x2;
 438       }
 439     }
 440   }
 441 }
 442
 443 double MeshPartition::getMinSurfaceStencilEdgeLength(vtkIdType id_node)
 444 {
 445   double l_min, l_max;
 446   computeMinAndMaxSurfaceStencilEdgeLengths(id_node, l_min, l_max);
 447   return l_min;
 448 }
 449
 450 double MeshPartition::getMaxSurfaceStencilEdgeLength(vtkIdType id_node)
 451 {
 452   double l_min, l_max;
 453   computeMinAndMaxSurfaceStencilEdgeLengths(id_node, l_min, l_max);
 454   return l_max;
 455 }
 456
 457 int MeshPartition::getNumberOfFeatureNeighbours(vtkIdType id_node)
 458 {
 459   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 460   int N = 0;
 461   for (int i = 0; i < n2nGSize(id_node); ++i) {
 462     char type = node_type->GetValue(n2nGG(id_node, i));
 463     if (type == EG_FEATURE_EDGE_VERTEX || type == EG_FEATURE_CORNER_VERTEX) {
 464       ++N;
 465     }
 466   }
 467   return N;
 468 }
 469
 470 int MeshPartition::getEdgeType(vtkIdType id_node1, vtkIdType id_node2)
 471 {
 472   QList <vtkIdType> edge_faces;
 473   getEdgeFaces(id_node1, id_node2, edge_faces);
 474   int edge = 0;
 475   if (edge_faces.size() == 1) {
 476     edge = 2;
 477   } else if (edge_faces.size() >= 2) {
 478     EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 479     QSet<int> bcs;
 480     foreach (vtkIdType id_face, edge_faces) {
 481       bcs.insert(cell_code->GetValue(id_face));
 482     }
 483     if (bcs.size() > 1) {
 484       edge = 2;
 485     } else {
 486       edge = 1;
 487     }
 488   }
 489   return edge;
 490 }
 491
 492 int MeshPartition::computeTopoDistance(vtkIdType id_node1, vtkIdType id_node2, int max_dist, int restriction_type)
 493 {
 494   int dist = 0;
 495   QSet<vtkIdType> candidates;
 496   candidates.insert(id_node1);
 497   while (dist < max_dist && !candidates.contains(id_node2)) {
 498     foreach (vtkIdType id_node, candidates) {
 499       for (int i = 0; i < n2nGSize(id_node); ++i) {
 500         vtkIdType id_neigh = n2nGG(id_node,i);
 501         if (!candidates.contains(id_neigh)) {
 502           if (getEdgeType(id_node, id_neigh) >= restriction_type) {
 503             candidates.insert(id_neigh);
 504           }
 505         }
 506       }
 507     }
 508     ++dist;
 509   }
 510   return max_dist;
 511 }