src/libengrid/vtkEgNormalExtrusion.cxx

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