src/libengrid/meshpartition.cpp

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