src/libengrid/surfaceoperation.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 #include "surfaceoperation.h"
  22
  23 #include "guimainwindow.h"
  24
  25 #include <vtkCharArray.h>
  26 #include <vtkMath.h>
  27 #include <vtkCellArray.h>
  28 #include <vtkPolygon.h>
  29
  30 #include "geometrytools.h"
  31 #include "meshqualityfaceorientation.h"
  32
  33 using namespace GeometryTools;
  34
  35 SurfaceOperation::SurfaceOperation() : Operation()
  36 {
  37   //default values for determining node types and for smoothing operations
  38   getSet("surface meshing", "edge angle to determine fixed vertices", 180, m_EdgeAngle);
  39   getSet("surface meshing", "feature angle", 30, m_FeatureAngle);
  40   getSet("surface meshing", "threshold for face orientation quality", 0.1, m_FaceOrientationThreshold);
  41   m_FeatureAngle = GeometryTools::deg2rad(m_FeatureAngle);
  42   m_EdgeAngle = GeometryTools::deg2rad(m_EdgeAngle);
  43   setEdgeAngle(m_EdgeAngle);
  44   m_StretchingFactor = 0;
  45   m_UniformSnapPoints = false;
  46
  47   EG_STOPDATE("2015-03-01");
  48   m_StrictFeatureSnap = true;
  49 }
  50
  51 void SurfaceOperation::operate()
  52 {
  53
  54 }
  55
  56 ostream& operator<<(ostream &out, stencil_t S)
  57 {
  58   out << "S.id_cell = " << S.id_cell << " ";
  59   out << "S.id_node = " << S.id_node << " ";
  60   out << "S.sameBC = " << S.sameBC << " ";
  61   out << "S.type = " << S.type_cell << " ";
  62   out << "S.p1 = " << S.p1 << " ";
  63   out << "S.p2 = " << S.p2 << " ";
  64   return(out);
  65 }
  66
  67 stencil_t SurfaceOperation::getStencil(vtkIdType id_cell1, int j1)
  68 {
  69   stencil_t S;
  70   {
  71     vtkIdType N_pts, *pts;
  72     m_Grid->GetCellPoints(id_cell1, N_pts, pts);
  73     S.p1 = pts[j1];
  74     S.p2 = pts[0];
  75     if (j1 < N_pts - 1) {
  76       S.p2 = pts[j1 + 1];
  77     }
  78   }
  79   QSet<vtkIdType> cells_p1;
  80   for (int i = 0; i < m_Part.n2cGSize(S.p1); ++i) {
  81     vtkIdType id_cell = m_Part.n2cGG(S.p1, i);
  82     if (id_cell != id_cell1) {
  83       cells_p1.insert(id_cell);
  84     }
  85   }
  86   QSet<vtkIdType> cells_p2;
  87   for (int i = 0; i < m_Part.n2cGSize(S.p2); ++i) {
  88     vtkIdType id_cell = m_Part.n2cGG(S.p2, i);
  89     if (id_cell != id_cell1) {
  90       cells_p2.insert(id_cell);
  91     }
  92   }
  93   QSet<vtkIdType> cells = cells_p1.intersect(cells_p2);
  94   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  95   S.sameBC = true;
  96   S.id_cell.resize(1);
  97   S.id_cell[0] = id_cell1;
  98   foreach (vtkIdType id_cell, cells) {
  99     if (isSurface(id_cell, m_Grid)) {
 100       S.id_cell.push_back(id_cell);
 101       if (cell_code->GetValue(id_cell) != cell_code->GetValue(id_cell1)) {
 102         S.sameBC = false;
 103       }
 104     }
 105   }
 106   S.id_node.resize(S.id_cell.size());
 107   S.type_cell.resize(S.id_cell.size());
 108   for (int i = 0; i < S.id_cell.size(); ++i) {
 109     vtkIdType N_pts, *pts;
 110     m_Grid->GetCellPoints(S.id_cell[i], N_pts, pts);
 111     S.type_cell[i] = m_Grid->GetCellType(S.id_cell[i]);
 112     for (int j = 0; j < N_pts; ++j) {
 113       if (pts[j] != S.p1 && pts[j] != S.p2) {
 114         S.id_node[i] = pts[j];
 115         break;
 116       }
 117     }
 118   }
 119   return S;
 120 }
 121
 122 int SurfaceOperation::UpdateCurrentMeshDensity()
 123 {
 124   if ( DebugLevel > 0 ) {
 125     cout << "===UpdateMeshDensity START===" << endl;
 126   }
 127   QVector<vtkIdType> cells;
 128   getAllSurfaceCells( cells, m_Grid );
 129   EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
 130   EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
 131   setGrid( m_Grid );
 132   setCells( cells );
 133   if ( DebugLevel > 5 ) {
 134     cout << "cells.size()=" << cells.size() << endl;
 135   }
 136   EG_VTKDCN( vtkDoubleArray, node_meshdensity_current, m_Grid, "node_meshdensity_current" );
 137   l2g_t nodes = getPartNodes();
 138   foreach( vtkIdType node, nodes ) {
 139     node_meshdensity_current->SetValue( node, CurrentMeshDensity( node ) );
 140   }
 141   if ( DebugLevel > 0 ) {
 142     cout << "===UpdateMeshDensity END===" << endl;
 143   }
 144   return( 0 ); ///\todo what for???
 145 }
 146
 147 void SurfaceOperation::readVMD()
 148 {
 149   QString buffer = GuiMainWindow::pointer()->getXmlSection("engrid/surface/table").replace("\n", " ");
 150   int row_count = 0;
 151   int column_count = 0;
 152   m_VMDvector.clear();
 153
 154   if(!buffer.isEmpty()) {
 155     QTextStream in(&buffer, QIODevice::ReadOnly);
 156     in >> row_count >> column_count;
 157     QVector<int> tmp_bcs;
 158     GuiMainWindow::pointer()->getAllBoundaryCodes(tmp_bcs);
 159     {
 160       m_VMDvector.fill(VertexMeshDensity(), row_count);
 161       for (int i = 0; i < row_count; ++i) {
 162         int row, column;
 163         QString formula;
 164         for (int j = 0; j < tmp_bcs.size(); ++j) {
 165           int bc = tmp_bcs[j];
 166           if (j < column_count - 3) {
 167             in >> row >> column >> formula;
 168           } else {
 169             formula = "1";
 170           }
 171           m_VMDvector[row].BCmap[bc] = formula.toInt();
 172         }
 173         in >> row >> column >> formula;
 174         m_VMDvector[row].type = Str2VertexType(formula);
 175         in >> row >> column >> formula;
 176         if (formula == "{{{empty}}}") {
 177           formula = "";
 178         }
 179         m_VMDvector[i].setNodes(formula);
 180         in >> row >> column >> formula;
 181         m_VMDvector[i].density = formula.toDouble();
 182       }
 183     }
 184   }
 185 }
 186
 187 /*
 188 void SurfaceOperation::updateNodeInfo()
 189 {
 190   setAllCells();
 191   readVMD();
 192   l2g_t nodes = getPartNodes();
 193   computeNormals();
 194   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 195   EG_VTKDCN(vtkIntArray, node_type_counter, m_Grid, "node_type_counter");
 196   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 197
 198   bool num_non_simple = 0;
 199   foreach (vtkIdType id_node, nodes) {
 200     if (node_type->GetValue(id_node) != EG_SIMPLE_VERTEX) {
 201       ++num_non_simple;
 202     }
 203   }
 204
 205   foreach (vtkIdType id_node, nodes) {
 206
 207     char old_type = node_type->GetValue(id_node);
 208
 209     if (!m_AngleFeatureDefinition) {
 210       char new_type = getNodeType(id_node, true);
 211       if (old_type == EG_FIXED_VERTEX && new_type != EG_FIXED_VERTEX) {
 212         //EG_BUG;
 213       }
 214
 215       if ((old_type != EG_FEATURE_CORNER_VERTEX && old_type != EG_FEATURE_EDGE_VERTEX) || m_BCodeFeatureDefinition) {
 216         node_type->SetValue(id_node, new_type);
 217       } else {
 218
 219         bool done = false;
 220
 221         // EG_FEATURE_CORNER_VERTEX -> any other type
 222         if (old_type == EG_FEATURE_CORNER_VERTEX && new_type != EG_FEATURE_CORNER_VERTEX) {
 223           if (node_type_counter->GetValue(id_node) > m_TypeProtectionCount) {
 224             node_type->SetValue(id_node, new_type);
 225           } else {
 226             node_type_counter->SetValue(id_node, node_type_counter->GetValue(id_node) + 1);
 227           }
 228           done = true;
 229         }
 230
 231         // EG_FEATURE_EDGE_VERTEX -> EG_SIMPLE_VERTEX
 232         if (old_type == EG_FEATURE_EDGE_VERTEX && new_type == EG_SIMPLE_VERTEX) {
 233           if (node_type_counter->GetValue(id_node) > m_TypeProtectionCount) {
 234             node_type->SetValue(id_node, new_type);
 235           } else {
 236             node_type_counter->SetValue(id_node, node_type_counter->GetValue(id_node) + 1);
 237           }
 238           done = true;
 239         }
 240
 241         if (!done) {
 242           node_type->SetValue(id_node, new_type);
 243           node_type_counter->SetValue(id_node, 0);
 244         }
 245
 246       }
 247     } else {
 248       if (old_type == EG_SIMPLE_VERTEX) {
 249         char new_type = getNodeType(id_node);
 250         if (new_type == EG_FEATURE_CORNER_VERTEX || new_type == EG_FEATURE_EDGE_VERTEX ) {
 251           if (num_non_simple == 0) {
 252             node_type->SetValue(id_node, new_type);
 253           }
 254         } else {
 255           if (num_non_simple == 0) {
 256             node_type->SetValue(id_node, new_type);
 257           }
 258         }
 259       }
 260     }
 261
 262     //density index from table
 263     EG_VTKDCN(vtkIntArray, node_specified_density, m_Grid, "node_specified_density");
 264
 265     VertexMeshDensity nodeVMD = getVMD(id_node);
 266     int idx = nodeVMD.findSmallestVMD(m_VMDvector);
 267     node_specified_density->SetValue(id_node, idx);
 268   }
 269
 270   // mesh quality
 271 //  if (!m_BCodeFeatureDefinition) {
 272 //    MeshQualityFaceOrientation mesh_quality;
 273 //    mesh_quality();
 274 //    EG_VTKDCN(vtkDoubleArray, node_mesh_quality, m_Grid, "node_mesh_quality");
 275 //    EG_FORALL_NODES(id_node, m_Grid) {
 276 //      if (node_mesh_quality->GetValue(id_node) < m_FaceOrientationThreshold) {
 277 //        node_type->SetValue(id_node, EG_SIMPLE_VERTEX);
 278 //      }
 279 //    }
 280 //  }
 281
 282   updatePotentialSnapPoints();
 283 }
 284 */
 285
 286 void SurfaceOperation::updateNodeInfo()
 287 {
 288   setAllCells();
 289   readVMD();
 290   l2g_t nodes = getPartNodes();
 291   computeNormals();
 292   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 293   EG_VTKDCN(vtkIntArray, node_type_counter, m_Grid, "node_type_counter");
 294   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 295
 296   // check if there are non simple vertices
 297   bool all_simple = true;
 298   foreach (vtkIdType id_node, nodes) {
 299     if (node_type->GetValue(id_node) != EG_SIMPLE_VERTEX) {
 300       all_simple = false;
 301       break;
 302     }
 303   }
 304
 305   foreach (vtkIdType id_node, nodes) {
 306
 307     if (all_simple) {
 308       node_type->SetValue(id_node, getNodeType(id_node, true));
 309     }
 310
 311     //density index from table
 312     EG_VTKDCN(vtkIntArray, node_specified_density, m_Grid, "node_specified_density");
 313
 314     VertexMeshDensity nodeVMD = getVMD(id_node);
 315     int idx = nodeVMD.findSmallestVMD(m_VMDvector);
 316     node_specified_density->SetValue(id_node, idx);
 317   }
 318
 319   /*
 320   foreach (vtkIdType id_node, nodes) {
 321
 322     char old_type = node_type->GetValue(id_node);
 323     char new_type = getNodeType(id_node, true);
 324
 325     if (new_type == EG_FIXED_VERTEX) {
 326       if (old_type == EG_SIMPLE_VERTEX) {
 327         node_type->SetValue(id_node, new_type);
 328       }
 329     }
 330
 331     if (new_type > old_type) {
 332       node_type->SetValue(id_node, new_type);
 333     } else {
 334     }
 335
 336     //density index from table
 337     EG_VTKDCN(vtkIntArray, node_specified_density, m_Grid, "node_specified_density");
 338
 339     VertexMeshDensity nodeVMD = getVMD(id_node);
 340     int idx = nodeVMD.findSmallestVMD(m_VMDvector);
 341     node_specified_density->SetValue(id_node, idx);
 342   }
 343   */
 344
 345   updatePotentialSnapPoints();
 346 }
 347
 348
 349 bool SurfaceOperation::checkSnapPointPairForBcMatch(vtkIdType id_node1, vtkIdType id_node2)
 350 {
 351   QSet<int> bcs1, bcs2;
 352   for (int i = 0; i < m_Part.n2bcGSize(id_node1); ++i) {
 353     bcs1.insert(m_Part.n2bcG(id_node1, i));
 354   }
 355   for (int i = 0; i < m_Part.n2bcGSize(id_node2); ++i) {
 356     bcs2.insert(m_Part.n2bcG(id_node2, i));
 357   }
 358   if (bcs2.contains(bcs1)) {
 359     return true;
 360   }
 361   return false;
 362 }
 363
 364 void SurfaceOperation::updatePotentialSnapPoints()
 365 {
 366   setAllSurfaceCells();
 367   l2g_t nodes  = getPartNodes();
 368
 369   m_PotentialSnapPoints.resize(m_Grid->GetNumberOfPoints());
 370
 371   if (m_UniformSnapPoints) {
 372     m_PotentialSnapPoints.resize(m_Grid->GetNumberOfPoints());
 373     foreach( vtkIdType id_node, nodes ) {
 374       for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 375         m_PotentialSnapPoints[id_node].append(m_Part.n2nGG(id_node, i));
 376       }
 377     }
 378     return;
 379   }
 380
 381   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 382
 383   foreach( vtkIdType id_node1, nodes ) {
 384     m_PotentialSnapPoints[id_node1].clear();
 385     char type1 = node_type->GetValue(id_node1);
 386     if (type1 != EG_FIXED_VERTEX) { // fixed vertices do not have any snap-points
 387       QSet<vtkIdType> exclude_nodes;
 388       if (type1 == EG_FEATURE_EDGE_VERTEX || type1 == EG_BOUNDARY_EDGE_VERTEX) {
 389         for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
 390           vtkIdType id_node2 = m_Part.n2nGG(id_node1, i);
 391           char type2 = node_type->GetValue(id_node2);
 392           if (type2 == EG_FEATURE_CORNER_VERTEX || type2 == EG_FIXED_VERTEX) {
 393             for (int j = 0; j < m_Part.n2nGSize(id_node2); ++j) {
 394               vtkIdType id_node3 = m_Part.n2nGG(id_node2, j);
 395               exclude_nodes.insert(id_node3);
 396             }
 397           }
 398         }
 399       }
 400       for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
 401         vtkIdType id_node2 = m_Part.n2nGG(id_node1, i);
 402         char type2 = node_type->GetValue(id_node2);
 403         /*
 404         if (m_StrictFeatureSnap) {
 405           if (   (type1 == EG_SIMPLE_VERTEX)
 406                  || (type1 == EG_FEATURE_EDGE_VERTEX && (type2 == EG_FEATURE_EDGE_VERTEX || type2 == EG_FEATURE_CORNER_VERTEX))
 407                  || (type1 == EG_FEATURE_CORNER_VERTEX && type2 == EG_FEATURE_CORNER_VERTEX)
 408                  || (type1 == EG_BOUNDARY_EDGE_VERTEX && (type2 == EG_BOUNDARY_EDGE_VERTEX || type2 == EG_FIXED_VERTEX)))
 409           {
 410             if (!exclude_nodes.contains(id_node2)) {
 411               if (checkSnapPointPairForBcMatch(id_node1, id_node2)) {
 412                 m_PotentialSnapPoints[id_node1].append(id_node2);
 413               }
 414             }
 415           }
 416         */
 417         if (isFeatureNode(id_node1)) {
 418           char type2 = node_type->GetValue(id_node2);
 419           if (type1 == EG_FEATURE_EDGE_VERTEX && type2 >= type1) {
 420             QVector<vtkIdType> edge_faces;
 421             m_Part.getEdgeFaces(id_node1, id_node2, edge_faces);
 422             if (edge_faces.size() == 2) {
 423               vec3_t n1 = cellNormal(m_Grid, edge_faces[0]);
 424               vec3_t n2 = cellNormal(m_Grid, edge_faces[1]);
 425               if (GeometryTools::angle(n1, n2) >= m_FeatureAngle) {
 426                 m_PotentialSnapPoints[id_node1].append(id_node2);
 427               }
 428             }
 429           }
 430         } else {
 431           if (   (type1 == EG_SIMPLE_VERTEX)
 432                  || (type1 == EG_FEATURE_EDGE_VERTEX)
 433                  || (type1 == EG_BOUNDARY_EDGE_VERTEX && (type2 == EG_BOUNDARY_EDGE_VERTEX || type2 == EG_FIXED_VERTEX)))
 434           {
 435             if (checkSnapPointPairForBcMatch(id_node1, id_node2)) {
 436               m_PotentialSnapPoints[id_node1].append(id_node2);
 437             }
 438           }
 439         }
 440       }
 441
 442       // make sure feature edge vertices have at least two snap points ...
 443       /*
 444     if (type1 == EG_FEATURE_EDGE_VERTEX) {
 445       if (m_PotentialSnapPoints[id_node1].size() < 2) {
 446         m_PotentialSnapPoints[id_node1].clear();
 447       }
 448     }
 449     */
 450     }
 451   }
 452 }
 453
 454 double SurfaceOperation::edgeAngle(vtkIdType id_node1, vtkIdType id_node2)
 455 {
 456   QVector<vtkIdType> faces;
 457   m_Part.getEdgeFaces(id_node1, id_node2, faces);
 458   if (faces.size() == 2) {
 459     vec3_t n1 = cellNormal(m_Grid, faces[0]);
 460     vec3_t n2 = cellNormal(m_Grid, faces[1]);
 461     return GeometryTools::angle(n1, n2);
 462   }
 463   return 0;
 464 }
 465
 466 /*
 467 char SurfaceOperation::getNodeType(vtkIdType id_node, bool fix_unselected)
 468 {
 469   char type = EG_SIMPLE_VERTEX;
 470
 471   if (m_Part.n2nGSize(id_node) == 0) {
 472     return EG_FIXED_VERTEX;
 473   }
 474
 475   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 476
 477   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 478     vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 479
 480     // fix all vertices that are part of a volume element, quad, or polygon
 481     if (m_Grid->GetCellType(id_cell) != VTK_TRIANGLE) {
 482       type = max(type, char(EG_FIXED_VERTEX));
 483     }
 484
 485     int bc = cell_code->GetValue(id_cell);
 486
 487     // fix nodes which belong to faces with unselected boundary codes
 488     if (!m_BoundaryCodes.contains(bc) && fix_unselected) {
 489       type = max(type, char(EG_FIXED_VERTEX));
 490     }
 491
 492   }
 493
 494   int num_bcs = m_Part.n2bcGSize(id_node);
 495   if (num_bcs >= 3 || num_bcs == 0) {
 496     type = max(type, char(EG_FIXED_VERTEX));
 497   }
 498   if (num_bcs == 2) {
 499     QList<vtkIdType> neigh;
 500     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 501       vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 502       if (getEdgeType(id_node, id_neigh) == EG_BOUNDARY_EDGE_VERTEX) {
 503         neigh << id_neigh;
 504       }
 505     }
 506     if (neigh.size() == 2) {
 507       if (M_PI - GeometryTools::angle(m_Grid, neigh[0], id_node, neigh[1]) >= m_EdgeAngle) {
 508         type = max(type, char(EG_FIXED_VERTEX));
 509       }
 510     }
 511     type = max(type, char(EG_BOUNDARY_EDGE_VERTEX));
 512   }
 513
 514   // features are defined via angles
 515   QList<vtkIdType> feature_neigh;
 516   for (int i_neigh = 0; i_neigh < m_Part.n2nGSize(id_node); ++i_neigh) {
 517     vtkIdType id_neigh = m_Part.n2nGG(id_node, i_neigh);
 518     QVector<vtkIdType> id_face;
 519     m_Part.getEdgeFaces(id_node, id_neigh, id_face);
 520     if (id_face.size() == 2) {
 521       vec3_t n1    = cellNormal(m_Grid, id_face[0]);
 522       vec3_t n2    = cellNormal(m_Grid, id_face[1]);
 523       double alpha = GeometryTools::angle(n1, n2);
 524       if (alpha >= m_FeatureAngle) {
 525         feature_neigh.append(id_neigh);
 526       }
 527     }
 528   }
 529   if (feature_neigh.size() == 0) {
 530     type = max(type, char(EG_SIMPLE_VERTEX));
 531   } else if (feature_neigh.size() == 2) {
 532
 533     vec3_t x1, x2, xc;
 534     m_Grid->GetPoint(id_node, xc.data());
 535     m_Grid->GetPoint(feature_neigh[0], x1.data());
 536     m_Grid->GetPoint(feature_neigh[1], x2.data());
 537     double alpha1 = edgeAngle(id_node, feature_neigh[0]);
 538     double alpha2 = edgeAngle(id_node, feature_neigh[1]);
 539     if (alpha1 >= m_FeatureAngle && alpha2 >= m_FeatureAngle) {
 540       if (GeometryTools::angle(xc - x1, x2 - xc) > m_EdgeAngle) {
 541         type = max(type, char(EG_FIXED_VERTEX));
 542       }
 543     }
 544     type = max(type, char(EG_FEATURE_EDGE_VERTEX));
 545   } else {
 546     type = max(type, char(EG_FIXED_VERTEX));
 547     //return EG_FIXED_VERTEX; //EG_FEATURE_CORNER_VERTEX;
 548   }
 549
 550   return type;
 551 }
 552 */
 553
 554 char SurfaceOperation::getNodeType(vtkIdType id_node, bool fix_unselected)
 555 {
 556   l2g_t  nodes = getPartNodes();
 557   g2l_t _nodes = getPartLocalNodes();
 558   l2l_t  n2n   = getPartN2N();
 559
 560   // fix all vertices that are part of a volume element or a quad
 561   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 562     if (m_Grid->GetCellType(m_Part.n2cGG(id_node, i)) != VTK_TRIANGLE) {
 563       return EG_FIXED_VERTEX;
 564     }
 565   }
 566
 567   //initialize default value
 568   char type = EG_SIMPLE_VERTEX;
 569
 570   //loop through edges around id_node
 571
 572   QVector <vtkIdType> edges;
 573
 574   double CosEdgeAngle = cos(this->m_EdgeAngle);
 575
 576   foreach (int i_node2, n2n[_nodes[id_node]]) {
 577     vtkIdType id_node2 = nodes[i_node2];
 578     //-----------------------
 579     //determine edge type
 580     char edge = getEdgeType(id_node2, id_node, fix_unselected);
 581
 582     //-----------------------
 583     //determine node type pre-processing (count nb of complex edges if the node is complex, otherwise, just count the nb of edges)
 584     if (edge && type == EG_SIMPLE_VERTEX) {
 585       edges.clear();
 586       edges.push_back( id_node2 );
 587       type = edge;
 588     } else if (( edge && type == EG_BOUNDARY_EDGE_VERTEX) ||
 589                ( edge && type == EG_FEATURE_EDGE_VERTEX) ||
 590                (!edge && type == EG_SIMPLE_VERTEX)) {
 591       edges.push_back(id_node2);
 592       if (type && edge == EG_BOUNDARY_EDGE_VERTEX) {
 593         type = EG_BOUNDARY_EDGE_VERTEX;//VTK_BOUNDARY_EDGE_VERTEX has priority over VTK_FEATURE_EDGE_VERTEX
 594       }
 595     }
 596   }
 597   //-----------------------
 598   //determine node type post-processing
 599   if (type == EG_FEATURE_EDGE_VERTEX || type == EG_BOUNDARY_EDGE_VERTEX) { //see how many edges; if two, what the angle is
 600
 601     if (edges.size() == 2) { //check angle between edges
 602       double x1[3], x2[3], x3[3], l1[3], l2[3];
 603       m_Grid->GetPoint( edges[0], x1 );
 604       m_Grid->GetPoint( id_node, x2 );
 605       m_Grid->GetPoint( edges[1], x3 );
 606       for ( int k = 0; k < 3; k++ ) {
 607         l1[k] = x2[k] - x1[k];
 608         l2[k] = x3[k] - x2[k];
 609       }
 610       if ( vtkMath::Normalize( l1 ) >= 0.0 &&
 611            vtkMath::Normalize( l2 ) >= 0.0 &&
 612            vtkMath::Dot( l1, l2 ) < CosEdgeAngle ) {
 613              type = EG_FIXED_VERTEX;
 614       }
 615     } else {
 616       type = EG_FIXED_VERTEX;
 617     }
 618
 619     //if along edge
 620   } //if edge vertex
 621
 622   return type;
 623 }
 624
 625
 626
 627 int SurfaceOperation::getEdgeCells(vtkIdType id_node1, vtkIdType id_node2, QVector <vtkIdType> &EdgeCells)
 628 {
 629   g2l_t _nodes = getPartLocalNodes();
 630   l2g_t cells  = getPartCells();
 631   l2l_t n2c    = getPartN2C();
 632
 633   QSet<vtkIdType> S1;
 634   foreach (int i, n2c[_nodes[id_node1]]) {
 635     S1.insert(cells[i]);
 636   }
 637
 638   QSet<vtkIdType> S2;
 639   foreach( int i, n2c[_nodes[id_node2]] ) {
 640     S2.insert(cells[i]);
 641   }
 642
 643   S2.intersect(S1);
 644   EdgeCells = set2Vector(S2, false);
 645   return EdgeCells.size();
 646 }
 647
 648 int SurfaceOperation::getEdgeCells( vtkIdType id_node1, vtkIdType id_node2, QSet <vtkIdType> &EdgeCells )
 649 {
 650   g2l_t _nodes = getPartLocalNodes();
 651   l2g_t cells  = getPartCells();
 652   l2l_t n2c    = getPartN2C();
 653
 654   QSet<vtkIdType> S1;
 655   foreach( int i, n2c[_nodes[id_node1]] ) {
 656     S1.insert( cells[i] );
 657   }
 658
 659   QSet<vtkIdType> S2;
 660   foreach( int i, n2c[_nodes[id_node2]] ) {
 661     S2.insert( cells[i] );
 662   }
 663
 664   EdgeCells = S2.intersect( S1 );
 665   return EdgeCells.size();
 666 }
 667
 668 char SurfaceOperation::getEdgeType(vtkIdType id_node1, vtkIdType id_node2, bool fix_unselected)
 669 {
 670   // compute number of cells around edge [a_node,p2] and put them into neighbour_cells
 671   QVector <vtkIdType> neighbour_cells;
 672   int num_neigh = getEdgeCells(id_node1, id_node2, neighbour_cells) - 1;
 673
 674   // set default value
 675   char edge = EG_SIMPLE_VERTEX;
 676
 677   if (num_neigh == 0) {
 678     edge = EG_BOUNDARY_EDGE_VERTEX;
 679   } else if (num_neigh >= 2) {
 680     edge = EG_BOUNDARY_EDGE_VERTEX;
 681   } else if (num_neigh == 1) {
 682
 683     /*
 684     if (m_Part.isFeatureEdge(id_node1, id_node2, m_FeatureAngle)) {
 685       edge = EG_FEATURE_EDGE_VERTEX;
 686     }
 687     */
 688     vec3_t n1 = cellNormal(m_Grid, neighbour_cells[0]);
 689     vec3_t n2 = cellNormal(m_Grid, neighbour_cells[1]);
 690     if (GeometryTools::angle(n1, n2) > m_FeatureAngle) {
 691       edge = EG_FEATURE_EDGE_VERTEX;
 692     }
 693
 694     // check the boundary codes
 695     EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
 696     int cell_code_0 = cell_code->GetValue( neighbour_cells[0] );
 697     int cell_code_1 = cell_code->GetValue( neighbour_cells[1] );
 698     if (cell_code_0 !=  cell_code_1) {
 699       edge = EG_BOUNDARY_EDGE_VERTEX;
 700     }
 701
 702     if (fix_unselected) {
 703       if (!m_BoundaryCodes.contains(cell_code_0) || !m_BoundaryCodes.contains(cell_code_1)) {
 704         // does not make sense, but should make the points of the edge fixed
 705         edge = EG_FIXED_VERTEX;
 706       }
 707     }
 708   }
 709
 710   return edge;
 711 }
 712
 713 VertexMeshDensity SurfaceOperation::getVMD( vtkIdType id_node )
 714 {
 715   g2l_t _nodes = getPartLocalNodes();
 716   l2g_t  cells = getPartCells();
 717   l2l_t  n2c   = getPartN2C();
 718
 719   EG_VTKDCN( vtkCharArray, node_type, m_Grid, "node_type" );
 720   EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
 721
 722   VertexMeshDensity VMD;
 723   VMD.type = node_type->GetValue( id_node );
 724   VMD.density = 0;
 725   VMD.CurrentNode = id_node;
 726
 727   foreach( int i_cell, n2c[_nodes[id_node]] ) {
 728     vtkIdType id_cell = cells[i_cell];
 729     VMD.BCmap[cell_code->GetValue( id_cell )] = 2;
 730   }
 731   return( VMD );
 732 }
 733
 734 //////////////////////////////////////////////
 735 double SurfaceOperation::currentVertexAvgDist( vtkIdType id_node )
 736 {
 737   l2g_t  nodes = getPartNodes();
 738   g2l_t _nodes = getPartLocalNodes();
 739   l2l_t  n2n   = getPartN2N();
 740
 741   double total_dist = 0;
 742   double avg_dist = 0;
 743   int N = n2n[_nodes[id_node]].size();
 744   vec3_t C;
 745   m_Grid->GetPoint( id_node, C.data() );
 746   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 747     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 748     vec3_t M;
 749     m_Grid->GetPoint( id_node_neighbour, M.data() );
 750     total_dist += ( M - C ).abs();
 751   }
 752   avg_dist = total_dist / ( double )N;
 753   return( avg_dist );
 754 }
 755
 756 double SurfaceOperation::CurrentMeshDensity( vtkIdType id_node )
 757 {
 758   return 1.0 / currentVertexAvgDist( id_node );
 759 }
 760
 761 ///\todo change meshdensity fields to edgelength fields since this is what is mostly used?
 762
 763 /// desired edge length for id_node
 764 double SurfaceOperation::desiredEdgeLength( vtkIdType id_node )
 765 {
 766   EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
 767   return( 1.0 / characteristic_length_desired->GetValue( id_node ) );
 768 }
 769
 770 /// mean desired edge length for id_cell
 771 double SurfaceOperation::meanDesiredEdgeLength( vtkIdType id_cell )
 772 {
 773   vtkIdType num_pts, *pts;
 774   m_Grid->GetCellPoints( id_cell, num_pts, pts );
 775   int total = 0;
 776   for ( int i = 0; i < num_pts; i++ ) {
 777     total += desiredEdgeLength( pts[i] );
 778   }
 779   return total / ( double )num_pts;
 780 }
 781
 782 QVector <vtkIdType> SurfaceOperation::getPotentialSnapPoints( vtkIdType id_node )
 783 {
 784   if ((id_node < 0) || (id_node >= m_PotentialSnapPoints.size())) {
 785     // UpdatePotentialSnapPoints should probably be called before using this function.
 786     EG_BUG;
 787   }
 788   return m_PotentialSnapPoints[id_node];
 789 }
 790
 791 bool SurfaceOperation::isCell(vtkIdType id_node1, vtkIdType id_node2, vtkIdType id_node3)
 792 {
 793   QVector <vtkIdType> EdgeCells_12;
 794   QVector <vtkIdType> EdgeCells_13;
 795   QVector <vtkIdType> inter;
 796
 797   getEdgeCells( id_node1, id_node2, EdgeCells_12 );
 798   getEdgeCells( id_node1, id_node3, EdgeCells_13 );
 799   qcontIntersection( EdgeCells_12, EdgeCells_13, inter );
 800   if(inter.size()>1) {
 801     qWarning()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
 802     qWarning()<<"EdgeCells_12="<<EdgeCells_12;
 803     qWarning()<<"EdgeCells_13="<<EdgeCells_13;
 804     qWarning()<<"inter="<<inter;
 805     writeGrid(m_Grid, "abort");
 806     EG_BUG;// multiple cells in the same place
 807   }
 808   if(DebugLevel>100) {
 809     qDebug()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
 810     qDebug()<<"EdgeCells_12="<<EdgeCells_12;
 811     qDebug()<<"EdgeCells_13="<<EdgeCells_13;
 812     qDebug()<<"inter="<<inter;
 813   }
 814   return(inter.size()>0);
 815 }
 816
 817 void SurfaceOperation::computeNormals()
 818 {
 819   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 820   m_NodeNormal.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
 821   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 822     if (m_Part.localNode(id_node) != -1) {
 823       QSet<int> bcs;
 824       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 825         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 826         if (isSurface(id_cell, m_Grid)) {
 827           int bc = cell_code->GetValue(id_cell);
 828           if (m_BoundaryCodes.contains(bc)) {
 829             bcs.insert(bc);
 830           }
 831         }
 832       }
 833       int num_bcs = bcs.size();
 834       QVector<vec3_t> normal(num_bcs, vec3_t(0,0,0));
 835       QMap<int,int> bcmap;
 836       int i_bc = 0;
 837       foreach (int bc, bcs) {
 838         bcmap[bc] = i_bc;
 839         ++i_bc;
 840       }
 841       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 842         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 843         if (isSurface(id_cell, m_Grid)) {
 844           int bc = cell_code->GetValue(id_cell);
 845           if (m_BoundaryCodes.contains(bc)) {
 846             vtkIdType N_pts, *pts;
 847             m_Grid->GetCellPoints(id_cell, N_pts, pts);
 848             vec3_t a, b, c;
 849             for (int j = 0; j < N_pts; ++j) {
 850               if (pts[j] == id_node) {
 851                 m_Grid->GetPoint(pts[j], a.data());
 852                 if (j > 0) {
 853                   m_Grid->GetPoint(pts[j-1], b.data());
 854                 } else {
 855                   m_Grid->GetPoint(pts[N_pts-1], b.data());
 856                 }
 857                 if (j < N_pts - 1) {
 858                   m_Grid->GetPoint(pts[j+1], c.data());
 859                 } else {
 860                   m_Grid->GetPoint(pts[0], c.data());
 861                 }
 862               }
 863             }
 864             vec3_t u = b - a;
 865             vec3_t v = c - a;
 866             double alpha = GeometryTools::angle(u, v);
 867             vec3_t n = u.cross(v);
 868             n.normalise();
 869             if (checkVector(n)) {
 870               normal[bcmap[bc]] -= alpha*n;
 871             }
 872           }
 873         }
 874       }
 875       for (int i = 0; i < num_bcs; ++i) {
 876         normal[i].normalise();
 877       }
 878       if (num_bcs > 0) {
 879         if (num_bcs > 1) {
 880           if (num_bcs == 3) {
 881             for (int i = 0; i < num_bcs; ++i) {
 882               for (int j = i + 1; j < num_bcs; ++j) {
 883                 vec3_t n = normal[i] + normal[j];
 884                 n.normalise();
 885                 m_NodeNormal[id_node] += n;
 886               }
 887             }
 888           } else {
 889             for (int i = 0; i < num_bcs; ++i) {
 890               m_NodeNormal[id_node] += normal[i];
 891             }
 892           }
 893         } else {
 894           m_NodeNormal[id_node] = normal[0];
 895         }
 896         m_NodeNormal[id_node].normalise();
 897       }
 898     }
 899   }
 900 }
 901
 902 double SurfaceOperation::normalIrregularity(vtkIdType id_node)
 903 {
 904   double nirr = 0;
 905   QVector<vec3_t> nc(m_Part.n2cGSize(id_node));
 906   for (int i = 0; i < nc.size(); ++i) {
 907     nc[i] = GeometryTools::cellNormal(m_Grid, m_Part.n2cGG(id_node, i));
 908     nc[i].normalise();
 909   }
 910   for (int i = 0; i < nc.size(); ++i) {
 911     for (int j = i + 1; j < nc.size(); ++j) {
 912       nirr += 1.0 - nc[i]*nc[j];
 913     }
 914   }
 915   return nirr;
 916 }
 917
 918 double SurfaceOperation::getSurfaceDeviation(vtkIdType id_node)
 919 {
 920   vec3_t x;
 921   m_Grid->GetPoint(id_node, x.data());
 922   double d = 0;
 923   for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 924     vec3_t xe;
 925     m_Grid->GetPoint(m_Part.n2nGG(id_node, i), xe.data());
 926     d += (x - xe)*m_NodeNormal[id_node];
 927   }
 928   if (m_Part.n2nGSize(id_node) > 0) {
 929     d *= 1.0/m_Part.n2nGSize(id_node);
 930   }
 931   return d;
 932 }
 933
 934 bool SurfaceOperation::isFeatureNode(vtkIdType id_node)
 935 {
 936   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 937   char type = node_type->GetValue(id_node);
 938   if (type == EG_FEATURE_CORNER_VERTEX || type == EG_FEATURE_EDGE_VERTEX) {
 939     return true;
 940   }
 941   return false;
 942 }