src/libengrid/sharpenedges.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2012 enGits GmbH                                     +
   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
  24 #include "sharpenedges.h"
  25 #include "guimainwindow.h"
  26 #include "laplacesmoother.h"
  27
  28 SharpenEdges::SharpenEdges()
  29 {
  30   EG_TYPENAME;
  31   m_PerformGeometricTests             = true;
  32   m_UseProjectionForSmoothing         = true;
  33   m_UseNormalCorrectionForSmoothing   = false;
  34   m_AllowFeatureEdgeSwapping          = false;
  35   m_RespectFeatureEdgesForDeleteNodes = false;
  36   m_FeatureAngle                      = deg2rad(180);
  37   m_EdgeAngle                         = deg2rad(180);
  38   m_FeatureAngleForDeleteNodes        = deg2rad(20);
  39   m_NumDelaunaySweeps                 = 5;
  40   m_NumSmoothSteps                    = 5;
  41 }
  42
  43
  44 void SharpenEdges::assignBCs()
  45 {
  46   QVector<int> bc_set(m_Grid->GetNumberOfCells(),0);
  47   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  48   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  49     bool certain_bc = true;
  50     vtkIdType N_pts, *pts;
  51     m_Grid->GetCellPoints(id_cell, N_pts, pts);
  52     for (int i = 0; i < N_pts; ++i) {
  53       vtkIdType id_node = pts[i];
  54       for (int j = 0; j < m_Part.n2cGSize(pts[i]); ++j) {
  55         if (cell_code->GetValue(id_cell) != cell_code->GetValue(m_Part.n2cGG(pts[i], j))) {
  56           certain_bc = false;
  57           break;
  58         }
  59       }
  60     }
  61     if (certain_bc) {
  62       bc_set[id_cell] = cell_code->GetValue(id_cell);
  63     }
  64   }
  65   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  66     if (bc_set[id_cell] == 0) {
  67       cell_code->SetValue(id_cell, 9999);
  68     }
  69   }
  70   GuiMainWindow::pointer()->storeSurfaceProjection();
  71   int N;
  72   int count = 0;
  73   double scal_par = 1.0;
  74   do {
  75     N = 0;
  76     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  77       if (bc_set[id_cell] == 0) {
  78         ++N;
  79         vec3_t n1 = cellNormal(m_Grid, id_cell);
  80         int bc = 9999;
  81         double scal_max = -3;
  82         vtkIdType id_copy = -1;
  83         for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
  84           vtkIdType id_neigh = m_Part.c2cGG(id_cell, i);
  85           vec3_t n2 = cellNormal(m_Grid, id_neigh);
  86           n1.normalise();
  87           n2.normalise();
  88           double scal = n1*n2;
  89           bool parallel = false;
  90           if (scal >= scal_par && cell_code->GetValue(id_neigh) < 9999) {
  91             scal = 1.0;
  92             parallel = true;
  93           }
  94           if (scal > scal_max || parallel) {
  95             bc = cell_code->GetValue(id_neigh);
  96             scal_max = scal;
  97           }
  98         }
  99         if (bc < 9999) {
 100           bc_set[id_cell] = bc;
 101           cell_code->SetValue(id_cell, bc);
 102         }
 103       }
 104     }
 105     scal_par *= 0.99;
 106     cout << scal_par << ',' << N << endl;
 107     ++count;
 108   } while (N > 0 && count < 1000);
 109 }
 110
 111
 112 /*
 113 void SharpenEdges::assignBCs()
 114 {
 115   QVector<int> bc_set(m_Grid->GetNumberOfCells(),0);
 116   QVector<int> bc_old(m_Grid->GetNumberOfCells(),0);
 117   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 118   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 119     bool certain_bc = true;
 120     vtkIdType N_pts, *pts;
 121     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 122     bc_old[id_cell] = cell_code->GetValue(id_cell);
 123     for (int i = 0; i < N_pts; ++i) {
 124       for (int j = 0; j < m_Part.n2cGSize(pts[i]); ++j) {
 125         if (cell_code->GetValue(id_cell) != cell_code->GetValue(m_Part.n2cGG(pts[i], j))) {
 126           certain_bc = false;
 127           break;
 128         }
 129       }
 130     }
 131     if (certain_bc) {
 132       bc_set[id_cell] = cell_code->GetValue(id_cell);
 133     }
 134   }
 135   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 136     if (bc_set[id_cell] == 0) {
 137       cell_code->SetValue(id_cell, 9999);
 138     }
 139   }
 140   GuiMainWindow::pointer()->storeSurfaceProjection();
 141
 142   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 143     if (cell_code->GetValue(id_cell) == 9999) {
 144       cell_code->SetValue(id_cell, bc_old[id_cell]);
 145     }
 146   }
 147
 148   int N;
 149   int count = 0;
 150   do {
 151     N = 0;
 152     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 153       bool ok = false;
 154       int bc;
 155       for (int i = 0; i < m_Part.c2cGSize(id_cell); ++i) {
 156         bc = cell_code->GetValue(m_Part.c2cGG(id_cell, i));
 157         if (bc == cell_code->GetValue(id_cell)) {
 158           ok = true;
 159           break;
 160         }
 161       }
 162       if (!ok) {
 163         cell_code->SetValue(id_cell, bc);
 164         ++N;
 165       }
 166     }
 167     ++count;
 168   } while (N > 0 && count < 1000);
 169 }
 170 */
 171
 172 void SharpenEdges::operate()
 173 {
 174   cout << "sharpening edges" << endl;
 175   assignBCs();
 176   return;
 177
 178
 179   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 180   QVector<QSet<int> > n2bc(m_Grid->GetNumberOfPoints());
 181   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 182     vtkIdType N_pts, *pts;
 183     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 184     int bc = cell_code->GetValue(id_cell);
 185     for (int i = 0; i < N_pts; ++i) {
 186       n2bc[pts[i]].insert(bc);
 187     }
 188   }
 189   int N = 0;
 190   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 191     if (n2bc[id_node].size() > 1) {
 192       ++N;
 193     }
 194   }
 195   int total_count = 0;
 196   int count = 0;
 197   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 198     if (n2bc[id_node].size() > 1) {
 199       ++count;
 200       ++total_count;
 201       vec3_t x;
 202       m_Grid->GetPoint(id_node, x.data());
 203       for (int i_proj_iter = 0; i_proj_iter < 20; ++i_proj_iter) {
 204         foreach (int bc, n2bc[id_node]) {
 205           x = GuiMainWindow::pointer()->getSurfProj(bc)->project(x, id_node);
 206         }
 207       }
 208       m_Grid->GetPoints()->SetPoint(id_node, x.data());
 209       if (count >= 100) {
 210         cout << total_count << " of " << N << " nodes corrected" << endl;
 211         count = 0;
 212       }
 213     }
 214   }
 215   GuiMainWindow::pointer()->storeSurfaceProjection();
 216
 217
 218 /*
 219   {
 220     LaplaceSmoother lap;
 221     lap.setGrid(m_Grid);
 222     QVector<vtkIdType> cls;
 223     m_BoundaryCodes = GuiMainWindow::pointer()->getAllBoundaryCodes();
 224     getSurfaceCells(m_BoundaryCodes, cls, m_Grid);
 225     lap.setCells(cls);
 226     lap.setNumberOfIterations(5);
 227     lap.setBoundaryCodes(m_BoundaryCodes);//IMPORTANT: so that unselected nodes become fixed when node types are updated!
 228     lap.setProjectionOn();
 229     lap.setNormalCorrectionOff();
 230     lap.setFreeProjectionForEdgesOn();
 231     lap();
 232   }
 233 */
 234
 235   UpdateCellIndex(m_Grid);
 236   GuiMainWindow::pointer()->updateBoundaryCodes(true);
 237   GuiMainWindow::pointer()->updateActors();
 238 }