src/fixcadgeometry.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008,2009 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 "fixcadgeometry.h"
  24 #include "surfacemesher.h"
  25 #include "vertexmeshdensity.h"
  26 #include "guimainwindow.h"
  27 #include "deletestraynodes.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_FeatureAngle = GeometryTools::deg2rad(0.5);
  45   m_AllowFeatureEdgeSwapping = false;
  46   m_AllowSmallAreaSwapping = true;
  47   m_NumDelaunaySweeps = 100;
  48 }
  49
  50 void FixCadGeometry::operate()
  51 {
  52   {
  53     DeleteStrayNodes del;
  54     del();
  55   }
  56
  57   setAllCells();
  58
  59   //prepare BCmap
  60   EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
  61   QSet <int> bcs;
  62   l2g_t cells = m_Part.getCells();
  63   foreach (vtkIdType id_cell, cells) {
  64     if (isSurface(id_cell, m_Grid)) {
  65       bcs.insert(bc->GetValue(id_cell));
  66     }
  67   }
  68   QMap <int,int> BCmap;
  69   foreach(int bc, bcs) {
  70     BCmap[bc] = 1;
  71   }
  72
  73   //set density infinite
  74   VertexMeshDensity VMD;
  75   VMD.density = 1e99;
  76   VMD.BCmap = BCmap;
  77   cout<<"VMD="<<VMD<<endl;
  78   qWarning()<<"VMD.BCmap="<<VMD.BCmap;
  79   m_VMDvector.push_back(VMD);
  80   cout<<"m_VMDvector="<<m_VMDvector<<endl;
  81
  82   customUpdateNodeInfo();
  83
  84   //call surface mesher
  85   setGrid(m_Grid);
  86   setBoundaryCodes(bcs);
  87   setVertexMeshDensityVector(m_VMDvector);
  88   setDesiredLength();
  89   mesher();
  90
  91   // finalise
  92   createIndices(m_Grid);
  93   customUpdateNodeInfo();
  94   setDesiredLength();
  95 }
  96
  97 void FixCadGeometry::mesher()
  98 {
  99   setDesiredLength();
 100   if (m_BoundaryCodes.size() == 0) {
 101     return;
 102   }
 103   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired");
 104   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 105     characteristic_length_desired->SetValue(id_node, 1e-6);
 106   }
 107   int num_deleted = 0;
 108   int iter = 0;
 109   bool done = false;
 110   swap();
 111   while (!done) {
 112     ++iter;
 113     cout << "CAD fix iteration " << iter << ":" << endl;
 114     setDesiredLength();
 115     customUpdateNodeInfo();
 116     setDesiredLength();
 117     int num_deleted = deleteNodes();
 118     cout << "  deleted nodes  : " << num_deleted << endl;
 119     swap();
 120     setDesiredLength();
 121     done = (iter >= 20) || (num_deleted == 0);
 122     cout << "  total nodes    : " << m_Grid->GetNumberOfPoints() << endl;
 123     cout << "  total cells    : " << m_Grid->GetNumberOfCells() << endl;
 124   }
 125   createIndices(m_Grid);
 126   customUpdateNodeInfo();
 127   setDesiredLength();
 128 }
 129
 130 void FixCadGeometry::setDesiredLength(double L)
 131 {
 132   setAllSurfaceCells();
 133   l2g_t  nodes = getPartNodes();
 134   g2l_t _nodes = getPartLocalNodes();
 135   l2g_t  cells = getPartCells();
 136   l2l_t  n2n   = getPartN2N();
 137   l2l_t  c2c   = getPartC2C();
 138
 139   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired,   m_Grid, "node_meshdensity_desired");
 140   EG_VTKDCN(vtkIntArray,    characteristic_length_specified, m_Grid, "node_specified_density");
 141
 142   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 143     vtkIdType id_node = nodes[i_nodes];
 144     characteristic_length_specified->SetValue(id_node, 0);
 145     characteristic_length_desired->SetValue(id_node, L);
 146   }
 147 }
 148
 149 void FixCadGeometry::customUpdateNodeInfo()
 150 {
 151   setAllCells();
 152   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 153   l2g_t cells = getPartCells();
 154   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 155     node_type->SetValue(id_node, VTK_FIXED_VERTEX);
 156   }
 157   foreach (vtkIdType id_cell, cells) {
 158     if (isSurface(id_cell, m_Grid)) {
 159       vtkIdType N_pts, *pts;
 160       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 161       if (N_pts == 3) {
 162         vec3_t n1 = cellNormal(m_Grid, id_cell);
 163         QVector<int> num_bad_edges(3, 0);
 164         for (int i = 0; i < N_pts; ++i) {
 165           int i1 = i;
 166           int i2 = 0;
 167           if (i < N_pts - 1) {
 168             i2= i+1;
 169           }
 170           vec3_t x1, x2;
 171           m_Grid->GetPoint(pts[i1], x1.data());
 172           m_Grid->GetPoint(pts[i2], x2.data());
 173           double L = (x1 - x2).abs();
 174           vtkIdType id_neigh_cell = m_Part.c2cGG(id_cell, i);
 175           if (id_neigh_cell != -1) {
 176             bool bad_edge = false;
 177             vec3_t n2 = cellNormal(m_Grid, id_neigh_cell);
 178             if (angle(n1, n2) > deg2rad(179.5)) {
 179               bad_edge = true;
 180             }
 181             if (bad_edge) {
 182               ++num_bad_edges[i1];
 183               ++num_bad_edges[i2];
 184             }
 185           }
 186         }
 187         for (int i = 0; i < N_pts; ++i) {
 188           if (num_bad_edges[i] >= 2) {
 189             node_type->SetValue(pts[i], VTK_SIMPLE_VERTEX);
 190           }
 191         }
 192       }
 193     }
 194   }
 195 }
 196