X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Hexa_3D.cxx;h=edfca5ae541f8aa0b87da28646bd615836a03dd8;hp=1ad459448f597df5017c01702262aefa2abc70ff;hb=6b45e1948462a8feb8618ea810ff2ec7fe00b6ea;hpb=6bac08c1a81f34d3f21c550bd92f83654b2546a5 diff --git a/src/StdMeshers/StdMeshers_Hexa_3D.cxx b/src/StdMeshers/StdMeshers_Hexa_3D.cxx index 1ad459448..edfca5ae5 100644 --- a/src/StdMeshers/StdMeshers_Hexa_3D.cxx +++ b/src/StdMeshers/StdMeshers_Hexa_3D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2019 CEA/DEN, EDF R&D, OPEN CASCADE // // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS @@ -20,7 +20,7 @@ // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -// SMESH SMESH : implementaion of SMESH idl descriptions +// SMESH SMESH : implementation of SMESH idl descriptions // File : StdMeshers_Hexa_3D.cxx // Moved here from SMESH_Hexa_3D.cxx // Author : Paul RASCLE, EDF @@ -70,14 +70,14 @@ static bool EvaluatePentahedralMesh(SMESH_Mesh &, const TopoDS_Shape &, */ //============================================================================= -StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, int studyId, SMESH_Gen * gen) - :SMESH_3D_Algo(hypId, studyId, gen) +StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, SMESH_Gen * gen) + :SMESH_3D_Algo(hypId, gen) { - MESSAGE("StdMeshers_Hexa_3D::StdMeshers_Hexa_3D"); _name = "Hexa_3D"; _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID); // 1 bit /shape type _requireShape = false; _compatibleHypothesis.push_back("ViscousLayers"); + _quadAlgo = new StdMeshers_Quadrangle_2D( gen->GetANewId(), _gen ); } //============================================================================= @@ -88,7 +88,8 @@ StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, int studyId, SMESH_Gen * gen) StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D() { - MESSAGE("StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D"); + delete _quadAlgo; + _quadAlgo = 0; } //============================================================================= @@ -144,6 +145,7 @@ namespace //============================================================================= typedef boost::shared_ptr< FaceQuadStruct > FaceQuadStructPtr; + typedef std::vector TXYZColumn; // symbolic names of box sides enum EBoxSides{ B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_NB_SIDES }; @@ -151,6 +153,8 @@ namespace // symbolic names of sides of quadrangle enum EQuadSides{ Q_BOTTOM=0, Q_RIGHT, Q_TOP, Q_LEFT, Q_NB_SIDES }; + enum EAxes{ COO_X=1, COO_Y, COO_Z }; + //============================================================================= /*! * \brief Container of nodes of structured mesh on a qudrangular geom FACE @@ -163,9 +167,12 @@ namespace // map of (node parameter on EDGE) to (column (vector) of nodes) TParam2ColumnMap _u2nodesMap; - // node column's taken form _u2nodesMap taking into account sub-shape orientation + // node column's taken from _u2nodesMap taking into account sub-shape orientation vector _columns; + // columns of normalized parameters of nodes within the unitary cube + vector _ijkColumns; + // geometry of a cube side TopoDS_Face _sideF; @@ -177,11 +184,15 @@ namespace { return SMESH_TNodeXYZ( GetNode( iCol, iRow )); } + gp_XYZ& GetIJK(int iCol, int iRow) + { + return _ijkColumns[iCol][iRow]; + } }; //================================================================================ /*! - * \brief Convertor of a pair of integers to a sole index + * \brief Converter of a pair of integers to a sole index */ struct _Indexer { @@ -223,7 +234,7 @@ namespace for ( int i = 1; i < 6; ++i ) { if ( !quad[i] ) continue; - for ( unsigned iS = 0; iS < quad[i]->side.size(); ++iS ) + for ( size_t iS = 0; iS < quad[i]->side.size(); ++iS ) { const StdMeshers_FaceSidePtr side2 = quad[i]->side[iS]; if (( side->FirstVertex().IsSame( side2->FirstVertex() ) || @@ -236,9 +247,9 @@ namespace if ( iS != Q_BOTTOM ) { vector< FaceQuadStruct::Side > newSides; - for ( unsigned j = iS; j < quad[i]->side.size(); ++j ) + for ( size_t j = iS; j < quad[i]->side.size(); ++j ) newSides.push_back( quad[i]->side[j] ); - for ( unsigned j = 0; j < iS; ++j ) + for ( size_t j = 0; j < iS; ++j ) newSides.push_back( quad[i]->side[j] ); quad[i]->side.swap( newSides ); } @@ -276,6 +287,56 @@ namespace } return ( n == n00 || n == n01 || n == n10 || n == n11 ); } + + //================================================================================ + /*! + * \brief Fill in _FaceGrid::_ijkColumns + * \param [in,out] fg - a _FaceGrid + * \param [in] i1 - coordinate index along _columns + * \param [in] i2 - coordinate index along _columns[i] + * \param [in] v3 - value of the constant parameter + */ + //================================================================================ + + void computeIJK( _FaceGrid& fg, int i1, int i2, double v3 ) + { + gp_XYZ ijk( v3, v3, v3 ); + const size_t nbCol = fg._columns.size(); + const size_t nbRow = fg._columns[0].size(); + + fg._ijkColumns.resize( nbCol ); + for ( size_t i = 0; i < nbCol; ++i ) + fg._ijkColumns[ i ].resize( nbRow, ijk ); + + vector< double > len( nbRow ); + len[0] = 0; + for ( size_t i = 0; i < nbCol; ++i ) + { + gp_Pnt pPrev = fg.GetXYZ( i, 0 ); + for ( size_t j = 1; j < nbRow; ++j ) + { + gp_Pnt p = fg.GetXYZ( i, j ); + len[ j ] = len[ j-1 ] + p.Distance( pPrev ); + pPrev = p; + } + for ( size_t j = 0; j < nbRow; ++j ) + fg.GetIJK( i, j ).SetCoord( i2, len[ j ]/len.back() ); + } + + len.resize( nbCol ); + for ( size_t j = 0; j < nbRow; ++j ) + { + gp_Pnt pPrev = fg.GetXYZ( 0, j ); + for ( size_t i = 1; i < nbCol; ++i ) + { + gp_Pnt p = fg.GetXYZ( i, j ); + len[ i ] = len[ i-1 ] + p.Distance( pPrev ); + pPrev = p; + } + for ( size_t i = 0; i < nbCol; ++i ) + fg.GetIJK( i, j ).SetCoord( i1, len[ i ]/len.back() ); + } + } } //============================================================================= @@ -293,7 +354,6 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, { // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside //Unexpect aCatch(SalomeException); - MESSAGE("StdMeshers_Hexa_3D::Compute"); SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); // Shape verification @@ -306,11 +366,11 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, if ( exp.Next(), exp.More() ) return error(COMPERR_BAD_SHAPE, "More than one SHELL in the geometry"); - TopTools_IndexedMapOfShape FF; + TopTools_IndexedMapOfShape FF, EE; TopExp::MapShapes( aShape, TopAbs_FACE, FF); if ( FF.Extent() != 6) { - static StdMeshers_CompositeHexa_3D compositeHexa(_gen->GetANewId(), 0, _gen); + static StdMeshers_CompositeHexa_3D compositeHexa(_gen->GetANewId(), _gen); if ( !compositeHexa.Compute( aMesh, aShape )) return error( compositeHexa.GetComputeError() ); return true; @@ -318,13 +378,23 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, // Find sides of a cube // --------------------- - + + // tool creating quadratic elements if needed + SMESH_MesherHelper helper (aMesh); + _quadraticMesh = helper.IsQuadraticSubMesh(aShape); + + TopExp::MapShapes( aShape, TopAbs_EDGE, EE ); + SMESH_MesherHelper* faceHelper = ( EE.Size() == 12 ) ? 0 : &helper; + FaceQuadStructPtr quad[ 6 ]; - StdMeshers_Quadrangle_2D quadAlgo( _gen->GetANewId(), GetStudyId(), _gen); for ( int i = 0; i < 6; ++i ) { - if ( !( quad[i] = FaceQuadStructPtr( quadAlgo.CheckNbEdges( aMesh, FF( i+1 ))))) - return error( quadAlgo.GetComputeError() ); + if ( faceHelper ) + faceHelper->SetSubShape( FF( i+1 )); + if ( !( quad[i] = FaceQuadStructPtr( _quadAlgo->CheckNbEdges( aMesh, FF( i+1 ), + /*considerMesh=*/true, + faceHelper)))) + return error( _quadAlgo->GetComputeError() ); if ( quad[i]->side.size() != 4 ) return error( COMPERR_BAD_SHAPE, "Not a quadrangular box side" ); } @@ -379,10 +449,6 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, // Check presence of regular grid mesh on FACEs of the cube // ------------------------------------------------------------ - // tool creating quadratic elements if needed - SMESH_MesherHelper helper (aMesh); - _quadraticMesh = helper.IsQuadraticSubMesh(aShape); - for ( int i = 0; i < 6; ++i ) { const TopoDS_Face& F = aCubeSide[i]._quad->face; @@ -444,19 +510,19 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, { aCubeSide[i]._columns.resize( aCubeSide[i]._u2nodesMap.size() ); - int iFwd = 0, iRev = aCubeSide[i]._columns.size()-1; - int* pi = isReverse[i] ? &iRev : &iFwd; + size_t iFwd = 0, iRev = aCubeSide[i]._columns.size()-1; + size_t* pi = isReverse[i] ? &iRev : &iFwd; TParam2ColumnMap::iterator u2nn = aCubeSide[i]._u2nodesMap.begin(); for ( ; iFwd < aCubeSide[i]._columns.size(); --iRev, ++iFwd, ++u2nn ) aCubeSide[i]._columns[ *pi ].swap( u2nn->second ); aCubeSide[i]._u2nodesMap.clear(); } - + if ( proxymesh ) for ( int i = 0; i < 6; ++i ) - for ( unsigned j = 0; j < aCubeSide[i]._columns.size(); ++j) - for ( unsigned k = 0; k < aCubeSide[i]._columns[j].size(); ++k) + for ( size_t j = 0; j < aCubeSide[i]._columns.size(); ++j) + for ( size_t k = 0; k < aCubeSide[i]._columns[j].size(); ++k) { const SMDS_MeshNode* & n = aCubeSide[i]._columns[j][k]; n = proxymesh->GetProxyNode( n ); @@ -477,9 +543,23 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, _FaceGrid* fBack = & aCubeSide[ B_BACK ]; // cube size measured in nb of nodes - int x, xSize = fBottom->_columns.size() , X = xSize - 1; - int y, ySize = fLeft->_columns.size() , Y = ySize - 1; - int z, zSize = fLeft->_columns[0].size(), Z = zSize - 1; + size_t x, xSize = fBottom->_columns.size() , X = xSize - 1; + size_t y, ySize = fLeft->_columns.size() , Y = ySize - 1; + size_t z, zSize = fLeft->_columns[0].size(), Z = zSize - 1; + + // check sharing of FACEs (IPAL54417) + if ( fFront ->_columns.size() != xSize || + fBack ->_columns.size() != xSize || + fTop ->_columns.size() != xSize || + + fRight ->_columns.size() != ySize || + fTop ->_columns[0].size() != ySize || + fBottom->_columns[0].size() != ySize || + + fRight ->_columns[0].size() != zSize || + fFront ->_columns[0].size() != zSize || + fBack ->_columns[0].size() != zSize ) + return error( COMPERR_BAD_SHAPE, "Not sewed faces" ); // columns of internal nodes "rising" from nodes of fBottom _Indexer colIndex( xSize, ySize ); @@ -517,6 +597,14 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, } } + // compute normalized parameters of nodes on sides (PAL23189) + computeIJK( *fBottom, COO_X, COO_Y, /*z=*/0. ); + computeIJK( *fRight, COO_Y, COO_Z, /*x=*/1. ); + computeIJK( *fTop, COO_X, COO_Y, /*z=*/1. ); + computeIJK( *fLeft, COO_Y, COO_Z, /*x=*/0. ); + computeIJK( *fFront, COO_X, COO_Z, /*y=*/0. ); + computeIJK( *fBack, COO_X, COO_Z, /*y=*/1. ); + // projection points of the internal node on cube sub-shapes by which // coordinates of the internal node are computed vector pointsOnShapes( SMESH_Block::ID_Shell ); @@ -531,13 +619,14 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, pointsOnShapes[ SMESH_Block::ID_V011 ] = fTop->GetXYZ( 0, Y ); pointsOnShapes[ SMESH_Block::ID_V111 ] = fTop->GetXYZ( X, Y ); + gp_XYZ params; // normalized parameters of an internal node within the unit box for ( x = 1; x < xSize-1; ++x ) { - gp_XYZ params; // normalized parameters of internal node within a unit box - params.SetCoord( 1, x / double(X) ); + const double rX = x / double(X); for ( y = 1; y < ySize-1; ++y ) { - params.SetCoord( 2, y / double(Y) ); + const double rY = y / double(Y); + // a column to fill in during z loop vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )]; // projection points on horizontal edges @@ -554,23 +643,36 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = fTop ->GetXYZ( x, y ); for ( z = 1; z < zSize-1; ++z ) // z loop { - params.SetCoord( 3, z / double(Z) ); + const double rZ = z / double(Z); + + const gp_XYZ& pBo = fBottom->GetIJK( x, y ); + const gp_XYZ& pTo = fTop ->GetIJK( x, y ); + const gp_XYZ& pFr = fFront ->GetIJK( x, z ); + const gp_XYZ& pBa = fBack ->GetIJK( x, z ); + const gp_XYZ& pLe = fLeft ->GetIJK( y, z ); + const gp_XYZ& pRi = fRight ->GetIJK( y, z ); + params.SetCoord( 1, 0.5 * ( pBo.X() * ( 1. - rZ ) + pTo.X() * rZ + + pFr.X() * ( 1. - rY ) + pBa.X() * rY )); + params.SetCoord( 2, 0.5 * ( pBo.Y() * ( 1. - rZ ) + pTo.Y() * rZ + + pLe.Y() * ( 1. - rX ) + pRi.Y() * rX )); + params.SetCoord( 3, 0.5 * ( pFr.Z() * ( 1. - rY ) + pBa.Z() * rY + + pLe.Z() * ( 1. - rX ) + pRi.Z() * rX )); + // projection points on vertical edges - pointsOnShapes[ SMESH_Block::ID_E00z ] = fFront->GetXYZ( 0, z ); - pointsOnShapes[ SMESH_Block::ID_E10z ] = fFront->GetXYZ( X, z ); - pointsOnShapes[ SMESH_Block::ID_E01z ] = fBack->GetXYZ( 0, z ); + pointsOnShapes[ SMESH_Block::ID_E00z ] = fFront->GetXYZ( 0, z ); + pointsOnShapes[ SMESH_Block::ID_E10z ] = fFront->GetXYZ( X, z ); + pointsOnShapes[ SMESH_Block::ID_E01z ] = fBack->GetXYZ( 0, z ); pointsOnShapes[ SMESH_Block::ID_E11z ] = fBack->GetXYZ( X, z ); // projection points on vertical faces - pointsOnShapes[ SMESH_Block::ID_Fx0z ] = fFront->GetXYZ( x, z ); - pointsOnShapes[ SMESH_Block::ID_Fx1z ] = fBack ->GetXYZ( x, z ); - pointsOnShapes[ SMESH_Block::ID_F0yz ] = fLeft ->GetXYZ( y, z ); + pointsOnShapes[ SMESH_Block::ID_Fx0z ] = fFront->GetXYZ( x, z ); + pointsOnShapes[ SMESH_Block::ID_Fx1z ] = fBack ->GetXYZ( x, z ); + pointsOnShapes[ SMESH_Block::ID_F0yz ] = fLeft ->GetXYZ( y, z ); pointsOnShapes[ SMESH_Block::ID_F1yz ] = fRight->GetXYZ( y, z ); // compute internal node coordinates gp_XYZ coords; SMESH_Block::ShellPoint( params, pointsOnShapes, coords ); column[ z ] = helper.AddNode( coords.X(), coords.Y(), coords.Z() ); - } } } @@ -619,7 +721,7 @@ bool StdMeshers_Hexa_3D::Evaluate(SMESH_Mesh & aMesh, } if (meshFaces.size() != 6) { //return error(COMPERR_BAD_SHAPE, TComm(meshFaces.size())<<" instead of 6 faces in a block"); - static StdMeshers_CompositeHexa_3D compositeHexa(-10, 0, aMesh.GetGen()); + static StdMeshers_CompositeHexa_3D compositeHexa(-10, aMesh.GetGen()); return compositeHexa.Evaluate(aMesh, aShape, aResMap); } @@ -729,7 +831,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper static StdMeshers_HexaFromSkin_3D * algo = 0; if ( !algo ) { SMESH_Gen* gen = aMesh.GetGen(); - algo = new StdMeshers_HexaFromSkin_3D( gen->GetANewId(), 0, gen ); + algo = new StdMeshers_HexaFromSkin_3D( gen->GetANewId(), gen ); } algo->InitComputeError(); algo->Compute( aMesh, aHelper ); @@ -797,7 +899,7 @@ SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh & aMesh, static StdMeshers_Prism_3D * aPrism3D = 0; if ( !aPrism3D ) { SMESH_Gen* gen = aMesh.GetGen(); - aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), 0, gen ); + aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), gen ); } SMESH_Hypothesis::Hypothesis_Status aStatus; if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) { @@ -828,7 +930,7 @@ bool EvaluatePentahedralMesh(SMESH_Mesh & aMesh, static StdMeshers_Prism_3D * aPrism3D = 0; if ( !aPrism3D ) { SMESH_Gen* gen = aMesh.GetGen(); - aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), 0, gen ); + aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), gen ); } SMESH_Hypothesis::Hypothesis_Status aStatus; if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {