src/libengrid/updatedesiredmeshdensity.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 "updatedesiredmeshdensity.h"
  22 #include "guimainwindow.h"
  23 #include "pointfinder.h"
  24 #include "cgaltricadinterface.h"
  25
  26 #include <vtkCharArray.h>
  27
  28 UpdateDesiredMeshDensity::UpdateDesiredMeshDensity() : SurfaceOperation()
  29 {
  30   EG_TYPENAME;
  31   m_MaxEdgeLength = 1e99;
  32   m_NodesPerQuarterCircle = 0;
  33   m_OnlySurfaceCells = true;
  34
  35   m_GrowthFactor = 0.0;
  36   m_MinEdgeLength = 0.0;
  37   m_FeatureResolution2D = 0;
  38   m_FeatureResolution3D = 0;
  39
  40   getSet("surface meshing", "threshold angle for feature resolution", 20, m_FeatureThresholdAngle);
  41   m_FeatureThresholdAngle = deg2rad(m_FeatureThresholdAngle);
  42   m_Relaxation = true;
  43
  44 }
  45
  46 double UpdateDesiredMeshDensity::computeSearchDistance(vtkIdType id_face)
  47 {
  48   vtkIdType N_pts, *pts;
  49   m_Grid->GetCellPoints(id_face, N_pts, pts);
  50   QVector<vec3_t> x(N_pts + 1);
  51   for (int i = 0; i < N_pts; ++i) {
  52     m_Grid->GetPoint(pts[i], x[i].data());
  53   }
  54   x[N_pts] = x[0];
  55   double L = 0;
  56   for (int i = 0; i < N_pts; ++i) {
  57     L = max(L, (x[i] - x[i+1]).abs());
  58   }
  59   return L;
  60 }
  61
  62 void UpdateDesiredMeshDensity::computeExistingLengths()
  63 {
  64   QSet<int> all_bcs = GuiMainWindow::pointer()->getAllBoundaryCodes();
  65   QSet<int> fixed_bcs = all_bcs - m_BoundaryCodes;
  66   QVector<double> edge_length(m_Grid->GetNumberOfPoints(), 1e99);
  67   QVector<int> edge_count(m_Grid->GetNumberOfPoints(), 0);
  68   m_Fixed.fill(false, m_Grid->GetNumberOfPoints());
  69   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  70   for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
  71     if (isSurface(id_cell, m_Grid)) {
  72       if (fixed_bcs.contains(cell_code->GetValue(id_cell))) {
  73         vtkIdType N_pts, *pts;
  74         m_Grid->GetCellPoints(id_cell, N_pts, pts);
  75         QVector<vec3_t> x(N_pts);
  76         for (int i = 0; i < N_pts; ++i) {
  77           m_Grid->GetPoint(pts[i], x[i].data());
  78           m_Fixed[pts[i]] = true;
  79         }
  80         for (int i = 0; i < N_pts; ++i) {
  81           int j = i + 1;
  82           if (j >= N_pts) {
  83             j = 0;
  84           }
  85           double L = (x[i] - x[j]).abs();
  86           edge_length[pts[i]] = min(edge_length[pts[i]], L);
  87           edge_length[pts[j]] = min(edge_length[pts[j]], L);
  88           ++edge_count[pts[i]];
  89           ++edge_count[pts[j]];
  90         }
  91       }
  92     }
  93   }
  94   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired,   m_Grid, "node_meshdensity_desired");
  95   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
  96     if (edge_count[id_node] > 0) {
  97       if (edge_length[id_node] > 1e98) {
  98         EG_BUG;
  99       }
 100       characteristic_length_desired->SetValue(id_node, edge_length[id_node]);
 101     }
 102   }
 103 }
 104
 105 void UpdateDesiredMeshDensity::computeFeature(const QList<point_t> points, QVector<double> &cl_pre, double res, int restriction_type)
 106 {
 107   int N = 0;
 108   QVector<vec3_t> pts(points.size());
 109   for (int i = 0; i < points.size(); ++i) {
 110     pts[i] = points[i].x;
 111   }
 112   PointFinder pfind;
 113   pfind.setMaxNumPoints(5000);
 114   pfind.setPoints(pts);
 115   for (int i = 0; i < points.size(); ++i) {
 116     double h = -1;
 117     foreach (int i_points, points[i].idx) {
 118       h = max(h, cl_pre[i_points]);
 119     }
 120     if (h < 0) {
 121       h = 1e99;
 122     }
 123     QVector<int> close_points;
 124     pfind.getClosePoints(points[i].x, close_points, res*points[i].L);
 125     foreach (int j, close_points) {
 126       if (i != j) {
 127         ++N;
 128         vec3_t x1 = points[i].x;
 129         vec3_t x2 = points[j].x;
 130         vec3_t n1 = points[i].n;
 131         vec3_t n2 = points[j].n;
 132         vec3_t v = x2 - x1;
 133         if (n1*n2 < 0) {
 134           if (n1*v > 0) {
 135             if (fabs(GeometryTools::angle(n1, (-1)*n2)) <= m_FeatureThresholdAngle) {
 136               double l = v.abs()/fabs(n1*n2);
 137               if (l/res < h) {
 138                 // check topological distance
 139                 int search_dist = int(ceil(res));
 140                 bool topo_dist_ok = true;
 141                 foreach (int k, points[i].idx) {
 142                   vtkIdType id_node1 = m_Part.globalNode(k);
 143                   foreach (int l, points[j].idx) {
 144                     vtkIdType id_node2 = m_Part.globalNode(l);
 145                     if (m_Part.computeTopoDistance(id_node1, id_node2, search_dist, restriction_type) < search_dist) {
 146                       topo_dist_ok = false;
 147                       break;
 148                     }
 149                     if (!topo_dist_ok) {
 150                       break;
 151                     }
 152                   }
 153                 }
 154                 if (topo_dist_ok) {
 155                   h = l/res;
 156                 }
 157               }
 158             }
 159           }
 160         }
 161       }
 162     }
 163     foreach (int i_points, points[i].idx) {
 164       cl_pre[i_points] = min(h, cl_pre[i_points]);
 165     }
 166   }
 167 }
 168
 169 void UpdateDesiredMeshDensity::computeFeature2D(QVector<double> &cl_pre)
 170 {
 171   if (m_FeatureResolution2D < 1e-3) {
 172     return;
 173   }
 174   QSet<int> bcs = getAllBoundaryCodes(m_Grid);
 175   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 176   foreach (int bc, bcs) {
 177     QList<point_t> points;
 178     for (vtkIdType id_face = 0; id_face < m_Grid->GetNumberOfCells(); ++id_face) {
 179       if (isSurface(id_face, m_Grid)) {
 180         if (cell_code->GetValue(id_face) == bc) {
 181           vec3_t xc = cellCentre(m_Grid, id_face);
 182           vtkIdType num_pts, *pts;
 183           m_Grid->GetCellPoints(id_face, num_pts, pts);
 184           QVector<vec3_t> xn(num_pts + 1);
 185           QVector<vtkIdType> idx(num_pts + 1);
 186           for (int i_pts = 0; i_pts < num_pts; ++i_pts) {
 187             m_Grid->GetPoint(pts[i_pts], xn[i_pts].data());
 188             idx[i_pts] = m_Part.localNode(pts[i_pts]);
 189           }
 190           xn[num_pts] = xn[0];
 191           idx[num_pts] = idx[0];
 192           for (int i_neigh = 0; i_neigh < m_Part.c2cGSize(id_face); ++i_neigh) {
 193             vtkIdType id_neigh = m_Part.c2cGG(id_face, i_neigh);
 194             if (id_neigh != -1) {
 195               if (cell_code->GetValue(id_neigh) != bc) {
 196                 point_t P;
 197                 P.x = 0.5*(xn[i_neigh] + xn[i_neigh + 1]);
 198                 P.n = xc - xn[i_neigh];
 199                 vec3_t v = xn[i_neigh + 1] - xn[i_neigh];
 200                 v.normalise();
 201                 P.n -= (P.n*v)*v;
 202                 P.n.normalise();
 203                 P.idx.append(idx[i_neigh]);
 204                 P.idx.append(idx[i_neigh+1]);
 205                 P.L = computeSearchDistance(id_face);
 206                 points.append(P);
 207               }
 208             }
 209           }
 210         }
 211       }
 212     }
 213     computeFeature(points, cl_pre, m_FeatureResolution2D, 2);
 214   }
 215 }
 216
 217 void UpdateDesiredMeshDensity::computeFeature3D(QVector<double> &cl_pre)
 218 {
 219   // add mesh density radiation to 3D feature resolution
 220   EG_STOPDATE("2015-06-01");
 221
 222   if (m_FeatureResolution3D < 1e-3) {
 223     return;
 224   }
 225   CgalTriCadInterface cad(m_Grid);
 226
 227   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 228     vec3_t n0 = m_Part.globalNormal(id_node);
 229     vec3_t x0;
 230     m_Grid->GetPoint(id_node, x0.data());
 231     QVector<QPair<vec3_t, vtkIdType> > inters1;
 232     cad.computeIntersections(x0, n0, inters1);
 233     QVector<QPair<vec3_t, vtkIdType> > inters2;
 234     cad.computeIntersections(x0, (-1)*n0, inters2);
 235     QVector<QPair<vec3_t, vtkIdType> > intersections = inters1 + inters2;
 236     double d_min = EG_LARGE_REAL;
 237     vtkIdType id_tri_min = -1;
 238     for (int i = 0; i < intersections.size(); ++i) {
 239       vec3_t x = intersections[i].first;
 240       vtkIdType id_tri = intersections[i].second;
 241       vec3_t n = GeometryTools::cellNormal(m_Grid, id_tri);
 242       n.normalise();
 243       vec3_t dx = x - x0;
 244       if (dx*n0 < 0 && n*n0 < 0) {
 245         double d = dx.abs();
 246         if (d < d_min) {
 247           d_min = d;
 248           id_tri_min = id_tri;
 249         }
 250       }
 251     }
 252     if (id_tri_min > -1) {
 253       double h = d_min/m_FeatureResolution3D;
 254       int max_topo_dist = int(ceil(m_FeatureResolution3D));
 255       bool topo_dist_ok = true;
 256       EG_GET_CELL(id_tri_min, m_Grid);
 257       for (int i = 0; i < num_pts; ++i) {
 258         int topo_dist = m_Part.computeTopoDistance(pts[i], id_node, max_topo_dist, 0);
 259         if (topo_dist < max_topo_dist) {
 260           topo_dist_ok = false;
 261           break;
 262         }
 263       }
 264       if (topo_dist_ok) {
 265         cl_pre[id_node] = h;
 266         for (int i = 0; i < num_pts; ++i) {
 267           cl_pre[pts[i]] = h;
 268         }
 269       }
 270     }
 271   }
 272 }
 273
 274 void UpdateDesiredMeshDensity::readSettings()
 275 {
 276   readVMD();
 277   QString buffer = GuiMainWindow::pointer()->getXmlSection("engrid/surface/settings").replace("\n", " ");
 278   QTextStream in(&buffer, QIODevice::ReadOnly);
 279   in >> m_MaxEdgeLength;
 280   in >> m_MinEdgeLength;
 281   in >> m_GrowthFactor;
 282   in >> m_NodesPerQuarterCircle;
 283   int num_bcs;
 284   in >> num_bcs;
 285   QVector<int> tmp_bcs;
 286   GuiMainWindow::pointer()->getAllBoundaryCodes(tmp_bcs);
 287   m_BoundaryCodes.clear();
 288   if (num_bcs == tmp_bcs.size()) {
 289     foreach (int bc, tmp_bcs) {
 290       int check_state;
 291       in >> check_state;
 292       if (check_state == 1) {
 293         m_BoundaryCodes.insert(bc);
 294       }
 295     }
 296   }
 297   if (!in.atEnd()) {
 298     in >> m_FeatureResolution2D;
 299     in >> m_FeatureResolution3D;
 300   }
 301 }
 302
 303 void UpdateDesiredMeshDensity::operate()
 304 {
 305   m_ELSManager.read();
 306   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 307
 308   if (m_OnlySurfaceCells) {
 309     setAllSurfaceCells();
 310   } else {
 311     setAllCells();
 312   }
 313   l2g_t  nodes = getPartNodes();
 314   l2g_t  cells = getPartCells();
 315   l2l_t  n2n   = getPartN2N();
 316
 317   EG_VTKDCN(vtkDoubleArray, characteristic_length_desired,   m_Grid, "node_meshdensity_desired");
 318   EG_VTKDCN(vtkIntArray,    characteristic_length_specified, m_Grid, "node_specified_density");
 319
 320   QVector<vec3_t> normals(cells.size(), vec3_t(0,0,0));
 321   QVector<vec3_t> centres(cells.size(), vec3_t(0,0,0));
 322   QVector<double> cl_pre(nodes.size(), 1e99);
 323
 324   computeExistingLengths();
 325   if (m_BoundaryCodes.size() == 0) {
 326     return;
 327   }
 328
 329   if (m_NodesPerQuarterCircle > 1e-3) {
 330     QVector<double> R(nodes.size(), 1e99);
 331
 332 //#pragma omp parallel for
 333     for (int i_cell = 0; i_cell < cells.size(); ++i_cell) {
 334       vtkIdType id_cell = cells[i_cell];
 335       vtkIdType N_pts, *pts;
 336       m_Grid->GetCellPoints(id_cell, N_pts, pts);
 337       int bc = cell_code->GetValue(id_cell);
 338       for (int i = 0; i < N_pts; ++i) {
 339         int i_nodes = m_Part.localNode(pts[i]);
 340         R[i_nodes] = min(R[i_nodes], fabs(GuiMainWindow::pointer()->getCadInterface(bc)->getRadius(pts[i])));
 341       }
 342     }
 343
 344     for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 345       if (cl_pre[i_nodes] == 0) {
 346         EG_BUG;
 347       }
 348       cl_pre[i_nodes] = max(m_MinEdgeLength, min(cl_pre[i_nodes], 0.5*R[i_nodes]*M_PI/m_NodesPerQuarterCircle));
 349     }
 350   }
 351
 352   // cells across branches
 353   computeFeature2D(cl_pre);
 354   computeFeature3D(cl_pre);
 355
 356   // set everything to desired mesh density and find maximal mesh-density
 357   double cl_min = 1e99;
 358   int i_nodes_min = -1;
 359   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 360     vtkIdType id_node = nodes[i_nodes];
 361     double cl = m_MaxEdgeLength;
 362     if (m_BoundaryCodes.size() > 0) {
 363       int idx = characteristic_length_specified->GetValue(id_node);
 364       if ((idx >= 0) && (idx < m_VMDvector.size())) {
 365         cl = m_VMDvector[idx].density;
 366       }
 367     }
 368     if (m_Fixed[id_node]) {
 369       cl = characteristic_length_desired->GetValue(id_node);
 370     }
 371     cl = min(cl_pre[i_nodes], cl);
 372     vec3_t x;
 373     m_Grid->GetPoint(id_node, x.data());
 374     double cl_src = m_ELSManager.minEdgeLength(x);
 375     if (cl_src > 0) {
 376       cl = min(cl, cl_src);
 377     }
 378
 379     cl = max(m_MinEdgeLength, cl);
 380
 381     if(cl == 0) {
 382       EG_BUG;
 383     }
 384     characteristic_length_desired->SetValue(id_node, cl);
 385
 386     if (cl < cl_min) {
 387       for (int i = 0; i < m_Part.n2bcGSize(id_node); ++i) {
 388         if (m_BoundaryCodes.contains(m_Part.n2bcG(id_node, i))) {
 389           cl_min = cl;
 390           i_nodes_min = i_nodes;
 391           break;
 392         }
 393       }
 394     }
 395   }
 396   if (i_nodes_min == -1) {
 397     EG_ERR_RETURN("There are no edges that need improving.")
 398   }
 399
 400
 401   if (m_Relaxation) {
 402
 403     // start from smallest characteristic length and loop as long as nodes are updated
 404     int num_updated = 0;
 405
 406     do {
 407       num_updated = 0;
 408       for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 409         double cli = characteristic_length_desired->GetValue(nodes[i_nodes]);
 410         if (cli <= cl_min) {
 411           vec3_t xi;
 412           m_Grid->GetPoint(nodes[i_nodes], xi.data());
 413           for (int j = 0; j < n2n[i_nodes].size(); ++j) {
 414             int j_nodes = n2n[i_nodes][j];
 415             double clj = characteristic_length_desired->GetValue(nodes[j_nodes]);
 416             if (clj > cli && clj > cl_min) {
 417               vec3_t xj;
 418               m_Grid->GetPoint(nodes[j_nodes], xj.data());
 419               ++num_updated;
 420               double L_new = min(m_MaxEdgeLength, cli * m_GrowthFactor);
 421               if (!m_Fixed[nodes[j_nodes]]) {
 422
 423                 double cl_min = min(characteristic_length_desired->GetValue(nodes[j_nodes]), L_new);
 424                 if(cl_min==0) {
 425                   qWarning()<<"m_MaxEdgeLength="<<m_MaxEdgeLength;
 426                   qWarning()<<"cli="<<cli;
 427                   qWarning()<<"m_GrowthFactor="<<m_GrowthFactor;
 428                   qWarning()<<"characteristic_length_desired->GetValue(nodes[j_nodes])="<<characteristic_length_desired->GetValue(nodes[j_nodes]);
 429                   qWarning()<<"L_new="<<L_new;
 430                   EG_BUG;
 431                 }
 432                 characteristic_length_desired->SetValue(nodes[j_nodes], cl_min);
 433
 434               }
 435             }
 436           }
 437         }
 438       }
 439       cl_min *= m_GrowthFactor;
 440     } while (num_updated > 0);
 441   }
 442 }