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