src/updatedesiredmeshdensity.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 "updatedesiredmeshdensity.h"
  24 #include "guimainwindow.h"
  25
  26 #include <vtkCharArray.h>
  27
  28 UpdateDesiredMeshDensity::UpdateDesiredMeshDensity() : SurfaceOperation()
  29 {
  30   EG_TYPENAME;
  31   m_MaxEdgeLength = 1e99;
  32   m_NodesPerQuarterCircle = 0;
  33 }
  34
  35
  36 void UpdateDesiredMeshDensity::computeExistingLengths()
  37 {
  38   QSet<int> all_bcs;
  39   GuiMainWindow::pointer()->getAllBoundaryCodes(all_bcs);
  40   QSet<int> fixed_bcs = all_bcs - m_BoundaryCodes;
  41   QVector<double> edge_length(m_Grid->GetNumberOfPoints(), 1e99);
  42   QVector<int> edge_count(m_Grid->GetNumberOfPoints(), 0);
  43   m_Fixed.fill(false, m_Grid->GetNumberOfPoints());
  44   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  45   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  46     if (isSurface(id_cell, m_Grid)) {
  47       if (fixed_bcs.contains(cell_code->GetValue(id_cell))) {
  48         vtkIdType N_pts, *pts;
  49         m_Grid->GetCellPoints(id_cell, N_pts, pts);
  50         vec3_t x[N_pts];
  51         for (int i = 0; i < N_pts; ++i) {
  52           m_Grid->GetPoint(pts[i], x[i].data());
  53           m_Fixed[pts[i]] = true;
  54         }
  55         for (int i = 0; i < N_pts; ++i) {
  56           int j = i + 1;
  57           if (j >= N_pts) {
  58             j = 0;
  59           }
  60           double L = (x[i] - x[j]).abs();
  61           edge_length[pts[i]] = min(edge_length[pts[i]], L);
  62           edge_length[pts[j]] = min(edge_length[pts[j]], L);
  63           ++edge_count[pts[i]];
  64           ++edge_count[pts[j]];
  65         }
  66       }
  67     }
  68   }
  69   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired,   m_Grid, "node_meshdensity_desired");
  70   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
  71     if (edge_count[id_node] > 0) {
  72       if (edge_length[id_node] > 1e98) {
  73         EG_BUG;
  74       }
  75       characteristic_length_desired->SetValue(id_node, edge_length[id_node]);
  76     }
  77   }
  78 }
  79
  80 void UpdateDesiredMeshDensity::operate()
  81 {
  82   m_ELSManager.read();
  83   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  84
  85   setAllSurfaceCells();
  86   l2g_t  nodes = getPartNodes();
  87   g2l_t _nodes = getPartLocalNodes();
  88   l2g_t  cells = getPartCells();
  89   l2l_t  n2n   = getPartN2N();
  90   l2l_t  c2c   = getPartC2C();
  91
  92   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired,   m_Grid, "node_meshdensity_desired");
  93   EG_VTKDCN(vtkIntArray,    characteristic_length_specified, m_Grid, "node_specified_density");
  94
  95   QVector<vec3_t> normals(cells.size(), vec3_t(0,0,0));
  96   QVector<vec3_t> centres(cells.size(), vec3_t(0,0,0));
  97   QVector<double> cl_radius(nodes.size(), 1e99);
  98
  99   computeExistingLengths();
 100   if (m_BoundaryCodes.size() == 0) {
 101     return;
 102   }
 103
 104   if (m_NodesPerQuarterCircle > 1e-3) {
 105
 106     // compute node normals
 107     for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
 108       normals[i_cells] = GeometryTools::cellNormal(m_Grid, cells[i_cells]);
 109       normals[i_cells].normalise();
 110       centres[i_cells] = cellCentre(m_Grid, cells[i_cells]);
 111     }
 112
 113     // compute characteristic length according to nodes per quarter circle
 114     for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
 115       vec3_t xi = centres[i_cells];
 116       vtkIdType N_pts, *pts;
 117       m_Grid->GetCellPoints(cells[i_cells], N_pts, pts);
 118       for (int j = 0; j < c2c[i_cells].size(); ++j) {
 119         int j_cells = c2c[i_cells][j];
 120         if (cell_code->GetValue(cells[i_cells]) == cell_code->GetValue(cells[j_cells])) {
 121           vec3_t xj = centres[j_cells];
 122           double cosa  = normals[i_cells]*normals[j_cells];
 123           double alpha = acos(cosa);
 124           if (alpha > 0.01*M_PI) {
 125             vec3_t va = xi - xj;
 126             vec3_t n1 = normals[i_cells] + normals[j_cells];
 127             vec3_t n2 = GeometryTools::orthogonalVector(n1);
 128             n2.normalise();
 129             va -= (va*n2)*n2;
 130             double a  = va.abs();
 131             double R  = 0.5*a/sin(alpha);
 132             double cl = 0.5*R*M_PI/m_NodesPerQuarterCircle;
 133             for (int k = 0; k < N_pts; ++k) {
 134               cl_radius[_nodes[pts[k]]] = min(cl_radius[_nodes[pts[k]]], cl);
 135             }
 136           }
 137         }
 138       }
 139     }
 140   }
 141
 142   // set everything to desired mesh density and find maximal mesh-density
 143   double cl_min = 1e99;
 144   int i_nodes_min = -1;
 145   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 146     vtkIdType id_node = nodes[i_nodes];
 147     double cl = 1e99;
 148     if (m_BoundaryCodes.size() > 0) {
 149       int idx = characteristic_length_specified->GetValue(id_node);
 150       if (idx != -1) {
 151         if (idx >= m_VMDvector.size()) {
 152           qWarning()<<"idx="<<idx;
 153           qWarning()<<"m_VMDvector.size()="<<m_VMDvector.size();
 154           EG_BUG;
 155         }
 156         cl = m_VMDvector[idx].density;
 157       }
 158     }
 159     if (m_Fixed[id_node]) {
 160       cl = characteristic_length_desired->GetValue(id_node);
 161     }
 162     cl = min(cl_radius[i_nodes], cl);
 163     vec3_t x;
 164     m_Grid->GetPoint(id_node, x.data());
 165     double cl_src = m_ELSManager.minEdgeLength(x);
 166     if (cl_src > 0) {
 167       cl = min(cl, cl_src);
 168     }
 169
 170     if(cl == 0) EG_BUG;
 171     characteristic_length_desired->SetValue(id_node, cl);
 172
 173     if (cl < cl_min) {
 174       cl_min = cl;
 175       i_nodes_min = i_nodes;
 176     }
 177   }
 178   if (i_nodes_min == -1) {
 179     EG_BUG;
 180   }
 181
 182   // start from smallest characteristic length and loop as long as nodes are updated
 183   int num_updated = 0;
 184
 185   do {
 186     num_updated = 0;
 187     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 188       double cli = characteristic_length_desired->GetValue(nodes[i_nodes]);
 189       if (cli <= cl_min) {
 190         vec3_t xi;
 191         m_Grid->GetPoint(nodes[i_nodes], xi.data());
 192         for (int j = 0; j < n2n[i_nodes].size(); ++j) {
 193           int j_nodes = n2n[i_nodes][j];
 194           double clj = characteristic_length_desired->GetValue(nodes[j_nodes]);
 195           if (clj > cli && clj > cl_min) {
 196             vec3_t xj;
 197             m_Grid->GetPoint(nodes[j_nodes], xj.data());
 198             ++num_updated;
 199             double L_new = min(m_MaxEdgeLength, cli * m_GrowthFactor);
 200             if (!m_Fixed[nodes[j_nodes]]) {
 201
 202               double cl_min = min(characteristic_length_desired->GetValue(nodes[j_nodes]), L_new);
 203               if(cl_min==0) {
 204                 qWarning()<<"m_MaxEdgeLength="<<m_MaxEdgeLength;
 205                 qWarning()<<"cli="<<cli;
 206                 qWarning()<<"m_GrowthFactor="<<m_GrowthFactor;
 207                 qWarning()<<"characteristic_length_desired->GetValue(nodes[j_nodes])="<<characteristic_length_desired->GetValue(nodes[j_nodes]);
 208                 qWarning()<<"L_new="<<L_new;
 209                 EG_BUG;
 210               }
 211               characteristic_length_desired->SetValue(nodes[j_nodes], cl_min);
 212
 213             }
 214           }
 215         }
 216       }
 217     }
 218     cl_min *= m_GrowthFactor;
 219   } while (num_updated > 0);
 220
 221   // do a simple averaging step
 222   /*
 223   QVector<double> cl_save(m_Grid->GetNumberOfPoints(), 0.0);
 224   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 225     cl_save[id_node] = characteristic_length_desired->GetValue(id_node);
 226   }
 227   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 228     double cl_new = 0;
 229     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 230       cl_new += cl_save[m_Part.n2nGG(id_node, i)];
 231     }
 232     if (m_Part.n2nGSize(id_node) > 0) {
 233       characteristic_length_desired->SetValue(id_node, cl_new/m_Part.n2nGSize(id_node));
 234     }
 235   }
 236   */
 237 }