src/libengrid/fixcadgeometry.cpp

   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 "fixcadgeometry.h"
  22 #include "surfacemesher.h"
  23 #include "vertexmeshdensity.h"
  24 #include "guimainwindow.h"
  25 #include "deletestraynodes.h"
  26 #include "correctsurfaceorientation.h"
  27 #include "swaptriangles.h"
  28
  29 #include <vtkMath.h>
  30
  31 #include "geometrytools.h"
  32 using namespace GeometryTools;
  33
  34 #include <QtDebug>
  35 #include <iostream>
  36 using namespace std;
  37
  38 FixCadGeometry::FixCadGeometry()
  39 {
  40   EG_TYPENAME;
  41   m_PerformGeometricTests = true;
  42   m_UseProjectionForSmoothing = false;
  43   m_UseNormalCorrectionForSmoothing = true;
  44   m_OriginalFeatureAngle = m_FeatureAngle;
  45   m_FeatureAngle = GeometryTools::deg2rad(0.5);
  46   m_AllowFeatureEdgeSwapping = false;
  47   m_AllowSmallAreaSwapping = true;
  48   m_NumDelaunaySweeps = 1;
  49 }
  50
  51 void FixCadGeometry::callMesher()
  52 {
  53   setDesiredLength();
  54   if (m_BoundaryCodes.size() == 0) {
  55     return;
  56   }
  57   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired");
  58   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
  59     characteristic_length_desired->SetValue(id_node, 1e-6);
  60   }
  61   int num_deleted = 0;
  62   int iter = 0;
  63   bool done = false;
  64   int count = 0;
  65   while (!done) {
  66     SwapTriangles swap;
  67     swap.setGrid(m_Grid);
  68     swap.setRespectBC(true);
  69     swap.setFeatureSwap(m_AllowFeatureEdgeSwapping);
  70     swap.setFeatureAngle(m_FeatureAngle);
  71     swap.setMaxNumLoops(1);
  72     swap.setSmallAreaSwap(m_AllowSmallAreaSwapping);
  73     swap.setDelaunayThreshold(1e6);
  74     swap.setVerboseOn();
  75     QSet<int> rest_bcs = GuiMainWindow::pointer()->getAllBoundaryCodes();
  76     rest_bcs -= m_BoundaryCodes;
  77     swap.setBoundaryCodes(rest_bcs);
  78     swap();
  79     count = 0;
  80     if (swap.getNumSwaps() == 0 || count >= 100) {
  81       done = true;
  82     }
  83   }
  84
  85   /*
  86   while (!done) {
  87     ++iter;
  88     cout << "CAD fix iteration " << iter << ":" << endl;
  89     setDesiredLength();
  90     customUpdateNodeInfo();
  91     setDesiredLength();
  92     int num_deleted = deleteNodes();
  93     cout << "  deleted nodes  : " << num_deleted << endl;
  94     swap();
  95     setDesiredLength();
  96     done = (iter >= 20) || (num_deleted == 0);
  97     cout << "  total nodes    : " << m_Grid->GetNumberOfPoints() << endl;
  98     cout << "  total cells    : " << m_Grid->GetNumberOfCells() << endl;
  99   }
 100   */
 101   createIndices(m_Grid);
 102   customUpdateNodeInfo();
 103   setDesiredLength();
 104 }
 105
 106 void FixCadGeometry::setDesiredLength(double L)
 107 {
 108   setAllSurfaceCells();
 109   l2g_t  nodes = getPartNodes();
 110
 111   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired,   m_Grid, "node_meshdensity_desired");
 112   EG_VTKDCN(vtkIntArray,    characteristic_length_specified, m_Grid, "node_specified_density");
 113
 114   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 115     vtkIdType id_node = nodes[i_nodes];
 116     characteristic_length_specified->SetValue(id_node, 0);
 117     characteristic_length_desired->SetValue(id_node, L);
 118   }
 119 }
 120
 121 void FixCadGeometry::customUpdateNodeInfo()
 122 {
 123   cout << "updating node information ..." << endl;
 124   setAllCells();
 125   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 126   l2g_t cells = getPartCells();
 127   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 128     node_type->SetValue(id_node, EG_FIXED_VERTEX);
 129   }
 130   foreach (vtkIdType id_cell, cells) {
 131     if (isSurface(id_cell, m_Grid)) {
 132       vtkIdType N_pts, *pts;
 133       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 134       if (N_pts == 3) {
 135         vec3_t n1 = cellNormal(m_Grid, id_cell);
 136         QVector<int> num_bad_edges(3, 0);
 137         for (int i = 0; i < N_pts; ++i) {
 138           int i1 = i;
 139           int i2 = 0;
 140           if (i < N_pts - 1) {
 141             i2= i+1;
 142           }
 143           vec3_t x1, x2;
 144           m_Grid->GetPoint(pts[i1], x1.data());
 145           m_Grid->GetPoint(pts[i2], x2.data());
 146           double L = (x1 - x2).abs();
 147           vtkIdType id_neigh_cell = m_Part.c2cGG(id_cell, i);
 148           if (id_neigh_cell != -1) {
 149             bool bad_edge = false;
 150             vec3_t n2 = cellNormal(m_Grid, id_neigh_cell);
 151             if (angle(n1, n2) > deg2rad(179.5)) {
 152               bad_edge = true;
 153             }
 154             if (bad_edge) {
 155               ++num_bad_edges[i1];
 156               ++num_bad_edges[i2];
 157             }
 158           }
 159         }
 160         for (int i = 0; i < N_pts; ++i) {
 161           if (num_bad_edges[i] >= 2) {
 162             node_type->SetValue(pts[i], EG_SIMPLE_VERTEX);
 163           }
 164         }
 165       }
 166     }
 167   }
 168   cout << "done." << endl;
 169 }
 170
 171 void FixCadGeometry::copyFaces(const QVector<bool> &copy_face)
 172 {
 173   int N_del = 0;
 174   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 175     if (!copy_face[id_cell]) {
 176       ++N_del;
 177     }
 178   }
 179   int N_cells = m_Grid->GetNumberOfCells() - N_del;
 180   int N_nodes = 0;
 181   QVector<bool> copy_node(m_Grid->GetNumberOfPoints(), false);
 182   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 183     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 184       if (copy_face[m_Part.n2cGG(id_node, i)]) {
 185         copy_node[id_node] = true;
 186         ++N_nodes;
 187         break;
 188       }
 189     }
 190   }
 191   EG_VTKSP(vtkUnstructuredGrid, new_grid);
 192   allocateGrid(new_grid, N_cells, N_nodes);
 193   QVector<vtkIdType> old2new(m_Grid->GetNumberOfPoints(), -1);
 194   vtkIdType id_new_node = 0;
 195   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 196     if (copy_node[id_node]) {
 197       old2new[id_node] = id_new_node;
 198       vec3_t x;
 199       m_Grid->GetPoint(id_node, x.data());
 200       new_grid->GetPoints()->SetPoint(id_new_node, x.data());
 201       copyNodeData(m_Grid, id_node, new_grid, id_new_node);
 202       ++id_new_node;
 203     }
 204   }
 205   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 206     if (copy_face[id_cell]) {
 207       vtkIdType N_pts, *pts;
 208       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 209       vtkIdType id_new_cell = new_grid->InsertNextCell(m_Grid->GetCellType(id_cell), N_pts, pts);
 210       copyCellData(m_Grid, id_cell, new_grid, id_new_cell);
 211     }
 212   }
 213   makeCopy(new_grid, m_Grid);
 214   cout << "  deleted " << N_del << " faces" << endl;
 215 }
 216
 217 void FixCadGeometry::fixNonManifold1()
 218 {
 219   cout << "fixing non-manifold edges\n  (pass 1) ..." << endl;
 220   QVector<bool> copy_face(m_Grid->GetNumberOfCells(), true);
 221   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 222     for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
 223       if (m_Part.c2cGG(id_cell, i) == -1) {
 224         copy_face[id_cell] = false;
 225         break;
 226       }
 227     }
 228   }
 229   copyFaces(copy_face);
 230   cout << "done." << endl;
 231 }
 232
 233 void FixCadGeometry::fixNonManifold2()
 234 {
 235   cout << "fixing non-manifold edges\n  (pass 2) ..." << endl;
 236   QVector<bool> copy_face(m_Grid->GetNumberOfCells(), true);
 237   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 238     vtkIdType N_pts, *pts;
 239     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 240     if (N_pts < 3) {
 241       copy_face[id_cell] = false;
 242       break;
 243     }
 244     QVector<QSet<vtkIdType> > n2c(N_pts);
 245     for (int i = 0; i < N_pts; ++i) {
 246       for (int j = 0; j < m_Part.n2cGSize(pts[i]); ++j) {
 247         n2c[i].insert(m_Part.n2cGG(pts[i], j));
 248       }
 249     }
 250     QSet<vtkIdType> faces = n2c[0];
 251     for (int i = 1; i < N_pts; ++i) {
 252       faces = faces.intersect(n2c[i]);
 253     }
 254     if (faces.size() > 1) {
 255       vtkIdType id_del = id_cell;
 256       foreach (vtkIdType id, faces) {
 257         id_del = min(id_del, id);
 258       }
 259       copy_face[id_del] = false;
 260     }
 261   }
 262   copyFaces(copy_face);
 263   cout << "done." << endl;
 264 }
 265
 266 void FixCadGeometry::markNonManifold()
 267 {
 268   cout << "marking non-manifold edges\n" << endl;
 269   int new_bc = 0;
 270   m_NumNonManifold = 0;
 271   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 272   QVector<bool> nm_face(m_Grid->GetNumberOfCells(), false);
 273   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 274     if (isSurface(id_cell, m_Grid)) {
 275       new_bc = max(new_bc, cell_code->GetValue(id_cell));
 276       vtkIdType N_pts, *pts;
 277       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 278       for (int i = 0; i < N_pts; ++i) {
 279         QSet <vtkIdType> edge_cells;
 280         int N = 0;
 281         if (i < N_pts - 1) {
 282           N = getEdgeCells(pts[i], pts[i+1], edge_cells);
 283         } else {
 284           N = getEdgeCells(pts[i], pts[0], edge_cells);
 285         }
 286         if (N != 2) {
 287           nm_face[id_cell] = true;
 288           ++m_NumNonManifold;
 289           break;
 290         }
 291       }
 292     }
 293   }
 294   m_BoundaryCodes = GuiMainWindow::pointer()->getAllBoundaryCodes();
 295   foreach (int bc, m_BoundaryCodes) {
 296     ++new_bc;
 297     QString bc_name = "NM_" + GuiMainWindow::pointer()->getBC(bc).getName();
 298     for (int level = 0; level < 2; ++level) {
 299       QVector<bool> new_nm_face = nm_face;
 300       for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 301         if (isSurface(id_cell, m_Grid)) {
 302           if ((cell_code->GetValue(id_cell)) == bc && nm_face[id_cell]) {
 303             vtkIdType N_pts, *pts;
 304             m_Grid->GetCellPoints(id_cell, N_pts, pts);
 305             for (int i = 0; i < N_pts; ++i) {
 306               for (int j = 0; j < m_Part.n2cGSize(pts[i]); ++j) {
 307                 vtkIdType id_neigh = m_Part.n2cGG(pts[i], j);
 308                 if (cell_code->GetValue(id_neigh) == bc ) {
 309                   new_nm_face[id_neigh] = true;
 310                 }
 311               }
 312             }
 313           }
 314         }
 315       }
 316       nm_face = new_nm_face;
 317     }
 318     GuiMainWindow::pointer()->setBC(new_bc, BoundaryCondition(bc_name, "patch", new_bc));
 319     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 320       if (isSurface(id_cell, m_Grid)) {
 321         if (cell_code->GetValue(id_cell) == bc && nm_face[id_cell]) {
 322           cell_code->SetValue(id_cell, new_bc);
 323         }
 324       }
 325     }
 326   }
 327   GuiMainWindow::pointer()->updateBoundaryCodes(true);
 328   cout << "done." << endl;
 329 }
 330
 331 void FixCadGeometry::operate()
 332 {
 333   {
 334     DeleteStrayNodes del;
 335     del();
 336   }
 337
 338   setAllCells();
 339   markNonManifold();
 340
 341   if (m_NumNonManifold == 0) {
 342
 343     //prepare BCmap
 344     EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
 345     QSet <int> bcs;
 346     l2g_t cells = m_Part.getCells();
 347     foreach (vtkIdType id_cell, cells) {
 348       if (isSurface(id_cell, m_Grid)) {
 349         bcs.insert(bc->GetValue(id_cell));
 350       }
 351     }
 352     QMap <int,int> BCmap;
 353     foreach(int bc, bcs) {
 354       BCmap[bc] = 1;
 355     }
 356
 357     //set density infinite
 358     VertexMeshDensity VMD;
 359     VMD.density = 1e99;
 360     VMD.BCmap = BCmap;
 361     qWarning()<<"VMD.BCmap="<<VMD.BCmap;
 362     m_VMDvector.push_back(VMD);
 363
 364     customUpdateNodeInfo();
 365
 366     // fix non manifold edges
 367     //fixNonManifold1();
 368     //fixNonManifold2();
 369
 370     //call surface mesher
 371     setGrid(m_Grid);
 372     setBoundaryCodes(bcs);
 373     setVertexMeshDensityVector(m_VMDvector);
 374     setDesiredLength();
 375
 376     callMesher();
 377
 378     /*
 379     // correct surface orientation
 380     CorrectSurfaceOrientation correct;
 381     correct.setGrid(m_Grid);
 382     correct.setAllCells();
 383     correct();
 384     */
 385   }
 386
 387   // finalise
 388   m_FeatureAngle = m_OriginalFeatureAngle;
 389   createIndices(m_Grid);
 390   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 391   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 392     node_type->SetValue(id_node, EG_SIMPLE_VERTEX);
 393   }
 394   updateNodeInfo();
 395   setDesiredLength();
 396 }
 397