src/libengrid/converttopolymesh.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 "converttopolymesh.h"
  22 #include "polymesh.h"
  23 #include "guimainwindow.h"
  24
  25 ConvertToPolyMesh::ConvertToPolyMesh()
  26 {
  27   m_Optimise = false;
  28   m_SplitCells = false;
  29   m_SplitFaces = false;
  30   m_PullInFactor = 0.0;
  31 }
  32
  33 void ConvertToPolyMesh::operate()
  34 {
  35   QList<VolumeDefinition> vols = mainWindow()->getAllVols();
  36   if (vols.size() > 1) {
  37     EG_ERR_RETURN("can only handle grids with a single volume at the moment");
  38   }
  39   vols.first().setVC(0);
  40   mainWindow()->setAllVols(vols);
  41
  42   PolyMesh pmesh(m_Grid, true, m_PullInFactor, m_Optimise, m_SplitFaces, m_SplitCells);
  43
  44   // count the number of boundary faces in the polygonal mesh
  45   int num_boundary_faces = 0;
  46   for (int i = 0; i < pmesh.numFaces(); ++i) {
  47     if (pmesh.boundaryCode(i) != 0) {
  48       if (pmesh.neighbour(i) != -1) {
  49         EG_BUG;
  50       }
  51       ++num_boundary_faces;
  52     }
  53   }
  54
  55   EG_VTKSP(vtkUnstructuredGrid, grid);
  56   EgVtkObject::allocateGrid(grid, pmesh.numPolyCells() + num_boundary_faces, pmesh.totalNumNodes());
  57
  58   // copy all nodes to the new grid
  59   for (vtkIdType id_node = 0; id_node < pmesh.totalNumNodes(); ++id_node) {
  60     vec3_t x = pmesh.nodeVector(id_node);
  61     grid->GetPoints()->SetPoint(id_node, x.data());
  62   }
  63
  64   // copy all polyhedral cells to the new grid
  65   EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
  66   for (int i = 0; i < pmesh.numCells(); ++i) {
  67     vtkIdType num_ids = 1;
  68     for (int j = 0; j < pmesh.numFacesOfPCell(i); ++j) {
  69       num_ids += 1 + pmesh.numNodes(pmesh.pcell2Face(i, j));
  70     }
  71     EG_VTKSP(vtkIdList, ptids);
  72     ptids->SetNumberOfIds(num_ids);
  73     vtkIdType i_id = 0;
  74     ptids->SetId(i_id++, pmesh.numFacesOfPCell(i));
  75     for (int j = 0; j < pmesh.numFacesOfPCell(i); ++j) {
  76       int i_face = pmesh.pcell2Face(i, j);
  77       ptids->SetId(i_id++, pmesh.numNodes(i_face));
  78       for (int k = 0; k < pmesh.numNodes(i_face); ++k) {
  79         if (pmesh.owner(i_face) == i) {
  80           ptids->SetId(i_id++, pmesh.nodeIndex(i_face, k));
  81         } else {
  82           ptids->SetId(i_id++, pmesh.nodeIndex(i_face, pmesh.numNodes(i_face) - 1 - k));
  83         }
  84       }
  85     }
  86     vtkIdType id_cell = grid->InsertNextCell(VTK_POLYHEDRON, ptids);
  87     cell_code->SetValue(id_cell, 1);
  88   }
  89
  90   // copy all boundary faces to the new grid
  91   for (int i = 0; i < pmesh.numFaces(); ++i) {
  92     if (pmesh.boundaryCode(i) != 0) {
  93       EG_VTKSP(vtkIdList, ptids);
  94       ptids->SetNumberOfIds(pmesh.numNodes(i));
  95       for (int j = 0; j < pmesh.numNodes(i); ++j) {
  96         ptids->SetId(j, pmesh.nodeIndex(i, j));
  97       }
  98       vtkIdType id_cell = grid->InsertNextCell(VTK_POLYGON, ptids);
  99       cell_code->SetValue(id_cell, pmesh.boundaryCode(i));
 100     }
 101   }
 102
 103   // compute surface integrals to detect open cells
 104   {
 105     MeshPartition part(grid, true);
 106     EG_VTKDCC(vtkDoubleArray, cell_mesh_quality, grid, "cell_mesh_quality");
 107     double max_err = 0;
 108     double min_dist =  1e99;
 109     double max_dist = -1e99;
 110     double average_dist = 0;
 111     int N = 0;
 112     for (vtkIdType id_cell = 0; id_cell < grid->GetNumberOfCells(); ++id_cell) {
 113       if (isVolume(id_cell, grid)) {
 114         vec3_t surf_int(0,0,0);
 115         double A = 0;
 116         vec3_t x = cellCentre(grid, id_cell);
 117         for (int i_face = 0; i_face < part.c2cGSize(id_cell); ++i_face) {
 118           vec3_t n = getNormalOfCell(grid, id_cell, i_face);
 119           surf_int += n;
 120           A += n.abs();
 121           vec3_t xf = getCentreOfCellFace(grid, id_cell, i_face);
 122           n.normalise();
 123           double dist = (xf - x)*n;
 124           max_dist = max(max_dist, dist);
 125           min_dist = min(min_dist, dist);
 126           average_dist += dist;
 127           ++N;
 128         }
 129         max_err = max(max_err, surf_int.abs()/A);
 130       }
 131     }
 132     average_dist /= N;
 133     cout << "maximal cell error: " << max_err << endl;
 134     cout << "maximal face distance: " << max_dist << endl;
 135     cout << "minimal face distance: " << min_dist << endl;
 136     cout << "average face distance: " << average_dist << endl;
 137   }
 138
 139   // replace the original mesh
 140   makeCopy(grid, m_Grid);
 141   m_Grid->Modified();
 142 }