1 // File : FrontTrack_NodesOnGeom.cxx
2 // Created : Tue Apr 25 20:48:23 2017
3 // Author : Edward AGAPOV (eap)
5 #include "FrontTrack_NodesOnGeom.hxx"
6 #include "FrontTrack_Utils.hxx"
8 #include <MEDCouplingMemArray.hxx>
18 * \brief Close a file at destruction
24 FileCloser( FILE * file ): _file( file ) {}
25 ~FileCloser() { if ( _file ) ::fclose( _file ); }
29 //================================================================================
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
39 //================================================================================
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 )
46 _nodeCoords = theNodeCoords;
48 FILE * file = ::fopen( theNodeFile.c_str(), "r" );
50 throw std::invalid_argument( "Can't open an input node file: " + theNodeFile );
52 FileCloser fileCloser( file );
54 // -------------------------------------
55 // get shape dimension by the file name
56 // -------------------------------------
58 // hope the file name is something like "fr2D.**"
59 int dimPos = theNodeFile.size() - 5;
60 if ( theNodeFile[ dimPos ] == '2' )
62 else if ( theNodeFile[ dimPos ] == '1' )
65 throw std::invalid_argument( "Can't define dimension by node file name " + theNodeFile );
67 // -------------------------------------
68 // read geom group names; several lines
69 // -------------------------------------
71 std::vector< std::string > geomNames;
73 const int maxLineLen = 256;
74 char line[ maxLineLen ];
76 long int pos = ::ftell( file );
77 while ( ::fgets( line, maxLineLen, file )) // read a line
80 return; // no nodes in the file
82 // check if the line describes node ids in format 3I10 (e.g. " 120 1 43\n")
83 size_t lineLen = strlen( line );
85 ::isdigit( line[9] ) &&
87 ::isdigit( line[19] ) &&
89 ::isdigit( line[29] ) &&
90 ::isspace( line[30] ))
93 geomNames.push_back( line + 1 ); // skip the 1st white space
95 pos = ::ftell( file ); // remember the position to return if the next line holds node ids
98 ::fseek( file, pos, SEEK_SET ); // return to the 1st line holding nodes ids
105 FT_NodeToMove nodeIds;
106 std::vector< int > ids;
108 const int nbNodes = theNodeCoords->getNumberOfTuples(); // to check validity of node IDs
110 while ( ::fgets( line, maxLineLen, file )) // read a line
112 // find node ids in the line
114 char *beg = line, *end = 0;
118 while (( id = ::strtol( beg, &end, 10 )) &&
123 throw std::invalid_argument( "Too large node ID: " + FT_Utils::toStr( id ));
127 if ( ids.size() >= 3 )
129 std::vector< int >::iterator i = ids.begin();
130 nodeIds._nodeToMove = *i;
131 nodeIds._neighborNodes.assign( ++i, ids.end() );
133 _nodes.push_back( nodeIds );
140 // -----------------------------------------------------------------
141 // try to find FT_Projector's to boundary sub-shapes by group names
142 // -----------------------------------------------------------------
144 _allProjectors = & theAllProjectorsByDim[ _shapeDim - 1 ];
146 _projectors.reserve( geomNames.size() );
147 std::vector< const FT_Projector* > projectors;
149 for ( size_t i = 0; i < geomNames.size(); ++i )
151 std::string & groupName = geomNames[i];
153 // remove trailing white spaces
154 for ( int iC = groupName.size() - 1; iC >= 0; --iC )
156 if ( ::isspace( groupName[iC] ) )
157 groupName.resize( iC );
161 if ( groupName.empty() )
164 _groupNames.push_back( groupName ); // keep _groupNames for easier debug :)
166 // get projectors by group name
167 theXaoGroups.getProjectors( groupName, _shapeDim,
168 theAllProjectorsByDim[ _shapeDim-1 ], projectors );
171 // ------------------------------
172 // check the found FT_Projector's
173 // ------------------------------
175 if ( projectors.size() == 1 )
177 _projectors.push_back( *projectors[ 0 ]);
182 for ( size_t i = 0; i < _nodes.size(); ++i )
183 nodesBox.Add( getPoint( _nodes[i]._nodeToMove ));
185 if ( projectors.size() > 1 )
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 ]);
194 if ( _projectors.empty() )
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 ]);
202 if ( _projectors.empty() && !_nodes.empty() )
203 throw std::runtime_error("No boundary shape found for nodes in file " + theNodeFile );
207 // prepare for projection - create real projectors
208 for ( size_t i = 0; i < _projectors.size(); ++i )
209 _projectors[ i ].prepareForProjection();
213 //================================================================================
215 * \brief Project nodes to the shapes and move them to new positions
217 //================================================================================
219 void FT_NodesOnGeom::projectAndMove()
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();
230 // set nodes in the order suitable for optimal projection
233 // project and move nodes
235 std::vector< FT_NodeToMove* > notProjectedNodes;
236 size_t iP, iProjector;
240 std::cout << ".. _projectors.size() = " << _projectors.size() << std::endl;
241 std::cout << ".. _nodesOrder.size() = " << _nodesOrder.size() << std::endl;
243 if ( _projectors.size() > 1 )
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
250 for ( size_t i = 0; i < _nodesOrder.size(); ++i )
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 ))
260 moveNode( nn._nodeToMove, newXyz );
262 else // a node is not on iProjector-th shape, find the shape it is on
264 for ( iP = 1; iP < _projectors.size(); ++iP ) // check _projectors other than iProjector
266 iProjector = ( iProjector + 1 ) % _projectors.size();
267 if ( _projectors[ iProjector ].projectAndClassify( xyz, maxDist2, newXyz,
268 nn._params, nn._nearParams ))
270 moveNode( nn._nodeToMove, newXyz );
274 if ( iP == _projectors.size() )
276 notProjectedNodes.push_back( &nn );
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;
289 for ( size_t i = 0; i < _nodesOrder.size(); ++i )
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] );
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 );
302 notProjectedNodes.push_back( &nn );
307 if ( !notProjectedNodes.empty() )
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
313 std::vector< FT_Projector > & projectors = *_allProjectors;
316 for ( size_t i = 0; i < notProjectedNodes.size(); ++i )
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;
326 for ( iP = 0; iP < projectors.size(); ++iP )
328 projectors[ iProjector ].prepareForProjection();
329 projectors[ iProjector ].tryWithoutPrevSolution( true );
331 if (( ok = projectors[ iProjector ].isOnShape( xyz1, tol2, nn._params, nn._nearParams )) &&
332 ( ok = projectors[ iProjector ].isOnShape( xyz2, tol2, nn._params, nn._params )))
334 if ( nn._neighborNodes.size() == 4 )
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 );
343 if ( ok && projectors[iProjector].project( xyz, maxDist2, newXyz, nn._params, nn._params ))
345 moveNode( nn._nodeToMove, newXyz );
348 iProjector = ( iProjector + 1 ) % projectors.size();
350 if ( iP == projectors.size() )
354 std::cerr << "Error: not projected node " << nn._nodeToMove << std::endl;
360 //================================================================================
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
365 //================================================================================
367 void FT_NodesOnGeom::putNodesInOrder()
369 if ( !_nodesOrder.empty() )
372 // check if any of projectors can use parameters of a previously projected node on a shape
373 // to speed up projection
375 bool isPrevSolutionUsed = false;
376 for ( size_t i = 0; i < _projectors.size() && !isPrevSolutionUsed; ++i )
377 isPrevSolutionUsed = _projectors[i].canUsePrevSolution();
379 if ( !isPrevSolutionUsed )
381 _nodesOrder.resize( _nodes.size() );
382 for ( size_t i = 0; i < _nodesOrder.size(); ++i )
383 _nodesOrder[ i ] = i;
387 // make a map to find a neighbor projected node
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;
394 int mapSize = ( _shapeDim == 1 ) ? _nodes.size() + 1 : _nodes.size() * 3;
396 neigborsMap.ReSize( mapSize );
398 std::vector< FT_NodeToMove* > linkVec, *linkVecPtr;
399 const int maxNbLinks = ( _shapeDim == 1 ) ? 2 : 6; // usual nb of links
401 for ( size_t i = 0; i < _nodes.size(); ++i )
403 FT_NodeToMove& nn = _nodes[i];
404 for ( size_t iN = 0; iN < nn._neighborNodes.size(); ++iN )
406 if ( !( linkVecPtr = neigborsMap.ChangeSeek( nn._neighborNodes[ iN ] )))
408 linkVecPtr = neigborsMap.Bound( nn._neighborNodes[ iN ], linkVec );
409 linkVecPtr->reserve( maxNbLinks );
411 linkVecPtr->push_back( & nn );
415 // fill in _nodesOrder
417 _nodesOrder.reserve( _nodes.size() );
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
423 while ( !queue.empty() )
425 FT_NodeToMove* nn = queue.front();
428 _nodesOrder.push_back( nn - & _nodes[0] );
430 // add neighbors to the queue and set their _nearParams = nn->_params
431 for ( size_t iN = 0; iN < nn->_neighborNodes.size(); ++iN )
433 std::vector< FT_NodeToMove* >& linkVec = neigborsMap( nn->_neighborNodes[ iN ]);
434 for ( size_t iL = 0; iL < linkVec.size(); ++iL )
436 FT_NodeToMove* nnn = linkVec[ iL ];
437 if ( nnn != nn && nnn->_nearParams == 0 )
439 nnn->_nearParams = nn->_params;
440 queue.push_back( nnn );
445 _nodes[0]._nearParams = 0; // reset
448 //================================================================================
450 * \brief Get node coordinates. Node IDs count from a unit
452 //================================================================================
454 gp_Pnt FT_NodesOnGeom::getPoint( const int nodeID )
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] );
461 //================================================================================
463 * \brief change node coordinates
465 //================================================================================
467 void FT_NodesOnGeom::moveNode( const int nodeID, const gp_Pnt& newXyz )
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] );