src/ADAPTFrontTrack/FrontTrack_NodesOnGeom.cxx

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