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
  22 #include "vtkEgNormalExtrusion.h"
  23 #include "geometrytools.h"
  24
  25 vtkStandardNewMacro(vtkEgNormalExtrusion)
  26
  27 vtkEgNormalExtrusion::vtkEgNormalExtrusion()
  28 {
  29   m_Mode = normal;
  30   m_ScaleX = 1;
  31   m_ScaleY = 1;
  32   m_ScaleZ = 1;
  33   m_RemoveInternalFaces = true;
  34   m_Curve1 = NULL;
  35   m_Curve2 = NULL;
  36   m_OrthoExtrusion = false;
  37   m_CustomBottomBc = -1;
  38   m_CustomSideBc = -1;
  39   m_CustomTopBc = -1;
  40 }
  41
  42 int vtkEgNormalExtrusion::getNewBc(MeshPartition *part, vtkIdType id_node1, vtkIdType id_node2)
  43 {
  44   QVector<vtkIdType> nodes(2);
  45   nodes[0] = id_node1;
  46   nodes[1] = id_node2;
  47   QSet<int> common_bcs = part->getCommonBcs(nodes);
  48   QSet<int> potential_bcs = common_bcs - m_BoundaryCodes;
  49   if (potential_bcs.size() > 0) {
  50     return potential_bcs.values().first();
  51   }
  52   return m_UnknownBc;
  53 }
  54
  55 void vtkEgNormalExtrusion::SetCurves(Curve *curve1, Curve *curve2)
  56 {
  57   m_Curve1 = curve1;
  58   m_Curve2 = curve2;
  59   if (m_Curve1 && m_Curve2) {
  60     m_Mode = curve;
  61   } else {
  62     m_Mode = normal;
  63   }
  64 }
  65
  66 bool vtkEgNormalExtrusion::PrepareCurveExtrusion(QVector<vtkIdType> nodes)
  67 {
  68   if (!m_Curve1 || !m_Curve2) {
  69     return false;
  70   }
  71   mat3_t A;
  72   if (m_OrthoExtrusion) {
  73     A = m_Curve1->computeOrthoBase(0, m_Curve2);
  74   } else {
  75     A = m_Curve2->computeBase(0, m_Curve2);
  76   }
  77   A = A.inverse();
  78   m_InitialBaseCoords.resize(nodes.size());
  79   vec3_t x0 = m_Curve1->position(0);
  80   for (int i = 0; i < nodes.size(); ++i) {
  81     vec3_t x;
  82     m_Input->GetPoint(nodes[i], x.data());
  83     x -= x0;
  84     m_InitialBaseCoords[i] = A*x;
  85   }
  86   return true;
  87 }
  88
  89 mat3_t vtkEgNormalExtrusion::ComputeBase(double l)
  90 {
  91   mat3_t A;
  92   if (m_OrthoExtrusion) {
  93     A = m_Curve1->computeOrthoBase(l, m_Curve2);
  94   } else {
  95     A = m_Curve2->computeBase(l, m_Curve2);
  96   }
  97   return A;
  98 }
  99
 100 void vtkEgNormalExtrusion::ExecuteEg()
 101 {
 102   QVector<vtkIdType> cells, nodes, n1, n2;
 103   QVector<vec3_t> cell_normals, node_normals;
 104   ExtractBoundary(cells, nodes, m_BoundaryCodes, m_Input);
 105   if (m_Mode == normal) {
 106     computeNormals(cell_normals, node_normals, cells, nodes,m_Input);
 107   } else if (m_Mode == cylinder) {
 108     axis.normalise();
 109     node_normals.resize(nodes.size());
 110     for (int i = 0; i < node_normals.size(); ++i) {
 111       vec3_t x;
 112       m_Input->GetPoint(nodes[i],x.data());
 113       vec3_t x0 = x - ((x-origin)*axis)*axis;
 114       node_normals[i] = x0 - origin;
 115       node_normals[i].normalise();
 116     }
 117   } else if ((m_Mode == fixed) || (m_Mode == planar)) {
 118     fixed_normal.normalise();
 119     node_normals.resize(nodes.size());
 120     for (int i = 0; i < node_normals.size(); ++i) {
 121       node_normals[i] = fixed_normal;
 122     }
 123   }
 124   for (int i = 0; i < node_normals.size(); ++ i) {
 125     node_normals[i][0] *= m_ScaleX;
 126     node_normals[i][1] *= m_ScaleY;
 127     node_normals[i][2] *= m_ScaleZ;
 128   }
 129   n1.resize(nodes.size());
 130   n2.resize(nodes.size());
 131   if (m_Mode == curve) {
 132     if (!PrepareCurveExtrusion(nodes)) {
 133       return;
 134     }
 135   }
 136
 137   // mapping
 138   QVector<int> _cells, _nodes;
 139   createNodeMapping(nodes, _nodes, m_Input);
 140   createCellMapping(cells, _cells, m_Input);
 141   QVector<QSet<int> > n2c;
 142   createNodeToCell(cells, nodes, _nodes, n2c, m_Input);
 143
 144   vtkIdType NnewNodes = m_Input->GetNumberOfPoints() + (layer_y.size()-1)*nodes.size();
 145   vtkIdType NnewCells = m_Input->GetNumberOfCells() + (layer_y.size()-1)*cells.size();
 146
 147   // count the number of new surface elements (side walls)
 148   for (int i_cell = 0; i_cell < cells.size(); ++i_cell) {
 149     vtkIdType *pts;
 150     vtkIdType  Npts;
 151     m_Input->GetCellPoints(cells[i_cell], Npts, pts);
 152     QVector<vtkIdType> surf_pts(Npts);
 153     for (int i_pts = 0; i_pts < Npts; ++i_pts) {
 154       surf_pts[i_pts] = _nodes[pts[i_pts]];
 155     }
 156     for (int i_surf_pts = 0; i_surf_pts < Npts; ++i_surf_pts) {
 157       vtkIdType p1 = surf_pts[i_surf_pts];
 158       vtkIdType p2;
 159       if (i_surf_pts < Npts - 1) {
 160         p2 = surf_pts[i_surf_pts + 1];
 161       } else {
 162         p2 = surf_pts[0];
 163       }
 164       bool add_bd = false;
 165       {
 166         int N = 0;
 167         int cell;
 168         foreach(cell, n2c[p1]) {
 169           if (n2c[p2].contains(cell)) ++N;
 170         }
 171         if (N == 1) add_bd = true;
 172       }
 173       if (add_bd) {
 174         NnewCells += layer_y.size();
 175       }
 176     }
 177   }
 178
 179   // count the number of new surface elements (base)
 180   QVector<bool> is_boundary;
 181   is_boundary.fill(false, cells.size());
 182   {
 183     QVector<int> nvol;
 184     nvol.fill(0, nodes.size());
 185     for (vtkIdType id_cell = 0; id_cell < m_Input->GetNumberOfCells(); ++id_cell) {
 186       if (isVolume(id_cell, m_Input)) {
 187         vtkIdType Npts, *pts;
 188         m_Input->GetCellPoints(id_cell, Npts, pts);
 189         for (int i = 0; i < Npts; ++i) {
 190           if (_nodes[pts[i]] >= 0) {
 191             ++nvol[_nodes[pts[i]]];
 192           }
 193         }
 194       }
 195     }
 196     for (int i_cell = 0; i_cell < cells.size(); ++i_cell) {
 197       vtkIdType id_cell = cells[i_cell];
 198       vtkIdType Npts, *pts;
 199       m_Input->GetCellPoints(id_cell, Npts, pts);
 200       for (int i = 0; i < Npts; ++i) {
 201         if (nvol[_nodes[pts[i]]] == 0) {
 202           is_boundary[i_cell] = true;
 203         }
 204       }
 205       if (is_boundary[i_cell] || !m_RemoveInternalFaces) {
 206         ++NnewCells;
 207       }
 208     }
 209   }
 210
 211   // allocate memory for the new grid
 212   allocateGrid(m_Output, NnewCells, NnewNodes);
 213
 214   // boundary conditions
 215   EG_VTKDCC(vtkIntArray, cell_code1, m_Input, "cell_code");
 216   EG_VTKDCC(vtkIntArray, cell_code2, m_Output, "cell_code");
 217   EG_VTKDCC(vtkIntArray, orgdir, m_Output, "cell_orgdir");
 218   EG_VTKDCC(vtkIntArray, voldir, m_Output, "cell_voldir");
 219   EG_VTKDCC(vtkIntArray, curdir, m_Output, "cell_curdir");
 220
 221   m_UnknownBc = 1;
 222   for (vtkIdType id_cell = 0; id_cell < m_Input->GetNumberOfCells(); ++id_cell) {
 223     if (isSurface(id_cell, m_Input)) {
 224       if (cell_code1->GetValue(id_cell) >= m_UnknownBc) {
 225         m_UnknownBc = cell_code1->GetValue(id_cell) + 1;
 226       }
 227     }
 228   }
 229
 230   QMap<QPair<int,int>,int> new_bcs;
 231
 232   MeshPartition part(m_Input, true);
 233
 234   for (int i = 0; i < nodes.size(); ++i) {
 235     n2[i] = nodes[i];
 236     vec3_t x;
 237     m_Input->GetPoint(nodes[i],x.data());
 238     m_Output->GetPoints()->SetPoint(n2[i],x.data());
 239   }
 240   double total_layers = layer_y[layer_y.size()-1] - layer_y[0];
 241   QVector<double> total_dist(nodes.size());
 242   double L_max = -1e99;
 243   int i_max = 0;
 244   if (m_Mode == planar) {
 245     for (int i = 0; i < nodes.size(); ++i) {
 246       total_dist[i] = total_layers;
 247       vec3_t x_origin, x_target;
 248       m_Input->GetPoint(nodes[i],x_origin.data());
 249       x_target = x_origin + total_layers*node_normals[i];
 250       double L = x_target*fixed_normal;
 251       if (L > L_max) {
 252         i_max = i;
 253         L_max = L;
 254       }
 255     }
 256     vec3_t x_far;
 257     m_Input->GetPoint(nodes[i_max],x_far.data());
 258     x_far += min_dist*fixed_normal;
 259     for (int i = 0; i < nodes.size(); ++i) {
 260       total_dist[i] = total_layers;
 261       if (m_Mode == planar) {
 262         vec3_t x_origin;
 263         m_Input->GetPoint(nodes[i],x_origin.data());
 264         total_dist[i] = (x_far-x_origin)*fixed_normal;
 265       }
 266     }
 267   }
 268   for (int i_layer = 0; i_layer < layer_y.size() - 1; ++i_layer) {
 269     mat3_t current_base;
 270     vec3_t current_origin;
 271     if (m_Mode == curve) {
 272       current_origin = m_Curve1->position(layer_y[i_layer + 1]);
 273       current_base = ComputeBase(layer_y[i_layer + 1]);
 274     }
 275     for (int i = 0; i < n1.size(); ++i) {
 276       n1[i] = n2[i];
 277       n2[i] = i_layer*nodes.size() + i + m_Input->GetNumberOfPoints();
 278       vec3_t x1, x2;
 279       m_Output->GetPoint(n1[i],x1.data());
 280       if (m_Mode == rotation) {
 281
 282         x2 = x1 - origin;
 283         double alpha = (layer_y[i_layer + 1] - layer_y[i_layer])*M_PI/180.0;
 284         x2 = GeometryTools::rotate(x2,axis,alpha);
 285         x2 += origin;
 286
 287       } else if (m_Mode == curve) {
 288
 289         x2 = current_origin + current_base*m_InitialBaseCoords[i];
 290
 291       } else {
 292
 293         double dist = (layer_y[i_layer + 1] - layer_y[i_layer]);
 294         if (m_Mode == planar) {
 295           dist *= total_dist[i]/total_layers;
 296         }
 297         x2 = x1 + dist*node_normals[i];
 298
 299       }
 300       m_Output->GetPoints()->SetPoint(n2[i],x2.data());
 301     }
 302
 303     for (int i_cell = 0; i_cell < cells.size(); ++i_cell) {
 304       vtkIdType *pts;
 305       vtkIdType  Npts;
 306       m_Input->GetCellPoints(cells[i_cell], Npts, pts);
 307       QVector<vtkIdType> surf_pts(Npts);
 308       for (int i_pts = 0; i_pts < Npts; ++i_pts) {
 309         surf_pts[i_pts] = _nodes[pts[i_pts]];
 310       }
 311       for (int i_surf_pts = 0; i_surf_pts < Npts; ++i_surf_pts) {
 312         vtkIdType p1 = surf_pts[i_surf_pts];
 313         vtkIdType p2;
 314         if (i_surf_pts < Npts - 1) {
 315           p2 = surf_pts[i_surf_pts + 1];
 316         } else {
 317           p2 = surf_pts[0];
 318         }
 319         bool add_bd = false;
 320         {
 321           int N = 0;
 322           int cell;
 323           foreach(cell, n2c[p1]) {
 324             if (n2c[p2].contains(cell)) ++N;
 325           }
 326           if (N == 1) add_bd = true;
 327         }
 328         if (add_bd) {
 329           int new_bc = new_bcs[QPair<int,int>(p1,p2)];
 330           if (new_bc == 0) {
 331             new_bc = getNewBc(&part, nodes[p1], nodes[p2]);
 332             new_bcs[QPair<int,int>(p1,p2)] = new_bc;
 333           }
 334           if (m_CustomSideBc > 0) {
 335             new_bc = m_CustomSideBc;
 336           }
 337           vtkIdType quad_pts[4];
 338           quad_pts[0] = n1[p1];
 339           quad_pts[1] = n1[p2];
 340           quad_pts[2] = n2[p2];
 341           quad_pts[3] = n2[p1];
 342           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_QUAD,4,quad_pts);
 343           cell_code2->SetValue(id_new_cell, new_bc);
 344           orgdir->SetValue(id_new_cell, 0);
 345           curdir->SetValue(id_new_cell, 0);
 346           voldir->SetValue(id_new_cell, 0);
 347         }
 348       }
 349       if (Npts == 3) {
 350         {
 351           vtkIdType pri_pts[6];
 352           pri_pts[0] = n1[surf_pts[0]];
 353           pri_pts[1] = n1[surf_pts[2]];
 354           pri_pts[2] = n1[surf_pts[1]];
 355           pri_pts[3] = n2[surf_pts[0]];
 356           pri_pts[4] = n2[surf_pts[2]];
 357           pri_pts[5] = n2[surf_pts[1]];
 358           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_WEDGE,6,pri_pts);
 359           cell_code2->SetValue(id_new_cell, 0);
 360           orgdir->SetValue(id_new_cell, 0);
 361           curdir->SetValue(id_new_cell, 0);
 362           voldir->SetValue(id_new_cell, 0);
 363         }
 364         if (i_layer == layer_y.size() - 2) {
 365           vtkIdType tri_pts[3];
 366           tri_pts[0] = n2[surf_pts[0]];
 367           tri_pts[1] = n2[surf_pts[1]];
 368           tri_pts[2] = n2[surf_pts[2]];
 369           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_TRIANGLE,3,tri_pts);
 370           if (m_CustomTopBc > 0) {
 371             cell_code2->SetValue(id_new_cell, m_CustomTopBc);
 372           } else {
 373             cell_code2->SetValue(id_new_cell, cell_code1->GetValue(cells[i_cell]));
 374           }
 375           orgdir->SetValue(id_new_cell, 0);
 376           curdir->SetValue(id_new_cell, 0);
 377           voldir->SetValue(id_new_cell, 0);
 378         }
 379       }
 380       if (Npts == 4) {
 381         {
 382           vtkIdType pri_pts[8];
 383           pri_pts[0] = n1[surf_pts[0]];
 384           pri_pts[1] = n1[surf_pts[1]];
 385           pri_pts[2] = n1[surf_pts[2]];
 386           pri_pts[3] = n1[surf_pts[3]];
 387           pri_pts[4] = n2[surf_pts[0]];
 388           pri_pts[5] = n2[surf_pts[1]];
 389           pri_pts[6] = n2[surf_pts[2]];
 390           pri_pts[7] = n2[surf_pts[3]];
 391           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_HEXAHEDRON,8,pri_pts);
 392           cell_code2->SetValue(id_new_cell, 0);
 393           orgdir->SetValue(id_new_cell, 0);
 394           curdir->SetValue(id_new_cell, 0);
 395           voldir->SetValue(id_new_cell, 0);
 396         }
 397         if (i_layer == layer_y.size() - 2) {
 398           vtkIdType quad_pts[4];
 399           quad_pts[0] = n2[surf_pts[0]];
 400           quad_pts[1] = n2[surf_pts[1]];
 401           quad_pts[2] = n2[surf_pts[2]];
 402           quad_pts[3] = n2[surf_pts[3]];
 403           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_QUAD,4,quad_pts);
 404           if (m_CustomTopBc > 0) {
 405             cell_code2->SetValue(id_new_cell, m_CustomTopBc);
 406           } else {
 407             cell_code2->SetValue(id_new_cell, cell_code1->GetValue(cells[i_cell]));
 408           }
 409           orgdir->SetValue(id_new_cell, 0);
 410           curdir->SetValue(id_new_cell, 0);
 411           voldir->SetValue(id_new_cell, 0);
 412         }
 413       }
 414       if (Npts > 4) {
 415         {
 416           EG_VTKSP(vtkIdList, stream);
 417           stream->SetNumberOfIds(1 + 2*(Npts + 1) + Npts*5);
 418           vtkIdType id = 0;
 419           stream->SetId(id++, 2 + Npts);
 420
 421           // bottom face
 422           stream->SetId(id++, Npts);
 423           for (int i = Npts - 1; i >= 0; --i) {
 424             stream->SetId(id++, n1[surf_pts[i]]);
 425           }
 426
 427           // top face
 428           stream->SetId(id++, Npts);
 429           for (int i = 0; i < Npts; ++i) {
 430             stream->SetId(id++, n2[surf_pts[i]]);
 431           }
 432
 433           // side faces
 434           for (int i = 0; i < Npts; ++i) {
 435             stream->SetId(id++, 4);
 436             vtkIdType p1 = n1[surf_pts[i]];
 437             vtkIdType p2 = n1[surf_pts[0]];
 438             vtkIdType p3 = n2[surf_pts[0]];
 439             vtkIdType p4 = n2[surf_pts[i]];
 440             if (i < Npts - 1) {
 441               p2 = n1[surf_pts[i+1]];
 442               p3 = n2[surf_pts[i+1]];
 443             }
 444             stream->SetId(id++, p1);
 445             stream->SetId(id++, p2);
 446             stream->SetId(id++, p3);
 447             stream->SetId(id++, p4);
 448           }
 449
 450           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_POLYHEDRON, stream);
 451           cell_code2->SetValue(id_new_cell, 0);
 452           orgdir->SetValue(id_new_cell, 0);
 453           curdir->SetValue(id_new_cell, 0);
 454           voldir->SetValue(id_new_cell, 0);
 455         }
 456         if (i_layer == layer_y.size() - 2) {
 457           EG_VTKSP(vtkIdList, poly_pts);
 458           poly_pts->SetNumberOfIds(Npts);
 459           for (vtkIdType id = 0; id < Npts; ++id) {
 460             poly_pts->SetId(id, n2[surf_pts[id]]);
 461           }
 462           vtkIdType id_new_cell = m_Output->InsertNextCell(VTK_POLYGON, poly_pts);
 463           if (m_CustomTopBc > 0) {
 464             cell_code2->SetValue(id_new_cell, m_CustomTopBc);
 465           } else {
 466             cell_code2->SetValue(id_new_cell, cell_code1->GetValue(cells[i_cell]));
 467           }
 468           orgdir->SetValue(id_new_cell, 0);
 469           curdir->SetValue(id_new_cell, 0);
 470           voldir->SetValue(id_new_cell, 0);
 471         }
 472       }
 473     }
 474   }
 475
 476
 477   for (vtkIdType nodeId = 0; nodeId < m_Input->GetNumberOfPoints(); ++nodeId) {
 478     vec3_t x;
 479     m_Input->GetPoints()->GetPoint(nodeId, x.data());
 480     m_Output->GetPoints()->SetPoint(nodeId, x.data());
 481   }
 482
 483   // copy all original cells that were not part of the extrusion
 484   for (vtkIdType id_cell = 0; id_cell < m_Input->GetNumberOfCells(); ++id_cell) {
 485     if (_cells[id_cell] == -1 || !m_RemoveInternalFaces) {
 486       vtkIdType id_new_cell = copyCell(m_Input, id_cell, m_Output);
 487       copyCellData(m_Input, id_cell, m_Output, id_new_cell);
 488     }
 489   }
 490
 491   // close boundary where no volume cells were present in the original grid
 492
 493   for (int i_cell = 0; i_cell < cells.size(); ++i_cell) {
 494     if (is_boundary[i_cell]) {
 495       vtkIdType id_cell = cells[i_cell];
 496       vtkIdType *pts;
 497       vtkIdType  Npts;
 498       m_Input->GetCellPoints(id_cell, Npts, pts);
 499       vtkIdType id_new_cell = m_Output->InsertNextCell(m_Input->GetCellType(id_cell), Npts, pts);
 500       m_Output->GetCellPoints(id_new_cell, Npts, pts);
 501       QVector<vtkIdType> nodes(Npts);
 502       for (vtkIdType j = 0; j < Npts; ++j) nodes[j]          = pts[j];
 503       for (vtkIdType j = 0; j < Npts; ++j) pts[Npts - j - 1] = nodes[j];
 504       copyCellData(m_Input, id_cell, m_Output, id_new_cell);
 505       if (m_CustomBottomBc > 0) {
 506         cell_code2->SetValue(id_new_cell, m_CustomBottomBc);
 507       }
 508     }
 509   }
 510
 511   updateCellIndex(m_Output);
 512 }
 513
 514 void vtkEgNormalExtrusion::SetLayers(const QVector<double> &y)
 515 {
 516   layer_y.resize(y.size());
 517   qCopy(y.begin(), y.end(), layer_y.begin());
 518 }