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
  25 #include <vtkCellArray.h>
  26
  27 void SwapTriangles::prepare()
  28 {
  29   getAllCellsOfType(VTK_TRIANGLE, cells, grid);
  30   QList<vtkIdType> ex_cells;
  31   EG_VTKDCC(vtkIntArray, bc, grid, "cell_code");
  32   foreach (vtkIdType id_cell, cells) {
  33     if (!boundary_codes.contains(bc->GetValue(id_cell))) {
  34       ex_cells.append(id_cell);
  35     };
  36     if (grid->GetCellType(id_cell) != VTK_TRIANGLE) {
  37       EG_BUG;
  38     };
  39   };
  40   cells.resize(ex_cells.size());
  41   qCopy(ex_cells.begin(), ex_cells.end(), cells.begin());
  42   createCellMapping(cells, _cells, grid);
  43   createCellToCell(cells, c2c, grid);
  44   marked.resize(cells.size());
  45 };
  46
  47 void SwapTriangles::operate()
  48 {
  49   cout << "swapping edges of boundary triangles (Delaunay)" << endl;
  50   using namespace GeometryTools;
  51   prepare();
  52   int N_swaps;
  53   int N_total = 0;
  54   int loop = 1;
  55   do {
  56     N_swaps = 0;
  57     createCellToCell(cells, c2c, grid);
  58     //NOTE: This for loop can eventually be removed because if undefined, it's probably false.
  59     for (int i_cells = 0; i_cells < cells.size(); ++i_cells) {
  60       marked[i_cells] = false;
  61     };
  62     foreach (vtkIdType id_cell, cells) {
  63       if (!marked[_cells[id_cell]]) {
  64         for (int j = 0; j < 3; ++j) {
  65           bool swap = false;
  66           stencil_t S = getStencil(id_cell, j);
  67           if (S.valid) {
  68             if (!marked[_cells[S.id_cell2]]) {
  69               vec3_t x3[4], x3_0(0,0,0);
  70               vec2_t x[4];
  71               for (int k = 0; k < 4; ++k) {
  72                 grid->GetPoints()->GetPoint(S.p[k], x3[k].data());
  73                 x3_0 += x3[k];
  74               };
  75               vec3_t n1 = triNormal(x3[0], x3[1], x3[3]);
  76               vec3_t n2 = triNormal(x3[1], x3[2], x3[3]);
  77               n1.normalise();
  78               n2.normalise();
  79               if ( (n1*n2) > 0.8) {
  80                 vec3_t n = n1 + n2;
  81                 n.normalise();
  82                 vec3_t ex = orthogonalVector(n);
  83                 vec3_t ey = ex.cross(n);
  84                 for (int k = 0; k < 4; ++k) {
  85                   x[k] = vec2_t(x3[k]*ex, x3[k]*ey);
  86                 };
  87                 vec2_t r1, r2, r3, u1, u2, u3;
  88                 r1 = 0.5*(x[0] + x[1]); u1 = turnLeft(x[1] - x[0]);
  89                 r2 = 0.5*(x[1] + x[2]); u2 = turnLeft(x[2] - x[1]);
  90                 r3 = 0.5*(x[1] + x[3]); u3 = turnLeft(x[3] - x[1]);
  91                 double k, l;
  92                 vec2_t xm1, xm2;
  93                 bool ok = true;
  94                 if (intersection(k, l, r1, u1, r3, u3)) {
  95                   xm1 = r1 + k*u1;
  96                   if (intersection(k, l, r2, u2, r3, u3)) {
  97                     xm2 = r2 + k*u2;
  98                   } else {
  99                     ok = false;
 100                   };
 101                 } else {
 102                   ok = false;
 103                   swap = true;
 104                 };
 105                 if (ok) {
 106                   if ((xm1 - x[2]).abs() < (xm1 - x[0]).abs()) {
 107                     swap = true;
 108                   };
 109                   if ((xm2 - x[0]).abs() < (xm2 - x[2]).abs()) {
 110                     swap = true;
 111                   };
 112                 };
 113               };
 114             };
 115           };
 116           if (swap) {
 117             marked[_cells[S.id_cell1]] = true;
 118             marked[_cells[S.id_cell2]] = true;
 119             if (c2c[_cells[S.id_cell1]][0] != -1) {
 120               marked[c2c[_cells[S.id_cell1]][0]] = true;
 121             } else if (c2c[_cells[S.id_cell1]][1] != -1) {
 122               marked[c2c[_cells[S.id_cell1]][1]] = true;
 123             } else if (c2c[_cells[S.id_cell1]][2] != -1) {
 124               marked[c2c[_cells[S.id_cell1]][2]] = true;
 125             };
 126             if (c2c[_cells[S.id_cell2]][0] != -1) {
 127               marked[c2c[_cells[S.id_cell2]][0]] = true;
 128             } else if (c2c[_cells[S.id_cell2]][1] != -1) {
 129               marked[c2c[_cells[S.id_cell2]][1]] = true;
 130             } else if (c2c[_cells[S.id_cell2]][2] != -1) {
 131               marked[c2c[_cells[S.id_cell2]][2]] = true;
 132             };
 133             vtkIdType new_pts1[3], new_pts2[3];
 134             new_pts1[0] = S.p[1];
 135             new_pts1[1] = S.p[2];
 136             new_pts1[2] = S.p[0];
 137             new_pts2[0] = S.p[2];
 138             new_pts2[1] = S.p[3];
 139             new_pts2[2] = S.p[0];
 140             grid->ReplaceCell(S.id_cell1, 3, new_pts1);
 141             grid->ReplaceCell(S.id_cell2, 3, new_pts2);
 142             ++N_swaps;
 143             ++N_total;
 144             break;
 145           };
 146         };
 147       };
 148     };
 149     ++loop;
 150   } while ((N_swaps > 0) && (loop <= 20));
 151   cout << N_total << " triangles have been swapped" << endl;
 152 };
 153