X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FControls%2FSMESH_Controls.cxx;h=1ef3a5f2f4018efa8b7c14ca2454928338ecaa43;hp=fc7ba629b818dc78abcd10be7cda752a37276553;hb=24825596b3595ab0c0b1988368a886864d1080ef;hpb=6472eab132825fec572beda8276947593f85ffa1 diff --git a/src/Controls/SMESH_Controls.cxx b/src/Controls/SMESH_Controls.cxx index fc7ba629b..1ef3a5f2f 100644 --- a/src/Controls/SMESH_Controls.cxx +++ b/src/Controls/SMESH_Controls.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 // // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS @@ -35,6 +35,7 @@ #include "SMESH_MeshAlgos.hxx" #include "SMESH_OctreeNode.hxx" +#include #include #include @@ -996,6 +997,102 @@ namespace{ return aHeight; } + //================================================================================ + /*! + * \brief Standard quality of a tetrahedron but not normalized + */ + //================================================================================ + + double tetQualityByHomardMethod( const gp_XYZ & p1, + const gp_XYZ & p2, + const gp_XYZ & p3, + const gp_XYZ & p4 ) + { + gp_XYZ edgeVec[6]; + edgeVec[0] = ( p1 - p2 ); + edgeVec[1] = ( p2 - p3 ); + edgeVec[2] = ( p3 - p1 ); + edgeVec[3] = ( p4 - p1 ); + edgeVec[4] = ( p4 - p2 ); + edgeVec[5] = ( p4 - p3 ); + + double maxEdgeLen2 = edgeVec[0].SquareModulus(); + maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[1].SquareModulus() ); + maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[2].SquareModulus() ); + maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[3].SquareModulus() ); + maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[4].SquareModulus() ); + maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[5].SquareModulus() ); + double maxEdgeLen = Sqrt( maxEdgeLen2 ); + + gp_XYZ cross01 = edgeVec[0] ^ edgeVec[1]; + double sumArea = ( cross01 ).Modulus(); // actually double area + sumArea += ( edgeVec[0] ^ edgeVec[3] ).Modulus(); + sumArea += ( edgeVec[1] ^ edgeVec[4] ).Modulus(); + sumArea += ( edgeVec[2] ^ edgeVec[5] ).Modulus(); + + double sixVolume = Abs( cross01 * edgeVec[4] ); // 6 * volume + double quality = maxEdgeLen * sumArea / sixVolume; // not normalized!!! + return quality; + } + + //================================================================================ + /*! + * \brief HOMARD method of hexahedron quality + * 1. Decompose the hexa into 24 tetra: each face is splitted into 4 triangles by + * adding the diagonals and every triangle is connected to the center of the hexa. + * 2. Compute the quality of every tetra with the same formula as for the standard quality, + * except that the factor for the normalization is not the same because the final goal + * is to have a quality equal to 1 for a perfect cube. So the formula is: + * qual = max(lengthes of 6 edges) * (sum of surfaces of 4 faces) / (7.6569*6*volume) + * 3. The quality of the hexa is the highest value of the qualities of the 24 tetra + */ + //================================================================================ + + double hexQualityByHomardMethod( const TSequenceOfXYZ& P ) + { + gp_XYZ quadCenter[6]; + quadCenter[0] = ( P(1) + P(2) + P(3) + P(4) ) / 4.; + quadCenter[1] = ( P(5) + P(6) + P(7) + P(8) ) / 4.; + quadCenter[2] = ( P(1) + P(2) + P(6) + P(5) ) / 4.; + quadCenter[3] = ( P(2) + P(3) + P(7) + P(6) ) / 4.; + quadCenter[4] = ( P(3) + P(4) + P(8) + P(7) ) / 4.; + quadCenter[5] = ( P(1) + P(4) + P(8) + P(5) ) / 4.; + + gp_XYZ hexCenter = ( P(1) + P(2) + P(3) + P(4) + P(5) + P(6) + P(7) + P(8) ) / 8.; + + // quad 1 ( 1 2 3 4 ) + double quality = tetQualityByHomardMethod( P(1), P(2), quadCenter[0], hexCenter ); + quality = Max( quality, tetQualityByHomardMethod( P(2), P(3), quadCenter[0], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(3), P(4), quadCenter[0], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(4), P(1), quadCenter[0], hexCenter )); + // quad 2 ( 5 6 7 8 ) + quality = Max( quality, tetQualityByHomardMethod( P(5), P(6), quadCenter[1], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(6), P(7), quadCenter[1], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(7), P(8), quadCenter[1], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(8), P(5), quadCenter[1], hexCenter )); + // quad 3 ( 1 2 6 5 ) + quality = Max( quality, tetQualityByHomardMethod( P(1), P(2), quadCenter[2], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(2), P(6), quadCenter[2], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(6), P(5), quadCenter[2], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(5), P(1), quadCenter[2], hexCenter )); + // quad 4 ( 2 3 7 6 ) + quality = Max( quality, tetQualityByHomardMethod( P(2), P(3), quadCenter[3], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(3), P(7), quadCenter[3], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(7), P(6), quadCenter[3], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(6), P(2), quadCenter[3], hexCenter )); + // quad 5 ( 3 4 8 7 ) + quality = Max( quality, tetQualityByHomardMethod( P(3), P(4), quadCenter[4], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(4), P(8), quadCenter[4], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(8), P(7), quadCenter[4], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(7), P(3), quadCenter[4], hexCenter )); + // quad 6 ( 1 4 8 5 ) + quality = Max( quality, tetQualityByHomardMethod( P(1), P(4), quadCenter[5], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(4), P(8), quadCenter[5], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(8), P(5), quadCenter[5], hexCenter )); + quality = Max( quality, tetQualityByHomardMethod( P(5), P(1), quadCenter[5], hexCenter )); + + return quality / 7.65685424949; + } } double AspectRatio3D::GetValue( long theId ) @@ -1081,7 +1178,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) case 5:{ { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )}; - aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )}; @@ -1100,7 +1197,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) case 6:{ { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )}; - aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )}; @@ -1125,9 +1222,13 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) break; } case 8:{ + + return hexQualityByHomardMethod( P ); // bos #23982 + + { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )}; - aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality); + aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])); } { gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )}; @@ -1262,7 +1363,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P ) case 12: { gp_XYZ aXYZ[8] = {P( 1 ),P( 2 ),P( 4 ),P( 5 ),P( 7 ),P( 8 ),P( 10 ),P( 11 )}; - aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality); + aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])); } { gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )}; @@ -2044,7 +2145,7 @@ double Deflection2D::GetBadRate( double Value, int /*nbNodes*/ ) const */ //================================================================================ -double MultiConnection::GetValue( const TSequenceOfXYZ& P ) +double MultiConnection::GetValue( const TSequenceOfXYZ& /*P*/ ) { return 0; } @@ -2071,7 +2172,7 @@ SMDSAbs_ElementType MultiConnection::GetType() const */ //================================================================================ -double MultiConnection2D::GetValue( const TSequenceOfXYZ& P ) +double MultiConnection2D::GetValue( const TSequenceOfXYZ& /*P*/ ) { return 0; } @@ -4262,9 +4363,11 @@ private: bool isOutOfFace (const gp_Pnt& p); bool isOutOfEdge (const gp_Pnt& p); bool isOutOfVertex(const gp_Pnt& p); - bool isOutOfNone (const gp_Pnt& p) { return true; } + bool isOutOfNone (const gp_Pnt& /*p*/) { return true; } bool isBox (const TopoDS_Shape& s); + TopoDS_Shape prepareSolid( const TopoDS_Shape& theSolid ); + bool (Classifier::* myIsOutFun)(const gp_Pnt& p); BRepClass3d_SolidClassifier* mySolidClfr; // ptr because of a run-time forbidden copy-constructor Bnd_B3d myBox; @@ -4524,11 +4627,17 @@ bool ElementsOnShape::IsSatisfy (const SMDS_MeshElement* elem) centerXYZ /= elem->NbNodes(); isSatisfy = false; if ( myOctree ) + { + myWorkClassifiers.clear(); + myOctree->GetClassifiersAtPoint( centerXYZ, myWorkClassifiers ); for ( size_t i = 0; i < myWorkClassifiers.size() && !isSatisfy; ++i ) isSatisfy = ! myWorkClassifiers[i]->IsOut( centerXYZ ); + } else + { for ( size_t i = 0; i < myClassifiers.size() && !isSatisfy; ++i ) isSatisfy = ! myClassifiers[i].IsOut( centerXYZ ); + } } return isSatisfy; @@ -4613,7 +4722,7 @@ void ElementsOnShape::Classifier::Init( const TopoDS_Shape& theShape, } else { - mySolidClfr = new BRepClass3d_SolidClassifier(theShape); + mySolidClfr = new BRepClass3d_SolidClassifier( prepareSolid( theShape )); myIsOutFun = & ElementsOnShape::Classifier::isOutOfSolid; } break; @@ -4692,6 +4801,25 @@ ElementsOnShape::Classifier::~Classifier() delete mySolidClfr; mySolidClfr = 0; } +TopoDS_Shape ElementsOnShape::Classifier::prepareSolid( const TopoDS_Shape& theSolid ) +{ + // try to limit tolerance of theSolid down to myTol (issue #19026) + + // check if tolerance of theSolid is more than myTol + bool tolIsOk = true; // max tolerance is at VERTEXes + for ( TopExp_Explorer exp( theSolid, TopAbs_VERTEX ); exp.More() && tolIsOk; exp.Next() ) + tolIsOk = ( myTol >= BRep_Tool::Tolerance( TopoDS::Vertex( exp.Current() ))); + if ( tolIsOk ) + return theSolid; + + // make a copy to prevent the original shape from changes + TopoDS_Shape resultShape = BRepBuilderAPI_Copy( theSolid ); + + if ( !GEOMUtils::FixShapeTolerance( resultShape, TopAbs_SHAPE, myTol )) + return theSolid; + return resultShape; +} + bool ElementsOnShape::Classifier::isOutOfSolid( const gp_Pnt& p ) { if ( isOutOfBox( p )) return true;