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