X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MesherHelper.cxx;h=242900a749d0877a3c4f9e25ace56a0c2fdf43eb;hp=ce4e90cc4bb21feeb8aaf93b691f02b79ade5444;hb=3a821b48ba986795abde72f524fa1a62a5802df2;hpb=e74c29b7867ca230102fc580d831dda6b81b3b4c diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index ce4e90cc4..242900a74 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -139,9 +139,9 @@ SMESHDS_Mesh* SMESH_MesherHelper::GetMeshDS() const //======================================================================= //function : IsQuadraticSubMesh -//purpose : Check submesh for given shape: if all elements on this shape +//purpose : Check sub-meshes of a given shape: if all elements on sub-shapes // are quadratic, quadratic elements will be created. -// Also fill myTLinkNodeMap +// Fill myTLinkNodeMap //======================================================================= bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) @@ -149,7 +149,6 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) SMESHDS_Mesh* meshDS = GetMeshDS(); // we can create quadratic elements only if all elements // created on sub-shapes of given shape are quadratic - // also we have to fill myTLinkNodeMap myCreateQuadratic = true; mySeamShapeIds.clear(); myDegenShapeIds.clear(); @@ -162,9 +161,6 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh) } SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge ); - - //int nbOldLinks = myTLinkNodeMap.size(); - if ( !myMesh->HasShapeToMesh() ) { if (( myCreateQuadratic = myMesh->NbFaces( ORDER_QUADRATIC ))) @@ -373,6 +369,27 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) } } +//======================================================================= +/*! + * \brief Copy shape information from another helper. Used to improve performance + * since SetSubShape() can be time consuming if there are many edges + */ +//======================================================================= + +void SMESH_MesherHelper::CopySubShapeInfo(const SMESH_MesherHelper& other) +{ + this->myShape = other.myShape; + this->myShapeID = other.myShapeID; + this->myDegenShapeIds = other.myDegenShapeIds; + this->mySeamShapeIds = other.mySeamShapeIds; + this->myPar1[0] = other.myPar1[0]; + this->myPar1[1] = other.myPar1[1]; + this->myPar2[0] = other.myPar2[0]; + this->myPar2[1] = other.myPar2[1]; + this->myParIndex = other.myParIndex; + this->myFace2Surface = other.myFace2Surface; +} + //======================================================================= //function : ShapeToIndex //purpose : Convert a shape to its index in the SMESHDS_Mesh @@ -807,7 +824,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face& F, MESSAGE( "SMESH_MesherHelper::CheckNodeUV() failed to project" ); return false; } - Quantity_Parameter U,V; + Standard_Real U,V; projector.LowerDistanceParameters(U,V); uv.SetCoord( U,V ); surfPnt = surface->Value( U, V ); @@ -1113,7 +1130,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge& E, MESSAGE( "SMESH_MesherHelper::CheckNodeU() failed to project" ); return false; } - Quantity_Parameter U = projector->LowerDistanceParameter(); + Standard_Real U = projector->LowerDistanceParameter(); u = double( U ); curvPnt = curve->Value( u ); dist = nodePnt.Distance( curvPnt ); @@ -2162,13 +2179,30 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1, const SMDS_MeshNode* n14 = GetMediumNode( n1, n4, force3d, TopAbs_SOLID ); const SMDS_MeshNode* n25 = GetMediumNode( n2, n5, force3d, TopAbs_SOLID ); const SMDS_MeshNode* n36 = GetMediumNode( n3, n6, force3d, TopAbs_SOLID ); + if ( myCreateBiQuadratic ) + { + const SMDS_MeshNode* n1245 = GetCentralNode( n1,n2,n4,n5,n12,n25,n45,n14,force3d ); + const SMDS_MeshNode* n1346 = GetCentralNode( n1,n3,n4,n6,n31,n36,n64,n14,force3d ); + const SMDS_MeshNode* n2356 = GetCentralNode( n2,n3,n6,n5,n23,n36,n56,n25,force3d ); - if(id) - elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, - n12, n23, n31, n45, n56, n64, n14, n25, n36, id); + if(id) + elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, + n12, n23, n31, n45, n56, n64, n14, n25, n36, + n1245, n2356, n1346, id); + else + elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, + n12, n23, n31, n45, n56, n64, n14, n25, n36, + n1245, n2356, n1346); + } else - elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, - n12, n23, n31, n45, n56, n64, n14, n25, n36); + { + if(id) + elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6, + n12, n23, n31, n45, n56, n64, n14, n25, n36, id); + else + elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6, + n12, n23, n31, n45, n56, n64, n14, n25, n36); + } } if ( mySetElemOnShape && myShapeID > 0 ) meshDS->SetMeshElementOnShape( elem, myShapeID ); @@ -2652,7 +2686,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 } // nb rows of nodes - size_t prevNbRows = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here + size_t prevNbRows = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here size_t expectNbRows = faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 ); // to be added // fill theParam2ColumnMap column by column by passing from nodes on @@ -2825,8 +2859,9 @@ bool SMESH_MesherHelper::IsStructured( SMESH_subMesh* faceSM ) //purpose : Return true if 2D mesh on FACE is ditorted //======================================================================= -bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, - bool checkUV) +bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, + bool checkUV, + SMESH_MesherHelper* faceHelper) { if ( !faceSM || faceSM->GetSubShape().ShapeType() != TopAbs_FACE ) return false; @@ -2834,12 +2869,26 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, bool haveBadFaces = false; SMESH_MesherHelper helper( *faceSM->GetFather() ); + if ( faceHelper ) + helper.CopySubShapeInfo( *faceHelper ); helper.SetSubShape( faceSM->GetSubShape() ); const TopoDS_Face& F = TopoDS::Face( faceSM->GetSubShape() ); SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( F ); if ( !smDS || smDS->NbElements() == 0 ) return false; + bool subIdsValid = true; // shape ID of nodes is OK + if ( helper.HasSeam() ) + { + // check if nodes are bound to seam edges + SMESH_subMeshIteratorPtr smIt = faceSM->getDependsOnIterator(/*includeSelf=*/false); + while ( smIt->more() && subIdsValid ) + { + SMESH_subMesh* sm = smIt->next(); + if ( helper.IsSeamShape( sm->GetId() ) && sm->IsEmpty() ) + subIdsValid = false; + } + } SMDS_ElemIteratorPtr faceIt = smDS->GetElements(); double prevArea = 0; vector< const SMDS_MeshNode* > nodes; @@ -2864,12 +2913,20 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, if ( isOnDegen ) continue; } - // prepare to getting UVs + // prepare for getting UVs const SMDS_MeshNode* inFaceNode = 0; if ( helper.HasSeam() ) { for ( size_t i = 0; ( i < nodes.size() && !inFaceNode ); ++i ) if ( !helper.IsSeamShape( nodes[ i ]->getshapeId() )) + { inFaceNode = nodes[ i ]; + if ( !subIdsValid ) + { + gp_XY uv = helper.GetNodeUV( F, inFaceNode ); + if ( helper.IsOnSeam( uv )) + inFaceNode = NULL; + } + } if ( !inFaceNode ) continue; } @@ -2878,6 +2935,14 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, for ( size_t i = 0; i < nodes.size(); ++i ) uv[ i ] = helper.GetNodeUV( F, nodes[ i ], inFaceNode, toCheckUV ); + if ( !subIdsValid ) // fix uv on seam + { + gp_XY uvInFace = helper.GetNodeUV( F, inFaceNode ); + for ( size_t i = 0; i < uv.size(); ++i ) + if ( helper.IsOnSeam( uv[i] )) + uv[i] = helper.getUVOnSeam( uv[i], uvInFace ).XY(); + } + // compare orientation of triangles double faceArea = 0; for ( int iT = 0, nbT = nodes.size()-2; iT < nbT; ++iT ) @@ -2898,7 +2963,7 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh* faceSM, * \brief Find out elements orientation on a geometrical face * \param theFace - The face correctly oriented in the shape being meshed * \retval bool - true if the face normal and the normal of first element - * in the correspoding submesh point in different directions + * in the corresponding submesh point in different directions */ //================================================================================ @@ -3087,7 +3152,7 @@ TopAbs_Orientation SMESH_MesherHelper::GetSubShapeOri(const TopoDS_Shape& shape, //======================================================================= //function : IsSubShape -//purpose : +//purpose : //======================================================================= bool SMESH_MesherHelper::IsSubShape( const TopoDS_Shape& shape, @@ -3395,6 +3460,43 @@ double SMESH_MesherHelper::GetOtherParam(const double param) const return fabs(param-myPar1[i]) < fabs(param-myPar2[i]) ? myPar2[i] : myPar1[i]; } +//======================================================================= +//function : NbRealSeam +//purpose : Return a number of real seam edges in the shape set through +// IsQuadraticSubMesh() or SetSubShape(). A real seam edge encounters twice in a wire +//======================================================================= + +size_t SMESH_MesherHelper::NbRealSeam() const +{ + size_t nb = 0; + + std::set< int >::const_iterator id = mySeamShapeIds.begin(); + for ( ; id != mySeamShapeIds.end(); ++id ) + if ( *id < 0 ) ++nb; + else break; + + return nb; +} + +//======================================================================= +//function : IsOnSeam +//purpose : Check if UV is on seam. Return 0 if not, 1 for U seam, 2 for V seam +//======================================================================= + +int SMESH_MesherHelper::IsOnSeam(const gp_XY& uv) const +{ + for ( int i = U_periodic; i <= V_periodic ; ++i ) + if ( myParIndex & i ) + { + double p = uv.Coord( i ); + double tol = ( myPar2[i-1] - myPar1[i-1] ) / 100.; + if ( Abs( p - myPar1[i-1] ) < tol || + Abs( p - myPar2[i-1] ) < tol ) + return i; + } + return 0; +} + namespace { //======================================================================= @@ -4041,9 +4143,9 @@ namespace { // Structures used by FixQuadraticElements() { // code is valid for convex faces only gp_XYZ gc(0,0,0); - for ( TIDSortedNodeSet::const_iterator n = begin(); n!=end(); ++n) - gc += XYZ( *n ) / size(); - for (unsigned i = 0; i < _sides.size(); ++i ) + for ( TIDSortedNodeSet::const_iterator n = begin(); n != end(); ++n ) + gc += XYZ( *n ) / double( size() ); + for ( size_t i = 0; i < _sides.size(); ++i ) { if ( _sides[i] == bentLink ) continue; gp_Vec linkNorm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2())); @@ -4059,12 +4161,12 @@ namespace { // Structures used by FixQuadraticElements() return true; } return false; - + } //================================================================================ /*! - * \brief Find pairs of continues faces + * \brief Find pairs of continues faces */ //================================================================================ @@ -4276,7 +4378,7 @@ namespace { // Structures used by FixQuadraticElements() vector< TChain> & resultChains, SMDS_TypeOfPosition pos ) { - // put links in the set and evalute number of result chains by number of boundary links + // put links in the set and evaluate number of result chains by number of boundary links TLinkSet linkSet; size_t nbBndLinks = 0; for ( TChain::iterator lnk = allLinks.begin(); lnk != allLinks.end(); ++lnk ) { @@ -4799,6 +4901,7 @@ namespace { // Structures used by FixQuadraticElements() void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, bool volumeOnly) { + //MESSAGE("FixQuadraticElements " << volumeOnly); // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion if ( getenv("NO_FixQuadraticElements") ) return; @@ -4840,6 +4943,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, nbfaces = faces.Extent(); /*avoid "unused varianbles": */ nbfaces++, nbfaces--; #endif for ( TopTools_MapIteratorOfMapOfShape fIt( faces ); fIt.More(); fIt.Next() ) { + MESSAGE("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key())); MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key())); SMESH_MesherHelper h(*myMesh); h.SetSubShape( fIt.Key() ); @@ -5040,7 +5144,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, MSG("Internal chain - ignore"); continue; } - // mesure chain length and compute link position along the chain + // measure chain length and compute link position along the chain double chainLen = 0; vector< double > linkPos; TChain savedChain; // backup @@ -5196,10 +5300,11 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, // 4. Move nodes // ------------- - TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa; + TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa, biQuadPenta; const bool toFixCentralNodes = ( myMesh->NbBiQuadQuadrangles() + myMesh->NbBiQuadTriangles() + - myMesh->NbTriQuadraticHexas() ); + myMesh->NbTriQuadraticHexas() + + myMesh->NbBiQuadPrisms()); double distXYZ[4]; faceHlp.ToFixNodeParameters( true ); @@ -5235,6 +5340,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, case SMDSEntity_BiQuad_Quadrangle: biQuadQuas.insert( e ); break; case SMDSEntity_BiQuad_Triangle: biQuadTris.insert( e ); break; case SMDSEntity_TriQuad_Hexa: triQuadHexa.insert( e ); break; + case SMDSEntity_BiQuad_Penta: biQuadPenta.insert( e ); break; default:; } } @@ -5385,6 +5491,16 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, nCenterCoords.X(), nCenterCoords.Y(), nCenterCoords.Z()); } } + // treat tri-quadratic hexahedra + { + SMDS_VolumeTool volExp; + TIDSortedElemSet::iterator pentIt = biQuadPenta.begin(); + for ( ; pentIt != biQuadPenta.end(); ++pentIt ) + { + MESSAGE("---"); + volExp.Set( *pentIt, /*ignoreCentralNodes=*/false ); + } + } #ifdef _DEBUG_ // avoid warning: defined but not used operator<<() SMESH_Comment() << *links.begin() << *faces.begin();