1 // Copyright (C) 2007-2020 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, or (at your option) any later version.
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 : implementation 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, SMESH_Gen * gen)
74 :SMESH_3D_Algo(hypId, gen)
77 _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID); // 1 bit /shape type
78 _requireShape = false;
79 _compatibleHypothesis.push_back("ViscousLayers");
80 _quadAlgo = new StdMeshers_Quadrangle_2D( gen->GetANewId(), _gen );
83 //=============================================================================
87 //=============================================================================
89 StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D()
95 //=============================================================================
97 * Retrieves defined hypotheses
99 //=============================================================================
101 bool StdMeshers_Hexa_3D::CheckHypothesis
103 const TopoDS_Shape& aShape,
104 SMESH_Hypothesis::Hypothesis_Status& aStatus)
106 // check nb of faces in the shape
108 aStatus = SMESH_Hypothesis::HYP_BAD_GEOMETRY;
110 for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next())
117 _viscousLayersHyp = NULL;
119 const list<const SMESHDS_Hypothesis*>& hyps =
120 GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
121 list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
122 if ( h == hyps.end())
124 aStatus = SMESH_Hypothesis::HYP_OK;
128 // only StdMeshers_ViscousLayers can be used
130 for ( ; h != hyps.end(); ++h )
132 if ( !(_viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h )))
135 if ( !_viscousLayersHyp )
136 aStatus = HYP_INCOMPATIBLE;
138 error( _viscousLayersHyp->CheckHypothesis( aMesh, aShape, aStatus ));
140 return aStatus == HYP_OK;
145 //=============================================================================
147 typedef boost::shared_ptr< FaceQuadStruct > FaceQuadStructPtr;
148 typedef std::vector<gp_XYZ> TXYZColumn;
150 // symbolic names of box sides
151 enum EBoxSides{ B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_NB_SIDES };
153 // symbolic names of sides of quadrangle
154 enum EQuadSides{ Q_BOTTOM=0, Q_RIGHT, Q_TOP, Q_LEFT, Q_NB_SIDES };
156 enum EAxes{ COO_X=1, COO_Y, COO_Z };
158 //=============================================================================
160 * \brief Container of nodes of structured mesh on a qudrangular geom FACE
165 FaceQuadStructPtr _quad;
167 // map of (node parameter on EDGE) to (column (vector) of nodes)
168 TParam2ColumnMap _u2nodesMap;
170 // node column's taken from _u2nodesMap taking into account sub-shape orientation
171 vector<TNodeColumn> _columns;
173 // columns of normalized parameters of nodes within the unitary cube
174 vector<TXYZColumn> _ijkColumns;
176 // geometry of a cube side
179 const SMDS_MeshNode* GetNode(int iCol, int iRow) const
181 return _columns[iCol][iRow];
183 gp_XYZ GetXYZ(int iCol, int iRow) const
185 return SMESH_TNodeXYZ( GetNode( iCol, iRow ));
187 gp_XYZ& GetIJK(int iCol, int iRow)
189 return _ijkColumns[iCol][iRow];
193 //================================================================================
195 * \brief Converter of a pair of integers to a sole index
200 _Indexer( int xSize, int ySize ): _xSize(xSize), _ySize(ySize) {}
201 int size() const { return _xSize * _ySize; }
202 int operator()(const int x, const int y) const { return y * _xSize + x; }
205 //================================================================================
207 * \brief Appends a range of node columns from a map to another map
209 template< class TMapIterator >
210 void append( TParam2ColumnMap& toMap, TMapIterator from, TMapIterator to )
212 const SMDS_MeshNode* lastNode = toMap.rbegin()->second[0];
213 const SMDS_MeshNode* firstNode = from->second[0];
214 if ( lastNode == firstNode )
216 double u = toMap.rbegin()->first;
217 for (; from != to; ++from )
220 TParam2ColumnMap::iterator u2nn = toMap.insert( toMap.end(), make_pair ( u, TNodeColumn()));
221 u2nn->second.swap( from->second );
225 //================================================================================
227 * \brief Finds FaceQuadStruct having a side equal to a given one and rearranges
228 * the found FaceQuadStruct::side to have the given side at a Q_BOTTOM place
230 FaceQuadStructPtr getQuadWithBottom( StdMeshers_FaceSidePtr side,
231 FaceQuadStructPtr quad[ 6 ])
233 FaceQuadStructPtr foundQuad;
234 for ( int i = 1; i < 6; ++i )
236 if ( !quad[i] ) continue;
237 for ( size_t iS = 0; iS < quad[i]->side.size(); ++iS )
239 const StdMeshers_FaceSidePtr side2 = quad[i]->side[iS];
240 if (( side->FirstVertex().IsSame( side2->FirstVertex() ) ||
241 side->FirstVertex().IsSame( side2->LastVertex() ))
243 ( side->LastVertex().IsSame( side2->FirstVertex() ) ||
244 side->LastVertex().IsSame( side2->LastVertex() ))
247 if ( iS != Q_BOTTOM )
249 vector< FaceQuadStruct::Side > newSides;
250 for ( size_t j = iS; j < quad[i]->side.size(); ++j )
251 newSides.push_back( quad[i]->side[j] );
252 for ( size_t j = 0; j < iS; ++j )
253 newSides.push_back( quad[i]->side[j] );
254 quad[i]->side.swap( newSides );
256 foundQuad.swap(quad[i]);
263 //================================================================================
265 * \brief Returns true if the 1st base node of sideGrid1 belongs to sideGrid2
267 //================================================================================
269 bool beginsAtSide( const _FaceGrid& sideGrid1,
270 const _FaceGrid& sideGrid2,
271 SMESH_ProxyMesh::Ptr proxymesh )
273 const TNodeColumn& col0 = sideGrid2._u2nodesMap.begin()->second;
274 const TNodeColumn& col1 = sideGrid2._u2nodesMap.rbegin()->second;
275 const SMDS_MeshNode* n00 = col0.front();
276 const SMDS_MeshNode* n01 = col0.back();
277 const SMDS_MeshNode* n10 = col1.front();
278 const SMDS_MeshNode* n11 = col1.back();
279 const SMDS_MeshNode* n = (sideGrid1._u2nodesMap.begin()->second)[0];
282 n00 = proxymesh->GetProxyNode( n00 );
283 n10 = proxymesh->GetProxyNode( n10 );
284 n01 = proxymesh->GetProxyNode( n01 );
285 n11 = proxymesh->GetProxyNode( n11 );
286 n = proxymesh->GetProxyNode( n );
288 return ( n == n00 || n == n01 || n == n10 || n == n11 );
291 //================================================================================
293 * \brief Fill in _FaceGrid::_ijkColumns
294 * \param [in,out] fg - a _FaceGrid
295 * \param [in] i1 - coordinate index along _columns
296 * \param [in] i2 - coordinate index along _columns[i]
297 * \param [in] v3 - value of the constant parameter
299 //================================================================================
301 void computeIJK( _FaceGrid& fg, int i1, int i2, double v3 )
303 gp_XYZ ijk( v3, v3, v3 );
304 const size_t nbCol = fg._columns.size();
305 const size_t nbRow = fg._columns[0].size();
307 fg._ijkColumns.resize( nbCol );
308 for ( size_t i = 0; i < nbCol; ++i )
309 fg._ijkColumns[ i ].resize( nbRow, ijk );
311 vector< double > len( nbRow );
313 for ( size_t i = 0; i < nbCol; ++i )
315 gp_Pnt pPrev = fg.GetXYZ( i, 0 );
316 for ( size_t j = 1; j < nbRow; ++j )
318 gp_Pnt p = fg.GetXYZ( i, j );
319 len[ j ] = len[ j-1 ] + p.Distance( pPrev );
322 for ( size_t j = 0; j < nbRow; ++j )
323 fg.GetIJK( i, j ).SetCoord( i2, len[ j ]/len.back() );
327 for ( size_t j = 0; j < nbRow; ++j )
329 gp_Pnt pPrev = fg.GetXYZ( 0, j );
330 for ( size_t i = 1; i < nbCol; ++i )
332 gp_Pnt p = fg.GetXYZ( i, j );
333 len[ i ] = len[ i-1 ] + p.Distance( pPrev );
336 for ( size_t i = 0; i < nbCol; ++i )
337 fg.GetIJK( i, j ).SetCoord( i1, len[ i ]/len.back() );
342 //=============================================================================
344 * Generates hexahedron mesh on hexaedron like form using algorithm from
345 * "Application de l'interpolation transfinie � la cr�ation de maillages
346 * C0 ou G1 continus sur des triangles, quadrangles, tetraedres, pentaedres
347 * et hexaedres d�form�s."
348 * Alain PERONNET - 8 janvier 1999
350 //=============================================================================
352 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
353 const TopoDS_Shape & aShape)
355 // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside
356 //Unexpect aCatch(SalomeException);
357 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
359 // Shape verification
360 // ----------------------
362 // shape must be a solid (or a shell) with 6 faces
363 TopExp_Explorer exp(aShape,TopAbs_SHELL);
365 return error(COMPERR_BAD_SHAPE, "No SHELL in the geometry");
366 if ( exp.Next(), exp.More() )
367 return error(COMPERR_BAD_SHAPE, "More than one SHELL in the geometry");
369 TopTools_IndexedMapOfShape FF, EE;
370 TopExp::MapShapes( aShape, TopAbs_FACE, FF);
371 if ( FF.Extent() != 6)
373 static StdMeshers_CompositeHexa_3D compositeHexa(_gen->GetANewId(), _gen);
374 if ( !compositeHexa.Compute( aMesh, aShape ))
375 return error( compositeHexa.GetComputeError() );
379 // Find sides of a cube
380 // ---------------------
382 // tool creating quadratic elements if needed
383 SMESH_MesherHelper helper (aMesh);
384 _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
386 TopExp::MapShapes( aShape, TopAbs_EDGE, EE );
387 SMESH_MesherHelper* faceHelper = ( EE.Size() == 12 ) ? 0 : &helper;
389 FaceQuadStructPtr quad[ 6 ];
390 for ( int i = 0; i < 6; ++i )
393 faceHelper->SetSubShape( FF( i+1 ));
394 if ( !( quad[i] = FaceQuadStructPtr( _quadAlgo->CheckNbEdges( aMesh, FF( i+1 ),
395 /*considerMesh=*/true,
397 return error( _quadAlgo->GetComputeError() );
398 if ( quad[i]->side.size() != 4 )
399 return error( COMPERR_BAD_SHAPE, "Not a quadrangular box side" );
402 _FaceGrid aCubeSide[ 6 ];
404 swap( aCubeSide[B_BOTTOM]._quad, quad[0] );
405 swap( aCubeSide[B_BOTTOM]._quad->side[ Q_RIGHT],// direct the normal of bottom quad inside cube
406 aCubeSide[B_BOTTOM]._quad->side[ Q_LEFT ] );
408 aCubeSide[B_FRONT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_BOTTOM], quad );
409 aCubeSide[B_RIGHT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_RIGHT ], quad );
410 aCubeSide[B_BACK ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_TOP ], quad );
411 aCubeSide[B_LEFT ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_LEFT ], quad );
412 if ( aCubeSide[B_FRONT ]._quad )
413 aCubeSide[B_TOP]._quad = getQuadWithBottom( aCubeSide[B_FRONT ]._quad->side[Q_TOP ], quad );
415 for ( int i = 1; i < 6; ++i )
416 if ( !aCubeSide[i]._quad )
417 return error( COMPERR_BAD_SHAPE );
419 // Make viscous layers
420 // --------------------
422 SMESH_ProxyMesh::Ptr proxymesh;
423 if ( _viscousLayersHyp )
425 proxymesh = _viscousLayersHyp->Compute( aMesh, aShape, /*makeN2NMap=*/ true );
430 // Check if there are triangles on cube sides
431 // -------------------------------------------
433 if ( aMesh.NbTriangles() > 0 )
435 for ( int i = 0; i < 6; ++i )
437 const TopoDS_Face& sideF = aCubeSide[i]._quad->face;
438 const SMESHDS_SubMesh* smDS =
439 proxymesh ? proxymesh->GetSubMesh( sideF ) : meshDS->MeshElements( sideF );
440 if ( !SMESH_MesherHelper::IsSameElemGeometry( smDS, SMDSGeom_QUADRANGLE,
441 /*nullSubMeshRes=*/false ))
443 SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
449 // Check presence of regular grid mesh on FACEs of the cube
450 // ------------------------------------------------------------
452 for ( int i = 0; i < 6; ++i )
454 const TopoDS_Face& F = aCubeSide[i]._quad->face;
455 StdMeshers_FaceSidePtr baseQuadSide = aCubeSide[i]._quad->side[ Q_BOTTOM ];
456 list<TopoDS_Edge> baseEdges( baseQuadSide->Edges().begin(), baseQuadSide->Edges().end() );
458 // assure correctness of node positions on baseE:
459 // helper.GetNodeU() will fix positions if they are wrong
460 helper.ToFixNodeParameters( true );
461 for ( int iE = 0; iE < baseQuadSide->NbEdges(); ++iE )
463 const TopoDS_Edge& baseE = baseQuadSide->Edge( iE );
464 if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( baseE ))
467 helper.SetSubShape( baseE );
468 SMDS_ElemIteratorPtr eIt = smDS->GetElements();
469 while ( eIt->more() )
471 const SMDS_MeshElement* e = eIt->next();
472 // expect problems on a composite side
473 try { helper.GetNodeU( baseE, e->GetNode(0), e->GetNode(1), &ok); }
475 try { helper.GetNodeU( baseE, e->GetNode(1), e->GetNode(0), &ok); }
483 helper.LoadNodeColumns( aCubeSide[i]._u2nodesMap, F, baseEdges, meshDS, proxymesh.get());
486 // check if the loaded grid corresponds to nb of quadrangles on the FACE
487 const SMESHDS_SubMesh* faceSubMesh =
488 proxymesh ? proxymesh->GetSubMesh( F ) : meshDS->MeshElements( F );
489 const int nbQuads = faceSubMesh->NbElements();
490 const int nbHor = aCubeSide[i]._u2nodesMap.size() - 1;
491 const int nbVer = aCubeSide[i]._u2nodesMap.begin()->second.size() - 1;
492 ok = ( nbQuads == nbHor * nbVer );
496 SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
501 // Orient loaded grids of cube sides along axis of the unitary cube coord system
503 isReverse[B_BOTTOM] = beginsAtSide( aCubeSide[B_BOTTOM], aCubeSide[B_RIGHT ], proxymesh );
504 isReverse[B_TOP ] = beginsAtSide( aCubeSide[B_TOP ], aCubeSide[B_RIGHT ], proxymesh );
505 isReverse[B_FRONT ] = beginsAtSide( aCubeSide[B_FRONT ], aCubeSide[B_RIGHT ], proxymesh );
506 isReverse[B_BACK ] = beginsAtSide( aCubeSide[B_BACK ], aCubeSide[B_RIGHT ], proxymesh );
507 isReverse[B_LEFT ] = beginsAtSide( aCubeSide[B_LEFT ], aCubeSide[B_BACK ], proxymesh );
508 isReverse[B_RIGHT ] = beginsAtSide( aCubeSide[B_RIGHT ], aCubeSide[B_BACK ], proxymesh );
509 for ( int i = 0; i < 6; ++i )
511 aCubeSide[i]._columns.resize( aCubeSide[i]._u2nodesMap.size() );
513 size_t iFwd = 0, iRev = aCubeSide[i]._columns.size()-1;
514 size_t* pi = isReverse[i] ? &iRev : &iFwd;
515 TParam2ColumnMap::iterator u2nn = aCubeSide[i]._u2nodesMap.begin();
516 for ( ; iFwd < aCubeSide[i]._columns.size(); --iRev, ++iFwd, ++u2nn )
517 aCubeSide[i]._columns[ *pi ].swap( u2nn->second );
519 aCubeSide[i]._u2nodesMap.clear();
523 for ( int i = 0; i < 6; ++i )
524 for ( size_t j = 0; j < aCubeSide[i]._columns.size(); ++j)
525 for ( size_t k = 0; k < aCubeSide[i]._columns[j].size(); ++k)
527 const SMDS_MeshNode* & n = aCubeSide[i]._columns[j][k];
528 n = proxymesh->GetProxyNode( n );
531 // 4) Create internal nodes of the cube
532 // -------------------------------------
534 helper.SetSubShape( aShape );
535 helper.SetElementsOnShape(true);
537 // shortcuts to sides
538 _FaceGrid* fBottom = & aCubeSide[ B_BOTTOM ];
539 _FaceGrid* fRight = & aCubeSide[ B_RIGHT ];
540 _FaceGrid* fTop = & aCubeSide[ B_TOP ];
541 _FaceGrid* fLeft = & aCubeSide[ B_LEFT ];
542 _FaceGrid* fFront = & aCubeSide[ B_FRONT ];
543 _FaceGrid* fBack = & aCubeSide[ B_BACK ];
545 // cube size measured in nb of nodes
546 size_t x, xSize = fBottom->_columns.size() , X = xSize - 1;
547 size_t y, ySize = fLeft->_columns.size() , Y = ySize - 1;
548 size_t z, zSize = fLeft->_columns[0].size(), Z = zSize - 1;
550 // check sharing of FACEs (IPAL54417)
551 if ( fFront ->_columns.size() != xSize ||
552 fBack ->_columns.size() != xSize ||
553 fTop ->_columns.size() != xSize ||
555 fRight ->_columns.size() != ySize ||
556 fTop ->_columns[0].size() != ySize ||
557 fBottom->_columns[0].size() != ySize ||
559 fRight ->_columns[0].size() != zSize ||
560 fFront ->_columns[0].size() != zSize ||
561 fBack ->_columns[0].size() != zSize )
562 return error( COMPERR_BAD_SHAPE, "Not sewed faces" );
564 // columns of internal nodes "rising" from nodes of fBottom
565 _Indexer colIndex( xSize, ySize );
566 vector< vector< const SMDS_MeshNode* > > columns( colIndex.size() );
568 // fill node columns by front and back box sides
569 for ( x = 0; x < xSize; ++x ) {
570 vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
571 vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
572 column0.resize( zSize );
573 column1.resize( zSize );
574 for ( z = 0; z < zSize; ++z ) {
575 column0[ z ] = fFront->GetNode( x, z );
576 column1[ z ] = fBack ->GetNode( x, z );
579 // fill node columns by left and right box sides
580 for ( y = 1; y < ySize-1; ++y ) {
581 vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
582 vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
583 column0.resize( zSize );
584 column1.resize( zSize );
585 for ( z = 0; z < zSize; ++z ) {
586 column0[ z ] = fLeft ->GetNode( y, z );
587 column1[ z ] = fRight->GetNode( y, z );
590 // get nodes from top and bottom box sides
591 for ( x = 1; x < xSize-1; ++x ) {
592 for ( y = 1; y < ySize-1; ++y ) {
593 vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
594 column.resize( zSize );
595 column.front() = fBottom->GetNode( x, y );
596 column.back() = fTop ->GetNode( x, y );
600 // compute normalized parameters of nodes on sides (PAL23189)
601 computeIJK( *fBottom, COO_X, COO_Y, /*z=*/0. );
602 computeIJK( *fRight, COO_Y, COO_Z, /*x=*/1. );
603 computeIJK( *fTop, COO_X, COO_Y, /*z=*/1. );
604 computeIJK( *fLeft, COO_Y, COO_Z, /*x=*/0. );
605 computeIJK( *fFront, COO_X, COO_Z, /*y=*/0. );
606 computeIJK( *fBack, COO_X, COO_Z, /*y=*/1. );
608 // projection points of the internal node on cube sub-shapes by which
609 // coordinates of the internal node are computed
610 vector<gp_XYZ> pointsOnShapes( SMESH_Block::ID_Shell );
612 // projections on vertices are constant
613 pointsOnShapes[ SMESH_Block::ID_V000 ] = fBottom->GetXYZ( 0, 0 );
614 pointsOnShapes[ SMESH_Block::ID_V100 ] = fBottom->GetXYZ( X, 0 );
615 pointsOnShapes[ SMESH_Block::ID_V010 ] = fBottom->GetXYZ( 0, Y );
616 pointsOnShapes[ SMESH_Block::ID_V110 ] = fBottom->GetXYZ( X, Y );
617 pointsOnShapes[ SMESH_Block::ID_V001 ] = fTop->GetXYZ( 0, 0 );
618 pointsOnShapes[ SMESH_Block::ID_V101 ] = fTop->GetXYZ( X, 0 );
619 pointsOnShapes[ SMESH_Block::ID_V011 ] = fTop->GetXYZ( 0, Y );
620 pointsOnShapes[ SMESH_Block::ID_V111 ] = fTop->GetXYZ( X, Y );
622 gp_XYZ params; // normalized parameters of an internal node within the unit box
623 for ( x = 1; x < xSize-1; ++x )
625 const double rX = x / double(X);
626 for ( y = 1; y < ySize-1; ++y )
628 const double rY = y / double(Y);
630 // a column to fill in during z loop
631 vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
632 // projection points on horizontal edges
633 pointsOnShapes[ SMESH_Block::ID_Ex00 ] = fBottom->GetXYZ( x, 0 );
634 pointsOnShapes[ SMESH_Block::ID_Ex10 ] = fBottom->GetXYZ( x, Y );
635 pointsOnShapes[ SMESH_Block::ID_E0y0 ] = fBottom->GetXYZ( 0, y );
636 pointsOnShapes[ SMESH_Block::ID_E1y0 ] = fBottom->GetXYZ( X, y );
637 pointsOnShapes[ SMESH_Block::ID_Ex01 ] = fTop->GetXYZ( x, 0 );
638 pointsOnShapes[ SMESH_Block::ID_Ex11 ] = fTop->GetXYZ( x, Y );
639 pointsOnShapes[ SMESH_Block::ID_E0y1 ] = fTop->GetXYZ( 0, y );
640 pointsOnShapes[ SMESH_Block::ID_E1y1 ] = fTop->GetXYZ( X, y );
641 // projection points on horizontal faces
642 pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = fBottom->GetXYZ( x, y );
643 pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = fTop ->GetXYZ( x, y );
644 for ( z = 1; z < zSize-1; ++z ) // z loop
646 const double rZ = z / double(Z);
648 const gp_XYZ& pBo = fBottom->GetIJK( x, y );
649 const gp_XYZ& pTo = fTop ->GetIJK( x, y );
650 const gp_XYZ& pFr = fFront ->GetIJK( x, z );
651 const gp_XYZ& pBa = fBack ->GetIJK( x, z );
652 const gp_XYZ& pLe = fLeft ->GetIJK( y, z );
653 const gp_XYZ& pRi = fRight ->GetIJK( y, z );
654 params.SetCoord( 1, 0.5 * ( pBo.X() * ( 1. - rZ ) + pTo.X() * rZ +
655 pFr.X() * ( 1. - rY ) + pBa.X() * rY ));
656 params.SetCoord( 2, 0.5 * ( pBo.Y() * ( 1. - rZ ) + pTo.Y() * rZ +
657 pLe.Y() * ( 1. - rX ) + pRi.Y() * rX ));
658 params.SetCoord( 3, 0.5 * ( pFr.Z() * ( 1. - rY ) + pBa.Z() * rY +
659 pLe.Z() * ( 1. - rX ) + pRi.Z() * rX ));
661 // projection points on vertical edges
662 pointsOnShapes[ SMESH_Block::ID_E00z ] = fFront->GetXYZ( 0, z );
663 pointsOnShapes[ SMESH_Block::ID_E10z ] = fFront->GetXYZ( X, z );
664 pointsOnShapes[ SMESH_Block::ID_E01z ] = fBack->GetXYZ( 0, z );
665 pointsOnShapes[ SMESH_Block::ID_E11z ] = fBack->GetXYZ( X, z );
666 // projection points on vertical faces
667 pointsOnShapes[ SMESH_Block::ID_Fx0z ] = fFront->GetXYZ( x, z );
668 pointsOnShapes[ SMESH_Block::ID_Fx1z ] = fBack ->GetXYZ( x, z );
669 pointsOnShapes[ SMESH_Block::ID_F0yz ] = fLeft ->GetXYZ( y, z );
670 pointsOnShapes[ SMESH_Block::ID_F1yz ] = fRight->GetXYZ( y, z );
672 // compute internal node coordinates
674 SMESH_Block::ShellPoint( params, pointsOnShapes, coords );
675 column[ z ] = helper.AddNode( coords.X(), coords.Y(), coords.Z() );
680 // side data no more needed, free memory
681 for ( int i = 0; i < 6; ++i )
682 aCubeSide[i]._columns.clear();
684 // 5) Create hexahedrons
685 // ---------------------
687 for ( x = 0; x < xSize-1; ++x ) {
688 for ( y = 0; y < ySize-1; ++y ) {
689 vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
690 vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
691 vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
692 vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
693 for ( z = 0; z < zSize-1; ++z )
695 // bottom face normal of a hexa mush point outside the volume
696 helper.AddVolume(col00[z], col01[z], col11[z], col10[z],
697 col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
704 //=============================================================================
708 //=============================================================================
710 bool StdMeshers_Hexa_3D::Evaluate(SMESH_Mesh & aMesh,
711 const TopoDS_Shape & aShape,
712 MapShapeNbElems& aResMap)
714 vector < SMESH_subMesh * >meshFaces;
715 TopTools_SequenceOfShape aFaces;
716 for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
717 aFaces.Append(exp.Current());
718 SMESH_subMesh *aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
720 meshFaces.push_back(aSubMesh);
722 if (meshFaces.size() != 6) {
723 //return error(COMPERR_BAD_SHAPE, TComm(meshFaces.size())<<" instead of 6 faces in a block");
724 static StdMeshers_CompositeHexa_3D compositeHexa(-10, aMesh.GetGen());
725 return compositeHexa.Evaluate(aMesh, aShape, aResMap);
730 //TopoDS_Shape aFace = meshFaces[i]->GetSubShape();
731 TopoDS_Shape aFace = aFaces.Value(i+1);
732 SMESH_Algo *algo = _gen->GetAlgo(aMesh, aFace);
734 std::vector<int> aResVec(SMDSEntity_Last);
735 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
736 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
737 aResMap.insert(std::make_pair(sm,aResVec));
738 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
739 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
742 string algoName = algo->GetName();
743 bool isAllQuad = false;
744 if (algoName == "Quadrangle_2D") {
745 MapShapeNbElemsItr anIt = aResMap.find(meshFaces[i]);
746 if( anIt == aResMap.end() ) continue;
747 std::vector<int> aVec = (*anIt).second;
748 int nbtri = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
753 return EvaluatePentahedralMesh(aMesh, aShape, aResMap);
757 // find number of 1d elems for 1 face
759 TopTools_MapOfShape Edges1;
760 bool IsQuadratic = false;
762 for (TopExp_Explorer exp(aFaces.Value(1), TopAbs_EDGE); exp.More(); exp.Next()) {
763 Edges1.Add(exp.Current());
764 SMESH_subMesh *sm = aMesh.GetSubMesh(exp.Current());
766 MapShapeNbElemsItr anIt = aResMap.find(sm);
767 if( anIt == aResMap.end() ) continue;
768 std::vector<int> aVec = (*anIt).second;
769 nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
771 IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
776 // find face opposite to 1 face
778 for(i=2; i<=6; i++) {
779 bool IsOpposite = true;
780 for(TopExp_Explorer exp(aFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) {
781 if( Edges1.Contains(exp.Current()) ) {
791 // find number of 2d elems on side faces
793 for(i=2; i<=6; i++) {
794 if( i == OppNum ) continue;
795 MapShapeNbElemsItr anIt = aResMap.find( meshFaces[i-1] );
796 if( anIt == aResMap.end() ) continue;
797 std::vector<int> aVec = (*anIt).second;
798 nb2d += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
801 MapShapeNbElemsItr anIt = aResMap.find( meshFaces[0] );
802 std::vector<int> aVec = (*anIt).second;
803 int nb2d_face0 = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
804 int nb0d_face0 = aVec[SMDSEntity_Node];
806 std::vector<int> aResVec(SMDSEntity_Last);
807 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
809 aResVec[SMDSEntity_Quad_Hexa] = nb2d_face0 * ( nb2d/nb1d );
810 int nb1d_face0_int = ( nb2d_face0*4 - nb1d ) / 2;
811 aResVec[SMDSEntity_Node] = nb0d_face0 * ( 2*nb2d/nb1d - 1 ) - nb1d_face0_int * nb2d/nb1d;
814 aResVec[SMDSEntity_Node] = nb0d_face0 * ( nb2d/nb1d - 1 );
815 aResVec[SMDSEntity_Hexa] = nb2d_face0 * ( nb2d/nb1d );
817 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
818 aResMap.insert(std::make_pair(sm,aResVec));
823 //================================================================================
825 * \brief Computes hexahedral mesh from 2D mesh of block
827 //================================================================================
829 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
831 static StdMeshers_HexaFromSkin_3D * algo = 0;
833 SMESH_Gen* gen = aMesh.GetGen();
834 algo = new StdMeshers_HexaFromSkin_3D( gen->GetANewId(), gen );
836 algo->InitComputeError();
837 algo->Compute( aMesh, aHelper );
838 return error( algo->GetComputeError());
841 //================================================================================
843 * \brief Return true if the algorithm can mesh this shape
844 * \param [in] aShape - shape to check
845 * \param [in] toCheckAll - if true, this check returns OK if all shapes are OK,
846 * else, returns OK if at least one shape is OK
848 //================================================================================
850 bool StdMeshers_Hexa_3D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
852 TopExp_Explorer exp0( aShape, TopAbs_SOLID );
853 if ( !exp0.More() ) return false;
855 for ( ; exp0.More(); exp0.Next() )
857 int nbFoundShells = 0;
858 TopExp_Explorer exp1( exp0.Current(), TopAbs_SHELL );
859 for ( ; exp1.More(); exp1.Next(), ++nbFoundShells)
860 if ( nbFoundShells == 2 ) break;
861 if ( nbFoundShells != 1 ) {
862 if ( toCheckAll ) return false;
865 exp1.Init( exp0.Current(), TopAbs_FACE );
866 int nbEdges = SMESH_MesherHelper::Count( exp1.Current(), TopAbs_EDGE, /*ignoreSame=*/true );
867 bool ok = ( nbEdges > 3 );
868 if ( toCheckAll && !ok ) return false;
869 if ( !toCheckAll && ok ) return true;
874 //=======================================================================
875 //function : ComputePentahedralMesh
877 //=======================================================================
879 SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh & aMesh,
880 const TopoDS_Shape & aShape,
881 SMESH_ProxyMesh* proxyMesh)
883 SMESH_ComputeErrorPtr err = SMESH_ComputeError::New();
886 err->myName = COMPERR_BAD_INPUT_MESH;
887 err->myComment = "Can't build pentahedral mesh on viscous layers";
891 StdMeshers_Penta_3D anAlgo;
893 bOK=anAlgo.Compute(aMesh, aShape);
895 err = anAlgo.GetComputeError();
897 if ( !bOK && anAlgo.ErrorStatus() == 5 )
899 static StdMeshers_Prism_3D * aPrism3D = 0;
901 SMESH_Gen* gen = aMesh.GetGen();
902 aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), gen );
904 SMESH_Hypothesis::Hypothesis_Status aStatus;
905 if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
906 aPrism3D->InitComputeError();
907 bOK = aPrism3D->Compute( aMesh, aShape );
908 err = aPrism3D->GetComputeError();
915 //=======================================================================
916 //function : EvaluatePentahedralMesh
918 //=======================================================================
920 bool EvaluatePentahedralMesh(SMESH_Mesh & aMesh,
921 const TopoDS_Shape & aShape,
922 MapShapeNbElems& aResMap)
924 StdMeshers_Penta_3D anAlgo;
925 bool bOK = anAlgo.Evaluate(aMesh, aShape, aResMap);
927 //err = anAlgo.GetComputeError();
928 //if ( !bOK && anAlgo.ErrorStatus() == 5 )
930 static StdMeshers_Prism_3D * aPrism3D = 0;
932 SMESH_Gen* gen = aMesh.GetGen();
933 aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), gen );
935 SMESH_Hypothesis::Hypothesis_Status aStatus;
936 if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
937 return aPrism3D->Evaluate(aMesh, aShape, aResMap);