Back to index

salome-smesh  6.5.0
StdMeshers_Hexa_3D.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 //  SMESH SMESH : implementaion of SMESH idl descriptions
00024 //  File   : StdMeshers_Hexa_3D.cxx
00025 //           Moved here from SMESH_Hexa_3D.cxx
00026 //  Author : Paul RASCLE, EDF
00027 //  Module : SMESH
00028 //
00029 #include "StdMeshers_Hexa_3D.hxx"
00030 
00031 #include "StdMeshers_CompositeHexa_3D.hxx"
00032 #include "StdMeshers_FaceSide.hxx"
00033 #include "StdMeshers_HexaFromSkin_3D.hxx"
00034 #include "StdMeshers_Penta_3D.hxx"
00035 #include "StdMeshers_Prism_3D.hxx"
00036 #include "StdMeshers_Quadrangle_2D.hxx"
00037 #include "StdMeshers_ViscousLayers.hxx"
00038 
00039 #include "SMESH_Comment.hxx"
00040 #include "SMESH_Gen.hxx"
00041 #include "SMESH_Mesh.hxx"
00042 #include "SMESH_MesherHelper.hxx"
00043 #include "SMESH_subMesh.hxx"
00044 
00045 #include "SMDS_MeshNode.hxx"
00046 
00047 #include <TopExp.hxx>
00048 #include <TopExp_Explorer.hxx>
00049 #include <TopTools_SequenceOfShape.hxx>
00050 #include <TopTools_MapOfShape.hxx>
00051 #include <TopoDS.hxx>
00052 
00053 #include "utilities.h"
00054 #include "Utils_ExceptHandlers.hxx"
00055 
00056 typedef SMESH_Comment TComm;
00057 
00058 using namespace std;
00059 
00060 static SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &,
00061                                                     const TopoDS_Shape &,
00062                                                     SMESH_ProxyMesh* proxyMesh=0);
00063 
00064 static bool EvaluatePentahedralMesh(SMESH_Mesh &, const TopoDS_Shape &,
00065                                     MapShapeNbElems &);
00066 
00067 //=============================================================================
00071 //=============================================================================
00072 
00073 StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, int studyId, SMESH_Gen * gen)
00074   :SMESH_3D_Algo(hypId, studyId, gen)
00075 {
00076   MESSAGE("StdMeshers_Hexa_3D::StdMeshers_Hexa_3D");
00077   _name = "Hexa_3D";
00078   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);       // 1 bit /shape type
00079   _requireShape = false;
00080   _compatibleHypothesis.push_back("ViscousLayers");
00081 }
00082 
00083 //=============================================================================
00087 //=============================================================================
00088 
00089 StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D()
00090 {
00091   MESSAGE("StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D");
00092 }
00093 
00094 //=============================================================================
00098 //=============================================================================
00099 
00100 bool StdMeshers_Hexa_3D::CheckHypothesis
00101                          (SMESH_Mesh&                          aMesh,
00102                           const TopoDS_Shape&                  aShape,
00103                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
00104 {
00105   // check nb of faces in the shape
00106 /*  PAL16229
00107   aStatus = SMESH_Hypothesis::HYP_BAD_GEOMETRY;
00108   int nbFaces = 0;
00109   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next())
00110     if ( ++nbFaces > 6 )
00111       break;
00112   if ( nbFaces != 6 )
00113     return false;
00114 */
00115 
00116   _viscousLayersHyp = NULL;
00117 
00118   const list<const SMESHDS_Hypothesis*>& hyps =
00119     GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
00120   list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
00121   if ( h == hyps.end())
00122   {
00123     aStatus = SMESH_Hypothesis::HYP_OK;
00124     return true;
00125   }
00126 
00127   aStatus = HYP_OK;
00128   for ( ; h != hyps.end(); ++h )
00129   {
00130     string hypName = (*h)->GetName();
00131     if ( find( _compatibleHypothesis.begin(),_compatibleHypothesis.end(),hypName )
00132          != _compatibleHypothesis.end() )
00133     {
00134       _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h );
00135     }
00136     else
00137     {
00138       aStatus = HYP_INCOMPATIBLE;
00139     }
00140   }
00141 
00142   if ( !_viscousLayersHyp )
00143     aStatus = HYP_INCOMPATIBLE;
00144 
00145   return aStatus == HYP_OK;
00146 }
00147 
00148 namespace
00149 {
00150   //=============================================================================
00151 
00152   typedef boost::shared_ptr< FaceQuadStruct > FaceQuadStructPtr;
00153 
00154   // symbolic names of box sides
00155   enum EBoxSides{ B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_NB_SIDES };
00156 
00157   // symbolic names of sides of quadrangle
00158   enum EQuadSides{ Q_BOTTOM=0, Q_RIGHT, Q_TOP, Q_LEFT, Q_NB_SIDES };
00159 
00160   //=============================================================================
00164   struct _FaceGrid
00165   {
00166     // face sides
00167     FaceQuadStructPtr _quad;
00168 
00169     // map of (node parameter on EDGE) to (column (vector) of nodes)
00170     TParam2ColumnMap _u2nodesMap;
00171 
00172     // node column's taken form _u2nodesMap taking into account sub-shape orientation
00173     vector<TNodeColumn> _columns;
00174 
00175     // geometry of a cube side
00176     TopoDS_Face _sideF;
00177 
00178     const SMDS_MeshNode* GetNode(int iCol, int iRow) const
00179     {
00180       return _columns[iCol][iRow];
00181     }
00182     gp_XYZ GetXYZ(int iCol, int iRow) const
00183     {
00184       return SMESH_TNodeXYZ( GetNode( iCol, iRow ));
00185     }
00186   };
00187 
00188   //================================================================================
00192   struct _Indexer
00193   {
00194     int _xSize, _ySize;
00195     _Indexer( int xSize, int ySize ): _xSize(xSize), _ySize(ySize) {}
00196     int size() const { return _xSize * _ySize; }
00197     int operator()(const int x, const int y) const { return y * _xSize + x; }
00198   };
00199 
00200   //================================================================================
00204   template< class TMapIterator >
00205   void append( TParam2ColumnMap& toMap, TMapIterator from, TMapIterator to )
00206   {
00207     const SMDS_MeshNode* lastNode = toMap.rbegin()->second[0];
00208     const SMDS_MeshNode* firstNode = from->second[0];
00209     if ( lastNode == firstNode )
00210       from++;
00211     double u = toMap.rbegin()->first;
00212     for (; from != to; ++from )
00213     {
00214       u += 1;
00215       TParam2ColumnMap::iterator u2nn = toMap.insert( toMap.end(), make_pair ( u, TNodeColumn()));
00216       u2nn->second.swap( from->second );
00217     }
00218   }
00219 
00220   //================================================================================
00225   FaceQuadStructPtr getQuadWithBottom( StdMeshers_FaceSide* side,
00226                                        FaceQuadStructPtr    quad[ 6 ])
00227   {
00228     FaceQuadStructPtr foundQuad;
00229     for ( int i = 1; i < 6; ++i )
00230     {
00231       if ( !quad[i] ) continue;
00232       for ( unsigned iS = 0; iS < quad[i]->side.size(); ++iS )
00233       {
00234         const StdMeshers_FaceSide* side2 = quad[i]->side[iS];
00235         if (( side->FirstVertex().IsSame( side2->FirstVertex() ) ||
00236               side->FirstVertex().IsSame( side2->LastVertex() ))
00237             &&
00238             ( side->LastVertex().IsSame( side2->FirstVertex() ) ||
00239               side->LastVertex().IsSame( side2->LastVertex() ))
00240             )
00241         {
00242           if ( iS != Q_BOTTOM )
00243           {
00244             vector< StdMeshers_FaceSide*> newSides;
00245             for ( unsigned j = iS; j < quad[i]->side.size(); ++j )
00246               newSides.push_back( quad[i]->side[j] );
00247             for ( unsigned j = 0; j < iS; ++j )
00248               newSides.push_back( quad[i]->side[j] );
00249             quad[i]->side.swap( newSides );
00250           }
00251           foundQuad.swap(quad[i]);
00252           return foundQuad;
00253         }
00254       }
00255     }
00256     return foundQuad;
00257   }
00258   //================================================================================
00262   //================================================================================
00263 
00264   bool beginsAtSide( const _FaceGrid&     sideGrid1,
00265                      const _FaceGrid&     sideGrid2,
00266                      SMESH_ProxyMesh::Ptr proxymesh )
00267   {
00268     const TNodeColumn& col0  = sideGrid2._u2nodesMap.begin()->second;
00269     const TNodeColumn& col1  = sideGrid2._u2nodesMap.rbegin()->second;
00270     const SMDS_MeshNode* n00 = col0.front();
00271     const SMDS_MeshNode* n01 = col0.back();
00272     const SMDS_MeshNode* n10 = col1.front();
00273     const SMDS_MeshNode* n11 = col1.back();
00274     const SMDS_MeshNode* n = (sideGrid1._u2nodesMap.begin()->second)[0];
00275     if ( proxymesh )
00276     {
00277       n00 = proxymesh->GetProxyNode( n00 );
00278       n10 = proxymesh->GetProxyNode( n10 );
00279       n01 = proxymesh->GetProxyNode( n01 );
00280       n11 = proxymesh->GetProxyNode( n11 );
00281       n   = proxymesh->GetProxyNode( n );
00282     }
00283     return ( n == n00 || n == n01 || n == n10 || n == n11 );
00284   }
00285 }
00286 
00287 //=============================================================================
00295 //=============================================================================
00296 
00297 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
00298                                  const TopoDS_Shape & aShape)// throw(SALOME_Exception)
00299 {
00300   // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside
00301   //Unexpect aCatch(SalomeException);
00302   MESSAGE("StdMeshers_Hexa_3D::Compute");
00303   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
00304 
00305   // Shape verification
00306   // ----------------------
00307 
00308   // shape must be a solid (or a shell) with 6 faces
00309   TopExp_Explorer exp(aShape,TopAbs_SHELL);
00310   if ( !exp.More() )
00311     return error(COMPERR_BAD_SHAPE, "No SHELL in the geometry");
00312   if ( exp.Next(), exp.More() )
00313     return error(COMPERR_BAD_SHAPE, "More than one SHELL in the geometry");
00314 
00315   TopTools_IndexedMapOfShape FF;
00316   TopExp::MapShapes( aShape, TopAbs_FACE, FF);
00317   if ( FF.Extent() != 6)
00318   {
00319     static StdMeshers_CompositeHexa_3D compositeHexa(_gen->GetANewId(), 0, _gen);
00320     if ( !compositeHexa.Compute( aMesh, aShape ))
00321       return error( compositeHexa.GetComputeError() );
00322     return true;
00323   }
00324 
00325   // Find sides of a cube
00326   // ---------------------
00327   
00328   FaceQuadStructPtr quad[ 6 ];
00329   StdMeshers_Quadrangle_2D quadAlgo( _gen->GetANewId(), GetStudyId(), _gen);
00330   for ( int i = 0; i < 6; ++i )
00331   {
00332     if ( !( quad[i] = FaceQuadStructPtr( quadAlgo.CheckNbEdges( aMesh, FF( i+1 )))))
00333       return error( quadAlgo.GetComputeError() );
00334     if ( quad[i]->side.size() != 4 )
00335       return error( COMPERR_BAD_SHAPE, "Not a quadrangular box side" );
00336   }
00337 
00338   _FaceGrid aCubeSide[ 6 ];
00339 
00340   swap( aCubeSide[B_BOTTOM]._quad, quad[0] );
00341   swap( aCubeSide[B_BOTTOM]._quad->side[ Q_RIGHT],// direct the normal of bottom quad inside cube
00342         aCubeSide[B_BOTTOM]._quad->side[ Q_LEFT ] );
00343 
00344   aCubeSide[B_FRONT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_BOTTOM], quad );
00345   aCubeSide[B_RIGHT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_RIGHT ], quad );
00346   aCubeSide[B_BACK ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_TOP   ], quad );
00347   aCubeSide[B_LEFT ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_LEFT  ], quad );
00348   if ( aCubeSide[B_FRONT ]._quad )
00349     aCubeSide[B_TOP]._quad = getQuadWithBottom( aCubeSide[B_FRONT ]._quad->side[Q_TOP ], quad );
00350 
00351   for ( int i = 1; i < 6; ++i )
00352     if ( !aCubeSide[i]._quad )
00353       return error( COMPERR_BAD_SHAPE );
00354 
00355   // Make viscous layers
00356   // --------------------
00357 
00358   SMESH_ProxyMesh::Ptr proxymesh;
00359   if ( _viscousLayersHyp )
00360   {
00361     proxymesh = _viscousLayersHyp->Compute( aMesh, aShape, /*makeN2NMap=*/ true );
00362     if ( !proxymesh )
00363       return false;
00364   }
00365 
00366   // Check if there are triangles on cube sides
00367   // -------------------------------------------
00368 
00369   if ( aMesh.NbTriangles() > 0 )
00370   {
00371     for ( int i = 0; i < 6; ++i )
00372     {
00373       const TopoDS_Face& sideF = aCubeSide[i]._quad->face;
00374       if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( sideF ))
00375       {
00376         bool isAllQuad = true;
00377         SMDS_ElemIteratorPtr fIt = smDS->GetElements();
00378         while ( fIt->more() && isAllQuad )
00379         {
00380           const SMDS_MeshElement* f = fIt->next();
00381           isAllQuad = ( f->NbCornerNodes() == 4 );
00382         }
00383         if ( !isAllQuad )
00384         {
00385           SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
00386           return error( err );
00387         }
00388       }
00389     }
00390   }
00391 
00392   // Check presence of regular grid mesh on FACEs of the cube
00393   // ------------------------------------------------------------
00394 
00395   // tool creating quadratic elements if needed
00396   SMESH_MesherHelper helper (aMesh);
00397   _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
00398 
00399   for ( int i = 0; i < 6; ++i )
00400   {
00401     const TopoDS_Face& F = aCubeSide[i]._quad->face;
00402     StdMeshers_FaceSide* baseQuadSide = aCubeSide[i]._quad->side[ Q_BOTTOM ];
00403     list<TopoDS_Edge> baseEdges( baseQuadSide->Edges().begin(), baseQuadSide->Edges().end() );
00404 
00405     // assure correctness of node positions on baseE:
00406     // helper.GetNodeU() will fix positions if they are wrong
00407     for ( int iE = 0; iE < baseQuadSide->NbEdges(); ++iE )
00408     {
00409       const TopoDS_Edge& baseE = baseQuadSide->Edge( iE );
00410       if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( baseE ))
00411       {
00412         bool ok;
00413         helper.SetSubShape( baseE );
00414         SMDS_ElemIteratorPtr eIt = smDS->GetElements();
00415         while ( eIt->more() )
00416         {
00417           const SMDS_MeshElement* e = eIt->next();
00418           // expect problems on a composite side
00419           try { helper.GetNodeU( baseE, e->GetNode(0), e->GetNode(1), &ok); }
00420           catch (...) {}
00421           try { helper.GetNodeU( baseE, e->GetNode(1), e->GetNode(0), &ok); }
00422           catch (...) {}
00423         }
00424       }
00425     }
00426 
00427     // load grid
00428     bool ok =
00429       helper.LoadNodeColumns( aCubeSide[i]._u2nodesMap, F, baseEdges, meshDS, proxymesh.get());
00430     if ( ok )
00431     {
00432       // check if the loaded grid corresponds to nb of quadrangles on the FACE
00433       const SMESHDS_SubMesh* faceSubMesh =
00434         proxymesh ? proxymesh->GetSubMesh( F ) : meshDS->MeshElements( F );
00435       const int nbQuads = faceSubMesh->NbElements();
00436       const int nbHor = aCubeSide[i]._u2nodesMap.size() - 1;
00437       const int nbVer = aCubeSide[i]._u2nodesMap.begin()->second.size() - 1;
00438       ok = ( nbQuads == nbHor * nbVer );
00439     }
00440     if ( !ok )
00441     {
00442       SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
00443       return error( err );
00444     }
00445   }
00446 
00447   // Orient loaded grids of cube sides along axis of the unitary cube coord system
00448   bool isReverse[6];
00449   isReverse[B_BOTTOM] = beginsAtSide( aCubeSide[B_BOTTOM], aCubeSide[B_RIGHT ], proxymesh );
00450   isReverse[B_TOP   ] = beginsAtSide( aCubeSide[B_TOP   ], aCubeSide[B_RIGHT ], proxymesh );
00451   isReverse[B_FRONT ] = beginsAtSide( aCubeSide[B_FRONT ], aCubeSide[B_RIGHT ], proxymesh );
00452   isReverse[B_BACK  ] = beginsAtSide( aCubeSide[B_BACK  ], aCubeSide[B_RIGHT ], proxymesh );
00453   isReverse[B_LEFT  ] = beginsAtSide( aCubeSide[B_LEFT  ], aCubeSide[B_BACK  ], proxymesh );
00454   isReverse[B_RIGHT ] = beginsAtSide( aCubeSide[B_RIGHT ], aCubeSide[B_BACK  ], proxymesh );
00455   for ( int i = 0; i < 6; ++i )
00456   {
00457     aCubeSide[i]._columns.resize( aCubeSide[i]._u2nodesMap.size() );
00458 
00459     int iFwd = 0, iRev = aCubeSide[i]._columns.size()-1;
00460     int* pi = isReverse[i] ? &iRev : &iFwd;
00461     TParam2ColumnMap::iterator u2nn = aCubeSide[i]._u2nodesMap.begin();
00462     for ( ; iFwd < aCubeSide[i]._columns.size(); --iRev, ++iFwd, ++u2nn )
00463       aCubeSide[i]._columns[ *pi ].swap( u2nn->second );
00464 
00465     aCubeSide[i]._u2nodesMap.clear();
00466   }
00467   
00468   if ( proxymesh )
00469     for ( int i = 0; i < 6; ++i )
00470       for ( unsigned j = 0; j < aCubeSide[i]._columns.size(); ++j)
00471         for ( unsigned k = 0; k < aCubeSide[i]._columns[j].size(); ++k)
00472         {
00473           const SMDS_MeshNode* & n = aCubeSide[i]._columns[j][k];
00474           n = proxymesh->GetProxyNode( n );
00475         }
00476 
00477   // 4) Create internal nodes of the cube
00478   // -------------------------------------
00479 
00480   helper.SetSubShape( aShape );
00481   helper.SetElementsOnShape(true);
00482 
00483   // shortcuts to sides
00484   _FaceGrid* fBottom = & aCubeSide[ B_BOTTOM ];
00485   _FaceGrid* fRight  = & aCubeSide[ B_RIGHT  ];
00486   _FaceGrid* fTop    = & aCubeSide[ B_TOP    ];
00487   _FaceGrid* fLeft   = & aCubeSide[ B_LEFT   ];
00488   _FaceGrid* fFront  = & aCubeSide[ B_FRONT  ];
00489   _FaceGrid* fBack   = & aCubeSide[ B_BACK   ];
00490 
00491   // cube size measured in nb of nodes
00492   int x, xSize = fBottom->_columns.size() , X = xSize - 1;
00493   int y, ySize = fLeft->_columns.size()   , Y = ySize - 1;
00494   int z, zSize = fLeft->_columns[0].size(), Z = zSize - 1;
00495 
00496   // columns of internal nodes "rising" from nodes of fBottom
00497   _Indexer colIndex( xSize, ySize );
00498   vector< vector< const SMDS_MeshNode* > > columns( colIndex.size() );
00499 
00500   // fill node columns by front and back box sides
00501   for ( x = 0; x < xSize; ++x ) {
00502     vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
00503     vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
00504     column0.resize( zSize );
00505     column1.resize( zSize );
00506     for ( z = 0; z < zSize; ++z ) {
00507       column0[ z ] = fFront->GetNode( x, z );
00508       column1[ z ] = fBack ->GetNode( x, z );
00509     }
00510   }
00511   // fill node columns by left and right box sides
00512   for ( y = 1; y < ySize-1; ++y ) {
00513     vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
00514     vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
00515     column0.resize( zSize );
00516     column1.resize( zSize );
00517     for ( z = 0; z < zSize; ++z ) {
00518       column0[ z ] = fLeft ->GetNode( y, z );
00519       column1[ z ] = fRight->GetNode( y, z );
00520     }
00521   }
00522   // get nodes from top and bottom box sides
00523   for ( x = 1; x < xSize-1; ++x ) {
00524     for ( y = 1; y < ySize-1; ++y ) {
00525       vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
00526       column.resize( zSize );
00527       column.front() = fBottom->GetNode( x, y );
00528       column.back()  = fTop   ->GetNode( x, y );
00529     }
00530   }
00531 
00532   // projection points of the internal node on cube sub-shapes by which
00533   // coordinates of the internal node are computed
00534   vector<gp_XYZ> pointsOnShapes( SMESH_Block::ID_Shell );
00535 
00536   // projections on vertices are constant
00537   pointsOnShapes[ SMESH_Block::ID_V000 ] = fBottom->GetXYZ( 0, 0 );
00538   pointsOnShapes[ SMESH_Block::ID_V100 ] = fBottom->GetXYZ( X, 0 );
00539   pointsOnShapes[ SMESH_Block::ID_V010 ] = fBottom->GetXYZ( 0, Y );
00540   pointsOnShapes[ SMESH_Block::ID_V110 ] = fBottom->GetXYZ( X, Y );
00541   pointsOnShapes[ SMESH_Block::ID_V001 ] = fTop->GetXYZ( 0, 0 );
00542   pointsOnShapes[ SMESH_Block::ID_V101 ] = fTop->GetXYZ( X, 0 );
00543   pointsOnShapes[ SMESH_Block::ID_V011 ] = fTop->GetXYZ( 0, Y );
00544   pointsOnShapes[ SMESH_Block::ID_V111 ] = fTop->GetXYZ( X, Y );
00545 
00546   for ( x = 1; x < xSize-1; ++x )
00547   {
00548     gp_XYZ params; // normalized parameters of internal node within a unit box
00549     params.SetCoord( 1, x / double(X) );
00550     for ( y = 1; y < ySize-1; ++y )
00551     {
00552       params.SetCoord( 2, y / double(Y) );
00553       // a column to fill in during z loop
00554       vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
00555       // projection points on horizontal edges
00556       pointsOnShapes[ SMESH_Block::ID_Ex00 ] = fBottom->GetXYZ( x, 0 );
00557       pointsOnShapes[ SMESH_Block::ID_Ex10 ] = fBottom->GetXYZ( x, Y );
00558       pointsOnShapes[ SMESH_Block::ID_E0y0 ] = fBottom->GetXYZ( 0, y );
00559       pointsOnShapes[ SMESH_Block::ID_E1y0 ] = fBottom->GetXYZ( X, y );
00560       pointsOnShapes[ SMESH_Block::ID_Ex01 ] = fTop->GetXYZ( x, 0 );
00561       pointsOnShapes[ SMESH_Block::ID_Ex11 ] = fTop->GetXYZ( x, Y );
00562       pointsOnShapes[ SMESH_Block::ID_E0y1 ] = fTop->GetXYZ( 0, y );
00563       pointsOnShapes[ SMESH_Block::ID_E1y1 ] = fTop->GetXYZ( X, y );
00564       // projection points on horizontal faces
00565       pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = fBottom->GetXYZ( x, y );
00566       pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = fTop   ->GetXYZ( x, y );
00567       for ( z = 1; z < zSize-1; ++z ) // z loop
00568       {
00569         params.SetCoord( 3, z / double(Z) );
00570         // projection points on vertical edges
00571         pointsOnShapes[ SMESH_Block::ID_E00z ] = fFront->GetXYZ( 0, z );    
00572         pointsOnShapes[ SMESH_Block::ID_E10z ] = fFront->GetXYZ( X, z );    
00573         pointsOnShapes[ SMESH_Block::ID_E01z ] = fBack->GetXYZ( 0, z );    
00574         pointsOnShapes[ SMESH_Block::ID_E11z ] = fBack->GetXYZ( X, z );
00575         // projection points on vertical faces
00576         pointsOnShapes[ SMESH_Block::ID_Fx0z ] = fFront->GetXYZ( x, z );    
00577         pointsOnShapes[ SMESH_Block::ID_Fx1z ] = fBack ->GetXYZ( x, z );    
00578         pointsOnShapes[ SMESH_Block::ID_F0yz ] = fLeft ->GetXYZ( y, z );    
00579         pointsOnShapes[ SMESH_Block::ID_F1yz ] = fRight->GetXYZ( y, z );
00580 
00581         // compute internal node coordinates
00582         gp_XYZ coords;
00583         SMESH_Block::ShellPoint( params, pointsOnShapes, coords );
00584         column[ z ] = helper.AddNode( coords.X(), coords.Y(), coords.Z() );
00585 
00586       }
00587     }
00588   }
00589 
00590   // side data no more needed, free memory
00591   for ( int i = 0; i < 6; ++i )
00592     aCubeSide[i]._columns.clear();
00593 
00594   // 5) Create hexahedrons
00595   // ---------------------
00596 
00597   for ( x = 0; x < xSize-1; ++x ) {
00598     for ( y = 0; y < ySize-1; ++y ) {
00599       vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
00600       vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
00601       vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
00602       vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
00603       for ( z = 0; z < zSize-1; ++z )
00604       {
00605         // bottom face normal of a hexa mush point outside the volume
00606         helper.AddVolume(col00[z],   col01[z],   col11[z],   col10[z],
00607                          col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
00608       }
00609     }
00610   }
00611   return true;
00612 }
00613 
00614 //=============================================================================
00618 //=============================================================================
00619 
00620 bool StdMeshers_Hexa_3D::Evaluate(SMESH_Mesh & aMesh,
00621                                   const TopoDS_Shape & aShape,
00622                                   MapShapeNbElems& aResMap)
00623 {
00624   vector < SMESH_subMesh * >meshFaces;
00625   TopTools_SequenceOfShape aFaces;
00626   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
00627     aFaces.Append(exp.Current());
00628     SMESH_subMesh *aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
00629     ASSERT(aSubMesh);
00630     meshFaces.push_back(aSubMesh);
00631   }
00632   if (meshFaces.size() != 6) {
00633     //return error(COMPERR_BAD_SHAPE, TComm(meshFaces.size())<<" instead of 6 faces in a block");
00634     static StdMeshers_CompositeHexa_3D compositeHexa(-10, 0, aMesh.GetGen());
00635     return compositeHexa.Evaluate(aMesh, aShape, aResMap);
00636   }
00637   
00638   int i = 0;
00639   for(; i<6; i++) {
00640     //TopoDS_Shape aFace = meshFaces[i]->GetSubShape();
00641     TopoDS_Shape aFace = aFaces.Value(i+1);
00642     SMESH_Algo *algo = _gen->GetAlgo(aMesh, aFace);
00643     if( !algo ) {
00644       std::vector<int> aResVec(SMDSEntity_Last);
00645       for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
00646       SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
00647       aResMap.insert(std::make_pair(sm,aResVec));
00648       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
00649       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
00650       return false;
00651     }
00652     string algoName = algo->GetName();
00653     bool isAllQuad = false;
00654     if (algoName == "Quadrangle_2D") {
00655       MapShapeNbElemsItr anIt = aResMap.find(meshFaces[i]);
00656       if( anIt == aResMap.end() ) continue;
00657       std::vector<int> aVec = (*anIt).second;
00658       int nbtri = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
00659       if( nbtri == 0 )
00660         isAllQuad = true;
00661     }
00662     if ( ! isAllQuad ) {
00663       return EvaluatePentahedralMesh(aMesh, aShape, aResMap);
00664     }
00665   }
00666   
00667   // find number of 1d elems for 1 face
00668   int nb1d = 0;
00669   TopTools_MapOfShape Edges1;
00670   bool IsQuadratic = false;
00671   bool IsFirst = true;
00672   for (TopExp_Explorer exp(aFaces.Value(1), TopAbs_EDGE); exp.More(); exp.Next()) {
00673     Edges1.Add(exp.Current());
00674     SMESH_subMesh *sm = aMesh.GetSubMesh(exp.Current());
00675     if( sm ) {
00676       MapShapeNbElemsItr anIt = aResMap.find(sm);
00677       if( anIt == aResMap.end() ) continue;
00678       std::vector<int> aVec = (*anIt).second;
00679       nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
00680       if(IsFirst) {
00681         IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
00682         IsFirst = false;
00683       }
00684     }
00685   }
00686   // find face opposite to 1 face
00687   int OppNum = 0;
00688   for(i=2; i<=6; i++) {
00689     bool IsOpposite = true;
00690     for(TopExp_Explorer exp(aFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) {
00691       if( Edges1.Contains(exp.Current()) ) {
00692         IsOpposite = false;
00693         break;
00694       }
00695     }
00696     if(IsOpposite) {
00697       OppNum = i;
00698       break;
00699     }
00700   }
00701   // find number of 2d elems on side faces
00702   int nb2d = 0;
00703   for(i=2; i<=6; i++) {
00704     if( i == OppNum ) continue;
00705     MapShapeNbElemsItr anIt = aResMap.find( meshFaces[i-1] );
00706     if( anIt == aResMap.end() ) continue;
00707     std::vector<int> aVec = (*anIt).second;
00708     nb2d += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
00709   }
00710   
00711   MapShapeNbElemsItr anIt = aResMap.find( meshFaces[0] );
00712   std::vector<int> aVec = (*anIt).second;
00713   int nb2d_face0 = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
00714   int nb0d_face0 = aVec[SMDSEntity_Node];
00715 
00716   std::vector<int> aResVec(SMDSEntity_Last);
00717   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
00718   if(IsQuadratic) {
00719     aResVec[SMDSEntity_Quad_Hexa] = nb2d_face0 * ( nb2d/nb1d );
00720     int nb1d_face0_int = ( nb2d_face0*4 - nb1d ) / 2;
00721     aResVec[SMDSEntity_Node] = nb0d_face0 * ( 2*nb2d/nb1d - 1 ) - nb1d_face0_int * nb2d/nb1d;
00722   }
00723   else {
00724     aResVec[SMDSEntity_Node] = nb0d_face0 * ( nb2d/nb1d - 1 );
00725     aResVec[SMDSEntity_Hexa] = nb2d_face0 * ( nb2d/nb1d );
00726   }
00727   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
00728   aResMap.insert(std::make_pair(sm,aResVec));
00729 
00730   return true;
00731 }
00732 
00733 //================================================================================
00737 //================================================================================
00738 
00739 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
00740 {
00741   static StdMeshers_HexaFromSkin_3D * algo = 0;
00742   if ( !algo ) {
00743     SMESH_Gen* gen = aMesh.GetGen();
00744     algo = new StdMeshers_HexaFromSkin_3D( gen->GetANewId(), 0, gen );
00745   }
00746   algo->InitComputeError();
00747   algo->Compute( aMesh, aHelper );
00748   return error( algo->GetComputeError());
00749 }
00750 
00751 //=======================================================================
00752 //function : ComputePentahedralMesh
00753 //purpose  : 
00754 //=======================================================================
00755 
00756 SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &          aMesh,
00757                                              const TopoDS_Shape &  aShape,
00758                                              SMESH_ProxyMesh*      proxyMesh)
00759 {
00760   SMESH_ComputeErrorPtr err = SMESH_ComputeError::New();
00761   if ( proxyMesh )
00762   {
00763     err->myName = COMPERR_BAD_INPUT_MESH;
00764     err->myComment = "Can't build pentahedral mesh on viscous layers";
00765     return err;
00766   }
00767   bool bOK;
00768   StdMeshers_Penta_3D anAlgo;
00769   //
00770   bOK=anAlgo.Compute(aMesh, aShape);
00771   //
00772   err = anAlgo.GetComputeError();
00773   //
00774   if ( !bOK && anAlgo.ErrorStatus() == 5 )
00775   {
00776     static StdMeshers_Prism_3D * aPrism3D = 0;
00777     if ( !aPrism3D ) {
00778       SMESH_Gen* gen = aMesh.GetGen();
00779       aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), 0, gen );
00780     }
00781     SMESH_Hypothesis::Hypothesis_Status aStatus;
00782     if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
00783       aPrism3D->InitComputeError();
00784       bOK = aPrism3D->Compute( aMesh, aShape );
00785       err = aPrism3D->GetComputeError();
00786     }
00787   }
00788   return err;
00789 }
00790 
00791 
00792 //=======================================================================
00793 //function : EvaluatePentahedralMesh
00794 //purpose  : 
00795 //=======================================================================
00796 
00797 bool EvaluatePentahedralMesh(SMESH_Mesh & aMesh,
00798                              const TopoDS_Shape & aShape,
00799                              MapShapeNbElems& aResMap)
00800 {
00801   StdMeshers_Penta_3D anAlgo;
00802   bool bOK = anAlgo.Evaluate(aMesh, aShape, aResMap);
00803 
00804   //err = anAlgo.GetComputeError();
00805   //if ( !bOK && anAlgo.ErrorStatus() == 5 )
00806   if( !bOK ) {
00807     static StdMeshers_Prism_3D * aPrism3D = 0;
00808     if ( !aPrism3D ) {
00809       SMESH_Gen* gen = aMesh.GetGen();
00810       aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), 0, gen );
00811     }
00812     SMESH_Hypothesis::Hypothesis_Status aStatus;
00813     if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
00814       return aPrism3D->Evaluate(aMesh, aShape, aResMap);
00815     }
00816   }
00817 
00818   return bOK;
00819 }