src/libengrid/swaptriangles.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2012 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
  30 using namespace GeometryTools;
  31
  32 SwapTriangles::SwapTriangles() : SurfaceOperation()
  33 {
  34   m_RespectBC = false;
  35   m_FeatureSwap = false;
  36   m_FeatureAngle = GeometryTools::deg2rad(30);
  37   m_MaxNumLoops = 20;
  38   m_SmallAreaSwap = false;
  39   m_SmallAreaRatio = 1e-3;
  40   m_Verbose = false;
  41   getSet("surface meshing", "small area ratio for edge-swapping", 1e-3, m_SmallAreaRatio);
  42   getSet("surface meshing", "threshold for surface errors", 10.0, m_SurfErrorThreshold);
  43 }
  44
  45 bool SwapTriangles::testOrientation(stencil_t S)
  46 {
  47   // old triangles
  48   vec3_t n1_old = triNormal(m_Grid, S.id_node[0], S.p1, S.p2);
  49   vec3_t n2_old = triNormal(m_Grid, S.id_node[1], S.p2, S.p1);
  50
  51   // new triangles
  52   vec3_t n1_new = triNormal(m_Grid, S.p1, S.id_node[1], S.id_node[0]);
  53   vec3_t n2_new = triNormal(m_Grid, S.p2, S.id_node[0], S.id_node[1]);
  54
  55   // top point
  56   vec3_t x_summit(0,0,0);
  57   vec3_t x[4];
  58   double l_max = 0;
  59   m_Grid->GetPoints()->GetPoint(S.id_node[0], x[0].data());
  60   m_Grid->GetPoints()->GetPoint(S.p1,         x[1].data());
  61   m_Grid->GetPoints()->GetPoint(S.id_node[1], x[2].data());
  62   m_Grid->GetPoints()->GetPoint(S.p2,         x[3].data());
  63   for (int k = 0; k < 4; ++k) {
  64     x_summit += x[k];
  65   }
  66   for (int k = 0; k < 4; ++k) {
  67     int i = k;
  68     int j = k + 1;
  69     if (j == 4) {
  70       j = 0;
  71     }
  72     l_max = max(l_max, (x[i]-x[j]).abs());
  73   }
  74   x_summit *= 0.25;
  75   vec3_t n = n1_old + n2_old;
  76   n.normalise();
  77   x_summit += 3*l_max*n;
  78
  79   // old volumes
  80   double V1_old = tetraVol(x[0], x[1], x[3], x_summit, true);
  81   double V2_old = tetraVol(x[2], x[3], x[1], x_summit, true);
  82   // new volumes
  83   double V1_new = tetraVol(x[1], x[2], x[0], x_summit, true);
  84   double V2_new = tetraVol(x[3], x[0], x[2], x_summit, true);
  85
  86   bool swap_ok = false;
  87   if (m_SmallAreaSwap) {
  88      swap_ok = V1_new>0 && V2_new>0;
  89   } else {
  90      swap_ok = V1_old>0 && V2_old>0 && V1_new>0 && V2_new>0;
  91   }
  92   return swap_ok;
  93 }
  94
  95 bool SwapTriangles::isEdge(vtkIdType id_node1, vtkIdType id_node2)
  96 {
  97   l2g_t  nodes = getPartNodes();
  98   g2l_t _nodes = getPartLocalNodes();
  99   l2l_t  n2n   = getPartN2N();
 100
 101   bool ret = false;
 102   foreach(int i_node, n2n[_nodes[id_node1]]) {
 103     vtkIdType id_node = nodes[i_node];
 104     if( id_node == id_node2 ) ret = true;
 105   }
 106   return(ret);
 107 }
 108
 109 void SwapTriangles::computeSurfaceErrors(const QVector<vec3_t> &x, int bc, double &err1, double &err2)
 110 {
 111   static int counter = 0;
 112   using namespace GeometryTools;
 113   err1 = 0;
 114   err2 = 0;
 115   if (m_SurfErrorThreshold <= 0) {
 116     return;
 117   }
 118   SurfaceProjection* proj = GuiMainWindow::pointer()->getSurfProj(bc);
 119   if (!proj) {
 120     return;
 121   }
 122
 123   vec3_t n = triNormal(x[0], x[1], x[3]) + triNormal(x[1], x[2], x[3]);
 124   n.normalise();
 125
 126   vec3_t xe1 = 0.5*(x[1] + x[3]);
 127   vec3_t xe2 = 0.5*(x[0] + x[2]);
 128
 129   vec3_t xe1_proj = proj->project(xe1, -1, true, n);
 130   vec3_t xe2_proj = proj->project(xe2, -1, true, n);
 131
 132   double L = 0;
 133   L = max(L, (x[0] - x[1]).abs());
 134   L = max(L, (x[0] - x[2]).abs());
 135   L = max(L, (x[0] - x[3]).abs());
 136   L = max(L, (x[1] - x[2]).abs());
 137   L = max(L, (x[1] - x[3]).abs());
 138   L = max(L, (x[2] - x[3]).abs());
 139
 140   err1 = (xe1 - xe1_proj).abs()/L;
 141   err2 = (xe2 - xe2_proj).abs()/L;
 142   ++counter;
 143 }
 144
 145 int SwapTriangles::swap()
 146 {
 147   int N_swaps = 0;
 148   setAllSurfaceCells();
 149   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 150   QVector<bool> marked(m_Grid->GetNumberOfCells(), false);
 151   for (int i = 0; i < m_Part.getNumberOfCells(); ++i) {
 152     vtkIdType id_cell = m_Part.globalCell(i);
 153     if (!m_BoundaryCodes.contains(cell_code->GetValue(id_cell)) && m_Grid->GetCellType(id_cell) == VTK_TRIANGLE) { //if it is a selected triangle
 154       if (!marked[id_cell] && !m_Swapped[id_cell]) {
 155         for (int j = 0; j < 3; ++j) {
 156           bool swap = false;
 157           stencil_t S = getStencil(id_cell, j);
 158           if (S.id_cell.size() == 2) {
 159             if ((S.id_cell[1] == 2171 && S.id_cell[0] == 2110) || (S.id_cell[0] == 2171 && S.id_cell[1] == 2110)) {
 160               cout << "break" << endl;
 161             }
 162           }
 163           if(S.id_cell.size() == 2 && S.sameBC) {
 164             if (S.type_cell[1] == VTK_TRIANGLE) {
 165               if(!isEdge(S.id_node[0], S.id_node[1]) ) {
 166                 if (!marked[S.id_cell[1]] && !m_Swapped[S.id_cell[1]]) {
 167                   QVector<vec3_t> x3(4);
 168                   vec2_t x[4];
 169
 170                   m_Grid->GetPoint(S.id_node[0], x3[0].data());
 171                   m_Grid->GetPoint(S.p1,         x3[1].data());
 172                   m_Grid->GetPoint(S.id_node[1], x3[2].data());
 173                   m_Grid->GetPoint(S.p2,         x3[3].data());
 174
 175                   vec3_t n1 = triNormal(x3[0], x3[1], x3[3]);
 176                   vec3_t n2 = triNormal(x3[1], x3[2], x3[3]);
 177
 178                   bool force_swap = false;
 179                   if (m_SmallAreaSwap) {
 180                     double A1 = n1.abs();
 181                     double A2 = n2.abs();
 182                     if (isnan(A1) || isnan(A2)) {
 183                       force_swap = true;
 184                     } else {
 185                       force_swap = A1 < m_SmallAreaRatio*A2 || A2 < m_SmallAreaRatio*A1;
 186                     }
 187                   }
 188                   if (m_FeatureSwap || GeometryTools::angle(n1, n2) < m_FeatureAngle || force_swap) {
 189
 190                     // surface errors
 191                     double se1, se2;
 192                     bool surf_block = false;
 193                     if (m_SurfErrorThreshold > 0) {
 194                       computeSurfaceErrors(x3, cell_code->GetValue(id_cell), se1, se2);
 195                       if (se2 > se1 && se2 > m_SurfErrorThreshold) {
 196                         surf_block = true;
 197                       } else {
 198                         if (se2 < se1 && se1 > m_SurfErrorThreshold) {
 199                           swap = true;
 200                         }
 201                       }
 202                     }
 203
 204                     // Delaunay
 205                     if (!swap && !surf_block) {
 206                       vec3_t n = n1 + n2;
 207                       n.normalise();
 208                       vec3_t ex = orthogonalVector(n);
 209                       vec3_t ey = ex.cross(n);
 210                       for (int k = 0; k < 4; ++k) {
 211                         x[k] = vec2_t(x3[k]*ex, x3[k]*ey);
 212                       }
 213                       vec2_t r1, r2, r3, u1, u2, u3;
 214                       r1 = 0.5*(x[0] + x[1]); u1 = turnLeft(x[1] - x[0]);
 215                       r2 = 0.5*(x[1] + x[2]); u2 = turnLeft(x[2] - x[1]);
 216                       r3 = 0.5*(x[1] + x[3]); u3 = turnLeft(x[3] - x[1]);
 217                       double k, l;
 218                       vec2_t xm1, xm2;
 219                       bool ok = true;
 220                       if (intersection(k, l, r1, u1, r3, u3)) {
 221                         xm1 = r1 + k*u1;
 222                         if (intersection(k, l, r2, u2, r3, u3)) {
 223                           xm2 = r2 + k*u2;
 224                         } else {
 225                           ok = false;
 226                         }
 227                       } else {
 228                         ok = false;
 229                         swap = true;
 230                       }
 231                       if (ok) {
 232                         if ((xm1 - x[2]).abs() < (xm1 - x[0]).abs()) {
 233                           swap = true;
 234                         }
 235                         if ((xm2 - x[0]).abs() < (xm2 - x[2]).abs()) {
 236                           swap = true;
 237                         }
 238                       }
 239                     }
 240                   //}// end of testswap
 241                   } //end of if feature angle
 242                 } //end of if l_marked
 243               } //end of if TestSwap
 244             }
 245           } //end of S valid
 246
 247           if (swap) {
 248             if (testOrientation(S)) {
 249               marked[S.id_cell[0]] = true;
 250               marked[S.id_cell[1]] = true;
 251               for (int k = 0; k < m_Part.n2cGSize(S.id_node[0]); ++k) {
 252                 vtkIdType id_neigh = m_Part.n2cGG(S.id_node[0], k);
 253                 marked[id_neigh] = true;
 254               }
 255               for (int k = 0; k < m_Part.n2cGSize(S.id_node[1]); ++k) {
 256                 vtkIdType id_neigh = m_Part.n2cGG(S.id_node[1], k);
 257                 marked[id_neigh] = true;
 258               }
 259               for (int k = 0; k < m_Part.n2cGSize(S.p1); ++k) {
 260                 vtkIdType id_neigh = m_Part.n2cGG(S.p1, k);
 261                 marked[id_neigh] = true;
 262               }
 263               for (int k = 0; k < m_Part.n2cGSize(S.p2); ++k) {
 264                 vtkIdType id_neigh = m_Part.n2cGG(S.p2, k);
 265                 marked[id_neigh] = true;
 266               }
 267               vtkIdType new_pts1[3], new_pts2[3];
 268               new_pts1[0] = S.p1;
 269               new_pts1[1] = S.id_node[1];
 270               new_pts1[2] = S.id_node[0];
 271               new_pts2[0] = S.id_node[1];
 272               new_pts2[1] = S.p2;
 273               new_pts2[2] = S.id_node[0];
 274               vtkIdType old_pts1[3], old_pts2[3];
 275               old_pts1[0] = S.id_node[0];
 276               old_pts1[1] = S.p1;
 277               old_pts1[2] = S.p2;
 278               old_pts2[0] = S.id_node[1];
 279               old_pts2[1] = S.p2;
 280               old_pts2[2] = S.p1;
 281               m_Grid->ReplaceCell(S.id_cell[0], 3, new_pts1);
 282               m_Grid->ReplaceCell(S.id_cell[1], 3, new_pts2);
 283               m_Swapped[S.id_cell[0]] = true;
 284               m_Swapped[S.id_cell[1]] = true;
 285               ++N_swaps;
 286               break;
 287             }
 288           } //end of if swap
 289         } //end of loop through sides
 290       } //end of if marked
 291     } //end of if selected triangle
 292   } //end of loop through cells
 293   return N_swaps;
 294 }
 295
 296 void SwapTriangles::operate()
 297 {
 298   if (m_Verbose) cout << "swapping edges for surface triangles ..." << endl;
 299   long int N_swaps      = 100000000;
 300   long int N_last_swaps = 100000001;
 301   int loop = 1;
 302   while ((N_swaps > 0) && (loop <= m_MaxNumLoops) && (N_swaps < N_last_swaps)) {
 303     N_last_swaps = N_swaps;
 304     if (m_Verbose) cout << "  loop " << loop << "/" << m_MaxNumLoops << endl;
 305     m_Swapped.fill(false, m_Grid->GetNumberOfCells());
 306     N_swaps = 0;
 307     int N;
 308     int sub_loop = 0;
 309     do {
 310       ++sub_loop;
 311       N = swap();
 312       if (m_Verbose) cout << "    sub-loop " << sub_loop << ": " << N << " swaps" << endl;
 313       N_swaps += N;
 314     } while (N > 0);
 315     ++loop;
 316   }
 317 }
 318