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 "fillplane.h"
  27 #include "vtkEgPolyDataToUnstructuredGridFilter.h"
  28 #include "guimainwindow.h"
  29 #include "correctsurfaceorientation.h"
  30
  31 FillPlane::FillPlane()
  32 {
  33   m_InverseDirection = false;
  34   m_DistTol = 0;
  35   m_AngleTol = deg2rad(10);
  36   m_BC = 0;
  37 }
  38
  39 vec3_t FillPlane::toPlane(vec3_t x)
  40 {
  41   x -= m_X0;
  42   x -= (x*m_N)*m_N;
  43   return vec3_t(x*m_G1, x*m_G2, 0);
  44 }
  45
  46 vec3_t FillPlane::fromPlane(vec3_t x)
  47 {
  48   return m_X0 + x[0]*m_G1 + x[1]*m_G2;
  49 }
  50
  51 bool FillPlane::isWithinTolerance(vec3_t x)
  52 {
  53   if (fabs((x - m_X0)*m_N) < m_DistTol) {
  54     return true;
  55   }
  56   return false;
  57 }
  58
  59 void FillPlane::createEdgesOnPlane(vtkUnstructuredGrid *edge_grid)
  60 {
  61   vtkIdType num_edges = 0;
  62   vtkIdType num_nodes = 0;
  63
  64   QVector<bool> is_edge_node(m_Grid->GetNumberOfPoints(), false);
  65   for (vtkIdType id_face = 0; id_face < m_Grid->GetNumberOfCells(); ++id_face) {
  66     vtkIdType num_pts, *pts;
  67     if (isSurface(id_face, m_Grid)) {
  68       m_Grid->GetCellPoints(id_face, num_pts, pts);
  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     vtkIdType num_pts, *pts;
 111     if (isSurface(id_face, m_Grid)) {
 112       m_Grid->GetCellPoints(id_face, num_pts, pts);
 113       for (int i = 0; i < num_pts; ++i) {
 114         if (m_Part.c2cGG(id_face, i) == -1) {
 115           vtkIdType id_node1 = pts[i];
 116           vtkIdType id_node2 = pts[0];
 117           if (i < num_pts - 1) {
 118             id_node2 = pts[i + 1];
 119           }
 120           if (is_edge_node[id_node1] && is_edge_node[id_node2]) {
 121             vtkIdType pts[2];
 122             pts[0] = node_map[id_node1];
 123             pts[1] = node_map[id_node2];
 124             edge_grid->InsertNextCell(VTK_LINE, 2, pts);
 125           }
 126         }
 127       }
 128     }
 129   }
 130 }
 131
 132 void FillPlane::closeLoops(vtkUnstructuredGrid *edge_grid)
 133 {
 134   bool done = false;
 135   while (!done) {
 136     QList<vtkIdType> end_nodes;
 137     QVector<int> count(edge_grid->GetNumberOfPoints(), 0);
 138     for (vtkIdType id_edge = 0; id_edge < edge_grid->GetNumberOfCells(); ++id_edge) {
 139       vtkIdType num_pts, *pts;
 140       edge_grid->GetCellPoints(id_edge, num_pts, pts);
 141       for (int i = 0; i < num_pts; ++i) {
 142         ++count[pts[i]];
 143       }
 144     }
 145     for (vtkIdType id_node = 0; id_node < edge_grid->GetNumberOfPoints(); ++id_node) {
 146       if (count[id_node] == 0) {
 147         EG_ERR_RETURN("unable to fill plane(s)");
 148       }
 149       if (count[id_node] > 2) {
 150         EG_ERR_RETURN("unable to fill plane(s)");
 151       }
 152       if (count[id_node] == 1) {
 153         end_nodes.append(id_node);
 154       }
 155     }
 156     if (end_nodes.size() % 2 != 0) {
 157       EG_ERR_RETURN("unable to fill plane(s)");
 158     }
 159     if (end_nodes.size() > 0) {
 160       double dist_min = EG_LARGE_REAL;
 161       vtkIdType id_fill1 = -1;
 162       vtkIdType id_fill2 = -1;
 163       foreach (vtkIdType id_node1, end_nodes) {
 164         vec3_t n1 = m_Part.globalNormal(m_NodeMap[id_node1]);
 165         vec3_t x1;
 166         edge_grid->GetPoint(id_node1, x1.data());
 167         foreach (vtkIdType id_node2, end_nodes) {
 168           if (id_node1 != id_node2) {
 169             vec3_t n2 = m_Part.globalNormal(m_NodeMap[id_node2]);
 170             vec3_t x2;
 171             edge_grid->GetPoint(id_node2, x2.data());
 172             double angle = GeometryTools::angle(n1, -1*n2);
 173             if (angle < m_AngleTol) {
 174               vec3_t d = x2 - x1;
 175               double dist = d.abs();
 176               d.normalise();
 177               if (d*n2 > 0.5 && d*n1 < -0.5) {
 178                 if (dist < dist_min) {
 179                   dist_min = dist;
 180                   id_fill1 = id_node1;
 181                   id_fill2 = id_node2;
 182                 }
 183               }
 184             }
 185           }
 186         }
 187       }
 188       if (id_fill1 == -1 || id_fill2 == -1) {
 189         EG_ERR_RETURN("unable to fill plane(s)");
 190       }
 191       vtkIdType pts[2];
 192       pts[0] = id_fill1;
 193       pts[1] = id_fill2;
 194       edge_grid->InsertNextCell(VTK_LINE, 2, pts);
 195     } else {
 196       done = true;
 197     }
 198   }
 199 }
 200
 201 void FillPlane::gridToPlane(vtkUnstructuredGrid *edge_grid)
 202 {
 203   for (vtkIdType id_node = 0; id_node < edge_grid->GetNumberOfPoints(); ++id_node) {
 204     vec3_t x;
 205     edge_grid->GetPoint(id_node, x.data());
 206     x = toPlane(x);
 207     edge_grid->GetPoints()->SetPoint(id_node, x.data());
 208   }
 209 }
 210
 211 void FillPlane::gridFromPlane(vtkUnstructuredGrid *edge_grid)
 212 {
 213   for (vtkIdType id_node = 0; id_node < edge_grid->GetNumberOfPoints(); ++id_node) {
 214     vec3_t x;
 215     edge_grid->GetPoint(id_node, x.data());
 216     x = fromPlane(x);
 217     edge_grid->GetPoints()->SetPoint(id_node, x.data());
 218   }
 219 }
 220
 221 void FillPlane::triangulate(vtkPolyData *edge_pdata, vtkUnstructuredGrid *tri_grid)
 222 {
 223   EG_VTKSP(vtkDelaunay2D, delaunay);
 224   EG_VTKSP(vtkEgPolyDataToUnstructuredGridFilter, pdata2grid);
 225   delaunay->SetInputData(edge_pdata);
 226   delaunay->SetSourceData(edge_pdata);
 227   delaunay->Update();
 228   pdata2grid->SetInputConnection(delaunay->GetOutputPort());
 229   pdata2grid->Update();
 230   makeCopy(pdata2grid->GetOutput(), tri_grid);
 231   createBasicFields(tri_grid, tri_grid->GetNumberOfCells(), tri_grid->GetNumberOfPoints());
 232   QSet<int> bcs = getAllBoundaryCodes(m_Grid);
 233   m_BC = 0;
 234   foreach (int bc, bcs) {
 235     m_BC = max(m_BC, bc);
 236   }
 237   ++m_BC;
 238   EG_VTKDCC(vtkIntArray, cell_code, tri_grid, "cell_code");
 239   for (vtkIdType id_face = 0; id_face < tri_grid->GetNumberOfCells(); ++id_face) {
 240     cell_code->SetValue(id_face, m_BC);
 241   }
 242 }
 243
 244 void FillPlane::order(vtkUnstructuredGrid *edge_grid, vtkPolyData *edge_pdata)
 245 {
 246   QVector<QVector<vtkIdType> > edges(edge_grid->GetNumberOfCells(), QVector<vtkIdType>(2));
 247   for (vtkIdType id_edge = 0; id_edge < edge_grid->GetNumberOfCells(); ++id_edge) {
 248     vtkIdType num_pts, *pts;
 249     edge_grid->GetCellPoints(id_edge, num_pts, pts);
 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::operate()
 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   EG_VTKSP(vtkUnstructuredGrid, edge_grid);
 357   EG_VTKSP(vtkUnstructuredGrid, tri_grid);
 358   EG_VTKSP(vtkPolyData, edge_pdata);
 359   createEdgesOnPlane(edge_grid);
 360   writeGrid(edge_grid, "filltest1");
 361   closeLoops(edge_grid);
 362   writeGrid(edge_grid, "filltest2");
 363   gridToPlane(edge_grid);
 364   order(edge_grid, edge_pdata);
 365   triangulate(edge_pdata, tri_grid);
 366   gridFromPlane(tri_grid);
 367   MeshPartition tri_part(tri_grid, true);
 368   m_Part.addPartition(tri_part, m_DistTol);
 369   m_Grid->Modified();
 370   CorrectSurfaceOrientation corr_surf;
 371   corr_surf();
 372   GuiMainWindow::pointer()->updateBoundaryCodes(true);
 373 }
 374
 375