#include <direct.h>
#endif
-////
-
-#include <MCAuto.hxx>
-#include <MEDCouplingMemArray.hxx>
-#include <MEDFileMesh.hxx>
-
-#include <XAO_Xao.hxx>
-#include <XAO_BrepGeometry.hxx>
-#include <XAO_Group.hxx>
-
-#include <stdexcept>
-#include <cstdio>
-#include <cstdlib>
-#include <list>
-#include <limits>
-
-#include <fcntl.h>
-#include <boost/filesystem.hpp>
-
-#include <OSD_Parallel.hxx>
-#include <BRepAdaptor_Curve.hxx>
-#include <BRepBndLib.hxx>
-#include <BRepTopAdaptor_FClass2d.hxx>
-#include <BRep_Tool.hxx>
-#include <Bnd_Box.hxx>
-#include <TopExp.hxx>
-#include <TopExp_Explorer.hxx>
-#include <TopTools_IndexedMapOfShape.hxx>
-#include <TopoDS.hxx>
-#include <TopoDS_Edge.hxx>
-#include <TopoDS_Face.hxx>
-#include <TopoDS_Vertex.hxx>
-#include <TopoDS_Shape.hxx>
-#include <ElCLib.hxx>
-#include <ElSLib.hxx>
-#include <GCPnts_UniformDeflection.hxx>
-#include <GeomAdaptor_Curve.hxx>
-#include <GeomLib_IsPlanarSurface.hxx>
-#include <ShapeAnalysis_Curve.hxx>
-#include <ShapeAnalysis_Surface.hxx>
-#include <gp_Circ.hxx>
-#include <gp_Cylinder.hxx>
-#include <gp_Dir.hxx>
-#include <gp_Pln.hxx>
-#include <gp_Pnt.hxx>
-#include <gp_Sphere.hxx>
-#include <gp_Vec.hxx>
-
-namespace boofs = boost::filesystem;
-
-////
-
namespace SMESHHOMARDImpl
{
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 #"<<i<<": " << theInputNodeFiles[i] << std::endl;
-#endif
- if ( !FT_Utils::fileExists( theInputNodeFiles[i] ))
- throw std::invalid_argument( "Input node file does not exist: " + theInputNodeFiles[i] );
- // the name of the groupe on line #1, then the numbers of nodes on line #>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 #"<<i<<": " << theNodeFiles[i] << std::endl; }
-#endif
-
-#ifdef _DEBUG_
- std::cout << "XAO file: " << theXaoFileName << std::endl;
-#endif
- if ( !FT_Utils::fileExists( theXaoFileName ))
- throw std::invalid_argument( "Input XAO file does not exist: " + theXaoFileName );
-
- // read a mesh
-
-#ifdef _DEBUG_
- std::cout << "Lecture du maillage" << std::endl;
-#endif
- MEDCoupling::MCAuto< MEDCoupling::MEDFileUMesh >
- 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::BrepGeometry*>( 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<XAO::BrepGeometry*>( 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<double>::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<double>::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<double>::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<int>::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/
