src/StdMeshers/StdMeshers_RadialPrism_3D.cxx

   1 // Copyright (C) 2007-2014  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 //  SMESH SMESH : implementaion of SMESH idl descriptions
  24 // File      : StdMeshers_RadialPrism_3D.cxx
  25 // Module    : SMESH
  26 // Created   : Fri Oct 20 11:37:07 2006
  27 // Author    : Edward AGAPOV (eap)
  28 //
  29
  30 #include "StdMeshers_RadialPrism_3D.hxx"
  31
  32 #include <Basics_OCCTVersion.hxx>
  33
  34 #include "StdMeshers_ProjectionUtils.hxx"
  35 #include "StdMeshers_NumberOfLayers.hxx"
  36 #include "StdMeshers_LayerDistribution.hxx"
  37 #include "StdMeshers_Prism_3D.hxx"
  38 #include "StdMeshers_Regular_1D.hxx"
  39
  40 #include "SMDS_MeshNode.hxx"
  41 #include "SMESHDS_SubMesh.hxx"
  42 #include "SMESH_Gen.hxx"
  43 #include "SMESH_Mesh.hxx"
  44 #include "SMESH_MesherHelper.hxx"
  45 #include "SMESH_subMesh.hxx"
  46
  47 #include "utilities.h"
  48
  49 #include <BRepAdaptor_Curve.hxx>
  50 #include <BRepBuilderAPI_MakeEdge.hxx>
  51 #include <BRep_Tool.hxx>
  52 #include <TopExp_Explorer.hxx>
  53 #include <TopTools_DataMapIteratorOfDataMapOfShapeShape.hxx>
  54 #include <TopTools_MapOfShape.hxx>
  55 #include <TopoDS.hxx>
  56 #include <TopoDS_Shell.hxx>
  57 #include <TopoDS_Solid.hxx>
  58 #include <gp.hxx>
  59 #include <gp_Pnt.hxx>
  60 #if OCC_VERSION_LARGE > 0x06050400
  61 #include <BRepClass3d.hxx>
  62 #else
  63 #include <BRepTools.hxx>
  64 #endif
  65
  66 using namespace std;
  67
  68 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
  69 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
  70
  71 namespace ProjectionUtils = StdMeshers_ProjectionUtils;
  72
  73 //=======================================================================
  74 //function : StdMeshers_RadialPrism_3D
  75 //purpose  :
  76 //=======================================================================
  77
  78 StdMeshers_RadialPrism_3D::StdMeshers_RadialPrism_3D(int hypId, int studyId, SMESH_Gen* gen)
  79   :SMESH_3D_Algo(hypId, studyId, gen)
  80 {
  81   _name = "RadialPrism_3D";
  82   _shapeType = (1 << TopAbs_SOLID);     // 1 bit per shape type
  83
  84   _compatibleHypothesis.push_back("LayerDistribution");
  85   _compatibleHypothesis.push_back("NumberOfLayers");
  86   myNbLayerHypo = 0;
  87   myDistributionHypo = 0;
  88 }
  89
  90 //================================================================================
  91 /*!
  92  * \brief Destructor
  93  */
  94 //================================================================================
  95
  96 StdMeshers_RadialPrism_3D::~StdMeshers_RadialPrism_3D()
  97 {}
  98
  99 //=======================================================================
 100 //function : CheckHypothesis
 101 //purpose  :
 102 //=======================================================================
 103
 104 bool StdMeshers_RadialPrism_3D::CheckHypothesis(SMESH_Mesh&                          aMesh,
 105                                                 const TopoDS_Shape&                  aShape,
 106                                                 SMESH_Hypothesis::Hypothesis_Status& aStatus)
 107 {
 108   // check aShape that must have 2 shells
 109 /*  PAL16229
 110   if ( ProjectionUtils::Count( aShape, TopAbs_SOLID, 0 ) != 1 ||
 111        ProjectionUtils::Count( aShape, TopAbs_SHELL, 0 ) != 2 )
 112   {
 113     aStatus = HYP_BAD_GEOMETRY;
 114     return false;
 115   }
 116 */
 117   myNbLayerHypo = 0;
 118   myDistributionHypo = 0;
 119
 120   list <const SMESHDS_Hypothesis * >::const_iterator itl;
 121
 122   const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
 123   if ( hyps.size() == 0 )
 124   {
 125     aStatus = SMESH_Hypothesis::HYP_MISSING;
 126     return false;  // can't work with no hypothesis
 127   }
 128
 129   if ( hyps.size() > 1 )
 130   {
 131     aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
 132     return false;
 133   }
 134
 135   const SMESHDS_Hypothesis *theHyp = hyps.front();
 136
 137   string hypName = theHyp->GetName();
 138
 139   if (hypName == "NumberOfLayers")
 140   {
 141     myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
 142     aStatus = SMESH_Hypothesis::HYP_OK;
 143     return true;
 144   }
 145   if (hypName == "LayerDistribution")
 146   {
 147     myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
 148     aStatus = SMESH_Hypothesis::HYP_OK;
 149     return true;
 150   }
 151   aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 152   return true;
 153 }
 154
 155 //=======================================================================
 156 //function : Compute
 157 //purpose  :
 158 //=======================================================================
 159
 160 bool StdMeshers_RadialPrism_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
 161 {
 162   TopExp_Explorer exp;
 163   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 164
 165   myHelper = new SMESH_MesherHelper( aMesh );
 166   myHelper->IsQuadraticSubMesh( aShape );
 167   // to delete helper at exit from Compute()
 168   std::auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
 169
 170   // get 2 shells
 171   TopoDS_Solid solid = TopoDS::Solid( aShape );
 172 #if OCC_VERSION_LARGE > 0x06050400
 173   TopoDS_Shell outerShell = BRepClass3d::OuterShell( solid );
 174 #else
 175   TopoDS_Shell outerShell = BRepTools::OuterShell( solid );
 176 #endif
 177   TopoDS_Shape innerShell;
 178   int nbShells = 0;
 179   for ( TopoDS_Iterator It (solid); It.More(); It.Next(), ++nbShells )
 180     if ( !outerShell.IsSame( It.Value() ))
 181       innerShell = It.Value();
 182   if ( nbShells != 2 )
 183     return error(COMPERR_BAD_SHAPE, SMESH_Comment("Must be 2 shells but not ")<<nbShells);
 184
 185   // ----------------------------------
 186   // Associate sub-shapes of the shells
 187   // ----------------------------------
 188
 189   ProjectionUtils::TShapeShapeMap shape2ShapeMaps[2];
 190   if ( !ProjectionUtils::FindSubShapeAssociation( innerShell, &aMesh,
 191                                                   outerShell, &aMesh,
 192                                                   shape2ShapeMaps[0])
 193        &&
 194        !ProjectionUtils::FindSubShapeAssociation( innerShell.Reversed(), &aMesh,
 195                                                   outerShell, &aMesh,
 196                                                   shape2ShapeMaps[1]))
 197     return error(COMPERR_BAD_SHAPE,"Topology of inner and outer shells seems different" );
 198
 199   int iMap;
 200   if ( shape2ShapeMaps[0].Extent() == shape2ShapeMaps[1].Extent() )
 201   {
 202     // choose an assiciation by shortest distance between VERTEXes
 203     double dist1 = 0, dist2 = 0;
 204     TopTools_DataMapIteratorOfDataMapOfShapeShape ssIt( shape2ShapeMaps[0]._map1to2 );
 205     for (; ssIt.More(); ssIt.Next() )
 206     {
 207       if ( ssIt.Key().ShapeType() != TopAbs_VERTEX ) continue;
 208       gp_Pnt pIn   = BRep_Tool::Pnt( TopoDS::Vertex( ssIt.Key() ));
 209       gp_Pnt pOut1 = BRep_Tool::Pnt( TopoDS::Vertex( ssIt.Value() ));
 210       gp_Pnt pOut2 = BRep_Tool::Pnt( TopoDS::Vertex( shape2ShapeMaps[1]( ssIt.Key() )));
 211       dist1 += pIn.SquareDistance( pOut1 );
 212       dist2 += pIn.SquareDistance( pOut2 );
 213     }
 214     iMap = ( dist1 < dist2 ) ? 0 : 1;
 215   }
 216   else
 217   {
 218     iMap = ( shape2ShapeMaps[0].Extent() > shape2ShapeMaps[1].Extent() ) ? 0 : 1;
 219   }
 220   ProjectionUtils::TShapeShapeMap& shape2ShapeMap = shape2ShapeMaps[iMap];
 221
 222   // ------------------
 223   // Make mesh
 224   // ------------------
 225
 226   TNode2ColumnMap node2columnMap;
 227   myLayerPositions.clear();
 228
 229   for ( exp.Init( outerShell, TopAbs_FACE ); exp.More(); exp.Next() )
 230   {
 231     // Corresponding sub-shapes
 232     TopoDS_Face outFace = TopoDS::Face( exp.Current() );
 233     TopoDS_Face inFace;
 234     if ( !shape2ShapeMap.IsBound( outFace, /*isOut=*/true )) {
 235       return error(SMESH_Comment("Corresponding inner face not found for face #" )
 236                    << meshDS->ShapeToIndex( outFace ));
 237     } else {
 238       inFace = TopoDS::Face( shape2ShapeMap( outFace, /*isOut=*/true ));
 239     }
 240
 241     // Find matching nodes of in and out faces
 242     TNodeNodeMap nodeIn2OutMap;
 243     if ( ! ProjectionUtils::FindMatchingNodesOnFaces( inFace, &aMesh, outFace, &aMesh,
 244                                                       shape2ShapeMap, nodeIn2OutMap ))
 245       return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Mesh on faces #")
 246                    << meshDS->ShapeToIndex( outFace ) << " and "
 247                    << meshDS->ShapeToIndex( inFace ) << " seems different" );
 248
 249     // Create volumes
 250
 251     SMDS_ElemIteratorPtr faceIt = meshDS->MeshElements( inFace )->GetElements();
 252     while ( faceIt->more() ) // loop on faces on inFace
 253     {
 254       const SMDS_MeshElement* face = faceIt->next();
 255       if ( !face || face->GetType() != SMDSAbs_Face )
 256         continue;
 257       int nbNodes = face->NbNodes();
 258       if ( face->IsQuadratic() )
 259         nbNodes /= 2;
 260
 261       // find node columns for each node
 262       vector< const TNodeColumn* > columns( nbNodes );
 263       for ( int i = 0; i < nbNodes; ++i )
 264       {
 265         const SMDS_MeshNode* nIn = face->GetNode( i );
 266         TNode2ColumnMap::iterator n_col = node2columnMap.find( nIn );
 267         if ( n_col != node2columnMap.end() ) {
 268           columns[ i ] = & n_col->second;
 269         }
 270         else {
 271           TNodeNodeMap::iterator nInOut = nodeIn2OutMap.find( nIn );
 272           if ( nInOut == nodeIn2OutMap.end() )
 273             RETURN_BAD_RESULT("No matching node for "<< nIn->GetID() <<
 274                               " in face "<< face->GetID());
 275           columns[ i ] = makeNodeColumn( node2columnMap, nIn, nInOut->second );
 276         }
 277       }
 278
 279       StdMeshers_Prism_3D::AddPrisms( columns, myHelper );
 280     }
 281   } // loop on faces of out shell
 282
 283   return true;
 284 }
 285
 286 //================================================================================
 287 /*!
 288  * \brief Create a column of nodes from outNode to inNode
 289   * \param n2ColMap - map of node columns to add a created column
 290   * \param outNode - botton node of a column
 291   * \param inNode - top node of a column
 292   * \retval const TNodeColumn* - a new column pointer
 293  */
 294 //================================================================================
 295
 296 TNodeColumn* StdMeshers_RadialPrism_3D::makeNodeColumn( TNode2ColumnMap&     n2ColMap,
 297                                                         const SMDS_MeshNode* outNode,
 298                                                         const SMDS_MeshNode* inNode)
 299 {
 300   SMESHDS_Mesh * meshDS = myHelper->GetMeshDS();
 301   int shapeID = myHelper->GetSubShapeID();
 302
 303   if ( myLayerPositions.empty() ) {
 304     gp_Pnt pIn = gpXYZ( inNode ), pOut = gpXYZ( outNode );
 305     computeLayerPositions( pIn, pOut );
 306   }
 307   int nbSegments = myLayerPositions.size() + 1;
 308
 309   TNode2ColumnMap::iterator n_col =
 310     n2ColMap.insert( make_pair( outNode, TNodeColumn() )).first;
 311   TNodeColumn & column = n_col->second;
 312   column.resize( nbSegments + 1 );
 313   column.front() = outNode;
 314   column.back() =  inNode;
 315
 316   gp_XYZ p1 = gpXYZ( outNode );
 317   gp_XYZ p2 = gpXYZ( inNode );
 318   for ( int z = 1; z < nbSegments; ++z )
 319   {
 320     double r = myLayerPositions[ z - 1 ];
 321     gp_XYZ p = ( 1 - r ) * p1 + r * p2;
 322     SMDS_MeshNode* n = meshDS->AddNode( p.X(), p.Y(), p.Z() );
 323     meshDS->SetNodeInVolume( n, shapeID );
 324     column[ z ] = n;
 325   }
 326
 327   return & column;
 328 }
 329
 330 //================================================================================
 331 //================================================================================
 332 /*!
 333  * \brief Class computing layers distribution using data of
 334  *        StdMeshers_LayerDistribution hypothesis
 335  */
 336 //================================================================================
 337 //================================================================================
 338
 339 class TNodeDistributor: public StdMeshers_Regular_1D
 340 {
 341   list <const SMESHDS_Hypothesis *> myUsedHyps;
 342 public:
 343   // -----------------------------------------------------------------------------
 344   static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
 345   {
 346     const int myID = -1000;
 347     TNodeDistributor* myHyp = dynamic_cast<TNodeDistributor*>( aMesh.GetHypothesis( myID ));
 348     if ( !myHyp )
 349       myHyp = new TNodeDistributor( myID, 0, aMesh.GetGen() );
 350     return myHyp;
 351   }
 352   // -----------------------------------------------------------------------------
 353   bool Compute( vector< double > &                  positions,
 354                 gp_Pnt                              pIn,
 355                 gp_Pnt                              pOut,
 356                 SMESH_Mesh&                         aMesh,
 357                 const StdMeshers_LayerDistribution* hyp)
 358   {
 359     double len = pIn.Distance( pOut );
 360     if ( len <= DBL_MIN ) return error("Too close points of inner and outer shells");
 361
 362     if ( !hyp || !hyp->GetLayerDistribution() )
 363       return error( "Invalid LayerDistribution hypothesis");
 364     myUsedHyps.clear();
 365     myUsedHyps.push_back( hyp->GetLayerDistribution() );
 366
 367     TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut );
 368     SMESH_Hypothesis::Hypothesis_Status aStatus;
 369     if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
 370       return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
 371                     "with LayerDistribution hypothesis");
 372
 373     BRepAdaptor_Curve C3D(edge);
 374     double f = C3D.FirstParameter(), l = C3D.LastParameter();
 375     list< double > params;
 376     if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false ))
 377       return error("StdMeshers_Regular_1D failed to compute layers distribution");
 378
 379     positions.clear();
 380     positions.reserve( params.size() );
 381     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
 382       positions.push_back( *itU / len );
 383     return true;
 384   }
 385 protected:
 386   // -----------------------------------------------------------------------------
 387   TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen)
 388     : StdMeshers_Regular_1D( hypId, studyId, gen)
 389   {
 390   }
 391   // -----------------------------------------------------------------------------
 392   virtual const list <const SMESHDS_Hypothesis *> &
 393     GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
 394   {
 395     return myUsedHyps;
 396   }
 397   // -----------------------------------------------------------------------------
 398 };
 399
 400 //================================================================================
 401 /*!
 402  * \brief Compute positions of nodes between the internal and the external surfaces
 403   * \retval bool - is a success
 404  */
 405 //================================================================================
 406
 407 bool StdMeshers_RadialPrism_3D::computeLayerPositions(const gp_Pnt& pIn,
 408                                                       const gp_Pnt& pOut)
 409 {
 410   if ( myNbLayerHypo )
 411   {
 412     int nbSegments = myNbLayerHypo->GetNumberOfLayers();
 413     myLayerPositions.resize( nbSegments - 1 );
 414     for ( int z = 1; z < nbSegments; ++z )
 415       myLayerPositions[ z - 1 ] = double( z )/ double( nbSegments );
 416     return true;
 417   }
 418   if ( myDistributionHypo ) {
 419     SMESH_Mesh * mesh = myHelper->GetMesh();
 420     if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions, pIn, pOut,
 421                                                             *mesh, myDistributionHypo ))
 422     {
 423       error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
 424       return false;
 425     }
 426   }
 427   RETURN_BAD_RESULT("Bad hypothesis");
 428 }
 429
 430
 431 //=======================================================================
 432 //function : Evaluate
 433 //purpose  :
 434 //=======================================================================
 435
 436 bool StdMeshers_RadialPrism_3D::Evaluate(SMESH_Mesh& aMesh,
 437                                          const TopoDS_Shape& aShape,
 438                                          MapShapeNbElems& aResMap)
 439 {
 440   // get 2 shells
 441   TopoDS_Solid solid = TopoDS::Solid( aShape );
 442 #if OCC_VERSION_LARGE > 0x06050400
 443   TopoDS_Shell outerShell = BRepClass3d::OuterShell( solid );
 444 #else
 445   TopoDS_Shell outerShell = BRepTools::OuterShell( solid );
 446 #endif
 447   TopoDS_Shape innerShell;
 448   int nbShells = 0;
 449   for ( TopoDS_Iterator It (solid); It.More(); It.Next(), ++nbShells )
 450     if ( !outerShell.IsSame( It.Value() ))
 451       innerShell = It.Value();
 452   if ( nbShells != 2 ) {
 453     std::vector<int> aResVec(SMDSEntity_Last);
 454     for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 455     SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 456     aResMap.insert(std::make_pair(sm,aResVec));
 457     SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 458     smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 459     return false;
 460   }
 461
 462   // Associate sub-shapes of the shells
 463   ProjectionUtils::TShapeShapeMap shape2ShapeMap;
 464   if ( !ProjectionUtils::FindSubShapeAssociation( outerShell, &aMesh,
 465                                                   innerShell, &aMesh,
 466                                                   shape2ShapeMap) ) {
 467     std::vector<int> aResVec(SMDSEntity_Last);
 468     for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 469     SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 470     aResMap.insert(std::make_pair(sm,aResVec));
 471     SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 472     smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 473     return false;
 474   }
 475
 476   // get info for outer shell
 477   int nb0d_Out=0, nb2d_3_Out=0, nb2d_4_Out=0;
 478   //TopTools_SequenceOfShape FacesOut;
 479   for (TopExp_Explorer exp(outerShell, TopAbs_FACE); exp.More(); exp.Next()) {
 480     //FacesOut.Append(exp.Current());
 481     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
 482     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
 483     std::vector<int> aVec = (*anIt).second;
 484     nb0d_Out += aVec[SMDSEntity_Node];
 485     nb2d_3_Out += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
 486     nb2d_4_Out += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 487   }
 488   int nb1d_Out = 0;
 489   TopTools_MapOfShape tmpMap;
 490   for (TopExp_Explorer exp(outerShell, TopAbs_EDGE); exp.More(); exp.Next()) {
 491     if( tmpMap.Contains( exp.Current() ) )
 492       continue;
 493     tmpMap.Add( exp.Current() );
 494     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
 495     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
 496     std::vector<int> aVec = (*anIt).second;
 497     nb0d_Out += aVec[SMDSEntity_Node];
 498     nb1d_Out += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
 499   }
 500   tmpMap.Clear();
 501   for (TopExp_Explorer exp(outerShell, TopAbs_VERTEX); exp.More(); exp.Next()) {
 502     if( tmpMap.Contains( exp.Current() ) )
 503       continue;
 504     tmpMap.Add( exp.Current() );
 505     nb0d_Out++;
 506   }
 507
 508   // get info for inner shell
 509   int nb0d_In=0, nb2d_3_In=0, nb2d_4_In=0;
 510   //TopTools_SequenceOfShape FacesIn;
 511   for (TopExp_Explorer exp(innerShell, TopAbs_FACE); exp.More(); exp.Next()) {
 512     //FacesIn.Append(exp.Current());
 513     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
 514     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
 515     std::vector<int> aVec = (*anIt).second;
 516     nb0d_In += aVec[SMDSEntity_Node];
 517     nb2d_3_In += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
 518     nb2d_4_In += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 519   }
 520   int nb1d_In = 0;
 521   tmpMap.Clear();
 522   bool IsQuadratic = false;
 523   bool IsFirst = true;
 524   for (TopExp_Explorer exp(innerShell, TopAbs_EDGE); exp.More(); exp.Next()) {
 525     if( tmpMap.Contains( exp.Current() ) )
 526       continue;
 527     tmpMap.Add( exp.Current() );
 528     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
 529     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
 530     std::vector<int> aVec = (*anIt).second;
 531     nb0d_In += aVec[SMDSEntity_Node];
 532     nb1d_In += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
 533     if(IsFirst) {
 534       IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
 535       IsFirst = false;
 536     }
 537   }
 538   tmpMap.Clear();
 539   for (TopExp_Explorer exp(innerShell, TopAbs_VERTEX); exp.More(); exp.Next()) {
 540     if( tmpMap.Contains( exp.Current() ) )
 541       continue;
 542     tmpMap.Add( exp.Current() );
 543     nb0d_In++;
 544   }
 545
 546   bool IsOK = (nb0d_Out==nb0d_In) && (nb1d_Out==nb1d_In) &&
 547               (nb2d_3_Out==nb2d_3_In) && (nb2d_4_Out==nb2d_4_In);
 548   if(!IsOK) {
 549     std::vector<int> aResVec(SMDSEntity_Last);
 550     for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 551     SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 552     aResMap.insert(std::make_pair(sm,aResVec));
 553     SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 554     smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 555     return false;
 556   }
 557
 558   int nbLayers = 0;
 559   if( myNbLayerHypo ) {
 560     nbLayers = myNbLayerHypo->GetNumberOfLayers();
 561   }
 562   if ( myDistributionHypo ) {
 563     if ( !myDistributionHypo->GetLayerDistribution() ) {
 564       std::vector<int> aResVec(SMDSEntity_Last);
 565       for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 566       SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 567       aResMap.insert(std::make_pair(sm,aResVec));
 568       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 569       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 570       return false;
 571     }
 572     TopExp_Explorer exp(outerShell, TopAbs_VERTEX);
 573     TopoDS_Vertex Vout = TopoDS::Vertex(exp.Current());
 574     TopoDS_Vertex Vin = TopoDS::Vertex( shape2ShapeMap(Vout) );
 575     if ( myLayerPositions.empty() ) {
 576       gp_Pnt pIn = BRep_Tool::Pnt(Vin);
 577       gp_Pnt pOut = BRep_Tool::Pnt(Vout);
 578       computeLayerPositions( pIn, pOut );
 579     }
 580     nbLayers = myLayerPositions.size() + 1;
 581   }
 582
 583   std::vector<int> aResVec(SMDSEntity_Last);
 584   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 585   if(IsQuadratic) {
 586     aResVec[SMDSEntity_Quad_Penta] = nb2d_3_Out * nbLayers;
 587     aResVec[SMDSEntity_Quad_Hexa] = nb2d_4_Out * nbLayers;
 588     int nb1d = ( nb2d_3_Out*3 + nb2d_4_Out*4 ) / 2;
 589     aResVec[SMDSEntity_Node] = nb0d_Out * ( 2*nbLayers - 1 ) - nb1d * nbLayers;
 590   }
 591   else {
 592     aResVec[SMDSEntity_Node] = nb0d_Out * ( nbLayers - 1 );
 593     aResVec[SMDSEntity_Penta] = nb2d_3_Out * nbLayers;
 594     aResVec[SMDSEntity_Hexa] = nb2d_4_Out * nbLayers;
 595   }
 596   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 597   aResMap.insert(std::make_pair(sm,aResVec));
 598
 599   return true;
 600 }
 601
 602 //================================================================================
 603 /*!
 604  * \brief Return true if the algorithm can mesh this shape
 605  *  \param [in] aShape - shape to check
 606  *  \param [in] toCheckAll - if true, this check returns OK if all shapes are OK,
 607  *              else, returns OK if all at least one shape is OK
 608  */
 609 //================================================================================
 610
 611 bool StdMeshers_RadialPrism_3D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
 612 {
 613   bool isCurShellApp;
 614   int nbFoundSolids = 0;
 615   for (TopExp_Explorer exp( aShape, TopAbs_SOLID ); exp.More(); exp.Next(), ++nbFoundSolids )
 616   {
 617     TopoDS_Shape shell[2];
 618     int nbShells = 0;
 619     for ( TopoDS_Iterator It (exp.Current()); It.More(); It.Next() )
 620     {
 621       nbShells++;
 622       if ( nbShells > 2 ) {
 623         if ( toCheckAll ) return false;
 624         break;
 625       }
 626       shell[ nbShells-1 ] = It.Value();
 627     }
 628     if ( nbShells != 2 ) continue;
 629
 630     int nbFaces1 = SMESH_MesherHelper:: Count( shell[0], TopAbs_FACE, 0 );
 631     int nbFaces2 = SMESH_MesherHelper:: Count( shell[1], TopAbs_FACE, 0 );
 632     if ( nbFaces1 != nbFaces2 ){
 633       if( toCheckAll ) return false;
 634       continue;
 635     }
 636     int nbEdges1 = SMESH_MesherHelper:: Count( shell[0], TopAbs_EDGE, 0 );
 637     int nbEdges2 = SMESH_MesherHelper:: Count( shell[1], TopAbs_EDGE, 0 );
 638     if ( nbEdges1 != nbEdges2 ){
 639       if( toCheckAll ) return false;
 640       continue;
 641     }
 642     int nbVertices1 = SMESH_MesherHelper:: Count( shell[0], TopAbs_VERTEX, 0 );
 643     int nbVertices2 = SMESH_MesherHelper:: Count( shell[1], TopAbs_VERTEX, 0 );
 644     if ( nbVertices1 != nbVertices2 ){
 645       if( toCheckAll ) return false;
 646       continue;
 647     }
 648     if ( !toCheckAll ) return true;
 649   }
 650   return ( toCheckAll && nbFoundSolids != 0);
 651 };