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