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