1 // Copyright (C) 2017-2022 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // File : FrontTrack_NodesOnGeom.cxx
20 // Created : Tue Apr 25 20:48:23 2017
21 // Author : Edward AGAPOV (eap)
23 #include "FrontTrack_NodesOnGeom.hxx"
24 #include "FrontTrack_Utils.hxx"
26 #include <MEDCouplingMemArray.hxx>
36 * \brief Close a file at destruction
42 FileCloser( FILE * file ): _file( file ) {}
43 ~FileCloser() { if ( _file ) ::fclose( _file ); }
47 //================================================================================
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
57 //================================================================================
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 )
64 _nodeCoords = theNodeCoords;
66 FILE * file = ::fopen( theNodeFile.c_str(), "r" );
68 throw std::invalid_argument( "Can't open an input node file: " + theNodeFile );
70 FileCloser fileCloser( file );
72 // -------------------------------------
73 // get shape dimension by the file name
74 // -------------------------------------
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' )
80 else if ( theNodeFile[ dimPos ] == '1' )
83 throw std::invalid_argument( "Can't define dimension by node file name " + theNodeFile );
85 std::cout << ". Dimension of the file " << theNodeFile << ": " << _shapeDim << std::endl;
88 // -------------------------------------
89 // read geom group names; several lines
90 // -------------------------------------
92 std::vector< std::string > geomNames;
94 const int maxLineLen = 256;
95 char line[ maxLineLen ];
97 long int pos = ::ftell( file );
98 while ( ::fgets( line, maxLineLen, file )) // read a line
102 return; // no nodes in the file
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] ) &&
110 ::isdigit( line[19] ) &&
112 ::isdigit( line[29] ) &&
113 ::isspace( line[30] ))
116 geomNames.push_back( line + 1 ); // skip the 1st white space
118 pos = ::ftell( file ); // remember the position to return if the next line holds node ids
121 ::fseek( file, pos, SEEK_SET ); // return to the 1st line holding nodes ids
128 FT_NodeToMove nodeIds;
129 std::vector< int > ids;
131 const int nbNodes = theNodeCoords->getNumberOfTuples(); // to check validity of node IDs
133 while ( ::fgets( line, maxLineLen, file )) // read a line
135 // find node ids in the line
137 char *beg = line, *end = 0;
141 while (( id = ::strtol( beg, &end, 10 )) &&
146 throw std::invalid_argument( "Too large node ID: " + FT_Utils::toStr( id ));
150 if ( ids.size() >= 3 )
152 std::vector< int >::iterator i = ids.begin();
153 nodeIds._nodeToMove = *i;
154 nodeIds._neighborNodes.assign( ++i, ids.end() );
156 _nodes.push_back( nodeIds );
163 // -----------------------------------------------------------------
164 // try to find FT_Projector's to boundary sub-shapes by group names
165 // -----------------------------------------------------------------
167 _allProjectors = & theAllProjectorsByDim[ _shapeDim - 1 ];
169 _projectors.reserve( geomNames.size() );
170 std::vector< const FT_Projector* > projectors;
172 for ( size_t i = 0; i < geomNames.size(); ++i )
174 std::string & groupName = geomNames[i];
176 std::cout << ". Group name: " << groupName << std::endl;
179 // remove trailing white spaces
180 for ( int iC = groupName.size() - 1; iC >= 0; --iC )
182 if ( ::isspace( groupName[iC] ) )
183 groupName.resize( iC );
187 if ( groupName.empty() )
190 _groupNames.push_back( groupName ); // keep _groupNames for easier debug :)
192 // get projectors by group name
193 theXaoGroups.getProjectors( groupName, _shapeDim,
194 theAllProjectorsByDim[ _shapeDim-1 ], projectors );
197 // ------------------------------
198 // check the found FT_Projector's
199 // ------------------------------
201 if ( projectors.size() == 1 )
203 _projectors.push_back( *projectors[ 0 ]);
208 for ( size_t i = 0; i < _nodes.size(); ++i )
209 nodesBox.Add( getPoint( _nodes[i]._nodeToMove ));
211 if ( projectors.size() > 1 )
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 ]);
220 if ( _projectors.empty() )
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 ]);
228 if ( _projectors.empty() && !_nodes.empty() )
229 throw std::runtime_error("No boundary shape found for nodes in file " + theNodeFile );
233 // prepare for projection - create real projectors
234 for ( size_t i = 0; i < _projectors.size(); ++i )
235 _projectors[ i ].prepareForProjection();
239 //================================================================================
241 * \brief Project nodes to the shapes and move them to new positions
243 //================================================================================
245 void FT_NodesOnGeom::projectAndMove()
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();
258 // set nodes in the order suitable for optimal projection
261 // project and move nodes
263 std::vector< FT_NodeToMove* > notProjectedNodes;
264 size_t iP, iProjector;
268 std::cout << ".. _projectors.size() = " << _projectors.size() << std::endl;
269 std::cout << ".. _nodesOrder.size() = " << _nodesOrder.size() << std::endl;
273 // 2.1. Avec plusieurs shapes
275 if ( _projectors.size() > 1 )
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
282 for ( size_t i = 0; i < _nodesOrder.size(); ++i )
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 ))
292 moveNode( nn._nodeToMove, newXyz );
294 else // a node is not on iProjector-th shape, find the shape it is on
296 for ( iP = 1; iP < _projectors.size(); ++iP ) // check _projectors other than iProjector
298 iProjector = ( iProjector + 1 ) % _projectors.size();
299 if ( _projectors[ iProjector ].projectAndClassify( xyz, maxDist2, newXyz,
300 nn._params, nn._nearParams ))
302 moveNode( nn._nodeToMove, newXyz );
306 if ( iP == _projectors.size() )
308 notProjectedNodes.push_back( &nn );
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;
320 // 2.2. Avec une seule shape
324 for ( size_t i = 0; i < _nodesOrder.size(); ++i )
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] );
331 // maxDist2 : le quart du carré de la distance entre les deux voisins du noeud à bouger
332 double maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.;
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;
336 if ( _projectors[ 0 ].project( xyz, maxDist2, newXyz,
337 nn._params, nn._nearParams ))
338 moveNode( nn._nodeToMove, newXyz );
340 notProjectedNodes.push_back( &nn );
346 if ( !notProjectedNodes.empty() )
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
352 std::vector< FT_Projector > & projectors = *_allProjectors;
355 for ( size_t i = 0; i < notProjectedNodes.size(); ++i )
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;
365 for ( iP = 0; iP < projectors.size(); ++iP )
367 projectors[ iProjector ].prepareForProjection();
368 projectors[ iProjector ].tryWithoutPrevSolution( true );
370 if (( ok = projectors[ iProjector ].isOnShape( xyz1, tol2, nn._params, nn._nearParams )) &&
371 ( ok = projectors[ iProjector ].isOnShape( xyz2, tol2, nn._params, nn._params )))
373 if ( nn._neighborNodes.size() == 4 )
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 );
382 if ( ok && projectors[iProjector].project( xyz, maxDist2, newXyz, nn._params, nn._params ))
384 moveNode( nn._nodeToMove, newXyz );
387 iProjector = ( iProjector + 1 ) % projectors.size();
389 if ( iP == projectors.size() )
393 std::cerr << "Error: not projected node " << nn._nodeToMove << std::endl;
399 //================================================================================
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
404 //================================================================================
406 void FT_NodesOnGeom::putNodesInOrder()
408 if ( !_nodesOrder.empty() )
411 // check if any of projectors can use parameters of a previously projected node on a shape
412 // to speed up projection
414 bool isPrevSolutionUsed = false;
415 for ( size_t i = 0; i < _projectors.size() && !isPrevSolutionUsed; ++i )
416 isPrevSolutionUsed = _projectors[i].canUsePrevSolution();
418 if ( !isPrevSolutionUsed )
420 _nodesOrder.resize( _nodes.size() );
421 for ( size_t i = 0; i < _nodesOrder.size(); ++i )
422 _nodesOrder[ i ] = i;
426 // make a map to find a neighbor projected node
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;
433 int mapSize = ( _shapeDim == 1 ) ? _nodes.size() + 1 : _nodes.size() * 3;
435 neigborsMap.ReSize( mapSize );
437 std::vector< FT_NodeToMove* > linkVec, *linkVecPtr;
438 const int maxNbLinks = ( _shapeDim == 1 ) ? 2 : 6; // usual nb of links
440 for ( size_t i = 0; i < _nodes.size(); ++i )
442 FT_NodeToMove& nn = _nodes[i];
443 for ( size_t iN = 0; iN < nn._neighborNodes.size(); ++iN )
445 if ( !( linkVecPtr = neigborsMap.ChangeSeek( nn._neighborNodes[ iN ] )))
447 linkVecPtr = neigborsMap.Bound( nn._neighborNodes[ iN ], linkVec );
448 linkVecPtr->reserve( maxNbLinks );
450 linkVecPtr->push_back( & nn );
454 // fill in _nodesOrder
456 _nodesOrder.reserve( _nodes.size() );
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
462 while ( !queue.empty() )
464 FT_NodeToMove* nn = queue.front();
467 _nodesOrder.push_back( nn - & _nodes[0] );
469 // add neighbors to the queue and set their _nearParams = nn->_params
470 for ( size_t iN = 0; iN < nn->_neighborNodes.size(); ++iN )
472 std::vector< FT_NodeToMove* >& linkVec = neigborsMap( nn->_neighborNodes[ iN ]);
473 for ( size_t iL = 0; iL < linkVec.size(); ++iL )
475 FT_NodeToMove* nnn = linkVec[ iL ];
476 if ( nnn != nn && nnn->_nearParams == 0 )
478 nnn->_nearParams = nn->_params;
479 queue.push_back( nnn );
484 _nodes[0]._nearParams = 0; // reset
487 //================================================================================
489 * \brief Get node coordinates. Node IDs count from a unit
491 //================================================================================
493 gp_Pnt FT_NodesOnGeom::getPoint( const int nodeID )
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] );
500 //================================================================================
502 * \brief change node coordinates
504 //================================================================================
506 void FT_NodesOnGeom::moveNode( const int nodeID, const gp_Pnt& newXyz )
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] );