Salome HOME
Merge branch 'master' into pre/penta18
[modules/smesh.git] / src / StdMeshers / StdMeshers_RadialPrism_3D.cxx
index b19f79aee8fc1fe57d993476923fab030db343c7..bbf339a8b124c022eda732e432df85d6ad46d563 100644 (file)
@@ -1,32 +1,36 @@
-//  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
 //
-//  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
-//  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+// Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
 //
-//  This library is free software; you can redistribute it and/or
-//  modify it under the terms of the GNU Lesser General Public
-//  License as published by the Free Software Foundation; either
-//  version 2.1 of the License.
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
 //
-//  This library is distributed in the hope that it will be useful,
-//  but WITHOUT ANY WARRANTY; without even the implied warranty of
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-//  Lesser General Public License for more details.
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
 //
-//  You should have received a copy of the GNU Lesser General Public
-//  License along with this library; if not, write to the Free Software
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
-//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
+
 //  SMESH SMESH : implementaion of SMESH idl descriptions
 // File      : StdMeshers_RadialPrism_3D.cxx
 // Module    : SMESH
 // Created   : Fri Oct 20 11:37:07 2006
 // Author    : Edward AGAPOV (eap)
 //
+
 #include "StdMeshers_RadialPrism_3D.hxx"
 
+#include <Basics_OCCTVersion.hxx>
+
 #include "StdMeshers_ProjectionUtils.hxx"
 #include "StdMeshers_NumberOfLayers.hxx"
 #include "StdMeshers_LayerDistribution.hxx"
@@ -37,7 +41,6 @@
 #include "SMESHDS_SubMesh.hxx"
 #include "SMESH_Gen.hxx"
 #include "SMESH_Mesh.hxx"
-#include "SMESH_MeshEditor.hxx"
 #include "SMESH_MesherHelper.hxx"
 #include "SMESH_subMesh.hxx"
 
 
 #include <BRepAdaptor_Curve.hxx>
 #include <BRepBuilderAPI_MakeEdge.hxx>
-#include <BRepTools.hxx>
+#include <BRep_Tool.hxx>
 #include <TopExp_Explorer.hxx>
+#include <TopTools_DataMapIteratorOfDataMapOfShapeShape.hxx>
+#include <TopTools_MapOfShape.hxx>
 #include <TopoDS.hxx>
 #include <TopoDS_Shell.hxx>
 #include <TopoDS_Solid.hxx>
 #include <gp.hxx>
 #include <gp_Pnt.hxx>
-
+#include <BRepClass3d.hxx>
 
 using namespace std;
 
 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
 
-typedef StdMeshers_ProjectionUtils TAssocTool;
+namespace ProjectionUtils = StdMeshers_ProjectionUtils;
 
 //=======================================================================
 //function : StdMeshers_RadialPrism_3D
@@ -70,7 +75,7 @@ StdMeshers_RadialPrism_3D::StdMeshers_RadialPrism_3D(int hypId, int studyId, SME
   :SMESH_3D_Algo(hypId, studyId, gen)
 {
   _name = "RadialPrism_3D";
-  _shapeType = (1 << TopAbs_SOLID);    // 1 bit per shape type
+  _shapeType = (1 << TopAbs_SOLID);     // 1 bit per shape type
 
   _compatibleHypothesis.push_back("LayerDistribution");
   _compatibleHypothesis.push_back("NumberOfLayers");
@@ -98,8 +103,8 @@ bool StdMeshers_RadialPrism_3D::CheckHypothesis(SMESH_Mesh&
 {
   // check aShape that must have 2 shells
 /*  PAL16229
-  if ( TAssocTool::Count( aShape, TopAbs_SOLID, 0 ) != 1 ||
-       TAssocTool::Count( aShape, TopAbs_SHELL, 0 ) != 2 )
+  if ( ProjectionUtils::Count( aShape, TopAbs_SOLID, 0 ) != 1 ||
+       ProjectionUtils::Count( aShape, TopAbs_SHELL, 0 ) != 2 )
   {
     aStatus = HYP_BAD_GEOMETRY;
     return false;
@@ -160,7 +165,7 @@ bool StdMeshers_RadialPrism_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& a
 
   // get 2 shells
   TopoDS_Solid solid = TopoDS::Solid( aShape );
-  TopoDS_Shell outerShell = BRepTools::OuterShell( solid );
+  TopoDS_Shell outerShell = BRepClass3d::OuterShell( solid );
   TopoDS_Shape innerShell;
   int nbShells = 0;
   for ( TopoDS_Iterator It (solid); It.More(); It.Next(), ++nbShells )
@@ -170,15 +175,42 @@ bool StdMeshers_RadialPrism_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& a
     return error(COMPERR_BAD_SHAPE, SMESH_Comment("Must be 2 shells but not ")<<nbShells);
 
   // ----------------------------------
-  // Associate subshapes of the shells
+  // Associate sub-shapes of the shells
   // ----------------------------------
 
-  TAssocTool::TShapeShapeMap shape2ShapeMap;
-  if ( !TAssocTool::FindSubShapeAssociation( outerShell, &aMesh,
-                                             innerShell, &aMesh,
-                                             shape2ShapeMap) )
+  ProjectionUtils::TShapeShapeMap shape2ShapeMaps[2];
+  bool mapOk1 = ProjectionUtils::FindSubShapeAssociation( innerShell, &aMesh,
+                                                          outerShell, &aMesh,
+                                                          shape2ShapeMaps[0]);
+  bool mapOk2 = ProjectionUtils::FindSubShapeAssociation( innerShell.Reversed(), &aMesh,
+                                                          outerShell, &aMesh,
+                                                          shape2ShapeMaps[1]);
+  if ( !mapOk1 && !mapOk2 )
     return error(COMPERR_BAD_SHAPE,"Topology of inner and outer shells seems different" );
 
+  int iMap;
+  if ( shape2ShapeMaps[0].Extent() == shape2ShapeMaps[1].Extent() )
+  {
+    // choose an assiciation by shortest distance between VERTEXes
+    double dist1 = 0, dist2 = 0;
+    TopTools_DataMapIteratorOfDataMapOfShapeShape ssIt( shape2ShapeMaps[0]._map1to2 );
+    for (; ssIt.More(); ssIt.Next() )
+    {
+      if ( ssIt.Key().ShapeType() != TopAbs_VERTEX ) continue;
+      gp_Pnt pIn   = BRep_Tool::Pnt( TopoDS::Vertex( ssIt.Key() ));
+      gp_Pnt pOut1 = BRep_Tool::Pnt( TopoDS::Vertex( ssIt.Value() ));
+      gp_Pnt pOut2 = BRep_Tool::Pnt( TopoDS::Vertex( shape2ShapeMaps[1]( ssIt.Key() )));
+      dist1 += pIn.SquareDistance( pOut1 );
+      dist2 += pIn.SquareDistance( pOut2 );
+    }
+    iMap = ( dist1 < dist2 ) ? 0 : 1;
+  }
+  else
+  {
+    iMap = ( shape2ShapeMaps[0].Extent() > shape2ShapeMaps[1].Extent() ) ? 0 : 1;
+  }
+  ProjectionUtils::TShapeShapeMap& shape2ShapeMap = shape2ShapeMaps[iMap];
+
   // ------------------
   // Make mesh
   // ------------------
@@ -188,20 +220,20 @@ bool StdMeshers_RadialPrism_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& a
 
   for ( exp.Init( outerShell, TopAbs_FACE ); exp.More(); exp.Next() )
   {
-    // Corresponding subshapes
+    // Corresponding sub-shapes
     TopoDS_Face outFace = TopoDS::Face( exp.Current() );
     TopoDS_Face inFace;
-    if ( !shape2ShapeMap.IsBound( outFace )) {
+    if ( !shape2ShapeMap.IsBound( outFace, /*isOut=*/true )) {
       return error(SMESH_Comment("Corresponding inner face not found for face #" )
                    << meshDS->ShapeToIndex( outFace ));
     } else {
-      inFace = TopoDS::Face( shape2ShapeMap( outFace ));
+      inFace = TopoDS::Face( shape2ShapeMap( outFace, /*isOut=*/true ));
     }
 
     // Find matching nodes of in and out faces
-    TNodeNodeMap nodeIn2OutMap;
-    if ( ! TAssocTool::FindMatchingNodesOnFaces( inFace, &aMesh, outFace, &aMesh,
-                                                 shape2ShapeMap, nodeIn2OutMap ))
+    ProjectionUtils::TNodeNodeMap nodeIn2OutMap;
+    if ( ! ProjectionUtils::FindMatchingNodesOnFaces( inFace, &aMesh, outFace, &aMesh,
+                                                      shape2ShapeMap, nodeIn2OutMap ))
       return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Mesh on faces #")
                    << meshDS->ShapeToIndex( outFace ) << " and "
                    << meshDS->ShapeToIndex( inFace ) << " seems different" );
@@ -304,11 +336,10 @@ public:
   static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
   {
     const int myID = -1000;
-    map < int, SMESH_1D_Algo * > & algoMap = aMesh.GetGen()->_map1D_Algo;
-    map < int, SMESH_1D_Algo * >::iterator id_algo = algoMap.find( myID );
-    if ( id_algo == algoMap.end() )
-      return new TNodeDistributor( myID, 0, aMesh.GetGen() );
-    return static_cast< TNodeDistributor* >( id_algo->second );
+    TNodeDistributor* myHyp = dynamic_cast<TNodeDistributor*>( aMesh.GetHypothesis( myID ));
+    if ( !myHyp )
+      myHyp = new TNodeDistributor( myID, 0, aMesh.GetGen() );
+    return myHyp;
   }
   // -----------------------------------------------------------------------------
   bool Compute( vector< double > &                  positions,
@@ -387,3 +418,224 @@ bool StdMeshers_RadialPrism_3D::computeLayerPositions(const gp_Pnt& pIn,
   }
   RETURN_BAD_RESULT("Bad hypothesis");
 }
+
+
+//=======================================================================
+//function : Evaluate
+//purpose  : 
+//=======================================================================
+
+bool StdMeshers_RadialPrism_3D::Evaluate(SMESH_Mesh& aMesh,
+                                         const TopoDS_Shape& aShape,
+                                         MapShapeNbElems& aResMap)
+{
+  // get 2 shells
+  TopoDS_Solid solid = TopoDS::Solid( aShape );
+  TopoDS_Shell outerShell = BRepClass3d::OuterShell( solid );
+  TopoDS_Shape innerShell;
+  int nbShells = 0;
+  for ( TopoDS_Iterator It (solid); It.More(); It.Next(), ++nbShells )
+    if ( !outerShell.IsSame( It.Value() ))
+      innerShell = It.Value();
+  if ( nbShells != 2 ) {
+    std::vector<int> aResVec(SMDSEntity_Last);
+    for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
+    SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+    aResMap.insert(std::make_pair(sm,aResVec));
+    SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+    smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
+    return false;
+  }
+
+  // Associate sub-shapes of the shells
+  ProjectionUtils::TShapeShapeMap shape2ShapeMap;
+  if ( !ProjectionUtils::FindSubShapeAssociation( outerShell, &aMesh,
+                                                  innerShell, &aMesh,
+                                                  shape2ShapeMap) ) {
+    std::vector<int> aResVec(SMDSEntity_Last);
+    for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
+    SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+    aResMap.insert(std::make_pair(sm,aResVec));
+    SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+    smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
+    return false;
+  }
+
+  // get info for outer shell
+  int nb0d_Out=0, nb2d_3_Out=0, nb2d_4_Out=0;
+  //TopTools_SequenceOfShape FacesOut;
+  for (TopExp_Explorer exp(outerShell, TopAbs_FACE); exp.More(); exp.Next()) {
+    //FacesOut.Append(exp.Current());
+    SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
+    MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
+    std::vector<int> aVec = (*anIt).second;
+    nb0d_Out += aVec[SMDSEntity_Node];
+    nb2d_3_Out += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
+    nb2d_4_Out += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
+  }
+  int nb1d_Out = 0;
+  TopTools_MapOfShape tmpMap;
+  for (TopExp_Explorer exp(outerShell, TopAbs_EDGE); exp.More(); exp.Next()) {
+    if( tmpMap.Contains( exp.Current() ) )
+      continue;
+    tmpMap.Add( exp.Current() );
+    SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
+    MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
+    std::vector<int> aVec = (*anIt).second;
+    nb0d_Out += aVec[SMDSEntity_Node];
+    nb1d_Out += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
+  }
+  tmpMap.Clear();
+  for (TopExp_Explorer exp(outerShell, TopAbs_VERTEX); exp.More(); exp.Next()) {
+    if( tmpMap.Contains( exp.Current() ) )
+      continue;
+    tmpMap.Add( exp.Current() );
+    nb0d_Out++;
+  }
+
+  // get info for inner shell
+  int nb0d_In=0, nb2d_3_In=0, nb2d_4_In=0;
+  //TopTools_SequenceOfShape FacesIn;
+  for (TopExp_Explorer exp(innerShell, TopAbs_FACE); exp.More(); exp.Next()) {
+    //FacesIn.Append(exp.Current());
+    SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
+    MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
+    std::vector<int> aVec = (*anIt).second;
+    nb0d_In += aVec[SMDSEntity_Node];
+    nb2d_3_In += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
+    nb2d_4_In += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
+  }
+  int nb1d_In = 0;
+  tmpMap.Clear();
+  bool IsQuadratic = false;
+  bool IsFirst = true;
+  for (TopExp_Explorer exp(innerShell, TopAbs_EDGE); exp.More(); exp.Next()) {
+    if( tmpMap.Contains( exp.Current() ) )
+      continue;
+    tmpMap.Add( exp.Current() );
+    SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
+    MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
+    std::vector<int> aVec = (*anIt).second;
+    nb0d_In += aVec[SMDSEntity_Node];
+    nb1d_In += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
+    if(IsFirst) {
+      IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
+      IsFirst = false;
+    }
+  }
+  tmpMap.Clear();
+  for (TopExp_Explorer exp(innerShell, TopAbs_VERTEX); exp.More(); exp.Next()) {
+    if( tmpMap.Contains( exp.Current() ) )
+      continue;
+    tmpMap.Add( exp.Current() );
+    nb0d_In++;
+  }
+
+  bool IsOK = (nb0d_Out==nb0d_In) && (nb1d_Out==nb1d_In) && 
+              (nb2d_3_Out==nb2d_3_In) && (nb2d_4_Out==nb2d_4_In);
+  if(!IsOK) {
+    std::vector<int> aResVec(SMDSEntity_Last);
+    for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
+    SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+    aResMap.insert(std::make_pair(sm,aResVec));
+    SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+    smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
+    return false;
+  }
+
+  int nbLayers = 0;
+  if( myNbLayerHypo ) {
+    nbLayers = myNbLayerHypo->GetNumberOfLayers();
+  }
+  if ( myDistributionHypo ) {
+    if ( !myDistributionHypo->GetLayerDistribution() ) {
+      std::vector<int> aResVec(SMDSEntity_Last);
+      for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
+      SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+      aResMap.insert(std::make_pair(sm,aResVec));
+      SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+      smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
+      return false;
+    }
+    TopExp_Explorer exp(outerShell, TopAbs_VERTEX);
+    TopoDS_Vertex Vout = TopoDS::Vertex(exp.Current());
+    TopoDS_Vertex Vin = TopoDS::Vertex( shape2ShapeMap(Vout) );
+    if ( myLayerPositions.empty() ) {
+      gp_Pnt pIn = BRep_Tool::Pnt(Vin);
+      gp_Pnt pOut = BRep_Tool::Pnt(Vout);
+      computeLayerPositions( pIn, pOut );
+    }
+    nbLayers = myLayerPositions.size() + 1;
+  }
+
+  std::vector<int> aResVec(SMDSEntity_Last);
+  for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
+  if(IsQuadratic) {
+    aResVec[SMDSEntity_Quad_Penta] = nb2d_3_Out * nbLayers;
+    aResVec[SMDSEntity_Quad_Hexa] = nb2d_4_Out * nbLayers;
+    int nb1d = ( nb2d_3_Out*3 + nb2d_4_Out*4 ) / 2;
+    aResVec[SMDSEntity_Node] = nb0d_Out * ( 2*nbLayers - 1 ) - nb1d * nbLayers;
+  }
+  else {
+    aResVec[SMDSEntity_Node] = nb0d_Out * ( nbLayers - 1 );
+    aResVec[SMDSEntity_Penta] = nb2d_3_Out * nbLayers;
+    aResVec[SMDSEntity_Hexa] = nb2d_4_Out * nbLayers;
+  }
+  SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+  aResMap.insert(std::make_pair(sm,aResVec));
+
+  return true;
+}
+
+//================================================================================
+/*!
+ * \brief Return true if the algorithm can mesh this shape
+ *  \param [in] aShape - shape to check
+ *  \param [in] toCheckAll - if true, this check returns OK if all shapes are OK,
+ *              else, returns OK if at least one shape is OK
+ */
+//================================================================================
+
+bool StdMeshers_RadialPrism_3D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
+{
+  int nbFoundSolids = 0;
+  for (TopExp_Explorer exp( aShape, TopAbs_SOLID ); exp.More(); exp.Next(), ++nbFoundSolids )
+  {
+    TopoDS_Shape shell[2];
+    int nbShells = 0;
+    for ( TopoDS_Iterator It (exp.Current()); It.More(); It.Next() )
+    {
+      nbShells++;
+      if ( nbShells > 2 ) {
+        if ( toCheckAll ) return false;
+        break;
+      }
+      shell[ nbShells-1 ] = It.Value();
+    }
+    if ( nbShells != 2 ) { 
+      if ( toCheckAll ) return false;  
+      continue; 
+    }
+
+    int nbFaces1 = SMESH_MesherHelper:: Count( shell[0], TopAbs_FACE, 0 );
+    int nbFaces2 = SMESH_MesherHelper:: Count( shell[1], TopAbs_FACE, 0 );
+    if ( nbFaces1 != nbFaces2 ){
+      if( toCheckAll ) return false;
+      continue;
+    }
+    int nbEdges1 = SMESH_MesherHelper:: Count( shell[0], TopAbs_EDGE, 0 );
+    int nbEdges2 = SMESH_MesherHelper:: Count( shell[1], TopAbs_EDGE, 0 );
+    if ( nbEdges1 != nbEdges2 ){
+      if( toCheckAll ) return false;
+      continue;
+    }
+    int nbVertices1 = SMESH_MesherHelper:: Count( shell[0], TopAbs_VERTEX, 0 );
+    int nbVertices2 = SMESH_MesherHelper:: Count( shell[1], TopAbs_VERTEX, 0 );
+    if ( nbVertices1 != nbVertices2 ){
+      if( toCheckAll ) return false;
+      continue;
+    }
+    if ( !toCheckAll ) return true;
+  }
+  return ( toCheckAll && nbFoundSolids != 0);
+};