src/libengrid/vtkEgNormalExtrusion.cxx

   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 #include "vtkEgNormalExtrusion.h"
  22
  23 vtkStandardNewMacro(vtkEgNormalExtrusion)
  24
  25 vtkEgNormalExtrusion::vtkEgNormalExtrusion()
  26 {
  27   mode = normal;
  28   m_ScaleX = 1;
  29   m_ScaleY = 1;
  30   m_ScaleZ = 1;
  31   m_RemoveInternalFaces = true;
  32 }
  33
  34 void vtkEgNormalExtrusion::ExecuteEg()
  35 {
  36   QVector<vtkIdType> cells, nodes, n1, n2;
  37   QVector<vec3_t> cell_normals, node_normals;
  38   ExtractBoundary(cells, nodes, m_BoundaryCodes, m_Input);
  39   if (mode == normal) {
  40     computeNormals(cell_normals, node_normals, cells, nodes,m_Input);
  41   } else if (mode == cylinder) {
  42     axis.normalise();
  43     node_normals.resize(nodes.size());
  44     for (int i = 0; i < node_normals.size(); ++i) {
  45       vec3_t x;
  46       m_Input->GetPoint(nodes[i],x.data());
  47       vec3_t x0 = x - ((x-origin)*axis)*axis;
  48       node_normals[i] = x0 - origin;
  49       node_normals[i].normalise();
  50     }
  51   } else if ((mode == fixed) || (mode == planar)) {
  52     fixed_normal.normalise();
  53     node_normals.resize(nodes.size());
  54     for (int i = 0; i < node_normals.size(); ++i) {
  55       node_normals[i] = fixed_normal;
  56     }
  57   }
  58   for (int i = 0; i < node_normals.size(); ++ i) {
  59     node_normals[i][0] *= m_ScaleX;
  60     node_normals[i][1] *= m_ScaleY;
  61     node_normals[i][2] *= m_ScaleZ;
  62   }
  63   n1.resize(nodes.size());
  64   n2.resize(nodes.size());
  65
  66   // mapping
  67   QVector<int> _cells, _nodes;
  68   createNodeMapping(nodes, _nodes, m_Input);
  69   createCellMapping(cells, _cells, m_Input);
  70   QVector<QSet<int> > n2c;
  71   createNodeToCell(cells, nodes, _nodes, n2c, m_Input);
  72
  73   vtkIdType NnewNodes = m_Input->GetNumberOfPoints() + (layer_y.size()-1)*nodes.size();
  74   vtkIdType NnewCells = m_Input->GetNumberOfCells() + (layer_y.size()-1)*cells.size();
  75
  76   // count the number of new surface elements (side walls)
  77   for (int i_cell = 0; i_cell < cells.size(); ++i_cell) {
  78     vtkIdType *pts;
  79     vtkIdType  Npts;
  80     m_Input->GetCellPoints(cells[i_cell], Npts, pts);
  81     QVector<vtkIdType> surf_pts(Npts);
  82     for (int i_pts = 0; i_pts < Npts; ++i_pts) {
  83       surf_pts[i_pts] = _nodes[pts[i_pts]];
  84     }
  85     for (int i_surf_pts = 0; i_surf_pts < Npts; ++i_surf_pts) {
  86       vtkIdType p1 = surf_pts[i_surf_pts];
  87       vtkIdType p2;
  88       if (i_surf_pts < Npts - 1) {
  89         p2 = surf_pts[i_surf_pts + 1];
  90       } else {
  91         p2 = surf_pts[0];
  92       }
  93       bool add_bd = false;
  94       {
  95         int N = 0;
  96         int cell;
  97         foreach(cell, n2c[p1]) {
  98           if (n2c[p2].contains(cell)) ++N;
  99         }
 100         if (N == 1) add_bd = true;
 101       }
 102       if (add_bd) {
 103         NnewCells += layer_y.size();
 104       }
 105     }
 106   }
 107
 108   // count the number of new surface elements (base)
 109   QVector<bool> is_boundary;
 110   is_boundary.fill(false, cells.size());
 111   {
 112     QVector<int> nvol;
 113     nvol.fill(0, nodes.size());
 114     for (vtkIdType id_cell = 0; id_cell < m_Input->GetNumberOfCells(); ++id_cell) {
 115       if (isVolume(id_cell, m_Input)) {
 116         vtkIdType Npts, *pts;
 117         m_Input->GetCellPoints(id_cell, Npts, pts);
 118         for (int i = 0; i < Npts; ++i) {
 119           if (_nodes[pts[i]] >= 0) {
 120             ++nvol[_nodes[pts[i]]];
 121           }
 122         }
 123       }
 124     }
 125     for (int i_cell = 0; i_cell < cells.size(); ++i_cell) {
 126       vtkIdType id_cell = cells[i_cell];
 127       vtkIdType Npts, *pts;
 128       m_Input->GetCellPoints(id_cell, Npts, pts);
 129       for (int i = 0; i < Npts; ++i) {
 130         if (nvol[_nodes[pts[i]]] == 0) {
 131           is_boundary[i_cell] = true;
 132         }
 133       }
 134       if (is_boundary[i_cell] || !m_RemoveInternalFaces) {
 135         ++NnewCells;
 136       }
 137     }
 138   }
 139
 140   // allocate memory for the new grid
 141   allocateGrid(m_Output, NnewCells, NnewNodes);
 142
 143   // boundary conditions
 144   EG_VTKDCC(vtkIntArray, cell_code1, m_Input, "cell_code");
 145   EG_VTKDCC(vtkIntArray, cell_code2, m_Output, "cell_code");
 146   EG_VTKDCC(vtkIntArray, orgdir, m_Output, "cell_orgdir");
 147   EG_VTKDCC(vtkIntArray, voldir, m_Output, "cell_voldir");
 148   EG_VTKDCC(vtkIntArray, curdir, m_Output, "cell_curdir");
 149
 150   int new_bc = 1;
 151   for (vtkIdType id_cell = 0; id_cell < m_Input->GetNumberOfCells(); ++id_cell) {
 152     if (isSurface(id_cell, m_Input)) {
 153       if (cell_code1->GetValue(id_cell) >= new_bc) {
 154         new_bc = cell_code1->GetValue(id_cell) + 1;
 155       }
 156     }
 157   }
 158
 159   for (int i = 0; i < nodes.size(); ++i) {
 160     n2[i] = nodes[i];
 161     vec3_t x;
 162     m_Input->GetPoint(nodes[i],x.data());
 163     m_Output->GetPoints()->SetPoint(n2[i],x.data());
 164   }
 165   double total_layers = layer_y[layer_y.size()-1] - layer_y[0];
 166   QVector<double> total_dist(nodes.size());
 167   double L_max = -1e99;
 168   int i_max = 0;
 169   if (mode == planar) {
 170     for (int i = 0; i < nodes.size(); ++i) {
 171       total_dist[i] = total_layers;
 172       vec3_t x_origin, x_target;
 173       m_Input->GetPoint(nodes[i],x_origin.data());
 174       x_target = x_origin + total_layers*node_normals[i];
 175       double L = x_target*fixed_normal;
 176       if (L > L_max) {
 177         i_max = i;
 178         L_max = L;
 179       }
 180     }
 181     vec3_t x_far;
 182     m_Input->GetPoint(nodes[i_max],x_far.data());
 183     x_far += min_dist*fixed_normal;
 184     for (int i = 0; i < nodes.size(); ++i) {
 185       total_dist[i] = total_layers;
 186       if (mode == planar) {
 187         vec3_t x_origin;
 188         m_Input->GetPoint(nodes[i],x_origin.data());
 189         total_dist[i] = (x_far-x_origin)*fixed_normal;
 190       }
 191     }
 192   }
 193   for (int i_layer = 0; i_layer < layer_y.size() - 1; ++i_layer) {
 194     for (int i = 0; i < n1.size(); ++i) {
 195       n1[i] = n2[i];
 196       n2[i] = i_layer*nodes.size() + i + m_Input->GetNumberOfPoints();
 197       vec3_t x1, x2;
 198       m_Output->GetPoint(n1[i],x1.data());
 199       if (mode == rotation) {
 200         x2 = x1 - origin;
 201         double alpha = (layer_y[i_layer + 1] - layer_y[i_layer])*M_PI/180.0;
 202         x2 = GeometryTools::rotate(x2,axis,alpha);
 203         x2 += origin;
 204       } else {
 205         double dist = (layer_y[i_layer + 1] - layer_y[i_layer]);
 206         if (mode == planar) {
 207           dist *= total_dist[i]/total_layers;
 208         }
 209         x2 = x1 + dist*node_normals[i];
 210       }
 211       m_Output->GetPoints()->SetPoint(n2[i],x2.data());
 212     }
 213
 214     for (int i_cell = 0; i_cell < cells.size(); ++i_cell) {
 215       vtkIdType *pts;
 216       vtkIdType  Npts;
 217       m_Input->GetCellPoints(cells[i_cell], Npts, pts);
 218       QVector<vtkIdType> surf_pts(Npts);
 219       for (int i_pts = 0; i_pts < Npts; ++i_pts) {
 220         surf_pts[i_pts] = _nodes[pts[i_pts]];
 221       }
 222       for (int i_surf_pts = 0; i_surf_pts < Npts; ++i_surf_pts) {
 223         vtkIdType p1 = surf_pts[i_surf_pts];
 224         vtkIdType p2;
 225         if (i_surf_pts < Npts - 1) {
 226           p2 = surf_pts[i_surf_pts + 1];
 227         } else {
 228           p2 = surf_pts[0];
 229         }
 230         bool add_bd = false;
 231         {
 232           int N = 0;
 233           int cell;
 234           foreach(cell, n2c[p1]) {
 235             if (n2c[p2].contains(cell)) ++N;
 236           }
 237           if (N == 1) add_bd = true;
 238         }
 239         if (add_bd) {
 240           vtkIdType quad_pts[4];
 241           quad_pts[0] = n1[p1];
 242           quad_pts[1] = n1[p2];
 243           quad_pts[2] = n2[p2];
 244           quad_pts[3] = n2[p1];
 245           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_QUAD,4,quad_pts);
 246           cell_code2->SetValue(id_new_cell, new_bc);
 247           orgdir->SetValue(id_new_cell, 0);
 248           curdir->SetValue(id_new_cell, 0);
 249           voldir->SetValue(id_new_cell, 0);
 250         }
 251       }
 252       if (Npts == 3) {
 253         {
 254           vtkIdType pri_pts[6];
 255           pri_pts[0] = n1[surf_pts[0]];
 256           pri_pts[1] = n1[surf_pts[2]];
 257           pri_pts[2] = n1[surf_pts[1]];
 258           pri_pts[3] = n2[surf_pts[0]];
 259           pri_pts[4] = n2[surf_pts[2]];
 260           pri_pts[5] = n2[surf_pts[1]];
 261           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_WEDGE,6,pri_pts);
 262           cell_code2->SetValue(id_new_cell, 0);
 263           orgdir->SetValue(id_new_cell, 0);
 264           curdir->SetValue(id_new_cell, 0);
 265           voldir->SetValue(id_new_cell, 0);
 266         }
 267         if (i_layer == layer_y.size() - 2) {
 268           vtkIdType tri_pts[3];
 269           tri_pts[0] = n2[surf_pts[0]];
 270           tri_pts[1] = n2[surf_pts[1]];
 271           tri_pts[2] = n2[surf_pts[2]];
 272           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_TRIANGLE,3,tri_pts);
 273           cell_code2->SetValue(id_new_cell, cell_code1->GetValue(cells[i_cell]));
 274           orgdir->SetValue(id_new_cell, 0);
 275           curdir->SetValue(id_new_cell, 0);
 276           voldir->SetValue(id_new_cell, 0);
 277         }
 278       }
 279       if (Npts == 4) {
 280         {
 281           vtkIdType pri_pts[8];
 282           pri_pts[0] = n1[surf_pts[0]];
 283           pri_pts[1] = n1[surf_pts[1]];
 284           pri_pts[2] = n1[surf_pts[2]];
 285           pri_pts[3] = n1[surf_pts[3]];
 286           pri_pts[4] = n2[surf_pts[0]];
 287           pri_pts[5] = n2[surf_pts[1]];
 288           pri_pts[6] = n2[surf_pts[2]];
 289           pri_pts[7] = n2[surf_pts[3]];
 290           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_HEXAHEDRON,8,pri_pts);
 291           cell_code2->SetValue(id_new_cell, 0);
 292           orgdir->SetValue(id_new_cell, 0);
 293           curdir->SetValue(id_new_cell, 0);
 294           voldir->SetValue(id_new_cell, 0);
 295         }
 296         if (i_layer == layer_y.size() - 2) {
 297           vtkIdType quad_pts[4];
 298           quad_pts[0] = n2[surf_pts[0]];
 299           quad_pts[1] = n2[surf_pts[1]];
 300           quad_pts[2] = n2[surf_pts[2]];
 301           quad_pts[3] = n2[surf_pts[3]];
 302           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_QUAD,4,quad_pts);
 303           cell_code2->SetValue(id_new_cell, cell_code1->GetValue(cells[i_cell]));
 304           orgdir->SetValue(id_new_cell, 0);
 305           curdir->SetValue(id_new_cell, 0);
 306           voldir->SetValue(id_new_cell, 0);
 307         }
 308       }
 309       if (Npts > 4) {
 310         {
 311           EG_VTKSP(vtkIdList, stream);
 312           stream->SetNumberOfIds(1 + 2*(Npts + 1) + Npts*5);
 313           vtkIdType id = 0;
 314           stream->SetId(id++, 2 + Npts);
 315
 316           // bottom face
 317           stream->SetId(id++, Npts);
 318           for (int i = Npts - 1; i >= 0; --i) {
 319             stream->SetId(id++, n1[surf_pts[i]]);
 320           }
 321
 322           // top face
 323           stream->SetId(id++, Npts);
 324           for (int i = 0; i < Npts; ++i) {
 325             stream->SetId(id++, n2[surf_pts[i]]);
 326           }
 327
 328           // side faces
 329           for (int i = 0; i < Npts; ++i) {
 330             stream->SetId(id++, 4);
 331             vtkIdType p1 = n1[surf_pts[i]];
 332             vtkIdType p2 = n1[surf_pts[0]];
 333             vtkIdType p3 = n2[surf_pts[0]];
 334             vtkIdType p4 = n2[surf_pts[i]];
 335             if (i < Npts - 1) {
 336               p2 = n1[surf_pts[i+1]];
 337               p3 = n2[surf_pts[i+1]];
 338             }
 339             stream->SetId(id++, p1);
 340             stream->SetId(id++, p2);
 341             stream->SetId(id++, p3);
 342             stream->SetId(id++, p4);
 343           }
 344
 345           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_POLYHEDRON, stream);
 346           cell_code2->SetValue(id_new_cell, 0);
 347           orgdir->SetValue(id_new_cell, 0);
 348           curdir->SetValue(id_new_cell, 0);
 349           voldir->SetValue(id_new_cell, 0);
 350         }
 351         if (i_layer == layer_y.size() - 2) {
 352           EG_VTKSP(vtkIdList, poly_pts);
 353           poly_pts->SetNumberOfIds(Npts);
 354           for (vtkIdType id = 0; id < Npts; ++id) {
 355             poly_pts->SetId(id, n2[surf_pts[id]]);
 356           }
 357           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_POLYGON, poly_pts);
 358           cell_code2->SetValue(id_new_cell, cell_code1->GetValue(cells[i_cell]));
 359           orgdir->SetValue(id_new_cell, 0);
 360           curdir->SetValue(id_new_cell, 0);
 361           voldir->SetValue(id_new_cell, 0);
 362         }
 363       }
 364     }
 365   }
 366
 367
 368   for (vtkIdType nodeId = 0; nodeId < m_Input->GetNumberOfPoints(); ++nodeId) {
 369     vec3_t x;
 370     m_Input->GetPoints()->GetPoint(nodeId, x.data());
 371     m_Output->GetPoints()->SetPoint(nodeId, x.data());
 372   }
 373
 374   // copy all original cells that were not part of the extrusion
 375   for (vtkIdType id_cell = 0; id_cell < m_Input->GetNumberOfCells(); ++id_cell) {
 376     if (_cells[id_cell] == -1 || !m_RemoveInternalFaces) {
 377       vtkIdType id_new_cell = copyCell(m_Input, id_cell, m_Output);
 378       copyCellData(m_Input, id_cell, m_Output, id_new_cell);
 379     }
 380   }
 381
 382   // close boundary where no volume cells were present in the original grid
 383
 384   for (int i_cell = 0; i_cell < cells.size(); ++i_cell) {
 385     if (is_boundary[i_cell]) {
 386       vtkIdType id_cell = cells[i_cell];
 387       vtkIdType *pts;
 388       vtkIdType  Npts;
 389       m_Input->GetCellPoints(id_cell, Npts, pts);
 390       vtkIdType id_new_cell = m_Output->InsertNextCell(m_Input->GetCellType(id_cell), Npts, pts);
 391       m_Output->GetCellPoints(id_new_cell, Npts, pts);
 392       QVector<vtkIdType> nodes(Npts);
 393       for (vtkIdType j = 0; j < Npts; ++j) nodes[j]          = pts[j];
 394       for (vtkIdType j = 0; j < Npts; ++j) pts[Npts - j - 1] = nodes[j];
 395       copyCellData(m_Input, id_cell, m_Output, id_new_cell);
 396     }
 397   }
 398
 399   UpdateCellIndex(m_Output);
 400 }
 401
 402 void vtkEgNormalExtrusion::SetLayers(const QVector<double> &y)
 403 {
 404   layer_y.resize(y.size());
 405   qCopy(y.begin(), y.end(), layer_y.begin());
 406 }
 407
 408