src/FrontTrack/FrontTrack_NodesOnGeom.cxx

   1 // File      : FrontTrack_NodesOnGeom.cxx
   2 // Created   : Tue Apr 25 20:48:23 2017
   3 // Author    : Edward AGAPOV (eap)
   4
   5 #include "FrontTrack_NodesOnGeom.hxx"
   6 #include "FrontTrack_Utils.hxx"
   7
   8 #include <MEDCouplingMemArray.hxx>
   9
  10 #include <cstdio>
  11 #include <cstdlib>
  12 #include <list>
  13 #include <stdexcept>
  14
  15 namespace
  16 {
  17   /*!
  18    * \brief Close a file at destruction
  19    */
  20   struct FileCloser
  21   {
  22     FILE * _file;
  23
  24     FileCloser( FILE * file ): _file( file ) {}
  25     ~FileCloser() { if ( _file ) ::fclose( _file ); }
  26   };
  27 }
  28
  29 //================================================================================
  30 /*!
  31  * \brief Read node ids from a file and find shapes for projection
  32  *  \param [in] theNodeFile - a name of file holding IDs of nodes that
  33  *         will be moved onto geometry
  34  *  \param [in] theXaoGroups - a tool returning FT_Projector's by XAO group name
  35  *  \param [inout] theNodeCoords - array of node coordinates
  36  *  \param [in] theAllProjectorsByDim - all projectors of 2 dimensions, ordered so that
  37  *         a vector index corresponds to a XAO sub-shape ID
  38  */
  39 //================================================================================
  40
  41 void FT_NodesOnGeom::read( const std::string&            theNodeFile,
  42                            const FT_Utils::XaoGroups&    theXaoGroups,
  43                            MEDCoupling::DataArrayDouble* theNodeCoords,
  44                            std::vector< FT_Projector > * theAllProjectorsByDim )
  45 {
  46   _nodeCoords = theNodeCoords;
  47
  48   FILE * file = ::fopen( theNodeFile.c_str(), "r" );
  49   if ( !file )
  50     throw std::invalid_argument( "Can't open an input node file: " + theNodeFile );
  51
  52   FileCloser fileCloser( file );
  53
  54   // -------------------------------------
  55   // get shape dimension by the file name
  56   // -------------------------------------
  57
  58   // hope the file name is something like "frnD.**" with n in (1,2)
  59   int dimPos = theNodeFile.size() - 5;
  60   if ( theNodeFile[ dimPos ] == '2' )
  61     _shapeDim = 2;
  62   else if ( theNodeFile[ dimPos ] == '1' )
  63     _shapeDim = 1;
  64   else
  65     throw std::invalid_argument( "Can't define dimension by node file name " + theNodeFile );
  66 #ifdef _DEBUG_
  67   std::cout << ". Dimension of the file " << theNodeFile << ": " << _shapeDim << std::endl;
  68 #endif
  69
  70   // -------------------------------------
  71   // read geom group names; several lines
  72   // -------------------------------------
  73
  74   std::vector< std::string > geomNames;
  75
  76   const int maxLineLen = 256;
  77   char line[ maxLineLen ];
  78
  79   long int pos = ::ftell( file );
  80   while ( ::fgets( line, maxLineLen, file )) // read a line
  81   {
  82     if ( ::feof( file ))
  83     {
  84       return; // no nodes in the file
  85     }
  86
  87     // check if the line describes node ids in format 3I10 (e.g. "       120         1        43\n")
  88     size_t lineLen = strlen( line );
  89     if ( lineLen  >= 31        &&
  90          ::isdigit( line[9] )  &&
  91          line[10] == ' '       &&
  92          ::isdigit( line[19] ) &&
  93          line[20] == ' '       &&
  94          ::isdigit( line[29] ) &&
  95          ::isspace( line[30] ))
  96       break;
  97
  98     geomNames.push_back( line + 1 ); // skip the 1st white space
  99
 100     pos = ::ftell( file ); // remember the position to return if the next line holds node ids
 101   }
 102
 103   ::fseek( file, pos, SEEK_SET ); // return to the 1st line holding nodes ids
 104
 105
 106   // --------------
 107   // read node ids
 108   // --------------
 109
 110   FT_NodeToMove nodeIds;
 111   std::vector< int > ids;
 112
 113   const int nbNodes = theNodeCoords->getNumberOfTuples(); // to check validity of node IDs
 114
 115   while ( ::fgets( line, maxLineLen, file )) // read a line
 116   {
 117     // find node ids in the line
 118
 119     char *beg = line, *end = 0;
 120     long int id;
 121
 122     ids.clear();
 123     while (( id = ::strtol( beg, &end, 10 )) &&
 124            ( beg != end ))
 125     {
 126       ids.push_back( id );
 127       if ( id > nbNodes )
 128         throw std::invalid_argument( "Too large node ID: " + FT_Utils::toStr( id ));
 129       beg = end;
 130     }
 131
 132     if ( ids.size() >= 3 )
 133     {
 134       std::vector< int >::iterator i = ids.begin();
 135       nodeIds._nodeToMove = *i;
 136       nodeIds._neighborNodes.assign( ++i, ids.end() );
 137
 138       _nodes.push_back( nodeIds );
 139     }
 140
 141     if ( ::feof( file ))
 142       break;
 143   }
 144
 145   // -----------------------------------------------------------------
 146   // try to find FT_Projector's to boundary sub-shapes by group names
 147   // -----------------------------------------------------------------
 148
 149   _allProjectors = & theAllProjectorsByDim[ _shapeDim - 1 ];
 150
 151   _projectors.reserve( geomNames.size() );
 152   std::vector< const FT_Projector* >  projectors;
 153
 154   for ( size_t i = 0; i < geomNames.size(); ++i )
 155   {
 156     std::string & groupName = geomNames[i];
 157 #ifdef _DEBUG_
 158     std::cout << ". Group name: " << groupName << std::endl;
 159 #endif
 160
 161     // remove trailing white spaces
 162     for ( int iC = groupName.size() - 1; iC >= 0; --iC )
 163     {
 164       if ( ::isspace( groupName[iC] ) )
 165         groupName.resize( iC );
 166       else
 167         break;
 168     }
 169     if ( groupName.empty() )
 170       continue;
 171
 172     _groupNames.push_back( groupName ); // keep _groupNames for easier debug :)
 173
 174     // get projectors by group name
 175     theXaoGroups.getProjectors( groupName, _shapeDim,
 176                                 theAllProjectorsByDim[ _shapeDim-1 ], projectors );
 177   }
 178
 179   // ------------------------------
 180   // check the found FT_Projector's
 181   // ------------------------------
 182
 183   if ( projectors.size() == 1 )
 184   {
 185     _projectors.push_back( *projectors[ 0 ]);
 186   }
 187   else
 188   {
 189     Bnd_Box nodesBox;
 190     for ( size_t i = 0; i < _nodes.size(); ++i )
 191       nodesBox.Add( getPoint( _nodes[i]._nodeToMove ));
 192
 193     if ( projectors.size() > 1 )
 194     {
 195       // more than one boundary shape;
 196       // try to filter off unnecessary projectors using a bounding box of nodes
 197       for ( size_t i = 0; i < projectors.size(); ++i )
 198         if ( !nodesBox.IsOut( projectors[ i ]->getBoundingBox() ))
 199           _projectors.push_back( *projectors[ i ]);
 200     }
 201
 202     if ( _projectors.empty() )
 203     {
 204       // select projectors using a bounding box of nodes
 205       std::vector< FT_Projector > & allProjectors = *_allProjectors;
 206       for ( size_t i = 0; i < allProjectors.size(); ++i )
 207         if ( !nodesBox.IsOut( allProjectors[ i ].getBoundingBox() ))
 208           _projectors.push_back( allProjectors[ i ]);
 209
 210       if ( _projectors.empty() && !_nodes.empty() )
 211         throw std::runtime_error("No boundary shape found for nodes in file " + theNodeFile );
 212     }
 213   }
 214
 215   // prepare for projection - create real projectors
 216   for ( size_t i = 0; i < _projectors.size(); ++i )
 217     _projectors[ i ].prepareForProjection();
 218
 219 }
 220
 221 //================================================================================
 222 /*!
 223  * \brief Project nodes to the shapes and move them to new positions
 224  */
 225 //================================================================================
 226
 227 void FT_NodesOnGeom::projectAndMove()
 228 {
 229   _OK = true;
 230 //
 231 // 1. Préalables
 232 //
 233   // check if all the shapes are planar
 234   bool isAllPlanar = true;
 235   for ( size_t i = 0; i < _projectors.size() &&  isAllPlanar; ++i )
 236     isAllPlanar = _projectors[i].isPlanarBoundary();
 237   if ( isAllPlanar )
 238     return;
 239
 240   // set nodes in the order suitable for optimal projection
 241   putNodesInOrder();
 242
 243   // project and move nodes
 244
 245   std::vector< FT_NodeToMove* > notProjectedNodes;
 246   size_t iP, iProjector;
 247   gp_Pnt newXyz;
 248
 249 #ifdef _DEBUG_
 250     std::cout << ".. _projectors.size() = " << _projectors.size() << std::endl;
 251     std::cout << ".. _nodesOrder.size() = " << _nodesOrder.size() << std::endl;
 252 #endif
 253 //
 254 // 2. Calculs
 255 // 2.1. Avec plusieurs shapes
 256 //
 257   if ( _projectors.size() > 1 )
 258   {
 259     // the nodes are to be projected onto several boundary shapes;
 260     // in addition to the projecting, classification on a shape is necessary
 261     // in order to find out on which of the shapes a node is to be projected
 262
 263     iProjector = 0;
 264     for ( size_t i = 0; i < _nodesOrder.size(); ++i )
 265     {
 266       FT_NodeToMove& nn = _nodes[ _nodesOrder[ i ]];
 267       gp_Pnt        xyz = getPoint( nn._nodeToMove );
 268       gp_Pnt       xyz1 = getPoint( nn._neighborNodes[0] );
 269       gp_Pnt       xyz2 = getPoint( nn._neighborNodes[1] );
 270       double   maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.;
 271       if ( _projectors[ iProjector ].projectAndClassify( xyz, maxDist2, newXyz,
 272                                                          nn._params, nn._nearParams ))
 273       {
 274         moveNode( nn._nodeToMove, newXyz );
 275       }
 276       else // a node is not on iProjector-th shape, find the shape it is on
 277       {
 278         for ( iP = 1; iP < _projectors.size(); ++iP ) // check _projectors other than iProjector
 279         {
 280           iProjector = ( iProjector + 1 ) % _projectors.size();
 281           if ( _projectors[ iProjector ].projectAndClassify( xyz, maxDist2, newXyz,
 282                                                              nn._params, nn._nearParams ))
 283           {
 284             moveNode( nn._nodeToMove, newXyz );
 285             break;
 286           }
 287         }
 288         if ( iP == _projectors.size() )
 289         {
 290           notProjectedNodes.push_back( &nn );
 291
 292 #ifdef _DEBUG_
 293           std::cerr << "Warning: no shape found for node " << nn._nodeToMove << std::endl;
 294           if ( !_groupNames.empty() )
 295             std::cerr << "Warning:    group -- " << _groupNames[0] << std::endl;
 296 #endif
 297         }
 298       }
 299     }
 300   }
 301 //
 302 // 2.2. Avec une seule shape
 303 //
 304   else // one shape
 305   {
 306     for ( size_t i = 0; i < _nodesOrder.size(); ++i )
 307     {
 308       FT_NodeToMove& nn = _nodes[ _nodesOrder[ i ]];
 309       gp_Pnt        xyz = getPoint( nn._nodeToMove );
 310       gp_Pnt       xyz1 = getPoint( nn._neighborNodes[0] );
 311       gp_Pnt       xyz2 = getPoint( nn._neighborNodes[1] );
 312
 313 // maxDist2 : le quart du carré de la distance entre les deux voisins du noeud à bouger
 314       double   maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.;
 315 #ifdef _DEBUG_
 316     std::cout << "\n.. maxDist2 = " << maxDist2 << " entre " << nn._neighborNodes[0] << " et " << nn._neighborNodes[1] << " - milieu " << nn._nodeToMove << " - d/2 = " << sqrt(maxDist2) << " - d = " << sqrt(xyz1.SquareDistance( xyz2 )) << std::endl;
 317 #endif
 318       if ( _projectors[ 0 ].project( xyz, maxDist2, newXyz,
 319                                      nn._params, nn._nearParams ))
 320         moveNode( nn._nodeToMove, newXyz );
 321       else
 322         notProjectedNodes.push_back( &nn );
 323     }
 324   }
 325 //
 326 // 3. Bilan
 327 //
 328   if ( !notProjectedNodes.empty() )
 329   {
 330     // project nodes that are not projected by any of _projectors;
 331     // a proper projector is selected by evaluation of a distance between neighbor nodes
 332     // and a shape
 333
 334     std::vector< FT_Projector > & projectors = *_allProjectors;
 335
 336     iProjector = 0;
 337     for ( size_t i = 0; i < notProjectedNodes.size(); ++i )
 338     {
 339       FT_NodeToMove& nn = *notProjectedNodes[ i ];
 340       gp_Pnt        xyz = getPoint( nn._nodeToMove );
 341       gp_Pnt       xyz1 = getPoint( nn._neighborNodes[0] );
 342       gp_Pnt       xyz2 = getPoint( nn._neighborNodes[1] );
 343       double   maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.;
 344       double       tol2 = 1e-6 * maxDist2;
 345
 346       bool ok;
 347       for ( iP = 0; iP < projectors.size(); ++iP )
 348       {
 349         projectors[ iProjector ].prepareForProjection();
 350         projectors[ iProjector ].tryWithoutPrevSolution( true );
 351
 352         if (( ok = projectors[ iProjector ].isOnShape( xyz1, tol2, nn._params, nn._nearParams )) &&
 353             ( ok = projectors[ iProjector ].isOnShape( xyz2, tol2, nn._params, nn._params )))
 354         {
 355           if ( nn._neighborNodes.size() == 4 )
 356           {
 357             gp_Pnt xyz1 = getPoint( nn._neighborNodes[2] );
 358             gp_Pnt xyz2 = getPoint( nn._neighborNodes[3] );
 359             if (( ok = projectors[ iProjector ].isOnShape( xyz1, tol2, nn._params, nn._params )))
 360               ok     = projectors[ iProjector ].isOnShape( xyz2, tol2, nn._params, nn._params );
 361           }
 362         }
 363
 364         if ( ok && projectors[iProjector].project( xyz, maxDist2, newXyz, nn._params, nn._params ))
 365         {
 366           moveNode( nn._nodeToMove, newXyz );
 367           break;
 368         }
 369         iProjector = ( iProjector + 1 ) % projectors.size();
 370       }
 371       if ( iP == projectors.size() )
 372       {
 373         _OK = false;
 374
 375         std::cerr << "Error: not projected node " << nn._nodeToMove << std::endl;
 376       }
 377     }
 378   }
 379 }
 380
 381 //================================================================================
 382 /*!
 383  * \brief Put nodes in the order for optimal projection and set FT_NodeToMove::_nearParams
 384  *        to point to a FT_NodeToMove::_params of a node that will be projected earlier
 385  */
 386 //================================================================================
 387
 388 void FT_NodesOnGeom::putNodesInOrder()
 389 {
 390   if ( !_nodesOrder.empty() )
 391     return;
 392
 393   // check if any of projectors can use parameters of a previously projected node on a shape
 394   // to speed up projection
 395
 396   bool isPrevSolutionUsed = false;
 397   for ( size_t i = 0; i < _projectors.size() &&  !isPrevSolutionUsed; ++i )
 398     isPrevSolutionUsed = _projectors[i].canUsePrevSolution();
 399
 400   if ( !isPrevSolutionUsed )
 401   {
 402     _nodesOrder.resize( _nodes.size() );
 403     for ( size_t i = 0; i < _nodesOrder.size(); ++i )
 404       _nodesOrder[ i ] = i;
 405     return;
 406   }
 407
 408   // make a map to find a neighbor projected node
 409
 410   // map of { FT_NodeToMove::_neighborNodes[i] } to { FT_NodeToMove* };
 411   // here we call FT_NodeToMove a 'link' as this data links a _neighborNodes[i] node to other nodes
 412   typedef NCollection_DataMap< int, std::vector< FT_NodeToMove* > > TNodeIDToLinksMap;
 413   TNodeIDToLinksMap neigborsMap;
 414
 415   int mapSize = ( _shapeDim == 1 ) ? _nodes.size() + 1 : _nodes.size() * 3;
 416   neigborsMap.Clear();
 417   neigborsMap.ReSize( mapSize );
 418
 419   std::vector< FT_NodeToMove* > linkVec, *linkVecPtr;
 420   const int maxNbLinks = ( _shapeDim == 1 ) ? 2 : 6; // usual nb of links
 421
 422   for ( size_t i = 0; i < _nodes.size(); ++i )
 423   {
 424     FT_NodeToMove& nn = _nodes[i];
 425     for ( size_t iN = 0; iN < nn._neighborNodes.size(); ++iN )
 426     {
 427       if ( !( linkVecPtr = neigborsMap.ChangeSeek( nn._neighborNodes[ iN ] )))
 428       {
 429         linkVecPtr = neigborsMap.Bound( nn._neighborNodes[ iN ], linkVec );
 430         linkVecPtr->reserve( maxNbLinks );
 431       }
 432       linkVecPtr->push_back( & nn );
 433     }
 434   }
 435
 436   // fill in _nodesOrder
 437
 438   _nodesOrder.reserve( _nodes.size() );
 439
 440   std::list< FT_NodeToMove* > queue;
 441   queue.push_back( &_nodes[0] );
 442   _nodes[0]._nearParams = _nodes[0]._params; // to avoid re-adding to the queue
 443
 444   while ( !queue.empty() )
 445   {
 446     FT_NodeToMove* nn = queue.front();
 447     queue.pop_front();
 448
 449     _nodesOrder.push_back( nn - & _nodes[0] );
 450
 451     // add neighbors to the queue and set their _nearParams = nn->_params
 452     for ( size_t iN = 0; iN < nn->_neighborNodes.size(); ++iN )
 453     {
 454       std::vector< FT_NodeToMove* >& linkVec = neigborsMap( nn->_neighborNodes[ iN ]);
 455       for ( size_t iL = 0; iL < linkVec.size(); ++iL )
 456       {
 457         FT_NodeToMove* nnn = linkVec[ iL ];
 458         if ( nnn != nn && nnn->_nearParams == 0 )
 459         {
 460           nnn->_nearParams = nn->_params;
 461           queue.push_back( nnn );
 462         }
 463       }
 464     }
 465   }
 466   _nodes[0]._nearParams = 0; // reset
 467 }
 468
 469 //================================================================================
 470 /*!
 471  * \brief Get node coordinates. Node IDs count from a unit
 472  */
 473 //================================================================================
 474
 475 gp_Pnt FT_NodesOnGeom::getPoint( const int nodeID )
 476 {
 477   const size_t dim = _nodeCoords->getNumberOfComponents();
 478   const double * xyz = _nodeCoords->getConstPointer() + ( dim * ( nodeID - 1 ));
 479   return gp_Pnt( xyz[0], xyz[1], dim == 2 ? 0 : xyz[2] );
 480 }
 481
 482 //================================================================================
 483 /*!
 484  * \brief change node coordinates
 485  */
 486 //================================================================================
 487
 488 void FT_NodesOnGeom::moveNode( const int nodeID, const gp_Pnt& newXyz )
 489 {
 490   const size_t dim = _nodeCoords->getNumberOfComponents();
 491   double z, *xyz = _nodeCoords->getPointer() + ( dim * ( nodeID - 1 ));
 492   newXyz.Coord( xyz[0], xyz[1], dim == 2 ? z : xyz[2] );
 493 }