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-08-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           if (type1 == EG_FEATURE_EDGE_VERTEX && isFeatureNode(id_node2)) {
 419             QVector<vtkIdType> edge_faces;
 420             m_Part.getEdgeFaces(id_node1, id_node2, edge_faces);
 421             if (edge_faces.size() == 2) {
 422               vec3_t n1 = cellNormal(m_Grid, edge_faces[0]);
 423               vec3_t n2 = cellNormal(m_Grid, edge_faces[1]);
 424               if (GeometryTools::angle(n1, n2) >= m_FeatureAngle) {
 425                 m_PotentialSnapPoints[id_node1].append(id_node2);
 426               }
 427             }
 428           }
 429         } else {
 430           if (   (type1 == EG_SIMPLE_VERTEX)
 431                  || (type1 == EG_FEATURE_EDGE_VERTEX)
 432                  || (type1 == EG_BOUNDARY_EDGE_VERTEX && (type2 == EG_BOUNDARY_EDGE_VERTEX || type2 == EG_FIXED_VERTEX)))
 433           {
 434             if (checkSnapPointPairForBcMatch(id_node1, id_node2)) {
 435               m_PotentialSnapPoints[id_node1].append(id_node2);
 436             }
 437           }
 438         }
 439       }
 440
 441       // make sure feature edge vertices have at least two snap points ...
 442       /*
 443     if (type1 == EG_FEATURE_EDGE_VERTEX) {
 444       if (m_PotentialSnapPoints[id_node1].size() < 2) {
 445         m_PotentialSnapPoints[id_node1].clear();
 446       }
 447     }
 448     */
 449     }
 450   }
 451 }
 452
 453 double SurfaceOperation::edgeAngle(vtkIdType id_node1, vtkIdType id_node2)
 454 {
 455   QVector<vtkIdType> faces;
 456   m_Part.getEdgeFaces(id_node1, id_node2, faces);
 457   if (faces.size() == 2) {
 458     vec3_t n1 = cellNormal(m_Grid, faces[0]);
 459     vec3_t n2 = cellNormal(m_Grid, faces[1]);
 460     return GeometryTools::angle(n1, n2);
 461   }
 462   return 0;
 463 }
 464
 465 /*
 466 char SurfaceOperation::getNodeType(vtkIdType id_node, bool fix_unselected)
 467 {
 468   char type = EG_SIMPLE_VERTEX;
 469
 470   if (m_Part.n2nGSize(id_node) == 0) {
 471     return EG_FIXED_VERTEX;
 472   }
 473
 474   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 475
 476   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 477     vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 478
 479     // fix all vertices that are part of a volume element, quad, or polygon
 480     if (m_Grid->GetCellType(id_cell) != VTK_TRIANGLE) {
 481       type = max(type, char(EG_FIXED_VERTEX));
 482     }
 483
 484     int bc = cell_code->GetValue(id_cell);
 485
 486     // fix nodes which belong to faces with unselected boundary codes
 487     if (!m_BoundaryCodes.contains(bc) && fix_unselected) {
 488       type = max(type, char(EG_FIXED_VERTEX));
 489     }
 490
 491   }
 492
 493   int num_bcs = m_Part.n2bcGSize(id_node);
 494   if (num_bcs >= 3 || num_bcs == 0) {
 495     type = max(type, char(EG_FIXED_VERTEX));
 496   }
 497   if (num_bcs == 2) {
 498     QList<vtkIdType> neigh;
 499     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 500       vtkIdType id_neigh = m_Part.n2nGG(id_node, i);
 501       if (getEdgeType(id_node, id_neigh) == EG_BOUNDARY_EDGE_VERTEX) {
 502         neigh << id_neigh;
 503       }
 504     }
 505     if (neigh.size() == 2) {
 506       if (M_PI - GeometryTools::angle(m_Grid, neigh[0], id_node, neigh[1]) >= m_EdgeAngle) {
 507         type = max(type, char(EG_FIXED_VERTEX));
 508       }
 509     }
 510     type = max(type, char(EG_BOUNDARY_EDGE_VERTEX));
 511   }
 512
 513   // features are defined via angles
 514   QList<vtkIdType> feature_neigh;
 515   for (int i_neigh = 0; i_neigh < m_Part.n2nGSize(id_node); ++i_neigh) {
 516     vtkIdType id_neigh = m_Part.n2nGG(id_node, i_neigh);
 517     QVector<vtkIdType> id_face;
 518     m_Part.getEdgeFaces(id_node, id_neigh, id_face);
 519     if (id_face.size() == 2) {
 520       vec3_t n1    = cellNormal(m_Grid, id_face[0]);
 521       vec3_t n2    = cellNormal(m_Grid, id_face[1]);
 522       double alpha = GeometryTools::angle(n1, n2);
 523       if (alpha >= m_FeatureAngle) {
 524         feature_neigh.append(id_neigh);
 525       }
 526     }
 527   }
 528   if (feature_neigh.size() == 0) {
 529     type = max(type, char(EG_SIMPLE_VERTEX));
 530   } else if (feature_neigh.size() == 2) {
 531
 532     vec3_t x1, x2, xc;
 533     m_Grid->GetPoint(id_node, xc.data());
 534     m_Grid->GetPoint(feature_neigh[0], x1.data());
 535     m_Grid->GetPoint(feature_neigh[1], x2.data());
 536     double alpha1 = edgeAngle(id_node, feature_neigh[0]);
 537     double alpha2 = edgeAngle(id_node, feature_neigh[1]);
 538     if (alpha1 >= m_FeatureAngle && alpha2 >= m_FeatureAngle) {
 539       if (GeometryTools::angle(xc - x1, x2 - xc) > m_EdgeAngle) {
 540         type = max(type, char(EG_FIXED_VERTEX));
 541       }
 542     }
 543     type = max(type, char(EG_FEATURE_EDGE_VERTEX));
 544   } else {
 545     type = max(type, char(EG_FIXED_VERTEX));
 546     //return EG_FIXED_VERTEX; //EG_FEATURE_CORNER_VERTEX;
 547   }
 548
 549   return type;
 550 }
 551 */
 552
 553 char SurfaceOperation::getNodeType(vtkIdType id_node, bool fix_unselected)
 554 {
 555   l2g_t  nodes = getPartNodes();
 556   g2l_t _nodes = getPartLocalNodes();
 557   l2l_t  n2n   = getPartN2N();
 558
 559   // fix all vertices that are part of a volume element or a quad
 560   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 561     if (m_Grid->GetCellType(m_Part.n2cGG(id_node, i)) != VTK_TRIANGLE) {
 562       return EG_FIXED_VERTEX;
 563     }
 564   }
 565
 566   //initialize default value
 567   char type = EG_SIMPLE_VERTEX;
 568
 569   //loop through edges around id_node
 570
 571   QVector <vtkIdType> edges;
 572
 573   double CosEdgeAngle = cos(this->m_EdgeAngle);
 574
 575   foreach( int i_node2, n2n[_nodes[id_node]] ) {
 576     vtkIdType id_node2 = nodes[i_node2];
 577     //-----------------------
 578     //determine edge type
 579     char edge = getEdgeType(id_node2, id_node, fix_unselected);
 580
 581     //-----------------------
 582     //determine node type pre-processing (count nb of complex edges if the node is complex, otherwise, just count the nb of edges)
 583     if ( edge && type == EG_SIMPLE_VERTEX ) {
 584       edges.clear();
 585       edges.push_back( id_node2 );
 586       type = edge;
 587     }
 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     } //if along edge
 616   } //if edge vertex
 617
 618   return type;
 619 }
 620
 621
 622
 623 int SurfaceOperation::getEdgeCells(vtkIdType id_node1, vtkIdType id_node2, QVector <vtkIdType> &EdgeCells)
 624 {
 625   g2l_t _nodes = getPartLocalNodes();
 626   l2g_t cells  = getPartCells();
 627   l2l_t n2c    = getPartN2C();
 628
 629   QSet<vtkIdType> S1;
 630   foreach (int i, n2c[_nodes[id_node1]]) {
 631     S1.insert(cells[i]);
 632   }
 633
 634   QSet<vtkIdType> S2;
 635   foreach( int i, n2c[_nodes[id_node2]] ) {
 636     S2.insert(cells[i]);
 637   }
 638
 639   S2.intersect(S1);
 640   EdgeCells = set2Vector(S2, false);
 641   return EdgeCells.size();
 642 }
 643
 644 int SurfaceOperation::getEdgeCells( vtkIdType id_node1, vtkIdType id_node2, QSet <vtkIdType> &EdgeCells )
 645 {
 646   g2l_t _nodes = getPartLocalNodes();
 647   l2g_t cells  = getPartCells();
 648   l2l_t n2c    = getPartN2C();
 649
 650   QSet<vtkIdType> S1;
 651   foreach( int i, n2c[_nodes[id_node1]] ) {
 652     S1.insert( cells[i] );
 653   }
 654
 655   QSet<vtkIdType> S2;
 656   foreach( int i, n2c[_nodes[id_node2]] ) {
 657     S2.insert( cells[i] );
 658   }
 659
 660   EdgeCells = S2.intersect( S1 );
 661   return EdgeCells.size();
 662 }
 663
 664 char SurfaceOperation::getEdgeType(vtkIdType id_node1, vtkIdType id_node2, bool fix_unselected)
 665 {
 666   // compute number of cells around edge [a_node,p2] and put them into neighbour_cells
 667   QVector <vtkIdType> neighbour_cells;
 668   int num_neigh = getEdgeCells(id_node1, id_node2, neighbour_cells) - 1;
 669
 670   // set default value
 671   char edge = EG_SIMPLE_VERTEX;
 672
 673   if (num_neigh == 0) {
 674     edge = EG_BOUNDARY_EDGE_VERTEX;
 675   } else if (num_neigh >= 2) {
 676     edge = EG_BOUNDARY_EDGE_VERTEX;
 677   } else if (num_neigh == 1) {
 678
 679     /*
 680     if (m_Part.isFeatureEdge(id_node1, id_node2, m_FeatureAngle)) {
 681       edge = EG_FEATURE_EDGE_VERTEX;
 682     }
 683     */
 684     vec3_t n1 = cellNormal(m_Grid, neighbour_cells[0]);
 685     vec3_t n2 = cellNormal(m_Grid, neighbour_cells[1]);
 686     if (GeometryTools::angle(n1, n2) > m_FeatureAngle) {
 687       edge = EG_FEATURE_EDGE_VERTEX;
 688     }
 689
 690     // check the boundary codes
 691     EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
 692     int cell_code_0 = cell_code->GetValue( neighbour_cells[0] );
 693     int cell_code_1 = cell_code->GetValue( neighbour_cells[1] );
 694     if (cell_code_0 !=  cell_code_1) {
 695       edge = EG_BOUNDARY_EDGE_VERTEX;
 696     }
 697
 698     if (fix_unselected) {
 699       if (!m_BoundaryCodes.contains(cell_code_0) || !m_BoundaryCodes.contains(cell_code_1)) {
 700         // does not make sense, but should make the points of the edge fixed
 701         edge = EG_FIXED_VERTEX;
 702       }
 703     }
 704   }
 705
 706   return edge;
 707 }
 708
 709 VertexMeshDensity SurfaceOperation::getVMD( vtkIdType id_node )
 710 {
 711   g2l_t _nodes = getPartLocalNodes();
 712   l2g_t  cells = getPartCells();
 713   l2l_t  n2c   = getPartN2C();
 714
 715   EG_VTKDCN( vtkCharArray, node_type, m_Grid, "node_type" );
 716   EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
 717
 718   VertexMeshDensity VMD;
 719   VMD.type = node_type->GetValue( id_node );
 720   VMD.density = 0;
 721   VMD.CurrentNode = id_node;
 722
 723   foreach( int i_cell, n2c[_nodes[id_node]] ) {
 724     vtkIdType id_cell = cells[i_cell];
 725     VMD.BCmap[cell_code->GetValue( id_cell )] = 2;
 726   }
 727   return( VMD );
 728 }
 729
 730 //////////////////////////////////////////////
 731 double SurfaceOperation::currentVertexAvgDist( vtkIdType id_node )
 732 {
 733   l2g_t  nodes = getPartNodes();
 734   g2l_t _nodes = getPartLocalNodes();
 735   l2l_t  n2n   = getPartN2N();
 736
 737   double total_dist = 0;
 738   double avg_dist = 0;
 739   int N = n2n[_nodes[id_node]].size();
 740   vec3_t C;
 741   m_Grid->GetPoint( id_node, C.data() );
 742   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 743     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 744     vec3_t M;
 745     m_Grid->GetPoint( id_node_neighbour, M.data() );
 746     total_dist += ( M - C ).abs();
 747   }
 748   avg_dist = total_dist / ( double )N;
 749   return( avg_dist );
 750 }
 751
 752 double SurfaceOperation::CurrentMeshDensity( vtkIdType id_node )
 753 {
 754   return 1.0 / currentVertexAvgDist( id_node );
 755 }
 756
 757 ///\todo change meshdensity fields to edgelength fields since this is what is mostly used?
 758
 759 /// desired edge length for id_node
 760 double SurfaceOperation::desiredEdgeLength( vtkIdType id_node )
 761 {
 762   EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
 763   return( 1.0 / characteristic_length_desired->GetValue( id_node ) );
 764 }
 765
 766 /// mean desired edge length for id_cell
 767 double SurfaceOperation::meanDesiredEdgeLength( vtkIdType id_cell )
 768 {
 769   vtkIdType num_pts, *pts;
 770   m_Grid->GetCellPoints( id_cell, num_pts, pts );
 771   int total = 0;
 772   for ( int i = 0; i < num_pts; i++ ) {
 773     total += desiredEdgeLength( pts[i] );
 774   }
 775   return total / ( double )num_pts;
 776 }
 777
 778 QVector <vtkIdType> SurfaceOperation::getPotentialSnapPoints( vtkIdType id_node )
 779 {
 780   if ((id_node < 0) || (id_node >= m_PotentialSnapPoints.size())) {
 781     // UpdatePotentialSnapPoints should probably be called before using this function.
 782     EG_BUG;
 783   }
 784   return m_PotentialSnapPoints[id_node];
 785 }
 786
 787 bool SurfaceOperation::isCell(vtkIdType id_node1, vtkIdType id_node2, vtkIdType id_node3)
 788 {
 789   QVector <vtkIdType> EdgeCells_12;
 790   QVector <vtkIdType> EdgeCells_13;
 791   QVector <vtkIdType> inter;
 792
 793   getEdgeCells( id_node1, id_node2, EdgeCells_12 );
 794   getEdgeCells( id_node1, id_node3, EdgeCells_13 );
 795   qcontIntersection( EdgeCells_12, EdgeCells_13, inter );
 796   if(inter.size()>1) {
 797     qWarning()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
 798     qWarning()<<"EdgeCells_12="<<EdgeCells_12;
 799     qWarning()<<"EdgeCells_13="<<EdgeCells_13;
 800     qWarning()<<"inter="<<inter;
 801     writeGrid(m_Grid, "abort");
 802     EG_BUG;// multiple cells in the same place
 803   }
 804   if(DebugLevel>100) {
 805     qDebug()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
 806     qDebug()<<"EdgeCells_12="<<EdgeCells_12;
 807     qDebug()<<"EdgeCells_13="<<EdgeCells_13;
 808     qDebug()<<"inter="<<inter;
 809   }
 810   return(inter.size()>0);
 811 }
 812
 813 void SurfaceOperation::computeNormals()
 814 {
 815   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 816   m_NodeNormal.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
 817   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 818     if (m_Part.localNode(id_node) != -1) {
 819       QSet<int> bcs;
 820       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 821         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 822         if (isSurface(id_cell, m_Grid)) {
 823           int bc = cell_code->GetValue(id_cell);
 824           if (m_BoundaryCodes.contains(bc)) {
 825             bcs.insert(bc);
 826           }
 827         }
 828       }
 829       int num_bcs = bcs.size();
 830       QVector<vec3_t> normal(num_bcs, vec3_t(0,0,0));
 831       QMap<int,int> bcmap;
 832       int i_bc = 0;
 833       foreach (int bc, bcs) {
 834         bcmap[bc] = i_bc;
 835         ++i_bc;
 836       }
 837       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 838         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 839         if (isSurface(id_cell, m_Grid)) {
 840           int bc = cell_code->GetValue(id_cell);
 841           if (m_BoundaryCodes.contains(bc)) {
 842             vtkIdType N_pts, *pts;
 843             m_Grid->GetCellPoints(id_cell, N_pts, pts);
 844             vec3_t a, b, c;
 845             for (int j = 0; j < N_pts; ++j) {
 846               if (pts[j] == id_node) {
 847                 m_Grid->GetPoint(pts[j], a.data());
 848                 if (j > 0) {
 849                   m_Grid->GetPoint(pts[j-1], b.data());
 850                 } else {
 851                   m_Grid->GetPoint(pts[N_pts-1], b.data());
 852                 }
 853                 if (j < N_pts - 1) {
 854                   m_Grid->GetPoint(pts[j+1], c.data());
 855                 } else {
 856                   m_Grid->GetPoint(pts[0], c.data());
 857                 }
 858               }
 859             }
 860             vec3_t u = b - a;
 861             vec3_t v = c - a;
 862             double alpha = GeometryTools::angle(u, v);
 863             vec3_t n = u.cross(v);
 864             n.normalise();
 865             if (checkVector(n)) {
 866               normal[bcmap[bc]] -= alpha*n;
 867             }
 868           }
 869         }
 870       }
 871       for (int i = 0; i < num_bcs; ++i) {
 872         normal[i].normalise();
 873       }
 874       if (num_bcs > 0) {
 875         if (num_bcs > 1) {
 876           if (num_bcs == 3) {
 877             for (int i = 0; i < num_bcs; ++i) {
 878               for (int j = i + 1; j < num_bcs; ++j) {
 879                 vec3_t n = normal[i] + normal[j];
 880                 n.normalise();
 881                 m_NodeNormal[id_node] += n;
 882               }
 883             }
 884           } else {
 885             for (int i = 0; i < num_bcs; ++i) {
 886               m_NodeNormal[id_node] += normal[i];
 887             }
 888           }
 889         } else {
 890           m_NodeNormal[id_node] = normal[0];
 891         }
 892         m_NodeNormal[id_node].normalise();
 893       }
 894     }
 895   }
 896 }
 897
 898 double SurfaceOperation::normalIrregularity(vtkIdType id_node)
 899 {
 900   double nirr = 0;
 901   QVector<vec3_t> nc(m_Part.n2cGSize(id_node));
 902   for (int i = 0; i < nc.size(); ++i) {
 903     nc[i] = GeometryTools::cellNormal(m_Grid, m_Part.n2cGG(id_node, i));
 904     nc[i].normalise();
 905   }
 906   for (int i = 0; i < nc.size(); ++i) {
 907     for (int j = i + 1; j < nc.size(); ++j) {
 908       nirr += 1.0 - nc[i]*nc[j];
 909     }
 910   }
 911   return nirr;
 912 }
 913
 914 bool SurfaceOperation::isConvexNode(vtkIdType id_node)
 915 {
 916   int N = m_Part.n2nGSize(id_node);
 917   if (N == 0) {
 918     return false;
 919   }
 920   vec3_t x1, x2(0,0,0);
 921   m_Grid->GetPoint(id_node, x1.data());
 922   for (int i = 0; i < N; ++i) {
 923     vec3_t x;
 924     m_Grid->GetPoint(m_Part.n2nGG(id_node, i), x.data());
 925     x2 += x;
 926   }
 927   x2 *= 1.0/N;
 928   if ((x1 - x2)*m_NodeNormal[id_node] > 0) {
 929     return true;
 930   }
 931   return false;
 932 }
 933
 934 double SurfaceOperation::getSurfaceDeviation(vtkIdType id_node)
 935 {
 936   vec3_t x;
 937   m_Grid->GetPoint(id_node, x.data());
 938   double d = 0;
 939   for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 940     vec3_t xe;
 941     m_Grid->GetPoint(m_Part.n2nGG(id_node, i), xe.data());
 942     d += (x - xe)*m_NodeNormal[id_node];
 943   }
 944   if (m_Part.n2nGSize(id_node) > 0) {
 945     d *= 1.0/m_Part.n2nGSize(id_node);
 946   }
 947   return d;
 948 }
 949
 950 bool SurfaceOperation::isFeatureNode(vtkIdType id_node)
 951 {
 952   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 953   char type = node_type->GetValue(id_node);
 954   if (type == EG_FEATURE_CORNER_VERTEX || type == EG_FEATURE_EDGE_VERTEX) {
 955     return true;
 956   }
 957   return false;
 958 }