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