src/createvolumemesh.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 "createvolumemesh.h"
  24 #include "deletetetras.h"
  25 #include "guimainwindow.h"
  26 #include <vtkXMLUnstructuredGridWriter.h>
  27
  28 CreateVolumeMesh::CreateVolumeMesh()
  29 {
  30   EG_TYPENAME;
  31   maxh     = 1e99;
  32   fineness = 0.0;
  33 }
  34
  35 void CreateVolumeMesh::setTraceCells(const QVector<vtkIdType> &cells)
  36 {
  37   trace_cells.resize(cells.size());
  38   qCopy(cells.begin(), cells.end(), trace_cells.begin());
  39 }
  40
  41 void CreateVolumeMesh::getTraceCells(QVector<vtkIdType> &cells)
  42 {
  43   cells.resize(trace_cells.size());
  44   qCopy(trace_cells.begin(), trace_cells.end(), cells.begin());
  45 }
  46
  47 void CreateVolumeMesh::prepare()
  48 {
  49   using namespace nglib;
  50   DeleteTetras del;
  51   del.setGrid(m_Grid);
  52   del.setAllCells();
  53   del();
  54   QVector<vtkIdType> cells, nodes;
  55   QVector<int>       _cells, _nodes;
  56   QVector<QVector< int > > c2c;
  57   QVector<QVector<int> > n2c;
  58   getAllCells(cells, m_Grid);
  59   createCellMapping(cells, _cells, m_Grid);
  60   getNodesFromCells(cells, nodes, m_Grid);
  61   createNodeMapping(nodes, _nodes, m_Grid);
  62   createCellToCell(cells, c2c, m_Grid);
  63   createNodeToCell(cells, nodes, _nodes, n2c, m_Grid);
  64   QList<QVector<vtkIdType> > ex_tri;
  65   int N1 = 0;
  66   int N2 = 0;
  67   int N3 = 0;
  68   int N4 = 0;
  69   foreach (vtkIdType id_cell, cells) {
  70     vtkIdType type_cell = m_Grid->GetCellType(id_cell);
  71     vtkIdType *pts, N_pts;
  72     m_Grid->GetCellPoints(id_cell, N_pts, pts);
  73     QVector<vtkIdType> T(3);
  74     if (type_cell == VTK_TRIANGLE) {
  75       if (findVolumeCell(m_Grid, id_cell, _nodes, cells, _cells, n2c) == -1) {
  76         T[0] = pts[0];
  77         T[1] = pts[1];
  78         T[2] = pts[2];
  79         ex_tri.append(T);
  80         ++N4;
  81       }
  82     } else if (type_cell == VTK_QUAD) {
  83       if (findVolumeCell(m_Grid, id_cell, _nodes, cells, _cells, n2c) == -1) {
  84         EG_BUG;
  85         T[0] = pts[0];
  86         T[1] = pts[1];
  87         T[2] = pts[2];
  88         ex_tri.append(T);
  89         T[0] = pts[2];
  90         T[1] = pts[3];
  91         T[2] = pts[0];
  92         ex_tri.append(T);
  93         ++N3;
  94       }
  95     } else if (type_cell == VTK_WEDGE) {
  96       if (c2c[id_cell][0] == -1) {
  97         T[0] = pts[0];
  98         T[1] = pts[2];
  99         T[2] = pts[1];
 100         ex_tri.append(T);
 101         ++N2;
 102       }
 103       if (c2c[id_cell][1] == -1) {
 104         T[0] = pts[3];
 105         T[1] = pts[4];
 106         T[2] = pts[5];
 107         ex_tri.append(T);
 108         ++N2;
 109       }
 110       if (c2c[id_cell][2] == -1) {
 111         T[0] = pts[0];
 112         T[1] = pts[1];
 113         T[2] = pts[4];
 114         ex_tri.append(T);
 115         T[0] = pts[0];
 116         T[1] = pts[4];
 117         T[2] = pts[3];
 118         ex_tri.append(T);
 119         ++N1;
 120       }
 121       if (c2c[id_cell][3] == -1) {
 122         T[0] = pts[4];
 123         T[1] = pts[1];
 124         T[2] = pts[2];
 125         ex_tri.append(T);
 126         T[0] = pts[4];
 127         T[1] = pts[2];
 128         T[2] = pts[5];
 129         ex_tri.append(T);
 130         ++N1;
 131       }
 132       if (c2c[id_cell][4] == -1) {
 133         T[0] = pts[0];
 134         T[1] = pts[3];
 135         T[2] = pts[2];
 136         ex_tri.append(T);
 137         T[0] = pts[3];
 138         T[1] = pts[5];
 139         T[2] = pts[2];
 140         ex_tri.append(T);
 141         ++N1;
 142       }
 143     } else {
 144       EG_BUG;
 145     }
 146   }
 147   cout << "*********************************************************************" << endl;
 148   cout << "prism quads     : " << N1 << endl;
 149   cout << "prism triangles : " << N2 << endl;
 150   cout << "stray quads     : " << N3 << endl;
 151   cout << "stray triangles : " << N4 << endl;
 152   cout << "*********************************************************************" << endl;
 153   tri.resize(ex_tri.size());
 154   qCopy(ex_tri.begin(), ex_tri.end(), tri.begin());
 155   add_to_ng.fill(false, m_Grid->GetNumberOfPoints());
 156   foreach (QVector<vtkIdType> T, tri) {
 157     add_to_ng[T[0]] = true;
 158     add_to_ng[T[1]] = true;
 159     add_to_ng[T[2]] = true;
 160   }
 161   num_nodes_to_add = 0;
 162   num_old_nodes = 0;
 163   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 164     if (add_to_ng[id_node]) {
 165       ++num_nodes_to_add;
 166     } else {
 167       ++num_old_nodes;
 168     }
 169   }
 170   old2tri.fill(-1, m_Grid->GetNumberOfPoints());
 171   m_NumTriangles = 0;
 172   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 173     if (add_to_ng[id_node]) {
 174       old2tri[id_node] = m_NumTriangles;
 175       ++m_NumTriangles;
 176     }
 177   }
 178   //writeDebugInfo();
 179 }
 180
 181 void CreateVolumeMesh::writeDebugInfo()
 182 {
 183   {
 184     EG_VTKSP(vtkUnstructuredGrid,tri_grid);
 185     allocateGrid(tri_grid, tri.size(), m_NumTriangles);
 186     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 187       if (add_to_ng[id_node]) {
 188         vec3_t x;
 189         m_Grid->GetPoint(id_node, x.data());
 190         tri_grid->GetPoints()->SetPoint(old2tri[id_node], x.data());
 191         copyNodeData(m_Grid, id_node, tri_grid, old2tri[id_node]);
 192       }
 193     }
 194     foreach (QVector<vtkIdType> T, tri) {
 195       vtkIdType pts[3];
 196       pts[0] = old2tri[T[0]];
 197       pts[1] = old2tri[T[1]];
 198       pts[2] = old2tri[T[2]];
 199       tri_grid->InsertNextCell(VTK_TRIANGLE, 3, pts);
 200     }
 201     writeGrid(tri_grid, "triangles");
 202   }
 203   {
 204     writeGrid(m_Grid, "last_grid");
 205   }
 206 }
 207
 208 void CreateVolumeMesh::computeMeshDensity()
 209 {
 210   using namespace nglib;
 211   m_ELSManager.read();
 212   QVector<vtkIdType>  cells;
 213   QVector<vtkIdType>  nodes;
 214   QVector<int>       _nodes;
 215   QVector<QVector<int> >  c2c;
 216   QVector<QSet<int> >     n2n;
 217   getAllCellsOfType(VTK_TETRA, cells, m_Grid);
 218   getNodesFromCells(cells, nodes, m_Grid);
 219   createNodeMapping(nodes, _nodes, m_Grid);
 220   createCellToCell(cells, c2c, m_Grid);
 221   createNodeToNode(cells, nodes, _nodes, n2n, m_Grid);
 222   QVector<bool> fixed(nodes.size(), false);
 223   for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
 224     vtkIdType id_cell = cells[i_cells];
 225     vtkIdType N_pts, *pts;
 226     m_Grid->GetCellPoints(id_cell, N_pts, pts);
 227     if (c2c[i_cells][0] == -1) {
 228       fixed[_nodes[pts[0]]] = true;
 229       fixed[_nodes[pts[1]]] = true;
 230       fixed[_nodes[pts[2]]] = true;
 231     }
 232     if (c2c[i_cells][1] == -1) {
 233       fixed[_nodes[pts[0]]] = true;
 234       fixed[_nodes[pts[1]]] = true;
 235       fixed[_nodes[pts[3]]] = true;
 236     }
 237     if (c2c[i_cells][2] == -1) {
 238       fixed[_nodes[pts[0]]] = true;
 239       fixed[_nodes[pts[2]]] = true;
 240       fixed[_nodes[pts[3]]] = true;
 241     }
 242     if (c2c[i_cells][3] == -1) {
 243       fixed[_nodes[pts[1]]] = true;
 244       fixed[_nodes[pts[2]]] = true;
 245       fixed[_nodes[pts[3]]] = true;
 246     }
 247   }
 248   QVector<double> H(nodes.size(), 0.0);
 249   double H_min = 1e99;
 250   vec3_t X1, X2;
 251   bool first_node = true;
 252   int N_non_fixed = 0;
 253   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 254     if (fixed[i_nodes]) {
 255       ++N_non_fixed;
 256       int N = 0;
 257       vec3_t xi;
 258       m_Grid->GetPoint(nodes[i_nodes], xi.data());
 259       if (first_node) {
 260         X1 = xi;
 261         X2 = xi;
 262         first_node = false;
 263       } else {
 264         for (int k = 0; k < 3; ++k) {
 265           X1[k] = min(xi[k], X1[k]);
 266           X2[k] = max(xi[k], X2[k]);
 267         }
 268       }
 269       foreach (int j_nodes, n2n[i_nodes]) {
 270         if (fixed[j_nodes]) {
 271           vec3_t xj;
 272           m_Grid->GetPoint(nodes[j_nodes], xj.data());
 273           H[i_nodes] += (xi-xj).abs();
 274           ++N;
 275         }
 276       }
 277       if (N < 2) {
 278         EG_BUG;
 279       }
 280       H[i_nodes] /= N;
 281       H_min = min(H[i_nodes], H_min);
 282     }
 283   }
 284   boxes.clear();
 285
 286   // pass 1
 287   {
 288     QString num = "0";
 289     cout << "relaxing mesh size (pass1): " << qPrintable(num) << "% done" << endl;
 290     if (N_non_fixed > 0) {
 291       double DH_max = 1e99;
 292       double DH_last;
 293       do {
 294         DH_last = DH_max;
 295         DH_max = 0;
 296         for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 297           if (!fixed[i_nodes]) {
 298             vec3_t x;
 299             m_Grid->GetPoint(nodes[i_nodes], x.data());
 300             double H0 = H[i_nodes];
 301             H[i_nodes] = 0.0;
 302             int N = 0;
 303             foreach (int j_nodes, n2n[i_nodes]) {
 304               H[i_nodes] += H[j_nodes];
 305               ++N;
 306             }
 307             if (N == 0) {
 308               EG_BUG;
 309             }
 310             H[i_nodes] /= N;
 311             double dH = 1.0*(H[i_nodes] - H0);
 312             H[i_nodes] = H0 + dH;
 313             DH_max = max(dH, DH_max);
 314           }
 315         }
 316         QString new_num;
 317         double e = min(1.0,max(0.0,-log10(DH_max/H_min)/3));
 318         new_num.setNum(100*e,'f',0);
 319         if (new_num != num) {
 320           num = new_num;
 321           cout << "relaxing mesh size (pass1): " << qPrintable(num) << "% done " << endl;
 322         }
 323       } while (DH_max > 1e-3*H_min && DH_max < DH_last);
 324     }
 325   }
 326
 327   // sources
 328   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 329     vec3_t x;
 330     m_Grid->GetPoint(nodes[i_nodes], x.data());
 331     double cl_src = m_ELSManager.minEdgeLength(x);
 332     if (cl_src > 0) {
 333       if (cl_src < H[i_nodes]) {
 334         H[i_nodes] = cl_src;
 335         fixed[i_nodes] = true;
 336         --N_non_fixed;
 337       }
 338     }
 339   }
 340
 341   // pass 2
 342   {
 343     QString num = "0";
 344     cout << "relaxing mesh size (pass2): " << qPrintable(num) << "% done" << endl;
 345     if (N_non_fixed > 0) {
 346       double DH_max = 1e99;
 347       double DH_last;
 348       do {
 349         DH_last = DH_max;
 350         DH_max = 0;
 351         for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 352           if (!fixed[i_nodes]) {
 353             bool relax_node = true;
 354             vec3_t x;
 355             m_Grid->GetPoint(nodes[i_nodes], x.data());
 356             double cl_src = m_ELSManager.minEdgeLength(x);
 357             if (cl_src > 0) {
 358               if (cl_src < H[i_nodes]) {
 359                 H[i_nodes] = cl_src;
 360                 relax_node = false;
 361               }
 362             }
 363             if (relax_node) {
 364               double H0 = H[i_nodes];
 365               H[i_nodes] = 0.0;
 366               int N = 0;
 367               foreach (int j_nodes, n2n[i_nodes]) {
 368                 H[i_nodes] += H[j_nodes];
 369                 ++N;
 370               }
 371               if (N == 0) {
 372                 EG_BUG;
 373               }
 374               H[i_nodes] /= N;
 375               double dH = 1.0*(H[i_nodes] - H0);
 376               H[i_nodes] = H0 + dH;
 377               DH_max = max(dH, DH_max);
 378             }
 379           }
 380         }
 381         QString new_num;
 382         double e = min(1.0,max(0.0,-log10(DH_max/H_min)/3));
 383         new_num.setNum(100*e,'f',0);
 384         if (new_num != num) {
 385           num = new_num;
 386           cout << "relaxing mesh size (pass2): " << qPrintable(num) << "% done " << endl;
 387         }
 388       } while (DH_max > 1e-3*H_min && DH_max < DH_last);
 389     }
 390   }
 391
 392   for (int i_nodes = 0; i_nodes < nodes.size(); ++i_nodes) {
 393     vec3_t x1, x2;
 394     m_Grid->GetPoint(nodes[i_nodes], x1.data());
 395     x2 = x1;
 396     foreach (int j_nodes, n2n[i_nodes]) {
 397       vec3_t xj;
 398       m_Grid->GetPoint(nodes[j_nodes], xj.data());
 399       for (int k = 0; k < 3; ++k) {
 400         x1[k] = min(xj[k], x1[k]);
 401         x2[k] = max(xj[k], x2[k]);
 402       }
 403     }
 404     box_t B;
 405     B.x1 =x1;
 406     B.x2 =x2;
 407     B.h = H[i_nodes];
 408     boxes.append(B);
 409   }
 410 }
 411
 412
 413 void CreateVolumeMesh::operate()
 414 {
 415   using namespace nglib;
 416   if (m_Grid->GetNumberOfCells() == 0) {
 417     EG_ERR_RETURN("The grid appears to be empty.");
 418   }
 419   nglib::Ng_Init();
 420   Ng_Meshing_Parameters mp;
 421   mp.maxh = maxh;
 422   mp.fineness = fineness;
 423   mp.secondorder = 0;
 424   Ng_Mesh *mesh = Ng_NewMesh();
 425   computeMeshDensity();
 426   prepare();
 427   QVector<vtkIdType> ng2eg(num_nodes_to_add+1);
 428   QVector<vtkIdType> eg2ng(m_Grid->GetNumberOfPoints());
 429   {
 430     int N = 1;
 431     for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 432       if (add_to_ng[id_node]) {
 433         vec3_t x;
 434         m_Grid->GetPoints()->GetPoint(id_node, x.data());
 435         Ng_AddPoint(mesh, x.data());
 436         ng2eg[N] = id_node;
 437         eg2ng[id_node] = N;
 438         ++N;
 439       }
 440     }
 441   }
 442
 443   foreach (QVector<vtkIdType> T, tri) {
 444     int trig[3];
 445     for (int i = 0; i < 3; ++i) {
 446       trig[i] = eg2ng[T[i]];
 447     }
 448     Ng_AddSurfaceElement(mesh, NG_TRIG, trig);
 449   }
 450   Ng_Result res;
 451   try {
 452     foreach (box_t B, boxes) Ng_RestrictMeshSizeBox(mesh, B.x1.data(), B.x2.data(), B.h);
 453     res = Ng_GenerateVolumeMesh (mesh, &mp);
 454   } catch (netgen::NgException ng_err) {
 455     writeDebugInfo();
 456     Error err;
 457     QString msg = "Netgen stopped with the following error:\n";
 458     msg += ng_err.What().c_str();
 459     msg += "\n\nDebug information has been saved to:\n" + mainWindow()->getCwd();
 460     err.setType(Error::ExitOperation);
 461     err.setText(msg);
 462     throw err;
 463   }
 464   if (res == NG_OK) {
 465     int Npoints_ng = Ng_GetNP(mesh);
 466     int Ncells_ng  = Ng_GetNE(mesh);
 467     int Nscells_ng = Ng_GetNSE(mesh);
 468     EG_VTKSP(vtkUnstructuredGrid,vol_grid);
 469     allocateGrid(vol_grid, m_Grid->GetNumberOfCells() + Ncells_ng, Npoints_ng + num_old_nodes);
 470     vtkIdType new_point = 0;
 471
 472     // copy existing points
 473     QVector<vtkIdType> old2new(m_Grid->GetNumberOfPoints(), -1);
 474     for (vtkIdType id_point = 0; id_point < m_Grid->GetNumberOfPoints(); ++id_point) {
 475       vec3_t x;
 476       m_Grid->GetPoints()->GetPoint(id_point, x.data());
 477       vol_grid->GetPoints()->SetPoint(new_point, x.data());
 478       copyNodeData(m_Grid, id_point, vol_grid, new_point);
 479       old2new[id_point] = new_point;
 480       ++new_point;
 481     }
 482
 483     // mark all surface nodes coming from NETGEN
 484     QVector<bool> ng_surf_node(Npoints_ng + 1, false);
 485     for (int i = 1; i <= Nscells_ng; ++i) {
 486       int pts[8];
 487       Ng_Surface_Element_Type ng_type;
 488       ng_type = Ng_GetSurfaceElement(mesh, i, pts);
 489       int N = 0;
 490       if (ng_type == NG_TRIG) {
 491         N = 3;
 492       } else if (ng_type == NG_QUAD) {
 493         N = 4;
 494       } else {
 495         EG_BUG;
 496       }
 497       for (int j = 0; j < N; ++j) {
 498         ng_surf_node[pts[j]] = true;
 499       }
 500     }
 501
 502     // add new points from NETGEN
 503     QVector<vtkIdType> ng2new(Npoints_ng+1, -1);
 504     for (int i = 1; i <= Npoints_ng; ++i) {
 505       if (!ng_surf_node[i]) {
 506         vec3_t x;
 507         Ng_GetPoint(mesh, i, x.data());
 508         vol_grid->GetPoints()->SetPoint(new_point, x.data());
 509         ng2new[i] = new_point;
 510         ++new_point;
 511       }
 512     }
 513
 514     // copy existing cells
 515     QVector<vtkIdType> old2new_cell(m_Grid->GetNumberOfCells(), -1);
 516     for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 517       bool ins_cell = false;
 518       vtkIdType type_cell = m_Grid->GetCellType(id_cell);
 519       if (type_cell == VTK_TRIANGLE) ins_cell = true;
 520       if (type_cell == VTK_QUAD)     ins_cell = true;
 521       if (type_cell == VTK_WEDGE)    ins_cell = true;
 522       if (ins_cell) {
 523         vtkIdType N_pts, *pts;
 524         m_Grid->GetCellPoints(id_cell, N_pts, pts);
 525         for (int i = 0; i < N_pts; ++i) {
 526           pts[i] = old2new[pts[i]];
 527           if (pts[i] == -1) {
 528             EG_BUG;
 529           }
 530         }
 531         vtkIdType id_new = vol_grid->InsertNextCell(type_cell, N_pts, pts);
 532         copyCellData(m_Grid, id_cell, vol_grid, id_new);
 533         old2new_cell[id_cell] = id_new;
 534       }
 535     }
 536
 537     // add new cells
 538     vtkIdType id_new_cell;
 539     for (vtkIdType cellId = 0; cellId < Ncells_ng; ++cellId) {
 540       int       pts[8];
 541       vtkIdType new_pts[4];
 542       for (int i = 0; i < 8; ++i) {
 543         pts[i] = 0;
 544       }
 545       Ng_Volume_Element_Type ng_type;
 546       ng_type = Ng_GetVolumeElement(mesh, cellId + 1, pts);
 547       if (ng_type != NG_TET) {
 548         EG_BUG;
 549       }
 550       for (int i = 0; i < 4; ++i) {
 551         if (!ng_surf_node[pts[i]]) {
 552           new_pts[i] = ng2new[pts[i]];
 553         } else {
 554           new_pts[i] = ng2eg[pts[i]];
 555         }
 556       }
 557       if (ng_type == NG_TET) {
 558         vtkIdType tet[4];
 559         tet[0] = new_pts[0];
 560         tet[1] = new_pts[1];
 561         tet[2] = new_pts[3];
 562         tet[3] = new_pts[2];
 563         id_new_cell = vol_grid->InsertNextCell(VTK_TETRA, 4, tet);
 564       } else if (ng_type == NG_PYRAMID) {
 565         EG_ERR_RETURN("pyramids cannot be handled yet");
 566       } else if (ng_type == NG_PRISM) {
 567         EG_ERR_RETURN("prisms cannot be handled yet");
 568       } else {
 569         EG_ERR_RETURN("bug encountered");
 570       }
 571     }
 572     makeCopy(vol_grid, m_Grid);
 573     for (int i = 0; i < trace_cells.size(); ++i) {
 574       if (old2new_cell[trace_cells[i]] == -1) {
 575         EG_BUG;
 576       }
 577       trace_cells[i] = old2new_cell[trace_cells[i]];
 578     }
 579   } else {
 580     Error err;
 581     QString msg = "NETGEN did not succeed.\nPlease check if the surface mesh is oriented correctly";
 582     msg += " (red edges on the outside of the domain)";
 583     err.setType(Error::ExitOperation);
 584     err.setText(msg);
 585     throw err;
 586   }
 587   Ng_DeleteMesh(mesh);
 588   Ng_Exit();
 589   cout << "\n\nNETGEN call finished" << endl;
 590   cout << endl;
 591 }
 592