src/surfaceoperation.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2010 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   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
 153   m_PotentialSnapPoints.resize(m_Grid->GetNumberOfPoints());
 154
 155   //initialize default values
 156   EG_VTKDCN( vtkCharArray, node_type, m_Grid, "node_type" );
 157   foreach( vtkIdType id_node, nodes ) {
 158     if ( update_node_types ) node_type->SetValue( id_node, VTK_SIMPLE_VERTEX );
 159     m_PotentialSnapPoints[id_node].clear();
 160   }
 161
 162   //cout<<"===pre-processing==="<<endl;
 163   int num_edges = 0;
 164   //We loop through edges
 165   foreach( vtkIdType id_cell, cells ) {
 166     vtkIdType *pts, Npts;
 167     m_Grid->GetCellPoints( id_cell, Npts, pts );
 168     for ( int i = 0; i < Npts; i++ ) {
 169
 170       QVector <vtkIdType> EdgeCells;
 171       int Ncells = getEdgeCells( pts[i], pts[(i+1)%Npts], EdgeCells );
 172
 173       bool already_visited_edge = false;
 174       for(int i_cell=0;i_cell<Ncells;i_cell++) {
 175         if(EdgeCells[i_cell]<id_cell) {
 176           already_visited_edge = true;
 177           break;
 178         }
 179       }
 180
 181       if ( already_visited_edge ) continue;//already visited edge
 182       num_edges++;
 183
 184       vtkIdType id_node1 = pts[i];
 185       vtkIdType id_node2 = pts[( i+1 )%Npts];
 186
 187       //-----------------------
 188       //determine edge type
 189       char edge = getEdgeType( id_node2, id_node1, fix_unselected );
 190       //-----------------------
 191       //determine node type pre-processing (count nb of complex edges if the node is complex, otherwise, just count the nb of edges)
 192       if ( edge && node_type->GetValue( id_node1 ) == VTK_SIMPLE_VERTEX ) {
 193         m_PotentialSnapPoints[id_node1].clear();
 194         m_PotentialSnapPoints[id_node1].push_back( id_node2 );
 195         if ( update_node_types ) node_type->SetValue( id_node1, edge );
 196       }
 197       else if (( edge && node_type->GetValue( id_node1 ) == VTK_BOUNDARY_EDGE_VERTEX ) ||
 198                ( edge && node_type->GetValue( id_node1 ) == VTK_FEATURE_EDGE_VERTEX ) ||
 199                ( !edge && node_type->GetValue( id_node1 ) == VTK_SIMPLE_VERTEX ) ) {
 200         m_PotentialSnapPoints[id_node1].push_back( id_node2 );
 201         if ( node_type->GetValue( id_node1 ) && edge == VTK_BOUNDARY_EDGE_VERTEX ) {
 202           if ( update_node_types ) node_type->SetValue( id_node1, VTK_BOUNDARY_EDGE_VERTEX );//VTK_BOUNDARY_EDGE_VERTEX has priority over VTK_FEATURE_EDGE_VERTEX
 203         }
 204       }
 205
 206       if ( edge && node_type->GetValue( id_node2 ) == VTK_SIMPLE_VERTEX ) {
 207         m_PotentialSnapPoints[id_node2].clear();
 208         m_PotentialSnapPoints[id_node2].push_back( id_node1 );
 209         if ( update_node_types ) node_type->SetValue( id_node2, edge );
 210       }
 211       else if (( edge && node_type->GetValue( id_node2 ) == VTK_BOUNDARY_EDGE_VERTEX ) ||
 212                ( edge && node_type->GetValue( id_node2 ) == VTK_FEATURE_EDGE_VERTEX ) ||
 213                ( !edge && node_type->GetValue( id_node2 ) == VTK_SIMPLE_VERTEX ) ) {
 214         m_PotentialSnapPoints[id_node2].push_back( id_node1 );
 215         if ( node_type->GetValue( id_node2 ) && edge == VTK_BOUNDARY_EDGE_VERTEX ) {
 216           if ( update_node_types ) node_type->SetValue( id_node2, VTK_BOUNDARY_EDGE_VERTEX );//VTK_BOUNDARY_EDGE_VERTEX has priority over VTK_FEATURE_EDGE_VERTEX
 217         }
 218       }
 219     }
 220   }
 221
 222   //cout<<"num_edges="<<num_edges<<endl;
 223
 224   //-----------------------
 225   //determine node type post-processing
 226   double CosEdgeAngle = cos(this->m_EdgeAngle);
 227   //cout<<"===post-processing==="<<endl;
 228   //This time, we loop through nodes
 229   foreach( vtkIdType id_node, nodes ) {
 230     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
 231
 232       if ( !this->m_BoundarySmoothing && node_type->GetValue( id_node ) == VTK_BOUNDARY_EDGE_VERTEX ) {
 233         if ( update_node_types ) node_type->SetValue( id_node, VTK_FIXED_VERTEX );
 234       } else if ( m_PotentialSnapPoints[id_node].size() != 2 ) {
 235         if ( update_node_types ) node_type->SetValue( id_node, VTK_FIXED_VERTEX );
 236       }
 237       else { //check angle between edges
 238         double x1[3], x2[3], x3[3], l1[3], l2[3];
 239         m_Grid->GetPoint( m_PotentialSnapPoints[id_node][0], x1 );
 240         m_Grid->GetPoint( id_node, x2 );
 241         m_Grid->GetPoint( m_PotentialSnapPoints[id_node][1], x3 );
 242         for ( int k = 0; k < 3; k++ ) {
 243           l1[k] = x2[k] - x1[k];
 244           l2[k] = x3[k] - x2[k];
 245         }
 246         if ( vtkMath::Normalize( l1 ) >= 0.0 &&
 247              vtkMath::Normalize( l2 ) >= 0.0 &&
 248              vtkMath::Dot( l1, l2 ) < CosEdgeAngle ) {
 249           if ( update_node_types ) node_type->SetValue( id_node, VTK_FIXED_VERTEX );
 250         }
 251       }//if along edge
 252     }//if edge vertex
 253   }
 254   //cout<<"m_PotentialSnapPoints.size()="<<m_PotentialSnapPoints.size()<<endl;
 255   //cout<<"=== UpdatePotentialSnapPoints END ==="<<endl;
 256   return( 0 );
 257 }
 258
 259 char SurfaceOperation::getNodeType( vtkIdType id_node, bool fix_unselected )
 260 {
 261   l2g_t  nodes = getPartNodes();
 262   g2l_t _nodes = getPartLocalNodes();
 263   l2l_t  n2n   = getPartN2N();
 264
 265   //initialize default value
 266   char type = VTK_SIMPLE_VERTEX;
 267
 268   //loop through edges around id_node
 269
 270   QVector <vtkIdType> edges;
 271
 272   double CosEdgeAngle = cos(this->m_EdgeAngle);
 273
 274   foreach( int i_node2, n2n[_nodes[id_node]] ) {
 275     vtkIdType id_node2 = nodes[i_node2];
 276     //-----------------------
 277     //determine edge type
 278     char edge = getEdgeType(id_node2, id_node, fix_unselected);
 279
 280     //-----------------------
 281     //determine node type pre-processing (count nb of complex edges if the node is complex, otherwise, just count the nb of edges)
 282     if ( edge && type == VTK_SIMPLE_VERTEX ) {
 283       edges.clear();
 284       edges.push_back( id_node2 );
 285       type = edge;
 286     }
 287     else if (( edge && type == VTK_BOUNDARY_EDGE_VERTEX ) ||
 288              ( edge && type == VTK_FEATURE_EDGE_VERTEX ) ||
 289              ( !edge && type == VTK_SIMPLE_VERTEX ) ) {
 290       edges.push_back( id_node2 );
 291       if ( type && edge == VTK_BOUNDARY_EDGE_VERTEX ) {
 292         type = VTK_BOUNDARY_EDGE_VERTEX;//VTK_BOUNDARY_EDGE_VERTEX has priority over VTK_FEATURE_EDGE_VERTEX
 293       }
 294     }
 295   }
 296   //-----------------------
 297   //determine node type post-processing
 298   if ( type == VTK_FEATURE_EDGE_VERTEX || type == VTK_BOUNDARY_EDGE_VERTEX ) { //see how many edges; if two, what the angle is
 299
 300     if ( !this->m_BoundarySmoothing && type == VTK_BOUNDARY_EDGE_VERTEX ) {
 301       type = VTK_FIXED_VERTEX;
 302     }
 303     else if ( edges.size() != 2 ) {
 304       type = VTK_FIXED_VERTEX;
 305     }
 306     else { //check angle between edges
 307       double x1[3], x2[3], x3[3], l1[3], l2[3];
 308       m_Grid->GetPoint( edges[0], x1 );
 309       m_Grid->GetPoint( id_node, x2 );
 310       m_Grid->GetPoint( edges[1], x3 );
 311       for ( int k = 0; k < 3; k++ ) {
 312         l1[k] = x2[k] - x1[k];
 313         l2[k] = x3[k] - x2[k];
 314       }
 315       if ( vtkMath::Normalize( l1 ) >= 0.0 &&
 316            vtkMath::Normalize( l2 ) >= 0.0 &&
 317            vtkMath::Dot( l1, l2 ) < CosEdgeAngle ) {
 318              type = VTK_FIXED_VERTEX;
 319       }
 320     }//if along edge
 321   }//if edge vertex
 322
 323   return( type );
 324 }
 325
 326 int SurfaceOperation::getEdgeCells( vtkIdType id_node1, vtkIdType id_node2, QVector <vtkIdType> &EdgeCells )
 327 {
 328   g2l_t _nodes = getPartLocalNodes();
 329   l2g_t cells  = getPartCells();
 330   l2l_t n2c    = getPartN2C();
 331
 332   QSet<vtkIdType> S1;
 333   foreach( int i, n2c[_nodes[id_node1]] ) {
 334     S1.insert( cells[i] );
 335   }
 336
 337   QSet<vtkIdType> S2;
 338   foreach( int i, n2c[_nodes[id_node2]] ) {
 339     S2.insert( cells[i] );
 340   }
 341
 342   S2.intersect( S1 );
 343   EdgeCells = Set2Vector( S2, false );
 344   return EdgeCells.size();
 345 }
 346
 347 int SurfaceOperation::getEdgeCells( vtkIdType id_node1, vtkIdType id_node2, QSet <vtkIdType> &EdgeCells )
 348 {
 349   g2l_t _nodes = getPartLocalNodes();
 350   l2g_t cells  = getPartCells();
 351   l2l_t n2c    = getPartN2C();
 352
 353   QSet<vtkIdType> S1;
 354   foreach( int i, n2c[_nodes[id_node1]] ) {
 355     S1.insert( cells[i] );
 356   }
 357
 358   QSet<vtkIdType> S2;
 359   foreach( int i, n2c[_nodes[id_node2]] ) {
 360     S2.insert( cells[i] );
 361   }
 362
 363   EdgeCells = S2.intersect( S1 );
 364   return EdgeCells.size();
 365 }
 366
 367 char SurfaceOperation::getEdgeType(vtkIdType a_node1, vtkIdType a_node2, bool fix_unselected)
 368 {
 369   double CosFeatureAngle = cos(this->m_FeatureAngle);
 370   bool feature_edges_disabled = m_FeatureAngle >= M_PI;
 371
 372   //compute number of cells around edge [a_node,p2] and put them into neighbour_cells
 373   QVector <vtkIdType> neighbour_cells;
 374   int numNei = getEdgeCells( a_node1, a_node2, neighbour_cells ) - 1;
 375
 376   //set default value
 377   char edge = VTK_SIMPLE_VERTEX;
 378
 379   if ( numNei == 0 ) {
 380     edge = VTK_BOUNDARY_EDGE_VERTEX;
 381   }
 382   else if ( numNei >= 2 ) {
 383     //qWarning() << "FATAL ERROR: edge belongs to more than 2 cells! This is not supported yet.";
 384     //EG_BUG;
 385     //edge = VTK_FEATURE_EDGE_VERTEX;
 386     edge = VTK_BOUNDARY_EDGE_VERTEX;
 387   }
 388   else if ( numNei == 1 ) {
 389     //check angle between cell1 and cell2 against FeatureAngle
 390     if (CosAngle(m_Grid, neighbour_cells[0], neighbour_cells[1] ) <= CosFeatureAngle && !feature_edges_disabled) {
 391       edge = VTK_FEATURE_EDGE_VERTEX;
 392     }
 393     //check the boundary codes
 394     EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
 395     int cell_code_0 = cell_code->GetValue( neighbour_cells[0] );
 396     int cell_code_1 = cell_code->GetValue( neighbour_cells[1] );
 397     if ( cell_code_0 !=  cell_code_1 ) {
 398       edge = VTK_BOUNDARY_EDGE_VERTEX;
 399     }
 400 //     qWarning()<<"m_BoundaryCodes="<<m_BoundaryCodes;
 401     if(m_BoundaryCodes.isEmpty()) {
 402       EG_ERR_RETURN("no boundary codes specified");
 403     }
 404     if(fix_unselected) {
 405       if( !m_BoundaryCodes.contains(cell_code_0) || !m_BoundaryCodes.contains(cell_code_1) ) {
 406         edge = VTK_FIXED_VERTEX;// does not make sense, but should make the points of the edge fixed
 407       }
 408     }
 409   }
 410
 411   return( edge );
 412 }
 413
 414 QSet <int> SurfaceOperation::getBCset( vtkIdType id_node )
 415 {
 416   g2l_t _nodes = getPartLocalNodes();
 417   l2g_t  cells = getPartCells();
 418   l2l_t  n2c   = getPartN2C();
 419
 420   EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
 421   QSet <int> bc;
 422   foreach( int i_cell, n2c[_nodes[id_node]] ) {
 423     vtkIdType id_cell = cells[i_cell];
 424     bc.insert( cell_code->GetValue( id_cell ) );
 425   }
 426   return( bc );
 427 }
 428
 429 VertexMeshDensity SurfaceOperation::getVMD( vtkIdType id_node )
 430 {
 431   g2l_t _nodes = getPartLocalNodes();
 432   l2g_t  cells = getPartCells();
 433   l2l_t  n2c   = getPartN2C();
 434
 435   EG_VTKDCN( vtkCharArray, node_type, m_Grid, "node_type" );
 436   EG_VTKDCC( vtkIntArray, cell_code, m_Grid, "cell_code" );
 437
 438   VertexMeshDensity VMD;
 439   VMD.type = node_type->GetValue( id_node );
 440   VMD.density = 0;
 441   VMD.CurrentNode = id_node;
 442
 443   foreach( int i_cell, n2c[_nodes[id_node]] ) {
 444     vtkIdType id_cell = cells[i_cell];
 445     VMD.BCmap[cell_code->GetValue( id_cell )] = 2;
 446   }
 447   return( VMD );
 448 }
 449
 450 //////////////////////////////////////////////
 451 double SurfaceOperation::CurrentVertexAvgDist( vtkIdType id_node )
 452 {
 453   l2g_t  nodes = getPartNodes();
 454   g2l_t _nodes = getPartLocalNodes();
 455   l2l_t  n2n   = getPartN2N();
 456
 457   double total_dist = 0;
 458   double avg_dist = 0;
 459   int N = n2n[_nodes[id_node]].size();
 460   vec3_t C;
 461   m_Grid->GetPoint( id_node, C.data() );
 462   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 463     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 464     vec3_t M;
 465     m_Grid->GetPoint( id_node_neighbour, M.data() );
 466     total_dist += ( M - C ).abs();
 467   }
 468   avg_dist = total_dist / ( double )N;
 469   return( avg_dist );
 470 }
 471
 472 double SurfaceOperation::CurrentMeshDensity( vtkIdType id_node )
 473 {
 474   return 1.0 / CurrentVertexAvgDist( id_node );
 475 }
 476
 477 double SurfaceOperation::DesiredVertexAvgDist( vtkIdType id_node )
 478 {
 479   l2g_t  nodes = getPartNodes();
 480   g2l_t _nodes = getPartLocalNodes();
 481   l2l_t  n2n   = getPartN2N();
 482
 483   double total_dist = 0;
 484   double avg_dist = 0;
 485   EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
 486   int N = n2n[_nodes[id_node]].size();
 487   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 488     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 489     total_dist += 1. / characteristic_length_desired->GetValue( id_node_neighbour );
 490   }
 491   avg_dist = total_dist / ( double )N;
 492   return( avg_dist );
 493 }
 494
 495 double SurfaceOperation::DesiredMeshDensity( vtkIdType id_node )
 496 {
 497   l2g_t  nodes = getPartNodes();
 498   g2l_t _nodes = getPartLocalNodes();
 499   l2l_t  n2n   = getPartN2N();
 500
 501   double total_density = 0;
 502   double avg_density = 0;
 503   EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
 504   int N = n2n[_nodes[id_node]].size();
 505   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 506     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 507     total_density += characteristic_length_desired->GetValue( id_node_neighbour );
 508   }
 509   avg_density = total_density / ( double )N;
 510   return( avg_density );
 511 }
 512
 513 ///////////////////////////////////////////
 514
 515 //---------------------------------------------------
 516 //Utility functions used in Roland's formulas
 517
 518 ///\todo change meshdensity fields to edgelength fields since this is what is mostly used?
 519
 520 /// desired edge length for id_node
 521 double SurfaceOperation::desiredEdgeLength( vtkIdType id_node )
 522 {
 523   EG_VTKDCN( vtkDoubleArray, characteristic_length_desired, m_Grid, "node_meshdensity_desired" );
 524   return( 1.0 / characteristic_length_desired->GetValue( id_node ) );
 525 }
 526
 527 //other functions
 528 ///perimeter
 529 double SurfaceOperation::perimeter( vtkIdType id_cell )
 530 {
 531   double ret = 0;
 532   vtkIdType num_pts, *pts;
 533   m_Grid->GetCellPoints( id_cell, num_pts, pts );
 534   for ( int i = 0; i < num_pts; i++ ) {
 535     vec3_t A, B;
 536     m_Grid->GetPoints()->GetPoint( pts[i], A.data() );
 537     m_Grid->GetPoints()->GetPoint( pts[( i+1 )%num_pts], B.data() );
 538     ret += ( B - A ).abs();
 539   }
 540   return( ret );
 541 }
 542
 543 /// mean desired edge length for id_cell
 544 double SurfaceOperation::meanDesiredEdgeLength( vtkIdType id_cell )
 545 {
 546   vtkIdType num_pts, *pts;
 547   m_Grid->GetCellPoints( id_cell, num_pts, pts );
 548   int total = 0;
 549   for ( int i = 0; i < num_pts; i++ ) {
 550     total += desiredEdgeLength( pts[i] );
 551   }
 552   return total / ( double )num_pts;
 553 }
 554
 555 ///\todo Should be renamed to be more explicit if possible
 556
 557 /// perimeter / sum of the desired edge lengths
 558 double SurfaceOperation::Q_L( vtkIdType id_cell )
 559 {
 560   double denom_sum = 0;
 561   vtkIdType num_pts, *pts;
 562   m_Grid->GetCellPoints( id_cell, num_pts, pts );
 563   for ( int i = 0; i < num_pts; i++ ) {
 564     denom_sum += desiredEdgeLength( pts[i] );
 565   }
 566   return( perimeter( id_cell ) / denom_sum );
 567 }
 568
 569 /// sum(2*edgelength,edges(id_node))/sum(desired edgelengths of each edgepoint,edges(id_node))
 570 double SurfaceOperation::Q_L1( vtkIdType id_node )
 571 {
 572   l2l_t n2n = getPartN2N();
 573   g2l_t _nodes = getPartLocalNodes();
 574   l2g_t nodes = getPartNodes();
 575
 576   double num_sum = 0;
 577   double denom_sum = 0;
 578   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 579     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 580     num_sum += 2 * distance( m_Grid, id_node_neighbour, id_node );
 581     denom_sum += desiredEdgeLength( id_node ) + desiredEdgeLength( id_node_neighbour );
 582   }
 583   return( num_sum / denom_sum );
 584 }
 585
 586 /// minimum of sum(2*edgelength)/sum(desired edgelengths of each edgepoint) for each edge of id_node
 587 double SurfaceOperation::Q_L2( vtkIdType id_node )
 588 {
 589   l2l_t n2n = getPartN2N();
 590   g2l_t _nodes = getPartLocalNodes();
 591   l2g_t nodes = getPartNodes();
 592
 593   QVector <double> V;
 594   double num, denom;
 595   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 596     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 597     num = 2 * distance( m_Grid, id_node_neighbour, id_node );
 598     denom = desiredEdgeLength( id_node ) + desiredEdgeLength( id_node_neighbour );
 599     V.push_back( num / denom );
 600   }
 601   qSort( V.begin(), V.end() );
 602   return( V[0] );
 603 }
 604
 605 /// Value to minimize for mesh smoothing. w allows putting more weight on the form or the area of triangles.
 606 double SurfaceOperation::T_min( int w )
 607 {
 608   l2g_t cells = getPartCells();
 609   double T = 0;
 610   foreach( vtkIdType id_cell, cells ) {
 611     T += areaOfCircumscribedCircle( m_Grid, id_cell ) / pow( cellVA( m_Grid, id_cell ), w ) * pow( meanDesiredEdgeLength( id_cell ), 2 * ( w - 1 ) );
 612   }
 613   return( T );
 614 }
 615
 616 //---------------------------------------------------
 617
 618 vtkIdType SurfaceOperation::getClosestNode( vtkIdType id_node )
 619 {
 620   l2l_t n2n = getPartN2N();
 621   g2l_t _nodes = getPartLocalNodes();
 622   l2g_t nodes = getPartNodes();
 623
 624   vec3_t C;
 625   m_Grid->GetPoint( id_node, C.data() );
 626   vtkIdType id_minlen = -1;
 627   double minlen = -1;
 628   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 629     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 630     vec3_t M;
 631     m_Grid->GetPoint( id_node_neighbour, M.data() );
 632     double len = ( M - C ).abs();
 633     if ( minlen < 0 or len < minlen ) {
 634       minlen = len;
 635       id_minlen = id_node_neighbour;
 636     }
 637   }
 638   return( id_minlen );
 639 }
 640
 641 vtkIdType SurfaceOperation::getFarthestNode( vtkIdType id_node )
 642 {
 643   l2l_t n2n = getPartN2N();
 644   g2l_t _nodes = getPartLocalNodes();
 645   l2g_t nodes = getPartNodes();
 646
 647   vec3_t C;
 648   m_Grid->GetPoint( id_node, C.data() );
 649   vtkIdType id_maxlen = -1;
 650   double maxlen = -1;
 651   foreach( int i_node_neighbour, n2n[_nodes[id_node]] ) {
 652     vtkIdType id_node_neighbour = nodes[i_node_neighbour];
 653     vec3_t M;
 654     m_Grid->GetPoint( id_node_neighbour, M.data() );
 655     double len = ( M - C ).abs();
 656     if ( maxlen < 0 or len > maxlen ) {
 657       maxlen = len;
 658       id_maxlen = id_node_neighbour;
 659     }
 660   }
 661   return( id_maxlen );
 662 }
 663
 664 QVector <vtkIdType> SurfaceOperation::getPotentialSnapPoints( vtkIdType id_node )
 665 {
 666   if ((id_node < 0) || (id_node >= m_PotentialSnapPoints.size())) {
 667     // UpdatePotentialSnapPoints should probably be called before using this function.
 668     EG_BUG;
 669   }
 670   return m_PotentialSnapPoints[id_node];
 671 }
 672
 673 bool SurfaceOperation::isCell(vtkIdType id_node1, vtkIdType id_node2, vtkIdType id_node3)
 674 {
 675   QVector <vtkIdType> EdgeCells_12;
 676   QVector <vtkIdType> EdgeCells_13;
 677   QVector <vtkIdType> inter;
 678
 679   getEdgeCells( id_node1, id_node2, EdgeCells_12 );
 680   getEdgeCells( id_node1, id_node3, EdgeCells_13 );
 681   qcontIntersection( EdgeCells_12, EdgeCells_13, inter );
 682   if(inter.size()>1) {
 683     qWarning()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
 684     qWarning()<<"EdgeCells_12="<<EdgeCells_12;
 685     qWarning()<<"EdgeCells_13="<<EdgeCells_13;
 686     qWarning()<<"inter="<<inter;
 687     writeGrid(m_Grid, "abort");
 688     EG_BUG;// multiple cells in the same place
 689   }
 690   if(DebugLevel>100) {
 691     qDebug()<<"(id_node1, id_node2, id_node3)="<<"("<<id_node1<<", "<<id_node2<<", "<<id_node3<<")";
 692     qDebug()<<"EdgeCells_12="<<EdgeCells_12;
 693     qDebug()<<"EdgeCells_13="<<EdgeCells_13;
 694     qDebug()<<"inter="<<inter;
 695   }
 696   return(inter.size()>0);
 697 }