# Ensure a proper linker behavior:
CMAKE_POLICY(SET CMP0003 NEW)
+IF(WIN32)
+ CMAKE_POLICY(SET CMP0020 OLD) # disable automatic linking to qtmain.lib
+ENDIF(WIN32)
# Versioning
# ===========
STRING(TOUPPER ${PROJECT_NAME} PROJECT_NAME_UC)
SET(${PROJECT_NAME_UC}_MAJOR_VERSION 8)
-SET(${PROJECT_NAME_UC}_MINOR_VERSION 3)
+SET(${PROJECT_NAME_UC}_MINOR_VERSION 4)
SET(${PROJECT_NAME_UC}_PATCH_VERSION 0)
SET(${PROJECT_NAME_UC}_VERSION
${${PROJECT_NAME_UC}_MAJOR_VERSION}.${${PROJECT_NAME_UC}_MINOR_VERSION}.${${PROJECT_NAME_UC}_PATCH_VERSION})
FIND_PACKAGE(SalomePythonInterp REQUIRED)
FIND_PACKAGE(SalomePythonLibs REQUIRED)
-# Other KERNEL optionals:
+# Other KERNEL optionals:
IF(SALOME_BUILD_TESTS)
ENABLE_TESTING()
FIND_PACKAGE(SalomeCppUnit)
ENDIF()
+# Advanced options:
+OPTION(FRONTTRACK_USE_TBB "Enable parallel computation using TBB" OFF)
+MARK_AS_ADVANCED(FRONTTRACK_USE_TBB)
+
+# Prerequisites
+# =============
+
+# MEDCoupling
+SET(MEDCOUPLING_ROOT_DIR $ENV{MEDCOUPLING_ROOT_DIR} CACHE PATH "Path to the MEDCoupling tool")
+IF(EXISTS ${MEDCOUPLING_ROOT_DIR})
+ LIST(APPEND CMAKE_MODULE_PATH "${MEDCOUPLING_ROOT_DIR}/cmake_files")
+ FIND_PACKAGE(SalomeMEDCoupling REQUIRED) # will reload HDF5, MEDFile, XDR, etc ...
+ELSE(EXISTS ${MEDCOUPLING_ROOT_DIR})
+ MESSAGE(FATAL_ERROR "We absolutely need the MEDCoupling tool, please define MEDCOUPLING_ROOT_DIR !")
+ENDIF(EXISTS ${MEDCOUPLING_ROOT_DIR})
+
##
## From KERNEL:
##
IF(EXISTS ${GEOM_ROOT_DIR})
LIST(APPEND CMAKE_MODULE_PATH "${GEOM_ROOT_DIR}/adm_local/cmake_files")
FIND_PACKAGE(SalomeGEOM REQUIRED)
+ ADD_DEFINITIONS(${GEOM_DEFINITIONS})
+ INCLUDE_DIRECTORIES(${GEOM_INCLUDE_DIRS})
ELSE(EXISTS ${GEOM_ROOT_DIR})
MESSAGE(FATAL_ERROR "We absolutely need a Salome GEOM, please define GEOM_ROOT_DIR")
ENDIF(EXISTS ${GEOM_ROOT_DIR})
# MedFile
FIND_PACKAGE(SalomeMEDFile REQUIRED)
+# OpenCASCADE
+FIND_PACKAGE(SalomeCAS REQUIRED)
+
+# TBB
+IF(FRONTTRACK_USE_TBB)
+ FIND_PACKAGE(SalomeTBB)
+ SALOME_LOG_OPTIONAL_PACKAGE(TBB FRONTTRACK_USE_TBB)
+ ADD_DEFINITIONS(-DHAVE_TBB)
+ENDIF(FRONTTRACK_USE_TBB)
+
+# Python
+FIND_PACKAGE(SalomePythonInterp REQUIRED)
+# SWIG
+FIND_PACKAGE(SalomeSWIG REQUIRED)
+# Boost
+FIND_PACKAGE(SalomeBoost REQUIRED)
+
# Detection summary:
SALOME_PACKAGE_REPORT_AND_CHECK()
For a reference to HOMARD, use:
+G. Nicolas, T. Fouquet, S. Geniaut, S. Cuvilliez, Improved Adaptive Mesh Refinement for Conformal Hexahedral Meshes, "*Advances in Engineering Software*", Vol. 102, pp. 14-28, 2016, doi:10.1016/j.advengsoft.2016.07.014
+
+`Available here <http://dx.doi.org/10.1016/j.advengsoft.2016.07.014>`_
+
G. Nicolas and T. Fouquet, Adaptive Mesh Refinement for Conformal Hexahedral Meshes, "*Finite Elements in Analysis and Design*", Vol. 67, pp. 1-12, 2013, doi:10.1016/j.finel.2012.11.008
`Available here <http://dx.doi.org/10.1016/j.finel.2012.11.008>`_
SALOME RELEASE 7.1:
- Correction of a bug in the filtering of the refinement by groups
- - Driving of the adaptation by a threshold that is based upon the mean and the standard deviation (see :ref:`tui_create_hypothese`)
- - TUI: the last time step for the driving instant can be selected by SetTimeStepRankLast (see :ref:`tui_create_iteration`)
+ - Driving of the adaptation by a threshold that is based upon the mean and the standard deviation (see :doc:`tui_create_hypothese`)
+ - TUI: the last time step for the driving instant can be selected by SetTimeStepRankLast (see :doc:`tui_create_iteration`)
- The destruction of the objects is available (GUI and TUI)
- The TUI functions are modified: they act onto the objects and no more by the names:
- ``homard.AssociateIterHypo(iter_name,hypo_name)`` becomes ``iter.AssociateHypo(hypo_name)``, ``homard.AssociateHypoZone(hypo_name,zone_name,type_use)`` becomes ``hypo.AddZone(zone_name,type_use)``, ``homard.CreateIteration(iter_name,iter_parent_name)`` becomes ``case.NextIteration(iter_name)`` ou ``iter.NextIteration(iter_name)``
+ ``homard.AssociateIterHypo(iter_name,hypo_name)`` becomes ``iter.AssociateHypo(hypo_name)``, ``homard.AssociateHypoZone(hypo_name,zone_name,type_use)`` becomes ``hypo.AddZone(zone_name,type_use)``, ``homard.CreateIteration(iter_name,iter_parent_name)`` becomes ``case.NextIteration(iter_name)`` or ``iter.NextIteration(iter_name)``
- Addition of a function:
- ``cas.LastIteration()`` : returns the last iteration of the case (see :ref:`tui_create_iteration`)
- - Mesh analysis (see :ref:`gui_mesh_info`)
+ ``cas.LastIteration()`` : returns the last iteration of the case (see :doc:`tui_create_iteration`)
+ - Mesh analysis (see :doc:`gui_mesh_info`)
SALOME RELEASE 7.2:
- - A cone can be used as a support for a 2D boundary (see :ref:`tui_create_boundary`)
- - The result can be published or not published in SMESH (see :ref:`tui_create_iteration`)
+ - A cone can be used as a support for a 2D boundary (see :doc:`tui_create_boundary`)
+ - The result can be published or not published in SMESH (see :doc:`tui_create_iteration`)
``iter.Compute(option)`` becomes ``iter.Compute(option1, option2)``
- - A serie of iterations that is stored in a directory can be followed (see :ref:`gui_create_case`)
+ - A serie of iterations that is stored in a directory can be followed (see :doc:`gui_create_case`)
SALOME RELEASE 7.3:
- English documentation
- - Options to track the computation (see :ref:`tui_create_iteration`): ``iter.SetInfoCompute(MessInfo)``
+ - Options to track the computation (see :doc:`tui_create_iteration`): ``iter.SetInfoCompute(MessInfo)``
- Automatic creation of YACS schema:
- ``case.WriteYACSSchema(ScriptFile, DirName, MeshFile)``: writes a schema YACS relative to the case (see :ref:`tui_create_case`)
+ ``case.WriteYACSSchema(ScriptFile, DirName, MeshFile)``: writes a schema YACS relative to the case (see :doc:`tui_create_case`)
- The initial mesh is not published in SMESH.
- The refinement of 2D meshes with quadrangles is improved to take into account 'staircases'.
- The fields overs elements can be interpolated with two ways: intensive or extensive.
- The time steps for the interpolation of the files can be choosen:
- ``SetFieldInterpTimeStep(FieldInterp, TimeStep)``: defines the time step for the interpolation of the field (see :ref:`tui_create_iteration`)
- ``SetFieldInterpTimeStepRank(FieldInterp, TimeStep, Rank)``: defines the time step and the rank for the interpolation of the field (see :ref:`tui_create_iteration`)
+ ``SetFieldInterpTimeStep(FieldInterp, TimeStep)``: defines the time step for the interpolation of the field (see :doc:`tui_create_iteration`)
+ ``SetFieldInterpTimeStepRank(FieldInterp, TimeStep, Rank)``: defines the time step and the rank for the interpolation of the field (see :doc:`tui_create_iteration`)
- If error occurs into the python instructions, the program stops.
-SALOME RELEASE 7.4 :
+SALOME RELEASE 7.4:
- Scroll bar for the large windows.
- The meshes are published only if requested.
- The torus is a new possibility for the analytical boundary.
- - The function ``SetAdapRefinUnRef(TypeAdap, TypeRaff, TypeDera)`` is suppressed. For an adaptation by a field or by zones, the mode refinement/unrefinement is automatically set from the other parameters of the hypothesis. For a uniform refinement, the function ``SetUnifRefinUnRef(TypeRaffDera)`` is now used (see :ref:`tui_create_hypothese`).
+ - The function ``SetAdapRefinUnRef(TypeAdap, TypeRaff, TypeDera)`` is suppressed. For an adaptation by a field or by zones, the mode refinement/unrefinement is automatically set from the other parameters of the hypothesis. For a uniform refinement, the function ``SetUnifRefinUnRef(TypeRaffDera)`` is now used (see :doc:`tui_create_hypothese`).
- If an error occurs, the name of the message file is displayed.
- Release 11.1 for the HOMARD binary. The conformal refinement of the meshes with hexaedra is improved: less meshes are produced.
-SALOME RELEASE 7.5 :
- - Additional extra output: diameters and qualities; the function ``SetExtraOutput(option)`` is now used (see :ref:`tui_create_hypothese`).
+SALOME RELEASE 7.5:
+ - Additional extra output: diameters and qualities; the function ``SetExtraOutput(option)`` is now used (see :doc:`tui_create_hypothese`).
- Release 11.2 for the HOMARD binary ; improvement for the 1D curve boundaries.
-SALOME RELEASE 7.6 :
+SALOME RELEASE 7.6:
- Interpolation of all the fields into the input file.
-SALOME RELEASE 8.1 :
+SALOME RELEASE 8.1:
- Interface with med 3.2.0
- New options within the choice conformal/non-conformal.
- Less meshes are involved in the case of a non conformal refinement of a mesh made of hexahedra.
- Correction of an error when refinement and unrefinement are coupled for a mesh made of hexahedra.
-SALOME RELEASE 8.3 :
+SALOME RELEASE 8.3:
- Interface with med 3.2.1
- Output of the list of the meshes that are connected with a non conformal link.
+SALOME RELEASE 8.4:
+ - Interface with med 3.3.0
+ - Correction for the conformal coarsening of a set of hexaedra
+ - Correction of the interpolation od a volumic field over prisms which are splitted into tetrahedra and pyramids
+ - Integer MED fields to glue the elements with different levels when non conformal refinement is used
#
"""
Python script for HOMARD
-Copyright EDF-R&D 2014
+Copyright EDF-R&D 2014, 2017
"""
-__revision__ = "V1.3"
+__revision__ = "V2.01"
import os
import sys
Pour une référence à HOMARD, utiliser :
+G. Nicolas, T. Fouquet, S. Geniaut, S. Cuvilliez, Improved Adaptive Mesh Refinement for Conformal Hexahedral Meshes, "*Advances in Engineering Software*", Vol. 102, pp. 14-28, 2016, doi:10.1016/j.advengsoft.2016.07.014
+
+`Accessible en cliquant ici <http://dx.doi.org/10.1016/j.advengsoft.2016.07.014>`_
+
G. Nicolas and T. Fouquet, Adaptive Mesh Refinement for Conformal Hexahedral Meshes, "*Finite Elements in Analysis and Design*", Vol. 67, pp. 1-12, 2013, doi:10.1016/j.finel.2012.11.008
`Accessible en cliquant ici <http://dx.doi.org/10.1016/j.finel.2012.11.008>`_
SALOME VERSION 7.1 :
- Correction d'une anomalie sur le filtrage du raffinement par les groupes
- - Pilotage de l'adaptation par un seuil basé sur la moyenne et l'écart-type (voir :ref:`tui_create_hypothese`)
- - En TUI, choix du dernier instant comme instant de pilotage par SetTimeStepRankLast (voir :ref:`tui_create_iteration`)
+ - Pilotage de l'adaptation par un seuil basé sur la moyenne et l'écart-type (voir :doc:`tui_create_hypothese`)
+ - En TUI, choix du dernier instant comme instant de pilotage par SetTimeStepRankLast (voir :doc:`tui_create_iteration`)
- Possibilité de détruire les objets (GUI et TUI)
- Modification des fonctions TUI pour qu'elles agissent sur les objets et non plus sur les noms :
``homard.AssociateIterHypo(iter_name,hypo_name)`` devient ``iter.AssociateHypo(hypo_name)``, ``homard.AssociateHypoZone(hypo_name,zone_name,type_use)`` devient ``hypo.AddZone(zone_name,type_use)``, ``homard.CreateIteration(iter_name,iter_parent_name)`` devient ``case.NextIteration(iter_name)`` ou ``iter.NextIteration(iter_name)``
- Ajout de fonctions :
- ``cas.LastIteration()`` : retourne la dernière itération de la descendance du cas (voir :ref:`tui_create_iteration`)
- - Analyse de maillages (voir :ref:`gui_mesh_info`)
+ ``cas.LastIteration()`` : retourne la dernière itération de la descendance du cas (voir :doc:`tui_create_iteration`)
+ - Analyse de maillages (voir :doc:`gui_mesh_info`)
SALOME VERSION 7.2 :
- - Possibilité d'utiliser un cone comme support de frontiere 2D (voir :ref:`tui_create_boundary`)
- - Choix de publier ou non le résultat dans SMESH (voir :ref:`tui_create_iteration`)
+ - Possibilité d'utiliser un cone comme support de frontiere 2D (voir :doc:`tui_create_boundary`)
+ - Choix de publier ou non le résultat dans SMESH (voir :doc:`tui_create_iteration`)
``iter.Compute(option)`` devient ``iter.Compute(option1, option2)``
- - Possibilité de poursuivre une suite d'itérations archivées dans un répertoire (voir :ref:`gui_create_case`)
+ - Possibilité de poursuivre une suite d'itérations archivées dans un répertoire (voir :doc:`gui_create_case`)
SALOME VERSION 7.3 :
- Documentation en anglais
- - Options pour suivre le calcul de l'itération (voir :ref:`tui_create_iteration`) : ``iter.SetInfoCompute(MessInfo)``
+ - Options pour suivre le calcul de l'itération (voir :doc:`tui_create_iteration`) : ``iter.SetInfoCompute(MessInfo)``
- Creation automatique de schéma YACS :
- ``case.WriteYACSSchema(ScriptFile, DirName, MeshFile)`` : écrit un schéma YACS correspondant au cas (voir :ref:`tui_create_case`)
+ ``case.WriteYACSSchema(ScriptFile, DirName, MeshFile)`` : écrit un schéma YACS correspondant au cas (voir :doc:`tui_create_case`)
- Le maillage initial n'est pas publié dans SMESH.
- Le raffinement de maillages 2D avec des quadrangles est amélioré pour prendre en compte des zones 'en escalier'.
- Les champs constants par mailles peuvent être interpolés selon les deux modes : intensif ou extensif.
- Possibilité de choix des pas de temps pour l'interpolation des champs :
- ``SetFieldInterpTimeStep(FieldInterp, TimeStep)`` : interpole le champ au pas de temps TimeStep (voir :ref:`tui_create_iteration`)
- ``SetFieldInterpTimeStepRank(FieldInterp, TimeStep, Rank)`` : interpole le champ au pas de temps TimeStep et au numéro d'ordre Rank (voir :ref:`tui_create_iteration`)
+ ``SetFieldInterpTimeStep(FieldInterp, TimeStep)`` : interpole le champ au pas de temps TimeStep (voir :doc:`tui_create_iteration`)
+ ``SetFieldInterpTimeStepRank(FieldInterp, TimeStep, Rank)`` : interpole le champ au pas de temps TimeStep et au numéro d'ordre Rank (voir :doc:`tui_create_iteration`)
- Arrêt en cas d'erreur dans les données des instructions python
SALOME VERSION 7.4 :
- Grandes fenêtres avec ascenseur.
- Publication des maillages à la demande.
- Suivi de frontière analytique torique.
- - La fonction ``SetAdapRefinUnRef(TypeAdap, TypeRaff, TypeDera)`` est supprimée. Pour une adaptation selon un champ ou des zones, le mode raffinement/déraffinement est automatiquement déduit des autres paramètres de l'hypothèse. Pour un raffinement uniforme, on utilisera ``SetUnifRefinUnRef(TypeRaffDera)`` (voir :ref:`tui_create_hypothese`).
+ - La fonction ``SetAdapRefinUnRef(TypeAdap, TypeRaff, TypeDera)`` est supprimée. Pour une adaptation selon un champ ou des zones, le mode raffinement/déraffinement est automatiquement déduit des autres paramètres de l'hypothèse. Pour un raffinement uniforme, on utilisera ``SetUnifRefinUnRef(TypeRaffDera)`` (voir :doc:`tui_create_hypothese`).
- En cas d'erreur dans l'adaptation, le nom du fichier de messages est affiché.
- Passage à la version 11.1 de l'exécutable HOMARD ; le raffinement conforme des maillages en hexaèdres est amélioré pour diminuer le nombre de mailles produites.
SALOME VERSION 7.5 :
- - Sorties supplémentaires optionnelles : diamètre et qualité ; on utilisera ``SetExtraOutput(option)`` (voir :ref:`tui_create_hypothese`).
+ - Sorties supplémentaires optionnelles : diamètre et qualité ; on utilisera ``SetExtraOutput(option)`` (voir :doc:`tui_create_hypothese`).
- Passage à la version 11.2 de l'exécutable HOMARD ; amélioration du suivi des frontières courbes 1D.
SALOME VERSION 7.6 :
- Interfaçage avec med 3.2.1
- Sortie des listes de mailles volumiques raccordées par non-conformité
+SALOME VERSION 8.4 :
+ - Interfaçage avec med 3.3.0
+ - Correction pour le déraffinement conforme d’un ensemble d’hexaèdres.
+ - Correction de l’interpolation d’un champ volumique sur les prismes découpés en tétraèdres et en pyramides.
+ - Utilisation de champs MED entiers pour les recollements non conformes de niveaux différents.
HOMARD_I
HOMARDGUI
HOMARD_SWIG
+ FrontTrack
+ FrontTrack_SWIG
tests
)
--- /dev/null
+# --- options ---
+# additional include directories
+INCLUDE_DIRECTORIES(
+ ${CAS_INCLUDE_DIRS}
+ ${Boost_INCLUDE_DIRS}
+ ${GEOM_INCLUDE_DIRS}
+ ${MEDFILE_INCLUDE_DIRS}
+ ${MEDCOUPLING_INCLUDE_DIRS}
+ ${TBB_INCLUDE_DIRS}
+)
+
+# additional preprocessor / compiler flags
+ADD_DEFINITIONS(
+ ${CAS_DEFINITIONS}
+ ${BOOST_DEFINITIONS}
+)
+
+IF(FRONTTRACK_USE_TBB)
+ SET(TBB_LIBS ${TBB_LIBRARIES})
+ENDIF(FRONTTRACK_USE_TBB)
+
+# libraries to link to
+SET(_link_LIBRARIES
+ ${CAS_TKShHealing}
+ ${CAS_TKernel}
+ ${CAS_TKBRep}
+ ${CAS_TKG3d}
+ ${CAS_TKTopAlgo}
+ ${CAS_TKGeomBase}
+ ${CAS_TKGeomAlgo}
+ ${GEOM_XAO}
+ ${MEDCoupling_medloader}
+ ${TBB_LIBS}
+ ${Boost_LIBRARIES}
+)
+
+# --- headers ---
+
+# header files
+SET(FRONTTRACK_HEADERS
+ FrontTrack.hxx
+)
+
+# --- sources ---
+
+# sources / static
+SET(FRONTTRACK_SOURCES
+ FrontTrack.cxx
+ FrontTrack_NodeGroups.cxx
+ FrontTrack_NodesOnGeom.cxx
+ FrontTrack_Projector.cxx
+ FrontTrack_Utils.cxx
+)
+
+# --- rules ---
+
+ADD_LIBRARY(FrontTrack ${FRONTTRACK_SOURCES})
+TARGET_LINK_LIBRARIES(FrontTrack ${_link_LIBRARIES} )
+INSTALL(TARGETS FrontTrack EXPORT ${PROJECT_NAME}TargetGroup DESTINATION ${SALOME_INSTALL_LIBS})
+
+INSTALL(FILES ${FRONTTRACK_HEADERS} DESTINATION ${SALOME_INSTALL_HEADERS})
--- /dev/null
+// 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] theNodeFiles - an array of names of files describing groups of nodes that
+ * will be moved onto the geometry
+ * \param [in] theXaoFileName - a path to a file in XAO format holding the geometry and
+ * the geometrical groups.
+ * \param [in] theIsParallel - if \c true, all processors are used to treat boundary shapes
+ * in parallel.
+ */
+void FrontTrack::track( const std::string& theInputMedFile,
+ const std::string& theOutputMedFile,
+ const std::vector< std::string > & theNodeFiles,
+ const std::string& theXaoFileName,
+ bool theIsParallel )
+{
+ // check arguments
+
+ if ( theNodeFiles.empty() )
+ return;
+
+#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 );
+
+ for ( size_t i = 0; i < theNodeFiles.size(); ++i )
+ {
+#ifdef _DEBUG_
+ std::cout << "Input node file:" << theNodeFiles[i] << std::endl;
+#endif
+ if ( !FT_Utils::fileExists( theNodeFiles[i] ))
+ throw std::invalid_argument( "Input node file does not exist: " + theNodeFiles[i] );
+ }
+
+#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");
+}
--- /dev/null
+// 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] inputMedFile - a MED file holding a mesh including nodes that will be
+ * moved onto the geometry
+ * \param [in] outputMedFile - a MED file to create, that will hold a modified mesh
+ * \param [in] nodeFiles - an array of names of files describing groups of nodes that
+ * will be moved onto the geometry
+ * \param [in] xaoFileName - a path to a file in XAO format holding the geometry and
+ * the geometrical groups.
+ * \param [in] isParallel - if \c true, all processors are used to treat boundary shapes
+ * in parallel.
+ */
+ void track( const std::string& inputMedFile,
+ const std::string& outputMedFile,
+ const std::vector< std::string > & nodeFiles,
+ const std::string& xaoFileName,
+ bool isParallel=true);
+
+};
+
+#endif
--- /dev/null
+// 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;
+ }
+}
--- /dev/null
+// 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
--- /dev/null
+// 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 "fr2D.**"
+ int dimPos = theNodeFile.size() - 5;
+ if ( theNodeFile[ dimPos ] == '2' )
+ _shapeDim = 2;
+ else if ( theNodeFile[ dimPos ] == '1' )
+ _shapeDim = 1;
+ else
+ throw std::invalid_argument( "Can't define dimension by node file name " + theNodeFile );
+
+ // -------------------------------------
+ // read geom group names; several lines
+ // -------------------------------------
+
+ std::vector< std::string > geomNames;
+
+ const int maxLineLen = 256;
+ char line[ maxLineLen ];
+
+ long int pos = ::ftell( file );
+ while ( ::fgets( line, maxLineLen, file )) // read a line
+ {
+ if ( ::feof( file ))
+ return; // no nodes in the file
+
+ // check if the line describes node ids in format 3I10 (e.g. " 120 1 43\n")
+ size_t lineLen = strlen( line );
+ if ( lineLen >= 31 &&
+ ::isdigit( line[9] ) &&
+ line[10] == ' ' &&
+ ::isdigit( line[19] ) &&
+ line[20] == ' ' &&
+ ::isdigit( line[29] ) &&
+ ::isspace( line[30] ))
+ break;
+
+ geomNames.push_back( line + 1 ); // skip the 1st white space
+
+ pos = ::ftell( file ); // remember the position to return if the next line holds node ids
+ }
+
+ ::fseek( file, pos, SEEK_SET ); // return to the 1st line holding nodes ids
+
+
+ // --------------
+ // read node ids
+ // --------------
+
+ FT_NodeToMove nodeIds;
+ std::vector< int > ids;
+
+ const int nbNodes = theNodeCoords->getNumberOfTuples(); // to check validity of node IDs
+
+ while ( ::fgets( line, maxLineLen, file )) // read a line
+ {
+ // find node ids in the line
+
+ char *beg = line, *end = 0;
+ long int id;
+
+ ids.clear();
+ while (( id = ::strtol( beg, &end, 10 )) &&
+ ( beg != end ))
+ {
+ ids.push_back( id );
+ if ( id > nbNodes )
+ throw std::invalid_argument( "Too large node ID: " + FT_Utils::toStr( id ));
+ beg = end;
+ }
+
+ if ( ids.size() >= 3 )
+ {
+ std::vector< int >::iterator i = ids.begin();
+ nodeIds._nodeToMove = *i;
+ nodeIds._neighborNodes.assign( ++i, ids.end() );
+
+ _nodes.push_back( nodeIds );
+ }
+
+ if ( ::feof( file ))
+ break;
+ }
+
+ // -----------------------------------------------------------------
+ // try to find FT_Projector's to boundary sub-shapes by group names
+ // -----------------------------------------------------------------
+
+ _allProjectors = & theAllProjectorsByDim[ _shapeDim - 1 ];
+
+ _projectors.reserve( geomNames.size() );
+ std::vector< const FT_Projector* > projectors;
+
+ for ( size_t i = 0; i < geomNames.size(); ++i )
+ {
+ std::string & groupName = geomNames[i];
+
+ // remove trailing white spaces
+ for ( int iC = groupName.size() - 1; iC >= 0; --iC )
+ {
+ if ( ::isspace( groupName[iC] ) )
+ groupName.resize( iC );
+ else
+ break;
+ }
+ if ( groupName.empty() )
+ continue;
+
+ _groupNames.push_back( groupName ); // keep _groupNames for easier debug :)
+
+ // get projectors by group name
+ theXaoGroups.getProjectors( groupName, _shapeDim,
+ theAllProjectorsByDim[ _shapeDim-1 ], projectors );
+ }
+
+ // ------------------------------
+ // check the found FT_Projector's
+ // ------------------------------
+
+ if ( projectors.size() == 1 )
+ {
+ _projectors.push_back( *projectors[ 0 ]);
+ }
+ else
+ {
+ Bnd_Box nodesBox;
+ for ( size_t i = 0; i < _nodes.size(); ++i )
+ nodesBox.Add( getPoint( _nodes[i]._nodeToMove ));
+
+ if ( projectors.size() > 1 )
+ {
+ // more than one boundary shape;
+ // try to filter off unnecessary projectors using a bounding box of nodes
+ for ( size_t i = 0; i < projectors.size(); ++i )
+ if ( !nodesBox.IsOut( projectors[ i ]->getBoundingBox() ))
+ _projectors.push_back( *projectors[ i ]);
+ }
+
+ if ( _projectors.empty() )
+ {
+ // select projectors using a bounding box of nodes
+ std::vector< FT_Projector > & allProjectors = *_allProjectors;
+ for ( size_t i = 0; i < allProjectors.size(); ++i )
+ if ( !nodesBox.IsOut( allProjectors[ i ].getBoundingBox() ))
+ _projectors.push_back( allProjectors[ i ]);
+
+ if ( _projectors.empty() && !_nodes.empty() )
+ throw std::runtime_error("No boundary shape found for nodes in file " + theNodeFile );
+ }
+ }
+
+ // prepare for projection - create real projectors
+ for ( size_t i = 0; i < _projectors.size(); ++i )
+ _projectors[ i ].prepareForProjection();
+
+}
+
+//================================================================================
+/*!
+ * \brief Project nodes to the shapes and move them to new positions
+ */
+//================================================================================
+
+void FT_NodesOnGeom::projectAndMove()
+{
+ _OK = true;
+
+ // check if all the shapes are planar
+ bool isAllPlanar = true;
+ for ( size_t i = 0; i < _projectors.size() && isAllPlanar; ++i )
+ isAllPlanar = _projectors[i].isPlanarBoundary();
+ if ( isAllPlanar )
+ return;
+
+ // set nodes in the order suitable for optimal projection
+ putNodesInOrder();
+
+ // project and move nodes
+
+ std::vector< FT_NodeToMove* > notProjectedNodes;
+ size_t iP, iProjector;
+ gp_Pnt newXyz;
+
+#ifdef _DEBUG_
+ std::cout << ".. _projectors.size() = " << _projectors.size() << std::endl;
+ std::cout << ".. _nodesOrder.size() = " << _nodesOrder.size() << std::endl;
+#endif
+ if ( _projectors.size() > 1 )
+ {
+ // the nodes are to be projected onto several boundary shapes;
+ // in addition to the projecting, classification on a shape is necessary
+ // in order to find out on which of the shapes a node is to be projected
+
+ iProjector = 0;
+ for ( size_t i = 0; i < _nodesOrder.size(); ++i )
+ {
+ FT_NodeToMove& nn = _nodes[ _nodesOrder[ i ]];
+ gp_Pnt xyz = getPoint( nn._nodeToMove );
+ gp_Pnt xyz1 = getPoint( nn._neighborNodes[0] );
+ gp_Pnt xyz2 = getPoint( nn._neighborNodes[1] );
+ double maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.;
+ if ( _projectors[ iProjector ].projectAndClassify( xyz, maxDist2, newXyz,
+ nn._params, nn._nearParams ))
+ {
+ moveNode( nn._nodeToMove, newXyz );
+ }
+ else // a node is not on iProjector-th shape, find the shape it is on
+ {
+ for ( iP = 1; iP < _projectors.size(); ++iP ) // check _projectors other than iProjector
+ {
+ iProjector = ( iProjector + 1 ) % _projectors.size();
+ if ( _projectors[ iProjector ].projectAndClassify( xyz, maxDist2, newXyz,
+ nn._params, nn._nearParams ))
+ {
+ moveNode( nn._nodeToMove, newXyz );
+ break;
+ }
+ }
+ if ( iP == _projectors.size() )
+ {
+ notProjectedNodes.push_back( &nn );
+
+#ifdef _DEBUG_
+ std::cerr << "Warning: no shape found for node " << nn._nodeToMove << std::endl;
+ if ( !_groupNames.empty() )
+ std::cerr << "Warning: group -- " << _groupNames[0] << std::endl;
+#endif
+ }
+ }
+ }
+ }
+ else // one shape
+ {
+ for ( size_t i = 0; i < _nodesOrder.size(); ++i )
+ {
+ FT_NodeToMove& nn = _nodes[ _nodesOrder[ i ]];
+ gp_Pnt xyz = getPoint( nn._nodeToMove );
+ gp_Pnt xyz1 = getPoint( nn._neighborNodes[0] );
+ gp_Pnt xyz2 = getPoint( nn._neighborNodes[1] );
+
+ // maxDist2 : le quart du carré de la distance entre les deux voisins du noued à bouger
+ double maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.;
+ if ( _projectors[ 0 ].project( xyz, maxDist2, newXyz,
+ nn._params, nn._nearParams ))
+ moveNode( nn._nodeToMove, newXyz );
+ else
+ notProjectedNodes.push_back( &nn );
+ }
+ }
+
+
+ if ( !notProjectedNodes.empty() )
+ {
+ // project nodes that are not projected by any of _projectors;
+ // a proper projector is selected by evaluation of a distance between neighbor nodes
+ // and a shape
+
+ std::vector< FT_Projector > & projectors = *_allProjectors;
+
+ iProjector = 0;
+ for ( size_t i = 0; i < notProjectedNodes.size(); ++i )
+ {
+ FT_NodeToMove& nn = *notProjectedNodes[ i ];
+ gp_Pnt xyz = getPoint( nn._nodeToMove );
+ gp_Pnt xyz1 = getPoint( nn._neighborNodes[0] );
+ gp_Pnt xyz2 = getPoint( nn._neighborNodes[1] );
+ double maxDist2 = xyz1.SquareDistance( xyz2 ) / 4.;
+ double tol2 = 1e-6 * maxDist2;
+
+ bool ok;
+ for ( iP = 0; iP < projectors.size(); ++iP )
+ {
+ projectors[ iProjector ].prepareForProjection();
+ 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] );
+}
--- /dev/null
+// 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
--- /dev/null
+// 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)
+ {
+ gp_Pnt proj;
+ Standard_Real param;
+
+ if ( prevSolution )
+ {
+ _dist = _projector.NextProject( prevSolution[0], _curve, P, _tol, proj, param );
+ }
+ else
+ {
+ _dist = _projector.Project( _curve, P, _tol, proj, param, false );
+ }
+ proj = _curve.Value( param );
+
+ newSolution[0] = param;
+
+ return proj;
+ }
+
+ //-----------------------------------------------------------------------------
+ // project a point to a curve and check if the projection is within the curve boundary
+ virtual bool projectAndClassify( const gp_Pnt& point,
+ const double maxDist2,
+ gp_Pnt& projection,
+ double* newSolution,
+ const double* prevSolution = 0)
+ {
+ projection = project( point, newSolution, prevSolution );
+ return ( _uRange[0] < newSolution[0] && newSolution[0] < _uRange[1] &&
+ _dist * _dist < maxDist2 );
+ }
+
+ //-----------------------------------------------------------------------------
+ // return true if a previously found solution can be used to speed up the projection
+ virtual bool canUsePrevSolution() const { return true; }
+ };
+
+ //================================================================================
+ /*!
+ * \brief Projector to a straight curve. Don't project, classify only
+ */
+ //================================================================================
+
+ struct LineProjector : public FT_RealProjector
+ {
+ gp_Pnt _p0, _p1;
+
+ //-----------------------------------------------------------------------------
+ LineProjector( TopoDS_Edge e )
+ {
+ e.Orientation( TopAbs_FORWARD );
+ _p0 = BRep_Tool::Pnt( TopExp::FirstVertex( e ));
+ _p1 = BRep_Tool::Pnt( TopExp::LastVertex ( e ));
+ }
+
+ //-----------------------------------------------------------------------------
+ // does nothing
+ virtual gp_Pnt project( const gp_Pnt& P,
+ double* newSolution,
+ const double* prevSolution = 0)
+ {
+ return P;
+ }
+ //-----------------------------------------------------------------------------
+ // check if a point lies within the line segment
+ virtual bool projectAndClassify( const gp_Pnt& point,
+ const double maxDist2,
+ gp_Pnt& projection,
+ double* newSolution,
+ const double* prevSolution = 0)
+ {
+ gp_Vec edge( _p0, _p1 );
+ gp_Vec p0p ( _p0, point );
+ double u = ( edge * p0p ) / edge.SquareMagnitude(); // param [0,1] on the edge
+ projection = ( 1. - u ) * _p0.XYZ() + u * _p1.XYZ(); // projection of the point on the edge
+ if ( u < 0 || 1 < u )
+ return false;
+
+ // check distance
+ return point.SquareDistance( projection ) < theEPS * theEPS;
+ }
+ };
+
+ //================================================================================
+ /*!
+ * \brief Projector to a circular edge
+ */
+ //================================================================================
+
+ struct CircleProjector : public FT_RealProjector
+ {
+ gp_Circ _circle;
+ double _uRange[2];
+
+ //-----------------------------------------------------------------------------
+ CircleProjector( const gp_Circ& c, const double f, const double l ):
+ _circle( c )
+ {
+ _uRange[0] = f;
+ _uRange[1] = l;
+ }
+
+ //-----------------------------------------------------------------------------
+ // project a point to the circle
+ virtual gp_Pnt project( const gp_Pnt& P,
+ double* newSolution,
+ const double* prevSolution = 0)
+ {
+ // assume that P is already on the the plane of circle, since
+ // it is in the middle of two points lying on the circle
+
+ // move P to the circle
+ const gp_Pnt& O = _circle.Location();
+ gp_Vec radiusVec( O, P );
+ double radius = radiusVec.Magnitude();
+ if ( radius < std::numeric_limits<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() );
+}
--- /dev/null
+// 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
--- /dev/null
+// 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;
+}
--- /dev/null
+// 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
--- /dev/null
+
+INCLUDE(${SWIG_USE_FILE})
+
+ADD_DEFINITIONS(${PYTHON_DEFINITIONS})
+
+SET_SOURCE_FILES_PROPERTIES(FrontTrack.i PROPERTIES CPLUSPLUS ON)
+SET_SOURCE_FILES_PROPERTIES(FrontTrack.i PROPERTIES SWIG_DEFINITIONS "-shadow")
+
+INCLUDE_DIRECTORIES(
+ ${PYTHON_INCLUDE_DIRS}
+ ${PTHREAD_INCLUDE_DIR} # pthread dependancy due to python2.7 library
+ ${CMAKE_CURRENT_SOURCE_DIR}
+ ${CMAKE_CURRENT_SOURCE_DIR}/../FrontTrack
+ )
+
+SWIG_ADD_MODULE(FrontTrack python FrontTrack.i)
+SWIG_LINK_LIBRARIES(FrontTrack ${PYTHON_LIBRARIES} ${PLATFORM_LIBS} FrontTrack)
+
+IF(WIN32)
+ SET_TARGET_PROPERTIES(_FrontTrack PROPERTIES DEBUG_OUTPUT_NAME _FrontTrack_d)
+ENDIF(WIN32)
+INSTALL(TARGETS ${SWIG_MODULE_FrontTrack_REAL_NAME} DESTINATION ${SALOME_INSTALL_PYTHON})
+
+SALOME_INSTALL_SCRIPTS(${CMAKE_CURRENT_BINARY_DIR}/FrontTrack.py ${SALOME_INSTALL_PYTHON})
--- /dev/null
+// 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"
<source>Output of the diameters</source>
<translation>直径の出力</translation>
</message>
+ <message>
+ <source>Output of the parents</source>
+ <translation>親直径の出力</translation>
+ </message>
+ <message>
+ <source>Output of the neighbours</source>
+ <translation>隣の直径の出力</translation>
+ </message>
<message>
<source>Create an iteration</source>
<translation>イテレーションの作成</translation>
Python script for HOMARD
Test test_1
"""
-__revision__ = "V3.1"
+__revision__ = "V4.01"
#========================================================================
TEST_NAME = "test_1"
laux = zones_1_et_2.GetZones()
nbzone = len(laux) // 2
jaux = 0
- for _ in range(nbzone) :
+ for iaux in range(nbzone) :
print(hyponame_2, " : ", dico[laux[jaux+1]], "sur la zone", laux[jaux])
jaux += 2
print(hyponame_2, " : champ utilisé :", zones_1_et_2.GetFieldName())
if ERROR :
raise Exception('Pb in homard_exec at iteration %d' %ERROR )
except Exception as eee:
- raise Exception('Pb in homard_exec: ' + str(eee))
+ raise Exception('Pb in homard_exec: '+str(eee.message))
#
# Test of the results
#
Python script for HOMARD
Test test_2
"""
-__revision__ = "V3.1"
+__revision__ = "V4.01"
#========================================================================
TEST_NAME = "test_2"
Python script for HOMARD
Test test_3
"""
-__revision__ = "V3.1"
+__revision__ = "V4.01"
#========================================================================
TEST_NAME = "test_3"
Python script for HOMARD
Test test_4
"""
-__revision__ = "V2.2"
+__revision__ = "V3.01"
#========================================================================
TEST_NAME = "test_4"
if ERROR :
raise Exception('Pb in geom_smesh_exec')
except Exception as eee:
- raise Exception('Pb in geom_smesh_exec: ' + str(eee))
+ raise Exception('Pb in geom_smesh_exec: '+str(eee.message))
HOMARD = salome.lcc.FindOrLoadComponent('FactoryServer', 'HOMARD')
assert HOMARD is not None, "Impossible to load homard engine"
if ERROR :
raise Exception('Pb in homard_exec at iteration %d' %ERROR )
except Exception as eee:
- raise Exception('Pb in homard_exec: ' + str(eee))
+ raise Exception('Pb in homard_exec: '+str(eee.message))
#
# Test of the results
#
Specific conditions for Code_Saturne
Test test_5
"""
-__revision__ = "V1.0"
+__revision__ = "V2.01"
#========================================================================
TEST_NAME = "test_5"
if ERROR :
raise Exception('Pb in mesh_exec')
except Exception as eee:
- raise Exception('Pb in mesh_exec: ' + str(eee))
+ raise Exception('Pb in mesh_exec: '+str(eee.message))
HOMARD = salome.lcc.FindOrLoadComponent('FactoryServer', 'HOMARD')
assert HOMARD is not None, "Impossible to load homard engine"
if ERROR :
raise Exception('Pb in homard_exec at iteration %d' %ERROR )
except Exception as eee:
- raise Exception('Pb in homard_exec: ' + str(eee))
+ raise Exception('Pb in homard_exec: '+str(eee.message))
#
# Test of the results
#
# -*- coding: utf-8 -*-
-# Copyright (C) 2011-2016 CEA/DEN, EDF R&D
+# Copyright (C) 2011-2017 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
Python script for HOMARD
Utilitaires pour les tests
"""
-__revision__ = "V1.4"
+__revision__ = "V3.01"
import os
+import MEDLoader as ml
+import shutil
#========================================================================
#========================================================================
def remove_dir(directory) :
#
#========================================================================
#========================================================================
+#
+def saveGeometry( xao_file, name, author="" ):
+ """
+ Save the geometry in a XAO file
+ """
+ import salome
+ from salome.geom import geomBuilder
+ geompy = geomBuilder.New(salome.myStudy)
+ error = 0
+
+ # find an object having groups in GEOM component
+ component = salome.myStudy.FindObjectByPath("/Geometry")
+ it = salome.myStudy.NewChildIterator( component )
+ geomObj = None
+ l_groups = list()
+ while it.More():
+ so = it.Value()
+ it.Next()
+ go = so.GetObject()
+ #
+ if go:
+ if ( go.GetName() == name ) :
+ subIt = salome.myStudy.NewChildIterator( so )
+ while subIt.More():
+ subSO = subIt.Value()
+ subIt.Next()
+ if not subSO.GetName(): continue # a reference
+ gr = subSO.GetObject()
+ if gr and geompy.ShapeIdToType( gr.GetType() ) == "GROUP":
+ l_groups.append( gr )
+ geomObj = go
+ break
+
+ if not geomObj:
+ raise RuntimeError("Cant find a geometry object in the SALOME study with name = '%s'" % name)
+
+ # save the geom object in a XAO file
+ l_fields = list()
+ ok = geompy.ExportXAO( geomObj, l_groups, l_fields, author, xao_file, "" )
+ error = not ok
+
+ return error
+#
+#========================================================================
+#========================================================================
+#
+def repositionnement (rep_calc, fic_med_brut, fic_med_new, xao_file, verbose=False) :
+#
+ """
+Pilote le repositionnement des noeuds qui ont bougé
+Entrées :
+ rep_calc : répertoire du calcul HOMARD qui est à traiter
+ fic_med_brut : fichier MED du calcul avec les coordonnées avant projection
+ fic_med_new : fichier MED du calcul avec les coordonnées après projection
+ xao_file : fichier XAO de la géométrie
+ """
+ if verbose :
+ ligne = "rep_calc = %s" % rep_calc
+ ligne += "\nfic_med_brut = %s" % fic_med_brut
+ ligne += "\nfic_med_new = %s" % fic_med_new
+ ligne += "\nxao_file = %s" % xao_file
+ print(ligne)
+
+ message = ""
+ erreur = 0
+ while not erreur :
+#
+# 1. l_fr = liste des fichiers des lignes/surfaces a suivre
+# Les fichiers des numéros de groupes par frontière sont renommés selon le support
+# à condition de ne pas être vide.
+#
+ fic_hom_med = None
+ laux = os.listdir(rep_calc)
+ #print laux
+ l_fr = list()
+ icpt_1D = 0
+ icpt_2D = 0
+ for fic in laux :
+ #print "\t" + fic
+ if ( fic[:5] == 'fort.' ) :
+ fic_fort = os.path.join(rep_calc, fic)
+ fichier = open (fic_fort, "r")
+ les_lignes = fichier.readlines()
+ fichier.close()
+ os.remove(fic_fort)
+ a_faire = False
+ for ligne in les_lignes[1:] :
+ laux1 = ligne.split()
+ if ( len(laux1) >= 3 ) :
+ a_faire = True
+ break
+ if a_faire :
+ if ( "1D" in les_lignes[0] ) :
+ nomfic_bis = "fr1D.%02d" % icpt_1D
+ icpt_1D += 1
+ else :
+ nomfic_bis = "fr2D.%02d" % icpt_2D
+ icpt_2D += 1
+ fic_1 = os.path.join(rep_calc, nomfic_bis)
+ fichier = open (fic_1, "w")
+ for ligne in les_lignes[1:] :
+ if ( ( "1D" not in ligne ) and ( "2D" not in ligne ) ) :
+ fichier.write(ligne)
+ fichier.close()
+ #print "\t\tajout de %s" % fic_1
+ l_fr.append(fic_1)
+ elif ( fic[-4:] == '.med' ) :
+ fic_hom_med = os.path.join(rep_calc, fic)
+ #print "\t\treperage de fic_hom_med =", fic_hom_med
+#
+# 2. Lancement du post-traitement si des noeuds sont concernés
+#
+ if l_fr :
+#
+ if verbose :
+ print("l_fr =", l_fr)
+ print("fic_hom_med =", fic_hom_med)
+#
+# 2.1. La fonction :
+# . prend le maillage brut dans le fichier fic_med_brut
+# . prend la liste des noeuds à bouger pour chaque groupe concerné
+# . prend la géométrie dans le fichier xao_file
+# . retourne le maillage transformé dans le fichier fic_med_new
+#
+ from FrontTrack import FrontTrack
+ ft = FrontTrack()
+ ft.track( fic_med_brut, fic_med_new, l_fr, xao_file )
+#
+# 2.2. Transfert des coordonnées dans le fichier HOMARD MED
+#
+ #if not fic_hom_med :
+ #message = "Impossible de trouver le fichier HOMARD MED dans %s" % rep_calc
+ #erreur = 12
+ #break
+##
+ #fic_coords = os.path.join(rep_calc, "coords")
+ #erreur, message = change_coords (fic_med_new, fic_coords, verbose)
+ ##erreur, message = change_coords_0 (fic_med_new, fic_hom_med, verbose)
+ #if erreur :
+ #break
+#
+# 2.3. Ménage de l'ancien fichier MED
+#
+ if ( fic_med_brut != fic_med_new ) :
+ print("Suppression du fichier %s" % fic_med_new)
+ os.remove(fic_med_brut)
+#
+# 3. Renommage du fichier si aucun noeud n'est concerné
+#
+ else :
+#
+ if ( fic_med_brut != fic_med_new ) :
+ os.rename(fic_med_brut, fic_med_new)
+#
+# 4. Mise à jour des coordonnées dans le fichier historique HOMARD/MED si des noeuds sont concernés
+#
+ if l_fr :
+#
+ erreur, message = maj_coords (rep_calc, fic_med_new, verbose)
+#
+ break
+#
+ return erreur, message
+#
+#========================================================================
+#========================================================================
+#
+def maj_coords (rep_calc, fic_med_calc, verbose=False) :
+#
+ """
+Met à jour les coordonnées du fichier de calcul vers le fichier HOMARD
+Entrées :
+ rep_calc : répertoire du calcul HOMARD qui est à traiter
+ fic_med_calc : fichier MED du calcul avec les coordonnées après projection
+ xao_file : fichier XAO de la géométrie
+ """
+ if verbose :
+ ligne = "rep_calc = %s" % rep_calc
+ ligne += "\nfic_med_calc = %s" % fic_med_calc
+ print(ligne)
+
+ message = ""
+ erreur = 0
+ while not erreur :
+#
+# 1. Recherche des inforamtions permanentes dans le fichier de configuration
+#
+ fic_conf = os.path.join(rep_calc, "HOMARD.Configuration")
+ fichier = open (fic_conf, "r")
+ les_lignes = fichier.readlines()
+ fichier.close()
+#
+ ligne0 = ""
+ icpt = 0
+ for ligne in les_lignes :
+ #print "\t" + ligne
+ for saux in ( "HOMai", "CCNoM" ) :
+ if ( saux+"NP1" in ligne ) :
+ iaux = ligne.index(saux)
+ ligne0 += saux + "N__" + ligne[iaux+8:]
+ icpt += 1
+ if ( "NumeIter" in ligne ) :
+ iaux = ligne.index("NumeIter")
+ saux = ligne[iaux+8:-1]
+ iaux = int(saux) + 1
+ s_iter = "%02d" % iaux
+ ligne0 += "NumeIter %s\n" % s_iter
+ icpt += 1
+#
+ if ( icpt != 3 ) :
+ message = "Erreur dans le décodage de %s\n" % fic_conf
+ message+= ligne0
+ erreur = 1
+ break
+#
+# 2. Création du fichier de configuration
+#
+ fic_conf_sv = os.path.join(rep_calc, "HOMARD.Configuration.majc")
+ fichier = open (fic_conf_sv, "w")
+#
+ ligne = ligne0
+ ligne += "ModeHOMA 5\n"
+ fic = os.path.join(rep_calc, "Liste.%s.maj_coords.log" % s_iter)
+ ligne += "ListeStd %s\n" % fic
+ ligne += "CCMaiN__ %s\n" % fic_med_calc
+ ligne += "RepeTrav %s\n" % rep_calc
+ ligne += "RepeInfo %s\n" % rep_calc
+ ligne += "Action homa\n"
+ ligne += "CCAssoci med\n"
+ ligne += "EcriFiHO N_SANS_FRONTIERE\n"
+ ligne += "MessInfo 10\n"
+#
+ fichier.write(ligne)
+ fichier.close()
+#
+# 3. Mise à jour
+# 3.1. Détermination de l'exécutable
+#
+ if "HOMARD_REP_EXE_PRIVATE" in os.environ :
+ HOMARD_REP_EXE = os.environ["HOMARD_REP_EXE_PRIVATE"]
+ else :
+ HOMARD_REP_EXE = os.environ["HOMARD_REP_EXE"]
+#
+ if "HOMARD_EXE_PRIVATE" in os.environ :
+ HOMARD_EXE = os.environ["HOMARD_EXE_PRIVATE"]
+ else :
+ HOMARD_EXE = os.environ["HOMARD_EXE"]
+#
+ homard_exe = os.path.join(HOMARD_REP_EXE, HOMARD_EXE)
+ if verbose :
+ ligne = "homard_exe = %s" % homard_exe
+ print(ligne)
+#
+ if not os.path.isfile(homard_exe) :
+ message = "homard_exe = %s" % homard_exe
+ message += "\nCe fichier executable n'existe pas."
+ erreur = 31
+ break
+#
+# 3.2. Lancement
+#
+ fic_conf = os.path.join(rep_calc, "HOMARD.Configuration")
+ shutil.copyfile(fic_conf_sv, fic_conf)
+#
+ os.chdir(rep_calc)
+ erreur = os.system (homard_exe)
+#
+ break
+#
+ return erreur, message
+#
+#========================================================================
+#========================================================================
+#
Python script for HOMARD
Test tutorial_1 associe au tutorial 1
"""
-__revision__ = "V3.1"
+__revision__ = "V4.01"
#========================================================================
TEST_NAME = "tutorial_1"
Python script for HOMARD
Test tutorial_2 associe au tutorial 2
"""
-__revision__ = "V3.1"
+__revision__ = "V4.01"
#========================================================================
TEST_NAME = "tutorial_2"
Python script for HOMARD
Test tutorial_3 associe au tutorial 3
"""
-__revision__ = "V3.1"
+__revision__ = "V4.01"
#========================================================================
TEST_NAME = "tutorial_3"
Python script for HOMARD
Test tutorial_4 associe au tutorial 4
"""
-__revision__ = "V3.1"
+__revision__ = "V4.01"
#========================================================================
TEST_NAME = "tutorial_4"
Python script for HOMARD
Test tutorial_5 associe au tutorial 5
"""
-__revision__ = "V3.1"
+__revision__ = "V4.01"
#========================================================================
TEST_NAME = "tutorial_5"
Maillage apres adaptation
MESH
- Date de creation : mardi 3 novembre 2015 a 9 h 1 mn 54 s
+ Date de creation : mardi 12 septembre 2017 a 8 h 47 mn 39 s
Dimension : 3
Degre : 1
C'est un maillage obtenu apres 3 adaptations.
************************************************************
* Noeuds *
************************************************************
- * Nombre total * 2503 *
+ * Nombre total * 2458 *
************************************************************
************************************************************
************************************************************
* Triangles *
************************************************************
- * Nombre total * 582 *
- * . dont triangles de regions 2D * 582 *
+ * Nombre total * 576 *
+ * . dont triangles de regions 2D * 576 *
* . dont triangles de bord * 0 *
* . dont triangles internes aux volumes * 0 *
************************************************************
* . du niveau 0 * 0 *
- * . du niveau 0.5 * 558 *
+ * . du niveau 0.5 * 552 *
* . du niveau 1 * 0 *
* . du niveau 1.5 * 24 *
************************************************************
************************************************************
* Quadrangles *
************************************************************
- * Nombre total * 455 *
- * . dont quadrangles de regions 2D * 81 *
- * . dont quadrangles de bord * 374 *
+ * Nombre total * 447 *
+ * . dont quadrangles de regions 2D * 84 *
+ * . dont quadrangles de bord * 363 *
* . dont quadrangles internes aux volumes * 0 *
************************************************************
- * . du niveau 0 * 322 *
- * . du niveau 0.5 * 81 *
- * . du niveau 1 * 52 *
+ * . du niveau 0 * 327 *
+ * . du niveau 0.5 * 84 *
+ * . du niveau 1 * 36 *
************************************************************
************************************************************
************************************************************
* Hexaedres *
************************************************************
- * Nombre total * 523 *
+ * Nombre total * 482 *
************************************************************
- * . du niveau 0 * 283 *
+ * . du niveau 0 * 290 *
* . du niveau 0.5 * 0 *
- * . du niveau 1 * 240 *
+ * . du niveau 1 * 192 *
************************************************************
************************************************************
* Pyramides *
************************************************************
- * Nombre total * 3559 *
+ * Nombre total * 3565 *
************************************************************
* . du niveau 0 * 0 *
- * . du niveau 0.5 * 3239 *
+ * . du niveau 0.5 * 3245 *
* . du niveau 1 * 0 *
* . du niveau 1.5 * 320 *
************************************************************
Maillage apres adaptation
MESH
- Date de creation : mardi 19 juillet 2016 a 11 h 51 mn 14 s
+ Date de creation : jeudi 14 septembre 2017 a 8 h 44 mn 8 s
Dimension : 3
Degre : 1
C'est un maillage obtenu apres 3 adaptations.
************************************************************
* Noeuds *
************************************************************
- * Nombre total * 2171 *
+ * Nombre total * 2008 *
************************************************************
************************************************************
* Quadrangles *
************************************************************
- * Nombre total * 664 *
- * . dont quadrangles de regions 2D * 0 *
- * . dont quadrangles de bord * 664 *
+ * Nombre total * 596 *
+ * . dont quadrangles de regions 2D * 12 *
+ * . dont quadrangles de bord * 584 *
* . dont quadrangles internes aux volumes * 0 *
************************************************************
* . du niveau 0 * 0 *
- * . du niveau 1 * 488 *
+ * . du niveau 1 * 420 *
* . du niveau 2 * 176 *
************************************************************
************************************************************
* Hexaedres *
************************************************************
- * Nombre total * 1567 *
+ * Nombre total * 1420 *
************************************************************
- * . du niveau 0 * 933 *
- * . du niveau 1 * 522 *
+ * . du niveau 0 * 954 *
+ * . du niveau 1 * 354 *
* . du niveau 2 * 112 *
************************************************************
