src/swaptriangles.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 "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 }
  37
  38 void SwapTriangles::operate()
  39 {
  40   //cout << "swapping edges of boundary triangles (Delaunay)" << endl;
  41
  42   static int nStatic_SwapTriangles;    // Value of nStatic_SwapTriangles is retained between each function call
  43   nStatic_SwapTriangles++;
  44   //cout << "nStatic_SwapTriangles is " << nStatic_SwapTriangles << endl;
  45
  46   int N_swaps;
  47   int N_total = 0;
  48   int loop = 1;
  49   do {
  50     N_swaps = 0;
  51     setAllSurfaceCells();
  52     l2g_t  cells = getPartCells();
  53     g2l_t _cells = getPartLocalCells();
  54     l2l_t  n2c   = getPartN2C();
  55     g2l_t _nodes = getPartLocalNodes();
  56     EG_VTKDCC(vtkIntArray, cell_code, grid, "cell_code");
  57     QVector<bool> l_marked(cells.size());
  58     foreach (vtkIdType id_cell, cells) {
  59       if (!boundary_codes.contains(cell_code->GetValue(id_cell)) && grid->GetCellType(id_cell) == VTK_TRIANGLE) { //if it is a selected triangle
  60         if (!l_marked[_cells[id_cell]]) {
  61           for (int j = 0; j < 3; ++j) {
  62             bool swap = false;
  63             stencil_t S = getStencil(id_cell, j);
  64             if(S.twocells && S.sameBC && (S.neighbour_type == VTK_TRIANGLE)) {
  65               if( ! isEdge(S.p[0],S.p[2]) ) {
  66                 if (!l_marked[_cells[S.id_cell2]]) {
  67                   vec3_t x3[4], x3_0(0,0,0);
  68                   vec2_t x[4];
  69                   for (int k = 0; k < 4; ++k) {
  70                     grid->GetPoints()->GetPoint(S.p[k], x3[k].data());
  71                     x3_0 += x3[k];
  72                   }
  73                   vec3_t n1 = triNormal(x3[0], x3[1], x3[3]);
  74                   vec3_t n2 = triNormal(x3[1], x3[2], x3[3]);
  75                   if ( m_FeatureSwap || (n1*n2) > 0.8*n1.abs()*n2.abs() ) {
  76                     if(testSwap(S)) {
  77                       vec3_t n = n1 + n2;
  78                       n.normalise();
  79                       vec3_t ex = orthogonalVector(n);
  80                       vec3_t ey = ex.cross(n);
  81                       for (int k = 0; k < 4; ++k) {
  82                         x[k] = vec2_t(x3[k]*ex, x3[k]*ey);
  83                       }
  84                       vec2_t r1, r2, r3, u1, u2, u3;
  85                       r1 = 0.5*(x[0] + x[1]); u1 = turnLeft(x[1] - x[0]);
  86                       r2 = 0.5*(x[1] + x[2]); u2 = turnLeft(x[2] - x[1]);
  87                       r3 = 0.5*(x[1] + x[3]); u3 = turnLeft(x[3] - x[1]);
  88                       double k, l;
  89                       vec2_t xm1, xm2;
  90                       bool ok = true;
  91                       if (intersection(k, l, r1, u1, r3, u3)) {
  92                         xm1 = r1 + k*u1;
  93                         if (intersection(k, l, r2, u2, r3, u3)) {
  94                           xm2 = r2 + k*u2;
  95                         } else {
  96                           ok = false;
  97                         }
  98                       } else {
  99                         ok = false;
 100                         swap = true;
 101                       }
 102                       if (ok) {
 103                         if ((xm1 - x[2]).abs() < (xm1 - x[0]).abs()) {
 104                           swap = true;
 105                         }
 106                         if ((xm2 - x[0]).abs() < (xm2 - x[2]).abs()) {
 107                           swap = true;
 108                         }
 109                       }
 110                     }// end of testswap
 111                   } //end of if feature angle
 112                 } //end of if l_marked
 113               } //end of if TestSwap
 114             } //end of S valid
 115
 116             if (swap) {
 117               l_marked[_cells[S.id_cell1]] = true;
 118               l_marked[_cells[S.id_cell2]] = true;
 119               for(int i = 0; i < 4; i++) {
 120                 foreach(int i_neighbour, n2c[_nodes[S.p[i]]]) {
 121                   l_marked[i_neighbour] = true;
 122                 }
 123               }
 124               vtkIdType new_pts1[3], new_pts2[3];
 125               new_pts1[0] = S.p[1];
 126               new_pts1[1] = S.p[2];
 127               new_pts1[2] = S.p[0];
 128               new_pts2[0] = S.p[2];
 129               new_pts2[1] = S.p[3];
 130               new_pts2[2] = S.p[0];
 131               vtkIdType old_pts1[3], old_pts2[3];
 132               old_pts1[0] = S.p[0];
 133               old_pts1[1] = S.p[1];
 134               old_pts1[2] = S.p[3];
 135               old_pts2[0] = S.p[2];
 136               old_pts2[1] = S.p[3];
 137               old_pts2[2] = S.p[1];
 138               grid->ReplaceCell(S.id_cell1, 3, new_pts1);
 139               grid->ReplaceCell(S.id_cell2, 3, new_pts2);
 140               ++N_swaps;
 141               ++N_total;
 142               break;
 143             } //end of if swap
 144           } //end of loop through sides
 145         } //end of if marked
 146       } //end of if selected triangle
 147     } //end of loop through cells
 148     ++loop;
 149   } while ((N_swaps > 0) && (loop <= 20));
 150   //cout << N_total << " triangles have been swapped" << endl;
 151 }
 152
 153 bool SwapTriangles::testSwap(stencil_t S)
 154 {
 155   //old triangles
 156   vec3_t n1_old=triNormal(grid,S.p[0],S.p[1],S.p[3]);
 157   vec3_t n2_old=triNormal(grid,S.p[2],S.p[3],S.p[1]);
 158
 159   //new triangles
 160   vec3_t n1_new=triNormal(grid,S.p[1],S.p[2],S.p[0]);
 161   vec3_t n2_new=triNormal(grid,S.p[3],S.p[0],S.p[2]);
 162
 163   //top point
 164   vec3_t Summit=n1_old+n2_old;
 165   vec3_t M[4];
 166   for (int k = 0; k < 4; ++k) {
 167     grid->GetPoints()->GetPoint(S.p[k], M[k].data());
 168   }
 169
 170   //old volumes
 171   double V1_old=tetraVol(M[0], Summit, M[1], M[3], true);
 172   double V2_old=tetraVol(M[2], Summit, M[3], M[1], true);
 173   //new volumes
 174   double V1_new=tetraVol(M[1], Summit, M[2], M[0], true);
 175   double V2_new=tetraVol(M[3], Summit, M[0], M[2], true);
 176
 177   return(V1_old>0 && V2_old>0 && V1_new>0 && V2_new>0 );
 178 }
 179
 180 bool SwapTriangles::isEdge(vtkIdType id_node1, vtkIdType id_node2)
 181 {
 182   l2g_t  nodes = getPartNodes();
 183   g2l_t _nodes = getPartLocalNodes();
 184   l2l_t  n2n   = getPartN2N();
 185
 186   bool ret = false;
 187   foreach(int i_node, n2n[_nodes[id_node1]]) {
 188     vtkIdType id_node = nodes[i_node];
 189     if( id_node == id_node2 ) ret = true;
 190   }
 191   return(ret);
 192 }