src/libengrid/swaptriangles.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2013 enGits GmbH                                      +
   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 "swaptriangles.h"
  24 #include "guimainwindow.h"
  25
  26 #include <vtkCellArray.h>
  27
  28 #include "checksurfaceintegrity.h"
  29 #include "statistics.h"
  30
  31 using namespace GeometryTools;
  32
  33 SwapTriangles::SwapTriangles() : SurfaceOperation()
  34 {
  35   m_RespectBC = false;
  36   m_FeatureSwap = false;
  37   m_FeatureAngle = GeometryTools::deg2rad(30);
  38   m_MaxNumLoops = 20;
  39   m_SmallAreaSwap = false;
  40   m_SmallAreaRatio = 1e-3;
  41   m_Verbose = false;
  42   getSet("surface meshing", "small area ratio for edge-swapping", 1e-3, m_SmallAreaRatio);
  43   getSet("surface meshing", "surface error threshold", 0.05, m_SurfErrorThreshold);
  44   m_SurfErrorThreshold = deg2rad(m_SurfErrorThreshold);
  45 }
  46
  47 bool SwapTriangles::testOrientation(stencil_t S)
  48 {
  49   // old triangles
  50   vec3_t n1_old = triNormal(m_Grid, S.id_node[0], S.p1, S.p2);
  51   vec3_t n2_old = triNormal(m_Grid, S.id_node[1], S.p2, S.p1);
  52
  53   // new triangles
  54   vec3_t n1_new = triNormal(m_Grid, S.p1, S.id_node[1], S.id_node[0]);
  55   vec3_t n2_new = triNormal(m_Grid, S.p2, S.id_node[0], S.id_node[1]);
  56
  57   // top point
  58   vec3_t x_summit(0,0,0);
  59   vec3_t x[4];
  60   double l_max = 0;
  61   m_Grid->GetPoints()->GetPoint(S.id_node[0], x[0].data());
  62   m_Grid->GetPoints()->GetPoint(S.p1,         x[1].data());
  63   m_Grid->GetPoints()->GetPoint(S.id_node[1], x[2].data());
  64   m_Grid->GetPoints()->GetPoint(S.p2,         x[3].data());
  65   for (int k = 0; k < 4; ++k) {
  66     x_summit += x[k];
  67   }
  68   for (int k = 0; k < 4; ++k) {
  69     int i = k;
  70     int j = k + 1;
  71     if (j == 4) {
  72       j = 0;
  73     }
  74     l_max = max(l_max, (x[i]-x[j]).abs());
  75   }
  76   x_summit *= 0.25;
  77   vec3_t n = n1_old + n2_old;
  78   n.normalise();
  79   x_summit += 3*l_max*n;
  80
  81   // old volumes
  82   double V1_old = tetraVol(x[0], x[1], x[3], x_summit, true);
  83   double V2_old = tetraVol(x[2], x[3], x[1], x_summit, true);
  84   // new volumes
  85   double V1_new = tetraVol(x[1], x[2], x[0], x_summit, true);
  86   double V2_new = tetraVol(x[3], x[0], x[2], x_summit, true);
  87
  88   bool swap_ok = false;
  89   if (m_SmallAreaSwap) {
  90      swap_ok = V1_new>0 && V2_new>0;
  91   } else {
  92      swap_ok = V1_old>0 && V2_old>0 && V1_new>0 && V2_new>0;
  93   }
  94   return swap_ok;
  95 }
  96
  97 bool SwapTriangles::isEdge(vtkIdType id_node1, vtkIdType id_node2)
  98 {
  99   l2g_t  nodes = getPartNodes();
 100   g2l_t _nodes = getPartLocalNodes();
 101   l2l_t  n2n   = getPartN2N();
 102
 103   bool ret = false;
 104   foreach(int i_node, n2n[_nodes[id_node1]]) {
 105     vtkIdType id_node = nodes[i_node];
 106     if( id_node == id_node2 ) ret = true;
 107   }
 108   return(ret);
 109 }
 110
 111 double SwapTriangles::computeSurfaceDistance(vec3_t x1, vec3_t x2, vec3_t x3, CadInterface *cad_interface)
 112 {
 113   static const double f13 = 1.0/3.0;
 114   vec3_t x = f13*(x1 + x2 + x3);
 115   vec3_t n = triNormal(x1, x2, x3);
 116   n.normalise();
 117   if (!checkVector(n)) {
 118     return 0;
 119   }
 120   vec3_t xp = cad_interface->snap(x, false);
 121   return (x - xp).abs();
 122 }
 123
 124 double SwapTriangles::edgeAngle(vec3_t x1, vec3_t x2, vec3_t x3, vec3_t x4)
 125 {
 126   static const double f13 = 1.0/3.0;
 127   vec3_t xf1 = f13*(x1 + x2 + x4);
 128   vec3_t xf2 = f13*(x2 + x3 + x4);
 129   vec3_t xe  = 0.5*(x2 + x4);
 130   vec3_t xfe = 0.5*(xf1 + xf2);
 131   vec3_t n1  = triNormal(x1, x2, x4);
 132   vec3_t n2  = triNormal(x2, x3, x4);
 133   vec3_t ve  = xe - xfe;
 134   double alpha = angle(n1, n2);
 135   if ((n1 + n2)*ve < 0) {
 136     return -alpha;
 137   }
 138   return alpha;
 139 }
 140
 141 vtkIdType SwapTriangles::neighbourNode(vtkIdType id_node0, vtkIdType id_node1, vtkIdType id_node2)
 142 {
 143   for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
 144     vtkIdType id_neigh = m_Part.n2nGG(id_node1, i);
 145     if (id_neigh != id_node0) {
 146       if (m_Part.hasNeighNode(id_node2, id_neigh)) {
 147         return id_neigh;
 148       }
 149     }
 150   }
 151   EG_BUG;
 152   return -1;
 153 }
 154
 155 bool SwapTriangles::swapDueToSurfaceNoise(stencil_t S)
 156 {
 157   vtkIdType p1  = S.id_node[0];
 158   vtkIdType p2  = S.p1;
 159   vtkIdType p3  = S.id_node[1];
 160   vtkIdType p4  = S.p2;
 161   vtkIdType p12 = neighbourNode(p4, p1, p2);
 162   vtkIdType p23 = neighbourNode(p4, p2, p3);
 163   vtkIdType p34 = neighbourNode(p2, p3, p4);
 164   vtkIdType p41 = neighbourNode(p2, p4, p1);
 165   vec3_t x1, x2, x3, x4, x12, x23, x34, x41;
 166   m_Grid->GetPoint(p1, x1.data());
 167   m_Grid->GetPoint(p2, x2.data());
 168   m_Grid->GetPoint(p3, x3.data());
 169   m_Grid->GetPoint(p4, x4.data());
 170   m_Grid->GetPoint(p12, x12.data());
 171   m_Grid->GetPoint(p23, x23.data());
 172   m_Grid->GetPoint(p34, x34.data());
 173   m_Grid->GetPoint(p41, x41.data());
 174
 175   double noise1, noise2;
 176   {
 177     double alpha_min =  M_PI;
 178     double alpha_max = -M_PI;
 179     double alpha;
 180
 181     alpha = edgeAngle(x4, x1, x12, x2);
 182     alpha_max = max(alpha, alpha_max);
 183     alpha_min = min(alpha, alpha_min);
 184
 185     alpha = edgeAngle(x4, x2, x23, x3);
 186     alpha_max = max(alpha, alpha_max);
 187     alpha_min = min(alpha, alpha_min);
 188
 189     alpha = edgeAngle(x2, x3, x34, x4);
 190     alpha_max = max(alpha, alpha_max);
 191     alpha_min = min(alpha, alpha_min);
 192
 193     alpha = edgeAngle(x2, x4, x41, x1);
 194     alpha_max = max(alpha, alpha_max);
 195     alpha_min = min(alpha, alpha_min);
 196
 197     alpha = edgeAngle(x1, x2, x3, x4);
 198     alpha_max = max(alpha, alpha_max);
 199     alpha_min = min(alpha, alpha_min);
 200
 201     noise1 = alpha_max - alpha_min;
 202   }
 203   {
 204     double alpha_min =  M_PI;
 205     double alpha_max = -M_PI;
 206     double alpha;
 207
 208     alpha = edgeAngle(x3, x1, x12, x2);
 209     alpha_max = max(alpha, alpha_max);
 210     alpha_min = min(alpha, alpha_min);
 211
 212     alpha = edgeAngle(x1, x2, x23, x3);
 213     alpha_max = max(alpha, alpha_max);
 214     alpha_min = min(alpha, alpha_min);
 215
 216     alpha = edgeAngle(x1, x3, x34, x4);
 217     alpha_max = max(alpha, alpha_max);
 218     alpha_min = min(alpha, alpha_min);
 219
 220     alpha = edgeAngle(x3, x4, x41, x1);
 221     alpha_max = max(alpha, alpha_max);
 222     alpha_min = min(alpha, alpha_min);
 223
 224     alpha = edgeAngle(x2, x3, x4, x1);
 225     alpha_max = max(alpha, alpha_max);
 226     alpha_min = min(alpha, alpha_min);
 227
 228     noise2 = alpha_max - alpha_min;
 229   }
 230
 231   if (noise1 > m_SurfErrorThreshold && noise1 > noise2) {
 232     return true;
 233   }
 234
 235   return false;
 236 }
 237
 238 void SwapTriangles::computeSurfaceErrors(const QVector<vec3_t> &x, int bc, double &err1, double &err2)
 239 {
 240   using namespace GeometryTools;
 241   err1 = 0;
 242   err2 = 1e10;
 243   CadInterface* cad_interface = GuiMainWindow::pointer()->getCadInterface(bc);
 244   if (!cad_interface) {
 245     return;
 246   }
 247
 248   double d013 = computeSurfaceDistance(x[0], x[1], x[3], cad_interface);
 249   double d123 = computeSurfaceDistance(x[1], x[2], x[3], cad_interface);
 250   double d012 = computeSurfaceDistance(x[0], x[1], x[2], cad_interface);
 251   double d023 = computeSurfaceDistance(x[0], x[2], x[3], cad_interface);
 252
 253   double L = max((x[1] - x[3]).abs(), (x[0] - x[2]).abs());
 254   err1 = max(d013, d123)/L;
 255   err2 = max(d012, d023)/L;
 256   if (err1 < m_SurfErrorThreshold) {
 257     err1 = 0;
 258   }
 259   if (err2 < m_SurfErrorThreshold) {
 260     err2 = 0;
 261   }
 262
 263 }
 264
 265 int SwapTriangles::swap()
 266 {
 267   int N_swaps = 0;
 268   setAllSurfaceCells();
 269   //computeAverageSurfaceError();
 270   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 271   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 272   QVector<bool> marked(m_Grid->GetNumberOfCells(), false);
 273   for (int i = 0; i < m_Part.getNumberOfCells(); ++i) {
 274     vtkIdType id_cell = m_Part.globalCell(i);
 275     if (!m_BoundaryCodes.contains(cell_code->GetValue(id_cell)) && m_Grid->GetCellType(id_cell) == VTK_TRIANGLE) { //if it is a selected triangle
 276       if (!marked[id_cell] && !m_Swapped[id_cell]) {
 277         for (int j = 0; j < 3; ++j) {
 278           bool swap = false;
 279           bool feature_improvement_swap = false;
 280           stencil_t S = getStencil(id_cell, j);
 281           if(S.id_cell.size() == 2 && S.sameBC) {
 282             if (S.type_cell[1] == VTK_TRIANGLE) {
 283               if(!isEdge(S.id_node[0], S.id_node[1]) ) {
 284                 if (!marked[S.id_cell[1]] && !m_Swapped[S.id_cell[1]]) {
 285                   QVector<vec3_t> x3(4);
 286                   vec2_t x[4];
 287
 288                   m_Grid->GetPoint(S.id_node[0], x3[0].data());
 289                   m_Grid->GetPoint(S.p1,         x3[1].data());
 290                   m_Grid->GetPoint(S.id_node[1], x3[2].data());
 291                   m_Grid->GetPoint(S.p2,         x3[3].data());
 292
 293                   vec3_t n1 = triNormal(x3[0], x3[1], x3[3]);
 294                   vec3_t n2 = triNormal(x3[1], x3[2], x3[3]);
 295
 296                   bool force_swap = false;
 297                   if (m_SmallAreaSwap) {
 298                     double A1 = n1.abs();
 299                     double A2 = n2.abs();
 300                     if (isnan(A1) || isnan(A2)) {
 301                       force_swap = true;
 302                     } else {
 303                       force_swap = A1 < m_SmallAreaRatio*A2 || A2 < m_SmallAreaRatio*A1;
 304                     }
 305                   }
 306                   if (!isFeatureNode(S.p1) && !isFeatureNode(S.p2)) {
 307
 308                     if (isFeatureNode(S.id_node[0]) && isFeatureNode(S.id_node[1])) {
 309                       //swap = true;
 310                     }
 311
 312                     // Delaunay
 313                     if (!swap) {
 314                       vec3_t n = n1 + n2;
 315                       n.normalise();
 316                       vec3_t ex = orthogonalVector(n);
 317                       vec3_t ey = ex.cross(n);
 318                       for (int k = 0; k < 4; ++k) {
 319                         x[k] = vec2_t(x3[k]*ex, x3[k]*ey);
 320                       }
 321                       vec2_t r1, r2, r3, u1, u2, u3;
 322                       r1 = 0.5*(x[0] + x[1]); u1 = turnLeft(x[1] - x[0]);
 323                       r2 = 0.5*(x[1] + x[2]); u2 = turnLeft(x[2] - x[1]);
 324                       r3 = 0.5*(x[1] + x[3]); u3 = turnLeft(x[3] - x[1]);
 325                       double k, l;
 326                       vec2_t xm1, xm2;
 327                       bool ok = true;
 328                       if (intersection(k, l, r1, u1, r3, u3)) {
 329                         xm1 = r1 + k*u1;
 330                         if (intersection(k, l, r2, u2, r3, u3)) {
 331                           xm2 = r2 + k*u2;
 332                         } else {
 333                           ok = false;
 334                         }
 335                       } else {
 336                         ok = false;
 337                         swap = true;
 338                       }
 339                       if (ok) {
 340                         if ((xm1 - x[2]).abs() < (xm1 - x[0]).abs()) {
 341                           swap = true;
 342                         }
 343                         if ((xm2 - x[0]).abs() < (xm2 - x[2]).abs()) {
 344                           swap = true;
 345                         }
 346                       }
 347                     }
 348                   } //end of if feature angle
 349                 } //end of if l_marked
 350               } //end of if TestSwap
 351             }
 352           } //end of S valid
 353
 354           if (swap) {
 355             if (testOrientation(S)) {
 356               marked[S.id_cell[0]] = true;
 357               marked[S.id_cell[1]] = true;
 358               for (int k = 0; k < m_Part.n2cGSize(S.id_node[0]); ++k) {
 359                 vtkIdType id_neigh = m_Part.n2cGG(S.id_node[0], k);
 360                 marked[id_neigh] = true;
 361               }
 362               for (int k = 0; k < m_Part.n2cGSize(S.id_node[1]); ++k) {
 363                 vtkIdType id_neigh = m_Part.n2cGG(S.id_node[1], k);
 364                 marked[id_neigh] = true;
 365               }
 366               for (int k = 0; k < m_Part.n2cGSize(S.p1); ++k) {
 367                 vtkIdType id_neigh = m_Part.n2cGG(S.p1, k);
 368                 marked[id_neigh] = true;
 369               }
 370               for (int k = 0; k < m_Part.n2cGSize(S.p2); ++k) {
 371                 vtkIdType id_neigh = m_Part.n2cGG(S.p2, k);
 372                 marked[id_neigh] = true;
 373               }
 374               vtkIdType new_pts1[3], new_pts2[3];
 375               new_pts1[0] = S.p1;
 376               new_pts1[1] = S.id_node[1];
 377               new_pts1[2] = S.id_node[0];
 378               new_pts2[0] = S.id_node[1];
 379               new_pts2[1] = S.p2;
 380               new_pts2[2] = S.id_node[0];
 381               vtkIdType old_pts1[3], old_pts2[3];
 382               old_pts1[0] = S.id_node[0];
 383               old_pts1[1] = S.p1;
 384               old_pts1[2] = S.p2;
 385               old_pts2[0] = S.id_node[1];
 386               old_pts2[1] = S.p2;
 387               old_pts2[2] = S.p1;
 388               m_Grid->ReplaceCell(S.id_cell[0], 3, new_pts1);
 389               m_Grid->ReplaceCell(S.id_cell[1], 3, new_pts2);
 390               m_Swapped[S.id_cell[0]] = true;
 391               m_Swapped[S.id_cell[1]] = true;
 392               ++N_swaps;
 393               break;
 394             }
 395           } //end of if swap
 396         } //end of loop through sides
 397       } //end of if marked
 398     } //end of if selected triangle
 399   } //end of loop through cells
 400   return N_swaps;
 401 }
 402
 403 void SwapTriangles::computeAverageSurfaceError()
 404 {
 405   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 406   QList<double> errors;
 407   EG_FORALL_NODES(id_node1, m_Grid) {
 408     vec3_t n1 = m_Part.globalNormal(id_node1);
 409     vec3_t x1;
 410     m_Grid->GetPoint(id_node1, x1.data());
 411     for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
 412       vtkIdType id_node2 = m_Part.n2nGG(id_node1, i);
 413       if (id_node1 < id_node2) {
 414         QVector<vtkIdType> faces;
 415         getEdgeCells(id_node1, id_node2, faces);
 416         if (faces.size() == 2) {
 417           int bc1 = cell_code->GetValue(faces[0]);
 418           int bc2 = cell_code->GetValue(faces[1]);
 419           if (bc1 == bc2) {
 420             CadInterface* cad_interface = GuiMainWindow::pointer()->getCadInterface(bc1);
 421             if (cad_interface) {
 422               vec3_t n2 = m_Part.globalNormal(id_node2);
 423               vec3_t n = 0.5*(n1 + n2);
 424               vec3_t x2;
 425               m_Grid->GetPoint(id_node2, x2.data());
 426               vec3_t x = 0.5*(x1 + x2);
 427               vec3_t xp = cad_interface->snap(x);
 428               double err = (x - xp).abs()/(x1 - x2).abs();
 429               errors.append(err);
 430             }
 431           }
 432         }
 433       }
 434     }
 435   }
 436
 437   m_AverageSurfaceError = Statistics::meanValue(errors);
 438   m_SurfaceErrorDeviation = Statistics::standardDeviation(errors, m_AverageSurfaceError);
 439   cout << "mean: " << m_AverageSurfaceError << "  sigma: " << m_SurfaceErrorDeviation << endl;
 440 }
 441
 442 void SwapTriangles::operate()
 443 {
 444   if (m_Verbose) cout << "swapping edges for surface triangles ..." << endl;
 445   long int N_swaps      = 100000000;
 446   long int N_last_swaps = 100000001;
 447   int loop = 1;
 448   while ((N_swaps > 0) && (loop <= m_MaxNumLoops) && (N_swaps < N_last_swaps)) {
 449     N_last_swaps = N_swaps;
 450     if (m_Verbose) cout << "  loop " << loop << "/" << m_MaxNumLoops << endl;
 451     m_Swapped.fill(false, m_Grid->GetNumberOfCells());
 452     N_swaps = 0;
 453     int N;
 454     int sub_loop = 0;
 455     do {
 456       ++sub_loop;
 457       N = swap();
 458       if (m_Verbose) cout << "    sub-loop " << sub_loop << ": " << N << " swaps" << endl;
 459       N_swaps += N;
 460     } while (N > 0);
 461     ++loop;
 462   }
 463 }
 464