X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_CompositeHexa_3D.cxx;h=d43ceadafcb83267eb9d65d47cac66703c7b34b8;hp=b0f10f8d0d174ec23abf4ef64c1425175b8cb0c5;hb=refs%2Ftags%2FV9_7_0b1;hpb=ecea056f6165eddf99403a8648657a9a573bcbe0 diff --git a/src/StdMeshers/StdMeshers_CompositeHexa_3D.cxx b/src/StdMeshers/StdMeshers_CompositeHexa_3D.cxx index b0f10f8d0..d43ceadaf 100644 --- a/src/StdMeshers/StdMeshers_CompositeHexa_3D.cxx +++ b/src/StdMeshers/StdMeshers_CompositeHexa_3D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2019 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2021 CEA/DEN, EDF R&D, OPEN CASCADE // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -38,10 +38,13 @@ #include "SMESH_MeshAlgos.hxx" #include "SMESH_MesherHelper.hxx" #include "SMESH_subMesh.hxx" +#include "StdMeshers_BlockRenumber.hxx" +#include "StdMeshers_FaceSide.hxx" #include "StdMeshers_ViscousLayers.hxx" #include #include +#include #include #include #include @@ -59,6 +62,7 @@ #include #include +#include #include #include @@ -120,7 +124,7 @@ public: bool Contain( const TopoDS_Vertex& vertex ) const; void AppendSide( const _FaceSide& side ); void SetBottomSide( int i ); - int GetNbSegments(SMESH_ProxyMesh& mesh, const SMESHDS_SubMesh* smToCheckEdges=0) const; + smIdType GetNbSegments(SMESH_ProxyMesh& mesh, const SMESHDS_SubMesh* smToCheckEdges=0) const; bool StoreNodes(SMESH_ProxyMesh& mesh, vector& myGrid, bool reverse, bool isProxy, const SMESHDS_SubMesh* smToCheckEdges=0 ); void SetID(EQuadSides id) { myID = id; } @@ -180,6 +184,8 @@ public: //** Methods to find and orient faces of 6 sides of the box **// TChildIterator GetChildren() const { return TChildIterator( myChildren.begin(), myChildren.end()); } + bool Contain( const TopoDS_Vertex& vertex ) const { return mySides.Contain( vertex ); } + public: //** Loading and access to mesh **// //!< Load nodes of a mesh @@ -271,7 +277,7 @@ private: //================================================================================ StdMeshers_CompositeHexa_3D::StdMeshers_CompositeHexa_3D(int hypId, SMESH_Gen* gen) - :SMESH_3D_Algo(hypId, gen) + :SMESH_3D_Algo(hypId, gen), _blockRenumberHyp( nullptr ) { _name = "CompositeHexa_3D"; _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID); // 1 bit /shape type @@ -283,10 +289,11 @@ StdMeshers_CompositeHexa_3D::StdMeshers_CompositeHexa_3D(int hypId, SMESH_Gen* g */ //================================================================================ -bool StdMeshers_CompositeHexa_3D::CheckHypothesis(SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape, +bool StdMeshers_CompositeHexa_3D::CheckHypothesis(SMESH_Mesh& /*aMesh*/, + const TopoDS_Shape& /*aShape*/, Hypothesis_Status& aStatus) { + _blockRenumberHyp = nullptr; aStatus = HYP_OK; return true; } @@ -420,9 +427,9 @@ namespace const TopTools_MapOfShape& cornerVV, TopTools_MapOfShape& internEE) { - TopTools_IndexedMapOfShape subEE; + TopTools_IndexedMapOfShape subEE, subFF; TopExp::MapShapes( shape, TopAbs_EDGE, subEE ); - //TopExp::MapShapes( shape, TopAbs_FACE, subFF ); + TopExp::MapShapes( shape, TopAbs_FACE, subFF ); TopoDS_Vertex VV[2]; TopTools_MapOfShape subChecked, ridgeEE; @@ -460,6 +467,8 @@ namespace { if ( !SMESH_MesherHelper::IsSubShape( ridgeE, *F )) continue; + if ( !subFF.Contains( *F )) + continue; if ( isContinuousMesh( ridgeE, TopoDS::Edge( *E ), TopoDS::Face( *F ), mesh )) { nextRidgeE = *E; @@ -543,6 +552,163 @@ namespace return faceFound; } + //================================================================================ + /*! + * \brief Rearrange block sides according to StdMeshers_BlockRenumber hypothesis + */ + //================================================================================ + + bool arrangeForRenumber( list< _QuadFaceGrid >& blockSides, + const TopTools_MapOfShape& cornerVertices, + SMESH_Mesh* mesh, + TopoDS_Vertex& v000, + TopoDS_Vertex& v001 ) + { + if ( v000.IsNull() ) + { + // block CS is not defined; + // renumber only if the block has an edge parallel to an axis of global CS + + v000 = StdMeshers_RenumberHelper::GetVertex000( cornerVertices ); + } + + Bnd_B3d bbox; + for ( auto it = cornerVertices.cbegin(); it != cornerVertices.cend(); ++it ) + bbox.Add( BRep_Tool::Pnt( TopoDS::Vertex( *it ))); + double tol = 1e-5 * Sqrt( bbox.SquareExtent() ); + + // get block edges starting at v000 + + std::vector< const _FaceSide* > edgesAtV000; + std::vector< gp_Vec > edgeDir; + std::vector< int > iParallel; // 0 - none, 1 - X, 2 - Y, 3 - Z + TopTools_MapOfShape lastVertices; + for ( _QuadFaceGrid & quad: blockSides ) + { + for ( int iS = 0; iS < 4 && edgesAtV000.size() < 3; ++iS ) + { + const _FaceSide* side = & quad.GetSide( iS ); + TopoDS_Vertex v1 = side->FirstVertex(), v2 = side->LastVertex(); + if (( v1.IsSame( v000 ) && !lastVertices.Contains( v2 )) || + ( v2.IsSame( v000 ) && !lastVertices.Contains( v1 ))) + { + bool reverse = v2.IsSame( v000 ); + if ( reverse ) + std::swap( v1, v2 ); + lastVertices.Add( v2 ); + + edgesAtV000.push_back( side ); + + gp_Pnt pf = BRep_Tool::Pnt( v1 ); + gp_Pnt pl = BRep_Tool::Pnt( v2 ); + gp_Vec vec( pf, pl ); + edgeDir.push_back( vec ); + + iParallel.push_back( 0 ); + if ( !v001.IsNull() ) + { + if ( quad.IsComplex() ) + for ( _QuadFaceGrid::TChildIterator chIt = quad.GetChildren(); chIt.more(); ) + { + const _QuadFaceGrid& child = chIt.next(); + if ( child.GetSide( iS ).Contain( v001 )) + { + iParallel.back() = 3; + break; + } + } + else if ( side->Contain( v001 )) + iParallel.back() = 3; + } + else + { + bool isStraight = true; + std::list< TopoDS_Edge > edges; + for ( int iE = 0; true; ++iE ) + { + TopoDS_Edge edge = side->Edge( iE ); + if ( edge.IsNull() ) + break; + edges.push_back( edge ); + if ( isStraight ) + isStraight = SMESH_Algo::IsStraight( edge ); + } + // is parallel to a GCS axis? + if ( isStraight ) + { + int nbDiff = (( Abs( vec.X() ) > tol ) + + ( Abs( vec.Y() ) > tol ) + + ( Abs( vec.Z() ) > tol ) ); + if ( nbDiff == 1 ) + iParallel.back() = ( Abs( vec.X() ) > tol ) ? 1 : ( Abs( vec.Y() ) > tol ) ? 2 : 3; + } + else + { + TopoDS_Face nullFace; + StdMeshers_FaceSide fSide( nullFace, edges, mesh, true, true ); + edgeDir.back() = gp_Vec( pf, fSide.Value3d( reverse ? 0.99 : 0.01 )); + } + } + } + } + } + if ( std::accumulate( iParallel.begin(), iParallel.end(), 0 ) == 0 ) + return false; + + // find edge OZ and edge OX + const _FaceSide* edgeOZ = nullptr, *edgeOY = nullptr, *edgeOX = nullptr; + auto iZIt = std::find( iParallel.begin(), iParallel.end(), 3 ); + if ( iZIt != iParallel.end() ) + { + int i = std::distance( iParallel.begin(), iZIt ); + edgeOZ = edgesAtV000[ i ]; + int iE1 = SMESH_MesherHelper::WrapIndex( i + 1, edgesAtV000.size() ); + int iE2 = SMESH_MesherHelper::WrapIndex( i + 2, edgesAtV000.size() ); + if (( edgeDir[ iE1 ] ^ edgeDir[ iE2 ] ) * edgeDir[ i ] < 0 ) + std::swap( iE1, iE2 ); + edgeOX = edgesAtV000[ iE1 ]; + edgeOY = edgesAtV000[ iE2 ]; + } + else + { + for ( size_t i = 0; i < edgesAtV000.size(); ++i ) + { + if ( !iParallel[ i ] ) + continue; + int iE1 = SMESH_MesherHelper::WrapIndex( i + 1, edgesAtV000.size() ); + int iE2 = SMESH_MesherHelper::WrapIndex( i + 2, edgesAtV000.size() ); + if (( edgeDir[ iE1 ] ^ edgeDir[ iE2 ] ) * edgeDir[ i ] < 0 ) + std::swap( iE1, iE2 ); + edgeOZ = edgesAtV000[ iParallel[i] == 1 ? iE2 : iE1 ]; + edgeOX = edgesAtV000[ iParallel[i] == 1 ? i : iE1 ]; + edgeOY = edgesAtV000[ iParallel[i] == 1 ? iE1 : i ]; + break; + } + } + + if ( !edgeOZ || !edgeOX || !edgeOY ) + return false; + + TopoDS_Vertex v100 = edgeOX->LastVertex(); + if ( v100.IsSame( v000 )) + v100 = edgeOX->FirstVertex(); + + // Find the left quad, one including v000 but not v100 + + for ( auto quad = blockSides.begin(); quad != blockSides.end(); ++quad ) + { + if ( quad->Contain( v000 ) && !quad->Contain( v100 )) // it's a left quad + { + if ( quad != blockSides.begin() ) + blockSides.splice( blockSides.begin(), blockSides, quad ); + blockSides.front().SetBottomSide( *edgeOZ ); // edgeOY + + return true; + } + } + return false; + } + } // namespace //================================================================================ @@ -555,6 +721,7 @@ bool StdMeshers_CompositeHexa_3D::findBoxFaces( const TopoDS_Shape& shape, list< _QuadFaceGrid >& boxFaces, SMESH_Mesh& mesh, SMESH_ProxyMesh& proxyMesh, + bool& toRenumber, _QuadFaceGrid * & fBottom, _QuadFaceGrid * & fTop, _QuadFaceGrid * & fFront, @@ -607,6 +774,22 @@ bool StdMeshers_CompositeHexa_3D::findBoxFaces( const TopoDS_Shape& shape, for ( exp.Init( shape, TopAbs_FACE); exp.More(); exp.Next(), ++boxFace ) boxFace->Init( TopoDS::Face( exp.Current() ), proxyMesh ); } + + toRenumber = _blockRenumberHyp; + if ( toRenumber ) + { + TopoDS_Vertex v000, v001; + _blockRenumberHyp->IsSolidIncluded( mesh, shape, v000, v001 ); + + toRenumber = arrangeForRenumber( boxFaces, cornerVertices, &mesh, v000, v001 ); + + if ( toRenumber ) + { + mesh.GetMeshDS()->Modified(); + mesh.GetMeshDS()->CompactMesh(); // remove numbering holes + } + } + // ---------------------------------------- // Find out position of faces within a box // ---------------------------------------- @@ -684,8 +867,9 @@ bool StdMeshers_CompositeHexa_3D::Compute(SMESH_Mesh& theMesh, // Try to find 6 side faces // ------------------------- list< _QuadFaceGrid > boxFaceContainer; + bool toRenumber = false; _QuadFaceGrid *fBottom, *fTop, *fFront, *fBack, *fLeft, *fRight; - if ( ! findBoxFaces( theShape, boxFaceContainer, theMesh, *proxyMesh, + if ( ! findBoxFaces( theShape, boxFaceContainer, theMesh, *proxyMesh, toRenumber, fBottom, fTop, fFront, fBack, fLeft, fRight)) return false; @@ -709,6 +893,8 @@ bool StdMeshers_CompositeHexa_3D::Compute(SMESH_Mesh& theMesh, fRight ->ComputeIJK( COO_Y, COO_Z, /*x=*/1. ); fTop ->ComputeIJK( COO_X, COO_Y, /*z=*/1. ); + StdMeshers_RenumberHelper renumHelper( theMesh, _blockRenumberHyp ); + int x, xSize = fBottom->GetNbHoriSegments(*proxyMesh) + 1, X = xSize - 1; int y, ySize = fBottom->GetNbVertSegments(*proxyMesh) + 1, Y = ySize - 1; int z, zSize = fFront ->GetNbVertSegments(*proxyMesh) + 1, Z = zSize - 1; @@ -766,8 +952,23 @@ bool StdMeshers_CompositeHexa_3D::Compute(SMESH_Mesh& theMesh, gp_XYZ params; // normalized parameters of an internal node within the unit box - for ( x = 1; x < xSize-1; ++x ) + if ( toRenumber ) + for ( y = 0; y < ySize; ++y ) + { + vector< const SMDS_MeshNode* >& columnXy = columns[ colIndex( X, y )]; + for ( z = 0; z < zSize; ++z ) + renumHelper.AddReplacingNode( columnXy[ z ] ); + } + + for ( x = X-1; x > 0; --x ) { + if ( toRenumber ) + { + vector< const SMDS_MeshNode* >& columnX0 = columns[ colIndex( x, 0 )]; + for ( z = 0; z < zSize; ++z ) + renumHelper.AddReplacingNode( columnX0[ z ] ); + } + const double rX = x / double(X); for ( y = 1; y < ySize-1; ++y ) { @@ -786,6 +987,10 @@ bool StdMeshers_CompositeHexa_3D::Compute(SMESH_Mesh& theMesh, // points projections on horizontal faces pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = fBottom->GetXYZ( x, y ); pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = fTop ->GetXYZ( x, y ); + + if ( toRenumber ) + renumHelper.AddReplacingNode( column[ 0 ] ); + for ( z = 1; z < zSize-1; ++z ) // z loop { // compute normalized parameters of an internal node within the unit box @@ -830,16 +1035,35 @@ bool StdMeshers_CompositeHexa_3D::Compute(SMESH_Mesh& theMesh, //cout << "Params: ( "<< params.X()<<", "<& columnXY = columns[ colIndex( x, Y )]; + for ( z = 0; z < zSize; ++z ) + renumHelper.AddReplacingNode( columnXY[ z ] ); } - } + } // for ( x = X-1; x > 0; --x ) + + if ( toRenumber ) + for ( y = 0; y < ySize; ++y ) + { + vector< const SMDS_MeshNode* >& column0Y = columns[ colIndex( 0, y )]; + for ( z = 0; z < zSize; ++z ) + renumHelper.AddReplacingNode( column0Y[ z ] ); + } + + // faces no more needed, free memory boxFaceContainer.clear(); // ---------------- // Add hexahedrons // ---------------- - for ( x = 0; x < xSize-1; ++x ) { + for ( x = xSize-2; true; --x ) { for ( y = 0; y < ySize-1; ++y ) { vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )]; vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )]; @@ -848,11 +1072,22 @@ bool StdMeshers_CompositeHexa_3D::Compute(SMESH_Mesh& theMesh, for ( z = 0; z < zSize-1; ++z ) { // bottom face normal of a hexa mush point outside the volume - helper.AddVolume(col00[z], col01[z], col11[z], col10[z], - col00[z+1], col01[z+1], col11[z+1], col10[z+1]); + helper.AddVolume(col10[z], col11[z], col11[z+1], col10[z+1], + col00[z], col01[z], col01[z+1], col00[z+1]); } } + if ( x == 0) + break; + } + if ( toRenumber ) + renumHelper.DoReplaceNodes(); + + if ( _blockRenumberHyp ) + { + return error( _blockRenumberHyp->CheckHypothesis( theMesh, theShape )); + } + return true; } @@ -873,7 +1108,8 @@ bool StdMeshers_CompositeHexa_3D::Evaluate(SMESH_Mesh& theMesh, // ------------------------- list< _QuadFaceGrid > boxFaceContainer; _QuadFaceGrid *fBottom, *fTop, *fFront, *fBack, *fLeft, *fRight; - if ( ! findBoxFaces( theShape, boxFaceContainer, theMesh, *proxyMesh, + bool toRenumber = false; + if ( ! findBoxFaces( theShape, boxFaceContainer, theMesh, *proxyMesh, toRenumber, fBottom, fTop, fFront, fBack, fLeft, fRight)) return false; @@ -885,32 +1121,32 @@ bool StdMeshers_CompositeHexa_3D::Evaluate(SMESH_Mesh& theMesh, lessComplexSide = & *face; // Get an 1D size of lessComplexSide - int nbSeg1 = 0; + smIdType nbSeg1 = 0; vector edges; if ( !lessComplexSide->GetHoriEdges(edges) ) return false; for ( size_t i = 0; i < edges.size(); ++i ) { - const vector& nbElems = aResMap[ theMesh.GetSubMesh( edges[i] )]; + const vector& nbElems = aResMap[ theMesh.GetSubMesh( edges[i] )]; if ( !nbElems.empty() ) - nbSeg1 += Max( nbElems[ SMDSEntity_Edge ], nbElems[ SMDSEntity_Quad_Edge ]); + nbSeg1 += std::max( nbElems[ SMDSEntity_Edge ], nbElems[ SMDSEntity_Quad_Edge ]); } // Get an 1D size of a box side orthogonal to lessComplexSide - int nbSeg2 = 0; + smIdType nbSeg2 = 0; _QuadFaceGrid* ortoSide = lessComplexSide->FindAdjacentForSide( Q_LEFT, boxFaceContainer, B_UNDEFINED ); edges.clear(); if ( !ortoSide || !ortoSide->GetHoriEdges(edges) ) return false; for ( size_t i = 0; i < edges.size(); ++i ) { - const vector& nbElems = aResMap[ theMesh.GetSubMesh( edges[i] )]; + const vector& nbElems = aResMap[ theMesh.GetSubMesh( edges[i] )]; if ( !nbElems.empty() ) - nbSeg2 += Max( nbElems[ SMDSEntity_Edge ], nbElems[ SMDSEntity_Quad_Edge ]); + nbSeg2 += std::max( nbElems[ SMDSEntity_Edge ], nbElems[ SMDSEntity_Quad_Edge ]); } // Get an 2D size of a box side orthogonal to lessComplexSide - int nbFaces = 0, nbQuadFace = 0; + smIdType nbFaces = 0, nbQuadFace = 0; list< TopoDS_Face > sideFaces; if ( ortoSide->IsComplex() ) for ( _QuadFaceGrid::TChildIterator child = ortoSide->GetChildren(); child.more(); ) @@ -921,7 +1157,7 @@ bool StdMeshers_CompositeHexa_3D::Evaluate(SMESH_Mesh& theMesh, list< TopoDS_Face >::iterator f = sideFaces.begin(); for ( ; f != sideFaces.end(); ++f ) { - const vector& nbElems = aResMap[ theMesh.GetSubMesh( *f )]; + const vector& nbElems = aResMap[ theMesh.GetSubMesh( *f )]; if ( !nbElems.empty() ) { nbFaces = nbElems[ SMDSEntity_Quadrangle ]; @@ -930,8 +1166,8 @@ bool StdMeshers_CompositeHexa_3D::Evaluate(SMESH_Mesh& theMesh, } // Fill nb of elements - vector aResVec(SMDSEntity_Last,0); - int nbSeg3 = ( nbFaces + nbQuadFace ) / nbSeg2; + vector aResVec(SMDSEntity_Last,0); + smIdType nbSeg3 = ( nbFaces + nbQuadFace ) / nbSeg2; aResVec[SMDSEntity_Node] = (nbSeg1-1) * (nbSeg2-1) * (nbSeg3-1); aResVec[SMDSEntity_Hexa] = nbSeg1 * nbFaces; aResVec[SMDSEntity_Quad_Hexa] = nbSeg1 * nbQuadFace; @@ -1244,7 +1480,7 @@ bool _QuadFaceGrid::LoadGrid( SMESH_ProxyMesh& mesh ) if ( fIt->next()->NbNodes() % 4 > 0 ) return error("Non-quadrangular mesh faces are not allowed on sides of a composite block"); - bool isProxy, isTmpElem; + bool isProxy = false, isTmpElem = false; if ( faceSubMesh && faceSubMesh->NbElements() > 0 ) { isProxy = dynamic_cast< const SMESH_ProxyMesh::SubMesh* >( faceSubMesh ); @@ -1660,7 +1896,7 @@ bool _QuadFaceGrid::GetNormal( const TopoDS_Vertex& v, gp_Vec& n ) const n = d1u.Crossed( d1v ); return true; } - catch (Standard_Failure) { + catch (Standard_Failure&) { return false; } } @@ -1993,9 +2229,9 @@ void _FaceSide::SetBottomSide( int i ) //purpose : //======================================================================= -int _FaceSide::GetNbSegments(SMESH_ProxyMesh& mesh, const SMESHDS_SubMesh* smToCheckEdges) const +smIdType _FaceSide::GetNbSegments(SMESH_ProxyMesh& mesh, const SMESHDS_SubMesh* smToCheckEdges) const { - int nb = 0; + smIdType nb = 0; if ( myChildren.empty() ) { nb = mesh.GetSubMesh(myEdge)->NbElements(); @@ -2092,7 +2328,7 @@ bool _FaceSide::StoreNodes(SMESH_ProxyMesh& mesh, smToCheckEdges, mesh )) break; else - nodes.erase( --( u_node1.base() )); + nodes.erase( --(( u_node2 = u_node1 ).base() )); } }