src/StdMeshers/StdMeshers_RadialPrism_3D.cxx

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