src/meshpartition.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
  24 #include "meshpartition.h"
  25 #include "guimainwindow.h"
  26
  27 #include <vtkKdTreePointLocator.h>
  28
  29 MeshPartition::MeshPartition()
  30 {
  31   m_Grid = NULL;
  32   resetTimeStamps();
  33 }
  34
  35 MeshPartition::MeshPartition(vtkUnstructuredGrid *grid, bool use_all_cells)
  36 {
  37   resetTimeStamps();
  38   m_Grid = grid;
  39   if (use_all_cells) {
  40     QVector<vtkIdType> cls(grid->GetNumberOfCells());
  41     for (vtkIdType id_cell = 0; id_cell < cls.size(); ++id_cell) {
  42       cls[id_cell] = id_cell;
  43     }
  44     setCells(cls);
  45   }
  46 }
  47
  48 MeshPartition::MeshPartition(QString volume_name)
  49 {
  50   resetTimeStamps();
  51   setVolume(volume_name);
  52 }
  53
  54 void MeshPartition::resetTimeStamps()
  55 {
  56   m_CellsStamp = 0;
  57   m_LCellsStamp = 0;
  58   m_NodesStamp = 0;
  59   m_LNodesStamp = 0;
  60   m_N2NStamp = 0;
  61   m_N2CStamp = 0;
  62   m_C2CStamp = 0;
  63 }
  64
  65 void MeshPartition::setVolume(QString volume_name)
  66 {
  67   m_Grid = GuiMainWindow::pointer()->getGrid();
  68   resetOrientation(m_Grid);
  69   VolumeDefinition V = GuiMainWindow::pointer()->getVol(volume_name);
  70   QList<vtkIdType> cls;
  71   EG_VTKDCC(vtkIntArray, cell_code,   m_Grid, "cell_code");
  72   EG_VTKDCC(vtkIntArray, cell_orgdir, m_Grid, "cell_orgdir");
  73   EG_VTKDCC(vtkIntArray, cell_curdir, m_Grid, "cell_curdir");
  74   EG_VTKDCC(vtkIntArray, cell_voldir, m_Grid, "cell_voldir");
  75   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  76     if (isSurface(id_cell, m_Grid)) {
  77       int bc = cell_code->GetValue(id_cell);
  78       cell_voldir->SetValue(id_cell, 0);
  79       if (V.getSign(bc) != 0) {
  80         cls.append(id_cell);
  81         if (V.getSign(bc) == -1) {
  82           cell_voldir->SetValue(id_cell, 1);
  83         }
  84       }
  85     } else {
  86       if (cell_code->GetValue(id_cell) == V.getVC()) {
  87         cls.append(id_cell);
  88       }
  89     }
  90   }
  91   setCells(cls);
  92 }
  93
  94 void MeshPartition::setVolumeOrientation()
  95 {
  96   EG_VTKDCC(vtkIntArray, cell_curdir, m_Grid, "cell_curdir");
  97   EG_VTKDCC(vtkIntArray, cell_voldir, m_Grid, "cell_voldir");
  98   foreach (vtkIdType id_cell, m_Cells) {
  99     if (isSurface(id_cell, m_Grid)) {
 100       if (cell_curdir->GetValue(id_cell) != cell_voldir->GetValue(id_cell)) {
 101         reorientateFace(m_Grid, id_cell);
 102       }
 103     }
 104   }
 105 }
 106
 107 void MeshPartition::setOriginalOrientation()
 108 {
 109   EG_VTKDCC(vtkIntArray, cell_curdir, m_Grid, "cell_curdir");
 110   EG_VTKDCC(vtkIntArray, cell_orgdir, m_Grid, "cell_orgdir");
 111   foreach (vtkIdType id_cell, m_Cells) {
 112     if (isSurface(id_cell, m_Grid)) {
 113       if (cell_curdir->GetValue(id_cell) != cell_orgdir->GetValue(id_cell)) {
 114         reorientateFace(m_Grid, id_cell);
 115       }
 116     }
 117   }
 118 }
 119
 120 void MeshPartition::setRemainder(const MeshPartition& part)
 121 {
 122   setGrid(part.getGrid());
 123   QVector<vtkIdType> rcells;
 124   getRestCells(m_Grid, part.m_Cells, rcells);
 125   setCells(rcells);
 126 }
 127
 128 void MeshPartition::extractToVtkGrid(vtkUnstructuredGrid *new_grid)
 129 {
 130   checkLNodes();
 131   allocateGrid(new_grid, m_Cells.size(), m_Nodes.size());
 132   for (int i_nodes = 0; i_nodes < m_Nodes.size(); ++i_nodes) {
 133     vec3_t x;
 134     m_Grid->GetPoints()->GetPoint(m_Nodes[i_nodes], x.data());
 135     new_grid->GetPoints()->SetPoint(i_nodes, x.data());
 136     copyNodeData(m_Grid, m_Nodes[i_nodes], new_grid, i_nodes);
 137   }
 138   foreach (vtkIdType id_cell, m_Cells) {
 139     vtkIdType N_pts, *pts;
 140     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 141     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 142     vtkIdType new_pts[N_pts];
 143     for (int i = 0; i < N_pts; ++i) {
 144       new_pts[i] = m_LNodes[pts[i]];
 145     }
 146     vtkIdType id_new_cell = new_grid->InsertNextCell(type_cell, N_pts, new_pts);
 147     copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 148   }
 149 }
 150
 151 void MeshPartition::addPartition(const MeshPartition& part)
 152 {
 153   if (m_Grid == part.m_Grid) {
 154     QVector<bool> cell_marked(m_Grid->GetNumberOfCells(), false);
 155     foreach (vtkIdType id_cell, m_Cells) {
 156       cell_marked[id_cell] = true;
 157     }
 158     foreach (vtkIdType id_cell, part.m_Cells) {
 159       cell_marked[id_cell] = true;
 160     }
 161     QList<vtkIdType> new_cells;
 162     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 163       if (cell_marked[id_cell]) {
 164         new_cells.append(id_cell);
 165       }
 166     }
 167     setCells(new_cells);
 168   } else {
 169     double tol = 1e-3*min(getSmallestEdgeLength(), part.getSmallestEdgeLength());
 170     EG_VTKSP(vtkUnstructuredGrid, new_grid);
 171     EG_VTKSP(vtkKdTreePointLocator,loc);
 172     loc->SetDataSet(m_Grid);
 173     loc->BuildLocator();
 174     QVector<vtkIdType> pnode2node(part.m_Grid->GetNumberOfPoints());
 175     vtkIdType N = m_Grid->GetNumberOfPoints();
 176     for (vtkIdType id_pnode = 0; id_pnode < part.m_Grid->GetNumberOfPoints(); ++id_pnode) {
 177       vec3_t xp, x;
 178       part.m_Grid->GetPoint(id_pnode, xp.data());
 179       vtkIdType id_node = loc->FindClosestPoint(xp.data());
 180       m_Grid->GetPoint(id_node, x.data());
 181       if ((x - xp).abs() < tol) {
 182         pnode2node[id_pnode] = id_node;
 183       } else {
 184         pnode2node[id_pnode] = N;
 185         ++N;
 186       }
 187     }
 188     allocateGrid(new_grid, m_Grid->GetNumberOfCells() + part.m_Cells.size(), N);
 189     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 190       vec3_t x;
 191       m_Grid->GetPoint(id_node, x.data());
 192       new_grid->GetPoints()->SetPoint(id_node, x.data());
 193       copyNodeData(m_Grid, id_node, new_grid, id_node);
 194     }
 195     QVector<vtkIdType> part_nodes;
 196     getNodesFromCells(part.m_Cells, part_nodes, part.m_Grid);
 197     foreach (vtkIdType id_pnode, part_nodes) {
 198       vec3_t x;
 199       part.m_Grid->GetPoint(id_pnode, x.data());
 200       new_grid->GetPoints()->SetPoint(pnode2node[id_pnode], x.data());
 201       copyNodeData(part.m_Grid, id_pnode, new_grid, pnode2node[id_pnode]);
 202     }
 203     QList<vtkIdType> new_cells;
 204     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 205       vtkIdType N_pts, *pts;
 206       vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 207       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 208       vtkIdType id_new_cell = new_grid->InsertNextCell(type_cell, N_pts, pts);
 209       copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 210       new_cells.append(id_new_cell);
 211     }
 212     foreach (vtkIdType id_pcell, part.m_Cells) {
 213       vtkIdType N_pts, *pts;
 214       vtkIdType type_cell = part.m_Grid->GetCellType(id_pcell);
 215       part.m_Grid->GetCellPoints(id_pcell, N_pts, pts);
 216       vtkIdType new_pts[N_pts];
 217       for (int i = 0; i < N_pts; ++i) {
 218         new_pts[i] = pnode2node[pts[i]];
 219       }
 220       vtkIdType id_new_cell = new_grid->InsertNextCell(type_cell, N_pts, new_pts);
 221       copyCellData(part.m_Grid, id_pcell, new_grid, id_new_cell);
 222       new_cells.append(id_new_cell);
 223     }
 224     makeCopy(new_grid, m_Grid);
 225   }
 226 }
 227
 228 double MeshPartition::getSmallestEdgeLength() const
 229 {
 230   double L = 1e99;
 231   foreach (vtkIdType id_cell, m_Cells) {
 232     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 233     vtkIdType N_pts, *pts;
 234     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 235     vec3_t x[N_pts];
 236     for (int i = 0; i < N_pts; ++i) {
 237       m_Grid->GetPoint(pts[i], x[i].data());
 238     }
 239     if        (type_cell == VTK_TRIANGLE) {
 240       L = min(L, (x[0] - x[1]).abs());
 241       L = min(L, (x[1] - x[2]).abs());
 242       L = min(L, (x[2] - x[0]).abs());
 243     } else if (type_cell == VTK_QUAD) {
 244       L = min(L, (x[0] - x[1]).abs());
 245       L = min(L, (x[1] - x[2]).abs());
 246       L = min(L, (x[2] - x[3]).abs());
 247       L = min(L, (x[3] - x[0]).abs());
 248     } else if (type_cell == VTK_TETRA) {
 249       L = min(L, (x[0] - x[1]).abs());
 250       L = min(L, (x[1] - x[2]).abs());
 251       L = min(L, (x[2] - x[0]).abs());
 252       L = min(L, (x[3] - x[0]).abs());
 253       L = min(L, (x[3] - x[1]).abs());
 254       L = min(L, (x[3] - x[2]).abs());
 255     } else if (type_cell == VTK_PYRAMID) {
 256       L = min(L, (x[0] - x[1]).abs());
 257       L = min(L, (x[1] - x[2]).abs());
 258       L = min(L, (x[2] - x[3]).abs());
 259       L = min(L, (x[3] - x[0]).abs());
 260       L = min(L, (x[4] - x[0]).abs());
 261       L = min(L, (x[4] - x[1]).abs());
 262       L = min(L, (x[4] - x[2]).abs());
 263       L = min(L, (x[4] - x[3]).abs());
 264     } else if (type_cell == VTK_WEDGE) {
 265       L = min(L, (x[0] - x[1]).abs());
 266       L = min(L, (x[1] - x[2]).abs());
 267       L = min(L, (x[2] - x[0]).abs());
 268       L = min(L, (x[3] - x[4]).abs());
 269       L = min(L, (x[4] - x[5]).abs());
 270       L = min(L, (x[5] - x[3]).abs());
 271       L = min(L, (x[0] - x[3]).abs());
 272       L = min(L, (x[1] - x[4]).abs());
 273       L = min(L, (x[2] - x[5]).abs());
 274     } else if (type_cell == VTK_HEXAHEDRON) {
 275       L = min(L, (x[0] - x[1]).abs());
 276       L = min(L, (x[1] - x[2]).abs());
 277       L = min(L, (x[2] - x[3]).abs());
 278       L = min(L, (x[3] - x[0]).abs());
 279       L = min(L, (x[4] - x[5]).abs());
 280       L = min(L, (x[5] - x[6]).abs());
 281       L = min(L, (x[6] - x[7]).abs());
 282       L = min(L, (x[7] - x[4]).abs());
 283       L = min(L, (x[0] - x[4]).abs());
 284       L = min(L, (x[1] - x[5]).abs());
 285       L = min(L, (x[2] - x[6]).abs());
 286       L = min(L, (x[3] - x[7]).abs());
 287     }
 288   }
 289   return L;
 290 }
 291