]> SALOME platform Git repositories - modules/smesh.git/commitdiff
Salome HOME
Suivi de frontière - départ
authorGérald NICOLAS <gerald.nicolas@edf.fr>
Wed, 20 May 2020 15:36:39 +0000 (17:36 +0200)
committerGérald NICOLAS <gerald.nicolas@edf.fr>
Wed, 20 May 2020 15:36:39 +0000 (17:36 +0200)
13 files changed:
src/ADAPTFrontTrack/CMakeLists.txt [new file with mode: 0644]
src/ADAPTFrontTrack/FrontTrack.cxx [new file with mode: 0644]
src/ADAPTFrontTrack/FrontTrack.hxx [new file with mode: 0644]
src/ADAPTFrontTrack/FrontTrack_NodeGroups.cxx [new file with mode: 0644]
src/ADAPTFrontTrack/FrontTrack_NodeGroups.hxx [new file with mode: 0644]
src/ADAPTFrontTrack/FrontTrack_NodesOnGeom.cxx [new file with mode: 0644]
src/ADAPTFrontTrack/FrontTrack_NodesOnGeom.hxx [new file with mode: 0644]
src/ADAPTFrontTrack/FrontTrack_Projector.cxx [new file with mode: 0644]
src/ADAPTFrontTrack/FrontTrack_Projector.hxx [new file with mode: 0644]
src/ADAPTFrontTrack/FrontTrack_Utils.cxx [new file with mode: 0644]
src/ADAPTFrontTrack/FrontTrack_Utils.hxx [new file with mode: 0644]
src/ADAPTFrontTrack_SWIG/CMakeLists.txt [new file with mode: 0644]
src/ADAPTFrontTrack_SWIG/FrontTrack.i [new file with mode: 0644]

diff --git a/src/ADAPTFrontTrack/CMakeLists.txt b/src/ADAPTFrontTrack/CMakeLists.txt
new file mode 100644 (file)
index 0000000..b19d068
--- /dev/null
@@ -0,0 +1,73 @@
+# Copyright (C) 2017-2020  CEA/DEN, EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+
+# --- options ---
+# additional include directories
+INCLUDE_DIRECTORIES(
+  ${OpenCASCADE_INCLUDE_DIR}
+  ${GEOM_INCLUDE_DIRS}
+  ${MEDFILE_INCLUDE_DIRS}
+  ${MEDCOUPLING_INCLUDE_DIRS}
+  ${TBB_INCLUDE_DIRS}
+)
+
+# additional preprocessor / compiler flags
+ADD_DEFINITIONS(
+  ${OpenCASCADE_DEFINITIONS}
+)
+
+IF(ADAPTFRONTTRACK_USE_TBB)
+  SET(TBB_LIBS ${TBB_LIBRARIES})
+ENDIF(ADAPTFRONTTRACK_USE_TBB)
+
+# libraries to link to
+SET(_link_LIBRARIES
+  ${OpenCASCADE_FoundationClasses_LIBRARIES}
+  ${OpenCASCADE_ModelingAlgorithms_LIBRARIES}
+  ${OpenCASCADE_ModelingData_LIBRARIES}
+  ${GEOM_XAO}
+  ${MEDCoupling_medloader}
+  ${TBB_LIBS}
+)
+
+# --- headers ---
+
+# header files
+SET(ADAPTFRONTTRACK_HEADERS
+  FrontTrack.hxx
+)
+
+# --- sources ---
+
+# sources / static
+SET(ADAPTFRONTTRACK_SOURCES
+  FrontTrack.cxx
+  FrontTrack_NodeGroups.cxx
+  FrontTrack_NodesOnGeom.cxx
+  FrontTrack_Projector.cxx
+  FrontTrack_Utils.cxx
+)
+
+# --- rules ---
+
+ADD_LIBRARY(ADAPTFrontTrack ${ADAPTFRONTTRACK_SOURCES})
+TARGET_LINK_LIBRARIES(ADAPTFrontTrack ${_link_LIBRARIES} )
+INSTALL(TARGETS ADAPTFrontTrack EXPORT ${PROJECT_NAME}TargetGroup DESTINATION ${SALOME_INSTALL_LIBS})
+
+INSTALL(FILES ${ADAPTFRONTTRACK_HEADERS} DESTINATION ${SALOME_INSTALL_HEADERS})
diff --git a/src/ADAPTFrontTrack/FrontTrack.cxx b/src/ADAPTFrontTrack/FrontTrack.cxx
new file mode 100644 (file)
index 0000000..91eb4b6
--- /dev/null
@@ -0,0 +1,164 @@
+// Copyright (C) 2017-2020  CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// 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 <MCAuto.hxx>
+#include <MEDCouplingMemArray.hxx>
+#include <MEDFileMesh.hxx>
+
+#include <XAO_Xao.hxx>
+#include <XAO_BrepGeometry.hxx>
+
+#include <stdexcept>
+
+#include <OSD_Parallel.hxx>
+
+/*!
+ * \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");
+}
diff --git a/src/ADAPTFrontTrack/FrontTrack.hxx b/src/ADAPTFrontTrack/FrontTrack.hxx
new file mode 100644 (file)
index 0000000..fa03316
--- /dev/null
@@ -0,0 +1,54 @@
+// Copyright (C) 2017-2020  CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File      : FrontTrack.hxx
+// Created   : Tue Apr 25 17:08:52 2017
+// Author    : Edward AGAPOV (eap)
+
+
+#ifndef __FrontTrack_HXX__
+#define __FrontTrack_HXX__
+
+#include <vector>
+#include <string>
+
+class FrontTrack
+{
+public:
+
+  /*!
+   * \brief Relocate nodes to lie on geometry
+   *  \param [in] theInputMedFile - a MED file holding a mesh including nodes that will be
+   *         moved onto the geometry
+   *  \param [in] theOutputMedFile - a MED file to create, that will hold a modified mesh
+   *  \param [in] theInputNodeFiles - an array of names of files describing groups of nodes that
+   *         will be moved onto the geometry
+   *  \param [in] theXaoFileName - a path to a file in XAO format  holding the geometry and
+   *         the geometrical groups.
+   *  \param [in] theIsParallel - if \c true, all processors are used to treat boundary shapes
+   *          in parallel.
+   */
+  void track( const std::string&                 theInputMedFile,
+              const std::string&                 theOutputMedFile,
+              const std::vector< std::string > & theInputNodeFiles,
+              const std::string&                 theXaoFileName,
+              bool                               theIsParallel=true);
+
+};
+
+#endif
diff --git a/src/ADAPTFrontTrack/FrontTrack_NodeGroups.cxx b/src/ADAPTFrontTrack/FrontTrack_NodeGroups.cxx
new file mode 100644 (file)
index 0000000..247def2
--- /dev/null
@@ -0,0 +1,134 @@
+// Copyright (C) 2017-2020  CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// 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 <MEDCouplingMemArray.hxx>
+#include <XAO_BrepGeometry.hxx>
+#include <XAO_Xao.hxx>
+
+#include <BRepBndLib.hxx>
+#include <Bnd_Box.hxx>
+#include <TopExp.hxx>
+#include <TopExp_Explorer.hxx>
+#include <TopTools_IndexedMapOfShape.hxx>
+#include <TopoDS_Shape.hxx>
+
+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 );
+#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;
+  }
+}
diff --git a/src/ADAPTFrontTrack/FrontTrack_NodeGroups.hxx b/src/ADAPTFrontTrack/FrontTrack_NodeGroups.hxx
new file mode 100644 (file)
index 0000000..9204ac8
--- /dev/null
@@ -0,0 +1,76 @@
+// Copyright (C) 2017-2020  CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// 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 <vector>
+#include <string>
+
+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/ADAPTFrontTrack/FrontTrack_NodesOnGeom.cxx b/src/ADAPTFrontTrack/FrontTrack_NodesOnGeom.cxx
new file mode 100644 (file)
index 0000000..1642234
--- /dev/null
@@ -0,0 +1,511 @@
+// Copyright (C) 2017-2020  CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// 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 <MEDCouplingMemArray.hxx>
+
+#include <cstdio>
+#include <cstdlib>
+#include <list>
+#include <stdexcept>
+
+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 "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] );
+}
diff --git a/src/ADAPTFrontTrack/FrontTrack_NodesOnGeom.hxx b/src/ADAPTFrontTrack/FrontTrack_NodesOnGeom.hxx
new file mode 100644 (file)
index 0000000..08db0e1
--- /dev/null
@@ -0,0 +1,117 @@
+// Copyright (C) 2017-2020  CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// 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 <Bnd_Box.hxx>
+#include <NCollection_DataMap.hxx>
+#include <TopoDS_Shape.hxx>
+#include <TColStd_DataMapOfIntegerInteger.hxx>
+
+#include <string>
+#include <vector>
+
+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<int>              _nodesOrder;
+
+};
+
+#endif
diff --git a/src/ADAPTFrontTrack/FrontTrack_Projector.cxx b/src/ADAPTFrontTrack/FrontTrack_Projector.cxx
new file mode 100644 (file)
index 0000000..a3bbfc6
--- /dev/null
@@ -0,0 +1,947 @@
+// Copyright (C) 2017-2020  CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File      : FrontTrack_Projector.cxx
+// Created   : Wed Apr 26 20:33:55 2017
+// Author    : Edward AGAPOV (eap)
+
+#include "FrontTrack_Projector.hxx"
+
+#include <BRepAdaptor_Curve.hxx>
+#include <BRepBndLib.hxx>
+#include <BRepTopAdaptor_FClass2d.hxx>
+#include <BRep_Tool.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 <TopExp.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopoDS_Vertex.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>
+
+#include <limits>
+
+//-----------------------------------------------------------------------------
+/*!
+ * \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() );
+}
diff --git a/src/ADAPTFrontTrack/FrontTrack_Projector.hxx b/src/ADAPTFrontTrack/FrontTrack_Projector.hxx
new file mode 100644 (file)
index 0000000..82578e2
--- /dev/null
@@ -0,0 +1,95 @@
+// Copyright (C) 2017-2020  CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// 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 <TopoDS_Shape.hxx>
+#include <Bnd_Box.hxx>
+
+struct FT_RealProjector;
+
+/*!
+ * \brief Projector of a point to a boundary shape. Wrapper of a real projection algo
+ */
+class FT_Projector
+{
+public:
+
+  FT_Projector(const TopoDS_Shape& shape = TopoDS_Shape());
+  FT_Projector(const FT_Projector& other);
+  ~FT_Projector();
+
+  // initialize with a boundary shape, compute the bounding box
+  void setBoundaryShape(const TopoDS_Shape& shape);
+
+  // return the boundary shape
+  const TopoDS_Shape& getShape() const { return _shape; }
+
+  // return the bounding box
+  const Bnd_Box getBoundingBox() const { return _bndBox; }
+
+
+  // create a real projector
+  void prepareForProjection();
+
+  // return true if a previously found solution can be used to speed up the projection
+  bool canUsePrevSolution() const;
+
+  // return true if projection is not needed
+  bool isPlanarBoundary() const;
+
+
+  // switch a mode of usage of prevSolution.
+  // If projection fails, to try to project without usage of prevSolution.
+  // By default this mode is off
+  void tryWithoutPrevSolution( bool toTry ) { _tryWOPrevSolution = toTry; }
+
+  // project a point to the boundary shape
+  bool project( const gp_Pnt& point,
+                const double  maxDist2,
+                gp_Pnt&       projection,
+                double*       newSolution,
+                const double* prevSolution = 0);
+
+  // project a point to the boundary shape and check if the projection is within the shape boundary
+  bool projectAndClassify( const gp_Pnt& point,
+                           const double  maxDist2,
+                           gp_Pnt&       projection,
+                           double*       newSolution,
+                           const double* prevSolution = 0);
+
+  // check if a point lies on the boundary shape
+  bool isOnShape( const gp_Pnt& point,
+                  const double  tol2,
+                  double*       newSolution,
+                  const double* prevSolution = 0);
+
+private:
+
+  FT_RealProjector* _realProjector;
+  Bnd_Box           _bndBox;
+  TopoDS_Shape      _shape;
+  bool              _tryWOPrevSolution;
+};
+
+#endif
diff --git a/src/ADAPTFrontTrack/FrontTrack_Utils.cxx b/src/ADAPTFrontTrack/FrontTrack_Utils.cxx
new file mode 100644 (file)
index 0000000..e74072a
--- /dev/null
@@ -0,0 +1,164 @@
+// Copyright (C) 2017-2020  CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File      : FrontTrack_Utils.cxx
+// Created   : Tue Apr 25 17:28:58 2017
+// Author    : Edward AGAPOV (eap)
+
+#include "FrontTrack_Utils.hxx"
+
+#include <XAO_Xao.hxx>
+#include <XAO_Group.hxx>
+
+#include <fcntl.h>
+#include <boost/filesystem.hpp>
+
+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<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;
+}
diff --git a/src/ADAPTFrontTrack/FrontTrack_Utils.hxx b/src/ADAPTFrontTrack/FrontTrack_Utils.hxx
new file mode 100644 (file)
index 0000000..38817cb
--- /dev/null
@@ -0,0 +1,72 @@
+// Copyright (C) 2017-2020  CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// 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 <vector>
+#include <string>
+#include <map>
+
+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/ADAPTFrontTrack_SWIG/CMakeLists.txt b/src/ADAPTFrontTrack_SWIG/CMakeLists.txt
new file mode 100644 (file)
index 0000000..1154382
--- /dev/null
@@ -0,0 +1,83 @@
+# Copyright (C) 2017-2020  CEA/DEN, EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+
+INCLUDE(${SWIG_USE_FILE})
+
+# --- options ---
+# additional include directories
+INCLUDE_DIRECTORIES(
+  ${PYTHON_INCLUDE_DIRS}
+  ${PTHREAD_INCLUDE_DIR} # pthread dependancy due to python2.7 library
+  ${CMAKE_CURRENT_SOURCE_DIR}
+  ${PROJECT_SOURCE_DIR}/src/ADAPTFrontTrack
+)
+
+# additional preprocessor / compiler flags
+ADD_DEFINITIONS(
+  ${PYTHON_DEFINITIONS}
+)
+
+# libraries to link to
+SET(_link_LIBRARIES
+  ${PYTHON_LIBRARIES}
+  ${PLATFORM_LIBS}
+  ADAPTFrontTrack
+)
+
+# --- headers ---
+
+# header files / to be processed by moc
+SET(ADAPTFrontTrack_Swig_HEADERS
+  FrontTrack_Swig.i
+)
+
+# --- sources ---
+
+# sources / static
+SET(ADAPTFrontTrack_Swig_SOURCES
+  ${ADAPTFrontTrack_Swig_HEADERS}
+)
+
+# swig flags
+SET_SOURCE_FILES_PROPERTIES(FrontTrack_Swig.i PROPERTIES CPLUSPLUS ON)
+SET_SOURCE_FILES_PROPERTIES(FrontTrack_Swig.i PROPERTIES SWIG_DEFINITIONS "-shadow")
+
+# --- scripts ---
+
+# scripts / swig wrappings
+SET(_swig_SCRIPTS
+  ${CMAKE_CURRENT_BINARY_DIR}/FrontTrack.py 
+)
+
+# --- rules ---
+IF(${CMAKE_VERSION} VERSION_LESS "3.8.0")
+  SWIG_ADD_MODULE(FrontTrack_Swig python ${ADAPTFrontTrack_Swig_SOURCES})
+ELSE()
+  SWIG_ADD_LIBRARY(FrontTrack_Swig LANGUAGE python SOURCES ${ADAPTFrontTrack_Swig_SOURCES})
+ENDIF()
+
+SWIG_LINK_LIBRARIES(FrontTrack_Swig ${_link_LIBRARIES})
+SWIG_CHECK_GENERATION(FrontTrack_Swig)
+IF(WIN32)
+  SET_TARGET_PROPERTIES(_FrontTrack_Swig PROPERTIES DEBUG_OUTPUT_NAME _FrontTrack_Swig_d)
+ENDIF(WIN32)
+
+INSTALL(TARGETS ${SWIG_MODULE_FrontTrack_Swig_REAL_NAME} DESTINATION ${SALOME_INSTALL_PYTHON})
+
+SALOME_INSTALL_SCRIPTS("${_swig_SCRIPTS}" ${SALOME_INSTALL_PYTHON} EXTRA_DPYS "${SWIG_MODULE_FrontTrack_Swig_REAL_NAME}")
diff --git a/src/ADAPTFrontTrack_SWIG/FrontTrack.i b/src/ADAPTFrontTrack_SWIG/FrontTrack.i
new file mode 100644 (file)
index 0000000..50c6c20
--- /dev/null
@@ -0,0 +1,67 @@
+// Copyright (C) 2017-2020  CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File      : FrontTrack.i
+// Created   : Fri Apr 28 17:36:20 2017
+// Author    : Edward AGAPOV (eap)
+
+%module FrontTrack
+
+%{
+#include "FrontTrack.hxx"
+#include <Standard_Failure.hxx>
+#include <Standard_ErrorHandler.hxx>
+#include <stdexcept>
+
+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 <std_string.i>
+%include <std_vector.i>
+
+%template(svec) std::vector<std::string>;
+
+//%feature("autodoc", "1");
+//%feature("docstring");
+
+%include "FrontTrack.hxx"