src/libengrid/fillplane.cpp

   1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   2 // +                                                                      +
   3 // + This file is part of enGrid.                                         +
   4 // +                                                                      +
   5 // + Copyright 2008-2014 enGits GmbH                                      +
   6 // +                                                                      +
   7 // + enGrid is free software: you can redistribute it and/or modify       +
   8 // + it under the terms of the GNU General Public License as published by +
   9 // + the Free Software Foundation, either version 3 of the License, or    +
  10 // + (at your option) any later version.                                  +
  11 // +                                                                      +
  12 // + enGrid is distributed in the hope that it will be useful,            +
  13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  15 // + GNU General Public License for more details.                         +
  16 // +                                                                      +
  17 // + You should have received a copy of the GNU General Public License    +
  18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  19 // +                                                                      +
  20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  21
  22 #include <vtkDelaunay2D.h>
  23 #include <vtkGeometryFilter.h>
  24 #include <vtkCellArray.h>
  25
  26 #include "engrid.h"
  27 #include "fillplane.h"
  28 #include "vtkEgPolyDataToUnstructuredGridFilter.h"
  29 #include "guimainwindow.h"
  30 #include "correctsurfaceorientation.h"
  31
  32 FillPlane::FillPlane()
  33 {
  34   m_InverseDirection = false;
  35   m_DistTol = 0;
  36   m_AngleTol = deg2rad(10);
  37   m_BC = 0;
  38 }
  39
  40 vec3_t FillPlane::toPlane(vec3_t x)
  41 {
  42   x -= m_X0;
  43   x -= (x*m_N)*m_N;
  44   return vec3_t(x*m_G1, x*m_G2, 0);
  45 }
  46
  47 vec3_t FillPlane::fromPlane(vec3_t x)
  48 {
  49   return m_X0 + x[0]*m_G1 + x[1]*m_G2;
  50 }
  51
  52 bool FillPlane::isWithinTolerance(vec3_t x)
  53 {
  54   if (fabs((x - m_X0)*m_N) < m_DistTol) {
  55     return true;
  56   }
  57   return false;
  58 }
  59
  60 void FillPlane::createEdgesOnPlane(vtkUnstructuredGrid *edge_grid)
  61 {
  62   vtkIdType num_edges = 0;
  63   vtkIdType num_nodes = 0;
  64
  65   QVector<bool> is_edge_node(m_Grid->GetNumberOfPoints(), false);
  66   for (vtkIdType id_face = 0; id_face < m_Grid->GetNumberOfCells(); ++id_face) {
  67     if (isSurface(id_face, m_Grid)) {
  68       EG_GET_CELL(id_face, m_Grid);
  69       for (int i = 0; i < num_pts; ++i) {
  70         if (m_Part.c2cGG(id_face, i) == -1) {
  71           vtkIdType id_node1 = pts[i];
  72           vtkIdType id_node2 = pts[0];
  73           if (i < num_pts - 1) {
  74             id_node2 = pts[i + 1];
  75           }
  76           vec3_t x1, x2;
  77           m_Grid->GetPoint(id_node1, x1.data());
  78           m_Grid->GetPoint(id_node2, x2.data());
  79           if (isWithinTolerance(x1) && isWithinTolerance(x2)) {
  80             is_edge_node[id_node1] = true;
  81             is_edge_node[id_node2] = true;
  82             ++num_edges;
  83           }
  84         }
  85       }
  86     }
  87   }
  88
  89
  90   QVector<vtkIdType> node_map(m_Grid->GetNumberOfPoints(), -1);
  91   for (vtkIdType id_node1 = 0; id_node1 < m_Grid->GetNumberOfPoints(); ++id_node1) {
  92     if (is_edge_node[id_node1]) {
  93       node_map[id_node1] = num_nodes;
  94       ++num_nodes;
  95     }
  96   }
  97   m_NodeMap.resize(num_nodes);
  98   allocateGrid(edge_grid, 2*num_edges, num_nodes, false);
  99   for (vtkIdType id_node1 = 0; id_node1 < m_Grid->GetNumberOfPoints(); ++id_node1) {
 100     if (node_map[id_node1] != -1) {
 101       m_NodeMap[node_map[id_node1]] = id_node1;
 102       vec3_t x;
 103       m_Grid->GetPoint(id_node1, x.data());
 104       edge_grid->GetPoints()->SetPoint(node_map[id_node1], x.data());
 105     }
 106   }
 107
 108
 109   for (vtkIdType id_face = 0; id_face < m_Grid->GetNumberOfCells(); ++id_face) {
 110     if (isSurface(id_face, m_Grid)) {
 111       EG_GET_CELL(id_face, m_Grid);
 112       for (int i = 0; i < num_pts; ++i) {
 113         if (m_Part.c2cGG(id_face, i) == -1) {
 114           vtkIdType id_node1 = pts[i];
 115           vtkIdType id_node2 = pts[0];
 116           if (i < num_pts - 1) {
 117             id_node2 = pts[i + 1];
 118           }
 119           if (is_edge_node[id_node1] && is_edge_node[id_node2]) {
 120             vtkIdType pts[2];
 121             pts[0] = node_map[id_node1];
 122             pts[1] = node_map[id_node2];
 123             edge_grid->InsertNextCell(VTK_LINE, 2, pts);
 124           }
 125         }
 126       }
 127     }
 128   }
 129 }
 130
 131 void FillPlane::closeLoops(vtkUnstructuredGrid *edge_grid)
 132 {
 133   bool done = false;
 134   while (!done) {
 135     QList<vtkIdType> end_nodes;
 136     QVector<int> count(edge_grid->GetNumberOfPoints(), 0);
 137     for (vtkIdType id_edge = 0; id_edge < edge_grid->GetNumberOfCells(); ++id_edge) {
 138       EG_GET_CELL(id_edge, edge_grid);
 139       for (int i = 0; i < num_pts; ++i) {
 140         ++count[pts[i]];
 141       }
 142     }
 143     for (vtkIdType id_node = 0; id_node < edge_grid->GetNumberOfPoints(); ++id_node) {
 144       if (count[id_node] == 0) {
 145         EG_ERR_RETURN("unable to fill plane(s)");
 146       }
 147       if (count[id_node] > 2) {
 148         EG_ERR_RETURN("unable to fill plane(s)");
 149       }
 150       if (count[id_node] == 1) {
 151         end_nodes.append(id_node);
 152       }
 153     }
 154     if (end_nodes.size() % 2 != 0) {
 155       EG_ERR_RETURN("unable to fill plane(s)");
 156     }
 157     if (end_nodes.size() > 0) {
 158       double dist_min = EG_LARGE_REAL;
 159       vtkIdType id_fill1 = -1;
 160       vtkIdType id_fill2 = -1;
 161       foreach (vtkIdType id_node1, end_nodes) {
 162         vec3_t n1 = m_Part.globalNormal(m_NodeMap[id_node1]);
 163         vec3_t x1;
 164         edge_grid->GetPoint(id_node1, x1.data());
 165         foreach (vtkIdType id_node2, end_nodes) {
 166           if (id_node1 != id_node2) {
 167             vec3_t n2 = m_Part.globalNormal(m_NodeMap[id_node2]);
 168             vec3_t x2;
 169             edge_grid->GetPoint(id_node2, x2.data());
 170             double angle = GeometryTools::angle(n1, -1*n2);
 171             if (angle < m_AngleTol) {
 172               vec3_t d = x2 - x1;
 173               double dist = d.abs();
 174               d.normalise();
 175               if (d*n2 > 0.5 && d*n1 < -0.5) {
 176                 if (dist < dist_min) {
 177                   dist_min = dist;
 178                   id_fill1 = id_node1;
 179                   id_fill2 = id_node2;
 180                 }
 181               }
 182             }
 183           }
 184         }
 185       }
 186       if (id_fill1 == -1 || id_fill2 == -1) {
 187         EG_ERR_RETURN("unable to fill plane(s)");
 188       }
 189       vtkIdType pts[2];
 190       pts[0] = id_fill1;
 191       pts[1] = id_fill2;
 192       edge_grid->InsertNextCell(VTK_LINE, 2, pts);
 193     } else {
 194       done = true;
 195     }
 196   }
 197 }
 198
 199 void FillPlane::gridToPlane(vtkUnstructuredGrid *edge_grid)
 200 {
 201   for (vtkIdType id_node = 0; id_node < edge_grid->GetNumberOfPoints(); ++id_node) {
 202     vec3_t x;
 203     edge_grid->GetPoint(id_node, x.data());
 204     x = toPlane(x);
 205     edge_grid->GetPoints()->SetPoint(id_node, x.data());
 206   }
 207 }
 208
 209 void FillPlane::gridFromPlane(vtkUnstructuredGrid *edge_grid)
 210 {
 211   for (vtkIdType id_node = 0; id_node < edge_grid->GetNumberOfPoints(); ++id_node) {
 212     vec3_t x;
 213     edge_grid->GetPoint(id_node, x.data());
 214     x = fromPlane(x);
 215     edge_grid->GetPoints()->SetPoint(id_node, x.data());
 216   }
 217 }
 218
 219 void FillPlane::triangulate(vtkPolyData *edge_pdata, vtkUnstructuredGrid *tri_grid, int bc)
 220 {
 221   EG_VTKSP(vtkDelaunay2D, delaunay);
 222   EG_VTKSP(vtkEgPolyDataToUnstructuredGridFilter, pdata2grid);
 223   delaunay->SetInputData(edge_pdata);
 224   delaunay->SetSourceData(edge_pdata);
 225   delaunay->Update();
 226   pdata2grid->SetInputConnection(delaunay->GetOutputPort());
 227   pdata2grid->Update();
 228   makeCopy(pdata2grid->GetOutput(), tri_grid);
 229   createBasicFields(tri_grid, tri_grid->GetNumberOfCells(), tri_grid->GetNumberOfPoints());
 230   if (bc < 0) {
 231     QSet<int> bcs = getAllBoundaryCodes(m_Grid);
 232     m_BC = 0;
 233     foreach (int bc, bcs) {
 234       m_BC = max(m_BC, bc);
 235     }
 236     ++m_BC;
 237     bc = m_BC;
 238   }
 239   EG_VTKDCC(vtkIntArray, cell_code, tri_grid, "cell_code");
 240   for (vtkIdType id_face = 0; id_face < tri_grid->GetNumberOfCells(); ++id_face) {
 241     cell_code->SetValue(id_face, bc);
 242   }
 243 }
 244
 245 void FillPlane::order(vtkUnstructuredGrid *edge_grid, vtkPolyData *edge_pdata)
 246 {
 247   QVector<QVector<vtkIdType> > edges(edge_grid->GetNumberOfCells(), QVector<vtkIdType>(2));
 248   for (vtkIdType id_edge = 0; id_edge < edge_grid->GetNumberOfCells(); ++id_edge) {
 249     EG_GET_CELL(id_edge, edge_grid);
 250     if (num_pts != 2) {
 251       EG_BUG;
 252     }
 253     edges[id_edge][0] = pts[0];
 254     edges[id_edge][1] = pts[1];
 255   }
 256   for (int i = 1; i < edges.size(); ++i) {
 257     for (int j = i; j < edges.size(); ++j) {
 258       QVector<vtkIdType> tmp_edge = edges[j];
 259       edges[j] = edges[i];
 260       if (tmp_edge[0] == edges[i-1][1]) {
 261         edges[i][0] = tmp_edge[0];
 262         edges[i][1] = tmp_edge[1];
 263         break;
 264       }
 265       if (tmp_edge[1] == edges[i-1][1]) {
 266         edges[i][0] = tmp_edge[1];
 267         edges[i][1] = tmp_edge[0];
 268         break;
 269       }
 270       edges[j] = tmp_edge;
 271     }
 272   }
 273   QList<vtkIdType> poly_nodes;
 274   for (vtkIdType id_edge = 0; id_edge < edge_grid->GetNumberOfCells(); ++id_edge) {
 275     poly_nodes.append(edges[id_edge][0]);
 276   }
 277   orderGeometrically(edge_grid, poly_nodes);
 278   EG_VTKSP(vtkPoints, points);
 279   EG_VTKSP(vtkCellArray, polys);
 280   foreach (vtkIdType id_node, poly_nodes) {
 281     vec3_t x;
 282     edge_grid->GetPoint(id_node, x.data());
 283     points->InsertNextPoint(x.data());
 284   }
 285   EG_VTKSP(vtkIdList, pts);
 286   pts->SetNumberOfIds(poly_nodes.size());
 287   for (vtkIdType i = 0; i < pts->GetNumberOfIds(); ++i) {
 288     pts->SetId(i, i);
 289   }
 290   polys->InsertNextCell(pts);
 291   edge_pdata->SetPoints(points);
 292   edge_pdata->SetPolys(polys);
 293 }
 294
 295 void FillPlane::orderGeometrically(vtkUnstructuredGrid* edge_grid, QList<vtkIdType>& poly_nodes)
 296 {
 297   if (poly_nodes.size() < 3) {
 298     return;
 299   }
 300   vec3_t x_centre(0,0,0);
 301   int N = 0;
 302   foreach (vtkIdType id_node, poly_nodes) {
 303     vec3_t x;
 304     edge_grid->GetPoint(id_node, x.data());
 305     x_centre += x;
 306     ++N;
 307   }
 308   x_centre *= 1.0/N;
 309   double scale_sum = 0;
 310   {
 311     bool first = true;
 312     vec3_t x1;
 313     foreach (vtkIdType id_node, poly_nodes) {
 314       vec3_t x;
 315       edge_grid->GetPoint(id_node, x.data());
 316       if (first) {
 317         x1 = x;
 318         first = false;
 319       } else {
 320         vec3_t x2 = x;
 321         vec3_t u = x2 - x1;
 322         vec3_t v = GeometryTools::rotate(u, m_N, deg2rad(90));
 323         vec3_t c = x_centre - 0.5*(x1 + x2);
 324         c.normalise();
 325         scale_sum += c*v;
 326         x1 = x;
 327       }
 328     }
 329   }
 330   bool invert = false;
 331   if (scale_sum < 0) {
 332     if (!m_InverseDirection) {
 333       invert = true;
 334     }
 335   } else {
 336     if (m_InverseDirection) {
 337       invert = true;
 338     }
 339   }
 340   if (invert) {
 341     QList<vtkIdType> reverted;
 342     foreach (vtkIdType id_node, poly_nodes) {
 343       reverted.prepend(id_node);
 344     }
 345     poly_nodes = reverted;
 346   }
 347 }
 348
 349 void FillPlane::setupTransformation()
 350 {
 351   m_G1 = GeometryTools::orthogonalVector(m_N);
 352   m_G2 = m_N.cross(m_G1);
 353   m_N.normalise();
 354   m_G1.normalise();
 355   m_G2.normalise();
 356 }
 357
 358 void FillPlane::operate()
 359 {
 360   setupTransformation();
 361   EG_VTKSP(vtkUnstructuredGrid, edge_grid);
 362   EG_VTKSP(vtkUnstructuredGrid, tri_grid);
 363   EG_VTKSP(vtkPolyData, edge_pdata);
 364   createEdgesOnPlane(edge_grid);
 365   writeGrid(edge_grid, "filltest1");
 366   closeLoops(edge_grid);
 367   writeGrid(edge_grid, "filltest2");
 368   gridToPlane(edge_grid);
 369   order(edge_grid, edge_pdata);
 370   triangulate(edge_pdata, tri_grid);
 371   gridFromPlane(tri_grid);
 372   MeshPartition tri_part(tri_grid, true);
 373   m_Part.addPartition(tri_part, m_DistTol);
 374   m_Grid->Modified();
 375   CorrectSurfaceOrientation corr_surf;
 376   corr_surf();
 377   GuiMainWindow::pointer()->updateBoundaryCodes(true);
 378 }
 379
 380