src/libengrid/cgaltricadinterface.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 "cgaltricadinterface.h"
  23 #include "meshpartition.h"
  24
  25 CgalTriCadInterface::CgalTriCadInterface(vtkUnstructuredGrid *grid)
  26 {
  27   m_BGrid = vtkUnstructuredGrid::New();
  28   makeCopy(grid, m_BGrid);
  29   m_BPart.setGrid(m_BGrid, true);
  30
  31   double feature_angle;
  32
  33   EG_STOPDATE("2014-08-01");
  34   //getXmlSetting("engrid/surface/settings", "feature_angle", feature_angle);
  35   getSet("surface meshing", "feature angle", 20, feature_angle);
  36   feature_angle = GeometryTools::deg2rad(feature_angle);
  37
  38   // build triangle tree
  39   {
  40     m_Triangles.clear();
  41     QVector<vtkIdType> tris;
  42     getAllCellsOfType(VTK_TRIANGLE, tris, m_BGrid);
  43     m_Triangles.fill(Triangle(), tris.size());
  44     m_Tri2Grid.fill(-1, tris.size());
  45     int i = 0;
  46     foreach (vtkIdType id_cell, tris) {
  47       EG_GET_CELL(id_cell, m_BGrid);
  48       vec3_t a, b, c;
  49       m_BGrid->GetPoint(pts[0], a.data());
  50       m_BGrid->GetPoint(pts[1], b.data());
  51       m_BGrid->GetPoint(pts[2], c.data());
  52       m_Triangles[i] = Triangle(Point(a[0], a[1], a[2]), Point(b[0], b[1], b[2]), Point(c[0], c[1], c[2]));
  53       m_Tri2Grid[i] = id_cell;
  54       ++i;
  55     }
  56     m_TriangleTree.rebuild(m_Triangles.begin(), m_Triangles.end());
  57     m_TriangleTree.accelerate_distance_queries();
  58   }
  59
  60   // build edge tree
  61   {
  62     m_Segments.clear();
  63     QList<Segment> segs;
  64     EG_VTKDCC(vtkIntArray, cell_code, m_BGrid, "cell_code");
  65     EG_FORALL_NODES(id_node1, m_BGrid) {
  66       for (int i = 0; i < m_BPart.n2nGSize(id_node1); ++i) {
  67         vtkIdType id_node2 = m_BPart.n2nGG(id_node1, i);
  68         if (id_node1 < id_node2) {
  69           QVector<vtkIdType> id_face;
  70           m_BPart.getEdgeFaces(id_node1, id_node2, id_face);
  71           if (id_face.size() == 2) {
  72             bool append_edge = false;
  73             if (cell_code->GetValue(id_face[0]) != cell_code->GetValue(id_face[1])) {
  74               append_edge = true;
  75             } else {
  76               vec3_t n1 = GeometryTools::cellNormal(m_BGrid, id_face[0]);
  77               vec3_t n2 = GeometryTools::cellNormal(m_BGrid, id_face[1]);
  78               double alpha = GeometryTools::angle(n1, n2);
  79               if (alpha >= feature_angle) {
  80                 append_edge = true;
  81               }
  82             }
  83             if (append_edge) {
  84               vec3_t x1, x2;
  85               m_BGrid->GetPoint(id_node1, x1.data());
  86               m_BGrid->GetPoint(id_node2, x2.data());
  87               segs.append(Segment(Point(x1[0], x1[1], x1[2]), Point(x2[0], x2[1], x2[2])));
  88             }
  89           }
  90         }
  91       }
  92     }
  93     m_Segments.resize(segs.size());
  94     qCopy(segs.begin(), segs.end(), m_Segments.begin());
  95     m_SegmentTree.rebuild(m_Segments.begin(), m_Segments.end());
  96     m_SegmentTree.accelerate_distance_queries();
  97   }
  98
  99   setName("CgalTriCadInterface");
 100   computeSurfaceCurvature();
 101   //m_ShootRayImplemented = true;
 102 }
 103
 104 void CgalTriCadInterface::computeSurfaceCurvature()
 105 {
 106   QVector<double> node_radius(m_BGrid->GetNumberOfPoints(), EG_LARGE_REAL);
 107   for (vtkIdType id_node = 0; id_node < m_BGrid->GetNumberOfPoints(); ++id_node) {
 108     vec3_t x1;
 109     m_BGrid->GetPoint(id_node, x1.data());
 110     vec3_t n1 = m_BPart.globalNormal(id_node);
 111     for (int i = 0; i < m_BPart.n2nGSize(id_node); ++i) {
 112       vtkIdType id_neigh = m_BPart.n2nGG(id_node, i);
 113       vec3_t n2 = m_BPart.globalNormal(id_neigh);
 114       double scal_prod = max(-1.0, min(1.0, n1*n2));
 115       double alpha = max(1e-3, acos(scal_prod));
 116       if (alpha > 1e-3) {
 117         vec3_t x2;
 118         m_BGrid->GetPoint(id_neigh, x2.data());
 119         double a = (x1 - x2).abs();
 120         node_radius[id_node] = min(node_radius[id_node], a/alpha);
 121       }
 122     }
 123   }
 124
 125   // compute weighted (distance) average of radii
 126   QVector<double> R_new(m_BGrid->GetNumberOfPoints(), 1e99);
 127   for (vtkIdType id_node = 0; id_node < m_BGrid->GetNumberOfPoints(); ++id_node) {
 128     vec3_t x1;
 129     m_BGrid->GetPoint(id_node, x1.data());
 130     int N = m_BPart.n2nGSize(id_node);
 131     QVector<double> L(N);
 132     QVector<double> R(N, -1);
 133     double Lmax = 0;
 134     bool average = false;
 135     for (int i = 0; i < N; ++i) {
 136       vtkIdType id_neigh = m_BPart.n2nGG(id_node, i);
 137       vec3_t x2;
 138       m_BGrid->GetPoint(id_neigh, x2.data());
 139       L[i] = (x2 - x1).abs();
 140       if (node_radius[id_neigh] < 1e90 && L[i] > 0) {
 141         R[i] = node_radius[id_neigh];
 142         Lmax = max(Lmax, L[i]);
 143         average = true;
 144       }
 145     }
 146     if (average) {
 147       R_new[id_node] = 0;
 148       double total_weight = 0;
 149       for (int i = 0; i < N; ++i) {
 150         if (R[i] > 0) {
 151           R_new[id_node] += Lmax/L[i] * R[i];
 152           total_weight += Lmax/L[i];
 153           //R_new[id_node] += R[i];
 154           //total_weight += 1.0;
 155         }
 156       }
 157       R_new[id_node] /= total_weight;
 158     }
 159   }
 160   node_radius = R_new;
 161
 162   m_Radius.fill(0, m_BGrid->GetNumberOfCells());
 163   EG_FORALL_CELLS(id_cell, m_BGrid) {
 164     EG_GET_CELL(id_cell, m_BGrid);
 165     for (int i = 0; i < num_pts; ++i) {
 166       m_Radius[id_cell] = max(m_Radius[id_cell], node_radius[pts[i]]);
 167     }
 168   }
 169 }
 170
 171 CgalTriCadInterface::Ray CgalTriCadInterface::createRay(vec3_t x1, vec3_t v)
 172 {
 173   Point p1(x1[0], x1[1], x1[2]);
 174   vec3_t x2 = x1 + v;
 175   Point p2(x2[0], x2[1], x2[2]);
 176   return Ray(p1, p2);
 177 }
 178
 179 vec3_t CgalTriCadInterface::snap(vec3_t x, bool)
 180 {
 181   Point p(x[0], x[1], x[2]);
 182   TrianglePointAndPrimitiveId result = m_TriangleTree.closest_point_and_primitive(p);
 183   Point cp = result.first;
 184   Triangle* T = result.second;
 185   int id = (T - m_Triangles.begin());
 186   vtkIdType id_face = m_Tri2Grid[id];
 187   m_LastNormal = GeometryTools::cellNormal(m_BGrid, id_face);
 188   m_LastNormal.normalise();
 189   m_LastRadius = m_Radius[id_face];
 190   vec3_t x_snap = vec3_t(cp[0], cp[1], cp[2]);
 191   return x_snap;
 192 }
 193
 194 vec3_t CgalTriCadInterface::snapNode(vtkIdType id_node, vec3_t x, bool correct_curvature)
 195 {
 196   vec3_t x_snap = snap(x, correct_curvature);
 197   vec3_t n_node = m_FPart.globalNormal(id_node);
 198   if (n_node*m_LastNormal < 0) {
 199     double dist_min = EG_LARGE_REAL;
 200     Ray ray = createRay(x, n_node);
 201     std::list<Intersection> intersections;
 202     m_TriangleTree.all_intersections(ray, std::back_inserter(intersections));
 203     for (std::list<Intersection>::iterator i = intersections.begin(); i != intersections.end(); ++i) {
 204       int id = (i->second - m_Triangles.begin());
 205       vtkIdType id_face = m_Tri2Grid[id];
 206       vec3_t n = GeometryTools::cellNormal(m_BGrid, id_face);
 207       n.normalise();
 208       if (n*n_node > 0) {
 209         if (i->first.type() == typeid(Point)) {
 210           Point p = boost::get<Point>(i->first);
 211           vec3_t xs(p[0], p[1], p[2]);
 212           double dist = (x - xs).abs();
 213           if (dist < dist_min) {
 214             dist_min = dist;
 215             x_snap = xs;
 216             m_LastNormal = n;
 217             m_LastRadius = m_Radius[id_face];
 218           }
 219         }
 220       }
 221     }
 222   }
 223   return x_snap;
 224 }
 225
 226 vec3_t CgalTriCadInterface::snapToEdge(vec3_t x)
 227 {
 228   Point p(x[0], x[1], x[2]);
 229   Point cp = m_SegmentTree.closest_point(p);
 230   vec3_t x_snap = vec3_t(cp[0], cp[1], cp[2]);
 231   return x_snap;
 232 }
 233
 234 vec3_t CgalTriCadInterface::snapToCorner(vec3_t)
 235 {
 236   notImplemented();
 237 }
 238
 239 CgalTriCadInterface::HitType CgalTriCadInterface::shootRay(vec3_t x, vec3_t v, vec3_t &x_hit, vec3_t &n_hit, double &r)
 240 {
 241   notImplemented();
 242 }
 243
 244 void CgalTriCadInterface::computeIntersections(vec3_t x, vec3_t v, QVector<QPair<vec3_t, vtkIdType> > &intersections)
 245 {
 246   Ray ray = createRay(x, v);
 247   std::list<Intersection> inters;
 248   m_TriangleTree.all_intersections(ray, std::back_inserter(inters));
 249   intersections.resize(inters.size());
 250   int i_intersections = 0;
 251   for (std::list<Intersection>::iterator i = inters.begin(); i != inters.end(); ++i) {
 252     int id = (i->second - m_Triangles.begin());
 253     if (i->first.type() == typeid(Point)) {
 254       intersections[i_intersections].second = m_Tri2Grid[id];
 255       Point p = boost::get<Point>(i->first);
 256       vec3_t xs(p[0], p[1], p[2]);
 257       intersections[i_intersections].first = xs;
 258     }
 259   }
 260 }