1 // Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
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.
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.
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
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 // SMESH SMESH : implementaion of SMESH idl descriptions
24 // File : StdMeshers_Hexa_3D.cxx
25 // Moved here from SMESH_Hexa_3D.cxx
26 // Author : Paul RASCLE, EDF
29 #include "StdMeshers_Hexa_3D.hxx"
31 #include "StdMeshers_CompositeHexa_3D.hxx"
32 #include "StdMeshers_FaceSide.hxx"
33 #include "StdMeshers_HexaFromSkin_3D.hxx"
34 #include "StdMeshers_Penta_3D.hxx"
35 #include "StdMeshers_Prism_3D.hxx"
36 #include "StdMeshers_Quadrangle_2D.hxx"
37 #include "StdMeshers_ViscousLayers.hxx"
39 #include "SMESH_Comment.hxx"
40 #include "SMESH_Gen.hxx"
41 #include "SMESH_Mesh.hxx"
42 #include "SMESH_MesherHelper.hxx"
43 #include "SMESH_subMesh.hxx"
45 #include "SMDS_MeshNode.hxx"
48 #include <TopExp_Explorer.hxx>
49 #include <TopTools_SequenceOfShape.hxx>
50 #include <TopTools_MapOfShape.hxx>
53 #include "utilities.h"
54 #include "Utils_ExceptHandlers.hxx"
56 typedef SMESH_Comment TComm;
60 static SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &,
62 SMESH_ProxyMesh* proxyMesh=0);
64 static bool EvaluatePentahedralMesh(SMESH_Mesh &, const TopoDS_Shape &,
67 //=============================================================================
71 //=============================================================================
73 StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, int studyId, SMESH_Gen * gen)
74 :SMESH_3D_Algo(hypId, studyId, gen)
76 MESSAGE("StdMeshers_Hexa_3D::StdMeshers_Hexa_3D");
78 _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID); // 1 bit /shape type
79 _requireShape = false;
80 _compatibleHypothesis.push_back("ViscousLayers");
83 //=============================================================================
87 //=============================================================================
89 StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D()
91 MESSAGE("StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D");
94 //=============================================================================
96 * Retrieves defined hypotheses
98 //=============================================================================
100 bool StdMeshers_Hexa_3D::CheckHypothesis
102 const TopoDS_Shape& aShape,
103 SMESH_Hypothesis::Hypothesis_Status& aStatus)
105 // check nb of faces in the shape
107 aStatus = SMESH_Hypothesis::HYP_BAD_GEOMETRY;
109 for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next())
116 _viscousLayersHyp = NULL;
118 const list<const SMESHDS_Hypothesis*>& hyps =
119 GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
120 list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
121 if ( h == hyps.end())
123 aStatus = SMESH_Hypothesis::HYP_OK;
128 for ( ; h != hyps.end(); ++h )
130 string hypName = (*h)->GetName();
131 if ( find( _compatibleHypothesis.begin(),_compatibleHypothesis.end(),hypName )
132 != _compatibleHypothesis.end() )
134 _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h );
138 aStatus = HYP_INCOMPATIBLE;
142 if ( !_viscousLayersHyp )
143 aStatus = HYP_INCOMPATIBLE;
145 return aStatus == HYP_OK;
150 //=============================================================================
152 typedef boost::shared_ptr< FaceQuadStruct > FaceQuadStructPtr;
154 // symbolic names of box sides
155 enum EBoxSides{ B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_NB_SIDES };
157 // symbolic names of sides of quadrangle
158 enum EQuadSides{ Q_BOTTOM=0, Q_RIGHT, Q_TOP, Q_LEFT, Q_NB_SIDES };
160 //=============================================================================
162 * \brief Container of nodes of structured mesh on a qudrangular geom FACE
167 FaceQuadStructPtr _quad;
169 // map of (node parameter on EDGE) to (column (vector) of nodes)
170 TParam2ColumnMap _u2nodesMap;
172 // node column's taken form _u2nodesMap taking into account sub-shape orientation
173 vector<TNodeColumn> _columns;
175 // geometry of a cube side
178 const SMDS_MeshNode* GetNode(int iCol, int iRow) const
180 return _columns[iCol][iRow];
182 gp_XYZ GetXYZ(int iCol, int iRow) const
184 return SMESH_TNodeXYZ( GetNode( iCol, iRow ));
188 //================================================================================
190 * \brief Convertor of a pair of integers to a sole index
195 _Indexer( int xSize, int ySize ): _xSize(xSize), _ySize(ySize) {}
196 int size() const { return _xSize * _ySize; }
197 int operator()(const int x, const int y) const { return y * _xSize + x; }
200 //================================================================================
202 * \brief Appends a range of node columns from a map to another map
204 template< class TMapIterator >
205 void append( TParam2ColumnMap& toMap, TMapIterator from, TMapIterator to )
207 const SMDS_MeshNode* lastNode = toMap.rbegin()->second[0];
208 const SMDS_MeshNode* firstNode = from->second[0];
209 if ( lastNode == firstNode )
211 double u = toMap.rbegin()->first;
212 for (; from != to; ++from )
215 TParam2ColumnMap::iterator u2nn = toMap.insert( toMap.end(), make_pair ( u, TNodeColumn()));
216 u2nn->second.swap( from->second );
220 //================================================================================
222 * \brief Finds FaceQuadStruct having a side equal to a given one and rearranges
223 * the found FaceQuadStruct::side to have the given side at a Q_BOTTOM place
225 FaceQuadStructPtr getQuadWithBottom( StdMeshers_FaceSide* side,
226 FaceQuadStructPtr quad[ 6 ])
228 FaceQuadStructPtr foundQuad;
229 for ( int i = 1; i < 6; ++i )
231 if ( !quad[i] ) continue;
232 for ( unsigned iS = 0; iS < quad[i]->side.size(); ++iS )
234 const StdMeshers_FaceSide* side2 = quad[i]->side[iS];
235 if (( side->FirstVertex().IsSame( side2->FirstVertex() ) ||
236 side->FirstVertex().IsSame( side2->LastVertex() ))
238 ( side->LastVertex().IsSame( side2->FirstVertex() ) ||
239 side->LastVertex().IsSame( side2->LastVertex() ))
242 if ( iS != Q_BOTTOM )
244 vector< StdMeshers_FaceSide*> newSides;
245 for ( unsigned j = iS; j < quad[i]->side.size(); ++j )
246 newSides.push_back( quad[i]->side[j] );
247 for ( unsigned j = 0; j < iS; ++j )
248 newSides.push_back( quad[i]->side[j] );
249 quad[i]->side.swap( newSides );
251 foundQuad.swap(quad[i]);
258 //================================================================================
260 * \brief Returns true if the 1st base node of sideGrid1 belongs to sideGrid2
262 //================================================================================
264 bool beginsAtSide( const _FaceGrid& sideGrid1,
265 const _FaceGrid& sideGrid2,
266 SMESH_ProxyMesh::Ptr proxymesh )
268 const TNodeColumn& col0 = sideGrid2._u2nodesMap.begin()->second;
269 const TNodeColumn& col1 = sideGrid2._u2nodesMap.rbegin()->second;
270 const SMDS_MeshNode* n00 = col0.front();
271 const SMDS_MeshNode* n01 = col0.back();
272 const SMDS_MeshNode* n10 = col1.front();
273 const SMDS_MeshNode* n11 = col1.back();
274 const SMDS_MeshNode* n = (sideGrid1._u2nodesMap.begin()->second)[0];
277 n00 = proxymesh->GetProxyNode( n00 );
278 n10 = proxymesh->GetProxyNode( n10 );
279 n01 = proxymesh->GetProxyNode( n01 );
280 n11 = proxymesh->GetProxyNode( n11 );
281 n = proxymesh->GetProxyNode( n );
283 return ( n == n00 || n == n01 || n == n10 || n == n11 );
287 //=============================================================================
289 * Generates hexahedron mesh on hexaedron like form using algorithm from
290 * "Application de l'interpolation transfinie à la création de maillages
291 * C0 ou G1 continus sur des triangles, quadrangles, tetraedres, pentaedres
292 * et hexaedres déformés."
293 * Alain PERONNET - 8 janvier 1999
295 //=============================================================================
297 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
298 const TopoDS_Shape & aShape)
300 // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside
301 //Unexpect aCatch(SalomeException);
302 MESSAGE("StdMeshers_Hexa_3D::Compute");
303 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
305 // Shape verification
306 // ----------------------
308 // shape must be a solid (or a shell) with 6 faces
309 TopExp_Explorer exp(aShape,TopAbs_SHELL);
311 return error(COMPERR_BAD_SHAPE, "No SHELL in the geometry");
312 if ( exp.Next(), exp.More() )
313 return error(COMPERR_BAD_SHAPE, "More than one SHELL in the geometry");
315 TopTools_IndexedMapOfShape FF;
316 TopExp::MapShapes( aShape, TopAbs_FACE, FF);
317 if ( FF.Extent() != 6)
319 static StdMeshers_CompositeHexa_3D compositeHexa(_gen->GetANewId(), 0, _gen);
320 if ( !compositeHexa.Compute( aMesh, aShape ))
321 return error( compositeHexa.GetComputeError() );
325 // Find sides of a cube
326 // ---------------------
328 FaceQuadStructPtr quad[ 6 ];
329 StdMeshers_Quadrangle_2D quadAlgo( _gen->GetANewId(), GetStudyId(), _gen);
330 for ( int i = 0; i < 6; ++i )
332 if ( !( quad[i] = FaceQuadStructPtr( quadAlgo.CheckNbEdges( aMesh, FF( i+1 )))))
333 return error( quadAlgo.GetComputeError() );
334 if ( quad[i]->side.size() != 4 )
335 return error( COMPERR_BAD_SHAPE, "Not a quadrangular box side" );
338 _FaceGrid aCubeSide[ 6 ];
340 swap( aCubeSide[B_BOTTOM]._quad, quad[0] );
341 swap( aCubeSide[B_BOTTOM]._quad->side[ Q_RIGHT],// direct the normal of bottom quad inside cube
342 aCubeSide[B_BOTTOM]._quad->side[ Q_LEFT ] );
344 aCubeSide[B_FRONT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_BOTTOM], quad );
345 aCubeSide[B_RIGHT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_RIGHT ], quad );
346 aCubeSide[B_BACK ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_TOP ], quad );
347 aCubeSide[B_LEFT ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_LEFT ], quad );
348 if ( aCubeSide[B_FRONT ]._quad )
349 aCubeSide[B_TOP]._quad = getQuadWithBottom( aCubeSide[B_FRONT ]._quad->side[Q_TOP ], quad );
351 for ( int i = 1; i < 6; ++i )
352 if ( !aCubeSide[i]._quad )
353 return error( COMPERR_BAD_SHAPE );
355 // Make viscous layers
356 // --------------------
358 SMESH_ProxyMesh::Ptr proxymesh;
359 if ( _viscousLayersHyp )
361 proxymesh = _viscousLayersHyp->Compute( aMesh, aShape, /*makeN2NMap=*/ true );
366 // Check if there are triangles on cube sides
367 // -------------------------------------------
369 if ( aMesh.NbTriangles() > 0 )
371 for ( int i = 0; i < 6; ++i )
373 const TopoDS_Face& sideF = aCubeSide[i]._quad->face;
374 if ( !SMESH_MesherHelper::IsSameElemGeometry( meshDS->MeshElements( sideF ),
376 /*nullSubMeshRes=*/false ))
378 SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
384 // Check presence of regular grid mesh on FACEs of the cube
385 // ------------------------------------------------------------
387 // tool creating quadratic elements if needed
388 SMESH_MesherHelper helper (aMesh);
389 _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
391 for ( int i = 0; i < 6; ++i )
393 const TopoDS_Face& F = aCubeSide[i]._quad->face;
394 StdMeshers_FaceSide* baseQuadSide = aCubeSide[i]._quad->side[ Q_BOTTOM ];
395 list<TopoDS_Edge> baseEdges( baseQuadSide->Edges().begin(), baseQuadSide->Edges().end() );
397 // assure correctness of node positions on baseE:
398 // helper.GetNodeU() will fix positions if they are wrong
399 helper.ToFixNodeParameters( true );
400 for ( int iE = 0; iE < baseQuadSide->NbEdges(); ++iE )
402 const TopoDS_Edge& baseE = baseQuadSide->Edge( iE );
403 if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( baseE ))
406 helper.SetSubShape( baseE );
407 SMDS_ElemIteratorPtr eIt = smDS->GetElements();
408 while ( eIt->more() )
410 const SMDS_MeshElement* e = eIt->next();
411 // expect problems on a composite side
412 try { helper.GetNodeU( baseE, e->GetNode(0), e->GetNode(1), &ok); }
414 try { helper.GetNodeU( baseE, e->GetNode(1), e->GetNode(0), &ok); }
422 helper.LoadNodeColumns( aCubeSide[i]._u2nodesMap, F, baseEdges, meshDS, proxymesh.get());
425 // check if the loaded grid corresponds to nb of quadrangles on the FACE
426 const SMESHDS_SubMesh* faceSubMesh =
427 proxymesh ? proxymesh->GetSubMesh( F ) : meshDS->MeshElements( F );
428 const int nbQuads = faceSubMesh->NbElements();
429 const int nbHor = aCubeSide[i]._u2nodesMap.size() - 1;
430 const int nbVer = aCubeSide[i]._u2nodesMap.begin()->second.size() - 1;
431 ok = ( nbQuads == nbHor * nbVer );
435 SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
440 // Orient loaded grids of cube sides along axis of the unitary cube coord system
442 isReverse[B_BOTTOM] = beginsAtSide( aCubeSide[B_BOTTOM], aCubeSide[B_RIGHT ], proxymesh );
443 isReverse[B_TOP ] = beginsAtSide( aCubeSide[B_TOP ], aCubeSide[B_RIGHT ], proxymesh );
444 isReverse[B_FRONT ] = beginsAtSide( aCubeSide[B_FRONT ], aCubeSide[B_RIGHT ], proxymesh );
445 isReverse[B_BACK ] = beginsAtSide( aCubeSide[B_BACK ], aCubeSide[B_RIGHT ], proxymesh );
446 isReverse[B_LEFT ] = beginsAtSide( aCubeSide[B_LEFT ], aCubeSide[B_BACK ], proxymesh );
447 isReverse[B_RIGHT ] = beginsAtSide( aCubeSide[B_RIGHT ], aCubeSide[B_BACK ], proxymesh );
448 for ( int i = 0; i < 6; ++i )
450 aCubeSide[i]._columns.resize( aCubeSide[i]._u2nodesMap.size() );
452 int iFwd = 0, iRev = aCubeSide[i]._columns.size()-1;
453 int* pi = isReverse[i] ? &iRev : &iFwd;
454 TParam2ColumnMap::iterator u2nn = aCubeSide[i]._u2nodesMap.begin();
455 for ( ; iFwd < aCubeSide[i]._columns.size(); --iRev, ++iFwd, ++u2nn )
456 aCubeSide[i]._columns[ *pi ].swap( u2nn->second );
458 aCubeSide[i]._u2nodesMap.clear();
462 for ( int i = 0; i < 6; ++i )
463 for ( unsigned j = 0; j < aCubeSide[i]._columns.size(); ++j)
464 for ( unsigned k = 0; k < aCubeSide[i]._columns[j].size(); ++k)
466 const SMDS_MeshNode* & n = aCubeSide[i]._columns[j][k];
467 n = proxymesh->GetProxyNode( n );
470 // 4) Create internal nodes of the cube
471 // -------------------------------------
473 helper.SetSubShape( aShape );
474 helper.SetElementsOnShape(true);
476 // shortcuts to sides
477 _FaceGrid* fBottom = & aCubeSide[ B_BOTTOM ];
478 _FaceGrid* fRight = & aCubeSide[ B_RIGHT ];
479 _FaceGrid* fTop = & aCubeSide[ B_TOP ];
480 _FaceGrid* fLeft = & aCubeSide[ B_LEFT ];
481 _FaceGrid* fFront = & aCubeSide[ B_FRONT ];
482 _FaceGrid* fBack = & aCubeSide[ B_BACK ];
484 // cube size measured in nb of nodes
485 int x, xSize = fBottom->_columns.size() , X = xSize - 1;
486 int y, ySize = fLeft->_columns.size() , Y = ySize - 1;
487 int z, zSize = fLeft->_columns[0].size(), Z = zSize - 1;
489 // columns of internal nodes "rising" from nodes of fBottom
490 _Indexer colIndex( xSize, ySize );
491 vector< vector< const SMDS_MeshNode* > > columns( colIndex.size() );
493 // fill node columns by front and back box sides
494 for ( x = 0; x < xSize; ++x ) {
495 vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
496 vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
497 column0.resize( zSize );
498 column1.resize( zSize );
499 for ( z = 0; z < zSize; ++z ) {
500 column0[ z ] = fFront->GetNode( x, z );
501 column1[ z ] = fBack ->GetNode( x, z );
504 // fill node columns by left and right box sides
505 for ( y = 1; y < ySize-1; ++y ) {
506 vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
507 vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
508 column0.resize( zSize );
509 column1.resize( zSize );
510 for ( z = 0; z < zSize; ++z ) {
511 column0[ z ] = fLeft ->GetNode( y, z );
512 column1[ z ] = fRight->GetNode( y, z );
515 // get nodes from top and bottom box sides
516 for ( x = 1; x < xSize-1; ++x ) {
517 for ( y = 1; y < ySize-1; ++y ) {
518 vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
519 column.resize( zSize );
520 column.front() = fBottom->GetNode( x, y );
521 column.back() = fTop ->GetNode( x, y );
525 // projection points of the internal node on cube sub-shapes by which
526 // coordinates of the internal node are computed
527 vector<gp_XYZ> pointsOnShapes( SMESH_Block::ID_Shell );
529 // projections on vertices are constant
530 pointsOnShapes[ SMESH_Block::ID_V000 ] = fBottom->GetXYZ( 0, 0 );
531 pointsOnShapes[ SMESH_Block::ID_V100 ] = fBottom->GetXYZ( X, 0 );
532 pointsOnShapes[ SMESH_Block::ID_V010 ] = fBottom->GetXYZ( 0, Y );
533 pointsOnShapes[ SMESH_Block::ID_V110 ] = fBottom->GetXYZ( X, Y );
534 pointsOnShapes[ SMESH_Block::ID_V001 ] = fTop->GetXYZ( 0, 0 );
535 pointsOnShapes[ SMESH_Block::ID_V101 ] = fTop->GetXYZ( X, 0 );
536 pointsOnShapes[ SMESH_Block::ID_V011 ] = fTop->GetXYZ( 0, Y );
537 pointsOnShapes[ SMESH_Block::ID_V111 ] = fTop->GetXYZ( X, Y );
539 for ( x = 1; x < xSize-1; ++x )
541 gp_XYZ params; // normalized parameters of internal node within a unit box
542 params.SetCoord( 1, x / double(X) );
543 for ( y = 1; y < ySize-1; ++y )
545 params.SetCoord( 2, y / double(Y) );
546 // a column to fill in during z loop
547 vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
548 // projection points on horizontal edges
549 pointsOnShapes[ SMESH_Block::ID_Ex00 ] = fBottom->GetXYZ( x, 0 );
550 pointsOnShapes[ SMESH_Block::ID_Ex10 ] = fBottom->GetXYZ( x, Y );
551 pointsOnShapes[ SMESH_Block::ID_E0y0 ] = fBottom->GetXYZ( 0, y );
552 pointsOnShapes[ SMESH_Block::ID_E1y0 ] = fBottom->GetXYZ( X, y );
553 pointsOnShapes[ SMESH_Block::ID_Ex01 ] = fTop->GetXYZ( x, 0 );
554 pointsOnShapes[ SMESH_Block::ID_Ex11 ] = fTop->GetXYZ( x, Y );
555 pointsOnShapes[ SMESH_Block::ID_E0y1 ] = fTop->GetXYZ( 0, y );
556 pointsOnShapes[ SMESH_Block::ID_E1y1 ] = fTop->GetXYZ( X, y );
557 // projection points on horizontal faces
558 pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = fBottom->GetXYZ( x, y );
559 pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = fTop ->GetXYZ( x, y );
560 for ( z = 1; z < zSize-1; ++z ) // z loop
562 params.SetCoord( 3, z / double(Z) );
563 // projection points on vertical edges
564 pointsOnShapes[ SMESH_Block::ID_E00z ] = fFront->GetXYZ( 0, z );
565 pointsOnShapes[ SMESH_Block::ID_E10z ] = fFront->GetXYZ( X, z );
566 pointsOnShapes[ SMESH_Block::ID_E01z ] = fBack->GetXYZ( 0, z );
567 pointsOnShapes[ SMESH_Block::ID_E11z ] = fBack->GetXYZ( X, z );
568 // projection points on vertical faces
569 pointsOnShapes[ SMESH_Block::ID_Fx0z ] = fFront->GetXYZ( x, z );
570 pointsOnShapes[ SMESH_Block::ID_Fx1z ] = fBack ->GetXYZ( x, z );
571 pointsOnShapes[ SMESH_Block::ID_F0yz ] = fLeft ->GetXYZ( y, z );
572 pointsOnShapes[ SMESH_Block::ID_F1yz ] = fRight->GetXYZ( y, z );
574 // compute internal node coordinates
576 SMESH_Block::ShellPoint( params, pointsOnShapes, coords );
577 column[ z ] = helper.AddNode( coords.X(), coords.Y(), coords.Z() );
583 // side data no more needed, free memory
584 for ( int i = 0; i < 6; ++i )
585 aCubeSide[i]._columns.clear();
587 // 5) Create hexahedrons
588 // ---------------------
590 for ( x = 0; x < xSize-1; ++x ) {
591 for ( y = 0; y < ySize-1; ++y ) {
592 vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
593 vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
594 vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
595 vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
596 for ( z = 0; z < zSize-1; ++z )
598 // bottom face normal of a hexa mush point outside the volume
599 helper.AddVolume(col00[z], col01[z], col11[z], col10[z],
600 col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
607 //=============================================================================
611 //=============================================================================
613 bool StdMeshers_Hexa_3D::Evaluate(SMESH_Mesh & aMesh,
614 const TopoDS_Shape & aShape,
615 MapShapeNbElems& aResMap)
617 vector < SMESH_subMesh * >meshFaces;
618 TopTools_SequenceOfShape aFaces;
619 for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
620 aFaces.Append(exp.Current());
621 SMESH_subMesh *aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
623 meshFaces.push_back(aSubMesh);
625 if (meshFaces.size() != 6) {
626 //return error(COMPERR_BAD_SHAPE, TComm(meshFaces.size())<<" instead of 6 faces in a block");
627 static StdMeshers_CompositeHexa_3D compositeHexa(-10, 0, aMesh.GetGen());
628 return compositeHexa.Evaluate(aMesh, aShape, aResMap);
633 //TopoDS_Shape aFace = meshFaces[i]->GetSubShape();
634 TopoDS_Shape aFace = aFaces.Value(i+1);
635 SMESH_Algo *algo = _gen->GetAlgo(aMesh, aFace);
637 std::vector<int> aResVec(SMDSEntity_Last);
638 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
639 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
640 aResMap.insert(std::make_pair(sm,aResVec));
641 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
642 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
645 string algoName = algo->GetName();
646 bool isAllQuad = false;
647 if (algoName == "Quadrangle_2D") {
648 MapShapeNbElemsItr anIt = aResMap.find(meshFaces[i]);
649 if( anIt == aResMap.end() ) continue;
650 std::vector<int> aVec = (*anIt).second;
651 int nbtri = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
656 return EvaluatePentahedralMesh(aMesh, aShape, aResMap);
660 // find number of 1d elems for 1 face
662 TopTools_MapOfShape Edges1;
663 bool IsQuadratic = false;
665 for (TopExp_Explorer exp(aFaces.Value(1), TopAbs_EDGE); exp.More(); exp.Next()) {
666 Edges1.Add(exp.Current());
667 SMESH_subMesh *sm = aMesh.GetSubMesh(exp.Current());
669 MapShapeNbElemsItr anIt = aResMap.find(sm);
670 if( anIt == aResMap.end() ) continue;
671 std::vector<int> aVec = (*anIt).second;
672 nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
674 IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
679 // find face opposite to 1 face
681 for(i=2; i<=6; i++) {
682 bool IsOpposite = true;
683 for(TopExp_Explorer exp(aFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) {
684 if( Edges1.Contains(exp.Current()) ) {
694 // find number of 2d elems on side faces
696 for(i=2; i<=6; i++) {
697 if( i == OppNum ) continue;
698 MapShapeNbElemsItr anIt = aResMap.find( meshFaces[i-1] );
699 if( anIt == aResMap.end() ) continue;
700 std::vector<int> aVec = (*anIt).second;
701 nb2d += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
704 MapShapeNbElemsItr anIt = aResMap.find( meshFaces[0] );
705 std::vector<int> aVec = (*anIt).second;
706 int nb2d_face0 = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
707 int nb0d_face0 = aVec[SMDSEntity_Node];
709 std::vector<int> aResVec(SMDSEntity_Last);
710 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
712 aResVec[SMDSEntity_Quad_Hexa] = nb2d_face0 * ( nb2d/nb1d );
713 int nb1d_face0_int = ( nb2d_face0*4 - nb1d ) / 2;
714 aResVec[SMDSEntity_Node] = nb0d_face0 * ( 2*nb2d/nb1d - 1 ) - nb1d_face0_int * nb2d/nb1d;
717 aResVec[SMDSEntity_Node] = nb0d_face0 * ( nb2d/nb1d - 1 );
718 aResVec[SMDSEntity_Hexa] = nb2d_face0 * ( nb2d/nb1d );
720 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
721 aResMap.insert(std::make_pair(sm,aResVec));
726 //================================================================================
728 * \brief Computes hexahedral mesh from 2D mesh of block
730 //================================================================================
732 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
734 static StdMeshers_HexaFromSkin_3D * algo = 0;
736 SMESH_Gen* gen = aMesh.GetGen();
737 algo = new StdMeshers_HexaFromSkin_3D( gen->GetANewId(), 0, gen );
739 algo->InitComputeError();
740 algo->Compute( aMesh, aHelper );
741 return error( algo->GetComputeError());
744 //=======================================================================
745 //function : ComputePentahedralMesh
747 //=======================================================================
749 SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh & aMesh,
750 const TopoDS_Shape & aShape,
751 SMESH_ProxyMesh* proxyMesh)
753 SMESH_ComputeErrorPtr err = SMESH_ComputeError::New();
756 err->myName = COMPERR_BAD_INPUT_MESH;
757 err->myComment = "Can't build pentahedral mesh on viscous layers";
761 StdMeshers_Penta_3D anAlgo;
763 bOK=anAlgo.Compute(aMesh, aShape);
765 err = anAlgo.GetComputeError();
767 if ( !bOK && anAlgo.ErrorStatus() == 5 )
769 static StdMeshers_Prism_3D * aPrism3D = 0;
771 SMESH_Gen* gen = aMesh.GetGen();
772 aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), 0, gen );
774 SMESH_Hypothesis::Hypothesis_Status aStatus;
775 if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
776 aPrism3D->InitComputeError();
777 bOK = aPrism3D->Compute( aMesh, aShape );
778 err = aPrism3D->GetComputeError();
785 //=======================================================================
786 //function : EvaluatePentahedralMesh
788 //=======================================================================
790 bool EvaluatePentahedralMesh(SMESH_Mesh & aMesh,
791 const TopoDS_Shape & aShape,
792 MapShapeNbElems& aResMap)
794 StdMeshers_Penta_3D anAlgo;
795 bool bOK = anAlgo.Evaluate(aMesh, aShape, aResMap);
797 //err = anAlgo.GetComputeError();
798 //if ( !bOK && anAlgo.ErrorStatus() == 5 )
800 static StdMeshers_Prism_3D * aPrism3D = 0;
802 SMESH_Gen* gen = aMesh.GetGen();
803 aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), 0, gen );
805 SMESH_Hypothesis::Hypothesis_Status aStatus;
806 if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
807 return aPrism3D->Evaluate(aMesh, aShape, aResMap);