From: Gerald NICOLAS Date: Wed, 31 May 2017 13:30:58 +0000 (+0200) Subject: Fronttrack X-Git-Tag: V8_4_0a2~2^2~1 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=9c4b695101f6dfd8b305c42f37449b25688bf5e1;p=modules%2Fhomard.git Fronttrack Introduction des programmes du suivi des frontières courbes. Travaux réalisés par OCC mai 2017 --- diff --git a/CMakeLists.txt b/CMakeLists.txt index 3d6a311e..15ea411e 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,6 +22,9 @@ PROJECT(SalomeHOMARD C CXX) # Ensure a proper linker behavior: CMAKE_POLICY(SET CMP0003 NEW) +IF(WIN32) + CMAKE_POLICY(SET CMP0020 OLD) # disable automatic linking to qtmain.lib +ENDIF(WIN32) # Versioning # =========== @@ -92,7 +95,7 @@ OPTION(SALOME_BUILD_TESTS "Build SALOME tests" ON) FIND_PACKAGE(SalomePythonInterp REQUIRED) FIND_PACKAGE(SalomePythonLibs REQUIRED) -# Other KERNEL optionals: +# Other KERNEL optionals: IF(SALOME_BUILD_TESTS) ENABLE_TESTING() FIND_PACKAGE(SalomeCppUnit) @@ -104,6 +107,22 @@ IF(SALOME_BUILD_DOC) ENDIF() +# Advanced options: +OPTION(FRONTTRACK_USE_TBB "Enable parallel computation using TBB" OFF) +MARK_AS_ADVANCED(FRONTTRACK_USE_TBB) + +# Prerequisites +# ============= + +# MEDCoupling +SET(MEDCOUPLING_ROOT_DIR $ENV{MEDCOUPLING_ROOT_DIR} CACHE PATH "Path to the MEDCoupling tool") +IF(EXISTS ${MEDCOUPLING_ROOT_DIR}) + LIST(APPEND CMAKE_MODULE_PATH "${MEDCOUPLING_ROOT_DIR}/cmake_files") + FIND_PACKAGE(SalomeMEDCoupling REQUIRED) # will reload HDF5, MEDFile, XDR, etc ... +ELSE(EXISTS ${MEDCOUPLING_ROOT_DIR}) + MESSAGE(FATAL_ERROR "We absolutely need the MEDCoupling tool, please define MEDCOUPLING_ROOT_DIR !") +ENDIF(EXISTS ${MEDCOUPLING_ROOT_DIR}) + ## ## From KERNEL: ## @@ -146,6 +165,8 @@ SET(GEOM_ROOT_DIR $ENV{GEOM_ROOT_DIR} CACHE PATH "Path to the Salome GEOM") IF(EXISTS ${GEOM_ROOT_DIR}) LIST(APPEND CMAKE_MODULE_PATH "${GEOM_ROOT_DIR}/adm_local/cmake_files") FIND_PACKAGE(SalomeGEOM REQUIRED) + ADD_DEFINITIONS(${GEOM_DEFINITIONS}) + INCLUDE_DIRECTORIES(${GEOM_INCLUDE_DIRS}) ELSE(EXISTS ${GEOM_ROOT_DIR}) MESSAGE(FATAL_ERROR "We absolutely need a Salome GEOM, please define GEOM_ROOT_DIR") ENDIF(EXISTS ${GEOM_ROOT_DIR}) @@ -170,6 +191,23 @@ ENDIF(EXISTS ${SMESH_ROOT_DIR}) # MedFile FIND_PACKAGE(SalomeMEDFile REQUIRED) +# OpenCASCADE +FIND_PACKAGE(SalomeCAS REQUIRED) + +# TBB +IF(FRONTTRACK_USE_TBB) + FIND_PACKAGE(SalomeTBB) + SALOME_LOG_OPTIONAL_PACKAGE(TBB FRONTTRACK_USE_TBB) + ADD_DEFINITIONS(-DHAVE_TBB) +ENDIF(FRONTTRACK_USE_TBB) + +# Python +FIND_PACKAGE(SalomePythonInterp REQUIRED) +# SWIG +FIND_PACKAGE(SalomeSWIG REQUIRED) +# Boost +FIND_PACKAGE(SalomeBoost REQUIRED) + # Detection summary: SALOME_PACKAGE_REPORT_AND_CHECK() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1226a565..92f55fa5 100755 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -25,6 +25,8 @@ SET(SUBDIRS_COMMON HOMARD_I HOMARDGUI HOMARD_SWIG + FrontTrack + FrontTrack_SWIG tests ) diff --git a/src/FrontTrack/CMakeLists.txt b/src/FrontTrack/CMakeLists.txt new file mode 100755 index 00000000..8f496978 --- /dev/null +++ b/src/FrontTrack/CMakeLists.txt @@ -0,0 +1,61 @@ +# --- options --- +# additional include directories +INCLUDE_DIRECTORIES( + ${CAS_INCLUDE_DIRS} + ${Boost_INCLUDE_DIRS} + ${GEOM_INCLUDE_DIRS} + ${MEDFILE_INCLUDE_DIRS} + ${MEDCOUPLING_INCLUDE_DIRS} + ${TBB_INCLUDE_DIRS} +) + +# additional preprocessor / compiler flags +ADD_DEFINITIONS( + ${CAS_DEFINITIONS} + ${BOOST_DEFINITIONS} +) + +IF(FRONTTRACK_USE_TBB) + SET(TBB_LIBS ${TBB_LIBRARIES}) +ENDIF(FRONTTRACK_USE_TBB) + +# libraries to link to +SET(_link_LIBRARIES + ${CAS_TKShHealing} + ${CAS_TKernel} + ${CAS_TKBRep} + ${CAS_TKG3d} + ${CAS_TKTopAlgo} + ${CAS_TKGeomBase} + ${CAS_TKGeomAlgo} + ${GEOM_XAO} + ${MEDCoupling_medloader} + ${TBB_LIBS} + ${Boost_LIBRARIES} +) + +# --- headers --- + +# header files +SET(FRONTTRACK_HEADERS + FrontTrack.hxx +) + +# --- sources --- + +# sources / static +SET(FRONTTRACK_SOURCES + FrontTrack.cxx + FrontTrack_NodeGroups.cxx + FrontTrack_NodesOnGeom.cxx + FrontTrack_Projector.cxx + FrontTrack_Utils.cxx +) + +# --- rules --- + +ADD_LIBRARY(FrontTrack ${FRONTTRACK_SOURCES}) +TARGET_LINK_LIBRARIES(FrontTrack ${_link_LIBRARIES} ) +INSTALL(TARGETS FrontTrack EXPORT ${PROJECT_NAME}TargetGroup DESTINATION ${SALOME_INSTALL_LIBS}) + +INSTALL(FILES ${FRONTTRACK_HEADERS} DESTINATION ${SALOME_INSTALL_HEADERS}) diff --git a/src/FrontTrack/FrontTrack.cxx b/src/FrontTrack/FrontTrack.cxx new file mode 100755 index 00000000..43ee23f1 --- /dev/null +++ b/src/FrontTrack/FrontTrack.cxx @@ -0,0 +1,97 @@ +// File : FrontTrack.cxx +// Created : Tue Apr 25 17:20:28 2017 +// Author : Edward AGAPOV (eap) + +#include "FrontTrack.hxx" +#include "FrontTrack_NodeGroups.hxx" +#include "FrontTrack_Utils.hxx" + +#include +#include +#include + +#include +#include + +#include + +#include + +/*! + * \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] theNodeFiles - 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 > & theNodeFiles, + const std::string& theXaoFileName, + bool theIsParallel ) +{ + // check arguments + + if ( theNodeFiles.empty() ) + return; + + if ( !FT_Utils::fileExists( theInputMedFile )) + throw std::invalid_argument( "Input MED file does not exist: " + theInputMedFile ); + + if ( !FT_Utils::canWrite( theOutputMedFile )) + throw std::invalid_argument( "Can't create the output MED file: " + theOutputMedFile ); + + for ( size_t i = 0; i < theNodeFiles.size(); ++i ) + if ( !FT_Utils::fileExists( theNodeFiles[i] )) + throw std::invalid_argument( "Input node file does not exist: " + theNodeFiles[i] ); + + if ( !FT_Utils::fileExists( theXaoFileName )) + throw std::invalid_argument( "Input XAO file does not exist: " + theXaoFileName ); + + + // read a mesh + + 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 + + XAO::Xao xao; + if ( !xao.importXAO( theXaoFileName ) || !xao.getGeometry() ) + throw std::invalid_argument( "Failed to read the XAO input file: " + theXaoFileName ); + + 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) + + FT_NodeGroups nodeGroups; + nodeGroups.read( theNodeFiles, &xao, nodeCoords ); + + + // project nodes to the boundary shapes and change their coordinates + + OSD_Parallel::For( 0, nodeGroups.nbOfGroups(), nodeGroups, !theIsParallel ); + + // save the modified mesh + + const int erase = 2; + mfMesh->write( theOutputMedFile, /*mode=*/erase ); + + if ( !nodeGroups.isOK() ) + throw std::runtime_error("Unable to project some nodes"); +} diff --git a/src/FrontTrack/FrontTrack.hxx b/src/FrontTrack/FrontTrack.hxx new file mode 100755 index 00000000..049cb4bd --- /dev/null +++ b/src/FrontTrack/FrontTrack.hxx @@ -0,0 +1,36 @@ +// File : FrontTrack.hxx +// Created : Tue Apr 25 17:08:52 2017 +// Author : Edward AGAPOV (eap) + + +#ifndef __FrontTrack_HXX__ +#define __FrontTrack_HXX__ + +#include +#include + +class FrontTrack +{ +public: + + /*! + * \brief Relocate nodes to lie on geometry + * \param [in] inputMedFile - a MED file holding a mesh including nodes that will be + * moved onto the geometry + * \param [in] outputMedFile - a MED file to create, that will hold a modified mesh + * \param [in] nodeFiles - an array of names of files describing groups of nodes that + * will be moved onto the geometry + * \param [in] xaoFileName - a path to a file in XAO format holding the geometry and + * the geometrical groups. + * \param [in] isParallel - if \c true, all processors are used to treat boundary shapes + * in parallel. + */ + void track( const std::string& inputMedFile, + const std::string& outputMedFile, + const std::vector< std::string > & nodeFiles, + const std::string& xaoFileName, + bool isParallel=true); + +}; + +#endif diff --git a/src/FrontTrack/FrontTrack_NodeGroups.cxx b/src/FrontTrack/FrontTrack_NodeGroups.cxx new file mode 100755 index 00000000..917a056a --- /dev/null +++ b/src/FrontTrack/FrontTrack_NodeGroups.cxx @@ -0,0 +1,113 @@ +// File : FrontTrack_NodeGroups.cxx +// Created : Tue Apr 25 19:17:47 2017 +// Author : Edward AGAPOV (eap) + +#include "FrontTrack_NodeGroups.hxx" +#include "FrontTrack_Projector.hxx" +#include "FrontTrack_Utils.hxx" + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace +{ + //================================================================================ + /*! + * \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 ); + + 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; + } +} diff --git a/src/FrontTrack/FrontTrack_NodeGroups.hxx b/src/FrontTrack/FrontTrack_NodeGroups.hxx new file mode 100755 index 00000000..ac03e73e --- /dev/null +++ b/src/FrontTrack/FrontTrack_NodeGroups.hxx @@ -0,0 +1,58 @@ +// File : FrontTrack_NodeGroups.hxx +// Created : Tue Apr 25 19:02:49 2017 +// Author : Edward AGAPOV (eap) + +#ifndef __FrontTrack_NodeGroups_HXX__ +#define __FrontTrack_NodeGroups_HXX__ + +#include "FrontTrack_NodesOnGeom.hxx" +#include "FrontTrack_Projector.hxx" + +#include +#include + +namespace MEDCoupling { + class DataArrayDouble; +} +namespace XAO { + class Xao; +} + +/*! + * \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 + +}; + +#endif diff --git a/src/FrontTrack/FrontTrack_NodesOnGeom.cxx b/src/FrontTrack/FrontTrack_NodesOnGeom.cxx new file mode 100755 index 00000000..d150c887 --- /dev/null +++ b/src/FrontTrack/FrontTrack_NodesOnGeom.cxx @@ -0,0 +1,465 @@ +// File : FrontTrack_NodesOnGeom.cxx +// Created : Tue Apr 25 20:48:23 2017 +// Author : Edward AGAPOV (eap) + +#include "FrontTrack_NodesOnGeom.hxx" +#include "FrontTrack_Utils.hxx" + +#include + +#include +#include +#include +#include + +namespace +{ + /*! + * \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 "fr2D.**" + 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 ); + + // ------------------------------------- + // 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]; + + // 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; + + // 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; + + 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 + } + } + } + } + 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] ); + double maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.; + if ( _projectors[ 0 ].project( xyz, maxDist2, newXyz, + nn._params, nn._nearParams )) + moveNode( nn._nodeToMove, newXyz ); + else + notProjectedNodes.push_back( &nn ); + } + } + + + 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(); + + 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] ); +} diff --git a/src/FrontTrack/FrontTrack_NodesOnGeom.hxx b/src/FrontTrack/FrontTrack_NodesOnGeom.hxx new file mode 100755 index 00000000..664d5849 --- /dev/null +++ b/src/FrontTrack/FrontTrack_NodesOnGeom.hxx @@ -0,0 +1,99 @@ +// File : FrontTrack_NodesOnGeom.hxx +// Created : Tue Apr 25 19:12:25 2017 +// Author : Edward AGAPOV (eap) + + +#ifndef __FrontTrack_NodesOnGeom_HXX__ +#define __FrontTrack_NodesOnGeom_HXX__ + +#include "FrontTrack_Projector.hxx" + +#include +#include +#include +#include + +#include +#include + +namespace FT_Utils { + struct XaoGroups; +} +namespace MEDCoupling { + class DataArrayDouble; +} +namespace XAO { + class BrepGeometry; +} + + //-------------------------------------------------------------------------------------------- +/*! + * \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; + +}; + +#endif diff --git a/src/FrontTrack/FrontTrack_Projector.cxx b/src/FrontTrack/FrontTrack_Projector.cxx new file mode 100755 index 00000000..ce24d4b0 --- /dev/null +++ b/src/FrontTrack/FrontTrack_Projector.cxx @@ -0,0 +1,905 @@ +// File : FrontTrack_Projector.cxx +// Created : Wed Apr 26 20:33:55 2017 +// Author : Edward AGAPOV (eap) + +#include "FrontTrack_Projector.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +//----------------------------------------------------------------------------- +/*! + * \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) + { + 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 ); + 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) + { + 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 ); +} + +//================================================================================ +/*! + * \brief Copy another projector + */ +//================================================================================ + +FT_Projector::FT_Projector(const FT_Projector& other) +{ + _realProjector = 0; + _shape = other._shape; + _bndBox = other._bndBox; +} + +//================================================================================ +/*! + * \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; + } + } +} + +//================================================================================ +/*! + * \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 ); + + return ( _realProjector->_dist * _realProjector->_dist < maxDist2 ); +} + +//================================================================================ +/*! + * \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; + + return _realProjector->projectAndClassify( point, maxDist2, projection, + newSolution, prevSolution ); +} + +//================================================================================ +/*! + * \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() ); +} diff --git a/src/FrontTrack/FrontTrack_Projector.hxx b/src/FrontTrack/FrontTrack_Projector.hxx new file mode 100755 index 00000000..f7b0aa76 --- /dev/null +++ b/src/FrontTrack/FrontTrack_Projector.hxx @@ -0,0 +1,70 @@ +// File : FrontTrack_Projector.hxx +// Created : Wed Apr 26 20:12:13 2017 +// Author : Edward AGAPOV (eap) + +#ifndef __FrontTrack_Projector_HXX__ +#define __FrontTrack_Projector_HXX__ + +#include +#include + +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; + + // 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; +}; + +#endif diff --git a/src/FrontTrack/FrontTrack_Utils.cxx b/src/FrontTrack/FrontTrack_Utils.cxx new file mode 100755 index 00000000..129d4f37 --- /dev/null +++ b/src/FrontTrack/FrontTrack_Utils.cxx @@ -0,0 +1,146 @@ +// File : FrontTrack_Utils.cxx +// Created : Tue Apr 25 17:28:58 2017 +// Author : Edward AGAPOV (eap) + +#include "FrontTrack_Utils.hxx" + +#include +#include + +#include +#include + +namespace boofs = boost::filesystem; + +//================================================================================ +/* + * \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 < allProjectors.size() ) + groupProjectors.push_back ( & allProjectors[ *id ]); + } + + return nbFound; +} diff --git a/src/FrontTrack/FrontTrack_Utils.hxx b/src/FrontTrack/FrontTrack_Utils.hxx new file mode 100755 index 00000000..b557ef55 --- /dev/null +++ b/src/FrontTrack/FrontTrack_Utils.hxx @@ -0,0 +1,54 @@ +// File : FrontTrack_Utils.hxx +// Created : Tue Apr 25 17:23:33 2017 +// Author : Edward AGAPOV (eap) + +#ifndef __FrontTrack_Utils_HXX__ +#define __FrontTrack_Utils_HXX__ + +#include "FrontTrack_Projector.hxx" + +#include +#include +#include + +namespace XAO { + class Xao; + class Group; +} + +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 + }; +} + +#endif diff --git a/src/FrontTrack_SWIG/CMakeLists.txt b/src/FrontTrack_SWIG/CMakeLists.txt new file mode 100755 index 00000000..735bdcfc --- /dev/null +++ b/src/FrontTrack_SWIG/CMakeLists.txt @@ -0,0 +1,24 @@ + +INCLUDE(${SWIG_USE_FILE}) + +ADD_DEFINITIONS(${PYTHON_DEFINITIONS}) + +SET_SOURCE_FILES_PROPERTIES(FrontTrack.i PROPERTIES CPLUSPLUS ON) +SET_SOURCE_FILES_PROPERTIES(FrontTrack.i PROPERTIES SWIG_DEFINITIONS "-shadow") + +INCLUDE_DIRECTORIES( + ${PYTHON_INCLUDE_DIRS} + ${PTHREAD_INCLUDE_DIR} # pthread dependancy due to python2.7 library + ${CMAKE_CURRENT_SOURCE_DIR} + ${CMAKE_CURRENT_SOURCE_DIR}/../FrontTrack + ) + +SWIG_ADD_MODULE(FrontTrack python FrontTrack.i) +SWIG_LINK_LIBRARIES(FrontTrack ${PYTHON_LIBRARIES} ${PLATFORM_LIBS} FrontTrack) + +IF(WIN32) + SET_TARGET_PROPERTIES(_FrontTrack PROPERTIES DEBUG_OUTPUT_NAME _FrontTrack_d) +ENDIF(WIN32) +INSTALL(TARGETS ${SWIG_MODULE_FrontTrack_REAL_NAME} DESTINATION ${SALOME_INSTALL_PYTHON}) + +SALOME_INSTALL_SCRIPTS(${CMAKE_CURRENT_BINARY_DIR}/FrontTrack.py ${SALOME_INSTALL_PYTHON}) diff --git a/src/FrontTrack_SWIG/FrontTrack.i b/src/FrontTrack_SWIG/FrontTrack.i new file mode 100755 index 00000000..8ec6b17f --- /dev/null +++ b/src/FrontTrack_SWIG/FrontTrack.i @@ -0,0 +1,49 @@ +// File : FrontTrack.i +// Created : Fri Apr 28 17:36:20 2017 +// Author : Edward AGAPOV (eap) + +%module FrontTrack + +%{ +#include "FrontTrack.hxx" +#include +#include +#include + +static PyObject* setOCCException(Standard_Failure& ex) +{ + std::string msg(ex.DynamicType()->Name()); + if ( ex.GetMessageString() && strlen( ex.GetMessageString() )) { + msg += ": "; + msg += ex.GetMessageString(); + } + PyErr_SetString(PyExc_Exception, msg.c_str() ); + return NULL; +} + +%} + + +%exception +{ + try { + OCC_CATCH_SIGNALS; + $action } + catch (Standard_Failure& ex) { + return setOCCException(ex); + } + catch (std::exception& ex) { + PyErr_SetString(PyExc_Exception, ex.what() ); + return NULL; + } +} + +%include +%include + +%template(svec) std::vector; + +//%feature("autodoc", "1"); +//%feature("docstring"); + +%include "FrontTrack.hxx"