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("2014-07-21");
  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   foreach (vtkIdType id_node, nodes) {
 296
 297     char old_type = node_type->GetValue(id_node);
 298     char new_type = getNodeType(id_node, true);
 299     if (new_type > old_type) {
 300       node_type->SetValue(id_node, new_type);
 301     }
 302
 303     //density index from table
 304     EG_VTKDCN(vtkIntArray, node_specified_density, m_Grid, "node_specified_density");
 305
 306     VertexMeshDensity nodeVMD = getVMD(id_node);
 307     int idx = nodeVMD.findSmallestVMD(m_VMDvector);
 308     node_specified_density->SetValue(id_node, idx);
 309   }
 310
 311   updatePotentialSnapPoints();
 312 }
 313
 314
 315 bool SurfaceOperation::checkSnapPointPairForBcMatch(vtkIdType id_node1, vtkIdType id_node2)
 316 {
 317   QSet<int> bcs1, bcs2;
 318   for (int i = 0; i < m_Part.n2bcGSize(id_node1); ++i) {
 319     bcs1.insert(m_Part.n2bcG(id_node1, i));
 320   }
 321   for (int i = 0; i < m_Part.n2bcGSize(id_node2); ++i) {
 322     bcs2.insert(m_Part.n2bcG(id_node2, i));
 323   }
 324   if (bcs2.contains(bcs1)) {
 325     return true;
 326   }
 327   return false;
 328 }
 329
 330 void SurfaceOperation::updatePotentialSnapPoints()
 331 {
 332   setAllSurfaceCells();
 333   l2g_t nodes  = getPartNodes();
 334
 335   m_PotentialSnapPoints.resize(m_Grid->GetNumberOfPoints());
 336
 337   if (m_UniformSnapPoints) {
 338     m_PotentialSnapPoints.resize(m_Grid->GetNumberOfPoints());
 339     foreach( vtkIdType id_node, nodes ) {
 340       for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 341         m_PotentialSnapPoints[id_node].append(m_Part.n2nGG(id_node, i));
 342       }
 343     }
 344     return;
 345   }
 346
 347   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 348
 349   foreach( vtkIdType id_node1, nodes ) {
 350     m_PotentialSnapPoints[id_node1].clear();
 351     char type1 = node_type->GetValue(id_node1);
 352     if (type1 != EG_FIXED_VERTEX) { // fixed vertices do not have any snap-points
 353       QSet<vtkIdType> exclude_nodes;
 354       if (type1 == EG_FEATURE_EDGE_VERTEX || type1 == EG_BOUNDARY_EDGE_VERTEX) {
 355         for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
 356           vtkIdType id_node2 = m_Part.n2nGG(id_node1, i);
 357           char type2 = node_type->GetValue(id_node2);
 358           if (type2 == EG_FEATURE_CORNER_VERTEX || type2 == EG_FIXED_VERTEX) {
 359             for (int j = 0; j < m_Part.n2nGSize(id_node2); ++j) {
 360               vtkIdType id_node3 = m_Part.n2nGG(id_node2, j);
 361               exclude_nodes.insert(id_node3);
 362             }
 363           }
 364         }
 365       }
 366       for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
 367         vtkIdType id_node2 = m_Part.n2nGG(id_node1, i);
 368         char type2 = node_type->GetValue(id_node2);
 369         /*
 370         if (m_StrictFeatureSnap) {
 371           if (   (type1 == EG_SIMPLE_VERTEX)
 372                  || (type1 == EG_FEATURE_EDGE_VERTEX && (type2 == EG_FEATURE_EDGE_VERTEX || type2 == EG_FEATURE_CORNER_VERTEX))
 373                  || (type1 == EG_FEATURE_CORNER_VERTEX && type2 == EG_FEATURE_CORNER_VERTEX)
 374                  || (type1 == EG_BOUNDARY_EDGE_VERTEX && (type2 == EG_BOUNDARY_EDGE_VERTEX || type2 == EG_FIXED_VERTEX)))
 375           {
 376             if (!exclude_nodes.contains(id_node2)) {
 377               if (checkSnapPointPairForBcMatch(id_node1, id_node2)) {
 378                 m_PotentialSnapPoints[id_node1].append(id_node2);
 379               }
 380             }
 381           }
 382         */
 383         if (isFeatureNode(id_node1)) {
 384           if (type1 == EG_FEATURE_EDGE_VERTEX && isFeatureNode(id_node2)) {
 385             QVector<vtkIdType> edge_faces;
 386             m_Part.getEdgeFaces(id_node1, id_node2, edge_faces);
 387             if (edge_faces.size() == 2) {
 388               vec3_t n1 = cellNormal(m_Grid, edge_faces[0]);
 389               vec3_t n2 = cellNormal(m_Grid, edge_faces[1]);
 390               if (GeometryTools::angle(n1, n2) >= m_FeatureAngle) {
 391                 m_PotentialSnapPoints[id_node1].append(id_node2);
 392               }
 393             }
 394           }
 395         } else {
 396           if (   (type1 == EG_SIMPLE_VERTEX)
 397                  || (type1 == EG_FEATURE_EDGE_VERTEX)
 398                  || (type1 == EG_BOUNDARY_EDGE_VERTEX && (type2 == EG_BOUNDARY_EDGE_VERTEX || type2 == EG_FIXED_VERTEX)))
 399           {
 400             if (checkSnapPointPairForBcMatch(id_node1, id_node2)) {
 401               m_PotentialSnapPoints[id_node1].append(id_node2);
 402             }
 403           }
 404         }
 405       }
 406
 407       // make sure feature edge vertices have at least two snap points ...
 408       /*
 409     if (type1 == EG_FEATURE_EDGE_VERTEX) {
 410       if (m_PotentialSnapPoints[id_node1].size() < 2) {
 411         m_PotentialSnapPoints[id_node1].clear();
 412       }
 413     }
 414     */
 415     }
 416   }
 417 }
 418
 419 double SurfaceOperation::edgeAngle(vtkIdType id_node1, vtkIdType id_node2)
 420 {
 421   QVector<vtkIdType> faces;
 422   m_Part.getEdgeFaces(id_node1, id_node2, faces);
 423   if (faces.size() == 2) {
 424     vec3_t n1 = cellNormal(m_Grid, faces[0]);
 425     vec3_t n2 = cellNormal(m_Grid, faces[1]);
 426     return GeometryTools::angle(n1, n2);
 427   }
 428   return 0;
 429 }
 430
 431 /*
 432 char SurfaceOperation::getNodeType(vtkIdType id_node, bool fix_unselected)
 433 {
 434   char type = EG_SIMPLE_VERTEX;
 435
 436   if (m_Part.n2nGSize(id_node) == 0) {
 437     return EG_FIXED_VERTEX;
 438   }
 439
 440   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 441
 442   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 443     vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 444
 445     // fix all vertices that are part of a volume element, quad, or polygon
 446     if (m_Grid->GetCellType(id_cell) != VTK_TRIANGLE) {
 447       type = max(type, char(EG_FIXED_VERTEX));
 448     }
 449
 450     int bc = cell_code->GetValue(id_cell);
 451
 452     // fix nodes which belong to faces with unselected boundary codes
 453     if (!m_BoundaryCodes.contains(bc) && fix_unselected) {
 454       type = max(type, char(EG_FIXED_VERTEX));
 455     }
 456
 457   }
 458
 459   int num_bcs = m_Part.n2bcGSize(id_node);
 460   if (num_bcs >= 3 || num_bcs == 0) {
 461     type = max(type, char(EG_FIXED_VERTEX));
 462   }
 463   if (num_bcs == 2) {
 464     QList<vtkIdType> neigh;
 465     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 466       vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 467       if (getEdgeType(id_node, id_neigh) == EG_BOUNDARY_EDGE_VERTEX) {
 468         neigh << id_neigh;
 469       }
 470     }
 471     if (neigh.size() == 2) {
 472       if (M_PI - GeometryTools::angle(m_Grid, neigh[0], id_node, neigh[1]) >= m_EdgeAngle) {
 473         type = max(type, char(EG_FIXED_VERTEX));
 474       }
 475     }
 476     type = max(type, char(EG_BOUNDARY_EDGE_VERTEX));
 477   }
 478
 479   // features are defined via angles
 480   QList<vtkIdType> feature_neigh;
 481   for (int i_neigh = 0; i_neigh < m_Part.n2nGSize(id_node); ++i_neigh) {
 482     vtkIdType id_neigh = m_Part.n2nGG(id_node, i_neigh);
 483     QVector<vtkIdType> id_face;
 484     m_Part.getEdgeFaces(id_node, id_neigh, id_face);
 485     if (id_face.size() == 2) {
 486       vec3_t n1    = cellNormal(m_Grid, id_face[0]);
 487       vec3_t n2    = cellNormal(m_Grid, id_face[1]);
 488       double alpha = GeometryTools::angle(n1, n2);
 489       if (alpha >= m_FeatureAngle) {
 490         feature_neigh.append(id_neigh);
 491       }
 492     }
 493   }
 494   if (feature_neigh.size() == 0) {
 495     type = max(type, char(EG_SIMPLE_VERTEX));
 496   } else if (feature_neigh.size() == 2) {
 497
 498     vec3_t x1, x2, xc;
 499     m_Grid->GetPoint(id_node, xc.data());
 500     m_Grid->GetPoint(feature_neigh[0], x1.data());
 501     m_Grid->GetPoint(feature_neigh[1], x2.data());
 502     double alpha1 = edgeAngle(id_node, feature_neigh[0]);
 503     double alpha2 = edgeAngle(id_node, feature_neigh[1]);
 504     if (alpha1 >= m_FeatureAngle && alpha2 >= m_FeatureAngle) {
 505       if (GeometryTools::angle(xc - x1, x2 - xc) > m_EdgeAngle) {
 506         type = max(type, char(EG_FIXED_VERTEX));
 507       }
 508     }
 509     type = max(type, char(EG_FEATURE_EDGE_VERTEX));
 510   } else {
 511     type = max(type, char(EG_FIXED_VERTEX));
 512     //return EG_FIXED_VERTEX; //EG_FEATURE_CORNER_VERTEX;
 513   }
 514
 515   return type;
 516 }
 517 */
 518
 519 char SurfaceOperation::getNodeType(vtkIdType id_node, bool fix_unselected)
 520 {
 521   l2g_t  nodes = getPartNodes();
 522   g2l_t _nodes = getPartLocalNodes();
 523   l2l_t  n2n   = getPartN2N();
 524
 525   // fix all vertices that are part of a volume element or a quad
 526   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 527     if (m_Grid->GetCellType(m_Part.n2cGG(id_node, i)) != VTK_TRIANGLE) {
 528       return EG_FIXED_VERTEX;
 529     }
 530   }
 531
 532   //initialize default value
 533   char type = EG_SIMPLE_VERTEX;
 534
 535   //loop through edges around id_node
 536
 537   QVector <vtkIdType> edges;
 538
 539   double CosEdgeAngle = cos(this->m_EdgeAngle);
 540
 541   foreach( int i_node2, n2n[_nodes[id_node]] ) {
 542     vtkIdType id_node2 = nodes[i_node2];
 543     //-----------------------
 544     //determine edge type
 545     char edge = getEdgeType(id_node2, id_node, fix_unselected);
 546
 547     //-----------------------
 548     //determine node type pre-processing (count nb of complex edges if the node is complex, otherwise, just count the nb of edges)
 549     if ( edge && type == EG_SIMPLE_VERTEX ) {
 550       edges.clear();
 551       edges.push_back( id_node2 );
 552       type = edge;
 553     }
 554     else if (( edge && type == EG_BOUNDARY_EDGE_VERTEX ) ||
 555              ( edge && type == EG_FEATURE_EDGE_VERTEX ) ||
 556              ( !edge && type == EG_SIMPLE_VERTEX ) ) {
 557       edges.push_back( id_node2 );
 558       if ( type && edge == EG_BOUNDARY_EDGE_VERTEX ) {
 559         type = EG_BOUNDARY_EDGE_VERTEX;//VTK_BOUNDARY_EDGE_VERTEX has priority over VTK_FEATURE_EDGE_VERTEX
 560       }
 561     }
 562   }
 563   //-----------------------
 564   //determine node type post-processing
 565   if ( type == EG_FEATURE_EDGE_VERTEX || type == EG_BOUNDARY_EDGE_VERTEX ) { //see how many edges; if two, what the angle is
 566
 567     if (edges.size() == 2) { //check angle between edges
 568       double x1[3], x2[3], x3[3], l1[3], l2[3];
 569       m_Grid->GetPoint( edges[0], x1 );
 570       m_Grid->GetPoint( id_node, x2 );
 571       m_Grid->GetPoint( edges[1], x3 );
 572       for ( int k = 0; k < 3; k++ ) {
 573         l1[k] = x2[k] - x1[k];
 574         l2[k] = x3[k] - x2[k];
 575       }
 576       if ( vtkMath::Normalize( l1 ) >= 0.0 &&
 577            vtkMath::Normalize( l2 ) >= 0.0 &&
 578            vtkMath::Dot( l1, l2 ) < CosEdgeAngle ) {
 579              type = EG_FIXED_VERTEX;
 580       }
 581     } //if along edge
 582   } //if edge vertex
 583
 584   return type;
 585 }
 586
 587
 588
 589 int SurfaceOperation::getEdgeCells(vtkIdType id_node1, vtkIdType id_node2, QVector <vtkIdType> &EdgeCells)
 590 {
 591   g2l_t _nodes = getPartLocalNodes();
 592   l2g_t cells  = getPartCells();
 593   l2l_t n2c    = getPartN2C();
 594
 595   QSet<vtkIdType> S1;
 596   foreach (int i, n2c[_nodes[id_node1]]) {
 597     S1.insert(cells[i]);
 598   }
 599
 600   QSet<vtkIdType> S2;
 601   foreach( int i, n2c[_nodes[id_node2]] ) {
 602     S2.insert(cells[i]);
 603   }
 604
 605   S2.intersect(S1);
 606   EdgeCells = set2Vector(S2, false);
 607   return EdgeCells.size();
 608 }
 609
 610 int SurfaceOperation::getEdgeCells( vtkIdType id_node1, vtkIdType id_node2, QSet <vtkIdType> &EdgeCells )
 611 {
 612   g2l_t _nodes = getPartLocalNodes();
 613   l2g_t cells  = getPartCells();
 614   l2l_t n2c    = getPartN2C();
 615
 616   QSet<vtkIdType> S1;
 617   foreach( int i, n2c[_nodes[id_node1]] ) {
 618     S1.insert( cells[i] );
 619   }
 620
 621   QSet<vtkIdType> S2;
 622   foreach( int i, n2c[_nodes[id_node2]] ) {
 623     S2.insert( cells[i] );
 624   }
 625
 626   EdgeCells = S2.intersect( S1 );
 627   return EdgeCells.size();
 628 }
 629
 630 char SurfaceOperation::getEdgeType(vtkIdType id_node1, vtkIdType id_node2, bool fix_unselected)
 631 {
 632   // compute number of cells around edge [a_node,p2] and put them into neighbour_cells
 633   QVector <vtkIdType> neighbour_cells;
 634   int num_neigh = getEdgeCells(id_node1, id_node2, neighbour_cells) - 1;
 635
 636   // set default value
 637   char edge = EG_SIMPLE_VERTEX;
 638
 639   if (num_neigh == 0) {
 640     edge = EG_BOUNDARY_EDGE_VERTEX;
 641   } else if (num_neigh >= 2) {
 642     edge = EG_BOUNDARY_EDGE_VERTEX;
 643   } else if (num_neigh == 1) {
 644
 645     if (m_Part.isFeatureEdge(id_node1, id_node2, m_FeatureAngle)) {
 646       edge = EG_FEATURE_EDGE_VERTEX;
 647     }
 648
 649     // check the boundary codes
 650     EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
 651     int cell_code_0 = cell_code->GetValue( neighbour_cells[0] );
 652     int cell_code_1 = cell_code->GetValue( neighbour_cells[1] );
 653     if (cell_code_0 !=  cell_code_1) {
 654       edge = EG_BOUNDARY_EDGE_VERTEX;
 655     }
 656
 657     if (fix_unselected) {
 658       if (!m_BoundaryCodes.contains(cell_code_0) || !m_BoundaryCodes.contains(cell_code_1)) {
 659         // does not make sense, but should make the points of the edge fixed
 660         edge = EG_FIXED_VERTEX;
 661       }
 662     }
 663   }
 664
 665   return edge;
 666 }
 667
 668 VertexMeshDensity SurfaceOperation::getVMD( vtkIdType id_node )
 669 {
 670   g2l_t _nodes = getPartLocalNodes();
 671   l2g_t  cells = getPartCells();
 672   l2l_t  n2c   = getPartN2C();
 673
 674   EG_VTKDCN( vtkCharArray, node_type, m_Grid, "node_type" );
 675   EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
 676
 677   VertexMeshDensity VMD;
 678   VMD.type = node_type->GetValue( id_node );
 679   VMD.density = 0;
 680   VMD.CurrentNode = id_node;
 681
 682   foreach( int i_cell, n2c[_nodes[id_node]] ) {
 683     vtkIdType id_cell = cells[i_cell];
 684     VMD.BCmap[cell_code->GetValue( id_cell )] = 2;
 685   }
 686   return( VMD );
 687 }
 688
 689 //////////////////////////////////////////////
 690 double SurfaceOperation::currentVertexAvgDist( vtkIdType id_node )
 691 {
 692   l2g_t  nodes = getPartNodes();
 693   g2l_t _nodes = getPartLocalNodes();
 694   l2l_t  n2n   = getPartN2N();
 695
 696   double total_dist = 0;
 697   double avg_dist = 0;
 698   int N = n2n[_nodes[id_node]].size();
 699   vec3_t C;
 700   m_Grid->GetPoint( id_node, C.data() );
 701   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 702     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 703     vec3_t M;
 704     m_Grid->GetPoint( id_node_neighbour, M.data() );
 705     total_dist += ( M - C ).abs();
 706   }
 707   avg_dist = total_dist / ( double )N;
 708   return( avg_dist );
 709 }
 710
 711 double SurfaceOperation::CurrentMeshDensity( vtkIdType id_node )
 712 {
 713   return 1.0 / currentVertexAvgDist( id_node );
 714 }
 715
 716 ///\todo change meshdensity fields to edgelength fields since this is what is mostly used?
 717
 718 /// desired edge length for id_node
 719 double SurfaceOperation::desiredEdgeLength( vtkIdType id_node )
 720 {
 721   EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
 722   return( 1.0 / characteristic_length_desired->GetValue( id_node ) );
 723 }
 724
 725 /// mean desired edge length for id_cell
 726 double SurfaceOperation::meanDesiredEdgeLength( vtkIdType id_cell )
 727 {
 728   vtkIdType num_pts, *pts;
 729   m_Grid->GetCellPoints( id_cell, num_pts, pts );
 730   int total = 0;
 731   for ( int i = 0; i < num_pts; i++ ) {
 732     total += desiredEdgeLength( pts[i] );
 733   }
 734   return total / ( double )num_pts;
 735 }
 736
 737 QVector <vtkIdType> SurfaceOperation::getPotentialSnapPoints( vtkIdType id_node )
 738 {
 739   if ((id_node < 0) || (id_node >= m_PotentialSnapPoints.size())) {
 740     // UpdatePotentialSnapPoints should probably be called before using this function.
 741     EG_BUG;
 742   }
 743   return m_PotentialSnapPoints[id_node];
 744 }
 745
 746 bool SurfaceOperation::isCell(vtkIdType id_node1, vtkIdType id_node2, vtkIdType id_node3)
 747 {
 748   QVector <vtkIdType> EdgeCells_12;
 749   QVector <vtkIdType> EdgeCells_13;
 750   QVector <vtkIdType> inter;
 751
 752   getEdgeCells( id_node1, id_node2, EdgeCells_12 );
 753   getEdgeCells( id_node1, id_node3, EdgeCells_13 );
 754   qcontIntersection( EdgeCells_12, EdgeCells_13, inter );
 755   if(inter.size()>1) {
 756     qWarning()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
 757     qWarning()<<"EdgeCells_12="<<EdgeCells_12;
 758     qWarning()<<"EdgeCells_13="<<EdgeCells_13;
 759     qWarning()<<"inter="<<inter;
 760     writeGrid(m_Grid, "abort");
 761     EG_BUG;// multiple cells in the same place
 762   }
 763   if(DebugLevel>100) {
 764     qDebug()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
 765     qDebug()<<"EdgeCells_12="<<EdgeCells_12;
 766     qDebug()<<"EdgeCells_13="<<EdgeCells_13;
 767     qDebug()<<"inter="<<inter;
 768   }
 769   return(inter.size()>0);
 770 }
 771
 772 void SurfaceOperation::computeNormals()
 773 {
 774   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 775   m_NodeNormal.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
 776   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 777     if (m_Part.localNode(id_node) != -1) {
 778       QSet<int> bcs;
 779       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 780         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 781         if (isSurface(id_cell, m_Grid)) {
 782           int bc = cell_code->GetValue(id_cell);
 783           if (m_BoundaryCodes.contains(bc)) {
 784             bcs.insert(bc);
 785           }
 786         }
 787       }
 788       int num_bcs = bcs.size();
 789       QVector<vec3_t> normal(num_bcs, vec3_t(0,0,0));
 790       QMap<int,int> bcmap;
 791       int i_bc = 0;
 792       foreach (int bc, bcs) {
 793         bcmap[bc] = i_bc;
 794         ++i_bc;
 795       }
 796       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 797         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 798         if (isSurface(id_cell, m_Grid)) {
 799           int bc = cell_code->GetValue(id_cell);
 800           if (m_BoundaryCodes.contains(bc)) {
 801             vtkIdType N_pts, *pts;
 802             m_Grid->GetCellPoints(id_cell, N_pts, pts);
 803             vec3_t a, b, c;
 804             for (int j = 0; j < N_pts; ++j) {
 805               if (pts[j] == id_node) {
 806                 m_Grid->GetPoint(pts[j], a.data());
 807                 if (j > 0) {
 808                   m_Grid->GetPoint(pts[j-1], b.data());
 809                 } else {
 810                   m_Grid->GetPoint(pts[N_pts-1], b.data());
 811                 }
 812                 if (j < N_pts - 1) {
 813                   m_Grid->GetPoint(pts[j+1], c.data());
 814                 } else {
 815                   m_Grid->GetPoint(pts[0], c.data());
 816                 }
 817               }
 818             }
 819             vec3_t u = b - a;
 820             vec3_t v = c - a;
 821             double alpha = GeometryTools::angle(u, v);
 822             vec3_t n = u.cross(v);
 823             n.normalise();
 824             if (checkVector(n)) {
 825               normal[bcmap[bc]] -= alpha*n;
 826             }
 827           }
 828         }
 829       }
 830       for (int i = 0; i < num_bcs; ++i) {
 831         normal[i].normalise();
 832       }
 833       if (num_bcs > 0) {
 834         if (num_bcs > 1) {
 835           if (num_bcs == 3) {
 836             for (int i = 0; i < num_bcs; ++i) {
 837               for (int j = i + 1; j < num_bcs; ++j) {
 838                 vec3_t n = normal[i] + normal[j];
 839                 n.normalise();
 840                 m_NodeNormal[id_node] += n;
 841               }
 842             }
 843           } else {
 844             for (int i = 0; i < num_bcs; ++i) {
 845               m_NodeNormal[id_node] += normal[i];
 846             }
 847           }
 848         } else {
 849           m_NodeNormal[id_node] = normal[0];
 850         }
 851         m_NodeNormal[id_node].normalise();
 852       }
 853     }
 854   }
 855 }
 856
 857 double SurfaceOperation::normalIrregularity(vtkIdType id_node)
 858 {
 859   double nirr = 0;
 860   QVector<vec3_t> nc(m_Part.n2cGSize(id_node));
 861   for (int i = 0; i < nc.size(); ++i) {
 862     nc[i] = GeometryTools::cellNormal(m_Grid, m_Part.n2cGG(id_node, i));
 863     nc[i].normalise();
 864   }
 865   for (int i = 0; i < nc.size(); ++i) {
 866     for (int j = i + 1; j < nc.size(); ++j) {
 867       nirr += 1.0 - nc[i]*nc[j];
 868     }
 869   }
 870   return nirr;
 871 }
 872
 873 bool SurfaceOperation::isConvexNode(vtkIdType id_node)
 874 {
 875   int N = m_Part.n2nGSize(id_node);
 876   if (N == 0) {
 877     return false;
 878   }
 879   vec3_t x1, x2(0,0,0);
 880   m_Grid->GetPoint(id_node, x1.data());
 881   for (int i = 0; i < N; ++i) {
 882     vec3_t x;
 883     m_Grid->GetPoint(m_Part.n2nGG(id_node, i), x.data());
 884     x2 += x;
 885   }
 886   x2 *= 1.0/N;
 887   if ((x1 - x2)*m_NodeNormal[id_node] > 0) {
 888     return true;
 889   }
 890   return false;
 891 }
 892
 893 double SurfaceOperation::getSurfaceDeviation(vtkIdType id_node)
 894 {
 895   vec3_t x;
 896   m_Grid->GetPoint(id_node, x.data());
 897   double d = 0;
 898   for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 899     vec3_t xe;
 900     m_Grid->GetPoint(m_Part.n2nGG(id_node, i), xe.data());
 901     d += (x - xe)*m_NodeNormal[id_node];
 902   }
 903   if (m_Part.n2nGSize(id_node) > 0) {
 904     d *= 1.0/m_Part.n2nGSize(id_node);
 905   }
 906   return d;
 907 }
 908
 909 bool SurfaceOperation::isFeatureNode(vtkIdType id_node)
 910 {
 911   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 912   char type = node_type->GetValue(id_node);
 913   if (type == EG_FEATURE_CORNER_VERTEX || type == EG_FEATURE_EDGE_VERTEX) {
 914     return true;
 915   }
 916   return false;
 917 }