src/libengrid/surfaceoperation.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2012 enGits GmbH                                     +
   7 // +                                                                      +
   8 // + enGrid is free software: you can redistribute it and/or modify       +
   9 // + it under the terms of the GNU General Public License as published by +
  10 // + the Free Software Foundation, either version 3 of the License, or    +
  11 // + (at your option) any later version.                                  +
  12 // +                                                                      +
  13 // + enGrid is distributed in the hope that it will be useful,            +
  14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  16 // + GNU General Public License for more details.                         +
  17 // +                                                                      +
  18 // + You should have received a copy of the GNU General Public License    +
  19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  20 // +                                                                      +
  21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  22 //
  23 #include "surfaceoperation.h"
  24
  25 #include "guimainwindow.h"
  26
  27 #include <vtkCharArray.h>
  28 #include <vtkMath.h>
  29 #include <vtkCellArray.h>
  30 #include <vtkPolygon.h>
  31
  32 #include "geometrytools.h"
  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", 180, m_FeatureAngle);
  40   m_FeatureAngle = GeometryTools::deg2rad(m_FeatureAngle);
  41   m_EdgeAngle = GeometryTools::deg2rad(m_EdgeAngle);
  42   setEdgeAngle(m_EdgeAngle);
  43   m_BoundarySmoothing = 1;
  44   m_StretchingFactor = 0;
  45   m_UniformSnapPoints = false;
  46 }
  47
  48 void SurfaceOperation::operate()
  49 {
  50
  51 }
  52
  53 ostream& operator<<(ostream &out, stencil_t S)
  54 {
  55   out << "S.id_cell = " << S.id_cell << " ";
  56   out << "S.id_node = " << S.id_node << " ";
  57   out << "S.sameBC = " << S.sameBC << " ";
  58   out << "S.type = " << S.type_cell << " ";
  59   out << "S.p1 = " << S.p1 << " ";
  60   out << "S.p2 = " << S.p2 << " ";
  61   return(out);
  62 }
  63
  64 stencil_t SurfaceOperation::getStencil(vtkIdType id_cell1, int j1)
  65 {
  66   stencil_t S;
  67   {
  68     vtkIdType N_pts, *pts;
  69     m_Grid->GetCellPoints(id_cell1, N_pts, pts);
  70     S.p1 = pts[j1];
  71     S.p2 = pts[0];
  72     if (j1 < N_pts - 1) {
  73       S.p2 = pts[j1 + 1];
  74     }
  75   }
  76   QSet<vtkIdType> cells_p1;
  77   for (int i = 0; i < m_Part.n2cGSize(S.p1); ++i) {
  78     vtkIdType id_cell = m_Part.n2cGG(S.p1, i);
  79     if (id_cell != id_cell1) {
  80       cells_p1.insert(id_cell);
  81     }
  82   }
  83   QSet<vtkIdType> cells_p2;
  84   for (int i = 0; i < m_Part.n2cGSize(S.p2); ++i) {
  85     vtkIdType id_cell = m_Part.n2cGG(S.p2, i);
  86     if (id_cell != id_cell1) {
  87       cells_p2.insert(id_cell);
  88     }
  89   }
  90   QSet<vtkIdType> cells = cells_p1.intersect(cells_p2);
  91   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
  92   S.sameBC = true;
  93   S.id_cell.resize(1);
  94   S.id_cell[0] = id_cell1;
  95   foreach (vtkIdType id_cell, cells) {
  96     if (isSurface(id_cell, m_Grid)) {
  97       S.id_cell.push_back(id_cell);
  98       if (cell_code->GetValue(id_cell) != cell_code->GetValue(id_cell1)) {
  99         S.sameBC = false;
 100       }
 101     }
 102   }
 103   S.id_node.resize(S.id_cell.size());
 104   S.type_cell.resize(S.id_cell.size());
 105   for (int i = 0; i < S.id_cell.size(); ++i) {
 106     vtkIdType N_pts, *pts;
 107     m_Grid->GetCellPoints(S.id_cell[i], N_pts, pts);
 108     S.type_cell[i] = m_Grid->GetCellType(S.id_cell[i]);
 109     for (int j = 0; j < N_pts; ++j) {
 110       if (pts[j] != S.p1 && pts[j] != S.p2) {
 111         S.id_node[i] = pts[j];
 112         break;
 113       }
 114     }
 115   }
 116   return S;
 117 }
 118
 119 int SurfaceOperation::UpdateCurrentMeshDensity()
 120 {
 121   if ( DebugLevel > 0 ) {
 122     cout << "===UpdateMeshDensity START===" << endl;
 123   }
 124   QVector<vtkIdType> cells;
 125   getAllSurfaceCells( cells, m_Grid );
 126   EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
 127   EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
 128   setGrid( m_Grid );
 129   setCells( cells );
 130   if ( DebugLevel > 5 ) {
 131     cout << "cells.size()=" << cells.size() << endl;
 132   }
 133   EG_VTKDCN( vtkDoubleArray, node_meshdensity_current, m_Grid, "node_meshdensity_current" );
 134   l2g_t nodes = getPartNodes();
 135   foreach( vtkIdType node, nodes ) {
 136     node_meshdensity_current->SetValue( node, CurrentMeshDensity( node ) );
 137   }
 138   if ( DebugLevel > 0 ) {
 139     cout << "===UpdateMeshDensity END===" << endl;
 140   }
 141   return( 0 ); ///\todo what for???
 142 }
 143
 144 void SurfaceOperation::readVMD()
 145 {
 146   QString buffer = GuiMainWindow::pointer()->getXmlSection("engrid/surface/table").replace("\n", " ");
 147   int row_count = 0;
 148   int column_count = 0;
 149   m_VMDvector.clear();
 150
 151   if(!buffer.isEmpty()) {
 152     QTextStream in(&buffer, QIODevice::ReadOnly);
 153     in >> row_count >> column_count;
 154     QVector<int> tmp_bcs;
 155     GuiMainWindow::pointer()->getAllBoundaryCodes(tmp_bcs);
 156     if (column_count == tmp_bcs.size() + 3) {
 157       m_VMDvector.fill(VertexMeshDensity(), row_count);
 158       for (int i = 0; i < row_count; ++i) {
 159         int row, column;
 160         QString formula;
 161         foreach (int bc, tmp_bcs) {
 162           in >> row >> column >> formula;
 163           m_VMDvector[row].BCmap[bc] = formula.toInt();
 164         }
 165         in >> row >> column >> formula;
 166         m_VMDvector[row].type = Str2VertexType(formula);
 167         in >> row >> column >> formula;
 168         if (formula == "{{{empty}}}") {
 169           formula = "";
 170         }
 171         m_VMDvector[i].setNodes(formula);
 172         in >> row >> column >> formula;
 173         m_VMDvector[i].density = formula.toDouble();
 174       }
 175     } else {
 176       EG_ERR_RETURN(QObject::tr("Mismatch of number of boundary codes!"));
 177     }
 178   }
 179 }
 180
 181 void SurfaceOperation::updateNodeInfo()
 182 {
 183   setAllCells();
 184   readVMD();
 185   l2g_t nodes = getPartNodes();
 186   computeNormals();
 187   foreach (vtkIdType id_node, nodes) {
 188     EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 189     node_type->SetValue(id_node, getNodeType(id_node, true));
 190
 191     //density index from table
 192     EG_VTKDCN(vtkIntArray, node_specified_density, m_Grid, "node_specified_density");
 193
 194     VertexMeshDensity nodeVMD = getVMD(id_node);
 195     int idx = nodeVMD.findSmallestVMD(m_VMDvector);
 196     node_specified_density->SetValue(id_node, idx);
 197   }
 198   updatePotentialSnapPoints();
 199 }
 200
 201 void SurfaceOperation::updatePotentialSnapPoints()
 202 {
 203   setAllSurfaceCells();
 204   m_PotentialSnapPoints.resize(m_Grid->GetNumberOfPoints());
 205
 206   if (m_UniformSnapPoints) {
 207     m_PotentialSnapPoints.resize(m_Grid->GetNumberOfPoints());
 208     EG_FORALL_NODES(id_node, m_Grid) {
 209       for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 210         m_PotentialSnapPoints[id_node].append(m_Part.n2nGG(id_node, i));
 211       }
 212     }
 213     return;
 214   }
 215
 216   EG_VTKDCN(vtkCharArray, node_type, m_Grid, "node_type");
 217   EG_FORALL_NODES(id_node1, m_Grid) {
 218     m_PotentialSnapPoints[id_node1].clear();
 219     char type1 = node_type->GetValue(id_node1);
 220     for (int i = 0; i < m_Part.n2nGSize(id_node1); ++i) {
 221       vtkIdType id_node2 = m_Part.n2nGG(id_node1, i);
 222       char type2 = node_type->GetValue(id_node2);
 223       if (   (type1 == EG_SIMPLE_VERTEX)
 224           || (type1 == EG_FEATURE_EDGE_VERTEX && (type2 == EG_FEATURE_EDGE_VERTEX || type2 == EG_FEATURE_CORNER_VERTEX))
 225           || (type1 == EG_BOUNDARY_EDGE_VERTEX && (type2 == EG_BOUNDARY_EDGE_VERTEX || type2 == EG_FIXED_VERTEX)))
 226       {
 227         m_PotentialSnapPoints[id_node1].append(id_node2);
 228       }
 229     }
 230   }
 231
 232 }
 233
 234 char SurfaceOperation::getNodeType(vtkIdType id_node, bool fix_unselected)
 235 {
 236   if (m_Part.n2nGSize(id_node) == 0) {
 237     return EG_FIXED_VERTEX;
 238   }
 239
 240   //char type = EG_SIMPLE_VERTEX;
 241   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 242
 243   QSet<int> bcs;
 244   for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 245     vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 246
 247     // fix all vertices that are part of a volume element or a quad
 248     if (m_Grid->GetCellType(id_cell) != VTK_TRIANGLE) {
 249       return EG_FIXED_VERTEX;
 250     }
 251
 252     int bc = cell_code->GetValue(id_cell);
 253
 254     // fix nodes which belong to faces with unselected boundary codes
 255     if (!m_BoundaryCodes.contains(bc) && fix_unselected) {
 256       return EG_FIXED_VERTEX;
 257     }
 258
 259     bcs.insert(bc);
 260   }
 261
 262   if (bcs.size() >= 3 || bcs.size() == 0) {
 263     return EG_FIXED_VERTEX;
 264   }
 265   if (bcs.size() == 2) {
 266     return EG_BOUNDARY_EDGE_VERTEX;
 267   }
 268
 269   SurfaceProjection* proj = GuiMainWindow::pointer()->getSurfProj(*bcs.begin(), true);
 270   if (proj) {
 271
 272     vec3_t x, x0;
 273     m_Grid->GetPoint(id_node, x0.data());
 274
 275     // compute average edge length
 276     double L = 0;
 277     for (int i = 0; i < m_Part.n2nGSize(id_node); ++i) {
 278       vec3_t xn;
 279       m_Grid->GetPoint(m_Part.n2nGG(id_node, i), xn.data());
 280       L += (x0 - xn).abs();
 281     }
 282     L /= m_Part.n2nGSize(id_node);
 283     L *= 0.1;
 284     bool convex = isConvexNode(id_node);
 285
 286     x0 = proj->project(x0, id_node, true, m_NodeNormal[id_node]);
 287     double radius = proj->lastProjRadius();
 288     if (convex) {
 289       x = x0 + L*m_NodeNormal[id_node];
 290     } else {
 291       x = x0 - L*m_NodeNormal[id_node];
 292     }
 293     vec3_t n = proj->lastProjNormal();
 294     if (GeometryTools::angle(n, m_NodeNormal[id_node]) > 0.5*m_FeatureAngle && !proj->lastProjFailed()) {
 295       double d = L/tan(0.5*m_FeatureAngle);
 296       static const int num_steps = 36;
 297       double D_alpha = 2*M_PI/num_steps;
 298       vec3_t v;
 299
 300       v = GeometryTools::orthogonalVector(m_NodeNormal[id_node]);
 301       int num_miss = 0;
 302       for (int i = 0; i < num_steps; ++i) {
 303         v = GeometryTools::rotate(v, m_NodeNormal[id_node], D_alpha);
 304         vec3_t xp = proj->project(x, id_node, true, v, true);
 305         if (proj->lastProjFailed()) {
 306           ++num_miss;
 307         } else {
 308           double l = (x - xp).abs();
 309           if (l >= d) {
 310             ++num_miss;
 311           }
 312         }
 313       }
 314       if (num_miss == 0) {
 315         return EG_FEATURE_CORNER_VERTEX;
 316       }
 317
 318       if (convex) {
 319         x = x0 + L*n;
 320       } else {
 321         x = x0 - L*n;
 322       }
 323       v = GeometryTools::orthogonalVector(n);
 324       int num_hit = 0;
 325       for (int i = 0; i < num_steps; ++i) {
 326         v = GeometryTools::rotate(v, n, D_alpha);
 327         vec3_t xp = proj->project(x, id_node, true, v, true);
 328         if (!proj->lastProjFailed()) {
 329           double l = (x - xp).abs();
 330           if (l < d) {
 331             ++num_hit;
 332           }
 333         }
 334       }
 335       if (num_hit > 0) {
 336         return EG_FEATURE_EDGE_VERTEX;
 337       }
 338     }
 339   }
 340
 341   return EG_SIMPLE_VERTEX;
 342 }
 343
 344 int SurfaceOperation::getEdgeCells(vtkIdType id_node1, vtkIdType id_node2, QVector <vtkIdType> &EdgeCells)
 345 {
 346   g2l_t _nodes = getPartLocalNodes();
 347   l2g_t cells  = getPartCells();
 348   l2l_t n2c    = getPartN2C();
 349
 350   QSet<vtkIdType> S1;
 351   foreach (int i, n2c[_nodes[id_node1]]) {
 352     S1.insert(cells[i]);
 353   }
 354
 355   QSet<vtkIdType> S2;
 356   foreach( int i, n2c[_nodes[id_node2]] ) {
 357     S2.insert(cells[i]);
 358   }
 359
 360   S2.intersect(S1);
 361   EdgeCells = Set2Vector(S2, false);
 362   return EdgeCells.size();
 363 }
 364
 365 int SurfaceOperation::getEdgeCells( vtkIdType id_node1, vtkIdType id_node2, QSet <vtkIdType> &EdgeCells )
 366 {
 367   g2l_t _nodes = getPartLocalNodes();
 368   l2g_t cells  = getPartCells();
 369   l2l_t n2c    = getPartN2C();
 370
 371   QSet<vtkIdType> S1;
 372   foreach( int i, n2c[_nodes[id_node1]] ) {
 373     S1.insert( cells[i] );
 374   }
 375
 376   QSet<vtkIdType> S2;
 377   foreach( int i, n2c[_nodes[id_node2]] ) {
 378     S2.insert( cells[i] );
 379   }
 380
 381   EdgeCells = S2.intersect( S1 );
 382   return EdgeCells.size();
 383 }
 384
 385 char SurfaceOperation::getEdgeType(vtkIdType a_node1, vtkIdType a_node2, bool fix_unselected)
 386 {
 387   double cos_feature_angle = cos(this->m_FeatureAngle);
 388   bool feature_edges_disabled = m_FeatureAngle >= M_PI;
 389
 390   // compute number of cells around edge [a_node,p2] and put them into neighbour_cells
 391   QVector <vtkIdType> neighbour_cells;
 392   int numNei = getEdgeCells(a_node1, a_node2, neighbour_cells) - 1;
 393
 394   // set default value
 395   char edge = EG_SIMPLE_VERTEX;
 396
 397   if (numNei == 0) {
 398     edge = EG_BOUNDARY_EDGE_VERTEX;
 399   } else if (numNei >= 2) {
 400     edge = EG_BOUNDARY_EDGE_VERTEX;
 401   } else if (numNei == 1) {
 402
 403     // check angle between cell1 and cell2 against FeatureAngle
 404     double cosa = cosAngle(m_Grid, neighbour_cells[0], neighbour_cells[1]);
 405     if (cosa <= cos_feature_angle && !feature_edges_disabled) {
 406       edge = EG_FEATURE_EDGE_VERTEX;
 407     }
 408
 409     // check the boundary codes
 410     EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
 411     int cell_code_0 = cell_code->GetValue( neighbour_cells[0] );
 412     int cell_code_1 = cell_code->GetValue( neighbour_cells[1] );
 413     if (cell_code_0 !=  cell_code_1) {
 414       edge = EG_BOUNDARY_EDGE_VERTEX;
 415     }
 416
 417     if (fix_unselected) {
 418       if (!m_BoundaryCodes.contains(cell_code_0) || !m_BoundaryCodes.contains(cell_code_1)) {
 419         // does not make sense, but should make the points of the edge fixed
 420         edge = EG_FIXED_VERTEX;
 421       }
 422     }
 423   }
 424
 425   return edge;
 426 }
 427
 428 VertexMeshDensity SurfaceOperation::getVMD( vtkIdType id_node )
 429 {
 430   g2l_t _nodes = getPartLocalNodes();
 431   l2g_t  cells = getPartCells();
 432   l2l_t  n2c   = getPartN2C();
 433
 434   EG_VTKDCN( vtkCharArray, node_type, m_Grid, "node_type" );
 435   EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
 436
 437   VertexMeshDensity VMD;
 438   VMD.type = node_type->GetValue( id_node );
 439   VMD.density = 0;
 440   VMD.CurrentNode = id_node;
 441
 442   foreach( int i_cell, n2c[_nodes[id_node]] ) {
 443     vtkIdType id_cell = cells[i_cell];
 444     VMD.BCmap[cell_code->GetValue( id_cell )] = 2;
 445   }
 446   return( VMD );
 447 }
 448
 449 //////////////////////////////////////////////
 450 double SurfaceOperation::currentVertexAvgDist( vtkIdType id_node )
 451 {
 452   l2g_t  nodes = getPartNodes();
 453   g2l_t _nodes = getPartLocalNodes();
 454   l2l_t  n2n   = getPartN2N();
 455
 456   double total_dist = 0;
 457   double avg_dist = 0;
 458   int N = n2n[_nodes[id_node]].size();
 459   vec3_t C;
 460   m_Grid->GetPoint( id_node, C.data() );
 461   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 462     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 463     vec3_t M;
 464     m_Grid->GetPoint( id_node_neighbour, M.data() );
 465     total_dist += ( M - C ).abs();
 466   }
 467   avg_dist = total_dist / ( double )N;
 468   return( avg_dist );
 469 }
 470
 471 double SurfaceOperation::CurrentMeshDensity( vtkIdType id_node )
 472 {
 473   return 1.0 / currentVertexAvgDist( id_node );
 474 }
 475
 476 ///\todo change meshdensity fields to edgelength fields since this is what is mostly used?
 477
 478 /// desired edge length for id_node
 479 double SurfaceOperation::desiredEdgeLength( vtkIdType id_node )
 480 {
 481   EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
 482   return( 1.0 / characteristic_length_desired->GetValue( id_node ) );
 483 }
 484
 485 /// mean desired edge length for id_cell
 486 double SurfaceOperation::meanDesiredEdgeLength( vtkIdType id_cell )
 487 {
 488   vtkIdType num_pts, *pts;
 489   m_Grid->GetCellPoints( id_cell, num_pts, pts );
 490   int total = 0;
 491   for ( int i = 0; i < num_pts; i++ ) {
 492     total += desiredEdgeLength( pts[i] );
 493   }
 494   return total / ( double )num_pts;
 495 }
 496
 497 QVector <vtkIdType> SurfaceOperation::getPotentialSnapPoints( vtkIdType id_node )
 498 {
 499   if ((id_node < 0) || (id_node >= m_PotentialSnapPoints.size())) {
 500     // UpdatePotentialSnapPoints should probably be called before using this function.
 501     EG_BUG;
 502   }
 503   return m_PotentialSnapPoints[id_node];
 504 }
 505
 506 bool SurfaceOperation::isCell(vtkIdType id_node1, vtkIdType id_node2, vtkIdType id_node3)
 507 {
 508   QVector <vtkIdType> EdgeCells_12;
 509   QVector <vtkIdType> EdgeCells_13;
 510   QVector <vtkIdType> inter;
 511
 512   getEdgeCells( id_node1, id_node2, EdgeCells_12 );
 513   getEdgeCells( id_node1, id_node3, EdgeCells_13 );
 514   qcontIntersection( EdgeCells_12, EdgeCells_13, inter );
 515   if(inter.size()>1) {
 516     qWarning()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
 517     qWarning()<<"EdgeCells_12="<<EdgeCells_12;
 518     qWarning()<<"EdgeCells_13="<<EdgeCells_13;
 519     qWarning()<<"inter="<<inter;
 520     writeGrid(m_Grid, "abort");
 521     EG_BUG;// multiple cells in the same place
 522   }
 523   if(DebugLevel>100) {
 524     qDebug()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
 525     qDebug()<<"EdgeCells_12="<<EdgeCells_12;
 526     qDebug()<<"EdgeCells_13="<<EdgeCells_13;
 527     qDebug()<<"inter="<<inter;
 528   }
 529   return(inter.size()>0);
 530 }
 531
 532 void SurfaceOperation::computeNormals()
 533 {
 534   EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
 535   m_NodeNormal.fill(vec3_t(0,0,0), m_Grid->GetNumberOfPoints());
 536   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
 537     if (m_Part.localNode(id_node) != -1) {
 538       QSet<int> bcs;
 539       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 540         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 541         if (isSurface(id_cell, m_Grid)) {
 542           int bc = cell_code->GetValue(id_cell);
 543           if (m_BoundaryCodes.contains(bc)) {
 544             bcs.insert(bc);
 545           }
 546         }
 547       }
 548       int num_bcs = bcs.size();
 549       QVector<vec3_t> normal(num_bcs, vec3_t(0,0,0));
 550       QMap<int,int> bcmap;
 551       int i_bc = 0;
 552       foreach (int bc, bcs) {
 553         bcmap[bc] = i_bc;
 554         ++i_bc;
 555       }
 556       for (int i = 0; i < m_Part.n2cGSize(id_node); ++i) {
 557         vtkIdType id_cell = m_Part.n2cGG(id_node, i);
 558         if (isSurface(id_cell, m_Grid)) {
 559           int bc = cell_code->GetValue(id_cell);
 560           if (m_BoundaryCodes.contains(bc)) {
 561             vtkIdType N_pts, *pts;
 562             m_Grid->GetCellPoints(id_cell, N_pts, pts);
 563             vec3_t a, b, c;
 564             for (int j = 0; j < N_pts; ++j) {
 565               if (pts[j] == id_node) {
 566                 m_Grid->GetPoint(pts[j], a.data());
 567                 if (j > 0) {
 568                   m_Grid->GetPoint(pts[j-1], b.data());
 569                 } else {
 570                   m_Grid->GetPoint(pts[N_pts-1], b.data());
 571                 }
 572                 if (j < N_pts - 1) {
 573                   m_Grid->GetPoint(pts[j+1], c.data());
 574                 } else {
 575                   m_Grid->GetPoint(pts[0], c.data());
 576                 }
 577               }
 578             }
 579             vec3_t u = b - a;
 580             vec3_t v = c - a;
 581             double alpha = GeometryTools::angle(u, v);
 582             vec3_t n = u.cross(v);
 583             n.normalise();
 584             if (checkVector(n)) {
 585               normal[bcmap[bc]] -= alpha*n;
 586             }
 587           }
 588         }
 589       }
 590       for (int i = 0; i < num_bcs; ++i) {
 591         normal[i].normalise();
 592       }
 593       if (num_bcs > 0) {
 594         if (num_bcs > 1) {
 595           if (num_bcs == 3) {
 596             for (int i = 0; i < num_bcs; ++i) {
 597               for (int j = i + 1; j < num_bcs; ++j) {
 598                 vec3_t n = normal[i] + normal[j];
 599                 n.normalise();
 600                 m_NodeNormal[id_node] += n;
 601               }
 602             }
 603           } else {
 604             for (int i = 0; i < num_bcs; ++i) {
 605               m_NodeNormal[id_node] += normal[i];
 606             }
 607           }
 608         } else {
 609           m_NodeNormal[id_node] = normal[0];
 610         }
 611         m_NodeNormal[id_node].normalise();
 612       }
 613     }
 614   }
 615 }
 616
 617 double SurfaceOperation::normalIrregularity(vtkIdType id_node)
 618 {
 619   double nirr = 0;
 620   QVector<vec3_t> nc(m_Part.n2cGSize(id_node));
 621   for (int i = 0; i < nc.size(); ++i) {
 622     nc[i] = GeometryTools::cellNormal(m_Grid, m_Part.n2cGG(id_node, i));
 623     nc[i].normalise();
 624   }
 625   for (int i = 0; i < nc.size(); ++i) {
 626     for (int j = i + 1; j < nc.size(); ++j) {
 627       nirr += 1.0 - nc[i]*nc[j];
 628     }
 629   }
 630   return nirr;
 631 }
 632
 633 bool SurfaceOperation::isConvexNode(vtkIdType id_node)
 634 {
 635   int N = m_Part.n2nGSize(id_node);
 636   if (N == 0) {
 637     return false;
 638   }
 639   vec3_t x1, x2(0,0,0);
 640   m_Grid->GetPoint(id_node, x1.data());
 641   for (int i = 0; i < N; ++i) {
 642     vec3_t x;
 643     m_Grid->GetPoint(m_Part.n2nGG(id_node, i), x.data());
 644     x2 += x;
 645   }
 646   x2 *= 1.0/N;
 647   if ((x1 - x2)*m_NodeNormal[id_node] < 0) {
 648     return true;
 649   }
 650   return false;
 651 }
 652
 653