src/surfaceprojection.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 "surfaceprojection.h"
  24
  25 #include "beziertriangle.h"
  26 #include <math.h>
  27
  28 #include <vtkUnstructuredGridWriter.h>
  29
  30 ///\todo check orientation of triangle on which we project to avoid projecting on the wrong side in case of close opposite sides of a surface
  31 ///\todo Delete those grids somewhere
  32 SurfaceProjection::SurfaceProjection() : SurfaceAlgorithm() {
  33   m_BGrid = vtkUnstructuredGrid::New();
  34
  35   m_InterpolationGrid = vtkUnstructuredGrid::New();
  36   m_BezierGrid = vtkUnstructuredGrid::New();
  37
  38   this->setGrid(m_BGrid);
  39
  40   getSet("surface meshing", "correct curvature (experimental)", false, m_correctCurvature);
  41
  42   m_ExactMode = 0;
  43 }
  44
  45 SurfaceProjection::~SurfaceProjection() {
  46   m_InterpolationGrid->Delete();
  47   m_BezierGrid->Delete();
  48 }
  49
  50 void SurfaceProjection::writeGridWithNormals(QString filename)
  51 {
  52   vtkDoubleArray *vectors_normals = vtkDoubleArray::New();
  53   vectors_normals->SetName("normals");
  54   vectors_normals->SetNumberOfComponents(3);
  55   vectors_normals->SetNumberOfTuples(m_BGrid->GetNumberOfPoints());
  56
  57   for (vtkIdType id_node = 0; id_node < m_BGrid->GetNumberOfPoints(); ++id_node) {
  58     vec3_t N = m_NodeNormals[id_node];
  59     double n[3];
  60     n[0] = N[0];
  61     n[1] = N[1];
  62     n[2] = N[2];
  63     vectors_normals->InsertTuple(id_node, n);
  64   }
  65
  66   m_BGrid->GetPointData()->SetVectors(vectors_normals);
  67
  68   vectors_normals->Delete();
  69
  70   saveGrid(m_BGrid, filename + "_BGrid_WithNormals");
  71 }
  72
  73 void debug_output(QVector < QPair<vec3_t, vec3_t> > points,   QVector < QPair<vec3_t, vec3_t> > lines)
  74 {
  75   QFile file("debug.vtk");
  76   if (!file.open(QIODevice::WriteOnly | QIODevice::Text))
  77     return;
  78
  79   QTextStream out(&file);
  80   out << "# vtk DataFile Version 2.0" << endl;
  81   out << "Unstructured Grid Example" << endl;
  82   out << "ASCII" << endl;
  83   out << "" << endl;
  84   out << "DATASET UNSTRUCTURED_GRID" << endl;
  85   out << "POINTS " << points.size() + 2*lines.size() << " double" << endl;
  86   for (int i = 0; i < points.size(); i++) {
  87     vec3_t P = points[i].first;
  88     out << P[0] << " " << P[1] << " " << P[2] << endl;
  89   }
  90   for (int i = 0; i < lines.size(); i++) {
  91     vec3_t P1 = lines[i].first;
  92     vec3_t P2 = lines[i].second;
  93     out << P1[0] << " " << P1[1] << " " << P1[2] << endl;
  94     out << P2[0] << " " << P2[1] << " " << P2[2] << endl;
  95   }
  96   out << "" << endl;
  97   out << "CELLS " << 1 + lines.size() << " " << 4 + 3*lines.size() << endl;
  98   out << "3 0 1 2" << endl;
  99   for (int i = 0; i < lines.size(); i++) {
 100     int P1 = points.size() + 2 * i;
 101     int P2 = P1 + 1;
 102     out << 2 << " " << P1 << " " << P2 << endl;
 103   }
 104   out << "" << endl;
 105   out << "CELL_TYPES " << 1 + lines.size() << endl;
 106   out << "5" << endl;
 107   for (int i = 0; i < lines.size(); i++) {
 108     out << "3" << endl;
 109   }
 110   out << "" << endl;
 111   out << "POINT_DATA " << points.size() + 2*lines.size() << endl;
 112   out << "VECTORS normals float" << endl;
 113   for (int i = 0; i < points.size(); i++) {
 114     vec3_t N = points[i].second;
 115     out << N[0] << " " << N[1] << " " << N[2] << endl;
 116   }
 117   for (int i = 0; i < lines.size(); i++) {
 118     out << 0 << " " << 0 << " " << 0 << endl;
 119     out << 0 << " " << 0 << " " << 0 << endl;
 120   }
 121 }
 122
 123 double interpolate(vec2_t M, vec2_t A, vec2_t nA, vec2_t B, vec2_t nB)
 124 {
 125   double x0 = A[0];
 126   double x1 = B[0];
 127   double x = M[0];
 128
 129   double f_x0 = A[1];
 130   double f_x1 = B[1];
 131   double fp_x0 = -nA[0] / nA[1];
 132   double fp_x1 = -nB[0] / nB[1];
 133
 134   mat4_t matrix;
 135   matrix[0][0] = 3 * pow(x0, 2); matrix[0][1] = 2 * x0;      matrix[0][2] = 1;  matrix[0][3] = 0;
 136   matrix[1][0] = 3 * pow(x1, 2); matrix[1][1] = 2 * x1;      matrix[1][2] = 1;  matrix[1][3] = 0;
 137   matrix[2][0] = pow(x0, 3);   matrix[2][1] = pow(x0, 2); matrix[2][2] = x0; matrix[2][3] = 1;
 138   matrix[3][0] = pow(x1, 3);   matrix[3][1] = pow(x1, 2); matrix[3][2] = x1; matrix[3][3] = 1;
 139
 140   mat4_t matrix_inv = matrix.inverse();
 141   vec4_t Y = vec4_t(fp_x0, fp_x1, f_x0, f_x1);
 142   vec4_t coeffs = matrix_inv * Y;
 143
 144   // f(x)=a*x^3 + b*x^2 + c*x + d
 145   double ret = coeffs[0] * pow(x, 3) + coeffs[1] * pow(x, 2) + coeffs[2] * x + coeffs[3];
 146   return ret;
 147 }
 148
 149 vec2_t interpolate_bezier(double t, vec2_t P0, vec2_t P1, vec2_t P2, vec2_t P3)
 150 {
 151   // B(t)=(1-t^3)*P0 + 3*(1-t)^2*t*P1 + 3*(1-t)*t^2*P2 + t^3*P3;
 152   return (1 - pow(t, 3))*P0 + 3*pow((1 - t), 2)*t*P1 + 3*(1 - t)*pow(t, 2)*P2 + pow(t, 3)*P3;
 153 }
 154
 155 vec3_t interpolate_bezier(double t, vec3_t P0, vec3_t P1, vec3_t P2, vec3_t P3)
 156 {
 157   // B(t)=(1-t^3)*P0 + 3*(1-t)^2*t*P1 + 3*(1-t)*t^2*P2 + t^3*P3;
 158   return (1 - pow(t, 3))*P0 + 3*pow((1 - t), 2)*t*P1 + 3*(1 - t)*pow(t, 2)*P2 + pow(t, 3)*P3;
 159 }
 160
 161 vec2_t interpolate_bezier_normals(double t, vec2_t P0, vec2_t N0, vec2_t P3, vec2_t N3)
 162 {
 163   vec2_t P1(N0[1], -N0[0]);
 164   vec2_t P2(-N3[1], N3[0]);
 165   return interpolate_bezier(t, P0, P1, P2, P3);
 166 }
 167
 168 vec3_t interpolate_bezier_normals(double t, vec3_t P0, vec3_t N0, vec3_t P3, vec3_t N3, vec3_t u, vec3_t v) {
 169   double N0x = N0 * u / u.abs2();
 170   double N0y = N0 * v / v.abs2();
 171   double N3x = N3 * u / u.abs2();
 172   double N3y = N3 * v / v.abs2();
 173
 174   vec3_t P0P1 = (N0y) * u + (-N0x) * v;
 175   vec3_t P3P2 = (-N3y) * u + (N3x) * v;
 176
 177   P0P1.normalise();
 178   P3P2.normalise();
 179
 180   vec3_t P1 = P0 + P0P1;
 181   vec3_t P2 = P3 + P3P2;
 182   return interpolate_bezier(t, P0, P1, P2, P3);
 183 }
 184
 185 // coordinate systems:
 186 // 3D:
 187 // global: X,Y,Z -> g_
 188 // local: g1,g2,g3 -> l_
 189 // 2D:
 190 // triangle: g1,g2 -> t_
 191 // middle planes
 192 // plane 1: AI1, g3 -> pm1_
 193 // plane 2: BI2, g3 -> pm2_
 194 // plane 3: CI3, g3 -> pm3_
 195 // orthogonal edge planes
 196 // plane 1: BC, g3 -> poe1_
 197 // plane 2: CA, g3 -> poe2_
 198 // plane 3: AB, g3 -> poe3_
 199 // non-orthogonal edge planes
 200 // plane 1: BC, nI3 -> pnoe1_
 201 // plane 2: CA, nI1 -> pnoe2_
 202 // plane 3: AB, nI2 -> pnoe3_
 203
 204 vec3_t SurfaceProjection::correctCurvature1(int i_tri, vec3_t g_M)
 205 {
 206   Triangle T = m_Triangles[i_tri];
 207   vec3_t g_A = T.m_a;
 208   vec3_t g_B = T.m_b;
 209   vec3_t g_C = T.m_c;
 210
 211   vec3_t l_A(0, 0, 0);
 212   vec3_t l_B(1, 0, 0);
 213   vec3_t l_C(0, 1, 0);
 214   vec3_t l_M = T.global3DToLocal3D(g_M);
 215
 216   vec2_t t_A(0, 0);
 217   vec2_t t_B(1, 0);
 218   vec2_t t_C(0, 1);
 219   vec2_t t_M(l_M[0], l_M[1]);
 220
 221   vec3_t g_nA = m_NodeNormals[T.m_id_a];
 222   vec3_t g_nB = m_NodeNormals[T.m_id_b];
 223   vec3_t g_nC = m_NodeNormals[T.m_id_c];
 224
 225   vec2_t pm1_A(0, 0);
 226   vec2_t pm1_I1(1, 0);
 227
 228   vec2_t pm2_B(0, 0);
 229   vec2_t pm2_I2(1, 0);
 230
 231   vec2_t pm3_C(0, 0);
 232   vec2_t pm3_I3(1, 0);
 233
 234   vec3_t g_g1 = T.m_g1;
 235   vec3_t g_g2 = T.m_g2;
 236   vec3_t g_g3 = T.m_g3;
 237
 238   vec3_t l_g1(1, 0, 0);
 239   vec3_t l_g2(0, 1, 0);
 240   vec3_t l_g3(0, 0, 1);
 241
 242   // calculating intersections
 243   double k1, k2;
 244   if (!intersection(k1, k2, t_A, t_M - t_A, t_B, t_C - t_B)) return(g_M);
 245   vec2_t t_I1 = t_A + k1 * (t_M - t_A);
 246   vec3_t g_nI1 = (1 - k2) * g_nB + k2 * g_nC;
 247   vec2_t pm1_M(1.0 / k1, 0);
 248   vec2_t poe1_I1(k2, 0);
 249   vec2_t pnoe1_I1(k2, 0);
 250
 251   if (!intersection(k1, k2, t_B, t_M - t_B, t_C, t_A - t_C)) return(g_M);
 252   vec2_t t_I2 = t_B + k1 * (t_M - t_B);
 253   vec3_t g_nI2 = (1 - k2) * g_nC + k2 * g_nA;
 254   vec2_t pm2_M(1.0 / k1, 0);
 255   vec2_t poe2_I2(k2, 0);
 256   vec2_t pnoe2_I2(k2, 0);
 257
 258   if (!intersection(k1, k2, t_C, t_M - t_C, t_A, t_B - t_A)) return(g_M);
 259   vec2_t t_I3 = t_C + k1 * (t_M - t_C);
 260   vec3_t g_nI3 = (1 - k2) * g_nA + k2 * g_nB;
 261   vec2_t pm3_M(1.0 / k1, 0);
 262   vec2_t poe3_I3(k2, 0);
 263   vec2_t pnoe3_I3(k2, 0);
 264
 265   // after knowing intersections
 266   vec3_t l_I1(t_I1[0], t_I1[1], 0);
 267   vec3_t l_I2(t_I2[0], t_I2[1], 0);
 268   vec3_t l_I3(t_I3[0], t_I3[1], 0);
 269
 270   vec3_t g_I1 = g_A + T.m_G * l_I1;
 271   vec3_t g_I2 = g_A + T.m_G * l_I2;
 272   vec3_t g_I3 = g_A + T.m_G * l_I3;
 273
 274   vec3_t l_nI1 = T.m_GI * g_nI1;
 275   vec3_t l_nI2 = T.m_GI * g_nI2;
 276   vec3_t l_nI3 = T.m_GI * g_nI3;
 277
 278   vec3_t l_nA = T.m_GI * g_nA;
 279   vec3_t l_nB = T.m_GI * g_nB;
 280   vec3_t l_nC = T.m_GI * g_nC;
 281
 282   vec3_t l_AI1 = l_I1 - l_A;
 283   vec3_t l_BI2 = l_I2 - l_B;
 284   vec3_t l_CI3 = l_I3 - l_C;
 285
 286   vec3_t g_AI1 = g_I1 - g_A;
 287   vec3_t g_BI2 = g_I2 - g_B;
 288   vec3_t g_CI3 = g_I3 - g_C;
 289
 290   vec2_t pm1_nA = projectVectorOnPlane(l_nA, l_AI1, l_g3);
 291   vec2_t pm1_nI1 = projectVectorOnPlane(l_nI1, l_AI1, l_g3);
 292
 293   vec2_t pm2_nB = projectVectorOnPlane(l_nB, l_BI2, l_g3);
 294   vec2_t pm2_nI2 = projectVectorOnPlane(l_nI2, l_BI2, l_g3);
 295
 296   vec2_t pm3_nC = projectVectorOnPlane(l_nC, l_CI3, l_g3);
 297   vec2_t pm3_nI3 = projectVectorOnPlane(l_nI3, l_CI3, l_g3);
 298
 299   // plane 1: BC, nI1 -> pnoe1_
 300   // plane 2: CA, nI2 -> pnoe2_
 301   // plane 3: AB, nI3 -> pnoe3_
 302
 303   vec3_t l_AB = l_B - l_A;
 304   vec3_t l_BC = l_C - l_B;
 305   vec3_t l_CA = l_A - l_C;
 306
 307   ///////////////
 308   vec2_t pnoe1_B(0, 0);
 309   vec2_t pnoe1_C(1, 0);
 310   vec2_t pnoe1_nB = projectVectorOnPlane(l_nB, l_BC, l_nI1);
 311   vec2_t pnoe1_nC = projectVectorOnPlane(l_nC, l_BC, l_nI1);
 312   vec2_t pnoe1_nI1(0, 1);
 313   vec2_t pnoe1_tB = turnRight(pnoe1_nB);
 314   vec2_t pnoe1_tC = turnRight(pnoe1_nC);
 315   ///////////////
 316   vec2_t pnoe2_C(0, 0);
 317   vec2_t pnoe2_A(1, 0);
 318   vec2_t pnoe2_nC = projectVectorOnPlane(l_nC, l_CA, l_nI2);
 319   vec2_t pnoe2_nA = projectVectorOnPlane(l_nA, l_CA, l_nI2);
 320   vec2_t pnoe2_nI2(0, 1);
 321   vec2_t pnoe2_tC = turnRight(pnoe2_nC);
 322   vec2_t pnoe2_tA = turnRight(pnoe2_nA);
 323   ///////////////
 324   vec2_t pnoe3_A(0, 0);
 325   vec2_t pnoe3_B(1, 0);
 326   vec2_t pnoe3_nA = projectVectorOnPlane(l_nA, l_AB, l_nI3);
 327   vec2_t pnoe3_nB = projectVectorOnPlane(l_nB, l_AB, l_nI3);
 328   vec2_t pnoe3_nI3(0, 1);
 329   vec2_t pnoe3_tA = turnRight(pnoe3_nA);
 330   vec2_t pnoe3_tB = turnRight(pnoe3_nB);
 331   ///////////////
 332   vec3_t g_J1;
 333   vec3_t g_J2;
 334   vec3_t g_J3;
 335
 336
 337   // interpolation attempts
 338   double z1 = interpolate(pm1_M, pm1_A, pm1_nA, pm1_I1, pm1_nI1);
 339   double z2 = interpolate(pm2_M, pm2_B, pm2_nB, pm2_I2, pm2_nI2);
 340   double z3 = interpolate(pm3_M, pm3_C, pm3_nC, pm3_I3, pm3_nI3);
 341   double z = z1;//(z1+z2+z3)/3.0;
 342
 343   vec3_t l_X = l_M + z * l_g3;
 344   vec3_t g_X = g_A + T.m_G * l_X;
 345
 346   // returning value
 347   return g_X;
 348 }// end of correctCurvature
 349
 350 vec3_t SurfaceProjection::cylinder(vec3_t center, double radius, vec3_t g_M) {
 351   vec3_t g_P;
 352   g_P[2] = g_M[2];
 353   double L = sqrt(g_M[0] * g_M[0] + g_M[1] * g_M[1]);
 354   g_P[0] = radius * g_M[0] / L;
 355   g_P[1] = radius * g_M[1] / L;
 356   return g_P;
 357 }
 358
 359 vec3_t SurfaceProjection::cylinder(vec3_t center, double radius, int i_tri, vec3_t r)
 360 {
 361   Triangle T = m_Triangles[i_tri];
 362   vec3_t g_A = T.m_a;
 363   vec3_t g_M = g_A + T.m_G * r;
 364   return cylinder(center, radius, g_M);
 365 }
 366
 367 vec3_t SurfaceProjection::project(vec3_t xp, vtkIdType id_node)
 368 {
 369   checkVector(xp);
 370   if (id_node < 0) {
 371     EG_BUG;
 372   }
 373
 374   vec3_t x_proj(1e99, 1e99, 1e99), r_proj(0, 0, 0);
 375   bool x_proj_set = false;
 376
 377   // initilizing booleans
 378   bool on_triangle = false;
 379   bool need_full_search = false;
 380
 381   if (DebugLevel > 90) {
 382     need_full_search = true;
 383   }
 384
 385   if (id_node >= m_ProjTriangles.size()) { //if there is no known triangle on which to project
 386     int old_size = m_ProjTriangles.size();
 387     m_ProjTriangles.resize(m_FGrid->GetNumberOfPoints());
 388     for (int i = old_size; i < m_ProjTriangles.size(); ++i) {
 389       m_ProjTriangles[i] = -1;
 390     }
 391     need_full_search = true;
 392   } else { //if there is a known triangle on which to project
 393     int i_triangles = m_ProjTriangles[id_node];
 394     if (i_triangles < 0) {
 395       need_full_search = true;
 396     } else {
 397       if (i_triangles >= m_Triangles.size()) {
 398         EG_BUG;
 399       }
 400       Triangle T = m_Triangles[i_triangles];
 401       vec3_t xi, ri;
 402       double d;
 403       int side;
 404       bool intersects = T.projectOnTriangle(xp, xi, ri, d, side, true);
 405       if (!intersects || (d > 0.1*T.m_smallest_length)) {
 406         need_full_search = true;
 407       } else {
 408         x_proj = xi; x_proj_set = true;
 409         if(DebugLevel>90) qWarning()<<"before full search T="<<T.m_a<<T.m_b<<T.m_c;
 410         if (x_proj[0] > 1e98) { // should never happen
 411           EG_BUG;
 412         }
 413         r_proj = ri;
 414         on_triangle = intersects;
 415         ++m_NumDirect;
 416       }
 417     }
 418   }
 419   if (DebugLevel > 90) {
 420     qWarning() << "before full search x_proj=" << x_proj;
 421   }
 422   if (need_full_search) {
 423     if (DebugLevel > 90) {
 424       qWarning() << "starting full search. m_Triangles.size()=" << m_Triangles.size();
 425     }
 426     ++m_NumFull;
 427     double d_min = 1e99;
 428     bool first = true;
 429     for (int i_triangles = 0; i_triangles < m_Triangles.size(); ++i_triangles) {
 430       Triangle T = m_Triangles[i_triangles];
 431       double d;
 432       vec3_t xi, ri;
 433       int side;
 434       bool intersects = T.projectOnTriangle(xp, xi, ri, d, side, true);
 435       if (DebugLevel > 90) {
 436         qWarning() << "T=" << T.m_a << T.m_b << T.m_c;
 437         qWarning() << "d=" << d;
 438       }
 439       if (d >= 1e99) {
 440         EG_BUG;
 441       }
 442       if (/*first ||*/ d < d_min) {
 443         x_proj = xi; x_proj_set = true;
 444         if (x_proj[0] > 1e98) { // should never happen
 445           EG_BUG;
 446         }
 447         r_proj = ri;
 448         d_min = d;
 449         m_ProjTriangles[id_node] = i_triangles;
 450         on_triangle = intersects;
 451         first = false;
 452         if (DebugLevel > 90) {
 453           qWarning() << "new triangle: T=" << T.m_a << T.m_b << T.m_c;
 454         }
 455       }
 456     }
 457     if (DebugLevel > 90) {
 458       qWarning() << "full search done";
 459     }
 460     if (!x_proj_set) { // should never happen
 461       checkVector(xp);
 462       qWarning() << "No projection found for point xp=" << xp[0] << xp[1] << xp[2] << endl;
 463       writeGrid(GuiMainWindow::pointer()->getGrid(), "griddump");
 464       EG_BUG;
 465     }
 466   }
 467   if(DebugLevel > 90) {
 468     qWarning() << "final x_proj=" << x_proj;
 469   }
 470   if (m_correctCurvature) {
 471     x_proj = correctCurvature2(m_ProjTriangles[id_node], xp);
 472   }
 473   return x_proj;
 474 }
 475
 476 int SurfaceProjection::getControlPoints_orthogonal(Triangle T, vec3_t& X_011, vec3_t& X_101, vec3_t& X_110, double Lmax) {
 477   vec3_t A = T.m_a;
 478   vec3_t B = T.m_b;
 479   vec3_t C = T.m_c;
 480   vec3_t nA = m_NodeNormals[T.m_id_a];
 481   vec3_t nB = m_NodeNormals[T.m_id_b];
 482   vec3_t nC = m_NodeNormals[T.m_id_c];
 483
 484   if (T.m_Valid) {
 485     X_011 = intersectionOnPlane(T.m_g3, B, nB, C, nC);
 486     X_101 = intersectionOnPlane(T.m_g3, C, nC, A, nA);
 487     X_110 = intersectionOnPlane(T.m_g3, A, nA, B, nB);
 488   } else {
 489     X_011 = 0.5 * (B + C);
 490     X_101 = 0.5 * (C + A);
 491     X_110 = 0.5 * (A + B);
 492   }
 493
 494   if (!checkVector(X_011)) EG_BUG;
 495   if (!checkVector(X_101)) EG_BUG;
 496   if (!checkVector(X_110)) EG_BUG;
 497
 498   limitControlPoints(T, X_011, X_101, X_110);
 499
 500   if (!checkVector(X_011)) EG_BUG;
 501   if (!checkVector(X_101)) EG_BUG;
 502   if (!checkVector(X_110)) EG_BUG;
 503
 504   return(0);
 505 }
 506
 507 int SurfaceProjection::getControlPoints_nonorthogonal(Triangle T, vec3_t& X_011, vec3_t& X_101, vec3_t& X_110, double Lmax) {
 508   vec3_t A = T.m_a;
 509   vec3_t B = T.m_b;
 510   vec3_t C = T.m_c;
 511   vec3_t nA = m_NodeNormals[T.m_id_a];
 512   vec3_t nB = m_NodeNormals[T.m_id_b];
 513   vec3_t nC = m_NodeNormals[T.m_id_c];
 514
 515 //   cout<<"A="<<A<<" B="<<B<<" C="<<C<<endl;
 516
 517   if ((0.5*(nB + nC)).abs2() == 0) EG_BUG;
 518   if ((0.5*(nC + nA)).abs2() == 0) EG_BUG;
 519   if ((0.5*(nA + nB)).abs2() == 0) EG_BUG;
 520
 521   if (T.m_Valid) {
 522     X_011 = intersectionOnPlane(0.5 * (nB + nC), B, nB, C, nC);
 523     X_101 = intersectionOnPlane(0.5 * (nC + nA), C, nC, A, nA);
 524     X_110 = intersectionOnPlane(0.5 * (nA + nB), A, nA, B, nB);
 525   } else {
 526     X_011 = 0.5 * (B + C);
 527     X_101 = 0.5 * (C + A);
 528     X_110 = 0.5 * (A + B);
 529   }
 530
 531   /// \todo make sure nBC,nCA,nAB are not null vectors!!!
 532   vec3_t nBC = 0.5 * (nB + nC);
 533   vec3_t nCA = 0.5 * (nC + nA);
 534   vec3_t nAB = 0.5 * (nA + nB);
 535
 536   if (!checkVector(X_011)) {
 537     qWarning() << X_011 << " = intersectionOnPlane(" << nBC << ", " << B << ", " << nB << ", " << C << ", " << nC << ")";
 538     EG_BUG;
 539   }
 540   if (!checkVector(X_101)) {
 541     qWarning() << X_101 << " = intersectionOnPlane(" << nCA << ", " << C << ", " << nC << ", " << A << ", " << nA << ")";
 542     EG_BUG;
 543   }
 544   if (!checkVector(X_110)) {
 545     qWarning() << X_110 << " = intersectionOnPlane(" << nAB << ", " << A << ", " << nA << ", " << B << ", " << nB << ")";
 546     EG_BUG;
 547   }
 548
 549   limitControlPoints(T, X_011, X_101, X_110);
 550
 551   if (!checkVector(X_011)) EG_BUG;
 552   if (!checkVector(X_101)) EG_BUG;
 553   if (!checkVector(X_110)) EG_BUG;
 554
 555   return(0);
 556 }
 557
 558 int SurfaceProjection::limitControlPoints(Triangle T, vec3_t& X_011, vec3_t& X_101, vec3_t& X_110) {
 559   vec3_t A = T.m_a;
 560   vec3_t B = T.m_b;
 561   vec3_t C = T.m_c;
 562   vec3_t nA = m_NodeNormals[T.m_id_a];
 563   vec3_t nB = m_NodeNormals[T.m_id_b];
 564   vec3_t nC = m_NodeNormals[T.m_id_c];
 565
 566   /*  vec3_t P_011 = projectPointOnEdge(X_011,B,C-B);
 567     vec3_t P_101 = projectPointOnEdge(X_101,C,A-C);
 568     vec3_t P_110 = projectPointOnEdge(X_110,A,B-A);*/
 569
 570   vec3_t P_011 = 0.5 * (B + C);
 571   vec3_t P_101 = 0.5 * (A + C);
 572   vec3_t P_110 = 0.5 * (A + B);
 573
 574 //   Lmax = 1.0*T.m_smallest_length;
 575
 576   double Lmax_011 = (C - B).abs();
 577   double Lmax_101 = (A - C).abs();
 578   double Lmax_110 = (B - A).abs();
 579
 580   double L_011 = (X_011 - P_011).abs();
 581   double L_101 = (X_101 - P_101).abs();
 582   double L_110 = (X_110 - P_110).abs();
 583
 584   checkVector(X_011);
 585   checkVector(P_011);
 586   checkVector(X_101);
 587   checkVector(P_101);
 588   checkVector(X_110);
 589   checkVector(P_110);
 590
 591   if (L_011 > Lmax_011) {
 592 //     qWarning()<<"WARNING: CONTROL POINT RESTRICTED: Lmax="<<Lmax;
 593 //     qWarning()<<"X_011="<<X_011<<"P_011="<<P_011<<"L_011="<<L_011;
 594     X_011 = P_011 + Lmax_011 / L_011 * (X_011 - P_011);
 595   }
 596   if (L_101 > Lmax_101) {
 597 //     qWarning()<<"WARNING: CONTROL POINT RESTRICTED: Lmax="<<Lmax;
 598 //     qWarning()<<"X_101="<<X_101<<"P_101="<<P_101<<"L_101="<<L_101;
 599     X_101 = P_101 + Lmax_101 / L_101 * (X_101 - P_101);
 600   }
 601   if (L_110 > Lmax_110) {
 602 //     qWarning()<<"WARNING: CONTROL POINT RESTRICTED: Lmax="<<Lmax;
 603 //     qWarning()<<"X_110="<<X_110<<"P_110="<<P_110<<"L_110="<<L_110;
 604     X_110 = P_110 + Lmax_110 / L_110 * (X_110 - P_110);
 605   }
 606   return(0);
 607 }
 608
 609 void SurfaceProjection::writeTriangleGrid(QString filename) {
 610   int N_cells = m_Triangles.size();
 611   int N_points = 3 * m_Triangles.size();
 612
 613   qDebug() << "N_cells=" << N_cells;
 614   qDebug() << "N_points=" << N_points;
 615
 616
 617   EG_VTKSP(vtkUnstructuredGrid, triangle_grid);
 618   allocateGrid(triangle_grid , N_cells, N_points);
 619
 620   vtkDoubleArray *vectors_hasneighbour0 = vtkDoubleArray::New();
 621   vectors_hasneighbour0->SetName("hasneighbour0");
 622   vectors_hasneighbour0->SetNumberOfComponents(3);
 623   vectors_hasneighbour0->SetNumberOfTuples(N_points);
 624
 625   vtkDoubleArray *vectors_hasneighbour1 = vtkDoubleArray::New();
 626   vectors_hasneighbour1->SetName("hasneighbour1");
 627   vectors_hasneighbour1->SetNumberOfComponents(3);
 628   vectors_hasneighbour1->SetNumberOfTuples(N_points);
 629
 630   vtkDoubleArray *vectors_hasneighbour2 = vtkDoubleArray::New();
 631   vectors_hasneighbour2->SetName("hasneighbour2");
 632   vectors_hasneighbour2->SetNumberOfComponents(3);
 633   vectors_hasneighbour2->SetNumberOfTuples(N_points);
 634
 635   vtkDoubleArray *vectors_neighbours = vtkDoubleArray::New();
 636   vectors_neighbours->SetName("neighbours");
 637   vectors_neighbours->SetNumberOfComponents(3);
 638   vectors_neighbours->SetNumberOfTuples(N_cells);
 639
 640   int node_count = 0;
 641   int cell_count = 0;
 642   for (int i = 0; i < m_Triangles.size(); i++) {
 643     vtkIdType pts[3];
 644     triangle_grid->GetPoints()->SetPoint(node_count, m_Triangles[i].m_a.data()); pts[0] = node_count; node_count++;
 645     triangle_grid->GetPoints()->SetPoint(node_count, m_Triangles[i].m_b.data()); pts[1] = node_count; node_count++;
 646     triangle_grid->GetPoints()->SetPoint(node_count, m_Triangles[i].m_c.data()); pts[2] = node_count; node_count++;
 647
 648     vec3_t v0, v1, v2;
 649     v0 = getEdgeNormal(m_Triangles[i].m_id_a, m_Triangles[i].m_id_b).normalise();
 650     v1 = getEdgeNormal(m_Triangles[i].m_id_b, m_Triangles[i].m_id_c).normalise();
 651     v2 = getEdgeNormal(m_Triangles[i].m_id_c, m_Triangles[i].m_id_a).normalise();
 652     if (!checkVector(v0) || !checkVector(v1) || !checkVector(v2)) {
 653       qWarning() << "v0=" << v0;
 654       qWarning() << "v1=" << v1;
 655       qWarning() << "v2=" << v2;
 656       EG_BUG;
 657     };
 658     vectors_hasneighbour0->InsertTuple(pts[0], v0.data());
 659     vectors_hasneighbour1->InsertTuple(pts[1], v1.data());
 660     vectors_hasneighbour2->InsertTuple(pts[2], v2.data());
 661
 662     vec3_t neighbours = m_Triangles[i].m_g3;
 663     if (!m_Triangles[i].m_has_neighbour[0]) neighbours += v0;
 664     if (!m_Triangles[i].m_has_neighbour[1]) neighbours += v1;
 665     if (!m_Triangles[i].m_has_neighbour[2]) neighbours += v2;
 666
 667     int cell = triangle_grid->InsertNextCell(VTK_TRIANGLE, 3, pts); cell_count++;
 668     vectors_neighbours->InsertTuple(cell, neighbours.data());
 669
 670   }
 671
 672   triangle_grid->GetPointData()->AddArray(vectors_hasneighbour0);
 673   triangle_grid->GetPointData()->AddArray(vectors_hasneighbour1);
 674   triangle_grid->GetPointData()->AddArray(vectors_hasneighbour2);
 675   triangle_grid->GetCellData()->AddArray(vectors_neighbours);
 676   vectors_hasneighbour0->Delete();
 677   vectors_hasneighbour1->Delete();
 678   vectors_hasneighbour2->Delete();
 679   vectors_neighbours->Delete();
 680
 681   saveGrid(triangle_grid, filename + "_TriangleGrid");
 682
 683 }
 684
 685 void SurfaceProjection::writeInterpolationGrid(QString filename) {
 686   int N_cells = m_BGrid->GetNumberOfCells() + 1 * m_BGrid->GetNumberOfCells();
 687   int N_points = m_BGrid->GetNumberOfPoints() + 3 * m_BGrid->GetNumberOfCells();
 688   allocateGrid(m_InterpolationGrid , N_cells, N_points);
 689   makeCopyNoAlloc(m_BGrid, m_InterpolationGrid);
 690
 691   int N = 10;
 692   int N_cells_per_triangle = (N - 1) * (N - 1);
 693   int N_points_per_triangle = (N * N + N) / 2;
 694   allocateGrid(m_BezierGrid, m_Triangles.size()*N_cells_per_triangle, m_Triangles.size()*N_points_per_triangle);
 695
 696   vtkIdType node_count = m_BGrid->GetNumberOfPoints();
 697   int cell_count = m_BGrid->GetNumberOfCells();
 698
 699   vtkIdType offset = 0;
 700
 701 /*  qWarning()<<"m_Triangles.size()="<<m_Triangles.size();
 702   qWarning()<<"m_BezierTriangles.size()="<<m_BezierTriangles.size();
 703   qWarning()<<"m_BGrid->GetNumberOfCells()="<<m_BGrid->GetNumberOfCells();
 704   qWarning()<<"m_BGrid->GetNumberOfPoints()="<<m_BGrid->GetNumberOfPoints();*/
 705
 706   for (int i_triangles = 0; i_triangles < m_Triangles.size(); ++i_triangles) {
 707     qDebug()<<"i_triangles="<<i_triangles;
 708
 709     BezierTriangle bezier_triangle = m_BezierTriangles[i_triangles];
 710     // add the local grid
 711     EG_VTKSP(vtkUnstructuredGrid, bezier);
 712     bezier_triangle.getBezierSurface(bezier, N);
 713     offset = addGrid(m_BezierGrid, bezier,offset);
 714
 715     vec3_t K1 = bezier_triangle.m_X_011;
 716     vec3_t K2 = bezier_triangle.m_X_101;
 717     vec3_t K3 = bezier_triangle.m_X_110;
 718
 719     vtkIdType idx_K1, idx_K2, idx_K3;
 720     m_InterpolationGrid->GetPoints()->SetPoint(node_count, K1.data()); idx_K1 = node_count; node_count++;
 721     m_InterpolationGrid->GetPoints()->SetPoint(node_count, K2.data()); idx_K2 = node_count; node_count++;
 722     m_InterpolationGrid->GetPoints()->SetPoint(node_count, K3.data()); idx_K3 = node_count; node_count++;
 723
 724     vtkIdType polyline_nonortho[7];
 725     polyline_nonortho[0] = bezier_triangle.m_id_a;
 726     polyline_nonortho[1] = idx_K3;
 727     polyline_nonortho[2] = bezier_triangle.m_id_b;
 728     polyline_nonortho[3] = idx_K1;
 729     polyline_nonortho[4] = bezier_triangle.m_id_c;
 730     polyline_nonortho[5] = idx_K2;
 731     polyline_nonortho[6] = bezier_triangle.m_id_a;
 732
 733 //     for(int i=0;i<7;i++) qWarning()<<"polyline_nonortho["<<i<<"]="<<polyline_nonortho[i];
 734
 735     m_InterpolationGrid->InsertNextCell(4, 7, polyline_nonortho); cell_count++;
 736
 737   }
 738
 739 /*  qWarning()<<"node_count="<<node_count;
 740   qWarning()<<"cell_count="<<cell_count;
 741
 742   qWarning()<<"m_InterpolationGrid->GetNumberOfPoints()="<<m_InterpolationGrid->GetNumberOfPoints();
 743   qWarning()<<"m_InterpolationGrid->GetNumberOfCells()="<<m_InterpolationGrid->GetNumberOfCells();*/
 744
 745   saveGrid(m_InterpolationGrid, filename + "_InterpolationGrid");
 746   saveGrid(m_BezierGrid, filename + "_BezierGrid");
 747   this->writeGrid(m_BGrid, filename + "_BGrid");
 748
 749 }
 750
 751 vec3_t SurfaceProjection::getEdgeNormal(vtkIdType id_node1, vtkIdType id_node2) {
 752   l2l_t  n2n   = getPartN2N();
 753   l2l_t  n2c = getPartN2C();
 754   l2g_t cells = getPartCells();
 755   g2l_t _cells = getPartLocalCells();
 756   l2g_t nodes = getPartNodes();
 757   g2l_t _nodes = getPartLocalNodes();
 758
 759   vec3_t x1, x2;
 760   m_Grid->GetPoints()->GetPoint(id_node1, x1.data());
 761   m_Grid->GetPoints()->GetPoint(id_node2, x2.data());
 762
 763   QVector <vtkIdType> edge_cells;
 764   getEdgeCells(id_node1, id_node2, edge_cells);
 765   vtkIdType id_cell = edge_cells[0];
 766   int side = getSide(id_cell, m_BGrid, id_node1, id_node2);
 767   vtkIdType *pts, N_pts;
 768   m_BGrid->GetCellPoints(id_cell, N_pts, pts);
 769   vec3_t t;
 770   if (pts[side] == id_node1) {
 771     t = x2 - x1;
 772   } else {
 773     t = x1 - x2;
 774   }
 775   vec3_t n = m_Triangles[id_cell].m_g3;
 776   vec3_t Nedge = t.cross(n);
 777 //   qDebug()<<"Nedge="<<Nedge;
 778   return Nedge;
 779 }
 780
 781 vec3_t SurfaceProjection::ellipsoid(vec3_t M) {
 782   vec3_t A = m_center;
 783   vec3_t AM = M - A;
 784   return(A + m_Rx.abs() * AM.normalise());
 785 }
 786
 787 vec3_t SurfaceProjection::ellipse(vec3_t M) {
 788   vec3_t A = m_center;
 789   vec3_t AM = M - A;
 790   double x = AM * m_Rx;
 791   double y = AM * m_Ry;
 792   double z = AM * m_Rz;
 793 //   qWarning()<<x<<y<<z;
 794   vec3_t P = A + x * m_Rx + y * m_Ry;
 795 //   qWarning()<<P;
 796   vec3_t AP = P - A;
 797 //   qWarning()<<AP;
 798 //   qWarning()<<AP.normalise();
 799 //   qWarning()<<m_Rx.abs();
 800 //   qWarning()<<A + m_Rx.abs() * AP.normalise();
 801   if (AP.abs() <= m_Rx.abs()) return P;
 802   else return(A + m_Rx.abs() * AP.normalise());
 803 }
 804
 805 vec3_t SurfaceProjection::rectangle(vec3_t M) {
 806   return vec3_t(0, 0, 0);
 807 }
 808
 809 vec3_t SurfaceProjection::cuboid(vec3_t M) {
 810   return vec3_t(0, 0, 0);
 811 }
 812
 813 vec3_t SurfaceProjection::cylinder(vec3_t M) {
 814   vec3_t A = m_center;
 815   vec3_t u = m_Rz;
 816   double k = (M - A) * u;
 817   vec3_t P = A + k * u;
 818   vec3_t PM = M - P;
 819   return P + m_Rx.abs() * PM.normalise();
 820   /*  A+ku=P;
 821     AP*PM=0;
 822     ku*(M-A-ku)=0;
 823     k*(u*(M-A))-k2*u2=0;
 824
 825     M-*/
 826 }
 827
 828 void SurfaceProjection::updateBackgroundGridInfo() {
 829   EG_VTKDCN(vtkCharArray, node_type, m_BGrid, "node_type");//node type
 830
 831   getAllCells(m_Cells, m_BGrid);
 832   getNodesFromCells(m_Cells, m_Nodes, m_BGrid);
 833
 834   setBoundaryCodes(GuiMainWindow::pointer()->getAllBoundaryCodes());
 835 //   qDebug()<<"getBoundaryCodes()="<<getBoundaryCodes();
 836
 837   setAllCells();
 838   readVMD();
 839
 840   UpdatePotentialSnapPoints(true, false);
 841   l2l_t  c2c   = getPartC2C();
 842   g2l_t _cells = getPartLocalCells();
 843
 844 //   qDebug()<<"getBoundaryCodes()="<<getBoundaryCodes();
 845
 846   QVector<int> m_LNodes(m_Nodes.size());
 847   for (int i = 0; i < m_LNodes.size(); ++i) {
 848     m_LNodes[i] = i;
 849   }
 850
 851   createNodeToNode(m_Cells, m_Nodes, m_LNodes, m_N2N, m_BGrid);
 852
 853   m_EdgeLength.fill(1e99, m_BGrid->GetNumberOfPoints());
 854   foreach(vtkIdType id_node, m_Nodes) {
 855     vec3_t x;
 856     m_BGrid->GetPoints()->GetPoint(id_node, x.data());
 857     foreach(vtkIdType id_neigh, m_N2N[id_node]) {
 858       vec3_t xn;
 859       m_BGrid->GetPoints()->GetPoint(id_neigh, xn.data());
 860       m_EdgeLength[id_node] = min(m_EdgeLength[id_node], (x - xn).abs());
 861     }
 862     if (m_N2N[id_node].size() < 2) {
 863       EG_BUG;
 864     }
 865   }
 866   // create m_Triangles
 867   m_Triangles.resize(m_BGrid->GetNumberOfCells());
 868   for (vtkIdType id_cell = 0; id_cell < m_BGrid->GetNumberOfCells(); ++id_cell) {
 869     m_Triangles[id_cell] = Triangle(m_BGrid, id_cell);
 870     // TODO: Store infos about neighbour cells for each triangle
 871
 872     for (int i = 0; i < 3; i++) {
 873       int i_cell = _cells[id_cell];
 874       if (c2c[i_cell][i] < 0) {
 875         m_Triangles[id_cell].m_has_neighbour[i] = false;
 876       } else {
 877         m_Triangles[id_cell].m_has_neighbour[i] = true;
 878       }
 879     }
 880
 881   }
 882
 883   // compute node normals
 884   m_NodeNormals.resize(m_BGrid->GetNumberOfPoints());
 885 //   QVector <double> NodeAngles(m_NodeNormals.size());
 886   for (vtkIdType id_node = 0; id_node < m_BGrid->GetNumberOfPoints(); ++id_node) {
 887     m_NodeNormals[id_node] = vec3_t(0, 0, 0);
 888 //     NodeAngles[id_node]=0;
 889   }
 890
 891   foreach(Triangle T, m_Triangles) {
 892     double angle_a = GeometryTools::angle(m_BGrid, T.m_id_c, T.m_id_a, T.m_id_b);
 893     double angle_b = GeometryTools::angle(m_BGrid, T.m_id_a, T.m_id_b, T.m_id_c);
 894     double angle_c = GeometryTools::angle(m_BGrid, T.m_id_b, T.m_id_c, T.m_id_a);
 895     if (isnan(angle_a) || isinf(angle_a)) EG_BUG;
 896     if (isnan(angle_b) || isinf(angle_b)) EG_BUG;
 897     if (isnan(angle_c) || isinf(angle_c)) EG_BUG;
 898     if (!checkVector(T.m_g3)) {
 899       qWarning() << "T.m_g3=" << T.m_g3;
 900       EG_BUG;
 901     }
 902     double total_angle = angle_a + angle_b + angle_c;
 903     m_NodeNormals[T.m_id_a] += angle_a * T.m_g3;
 904     m_NodeNormals[T.m_id_b] += angle_b * T.m_g3;
 905     m_NodeNormals[T.m_id_c] += angle_c * T.m_g3;
 906     if (!checkVector(m_NodeNormals[T.m_id_a])) EG_BUG;
 907     if (!checkVector(m_NodeNormals[T.m_id_b])) EG_BUG;
 908     if (!checkVector(m_NodeNormals[T.m_id_c])) EG_BUG;
 909   }
 910
 911 //   qDebug()<<"===STARTING NORMAL CALCULATION===";
 912
 913   for (vtkIdType id_node = 0; id_node < m_BGrid->GetNumberOfPoints(); ++id_node) {
 914 //     qDebug()<<"id_node="<<id_node<<" and node_type="<< VertexType2Str(node_type->GetValue(id_node));
 915 //     qDebug()<<"n2n["<<id_node<<"]="<<n2n[id_node];
 916
 917     // take into account curvature at boundaries
 918     if (false && node_type->GetValue(id_node) == VTK_BOUNDARY_EDGE_VERTEX) {
 919 //       qDebug()<<"looking for edges...";
 920       QVector <vtkIdType> id_snappoints = getPotentialSnapPoints(id_node);
 921 //       qDebug()<<"id_snappoints.size()="<<id_snappoints.size();
 922 //       qDebug()<<"id_snappoints[0]="<<id_snappoints[0];
 923 //       qDebug()<<"id_snappoints[1]="<<id_snappoints[1];
 924       vec3_t x0, x1, x2;
 925       m_Grid->GetPoints()->GetPoint(id_node, x0.data());
 926       m_Grid->GetPoints()->GetPoint(id_snappoints[0], x1.data());
 927       m_Grid->GetPoints()->GetPoint(id_snappoints[1], x2.data());
 928 //       qDebug()<<"x0="<<x0<<" x1="<<x1<<" x="<<x2;
 929
 930 //       t1.cross(n1);
 931 //       t2.cross(n2);
 932
 933       /*      QVector <vtkIdType> edge_cells_1;
 934             getEdgeCells(id_node,id_snappoints[0],edge_cells_1);
 935             vtkIdType id_cell = edge_cells_1[0];
 936             int side = getSide(id_cell,m_BGrid,id_node,id_snappoints[0]);
 937             vtkIdType *pts, N_pts;
 938             m_BGrid->GetCellPoints(id_cell, N_pts, pts);
 939             vec3_t t1;
 940             if(pts[side]==id_node) {
 941               t1 = x1-x0;
 942             }
 943             else {
 944               t1 = x0-x1;
 945             }
 946             vec3_t n1 = m_Triangles[id_cell].m_g3;
 947             vec3_t Nedge1 = t1.cross(n1);
 948             qDebug()<<"Nedge1="<<Nedge1;*/
 949
 950 //       QVector <int> i_cell_vector = n2c[_nodes[id_node]];
 951 //       vtkIdType id_cell = cells[i_cell_vector[0]];
 952
 953       vec3_t Nedge1 = getEdgeNormal(id_node, id_snappoints[0]);
 954       vec3_t Nedge2 = getEdgeNormal(id_node, id_snappoints[1]);
 955       vec3_t N = Nedge1 + Nedge2;//(x0-x1) + (x0-x2);
 956
 957       m_NodeNormals[id_node] = N;
 958 //       qDebug()<<"x0="<<x0[0]<<x0[1]<<x0[2];
 959 //       qDebug()<<"x1="<<x1[0]<<x1[1]<<x1[2];
 960 //       qDebug()<<"x2="<<x2[0]<<x2[1]<<x2[2];
 961     }
 962     m_NodeNormals[id_node].normalise();
 963 //     qDebug()<<"m_NodeNormals["<<id_node<<"]="<<m_NodeNormals[id_node];
 964   }
 965
 966   // get the control points
 967   ///------------------------------
 968   /// UNDER CONSTRUCTION
 969   ///------------------------------
 970
 971   if(false) {
 972     MeshPartition m_BGrid_partition(m_BGrid, true);
 973
 974     for (vtkIdType id_cell = 0; id_cell < m_BGrid->GetNumberOfCells(); ++id_cell) {
 975       qDebug()<<"id_cell="<<id_cell;
 976       Triangle T = m_Triangles[id_cell];
 977       for (int i_edge = 0; i_edge < 3; i_edge++) {
 978
 979         // preparations
 980         vtkIdType id_cell1 = id_cell;
 981         vtkIdType id_cell2 = m_BGrid_partition.c2cGG(id_cell, i_edge);
 982
 983         vtkIdType N_pts1, *pts1;
 984         if (id_cell1 != -1) {
 985           m_BGrid->GetCellPoints(id_cell1, N_pts1, pts1);
 986         }
 987         vtkIdType N_pts2, *pts2;
 988         if (id_cell2 != -1) {
 989           m_BGrid->GetCellPoints(id_cell2, N_pts2, pts2);
 990         }
 991
 992         vtkIdType p1 = pts1[i_edge];
 993         vtkIdType p2 = pts1[(i_edge+1)%N_pts1];
 994
 995         if (!m_ControlPoints.contains(OrderedPair(p1, p2))) {
 996
 997           // calculate control point
 998           vec3_t control_point;
 999           vec3_t X_011, X_101, X_110;
1000           getControlPoints_nonorthogonal(T, X_011, X_101, X_110, 1e99);
1001           if (i_edge == 0) {
1002             control_point = X_110;
1003           } else if (i_edge == 1) {
1004             control_point = X_011;
1005           } else {
1006             control_point = X_101;
1007           }
1008
1009           //id_cell1
1010           if (id_cell1 != -1) {
1011             for (int j_edge = 0; j_edge < 3; j_edge++) {
1012               vtkIdType a = pts1[j_edge];
1013               vtkIdType b = pts1[(j_edge+1)%N_pts1];
1014               vtkIdType c = pts1[(j_edge+2)%N_pts1];
1015               if (a == p1 || a == p2) {
1016                 if (!m_ControlPoints.contains(OrderedPair(b, c))) {
1017                   vec3_t other_control_point = m_ControlPoints[OrderedPair(b, c)];
1018
1019                   vec3_t segment_A, segment_B;
1020                   vec3_t triangle_A, triangle_B, triangle_C;
1021                   m_BGrid->GetPoints()->GetPoint(a, segment_A.data());
1022                   m_BGrid->GetPoints()->GetPoint(b, triangle_A.data());
1023                   m_BGrid->GetPoints()->GetPoint(c, triangle_B.data());
1024                   segment_B = control_point;
1025                   triangle_C = other_control_point;
1026
1027                   vec3_t xi, ri;
1028                   if (intersectEdgeAndTriangle(triangle_A, triangle_B, triangle_C, segment_A, segment_B, xi, ri)) {
1029                     control_point = xi;
1030                   }
1031
1032                 }
1033               }
1034             }
1035           }
1036
1037           //id_cell2
1038           if (id_cell2 != -1) {
1039             for (int j_edge = 0; j_edge < 3; j_edge++) {
1040               vtkIdType a = pts2[j_edge];
1041               vtkIdType b = pts2[(j_edge+1)%N_pts2];
1042               vtkIdType c = pts2[(j_edge+2)%N_pts2];
1043               if (a == p1 || a == p2) {
1044                 if (!m_ControlPoints.contains(OrderedPair(b, c))) {
1045                   vec3_t other_control_point = m_ControlPoints[OrderedPair(b, c)];
1046
1047                   vec3_t segment_A, segment_B;
1048                   vec3_t triangle_A, triangle_B, triangle_C;
1049                   m_BGrid->GetPoints()->GetPoint(a, segment_A.data());
1050                   m_BGrid->GetPoints()->GetPoint(b, triangle_A.data());
1051                   m_BGrid->GetPoints()->GetPoint(c, triangle_B.data());
1052                   segment_B = control_point;
1053                   triangle_C = other_control_point;
1054
1055                   vec3_t xi, ri;
1056                   if (intersectEdgeAndTriangle(triangle_A, triangle_B, triangle_C, segment_A, segment_B, xi, ri)) {
1057                     control_point = xi;
1058                   }
1059
1060                 }
1061               }
1062             }
1063           }
1064
1065           // set final point
1066           m_ControlPoints[OrderedPair(p1, p2)] = control_point;
1067
1068         }// end of if(!m_ControlPoints.contains(OrderedPair(p1,p2)))
1069       }// end of loop through edges
1070     }// end of loop through cells
1071
1072 //     qDebug()<<"=== CONTROL POINTS READY ===";
1073
1074   } else {
1075     for (int i_tri = 0; i_tri < m_Triangles.size(); i_tri++) {
1076       m_Triangles[i_tri].m_Normal_a = m_NodeNormals[m_Triangles[i_tri].m_id_a];
1077       m_Triangles[i_tri].m_Normal_b = m_NodeNormals[m_Triangles[i_tri].m_id_b];
1078       m_Triangles[i_tri].m_Normal_c = m_NodeNormals[m_Triangles[i_tri].m_id_c];
1079
1080       Triangle T = m_Triangles[i_tri];
1081       vec3_t X_200 = T.m_a;
1082       vec3_t X_020 = T.m_b;
1083       vec3_t X_002 = T.m_c;
1084       vec3_t X_011, X_101, X_110;
1085       getControlPoints_nonorthogonal(T, X_011, X_101, X_110, 1e99);
1086       m_ControlPoints[OrderedPair(T.m_id_b, T.m_id_c)] = X_011;
1087       m_ControlPoints[OrderedPair(T.m_id_c, T.m_id_a)] = X_101;
1088       m_ControlPoints[OrderedPair(T.m_id_a, T.m_id_b)] = X_110;
1089     }
1090   }
1091   ///------------------------------
1092
1093   // store the bezier triangles
1094   m_BezierTriangles.resize(m_Triangles.size());
1095   for (int i_tri = 0; i_tri < m_Triangles.size(); i_tri++) {
1096     Triangle T = m_Triangles[i_tri];
1097     vec3_t X_200 = T.m_a;
1098     vec3_t X_020 = T.m_b;
1099     vec3_t X_002 = T.m_c;
1100
1101     if (!m_ControlPoints.contains(OrderedPair(T.m_id_b, T.m_id_c))) EG_BUG;
1102     if (!m_ControlPoints.contains(OrderedPair(T.m_id_c, T.m_id_a))) EG_BUG;
1103     if (!m_ControlPoints.contains(OrderedPair(T.m_id_a, T.m_id_b))) EG_BUG;
1104
1105     vec3_t X_011 = m_ControlPoints[OrderedPair(T.m_id_b, T.m_id_c)];
1106     vec3_t X_101 = m_ControlPoints[OrderedPair(T.m_id_c, T.m_id_a)];
1107     vec3_t X_110 = m_ControlPoints[OrderedPair(T.m_id_a, T.m_id_b)];
1108
1109     m_BezierTriangles[i_tri] = BezierTriangle(X_200, X_020, X_002, X_011, X_101, X_110);
1110     m_BezierTriangles[i_tri].m_id_a = T.m_id_a;
1111     m_BezierTriangles[i_tri].m_id_b = T.m_id_b;
1112     m_BezierTriangles[i_tri].m_id_c = T.m_id_c;
1113     m_BezierTriangles[i_tri].m_has_neighbour = m_Triangles[i_tri].m_has_neighbour;
1114   }
1115
1116 //   qDebug()<<"=== BEZIER TRIANGLES READY ===";
1117
1118   // compute maximum angle per node
1119   QVector<double> min_cos(m_BGrid->GetNumberOfPoints(), 1.0);
1120   foreach(Triangle T, m_Triangles) {
1121     double cosa = T.m_g3 * m_NodeNormals[T.m_id_a];
1122     double cosb = T.m_g3 * m_NodeNormals[T.m_id_b];
1123     double cosc = T.m_g3 * m_NodeNormals[T.m_id_c];
1124     min_cos[T.m_id_a] = min(cosa, min_cos[T.m_id_a]);
1125     min_cos[T.m_id_b] = min(cosb, min_cos[T.m_id_b]);
1126     min_cos[T.m_id_c] = min(cosc, min_cos[T.m_id_c]);
1127   }
1128   for (vtkIdType id_node = 0; id_node < m_BGrid->GetNumberOfPoints(); ++id_node) {
1129     double s = sqrt(1.0 - sqr(min(1 - 1e-20, min_cos[id_node])));
1130     m_EdgeLength[id_node] *= m_RadiusFactor * min_cos[id_node] / s;
1131   }
1132
1133 //   qDebug()<<"=== UPDATE BGRID DONE ===";
1134 }
1135
1136 vec3_t SurfaceProjection::correctCurvature2(int i_tri, vec3_t g_M) {
1137   return m_BezierTriangles[i_tri].projectOnQuadraticBezierTriangle(g_M);
1138 }// end of correctCurvature2