src/libengrid/createvolumemesh.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
  22 #include "createvolumemesh.h"
  23 #include "deletetetras.h"
  24 #include "guimainwindow.h"
  25 #include "updatedesiredmeshdensity.h"
  26 #include <vtkXMLUnstructuredGridWriter.h>
  27
  28 CreateVolumeMesh::CreateVolumeMesh()
  29 {
  30   EG_TYPENAME;
  31 }
  32
  33 void CreateVolumeMesh::setTraceCells(const QVector<vtkIdType> &cells)
  34 {
  35   m_TraceCells.resize(cells.size());
  36   qCopy(cells.begin(), cells.end(), m_TraceCells.begin());
  37 }
  38
  39 void CreateVolumeMesh::getTraceCells(QVector<vtkIdType> &cells)
  40 {
  41   cells.resize(m_TraceCells.size());
  42   qCopy(m_TraceCells.begin(), m_TraceCells.end(), cells.begin());
  43 }
  44
  45 void CreateVolumeMesh::computeMeshDensity()
  46 {
  47   /*
  48   QString buffer = GuiMainWindow::pointer()->getXmlSection("engrid/surface/settings").replace("\n", " ");
  49   if (!buffer.isEmpty()) {
  50     QTextStream in(&buffer, QIODevice::ReadOnly);
  51     in >> m_MaxEdgeLength;
  52     in >> m_MinEdgeLength;
  53     in >> m_GrowthFactor;
  54   } else {
  55     m_MaxEdgeLength = 1000.0;
  56     m_MinEdgeLength = 0.0;
  57     m_GrowthFactor = 1.5;
  58   }
  59   m_ELSManager.read();
  60   QVector<double> H(m_Grid->GetNumberOfPoints(), m_MaxEdgeLength);
  61
  62   QVector<bool> fixed(m_Grid->GetNumberOfPoints(), false);
  63   double H_min = 1e99;
  64   vtkIdType id_min = -1;
  65   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
  66     bool volume_only = true;
  67     for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
  68       if (isSurface(m_Part.n2cGG(id_node, i), m_Grid)) {
  69         volume_only = false;
  70       }
  71       if (!volume_only) {
  72         fixed[id_node] = true;
  73         H[id_node] = 0;
  74         int N = 0;
  75         vec3_t xi;
  76         m_Grid->GetPoint(id_node, xi.data());
  77         for (int j = 0; j < m_Part.n2nGSize(id_node); ++j) {
  78           if (m_Part.n2nGG(id_node, j)) {
  79             vec3_t xj;
  80             m_Grid->GetPoint(m_Part.n2nGG(id_node, j), xj.data());
  81             H[id_node] += (xi-xj).abs();
  82             ++N;
  83           }
  84         }
  85         if (N < 2) {
  86           EG_BUG;
  87         }
  88         H[id_node] /= N;
  89                 if (H[id_node] < 0) {
  90               EG_BUG;
  91                 }
  92         if (H[id_node] < H_min) {
  93           id_min = id_node;
  94           H_min = H[id_node];
  95         }
  96       }
  97     }
  98   }
  99   if (id_min < 0) {
 100     EG_BUG;
 101   }
 102
 103   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 104     vec3_t x;
 105     m_Grid->GetPoint(id_node, x.data());
 106     double cl_src = m_ELSManager.minEdgeLength(x);
 107     if (cl_src > 0) {
 108       if (cl_src < H[id_node]) {
 109         H[id_node] = cl_src;
 110       }
 111     }
 112   }
 113
 114   QVector<bool> marked(m_Grid->GetNumberOfPoints(), false);
 115   marked[id_min] = true;
 116   bool done = false;
 117   while (!done) {
 118     done = true;
 119     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 120       if (marked[id_node] && H[id_node] <= H_min) {
 121         vec3_t x1;
 122         m_Grid->GetPoint(id_node, x1.data());
 123         for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 124           vtkIdType id_neigh = m_Part.n2nGG(id_node,i);
 125           if (!marked[id_neigh]) {
 126             vec3_t x2;
 127             m_Grid->GetPoint(id_neigh, x2.data());
 128             double dist = (x1 - x2).abs();
 129             double h = H[id_node] + (m_GrowthFactor - 1)*dist;
 130             if (h < 0) {
 131               EG_BUG;
 132             }
 133             H[id_neigh] = min(H[id_neigh], h);
 134             marked[id_neigh] = true;
 135             //H[id_neigh] += 1.0*H[id_node];
 136             done = false;
 137           }
 138         }
 139       }
 140     }
 141     H_min *= m_GrowthFactor;
 142   }
 143
 144   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 145     vec3_t x;
 146     m_Grid->GetPoint(id_node, x.data());
 147     double cl_src = m_ELSManager.minEdgeLength(x);
 148     if (cl_src > 0) {
 149       if (cl_src < H[id_node]) {
 150         H[id_node] = cl_src;
 151       }
 152     }
 153   }
 154
 155   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 156     vec3_t x1, x2;
 157     m_Grid->GetPoint(id_node, x1.data());
 158     x2 = x1;
 159     for (int j = 0; j < m_Part.n2nGSize(id_node); ++j) {
 160       vec3_t xj;
 161       m_Grid->GetPoint(m_Part.n2nGG(id_node, j), xj.data());
 162       for (int k = 0; k < 3; ++k) {
 163         x1[k] = min(xj[k], x1[k]);
 164         x2[k] = max(xj[k], x2[k]);
 165       }
 166     }
 167
 168     // Do something here!!
 169     EG_BUG;
 170
 171   }
 172   */
 173 }
 174
 175
 176 void CreateVolumeMesh::operate()
 177 {
 178   readSettings();
 179   double a = m_MaximalEdgeLength;
 180   double V = a*a*a/(6*sqrt(2.0));
 181   QString flags;
 182   flags.setNum(V);
 183   flags = QString("pq1.4a") + flags;
 184   cout << qPrintable(flags) << endl;
 185   tetgen(flags);
 186 }
 187