src/surfaceoperation.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008,2009 Oliver Gloth                                     +
   7 // +                                                                      +
   8 // + enGrid is free software: you can redistribute it and/or modify       +
   9 // + it under the terms of the GNU General Public License as published by +
  10 // + the Free Software Foundation, either version 3 of the License, or    +
  11 // + (at your option) any later version.                                  +
  12 // +                                                                      +
  13 // + enGrid is distributed in the hope that it will be useful,            +
  14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  16 // + GNU General Public License for more details.                         +
  17 // +                                                                      +
  18 // + You should have received a copy of the GNU General Public License    +
  19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  20 // +                                                                      +
  21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  22 //
  23 #include "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()
  36     : Operation()
  37 {
  38   //default values for determining node types and for smoothing operations
  39   Convergence = 0;
  40   NumberOfIterations = 20;
  41   RelaxationFactor = 0.01;
  42   m_AllowFeatureEdgeVertices = 1;//0 by default in VTK, but we need 1 to avoid the "potatoe effect" ^^
  43   FeatureAngle = 45;
  44   EdgeAngle = 15;
  45   BoundarySmoothing = 1;
  46 }
  47
  48 void SurfaceOperation::operate()
  49 {
  50
  51 }
  52
  53 ostream& operator<<( ostream &out, stencil_t S )
  54 {
  55   out << "S.id_cell1=" << S.id_cell1 << " ";
  56   out << "S.id_cell2=" << S.id_cell2 << " ";
  57   out << "S.sameBC=" << S.sameBC << " ";
  58   out << "S.twocells=" << S.twocells << " ";
  59   out << "S.neighbour_type=" << S.neighbour_type << " ";
  60   out << "[";
  61   for ( int i = 0; i < 4; i++ ) {
  62     out << S.p[i];
  63     if ( i != 3 ) out << ",";
  64   }
  65   out << "]";
  66   return( out );
  67 }
  68
  69 stencil_t SurfaceOperation::getStencil( vtkIdType id_cell1, int j1 )
  70 {
  71   l2g_t  cells = getPartCells();
  72   g2l_t _cells = getPartLocalCells();
  73   l2l_t  c2c   = getPartC2C();
  74
  75   if ( grid->GetCellType( id_cell1 ) != VTK_TRIANGLE ) {
  76     cout << "CELL IS NOT A TRIANGLE" << endl;
  77     EG_BUG;
  78   }
  79
  80   //return variable
  81   stencil_t S;
  82
  83   //default values:
  84   S.sameBC = false;
  85   S.twocells = false;
  86   S.neighbour_type = -1;
  87
  88   //initialize first cell
  89   S.id_cell1 = id_cell1;
  90   vtkIdType N1, *pts1;
  91   grid->GetCellPoints( S.id_cell1, N1, pts1 );
  92   //place points 0,1,3
  93   if ( j1 == 0 ) { S.p[0] = pts1[2]; S.p[1] = pts1[0]; S.p[3] = pts1[1]; }
  94   else if ( j1 == 1 ) { S.p[0] = pts1[0]; S.p[1] = pts1[1]; S.p[3] = pts1[2]; }
  95   else if ( j1 == 2 ) { S.p[0] = pts1[1]; S.p[1] = pts1[2]; S.p[3] = pts1[0]; };
  96
  97   //initialize second cell
  98   S.id_cell2 = -1;
  99   S.p[2] = -1;
 100
 101   //twocells
 102   if ( c2c[_cells[id_cell1]][j1] != -1 ) { //if neighbour cell
 103
 104     //twocells
 105     S.twocells = true;
 106     S.id_cell2 = cells[c2c[_cells[id_cell1]][j1]];
 107
 108     //sameBC
 109     EG_VTKDCC( vtkIntArray, cell_code, grid, "cell_code" );
 110     if ( cell_code->GetValue( S.id_cell1 ) == cell_code->GetValue( S.id_cell2 ) ) S.sameBC = true;
 111
 112     //neighbour_type
 113     S.neighbour_type = grid->GetCellType( S.id_cell2 );
 114     if ( S.neighbour_type == VTK_TRIANGLE ) {//if neighbour cell is a triangle
 115       vtkIdType N2, *pts2;
 116       grid->GetCellPoints( S.id_cell2, N2, pts2 );
 117
 118       //place point 2
 119       bool p2 = false;
 120       if ( c2c[_cells[S.id_cell2]][0] != -1 ) {
 121         if ( cells[c2c[_cells[S.id_cell2]][0]] == S.id_cell1 ) {
 122           S.p[2] = pts2[2];
 123           p2 = true;
 124         }
 125       }
 126       if ( c2c[_cells[S.id_cell2]][1] != -1 ) {
 127         if ( cells[c2c[_cells[S.id_cell2]][1]] == S.id_cell1 ) {
 128           S.p[2] = pts2[0];
 129           p2 = true;
 130         }
 131       }
 132       if ( c2c[_cells[S.id_cell2]][2] != -1 ) {
 133         if ( cells[c2c[_cells[S.id_cell2]][2]] == S.id_cell1 ) {
 134           S.p[2] = pts2[1];
 135           p2 = true;
 136         }
 137       }
 138
 139       if ( !p2 ) {//failed to place point 2, appears when cell1 is linked to cell2, but cell2 not to cell1
 140         cout << "S.id_cell1=" << S.id_cell1 << endl;
 141         cout << "S.id_cell2=" << S.id_cell2 << endl;
 142         GuiMainWindow::pointer()->saveAs( GuiMainWindow::pointer()->getFilePath() + "abort.egc", false );
 143         EG_BUG;
 144       }
 145     }
 146   }//end of if neighbour cell
 147   return S;
 148 }
 149
 150 int SurfaceOperation::UpdateCurrentMeshDensity()
 151 {
 152   if ( DebugLevel > 0 ) {
 153     cout << "===UpdateMeshDensity START===" << endl;
 154   }
 155   QVector<vtkIdType> cells;
 156   getAllSurfaceCells( cells, grid );
 157   EG_VTKDCC( vtkIntArray, cell_code, grid, "cell_code" );
 158   EG_VTKDCN( vtkDoubleArray, node_meshdensity_desired, grid, "node_meshdensity_desired" );
 159   setGrid( grid );
 160   setCells( cells );
 161   if ( DebugLevel > 5 ) {
 162     cout << "cells.size()=" << cells.size() << endl;
 163   }
 164   EG_VTKDCN( vtkDoubleArray, node_meshdensity_current, grid, "node_meshdensity_current" );
 165   l2g_t nodes = getPartNodes();
 166   foreach( vtkIdType node, nodes ) {
 167     node_meshdensity_current->SetValue( node, CurrentMeshDensity( node ) );
 168   }
 169   if ( DebugLevel > 0 ) {
 170     cout << "===UpdateMeshDensity END===" << endl;
 171   }
 172   return( 0 ); ///@@@ what for???
 173 }
 174
 175 int SurfaceOperation::UpdatePotentialSnapPoints( bool update_node_types, bool allow_feature_edge_vertices )
 176 {
 177   l2g_t nodes  = getPartNodes();
 178   l2g_t cells  = getPartCells();
 179   g2l_t _cells = getPartLocalCells();
 180   l2l_t c2c    = getPartC2C();
 181
 182   //cout<<"=== UpdatePotentialSnapPoints START ==="<<endl;
 183   //prepare
 184   setAllSurfaceCells();
 185
 186   m_PotentialSnapPoints.resize( nodes.size() );
 187
 188   //initialize default values
 189   EG_VTKDCN( vtkCharArray, node_type, grid, "node_type" );
 190   foreach( vtkIdType id_node, nodes ) {
 191     if ( update_node_types ) node_type->SetValue( id_node, VTK_SIMPLE_VERTEX );
 192     m_PotentialSnapPoints[id_node].clear();
 193   }
 194
 195   //cout<<"===pre-processing==="<<endl;
 196   int num_edges = 0;
 197   //We loop through edges
 198   foreach( vtkIdType id_cell, cells ) {
 199     vtkIdType *pts, Npts;
 200     grid->GetCellPoints( id_cell, Npts, pts );
 201     for ( int i = 0; i < Npts; i++ ) {
 202
 203       int i_neighbour_cell = c2c[_cells[id_cell]][i];
 204       if ( i_neighbour_cell >= 0 && cells[i_neighbour_cell] < id_cell ) continue;//already visited edge
 205       num_edges++;
 206
 207       vtkIdType id_node1 = pts[i];
 208       vtkIdType id_node2 = pts[( i+1 )%Npts];
 209
 210       //-----------------------
 211       //determine edge type
 212       char edge = getEdgeType( id_node2, id_node1, allow_feature_edge_vertices );
 213       //-----------------------
 214       //determine node type pre-processing (count nb of complex edges if the node is complex, otherwise, just count the nb of edges)
 215       if ( edge && node_type->GetValue( id_node1 ) == VTK_SIMPLE_VERTEX ) {
 216         m_PotentialSnapPoints[id_node1].clear();
 217         m_PotentialSnapPoints[id_node1].push_back( id_node2 );
 218         if ( update_node_types ) node_type->SetValue( id_node1, edge );
 219       }
 220       else if (( edge && node_type->GetValue( id_node1 ) == VTK_BOUNDARY_EDGE_VERTEX ) ||
 221                ( edge && node_type->GetValue( id_node1 ) == VTK_FEATURE_EDGE_VERTEX ) ||
 222                ( !edge && node_type->GetValue( id_node1 ) == VTK_SIMPLE_VERTEX ) ) {
 223         m_PotentialSnapPoints[id_node1].push_back( id_node2 );
 224         if ( node_type->GetValue( id_node1 ) && edge == VTK_BOUNDARY_EDGE_VERTEX ) {
 225           if ( update_node_types ) node_type->SetValue( id_node1, VTK_BOUNDARY_EDGE_VERTEX );//VTK_BOUNDARY_EDGE_VERTEX has priority over VTK_FEATURE_EDGE_VERTEX
 226         }
 227       }
 228
 229       if ( edge && node_type->GetValue( id_node2 ) == VTK_SIMPLE_VERTEX ) {
 230         m_PotentialSnapPoints[id_node2].clear();
 231         m_PotentialSnapPoints[id_node2].push_back( id_node1 );
 232         if ( update_node_types ) node_type->SetValue( id_node2, edge );
 233       }
 234       else if (( edge && node_type->GetValue( id_node2 ) == VTK_BOUNDARY_EDGE_VERTEX ) ||
 235                ( edge && node_type->GetValue( id_node2 ) == VTK_FEATURE_EDGE_VERTEX ) ||
 236                ( !edge && node_type->GetValue( id_node2 ) == VTK_SIMPLE_VERTEX ) ) {
 237         m_PotentialSnapPoints[id_node2].push_back( id_node1 );
 238         if ( node_type->GetValue( id_node2 ) && edge == VTK_BOUNDARY_EDGE_VERTEX ) {
 239           if ( update_node_types ) node_type->SetValue( id_node2, VTK_BOUNDARY_EDGE_VERTEX );//VTK_BOUNDARY_EDGE_VERTEX has priority over VTK_FEATURE_EDGE_VERTEX
 240         }
 241       }
 242     }
 243   }
 244
 245   //cout<<"num_edges="<<num_edges<<endl;
 246
 247   //-----------------------
 248   //determine node type post-processing
 249   double CosEdgeAngle = cos(( double ) vtkMath::RadiansFromDegrees( this->EdgeAngle ) );
 250   //cout<<"===post-processing==="<<endl;
 251   //This time, we loop through nodes
 252   foreach( vtkIdType id_node, nodes ) {
 253     if ( node_type->GetValue( id_node ) == VTK_FEATURE_EDGE_VERTEX || node_type->GetValue( id_node ) == VTK_BOUNDARY_EDGE_VERTEX ) { //see how many edges; if two, what the angle is
 254
 255       if ( !this->BoundarySmoothing &&
 256            node_type->GetValue( id_node ) == VTK_BOUNDARY_EDGE_VERTEX ) {
 257         if ( update_node_types ) node_type->SetValue( id_node, VTK_FIXED_VERTEX );
 258       }
 259       else if ( m_PotentialSnapPoints[id_node].size() != 2 ) {
 260         if ( update_node_types ) node_type->SetValue( id_node, VTK_FIXED_VERTEX );
 261       }
 262       else { //check angle between edges
 263         double x1[3], x2[3], x3[3], l1[3], l2[3];
 264         grid->GetPoint( m_PotentialSnapPoints[id_node][0], x1 );
 265         grid->GetPoint( id_node, x2 );
 266         grid->GetPoint( m_PotentialSnapPoints[id_node][1], x3 );
 267         for ( int k = 0; k < 3; k++ ) {
 268           l1[k] = x2[k] - x1[k];
 269           l2[k] = x3[k] - x2[k];
 270         }
 271         if ( vtkMath::Normalize( l1 ) >= 0.0 &&
 272              vtkMath::Normalize( l2 ) >= 0.0 &&
 273              vtkMath::Dot( l1, l2 ) < CosEdgeAngle ) {
 274           if ( update_node_types ) node_type->SetValue( id_node, VTK_FIXED_VERTEX );
 275         }
 276       }//if along edge
 277     }//if edge vertex
 278   }
 279   //cout<<"m_PotentialSnapPoints.size()="<<m_PotentialSnapPoints.size()<<endl;
 280   //cout<<"=== UpdatePotentialSnapPoints END ==="<<endl;
 281   return( 0 );
 282 }
 283
 284 char SurfaceOperation::getNodeType( vtkIdType id_node, bool allow_feature_edge_vertices )
 285 {
 286   l2g_t  nodes = getPartNodes();
 287   g2l_t _nodes = getPartLocalNodes();
 288   l2l_t  n2n   = getPartN2N();
 289
 290   //initialize default value
 291   char type = VTK_SIMPLE_VERTEX;
 292
 293   //loop through edges around id_node
 294
 295   QVector <vtkIdType> edges;
 296
 297   double CosEdgeAngle = cos(( double ) vtkMath::RadiansFromDegrees( this->EdgeAngle ) );
 298
 299   foreach( int i_node2, n2n[_nodes[id_node]] ) {
 300     vtkIdType id_node2 = nodes[i_node2];
 301     //-----------------------
 302     //determine edge type
 303     char edge = getEdgeType( id_node2, id_node, allow_feature_edge_vertices );
 304
 305     //-----------------------
 306     //determine node type pre-processing (count nb of complex edges if the node is complex, otherwise, just count the nb of edges)
 307     if ( edge && type == VTK_SIMPLE_VERTEX ) {
 308       edges.clear();
 309       edges.push_back( id_node2 );
 310       type = edge;
 311     }
 312     else if (( edge && type == VTK_BOUNDARY_EDGE_VERTEX ) ||
 313              ( edge && type == VTK_FEATURE_EDGE_VERTEX ) ||
 314              ( !edge && type == VTK_SIMPLE_VERTEX ) ) {
 315       edges.push_back( id_node2 );
 316       if ( type && edge == VTK_BOUNDARY_EDGE_VERTEX ) {
 317         type = VTK_BOUNDARY_EDGE_VERTEX;//VTK_BOUNDARY_EDGE_VERTEX has priority over VTK_FEATURE_EDGE_VERTEX
 318       }
 319     }
 320   }
 321   //-----------------------
 322   //determine node type post-processing
 323   if ( type == VTK_FEATURE_EDGE_VERTEX || type == VTK_BOUNDARY_EDGE_VERTEX ) { //see how many edges; if two, what the angle is
 324
 325     if ( !this->BoundarySmoothing &&
 326          type == VTK_BOUNDARY_EDGE_VERTEX ) {
 327       type = VTK_FIXED_VERTEX;
 328     }
 329     else if ( edges.size() != 2 ) {
 330       type = VTK_FIXED_VERTEX;
 331     }
 332     else { //check angle between edges
 333       double x1[3], x2[3], x3[3], l1[3], l2[3];
 334       grid->GetPoint( edges[0], x1 );
 335       grid->GetPoint( id_node, x2 );
 336       grid->GetPoint( edges[1], x3 );
 337       for ( int k = 0; k < 3; k++ ) {
 338         l1[k] = x2[k] - x1[k];
 339         l2[k] = x3[k] - x2[k];
 340       }
 341       if ( vtkMath::Normalize( l1 ) >= 0.0 &&
 342            vtkMath::Normalize( l2 ) >= 0.0 &&
 343            vtkMath::Dot( l1, l2 ) < CosEdgeAngle ) {
 344         type = VTK_FIXED_VERTEX;
 345       }
 346     }//if along edge
 347   }//if edge vertex
 348
 349   if ( !allow_feature_edge_vertices && type == VTK_FEATURE_EDGE_VERTEX ) EG_BUG;
 350   return( type );
 351 }
 352
 353 int SurfaceOperation::getEdgeCells( vtkIdType id_node1, vtkIdType id_node2, QVector <vtkIdType> &EdgeCells )
 354 {
 355   g2l_t _nodes = getPartLocalNodes();
 356   l2g_t cells  = getPartCells();
 357   l2l_t n2c    = getPartN2C();
 358
 359   QSet<vtkIdType> S1;
 360   foreach( int i, n2c[_nodes[id_node1]] ) {
 361     S1.insert( cells[i] );
 362   }
 363
 364   QSet<vtkIdType> S2;
 365   foreach( int i, n2c[_nodes[id_node2]] ) {
 366     S2.insert( cells[i] );
 367   }
 368
 369   S2.intersect( S1 );
 370   EdgeCells = Set2Vector( S2, false );
 371   return EdgeCells.size();
 372 }
 373
 374 int SurfaceOperation::getEdgeCells( vtkIdType id_node1, vtkIdType id_node2, QSet <vtkIdType> &EdgeCells )
 375 {
 376   g2l_t _nodes = getPartLocalNodes();
 377   l2g_t cells  = getPartCells();
 378   l2l_t n2c    = getPartN2C();
 379
 380   QSet<vtkIdType> S1;
 381   foreach( int i, n2c[_nodes[id_node1]] ) {
 382     S1.insert( cells[i] );
 383   }
 384
 385   QSet<vtkIdType> S2;
 386   foreach( int i, n2c[_nodes[id_node2]] ) {
 387     S2.insert( cells[i] );
 388   }
 389
 390   EdgeCells = S2.intersect( S1 );
 391   return EdgeCells.size();
 392 }
 393
 394 char SurfaceOperation::getEdgeType( vtkIdType a_node1, vtkIdType a_node2, bool allow_feature_edge_vertices )
 395 {
 396   double CosFeatureAngle = cos(( double ) vtkMath::RadiansFromDegrees( this->FeatureAngle ) );
 397
 398   //compute number of cells around edge [a_node,p2] and put them into neighbour_cells
 399   QVector <vtkIdType> neighbour_cells;
 400   int numNei = getEdgeCells( a_node1, a_node2, neighbour_cells ) - 1;
 401
 402   //set default value
 403   char edge = VTK_SIMPLE_VERTEX;
 404
 405   if ( numNei == 0 ) {
 406     edge = VTK_BOUNDARY_EDGE_VERTEX;
 407   }
 408   else if ( numNei >= 2 ) {
 409     qWarning() << "FATAL ERROR: edge belongs to more than 2 cells! This is not supported yet.";
 410     EG_BUG;
 411     edge = VTK_FEATURE_EDGE_VERTEX;
 412   }
 413   else if ( numNei == 1 ) {
 414     //check angle between cell1 and cell2 against FeatureAngle
 415     if ( allow_feature_edge_vertices && this->m_AllowFeatureEdgeVertices && CosAngle( grid, neighbour_cells[0], neighbour_cells[1] ) <= CosFeatureAngle ) {
 416       edge = VTK_FEATURE_EDGE_VERTEX;
 417     }
 418     //check the boundary codes
 419     EG_VTKDCC( vtkIntArray, cell_code, grid, "cell_code" );
 420     if ( cell_code->GetValue( neighbour_cells[0] ) !=  cell_code->GetValue( neighbour_cells[1] ) ) {
 421       edge = VTK_BOUNDARY_EDGE_VERTEX;
 422     }
 423   }
 424
 425   return( edge );
 426 }
 427
 428 QSet <int> SurfaceOperation::getBCset( vtkIdType id_node )
 429 {
 430   g2l_t _nodes = getPartLocalNodes();
 431   l2g_t  cells = getPartCells();
 432   l2l_t  n2c   = getPartN2C();
 433
 434   EG_VTKDCC( vtkIntArray, cell_code, grid, "cell_code" );
 435   QSet <int> bc;
 436   foreach( int i_cell, n2c[_nodes[id_node]] ) {
 437     vtkIdType id_cell = cells[i_cell];
 438     bc.insert( cell_code->GetValue( id_cell ) );
 439   }
 440   return( bc );
 441 }
 442
 443 VertexMeshDensity SurfaceOperation::getVMD( vtkIdType id_node )
 444 {
 445   g2l_t _nodes = getPartLocalNodes();
 446   l2g_t  cells = getPartCells();
 447   l2l_t  n2c   = getPartN2C();
 448
 449   EG_VTKDCN( vtkCharArray, node_type, grid, "node_type" );
 450   EG_VTKDCC( vtkIntArray, cell_code, grid, "cell_code" );
 451
 452   VertexMeshDensity VMD;
 453   VMD.type = node_type->GetValue( id_node );
 454   VMD.density = 0;
 455   VMD.CurrentNode = id_node;
 456
 457   foreach( int i_cell, n2c[_nodes[id_node]] ) {
 458     vtkIdType id_cell = cells[i_cell];
 459     VMD.BCmap[cell_code->GetValue( id_cell )] = 2;
 460   }
 461   return( VMD );
 462 }
 463
 464 //////////////////////////////////////////////
 465 double SurfaceOperation::CurrentVertexAvgDist( vtkIdType id_node )
 466 {
 467   l2g_t  nodes = getPartNodes();
 468   g2l_t _nodes = getPartLocalNodes();
 469   l2l_t  n2n   = getPartN2N();
 470
 471   double total_dist = 0;
 472   double avg_dist = 0;
 473   int N = n2n[_nodes[id_node]].size();
 474   vec3_t C;
 475   grid->GetPoint( id_node, C.data() );
 476   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 477     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 478     vec3_t M;
 479     grid->GetPoint( id_node_neighbour, M.data() );
 480     total_dist += ( M - C ).abs();
 481   }
 482   avg_dist = total_dist / ( double )N;
 483   return( avg_dist );
 484 }
 485
 486 double SurfaceOperation::CurrentMeshDensity( vtkIdType id_node )
 487 {
 488   return 1.0 / CurrentVertexAvgDist( id_node );
 489 }
 490
 491 double SurfaceOperation::DesiredVertexAvgDist( vtkIdType id_node )
 492 {
 493   l2g_t  nodes = getPartNodes();
 494   g2l_t _nodes = getPartLocalNodes();
 495   l2l_t  n2n   = getPartN2N();
 496
 497   double total_dist = 0;
 498   double avg_dist = 0;
 499   EG_VTKDCN( vtkDoubleArray, node_meshdensity_desired, grid, "node_meshdensity_desired" );
 500   int N = n2n[_nodes[id_node]].size();
 501   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 502     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 503     total_dist += 1. / node_meshdensity_desired->GetValue( id_node_neighbour );
 504   }
 505   avg_dist = total_dist / ( double )N;
 506   return( avg_dist );
 507 }
 508
 509 double SurfaceOperation::DesiredMeshDensity( vtkIdType id_node )
 510 {
 511   l2g_t  nodes = getPartNodes();
 512   g2l_t _nodes = getPartLocalNodes();
 513   l2l_t  n2n   = getPartN2N();
 514
 515   double total_density = 0;
 516   double avg_density = 0;
 517   EG_VTKDCN( vtkDoubleArray, node_meshdensity_desired, grid, "node_meshdensity_desired" );
 518   int N = n2n[_nodes[id_node]].size();
 519   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 520     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 521     total_density += node_meshdensity_desired->GetValue( id_node_neighbour );
 522   }
 523   avg_density = total_density / ( double )N;
 524   return( avg_density );
 525 }
 526
 527 ///////////////////////////////////////////
 528
 529 //---------------------------------------------------
 530 //Utility functions used in Roland's formulas
 531
 532 ///@@@ TODO: change meshdensity fields to edgelength fields since this is what is mostly used?
 533
 534 /// desired edge length for id_node
 535 double SurfaceOperation::desiredEdgeLength( vtkIdType id_node )
 536 {
 537   EG_VTKDCN( vtkDoubleArray, node_meshdensity_desired, grid, "node_meshdensity_desired" );
 538   return( 1.0 / node_meshdensity_desired->GetValue( id_node ) );
 539 }
 540
 541 //other functions
 542 ///perimeter
 543 double SurfaceOperation::perimeter( vtkIdType id_cell )
 544 {
 545   double ret = 0;
 546   vtkIdType num_pts, *pts;
 547   grid->GetCellPoints( id_cell, num_pts, pts );
 548   for ( int i = 0; i < num_pts; i++ ) {
 549     vec3_t A, B;
 550     grid->GetPoints()->GetPoint( pts[i], A.data() );
 551     grid->GetPoints()->GetPoint( pts[( i+1 )%num_pts], B.data() );
 552     ret += ( B - A ).abs();
 553   }
 554   return( ret );
 555 }
 556
 557 /// mean desired edge length for id_cell
 558 double SurfaceOperation::meanDesiredEdgeLength( vtkIdType id_cell )
 559 {
 560   vtkIdType num_pts, *pts;
 561   grid->GetCellPoints( id_cell, num_pts, pts );
 562   int total = 0;
 563   for ( int i = 0; i < num_pts; i++ ) {
 564     total += desiredEdgeLength( pts[i] );
 565   }
 566   return total / ( double )num_pts;
 567 }
 568
 569 ///@@@ TODO: Should be renamed to be more explicit if possible
 570
 571 /// perimeter / sum of the desired edge lengths
 572 double SurfaceOperation::Q_L( vtkIdType id_cell )
 573 {
 574   double denom_sum = 0;
 575   vtkIdType num_pts, *pts;
 576   grid->GetCellPoints( id_cell, num_pts, pts );
 577   for ( int i = 0; i < num_pts; i++ ) {
 578     denom_sum += desiredEdgeLength( pts[i] );
 579   }
 580   return( perimeter( id_cell ) / denom_sum );
 581 }
 582
 583 /// sum(2*edgelength,edges(id_node))/sum(desired edgelengths of each edgepoint,edges(id_node))
 584 double SurfaceOperation::Q_L1( vtkIdType id_node )
 585 {
 586   l2l_t n2n = getPartN2N();
 587   g2l_t _nodes = getPartLocalNodes();
 588   l2g_t nodes = getPartNodes();
 589
 590   double num_sum = 0;
 591   double denom_sum = 0;
 592   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 593     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 594     num_sum += 2 * distance( grid, id_node_neighbour, id_node );
 595     denom_sum += desiredEdgeLength( id_node ) + desiredEdgeLength( id_node_neighbour );
 596   }
 597   return( num_sum / denom_sum );
 598 }
 599
 600 /// minimum of sum(2*edgelength)/sum(desired edgelengths of each edgepoint) for each edge of id_node
 601 double SurfaceOperation::Q_L2( vtkIdType id_node )
 602 {
 603   l2l_t n2n = getPartN2N();
 604   g2l_t _nodes = getPartLocalNodes();
 605   l2g_t nodes = getPartNodes();
 606
 607   QVector <double> V;
 608   double num, denom;
 609   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 610     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 611     num = 2 * distance( grid, id_node_neighbour, id_node );
 612     denom = desiredEdgeLength( id_node ) + desiredEdgeLength( id_node_neighbour );
 613     V.push_back( num / denom );
 614   }
 615   qSort( V.begin(), V.end() );
 616   return( V[0] );
 617 }
 618
 619 /// Value to minimize for mesh smoothing. w allows putting more weight on the form or the area of triangles.
 620 double SurfaceOperation::T_min( int w )
 621 {
 622   l2g_t cells = getPartCells();
 623   double T = 0;
 624   foreach( vtkIdType id_cell, cells ) {
 625     T += areaOfCircumscribedCircle( grid, id_cell ) / pow( cellVA( grid, id_cell ), w ) * pow( meanDesiredEdgeLength( id_cell ), 2 * ( w - 1 ) );
 626   }
 627   return( T );
 628 }
 629
 630 //---------------------------------------------------
 631
 632 vtkIdType SurfaceOperation::getClosestNode( vtkIdType id_node )
 633 {
 634   l2l_t n2n = getPartN2N();
 635   g2l_t _nodes = getPartLocalNodes();
 636   l2g_t nodes = getPartNodes();
 637
 638   vec3_t C;
 639   grid->GetPoint( id_node, C.data() );
 640   vtkIdType id_minlen = -1;
 641   double minlen = -1;
 642   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 643     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 644     vec3_t M;
 645     grid->GetPoint( id_node_neighbour, M.data() );
 646     double len = ( M - C ).abs();
 647     if ( minlen < 0 or len < minlen ) {
 648       minlen = len;
 649       id_minlen = id_node_neighbour;
 650     }
 651   }
 652   return( id_minlen );
 653 }
 654
 655 vtkIdType SurfaceOperation::getFarthestNode( vtkIdType id_node )
 656 {
 657   l2l_t n2n = getPartN2N();
 658   g2l_t _nodes = getPartLocalNodes();
 659   l2g_t nodes = getPartNodes();
 660
 661   vec3_t C;
 662   grid->GetPoint( id_node, C.data() );
 663   vtkIdType id_maxlen = -1;
 664   double maxlen = -1;
 665   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 666     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 667     vec3_t M;
 668     grid->GetPoint( id_node_neighbour, M.data() );
 669     double len = ( M - C ).abs();
 670     if ( maxlen < 0 or len > maxlen ) {
 671       maxlen = len;
 672       id_maxlen = id_node_neighbour;
 673     }
 674   }
 675   return( id_maxlen );
 676 }
 677
 678 QVector <vtkIdType> SurfaceOperation::getPotentialSnapPoints( vtkIdType id_node )
 679 {
 680   if ( id_node < 0 || id_node >= m_PotentialSnapPoints.size() ) EG_BUG;
 681   return m_PotentialSnapPoints[id_node];
 682 }
 683
 684 // bool SurfaceOperation::DeletePoint( vtkIdType DeadNode, int& num_newpoints, int& num_newcells )
 685 // {
 686 //   QSet <vtkIdType> DeadNodes;
 687 //   DeadNodes.insert( DeadNode );
 688 //   bool ret = DeleteSetOfPoints( DeadNodes, num_newpoints, num_newcells );
 689 //   return( ret );
 690 // }
 691 // //End of DeletePoint