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