From: jfa Date: Tue, 21 Dec 2021 07:01:40 +0000 (+0300) Subject: Enable FrontTrack for CAO boundaries X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=010b8a9efdf1b9df0bf16373e7968e11a2bfa17a;p=modules%2Fsmesh.git Enable FrontTrack for CAO boundaries --- diff --git a/src/SMESH/SMESH_Homard.cxx b/src/SMESH/SMESH_Homard.cxx index f9567d554..9351216ac 100644 --- a/src/SMESH/SMESH_Homard.cxx +++ b/src/SMESH/SMESH_Homard.cxx @@ -1,4 +1,4 @@ -// HOMARD HOMARD : implementation of HOMARD idl descriptions +// SMESH HOMARD : implementation of SMESHHOMARD idl descriptions // // Copyright (C) 2011-2021 CEA/DEN, EDF R&D // @@ -18,18 +18,6 @@ // // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -// File : HOMARD_Boundary.cxx -// Author : Gerald NICOLAS, EDF -// Module : HOMARD -// -// Remarques : -// L'ordre de description des fonctions est le meme dans tous les fichiers -// HOMARD_aaaa.idl, HOMARD_aaaa.hxx, HOMARD_aaaa.cxx, HOMARD_aaaa_i.hxx, HOMARD_aaaa_i.cxx : -// 1. Les generalites : Name, Delete, DumpPython, Dump, Restore -// 2. Les caracteristiques -// 3. Le lien avec les autres structures -// -// Quand les 2 fonctions Setxxx et Getxxx sont presentes, Setxxx est decrit en premier #include "SMESH_Homard.hxx" @@ -47,6 +35,58 @@ #include #endif +//// + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boofs = boost::filesystem; + +//// + namespace SMESHHOMARDImpl { @@ -129,7 +169,7 @@ namespace SMESHHOMARDImpl std::vector coor = cas.GetBoundingBox(); os << separator() << coor.size(); - for ( int i = 0; i < coor.size(); i++ ) + for ( unsigned int i = 0; i < coor.size(); i++ ) os << separator() << coor[i]; std::list ListString = cas.GetIterations(); @@ -291,10 +331,10 @@ namespace SMESHHOMARDImpl } else { std::vector coor = boundary.GetCoords() ; - for ( int i = 0; i < coor.size(); i++ ) + for ( unsigned int i = 0; i < coor.size(); i++ ) os << separator() << coor[i]; std::vector limit = boundary.GetLimit(); - for ( int i = 0; i < limit.size(); i++ ) + for ( unsigned int i = 0; i < limit.size(); i++ ) os << separator() << limit[i]; } @@ -1084,7 +1124,7 @@ void HOMARD_Cas::SetBoundingBox( const std::vector& extremas ) { _Boite.clear(); _Boite.resize( extremas.size() ); - for ( int i = 0; i < extremas.size(); i++ ) + for ( unsigned int i = 0; i < extremas.size(); i++ ) _Boite[i] = extremas[i]; } //============================================================================= @@ -2937,4 +2977,1736 @@ int HOMARD_Iteration::GetInfoCompute() const return _MessInfo; } + +/*! + * \brief Relocate nodes to lie on geometry + * \param [in] theInputMedFile - a MED file holding a mesh including nodes that will be + * moved onto the geometry + * \param [in] theOutputMedFile - a MED file to create, that will hold a modified mesh + * \param [in] theInputNodeFiles - an array of names of files describing groups of nodes that + * will be moved onto the geometry + * \param [in] theXaoFileName - a path to a file in XAO format holding the geometry and + * the geometrical groups. + * \param [in] theIsParallel - if \c true, all processors are used to treat boundary shapes + * in parallel. + */ +void FrontTrack::track( const std::string& theInputMedFile, + const std::string& theOutputMedFile, + const std::vector< std::string > & theInputNodeFiles, + const std::string& theXaoFileName, + bool theIsParallel ) +{ + // check arguments +#ifdef _DEBUG_ + std::cout << "FrontTrack::track" << std::endl; +#endif + + if ( theInputNodeFiles.empty() ) + return; + +#ifdef _DEBUG_ + std::cout << "Input MED file: " << theInputMedFile << std::endl; +#endif + if ( !FT_Utils::fileExists( theInputMedFile )) + throw std::invalid_argument( "Input MED file does not exist: " + theInputMedFile ); + +#ifdef _DEBUG_ + std::cout << "Output MED file: " << theOutputMedFile << std::endl; +#endif + if ( !FT_Utils::canWrite( theOutputMedFile )) + throw std::invalid_argument( "Can't create the output MED file: " + theOutputMedFile ); + + std::vector< std::string > theNodeFiles ; + for ( size_t i = 0; i < theInputNodeFiles.size(); ++i ) + { +#ifdef _DEBUG_ + std::cout << "Initial input node file #"<1 + // keep only files with more than 1 line: + std::ifstream fichier(theInputNodeFiles[i].c_str()); + std::string s; + unsigned int nb_lines = 0; + while(std::getline(fichier,s)) ++nb_lines; +// std::cout << ". nb_lines: " << nb_lines << std::endl; + if ( nb_lines >= 2 ) { theNodeFiles.push_back( theInputNodeFiles[i] ); } + } +#ifdef _DEBUG_ + for ( size_t i = 0; i < theNodeFiles.size(); ++i ) + { std::cout << "Valid input node file #"< + mfMesh( MEDCoupling::MEDFileUMesh::New( theInputMedFile )); + if ( mfMesh.isNull() ) + throw std::invalid_argument( "Failed to read the input MED file: " + theInputMedFile ); + + MEDCoupling::DataArrayDouble * nodeCoords = mfMesh->getCoords(); + if ( !nodeCoords || nodeCoords->empty() ) + throw std::invalid_argument( "No nodes in the input mesh" ); + + + // read a geometry + +#ifdef _DEBUG_ + std::cout << "Lecture de la geometrie" << std::endl; +#endif + XAO::Xao xao; + if ( !xao.importXAO( theXaoFileName ) || !xao.getGeometry() ) + throw std::invalid_argument( "Failed to read the XAO input file: " + theXaoFileName ); + +#ifdef _DEBUG_ + std::cout << "Conversion en BREP" << std::endl; +#endif + XAO::BrepGeometry* xaoGeom = dynamic_cast( xao.getGeometry() ); + if ( !xaoGeom || xaoGeom->getTopoDS_Shape().IsNull() ) + throw std::invalid_argument( "Failed to get a BREP shape from the XAO input file" ); + + + // read groups of nodes and associate them with boundary shapes using names (no projection so far) + +#ifdef _DEBUG_ + std::cout << "Lecture des groupes" << std::endl; +#endif + FT_NodeGroups nodeGroups; + nodeGroups.read( theNodeFiles, &xao, nodeCoords ); +#ifdef _DEBUG_ + std::cout << "Nombre de groupes : " << nodeGroups.nbOfGroups() << std::endl; +#endif + + // project nodes to the boundary shapes and change their coordinates + +#ifdef _DEBUG_ + std::cout << "Projection des noeuds, theIsParallel=" << theIsParallel << std::endl; +#endif + OSD_Parallel::For( 0, nodeGroups.nbOfGroups(), nodeGroups, !theIsParallel ); + + // save the modified mesh + +#ifdef _DEBUG_ + std::cout << "Ecriture du maillage" << std::endl; +#endif + const int erase = 2; + mfMesh->write( theOutputMedFile, /*mode=*/erase ); + + if ( !nodeGroups.isOK() ) + throw std::runtime_error("Unable to project some nodes"); +} + + //================================================================================ + /*! + * \brief Initialize FT_Projector's with all sub-shapes of given type + * \param [in] theMainShape - the shape to explore + * \param [in] theSubType - the type of sub-shapes + * \param [out] theProjectors - the projectors + */ + //================================================================================ + + void getProjectors( const TopoDS_Shape& theMainShape, + const TopAbs_ShapeEnum theSubType, + std::vector< FT_Projector > & theProjectors ) + { + TopTools_IndexedMapOfShape subShapes; + TopExp::MapShapes( theMainShape, theSubType, subShapes ); +#ifdef _DEBUG_ + std::cout << ". Nombre de subShapes : " << subShapes.Size() << std::endl; +#endif + + theProjectors.resize( subShapes.Size() ); + for ( int i = 1; i <= subShapes.Size(); ++i ) + theProjectors[ i-1 ].setBoundaryShape( subShapes( i )); + } + +//================================================================================ +/*! + * \brief Load node groups from files + * \param [in] theNodeFiles - an array of names of files describing groups of nodes that + * will be moved onto geometry + * \param [in] theXaoGeom - the whole geometry to project on + * \param [inout] theNodeCoords - array of node coordinates + */ +//================================================================================ + +void FT_NodeGroups::read( const std::vector< std::string >& theNodeFiles, + const XAO::Xao* theXao, + MEDCoupling::DataArrayDouble* theNodeCoords ) +{ + // get projectors for all boundary sub-shapes; + // index of a projector in the vector corresponds to a XAO index of a sub-shape + XAO::BrepGeometry* xaoGeom = dynamic_cast( theXao->getGeometry() ); + getProjectors( xaoGeom->getTopoDS_Shape(), TopAbs_EDGE, _projectors[0] ); + getProjectors( xaoGeom->getTopoDS_Shape(), TopAbs_FACE, _projectors[1] ); + + _nodesOnGeom.resize( theNodeFiles.size() ); + + // read node IDs and look for projectors to boundary sub-shapes by group name + FT_Utils::XaoGroups xaoGroups( theXao ); + for ( size_t i = 0; i < theNodeFiles.size(); ++i ) + { + _nodesOnGeom[i].read( theNodeFiles[i], xaoGroups, theNodeCoords, _projectors ); + } +} + +//================================================================================ +/*! + * \brief Project and move nodes of a given group of nodes + */ +//================================================================================ + +void FT_NodeGroups::projectAndMove( const int groupIndex ) +{ + _nodesOnGeom[ groupIndex ].projectAndMove(); +} + +//================================================================================ +/*! + * \brief Return true if all nodes were successfully relocated + */ +//================================================================================ + +bool FT_NodeGroups::isOK() const +{ + for ( size_t i = 0; i < _nodesOnGeom.size(); ++i ) + if ( ! _nodesOnGeom[ i ].isOK() ) + return false; + + return true; +} + +//================================================================================ +/*! + * \brief Print some statistics on node groups + */ +//================================================================================ + +void FT_NodeGroups::dumpStat() const +{ + for ( size_t i = 0; i < _nodesOnGeom.size(); ++i ) + { + std::cout << _nodesOnGeom[i].getShapeDim() << "D " + << _nodesOnGeom[i].nbNodes() << " nodes" << std::endl; + } +} + + /*! + * \brief Close a file at destruction + */ + struct FileCloser + { + FILE * _file; + + FileCloser( FILE * file ): _file( file ) {} + ~FileCloser() { if ( _file ) ::fclose( _file ); } + }; + +//================================================================================ +/*! + * \brief Read node ids from a file and find shapes for projection + * \param [in] theNodeFile - a name of file holding IDs of nodes that + * will be moved onto geometry + * \param [in] theXaoGroups - a tool returning FT_Projector's by XAO group name + * \param [inout] theNodeCoords - array of node coordinates + * \param [in] theAllProjectorsByDim - all projectors of 2 dimensions, ordered so that + * a vector index corresponds to a XAO sub-shape ID + */ +//================================================================================ + +void FT_NodesOnGeom::read( const std::string& theNodeFile, + const FT_Utils::XaoGroups& theXaoGroups, + MEDCoupling::DataArrayDouble* theNodeCoords, + std::vector< FT_Projector > * theAllProjectorsByDim ) +{ + _nodeCoords = theNodeCoords; + + FILE * file = ::fopen( theNodeFile.c_str(), "r" ); + if ( !file ) + throw std::invalid_argument( "Can't open an input node file: " + theNodeFile ); + + FileCloser fileCloser( file ); + + // ------------------------------------- + // get shape dimension by the file name + // ------------------------------------- + + // hope the file name is something like "frnD.**" with n in (1,2) + int dimPos = theNodeFile.size() - 5; + if ( theNodeFile[ dimPos ] == '2' ) + _shapeDim = 2; + else if ( theNodeFile[ dimPos ] == '1' ) + _shapeDim = 1; + else + throw std::invalid_argument( "Can't define dimension by node file name " + theNodeFile ); +#ifdef _DEBUG_ + std::cout << ". Dimension of the file " << theNodeFile << ": " << _shapeDim << std::endl; +#endif + + // ------------------------------------- + // read geom group names; several lines + // ------------------------------------- + + std::vector< std::string > geomNames; + + const int maxLineLen = 256; + char line[ maxLineLen ]; + + long int pos = ::ftell( file ); + while ( ::fgets( line, maxLineLen, file )) // read a line + { + if ( ::feof( file )) + { + return; // no nodes in the file + } + + // check if the line describes node ids in format 3I10 (e.g. " 120 1 43\n") + size_t lineLen = strlen( line ); + if ( lineLen >= 31 && + ::isdigit( line[9] ) && + line[10] == ' ' && + ::isdigit( line[19] ) && + line[20] == ' ' && + ::isdigit( line[29] ) && + ::isspace( line[30] )) + break; + + geomNames.push_back( line + 1 ); // skip the 1st white space + + pos = ::ftell( file ); // remember the position to return if the next line holds node ids + } + + ::fseek( file, pos, SEEK_SET ); // return to the 1st line holding nodes ids + + + // -------------- + // read node ids + // -------------- + + FT_NodeToMove nodeIds; + std::vector< int > ids; + + const int nbNodes = theNodeCoords->getNumberOfTuples(); // to check validity of node IDs + + while ( ::fgets( line, maxLineLen, file )) // read a line + { + // find node ids in the line + + char *beg = line, *end = 0; + long int id; + + ids.clear(); + while (( id = ::strtol( beg, &end, 10 )) && + ( beg != end )) + { + ids.push_back( id ); + if ( id > nbNodes ) + throw std::invalid_argument( "Too large node ID: " + FT_Utils::toStr( id )); + beg = end; + } + + if ( ids.size() >= 3 ) + { + std::vector< int >::iterator i = ids.begin(); + nodeIds._nodeToMove = *i; + nodeIds._neighborNodes.assign( ++i, ids.end() ); + + _nodes.push_back( nodeIds ); + } + + if ( ::feof( file )) + break; + } + + // ----------------------------------------------------------------- + // try to find FT_Projector's to boundary sub-shapes by group names + // ----------------------------------------------------------------- + + _allProjectors = & theAllProjectorsByDim[ _shapeDim - 1 ]; + + _projectors.reserve( geomNames.size() ); + std::vector< const FT_Projector* > projectors; + + for ( size_t i = 0; i < geomNames.size(); ++i ) + { + std::string & groupName = geomNames[i]; +#ifdef _DEBUG_ + std::cout << ". Group name: " << groupName << std::endl; +#endif + + // remove trailing white spaces + for ( int iC = groupName.size() - 1; iC >= 0; --iC ) + { + if ( ::isspace( groupName[iC] ) ) + groupName.resize( iC ); + else + break; + } + if ( groupName.empty() ) + continue; + + _groupNames.push_back( groupName ); // keep _groupNames for easier debug :) + + // get projectors by group name + theXaoGroups.getProjectors( groupName, _shapeDim, + theAllProjectorsByDim[ _shapeDim-1 ], projectors ); + } + + // ------------------------------ + // check the found FT_Projector's + // ------------------------------ + + if ( projectors.size() == 1 ) + { + _projectors.push_back( *projectors[ 0 ]); + } + else + { + Bnd_Box nodesBox; + for ( size_t i = 0; i < _nodes.size(); ++i ) + nodesBox.Add( getPoint( _nodes[i]._nodeToMove )); + + if ( projectors.size() > 1 ) + { + // more than one boundary shape; + // try to filter off unnecessary projectors using a bounding box of nodes + for ( size_t i = 0; i < projectors.size(); ++i ) + if ( !nodesBox.IsOut( projectors[ i ]->getBoundingBox() )) + _projectors.push_back( *projectors[ i ]); + } + + if ( _projectors.empty() ) + { + // select projectors using a bounding box of nodes + std::vector< FT_Projector > & allProjectors = *_allProjectors; + for ( size_t i = 0; i < allProjectors.size(); ++i ) + if ( !nodesBox.IsOut( allProjectors[ i ].getBoundingBox() )) + _projectors.push_back( allProjectors[ i ]); + + if ( _projectors.empty() && !_nodes.empty() ) + throw std::runtime_error("No boundary shape found for nodes in file " + theNodeFile ); + } + } + + // prepare for projection - create real projectors + for ( size_t i = 0; i < _projectors.size(); ++i ) + _projectors[ i ].prepareForProjection(); + +} + +//================================================================================ +/*! + * \brief Project nodes to the shapes and move them to new positions + */ +//================================================================================ + +void FT_NodesOnGeom::projectAndMove() +{ + _OK = true; +// +// 1. Préalables +// + // check if all the shapes are planar + bool isAllPlanar = true; + for ( size_t i = 0; i < _projectors.size() && isAllPlanar; ++i ) + isAllPlanar = _projectors[i].isPlanarBoundary(); + if ( isAllPlanar ) + return; + + // set nodes in the order suitable for optimal projection + putNodesInOrder(); + + // project and move nodes + + std::vector< FT_NodeToMove* > notProjectedNodes; + size_t iP, iProjector; + gp_Pnt newXyz; + +#ifdef _DEBUG_ + std::cout << ".. _projectors.size() = " << _projectors.size() << std::endl; + std::cout << ".. _nodesOrder.size() = " << _nodesOrder.size() << std::endl; +#endif +// +// 2. Calculs +// 2.1. Avec plusieurs shapes +// + if ( _projectors.size() > 1 ) + { + // the nodes are to be projected onto several boundary shapes; + // in addition to the projecting, classification on a shape is necessary + // in order to find out on which of the shapes a node is to be projected + + iProjector = 0; + for ( size_t i = 0; i < _nodesOrder.size(); ++i ) + { + FT_NodeToMove& nn = _nodes[ _nodesOrder[ i ]]; + gp_Pnt xyz = getPoint( nn._nodeToMove ); + gp_Pnt xyz1 = getPoint( nn._neighborNodes[0] ); + gp_Pnt xyz2 = getPoint( nn._neighborNodes[1] ); + double maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.; + if ( _projectors[ iProjector ].projectAndClassify( xyz, maxDist2, newXyz, + nn._params, nn._nearParams )) + { + moveNode( nn._nodeToMove, newXyz ); + } + else // a node is not on iProjector-th shape, find the shape it is on + { + for ( iP = 1; iP < _projectors.size(); ++iP ) // check _projectors other than iProjector + { + iProjector = ( iProjector + 1 ) % _projectors.size(); + if ( _projectors[ iProjector ].projectAndClassify( xyz, maxDist2, newXyz, + nn._params, nn._nearParams )) + { + moveNode( nn._nodeToMove, newXyz ); + break; + } + } + if ( iP == _projectors.size() ) + { + notProjectedNodes.push_back( &nn ); + +#ifdef _DEBUG_ + std::cerr << "Warning: no shape found for node " << nn._nodeToMove << std::endl; + if ( !_groupNames.empty() ) + std::cerr << "Warning: group -- " << _groupNames[0] << std::endl; +#endif + } + } + } + } +// +// 2.2. Avec une seule shape +// + else // one shape + { + for ( size_t i = 0; i < _nodesOrder.size(); ++i ) + { + FT_NodeToMove& nn = _nodes[ _nodesOrder[ i ]]; + gp_Pnt xyz = getPoint( nn._nodeToMove ); + gp_Pnt xyz1 = getPoint( nn._neighborNodes[0] ); + gp_Pnt xyz2 = getPoint( nn._neighborNodes[1] ); + +// maxDist2 : le quart du carré de la distance entre les deux voisins du noeud à bouger + double maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.; +#ifdef _DEBUG_ + 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; +#endif + if ( _projectors[ 0 ].project( xyz, maxDist2, newXyz, + nn._params, nn._nearParams )) + moveNode( nn._nodeToMove, newXyz ); + else + notProjectedNodes.push_back( &nn ); + } + } +// +// 3. Bilan +// + if ( !notProjectedNodes.empty() ) + { + // project nodes that are not projected by any of _projectors; + // a proper projector is selected by evaluation of a distance between neighbor nodes + // and a shape + + std::vector< FT_Projector > & projectors = *_allProjectors; + + iProjector = 0; + for ( size_t i = 0; i < notProjectedNodes.size(); ++i ) + { + FT_NodeToMove& nn = *notProjectedNodes[ i ]; + gp_Pnt xyz = getPoint( nn._nodeToMove ); + gp_Pnt xyz1 = getPoint( nn._neighborNodes[0] ); + gp_Pnt xyz2 = getPoint( nn._neighborNodes[1] ); + double maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.; + double tol2 = 1e-6 * maxDist2; + + bool ok; + for ( iP = 0; iP < projectors.size(); ++iP ) + { + projectors[ iProjector ].prepareForProjection(); + projectors[ iProjector ].tryWithoutPrevSolution( true ); + + if (( ok = projectors[ iProjector ].isOnShape( xyz1, tol2, nn._params, nn._nearParams )) && + ( ok = projectors[ iProjector ].isOnShape( xyz2, tol2, nn._params, nn._params ))) + { + if ( nn._neighborNodes.size() == 4 ) + { + gp_Pnt xyz1 = getPoint( nn._neighborNodes[2] ); + gp_Pnt xyz2 = getPoint( nn._neighborNodes[3] ); + if (( ok = projectors[ iProjector ].isOnShape( xyz1, tol2, nn._params, nn._params ))) + ok = projectors[ iProjector ].isOnShape( xyz2, tol2, nn._params, nn._params ); + } + } + + if ( ok && projectors[iProjector].project( xyz, maxDist2, newXyz, nn._params, nn._params )) + { + moveNode( nn._nodeToMove, newXyz ); + break; + } + iProjector = ( iProjector + 1 ) % projectors.size(); + } + if ( iP == projectors.size() ) + { + _OK = false; + + std::cerr << "Error: not projected node " << nn._nodeToMove << std::endl; + } + } + } +} + +//================================================================================ +/*! + * \brief Put nodes in the order for optimal projection and set FT_NodeToMove::_nearParams + * to point to a FT_NodeToMove::_params of a node that will be projected earlier + */ +//================================================================================ + +void FT_NodesOnGeom::putNodesInOrder() +{ + if ( !_nodesOrder.empty() ) + return; + + // check if any of projectors can use parameters of a previously projected node on a shape + // to speed up projection + + bool isPrevSolutionUsed = false; + for ( size_t i = 0; i < _projectors.size() && !isPrevSolutionUsed; ++i ) + isPrevSolutionUsed = _projectors[i].canUsePrevSolution(); + + if ( !isPrevSolutionUsed ) + { + _nodesOrder.resize( _nodes.size() ); + for ( size_t i = 0; i < _nodesOrder.size(); ++i ) + _nodesOrder[ i ] = i; + return; + } + + // make a map to find a neighbor projected node + + // map of { FT_NodeToMove::_neighborNodes[i] } to { FT_NodeToMove* }; + // here we call FT_NodeToMove a 'link' as this data links a _neighborNodes[i] node to other nodes + typedef NCollection_DataMap< int, std::vector< FT_NodeToMove* > > TNodeIDToLinksMap; + TNodeIDToLinksMap neigborsMap; + + int mapSize = ( _shapeDim == 1 ) ? _nodes.size() + 1 : _nodes.size() * 3; + neigborsMap.Clear(); + neigborsMap.ReSize( mapSize ); + + std::vector< FT_NodeToMove* > linkVec, *linkVecPtr; + const int maxNbLinks = ( _shapeDim == 1 ) ? 2 : 6; // usual nb of links + + for ( size_t i = 0; i < _nodes.size(); ++i ) + { + FT_NodeToMove& nn = _nodes[i]; + for ( size_t iN = 0; iN < nn._neighborNodes.size(); ++iN ) + { + if ( !( linkVecPtr = neigborsMap.ChangeSeek( nn._neighborNodes[ iN ] ))) + { + linkVecPtr = neigborsMap.Bound( nn._neighborNodes[ iN ], linkVec ); + linkVecPtr->reserve( maxNbLinks ); + } + linkVecPtr->push_back( & nn ); + } + } + + // fill in _nodesOrder + + _nodesOrder.reserve( _nodes.size() ); + + std::list< FT_NodeToMove* > queue; + queue.push_back( &_nodes[0] ); + _nodes[0]._nearParams = _nodes[0]._params; // to avoid re-adding to the queue + + while ( !queue.empty() ) + { + FT_NodeToMove* nn = queue.front(); + queue.pop_front(); + + _nodesOrder.push_back( nn - & _nodes[0] ); + + // add neighbors to the queue and set their _nearParams = nn->_params + for ( size_t iN = 0; iN < nn->_neighborNodes.size(); ++iN ) + { + std::vector< FT_NodeToMove* >& linkVec = neigborsMap( nn->_neighborNodes[ iN ]); + for ( size_t iL = 0; iL < linkVec.size(); ++iL ) + { + FT_NodeToMove* nnn = linkVec[ iL ]; + if ( nnn != nn && nnn->_nearParams == 0 ) + { + nnn->_nearParams = nn->_params; + queue.push_back( nnn ); + } + } + } + } + _nodes[0]._nearParams = 0; // reset +} + +//================================================================================ +/*! + * \brief Get node coordinates. Node IDs count from a unit + */ +//================================================================================ + +gp_Pnt FT_NodesOnGeom::getPoint( const int nodeID ) +{ + const size_t dim = _nodeCoords->getNumberOfComponents(); + const double * xyz = _nodeCoords->getConstPointer() + ( dim * ( nodeID - 1 )); + return gp_Pnt( xyz[0], xyz[1], dim == 2 ? 0 : xyz[2] ); +} + +//================================================================================ +/*! + * \brief change node coordinates + */ +//================================================================================ + +void FT_NodesOnGeom::moveNode( const int nodeID, const gp_Pnt& newXyz ) +{ + const size_t dim = _nodeCoords->getNumberOfComponents(); + double z, *xyz = _nodeCoords->getPointer() + ( dim * ( nodeID - 1 )); + newXyz.Coord( xyz[0], xyz[1], dim == 2 ? z : xyz[2] ); +} + +//----------------------------------------------------------------------------- +/*! + * \brief Root class of a projector of a point to a boundary shape + */ +struct FT_RealProjector +{ + virtual ~FT_RealProjector() {} + + /*! + * \brief Project a point to a boundary shape + * \param [in] point - the point to project + * \param [out] newSolution - position on the shape (U or UV) found during the projection + * \param [in] prevSolution - position already found during the projection of a neighbor point + * \return gp_Pnt - the projection point + */ + virtual gp_Pnt project( const gp_Pnt& point, + double* newSolution, + const double* prevSolution = 0) = 0; + + /*! + * \brief Project a point to a boundary shape and check if the projection is within + * the shape boundary + * \param [in] point - the point to project + * \param [in] maxDist2 - the maximal allowed square distance between point and projection + * \param [out] projection - the projection point + * \param [out] newSolution - position on the shape (U or UV) found during the projection + * \param [in] prevSolution - position already found during the projection of a neighbor point + * \return bool - false if the projection point lies out of the shape boundary or + the distance the point and the projection is more than sqrt(maxDist2) + */ + virtual bool projectAndClassify( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution = 0) = 0; + + // return true if a previously found solution can be used to speed up the projection + + virtual bool canUsePrevSolution() const { return false; } + + + double _dist; // distance between the point being projected and its projection +}; + +namespace // actual projection algorithms +{ + const double theEPS = 1e-12; + + //================================================================================ + /*! + * \brief Projector to any curve + */ + //================================================================================ + + struct CurveProjector : public FT_RealProjector + { + BRepAdaptor_Curve _curve; + double _tol; + ShapeAnalysis_Curve _projector; + double _uRange[2]; + + //----------------------------------------------------------------------------- + CurveProjector( const TopoDS_Edge& e, const double tol ): + _curve( e ), _tol( tol ) + { + BRep_Tool::Range( e, _uRange[0], _uRange[1] ); + } + + //----------------------------------------------------------------------------- + // project a point to the curve + virtual gp_Pnt project( const gp_Pnt& P, + double* newSolution, + const double* prevSolution = 0) + { +#ifdef _DEBUG_ + std::cout << ".. project a point to the curve prevSolution = " << prevSolution << std::endl; +#endif + gp_Pnt proj; + Standard_Real param; + + if ( prevSolution ) + { + _dist = _projector.NextProject( prevSolution[0], _curve, P, _tol, proj, param ); + } + else + { + _dist = _projector.Project( _curve, P, _tol, proj, param, false ); + } +#ifdef _DEBUG_ + std::cout << ".. _dist : " << _dist << std::endl; +#endif + proj = _curve.Value( param ); + + newSolution[0] = param; + + return proj; + } + + //----------------------------------------------------------------------------- + // project a point to a curve and check if the projection is within the curve boundary + virtual bool projectAndClassify( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution = 0) + { +#ifdef _DEBUG_ + std::cout << ".. project a point to a curve and check " << std::endl; +#endif + projection = project( point, newSolution, prevSolution ); + return ( _uRange[0] < newSolution[0] && newSolution[0] < _uRange[1] && + _dist * _dist < maxDist2 ); + } + + //----------------------------------------------------------------------------- + // return true if a previously found solution can be used to speed up the projection + virtual bool canUsePrevSolution() const { return true; } + }; + + //================================================================================ + /*! + * \brief Projector to a straight curve. Don't project, classify only + */ + //================================================================================ + + struct LineProjector : public FT_RealProjector + { + gp_Pnt _p0, _p1; + + //----------------------------------------------------------------------------- + LineProjector( TopoDS_Edge e ) + { + e.Orientation( TopAbs_FORWARD ); + _p0 = BRep_Tool::Pnt( TopExp::FirstVertex( e )); + _p1 = BRep_Tool::Pnt( TopExp::LastVertex ( e )); + } + + //----------------------------------------------------------------------------- + // does nothing + virtual gp_Pnt project( const gp_Pnt& P, + double* newSolution, + const double* prevSolution = 0) + { + return P; + } + //----------------------------------------------------------------------------- + // check if a point lies within the line segment + virtual bool projectAndClassify( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution = 0) + { + gp_Vec edge( _p0, _p1 ); + gp_Vec p0p ( _p0, point ); + double u = ( edge * p0p ) / edge.SquareMagnitude(); // param [0,1] on the edge + projection = ( 1. - u ) * _p0.XYZ() + u * _p1.XYZ(); // projection of the point on the edge + if ( u < 0 || 1 < u ) + return false; + + // check distance + return point.SquareDistance( projection ) < theEPS * theEPS; + } + }; + + //================================================================================ + /*! + * \brief Projector to a circular edge + */ + //================================================================================ + + struct CircleProjector : public FT_RealProjector + { + gp_Circ _circle; + double _uRange[2]; + + //----------------------------------------------------------------------------- + CircleProjector( const gp_Circ& c, const double f, const double l ): + _circle( c ) + { + _uRange[0] = f; + _uRange[1] = l; + } + + //----------------------------------------------------------------------------- + // project a point to the circle + virtual gp_Pnt project( const gp_Pnt& P, + double* newSolution, + const double* prevSolution = 0) + { + // assume that P is already on the the plane of circle, since + // it is in the middle of two points lying on the circle + + // move P to the circle + const gp_Pnt& O = _circle.Location(); + gp_Vec radiusVec( O, P ); + double radius = radiusVec.Magnitude(); + if ( radius < std::numeric_limits::min() ) + return P; // P in on the axe + + gp_Pnt proj = O.Translated( radiusVec.Multiplied( _circle.Radius() / radius )); + + _dist = _circle.Radius() - radius; + + return proj; + } + + //----------------------------------------------------------------------------- + // project and check if a projection lies within the circular edge + virtual bool projectAndClassify( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution = 0) + { + _dist = -1; + projection = project( point, newSolution ); + if ( _dist < 0 || // ? + _dist * _dist > maxDist2 ) + return false; + + newSolution[0] = ElCLib::Parameter( _circle, projection ); + return ( _uRange[0] < newSolution[0] && newSolution[0] < _uRange[1] ); + } + }; + + //================================================================================ + /*! + * \brief Projector to any surface + */ + //================================================================================ + + struct SurfaceProjector : public FT_RealProjector + { + ShapeAnalysis_Surface _projector; + double _tol; + BRepTopAdaptor_FClass2d* _classifier; + + //----------------------------------------------------------------------------- + SurfaceProjector( const TopoDS_Face& face, const double tol, BRepTopAdaptor_FClass2d* cls ): + _projector( BRep_Tool::Surface( face )), + _tol( tol ), + _classifier( cls ) + { + } + //----------------------------------------------------------------------------- + // delete _classifier + ~SurfaceProjector() + { + delete _classifier; + } + + //----------------------------------------------------------------------------- + // project a point to a surface + virtual gp_Pnt project( const gp_Pnt& P, + double* newSolution, + const double* prevSolution = 0) + { + gp_Pnt2d uv; + + if ( prevSolution ) + { + gp_Pnt2d prevUV( prevSolution[0], prevSolution[1] ); + uv = _projector.NextValueOfUV( prevUV, P, _tol ); + } + else + { + uv = _projector.ValueOfUV( P, _tol ); + } + + uv.Coord( newSolution[0], newSolution[1] ); + + gp_Pnt proj = _projector.Value( uv ); + + _dist = _projector.Gap(); + + return proj; + } + + //----------------------------------------------------------------------------- + // project a point to a surface and check if the projection is within the surface boundary + virtual bool projectAndClassify( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution = 0) + { + projection = project( point, newSolution, prevSolution ); + return ( _dist * _dist < maxDist2 ) && classify( newSolution ); + } + + //----------------------------------------------------------------------------- + // check if the projection is within the shape boundary + bool classify( const double* newSolution ) + { + TopAbs_State state = _classifier->Perform( gp_Pnt2d( newSolution[0], newSolution[1]) ); + return ( state != TopAbs_OUT ); + } + + //----------------------------------------------------------------------------- + // return true if a previously found solution can be used to speed up the projection + virtual bool canUsePrevSolution() const { return true; } + }; + + //================================================================================ + /*! + * \brief Projector to a plane. Don't project, classify only + */ + //================================================================================ + + struct PlaneProjector : public SurfaceProjector + { + gp_Pln _plane; + bool _isRealPlane; // false means that a surface is planar but parametrization is different + + //----------------------------------------------------------------------------- + PlaneProjector( const gp_Pln& pln, + const TopoDS_Face& face, + BRepTopAdaptor_FClass2d* cls, + bool isRealPlane=true): + SurfaceProjector( face, 0, cls ), + _plane( pln ), + _isRealPlane( isRealPlane ) + {} + + //----------------------------------------------------------------------------- + // does nothing + virtual gp_Pnt project( const gp_Pnt& P, + double* newSolution, + const double* prevSolution = 0) + { + return P; + } + //----------------------------------------------------------------------------- + // check if a point lies within the boundry of the planar face + virtual bool projectAndClassify( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution = 0) + { + if ( _isRealPlane ) + { + ElSLib::PlaneParameters( _plane.Position(), point, newSolution[0], newSolution[1]); + projection = ElSLib::PlaneValue ( newSolution[0], newSolution[1], _plane.Position() ); + if ( projection.SquareDistance( point ) > theEPS * theEPS ) + return false; + + return SurfaceProjector::classify( newSolution ); + } + else + { + return SurfaceProjector::projectAndClassify( point, maxDist2, projection, + newSolution, prevSolution ); + } + } + //----------------------------------------------------------------------------- + // return true if a previously found solution can be used to speed up the projection + virtual bool canUsePrevSolution() const { return false; } + }; + + //================================================================================ + /*! + * \brief Projector to a cylinder + */ + //================================================================================ + + struct CylinderProjector : public SurfaceProjector + { + gp_Cylinder _cylinder; + + //----------------------------------------------------------------------------- + CylinderProjector( const gp_Cylinder& c, + const TopoDS_Face& face, + BRepTopAdaptor_FClass2d* cls ): + SurfaceProjector( face, 0, cls ), + _cylinder( c ) + {} + + //----------------------------------------------------------------------------- + // project a point to the cylinder + virtual gp_Pnt project( const gp_Pnt& P, + double* newSolution, + const double* prevSolution = 0) + { + // project the point P to the cylinder axis -> Pp + const gp_Pnt& O = _cylinder.Position().Location(); + const gp_Dir& axe = _cylinder.Position().Direction(); + gp_Vec trsl = gp_Vec( axe ).Multiplied( gp_Vec( O, P ).Dot( axe )); + gp_Pnt Pp = O.Translated( trsl ); + + // move Pp to the cylinder + gp_Vec radiusVec( Pp, P ); + double radius = radiusVec.Magnitude(); + if ( radius < std::numeric_limits::min() ) + return P; // P in on the axe + + gp_Pnt proj = Pp.Translated( radiusVec.Multiplied( _cylinder.Radius() / radius )); + + _dist = _cylinder.Radius() - radius; + + return proj; + } + //----------------------------------------------------------------------------- + // project a point to the cylinder and check if the projection is within the surface boundary + virtual bool projectAndClassify( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution = 0) + { + ElSLib::CylinderParameters( _cylinder.Position(), _cylinder.Radius(), point, + newSolution[0], newSolution[1]); + projection = ElSLib::CylinderValue( newSolution[0], newSolution[1], + _cylinder.Position(), _cylinder.Radius() ); + + return ( _dist * _dist < maxDist2 ) && SurfaceProjector::classify( newSolution ); + } + //----------------------------------------------------------------------------- + // return true if a previously found solution can be used to speed up the projection + virtual bool canUsePrevSolution() const { return false; } + }; + + //================================================================================ + /*! + * \brief Projector to a cone + */ + //================================================================================ + + struct ConeProjector : public SurfaceProjector + { + gp_Cone _cone; + + //----------------------------------------------------------------------------- + ConeProjector( const gp_Cone& c, + const TopoDS_Face& face, + BRepTopAdaptor_FClass2d* cls ): + SurfaceProjector( face, 0, cls ), + _cone( c ) + {} + + //----------------------------------------------------------------------------- + // project a point to the cone + virtual gp_Pnt project( const gp_Pnt& point, + double* newSolution, + const double* prevSolution = 0) + { + ElSLib::ConeParameters( _cone.Position(), _cone.RefRadius(), _cone.SemiAngle(), + point, newSolution[0], newSolution[1]); + gp_Pnt proj = ElSLib::ConeValue( newSolution[0], newSolution[1], + _cone.Position(), _cone.RefRadius(), _cone.SemiAngle() ); + _dist = point.Distance( proj ); + + return proj; + } + + //----------------------------------------------------------------------------- + // project a point to the cone and check if the projection is within the surface boundary + virtual bool projectAndClassify( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution = 0) + { + projection = project( point, newSolution, prevSolution ); + + return ( _dist * _dist < maxDist2 ) && SurfaceProjector::classify( newSolution ); + } + //----------------------------------------------------------------------------- + // return true if a previously found solution can be used to speed up the projection + virtual bool canUsePrevSolution() const { return false; } + }; + + //================================================================================ + /*! + * \brief Projector to a sphere + */ + //================================================================================ + + struct SphereProjector : public SurfaceProjector + { + gp_Sphere _sphere; + + //----------------------------------------------------------------------------- + SphereProjector( const gp_Sphere& s, + const TopoDS_Face& face, + BRepTopAdaptor_FClass2d* cls ): + SurfaceProjector( face, 0, cls ), + _sphere( s ) + {} + + //----------------------------------------------------------------------------- + // project a point to the sphere + virtual gp_Pnt project( const gp_Pnt& P, + double* newSolution, + const double* prevSolution = 0) + { + // move Pp to the Sphere + const gp_Pnt& O = _sphere.Location(); + gp_Vec radiusVec( O, P ); + double radius = radiusVec.Magnitude(); + if ( radius < std::numeric_limits::min() ) + return P; // P is on O + + gp_Pnt proj = O.Translated( radiusVec.Multiplied( _sphere.Radius() / radius )); + + _dist = _sphere.Radius() - radius; + + return proj; + } + + //----------------------------------------------------------------------------- + // project a point to the sphere and check if the projection is within the surface boundary + virtual bool projectAndClassify( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution = 0) + { + ElSLib::SphereParameters( _sphere.Position(), _sphere.Radius(), point, + newSolution[0], newSolution[1]); + projection = ElSLib::SphereValue( newSolution[0], newSolution[1], + _sphere.Position(), _sphere.Radius() ); + + return ( _dist * _dist < maxDist2 ) && SurfaceProjector::classify( newSolution ); + } + //----------------------------------------------------------------------------- + // return true if a previously found solution can be used to speed up the projection + virtual bool canUsePrevSolution() const { return false; } + }; + + //================================================================================ + /*! + * \brief Projector to a torus + */ + //================================================================================ + + struct TorusProjector : public SurfaceProjector + { + gp_Torus _torus; + + //----------------------------------------------------------------------------- + TorusProjector( const gp_Torus& t, + const TopoDS_Face& face, + BRepTopAdaptor_FClass2d* cls ): + SurfaceProjector( face, 0, cls ), + _torus( t ) + {} + + //----------------------------------------------------------------------------- + // project a point to the torus + virtual gp_Pnt project( const gp_Pnt& point, + double* newSolution, + const double* prevSolution = 0) + { + ElSLib::TorusParameters( _torus.Position(), _torus.MajorRadius(), _torus.MinorRadius(), + point, newSolution[0], newSolution[1]); + gp_Pnt proj = ElSLib::TorusValue( newSolution[0], newSolution[1], + _torus.Position(), _torus.MajorRadius(), _torus.MinorRadius() ); + _dist = point.Distance( proj ); + + return proj; + } + + //----------------------------------------------------------------------------- + // project a point to the torus and check if the projection is within the surface boundary + virtual bool projectAndClassify( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution = 0) + { + projection = project( point, newSolution, prevSolution ); + + return ( _dist * _dist < maxDist2 ) && SurfaceProjector::classify( newSolution ); + } + //----------------------------------------------------------------------------- + // return true if a previously found solution can be used to speed up the projection + virtual bool canUsePrevSolution() const { return false; } + }; + + //================================================================================ + /*! + * \brief Check if a curve can be considered straight + */ + //================================================================================ + + bool isStraight( const GeomAdaptor_Curve& curve, const double tol ) + { + // rough check: evaluate how far from a straight line connecting the curve ends + // stand several internal points of the curve + + const double f = curve.FirstParameter(); + const double l = curve.LastParameter(); + const gp_Pnt pf = curve.Value( f ); + const gp_Pnt pl = curve.Value( l ); + const gp_Vec lineVec( pf, pl ); + const double lineLen2 = lineVec.SquareMagnitude(); + if ( lineLen2 < std::numeric_limits< double >::min() ) + return false; // E seems closed + + const double nbSamples = 7; + for ( int i = 0; i < nbSamples; ++i ) + { + const double r = ( i + 1 ) / nbSamples; + const gp_Pnt pi = curve.Value( f * r + l * ( 1 - r )); + const gp_Vec vi( pf, pi ); + const double h2 = lineVec.Crossed( vi ).SquareMagnitude() / lineLen2; + if ( h2 > tol * tol ) + return false; + } + + // thorough check + GCPnts_UniformDeflection divider( curve, tol ); + return ( divider.IsDone() && divider.NbPoints() < 3 ); + } +} + +//================================================================================ +/*! + * \brief Initialize with a boundary shape + */ +//================================================================================ + +FT_Projector::FT_Projector(const TopoDS_Shape& shape) +{ + _realProjector = 0; + setBoundaryShape( shape ); + _tryWOPrevSolution = false; +} + +//================================================================================ +/*! + * \brief Copy another projector + */ +//================================================================================ + +FT_Projector::FT_Projector(const FT_Projector& other) +{ + _realProjector = 0; + _shape = other._shape; + _bndBox = other._bndBox; + _tryWOPrevSolution = false; +} + +//================================================================================ +/*! + * \brief Destructor. Delete _realProjector + */ +//================================================================================ + +FT_Projector::~FT_Projector() +{ + delete _realProjector; +} + +//================================================================================ +/*! + * \brief Initialize with a boundary shape. Compute the bounding box + */ +//================================================================================ + +void FT_Projector::setBoundaryShape(const TopoDS_Shape& shape) +{ + delete _realProjector; _realProjector = 0; + _shape = shape; + if ( shape.IsNull() ) + return; + + BRepBndLib::Add( shape, _bndBox ); + _bndBox.Enlarge( 1e-5 * sqrt( _bndBox.SquareExtent() )); +} + +//================================================================================ +/*! + * \brief Create a real projector + */ +//================================================================================ + +void FT_Projector::prepareForProjection() +{ + if ( _shape.IsNull() || _realProjector ) + return; + + if ( _shape.ShapeType() == TopAbs_EDGE ) + { + const TopoDS_Edge& edge = TopoDS::Edge( _shape ); + + double tol = 1e-6 * sqrt( _bndBox.SquareExtent() ); + + double f,l; + Handle(Geom_Curve) curve = BRep_Tool::Curve( edge, f,l ); + if ( curve.IsNull() ) + return; // degenerated edge + + GeomAdaptor_Curve acurve( curve, f, l ); + switch ( acurve.GetType() ) + { + case GeomAbs_Line: + _realProjector = new LineProjector( edge ); + break; + case GeomAbs_Circle: + _realProjector = new CircleProjector( acurve.Circle(), f, l ); + break; + case GeomAbs_BezierCurve: + case GeomAbs_BSplineCurve: + case GeomAbs_OffsetCurve: + case GeomAbs_OtherCurve: + if ( isStraight( acurve, tol )) + { + _realProjector = new LineProjector( edge ); + break; + } + case GeomAbs_Ellipse: + case GeomAbs_Hyperbola: + case GeomAbs_Parabola: + _realProjector = new CurveProjector( edge, tol ); + } + } + else if ( _shape.ShapeType() == TopAbs_FACE ) + { + TopoDS_Face face = TopoDS::Face( _shape ); + + Handle(Geom_Surface) surface = BRep_Tool::Surface( face ); + if ( surface.IsNull() ) + return; + + GeomAdaptor_Surface asurface( surface ); + Standard_Real tol = BRep_Tool::Tolerance( face ); + Standard_Real toluv = Min( asurface.UResolution( tol ), asurface.VResolution( tol )); + BRepTopAdaptor_FClass2d* classifier = new BRepTopAdaptor_FClass2d( face, toluv ); + + switch ( asurface.GetType() ) + { + case GeomAbs_Plane: + _realProjector = new PlaneProjector( asurface.Plane(), face, classifier ); + break; + case GeomAbs_Cylinder: + _realProjector = new CylinderProjector( asurface.Cylinder(), face, classifier ); + break; + case GeomAbs_Sphere: + _realProjector = new SphereProjector( asurface.Sphere(), face, classifier ); + break; + case GeomAbs_Cone: + _realProjector = new ConeProjector( asurface.Cone(), face, classifier ); + break; + case GeomAbs_Torus: + _realProjector = new TorusProjector( asurface.Torus(), face, classifier ); + break; + case GeomAbs_BezierSurface: + case GeomAbs_BSplineSurface: + case GeomAbs_SurfaceOfRevolution: + case GeomAbs_SurfaceOfExtrusion: + case GeomAbs_OffsetSurface: + case GeomAbs_OtherSurface: + GeomLib_IsPlanarSurface isPlaneCheck( surface, tol ); + if ( isPlaneCheck.IsPlanar() ) + { + _realProjector = new PlaneProjector( isPlaneCheck.Plan(), face, classifier, + /*isRealPlane=*/false); + } + else + { + _realProjector = new SurfaceProjector( face, tol, classifier ); + } + break; + } + + if ( !_realProjector ) + delete classifier; + } +} + +//================================================================================ +/*! + * \brief Return true if projection is not needed + */ +//================================================================================ + +bool FT_Projector::isPlanarBoundary() const +{ + return ( dynamic_cast< LineProjector* >( _realProjector ) || + dynamic_cast< PlaneProjector* >( _realProjector ) ); +} + +//================================================================================ +/*! + * \brief Check if a point lies on the boundary shape + * \param [in] point - the point to check + * \param [in] tol2 - a square tolerance allowing to decide whether a point is on the shape + * \param [in] newSolution - position on the shape (U or UV) of the point found + * during projecting + * \param [in] prevSolution - position on the shape (U or UV) of a neighbor point + * \return bool - \c true if the point lies on the boundary shape + * + * This method is used to select a shape by checking if all neighbor nodes of a node to move + * lie on a shape. + */ +//================================================================================ + +bool FT_Projector::isOnShape( const gp_Pnt& point, + const double tol2, + double* newSolution, + const double* prevSolution) +{ + if ( _bndBox.IsOut( point ) || !_realProjector ) + return false; + + gp_Pnt proj; + if ( isPlanarBoundary() ) + return projectAndClassify( point, tol2, proj, newSolution, prevSolution ); + + return project( point, tol2, proj, newSolution, prevSolution ); +} + +//================================================================================ +/*! + * \brief Project a point to the boundary shape + * \param [in] point - the point to project + * \param [in] maxDist2 - the maximal square distance between the point and the projection + * \param [out] projection - the projection + * \param [out] newSolution - position on the shape (U or UV) of the point found + * during projecting + * \param [in] prevSolution - already found position on the shape (U or UV) of a neighbor point + * \return bool - false if the distance between the point and the projection + * is more than sqrt(maxDist2) + * + * This method is used to project a node in the case where only one shape is found by name + */ +//================================================================================ + +bool FT_Projector::project( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution) +{ + if ( !_realProjector ) + return false; + + _realProjector->_dist = 1e100; + projection = _realProjector->project( point, newSolution, prevSolution ); + + bool ok = ( _realProjector->_dist * _realProjector->_dist < maxDist2 ); + if ( !ok && _tryWOPrevSolution && prevSolution ) + { + projection = _realProjector->project( point, newSolution ); + ok = ( _realProjector->_dist * _realProjector->_dist < maxDist2 ); + } + return ok; +} + +//================================================================================ +/*! + * \brief Project a point to the boundary shape and check if the projection lies within + * the shape boundary + * \param [in] point - the point to project + * \param [in] maxDist2 - the maximal square distance between the point and the projection + * \param [out] projection - the projection + * \param [out] newSolution - position on the shape (U or UV) of the point found + * during projecting + * \param [in] prevSolution - already found position on the shape (U or UV) of a neighbor point + * \return bool - false if the projection point lies out of the shape boundary or + * the distance between the point and the projection is more than sqrt(maxDist2) + * + * This method is used to project a node in the case where several shapes are selected for + * projection of a node group + */ +//================================================================================ + +bool FT_Projector::projectAndClassify( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution) +{ + if ( _bndBox.IsOut( point ) || !_realProjector ) + return false; + + bool ok = _realProjector->projectAndClassify( point, maxDist2, projection, + newSolution, prevSolution ); + if ( !ok && _tryWOPrevSolution && prevSolution ) + ok = _realProjector->projectAndClassify( point, maxDist2, projection, newSolution ); + + return ok; +} + +//================================================================================ +/*! + * \brief Return true if a previously found solution can be used to speed up the projection + */ +//================================================================================ + +bool FT_Projector::canUsePrevSolution() const +{ + return ( _realProjector && _realProjector->canUsePrevSolution() ); +} + +//================================================================================ +/* + * \brief Check if a file exists + */ +//================================================================================ + +bool FT_Utils::fileExists( const std::string& path ) +{ + if ( path.empty() ) + return false; + + boost::system::error_code err; + bool res = boofs::exists( path, err ); + + return err ? false : res; +} + +//================================================================================ +/*! + * \brief Check if a file can be created/overwritten + */ +//================================================================================ + +bool FT_Utils::canWrite( const std::string& path ) +{ + if ( path.empty() ) + return false; + + bool can = false; +#ifdef WIN32 + + HANDLE file = CreateFile( path.c_str(), // name of the write + GENERIC_WRITE, // open for writing + 0, // do not share + NULL, // default security + OPEN_ALWAYS, // CREATE NEW or OPEN EXISTING + FILE_ATTRIBUTE_NORMAL, // normal file + NULL); // no attr. template + can = ( file != INVALID_HANDLE_VALUE ); + CloseHandle( file ); + +#else + + int file = ::open( path.c_str(), + O_WRONLY | O_CREAT, + S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH ); // rw-r--r-- + can = ( file >= 0 ); + +#endif + + return can; +} + +//================================================================================ +/*! + * \brief Make a map of XAO groups + */ +//================================================================================ + +FT_Utils::XaoGroups::XaoGroups( const XAO::Xao* theXao ) +{ + XAO::Xao* xao = const_cast< XAO::Xao* >( theXao ); + + for ( int iG = 0; iG < theXao->countGroups(); ++iG ) + { + XAO::Group* group = xao->getGroup( iG ); + + if ( group->getDimension() == 1 ) + + _xaoGroups[ 0 ].insert( std::make_pair( group->getName(), group )); + + else if ( group->getDimension() == 2 ) + + _xaoGroups[ 1 ].insert( std::make_pair( group->getName(), group )); + } +} + +//================================================================================ +/*! + * \brief Return FT_Projector's by a group name + * \param [in] groupName - the group name + * \param [in] dim - the group dimension + * \param [in] allProjectors - the projector of all shapes of \a dim dimension + * \param [out] groupProjectors - projectors to shapes of the group + * \return int - number of found shapes + */ +//================================================================================ + +int FT_Utils::XaoGroups::getProjectors( const std::string& groupName, + const int dim, + const std::vector< FT_Projector > & allProjectors, + std::vector< const FT_Projector* > & groupProjectors) const +{ + // get namesake groups + + const TGroupByNameMap* groupMap = 0; + if ( dim == 1 ) + groupMap = &_xaoGroups[ 0 ]; + else if ( dim == 2 ) + groupMap = &_xaoGroups[ 1 ]; + else + return 0; + + TGroupByNameMap::const_iterator name2gr = groupMap->find( groupName ); + if ( name2gr == groupMap->end() ) + return 0; + + std::vector< XAO::Group* > groups; + groups.push_back( name2gr->second ); + + for ( ++name2gr; name2gr != groupMap->end(); ++name2gr ) + { + if ( name2gr->second->getName() == groupName ) + groups.push_back( name2gr->second ); + else + break; + } + + // get projectors + + int nbFound = 0; + for ( size_t i = 0; i < groups.size(); ++i ) + { + // IDs in XAO correspond to indices of allProjectors + std::set::iterator id = groups[i]->begin(), end = groups[i]->end(); + for ( ; id != end; ++id, ++nbFound ) + if ( *id < (int) allProjectors.size() ) + groupProjectors.push_back ( & allProjectors[ *id ]); + } + + return nbFound; +} + } // namespace SMESHHOMARDImpl /end/ diff --git a/src/SMESH/SMESH_Homard.hxx b/src/SMESH/SMESH_Homard.hxx index 45be48424..416fe3937 100644 --- a/src/SMESH/SMESH_Homard.hxx +++ b/src/SMESH/SMESH_Homard.hxx @@ -53,14 +53,29 @@ #include #include #include +#include #include #include +#include +#include +#include +#include + #if defined WIN32 #pragma warning ( disable: 4251 ) #endif +namespace MEDCoupling { + class DataArrayDouble; +} +namespace XAO { + class Xao; + class Group; + class BrepGeometry; +} + namespace SMESHHOMARDImpl { @@ -472,6 +487,240 @@ private: std::list _ListFieldInterpTSR; }; +// HOMARD/FrontTrack + +class FrontTrack +{ +public: + + /*! + * \brief Relocate nodes to lie on geometry + * \param [in] theInputMedFile - a MED file holding a mesh including nodes that will be + * moved onto the geometry + * \param [in] theOutputMedFile - a MED file to create, that will hold a modified mesh + * \param [in] theInputNodeFiles - an array of names of files describing groups of nodes that + * will be moved onto the geometry + * \param [in] theXaoFileName - a path to a file in XAO format holding the geometry and + * the geometrical groups. + * \param [in] theIsParallel - if \c true, all processors are used to treat boundary shapes + * in parallel. + */ + void track( const std::string& theInputMedFile, + const std::string& theOutputMedFile, + const std::vector< std::string > & theInputNodeFiles, + const std::string& theXaoFileName, + bool theIsParallel=true); + +}; + + +struct FT_RealProjector; + +/*! + * \brief Projector of a point to a boundary shape. Wrapper of a real projection algo + */ +class FT_Projector +{ +public: + + FT_Projector(const TopoDS_Shape& shape = TopoDS_Shape()); + FT_Projector(const FT_Projector& other); + ~FT_Projector(); + + // initialize with a boundary shape, compute the bounding box + void setBoundaryShape(const TopoDS_Shape& shape); + + // return the boundary shape + const TopoDS_Shape& getShape() const { return _shape; } + + // return the bounding box + const Bnd_Box getBoundingBox() const { return _bndBox; } + + + // create a real projector + void prepareForProjection(); + + // return true if a previously found solution can be used to speed up the projection + bool canUsePrevSolution() const; + + // return true if projection is not needed + bool isPlanarBoundary() const; + + + // switch a mode of usage of prevSolution. + // If projection fails, to try to project without usage of prevSolution. + // By default this mode is off + void tryWithoutPrevSolution( bool toTry ) { _tryWOPrevSolution = toTry; } + + // project a point to the boundary shape + bool project( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution = 0); + + // project a point to the boundary shape and check if the projection is within the shape boundary + bool projectAndClassify( const gp_Pnt& point, + const double maxDist2, + gp_Pnt& projection, + double* newSolution, + const double* prevSolution = 0); + + // check if a point lies on the boundary shape + bool isOnShape( const gp_Pnt& point, + const double tol2, + double* newSolution, + const double* prevSolution = 0); + +private: + + FT_RealProjector* _realProjector; + Bnd_Box _bndBox; + TopoDS_Shape _shape; + bool _tryWOPrevSolution; +}; + +namespace FT_Utils +{ + // Check if a file exists + bool fileExists( const std::string& path ); + + // Check if a file can be created/overwritten + bool canWrite( const std::string& path ); + + // Transform anything printable to a string + template< typename T> std::string toStr( const T& t ) + { + std::ostringstream s; + s << t; + return s.str(); + } + + //-------------------------------------------------------------------------------------------- + /*! + * \brief Return projectors by group name + */ + struct XaoGroups + { + XaoGroups( const XAO::Xao* xao ); + + int getProjectors( const std::string& groupName, + const int dim, + const std::vector< FT_Projector > & allProjectors, + std::vector< const FT_Projector* > & groupProjectors ) const; + private: + + typedef std::multimap< std::string, XAO::Group* > TGroupByNameMap; + TGroupByNameMap _xaoGroups[ 2 ]; // by dim + }; +} // namespace FT_Utils + +/*! + * \brief Node group and geometry to project onto + */ +class FT_NodesOnGeom +{ +public: + + // read node IDs form a file and try to find a boundary sub-shape by name + void read( const std::string& nodesFile, + const FT_Utils::XaoGroups& xaoGroups, + MEDCoupling::DataArrayDouble* nodeCoords, + std::vector< FT_Projector > * allProjectorsByDim); + + // chose boundary shapes by evaluating distance between nodes and shapes + //void choseShape( const std::vector< FT_Utils::ShapeAndBndBox >& shapeAndBoxList ); + + // project nodes to the shapes and move them to new positions + void projectAndMove(); + + // return true if all nodes were successfully relocated + bool isOK() const { return _OK; } + + // return dimension of boundary shapes + int getShapeDim() const { return _shapeDim; } + + // return nb of nodes to move + int nbNodes() const { return _nodes.size(); } + + +private: + + // put nodes in the order for optimal projection + void putNodesInOrder(); + + // get node coordinates + gp_Pnt getPoint( const int nodeID ); + + // change node coordinates + void moveNode( const int nodeID, const gp_Pnt& xyz ); + + + // Ids of a node to move and its 2 or 4 neighbors + struct FT_NodeToMove + { + int _nodeToMove; + std::vector< int > _neighborNodes; + + double _params[2]; // parameters on shape (U or UV) found by projection + double *_nearParams; // _params of a neighbor already projected node + + FT_NodeToMove(): _nearParams(0) {} + }; + + std::vector< std::string > _groupNames; + int _shapeDim; // dimension of boundary shapes + std::vector< FT_NodeToMove > _nodes; // ids of nodes to move and their neighbors + std::vector< FT_Projector > _projectors; // FT_Projector's initialized with boundary shapes + std::vector< FT_Projector > * _allProjectors; // FT_Projector's for all shapes of _shapeDim + MEDCoupling::DataArrayDouble* _nodeCoords; + bool _OK; // projecting is successful + + // map of { FT_NodeToMove::_neighborNodes[i] } to { FT_NodeToMove* } + // this map is used to find neighbor nodes + typedef NCollection_DataMap< int, std::vector< FT_NodeToMove* > > TNodeIDToLinksMap; + TNodeIDToLinksMap _neigborsMap; + std::vector _nodesOrder; + +}; + +/*! + * \brief Container of node groups. + */ +class FT_NodeGroups +{ +public: + + // Load node groups from files + void read( const std::vector< std::string >& nodeFiles, + const XAO::Xao* xaoGeom, + MEDCoupling::DataArrayDouble* nodeCoords ); + + // return number of groups of nodes to move + int nbOfGroups() const { return _nodesOnGeom.size(); } + + // Move nodes of a group in parallel mode + void operator() ( const int groupIndex ) const + { + const_cast< FT_NodeGroups* >( this )->projectAndMove( groupIndex ); + } + + // Project and move nodes of a given group of nodes + void projectAndMove( const int groupIndex ); + + // return true if all nodes were successfully relocated + bool isOK() const; + + // print some statistics on node groups + void dumpStat() const; + +private: + + std::vector< FT_NodesOnGeom > _nodesOnGeom; + std::vector< FT_Projector > _projectors[2]; // curves and surfaces separately + +}; + }; // namespace SMESHHOMARDImpl #endif diff --git a/src/SMESHGUI/SMESHGUI_HomardAdaptDlg.cxx b/src/SMESHGUI/SMESHGUI_HomardAdaptDlg.cxx index 40fe30e26..e1507e9db 100644 --- a/src/SMESHGUI/SMESHGUI_HomardAdaptDlg.cxx +++ b/src/SMESHGUI/SMESHGUI_HomardAdaptDlg.cxx @@ -472,6 +472,7 @@ bool SMESHGUI_HomardAdaptDlg::PushOnApply() // Compute and publish bool isSuccess = true; try { + SUIT_OverrideCursor aWaitCursor; isSuccess = myHomardGen->Compute() == 0; } catch( SALOME::SALOME_Exception& S_ex ) { @@ -540,6 +541,7 @@ void SMESHGUI_HomardAdaptDlg::selectionChanged() if (!myArgs->myInBrowserRadio->isChecked()) return; + //SUIT_OverrideCursor aWaitCursor; LightApp_SelectionMgr *selMgr = SMESHGUI::selectionMgr(); // get selected mesh @@ -580,6 +582,10 @@ void SMESHGUI_HomardAdaptDlg::selectionChanged() myArgs->mySelectOutMedFileLineEdit->setText(aFileInfo.absoluteFilePath()); } //} + + // Check data + if (!aMeshName.isEmpty()) + CheckCase(false); } void SMESHGUI_HomardAdaptDlg::SetFileName() @@ -587,6 +593,7 @@ void SMESHGUI_HomardAdaptDlg::SetFileName() // Input med file QString fileName0 = myArgs->mySelectInMedFileLineEdit->text().trimmed(); QString fileName = SMESH_HOMARD_QT_COMMUN::PushNomFichier(false, QString("med")); + //SUIT_OverrideCursor aWaitCursor; if (fileName.isEmpty()) { fileName = fileName0; if (fileName.isEmpty()) return; diff --git a/src/SMESH_I/SMESH_Homard_i.cxx b/src/SMESH_I/SMESH_Homard_i.cxx index 0c80e0f0f..19bfe7a68 100644 --- a/src/SMESH_I/SMESH_Homard_i.cxx +++ b/src/SMESH_I/SMESH_Homard_i.cxx @@ -25,9 +25,6 @@ #include "SMESH_File.hxx" -// TODO? -//#include "FrontTrack.hxx" - #include "utilities.h" #include "Basics_Utils.hxx" #include "Basics_DirUtils.hxx" @@ -1008,6 +1005,7 @@ CORBA::Long HOMARD_Gen_i::DeleteIteration(int numIter) MESSAGE ("DeleteIteration : numIter = " << numIter); if (numIter == 0) { + if (!CORBA::is_nil(myIteration1)) DeleteIteration(1); myIteration0 = SMESHHOMARD::HOMARD_Iteration::_nil(); } else { @@ -2070,16 +2068,18 @@ CORBA::Long HOMARD_Gen_i::ComputeCAO(SMESHHOMARD::HOMARD_Cas_var myCase, MESSAGE (". theXaoFileName = " << theXaoFileName); // B.5. Parallélisme - //bool theIsParallel = false; + bool theIsParallel = false; // C. Lancement des projections MESSAGE (". Lancement des projections"); // TODO? - //FrontTrack* myFrontTrack = new FrontTrack(); - //myFrontTrack->track(theInputMedFile, theOutputMedFile, theInputNodeFiles, theXaoFileName, theIsParallel); + SMESHHOMARDImpl::FrontTrack* myFrontTrack = new SMESHHOMARDImpl::FrontTrack(); + myFrontTrack->track(theInputMedFile, theOutputMedFile, + theInputNodeFiles, theXaoFileName, theIsParallel); // D. Transfert des coordonnées modifiées dans le fichier historique de HOMARD - // On lance une exécution spéciale de HOMARD en attendant de savoir le faire avec MEDCoupling + // On lance une exécution spéciale de HOMARD en attendant + // de savoir le faire avec MEDCoupling MESSAGE (". Transfert des coordonnées"); codret = ComputeCAObis(myIteration);