src/StdMeshers/StdMeshers_RadialPrism_3D.cxx

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