X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_QuadToTriaAdaptor.cxx;h=2e4ea50b8fb3effa05cc7ff8904de0bbbf3d0df4;hp=f4c519ae7a63df1b8f2729542efeb14219df71b9;hb=30ce546b0c5099ad1112929e2db94810e683e54b;hpb=f6e2eed4240c426f1e65b40d1bd7e8d109a4d4b5 diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx index f4c519ae7..2e4ea50b8 100644 --- a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx +++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2016 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 @@ -24,12 +24,16 @@ #include "StdMeshers_QuadToTriaAdaptor.hxx" +#include "SMDS_IteratorOnIterators.hxx" #include "SMDS_SetIterator.hxx" #include "SMESHDS_GroupBase.hxx" +#include "SMESHDS_Mesh.hxx" #include "SMESH_Algo.hxx" #include "SMESH_Group.hxx" +#include "SMESH_Mesh.hxx" #include "SMESH_MeshAlgos.hxx" #include "SMESH_MesherHelper.hxx" +#include "SMESH_subMesh.hxx" #include #include @@ -38,6 +42,7 @@ #include #include #include +#include #include #include @@ -115,20 +120,20 @@ namespace gp_Vec nJ = baseVec.Crossed( baJ ); // Check angle between normals - double angle = nI.Angle( nJ ); + double angle = nI.Angle( nJ ); bool tooClose = ( angle < 15. * M_PI / 180. ); // Check if pyramids collide if ( !tooClose && ( baI * baJ > 0 ) && ( nI * nJ > 0 )) { // find out if nI points outside of PrmI or inside - int dInd = baseNodesIndI[1] - baseNodesIndI[0]; + int dInd = baseNodesIndI[1] - baseNodesIndI[0]; bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0; - // find out sign of projection of nJ to baI + // find out sign of projection of baI to nJ double proj = baI * nJ; - tooClose = isOutI ? proj > 0 : proj < 0; + tooClose = ( isOutI ? proj > 0 : proj < 0 ); } // Check if PrmI and PrmJ are in same domain @@ -170,8 +175,9 @@ namespace continue; // f is a base quadrangle // check projections of face direction (baOFN) to triange normals (nI and nJ) - gp_Vec baOFN( base1, SMESH_TNodeXYZ( otherFaceNode )); - if ( nI * baOFN > 0 && nJ * baOFN > 0 ) + gp_Vec baOFN( base2, SMESH_TNodeXYZ( otherFaceNode )); + if ( nI * baOFN > 0 && nJ * baOFN > 0 && + baI* baOFN > 0 && baJ* baOFN > 0 ) // issue 0023212 { tooClose = false; // f is between pyramids break; @@ -254,6 +260,51 @@ namespace } } + //================================================================================ + /*! + * \brief Store an error about overlapping faces + */ + //================================================================================ + + bool overlapError( SMESH_Mesh& mesh, + const SMDS_MeshElement* face1, + const SMDS_MeshElement* face2, + const TopoDS_Shape& shape = TopoDS_Shape()) + { + if ( !face1 || !face2 ) return false; + + SMESH_Comment msg; + msg << "face " << face1->GetID() << " overlaps face " << face2->GetID(); + + SMESH_subMesh * sm = 0; + if ( shape.IsNull() ) + { + sm = mesh.GetSubMesh( mesh.GetShapeToMesh() ); + } + else if ( shape.ShapeType() >= TopAbs_SOLID ) + { + sm = mesh.GetSubMesh( shape ); + } + else + { + TopoDS_Iterator it ( shape ); + if ( it.More() ) + sm = mesh.GetSubMesh( it.Value() ); + } + if ( sm ) + { + SMESH_ComputeErrorPtr& err = sm->GetComputeError(); + if ( !err || err->IsOK() ) + { + err = SMESH_ComputeError::New( COMPERR_BAD_INPUT_MESH, msg, sm->GetAlgo() ); + err->myBadElements.push_back( face1 ); + err->myBadElements.push_back( face2 ); + } + } + //throw SALOME_Exception( msg.c_str() ); + + return false; // == "algo fails" + } } //================================================================================ @@ -266,6 +317,8 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr const SMDS_MeshElement* PrmJ, set & nodesToMove) { + // cout << endl << "Merge " << PrmI->GetID() << " " << PrmJ->GetID() << " " + // << PrmI->GetNode(4) << PrmJ->GetNode(4) << endl; const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove //int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume ); SMESH_TNodeXYZ Pj( Nrem ); @@ -288,7 +341,7 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr vector< const SMDS_MeshElement* > inverseElems // copy inverse elements to avoid iteration on changing container ( TStdElemIterator( CommonNode->GetInverseElementIterator(SMDSAbs_Face)), itEnd); - for ( unsigned i = 0; i < inverseElems.size(); ++i ) + for ( size_t i = 0; i < inverseElems.size(); ++i ) { const SMDS_MeshElement* FI = inverseElems[i]; const SMDS_MeshElement* FJEqual = 0; @@ -309,11 +362,12 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr } // set the common apex node to pyramids and triangles merged with J + vector< const SMDS_MeshNode* > nodes; inverseElems.assign( TStdElemIterator( Nrem->GetInverseElementIterator()), itEnd ); - for ( unsigned i = 0; i < inverseElems.size(); ++i ) + for ( size_t i = 0; i < inverseElems.size(); ++i ) { const SMDS_MeshElement* elem = inverseElems[i]; - vector< const SMDS_MeshNode* > nodes( elem->begin_nodes(), elem->end_nodes() ); + nodes.assign( elem->begin_nodes(), elem->end_nodes() ); nodes[ elem->GetType() == SMDSAbs_Volume ? PYRAM_APEX : TRIA_APEX ] = CommonNode; GetMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size()); } @@ -330,33 +384,34 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr //================================================================================ void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement* PrmI, - set& nodesToMove) + set& nodesToMove, + const bool isRecursion) { TIDSortedElemSet adjacentPyrams; bool mergedPyrams = false; - for(int k=0; k<4; k++) // loop on 4 base nodes of PrmI + for ( int k=0; k<4; k++ ) // loop on 4 base nodes of PrmI { - const SMDS_MeshNode* n = PrmI->GetNode(k); + const SMDS_MeshNode* n = PrmI->GetNode(k); SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); while ( vIt->more() ) { const SMDS_MeshElement* PrmJ = vIt->next(); - if ( PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second ) + if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second ) continue; - if ( PrmI != PrmJ && TooCloseAdjacent( PrmI, PrmJ, GetMesh()->HasShapeToMesh() )) + if ( TooCloseAdjacent( PrmI, PrmJ, GetMesh()->HasShapeToMesh() )) { MergePiramids( PrmI, PrmJ, nodesToMove ); mergedPyrams = true; // container of inverse elements can change - vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); + // vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); -- iterator re-implemented } } } - if ( mergedPyrams ) + if ( mergedPyrams && !isRecursion ) { TIDSortedElemSet::iterator prm; for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm) - MergeAdjacent( *prm, nodesToMove ); + MergeAdjacent( *prm, nodesToMove, true ); } } @@ -395,17 +450,17 @@ static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& PC, const gp_Vec& V) { gp_Pnt Pbest = PC; - const double a = P1.Distance(P2); - const double b = P1.Distance(PC); - const double c = P2.Distance(PC); - if( a < (b+c)/2 ) + const double a2 = P1.SquareDistance(P2); + const double b2 = P1.SquareDistance(PC); + const double c2 = P2.SquareDistance(PC); + if ( a2 < ( b2 + Sqrt( 4 * b2 * c2 ) + c2 ) / 4 ) // ( a < (b+c)/2 ) return Pbest; else { // find shift along V in order a to became equal to (b+c)/2 const double Vsize = V.Magnitude(); if ( fabs( Vsize ) > std::numeric_limits::min() ) { - const double shift = sqrt( a*a + (b*b-c*c)*(b*b-c*c)/16/a/a - (b*b+c*c)/2 ); + const double shift = sqrt( a2 + (b2-c2)*(b2-c2)/16/a2 - (b2+c2)/2 ); Pbest.ChangeCoord() += shift * V.XYZ() / Vsize; } } @@ -414,79 +469,59 @@ static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2, //======================================================================= //function : HasIntersection3 -//purpose : Auxilare for HasIntersection() -// find intersection point between triangle (P1,P2,P3) -// and segment [PC,P] +//purpose : Find intersection point between a triangle (P1,P2,P3) +// and a segment [PC,P] //======================================================================= static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3) { - //cout<<"HasIntersection3"< preci ) || - ( (PC.Y()-PIn.Y())*(P.Y()-PIn.Y()) > preci ) || - ( (PC.Z()-PIn.Z())*(P.Z()-PIn.Z()) > preci ); - if(IsExternal) { - return false; - } - // check if this point is internal for triangle (P1,P2,P3) - gp_Vec V1(PIn,P1); - gp_Vec V2(PIn,P2); - gp_Vec V3(PIn,P3); - if( V1.Magnitude() -ANGL_EPSILON && det < ANGL_EPSILON ) + return false; + + /* calculate distance from vert0 to ray origin */ + gp_XYZ tvec = orig - vert0; + + /* calculate U parameter and test bounds */ + double u = ( tvec * pvec ) / det; + //if (u < 0.0 || u > 1.0) + if (u < -EPSILON || u > 1.0 + EPSILON) + return false; + + /* prepare to test V parameter */ + gp_XYZ qvec = tvec ^ edge1; + + /* calculate V parameter and test bounds */ + double v = (dir * qvec) / det; + //if ( v < 0.0 || u + v > 1.0 ) + if ( v < -EPSILON || u + v > 1.0 + EPSILON) + return false; + + /* calculate t, ray intersects triangle */ + double t = (edge2 * qvec) / det; + + Pint = orig + dir * t; + + return ( t > 0. && t < segLen ); } //======================================================================= @@ -502,15 +537,15 @@ static bool HasIntersection(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, } else { bool check = false; - if( (aContour(1).Distance(aContour(2)) > 1.e-6) && - (aContour(1).Distance(aContour(3)) > 1.e-6) && - (aContour(2).Distance(aContour(3)) > 1.e-6) ) { + if( (aContour(1).SquareDistance(aContour(2)) > 1.e-12) && + (aContour(1).SquareDistance(aContour(3)) > 1.e-12) && + (aContour(2).SquareDistance(aContour(3)) > 1.e-12) ) { check = HasIntersection3( P, PC, Pint, aContour(1), aContour(2), aContour(3) ); } if(check) return true; - if( (aContour(1).Distance(aContour(4)) > 1.e-6) && - (aContour(1).Distance(aContour(3)) > 1.e-6) && - (aContour(4).Distance(aContour(3)) > 1.e-6) ) { + if( (aContour(1).SquareDistance(aContour(4)) > 1.e-12) && + (aContour(1).SquareDistance(aContour(3)) > 1.e-12) && + (aContour(4).SquareDistance(aContour(3)) > 1.e-12) ) { check = HasIntersection3( P, PC, Pint, aContour(1), aContour(3), aContour(4) ); } if(check) return true; @@ -521,54 +556,112 @@ static bool HasIntersection(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, //================================================================================ /*! - * \brief Checks if a line segment (P,PC) intersects any mesh face. - * \param P - first segment end - * \param PC - second segment end (it is a gravity center of quadrangle) - * \param Pint - (out) intersection point + * \brief Return allowed height of a pyramid + * \param Papex - optimal pyramid apex + * \param PC - gravity center of a quadrangle + * \param PN - four nodes of the quadrangle * \param aMesh - mesh - * \param aShape - shape to check faces on - * \param NotCheckedFace - mesh face not to check - * \retval bool - true if there is an intersection + * \param NotCheckedFace - the quadrangle face + * \param Shape - the shape being meshed + * \retval false if mesh invalidity detected */ //================================================================================ -bool StdMeshers_QuadToTriaAdaptor::CheckIntersection (const gp_Pnt& P, - const gp_Pnt& PC, - gp_Pnt& Pint, - SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape, - const SMDS_MeshElement* NotCheckedFace) +bool StdMeshers_QuadToTriaAdaptor::LimitHeight (gp_Pnt& Papex, + const gp_Pnt& PC, + const TColgp_Array1OfPnt& PN, + const vector& FNodes, + SMESH_Mesh& aMesh, + const SMDS_MeshElement* NotCheckedFace, + const bool UseApexRay, + const TopoDS_Shape& Shape) { if ( !myElemSearcher ) myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *aMesh.GetMeshDS() ); SMESH_ElementSearcher* searcher = const_cast(myElemSearcher); - bool res = false; - double dist = RealLast(); // find intersection closest to PC - gp_Pnt Pres; - - gp_Ax1 line( P, gp_Vec(P,PC)); - vector< const SMDS_MeshElement* > suspectElems; - searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems); + // Find intersection of faces with (P,PC) segment elongated 3 times + double height = Papex.Distance( PC ); + gp_Ax1 line( PC, gp_Vec( PC, Papex )); + gp_Pnt Pint, Ptest; + vector< const SMDS_MeshElement* > suspectFaces; TColgp_SequenceOfPnt aContour; - for ( size_t iF = 0; iF < suspectElems.size(); ++iF ) + + if ( UseApexRay ) { - const SMDS_MeshElement* face = suspectElems[iF]; + double idealHeight = height; + const SMDS_MeshElement* intFace = 0; + + // find intersection closest to PC + Ptest = PC.XYZ() + line.Direction().XYZ() * height * 3; + + searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces ); + for ( size_t iF = 0; iF < suspectFaces.size(); ++iF ) + { + const SMDS_MeshElement* face = suspectFaces[iF]; + if ( face == NotCheckedFace ) continue; + + aContour.Clear(); + for ( int i = 0, nb = face->NbCornerNodes(); i < nb; ++i ) + aContour.Append( SMESH_TNodeXYZ( face->GetNode(i) )); + + if ( HasIntersection( Ptest, PC, Pint, aContour )) + { + double dInt = PC.Distance( Pint ) / 3.; + if ( dInt < height ) + { + height = dInt; + intFace = face; + } + } + } + if ( height < 1e-2 * idealHeight && intFace ) + return overlapError( aMesh, NotCheckedFace, intFace, Shape ); + } + + // Find faces intersecting triangular facets of the pyramid (issue 23212) + + gp_XYZ center = PC.XYZ() + line.Direction().XYZ() * height * 0.5; + double diameter = Max( PN(1).Distance(PN(3)), PN(2).Distance(PN(4))); + suspectFaces.clear(); + searcher->GetElementsInSphere( center, diameter * 0.6, SMDSAbs_Face, suspectFaces); + + const double upShift = 1.5; + Ptest = PC.XYZ() + line.Direction().XYZ() * height * upShift; // tmp apex + + for ( size_t iF = 0; iF < suspectFaces.size(); ++iF ) + { + const SMDS_MeshElement* face = suspectFaces[iF]; if ( face == NotCheckedFace ) continue; - aContour.Clear(); - for ( int i = 0; i < face->NbCornerNodes(); ++i ) - aContour.Append( SMESH_TNodeXYZ( face->GetNode(i) )); - if ( HasIntersection(P, PC, Pres, aContour)) { - res = true; - double tmp = PC.Distance(Pres); - if ( tmp < dist ) { - Pint = Pres; - dist = tmp; + if ( face->GetNodeIndex( FNodes[0] ) >= 0 || + face->GetNodeIndex( FNodes[1] ) >= 0 || + face->GetNodeIndex( FNodes[2] ) >= 0 || + face->GetNodeIndex( FNodes[3] ) >= 0 ) + continue; // neighbor face of the quadrangle + + // limit height using points of intersection of face links with pyramid facets + int nbN = face->NbCornerNodes(); + gp_Pnt P1 = SMESH_TNodeXYZ( face->GetNode( nbN-1 )); // 1st link end + for ( int i = 0; i < nbN; ++i ) + { + gp_Pnt P2 = SMESH_TNodeXYZ( face->GetNode(i) ); // 2nd link end + + for ( int iN = 1; iN <= 4; ++iN ) // loop on pyramid facets + { + if ( HasIntersection3( P1, P2, Pint, PN(iN), PN(iN+1), Ptest )) + { + height = Min( height, gp_Vec( PC, Pint ) * line.Direction() ); + //Ptest = PC.XYZ() + line.Direction().XYZ() * height * upShift; // new tmp apex + } } + P1 = P2; } } - return res; + + Papex = PC.XYZ() + line.Direction().XYZ() * height; + + return true; } //================================================================================ @@ -720,25 +813,36 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, vector myPyramids; + const SMESHDS_SubMesh * aSubMeshDSFace; SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); SMESH_MesherHelper helper(aMesh); helper.IsQuadraticSubMesh(aShape); helper.SetElementsOnShape( true ); if ( myElemSearcher ) delete myElemSearcher; + vector< SMDS_ElemIteratorPtr > itVec; if ( aProxyMesh ) - myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, aProxyMesh->GetFaces(aShape)); + { + itVec.push_back( aProxyMesh->GetFaces( aShape )); + } else - myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS ); + { + for ( TopExp_Explorer exp(aShape,TopAbs_FACE); exp.More(); exp.Next() ) + if (( aSubMeshDSFace = meshDS->MeshElements( exp.Current() ))) + itVec.push_back( aSubMeshDSFace->GetElements() ); + } + typedef + SMDS_IteratorOnIterators< const SMDS_MeshElement*, vector< SMDS_ElemIteratorPtr > > TIter; + SMDS_ElemIteratorPtr faceIt( new TIter( itVec )); + myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, faceIt ); - const SMESHDS_SubMesh * aSubMeshDSFace; TColgp_Array1OfPnt PN(1,5); TColgp_Array1OfVec VN(1,4); vector FNodes(5); gp_Pnt PC; gp_Vec VNorm; - for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) + for ( TopExp_Explorer exp(aShape,TopAbs_FACE); exp.More(); exp.Next() ) { const TopoDS_Shape& aShapeFace = exp.Current(); if ( aProxyMesh ) @@ -809,17 +913,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, } else { // check possible intersection with other faces - gp_Pnt Pint; - gp_Vec VB(PC,PCbest); - gp_Pnt PCbestTmp = PC.XYZ() + VB.XYZ() * 3.0; - bool hasInters = CheckIntersection(PCbestTmp, PC, Pint, aMesh, aShape, face); - if ( hasInters ) { - double dist = PC.Distance(Pint)/3.; - if ( dist < height ) { - gp_Dir aDir( VB ); - PCbest = PC.XYZ() + aDir.XYZ() * dist; - } - } + if ( !LimitHeight( PCbest, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/true, aShape )) + return false; } // create node for PCbest SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() ); @@ -894,30 +989,31 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) SMESHDS_GroupBase* groupDS = 0; SMESH_Mesh::GroupIteratorPtr groupIt = aMesh.GetGroups(); while ( groupIt->more() ) + { + groupDS = 0; + SMESH_Group * group = groupIt->next(); + if ( !group ) continue; + groupDS = group->GetGroupDS(); + if ( !groupDS || groupDS->IsEmpty() ) { groupDS = 0; - SMESH_Group * group = groupIt->next(); - if ( !group ) continue; - groupDS = group->GetGroupDS(); - if ( !groupDS || groupDS->IsEmpty() ) - { - groupDS = 0; - continue; - } - if (groupDS->GetType() != SMDSAbs_Face) - { - groupDS = 0; - continue; - } - std::string grpName = group->GetName(); - if (grpName == groupName) - { - MESSAGE("group skinFaces provided"); - break; - } - else - groupDS = 0; + continue; } + if (groupDS->GetType() != SMDSAbs_Face) + { + groupDS = 0; + continue; + } + std::string grpName = group->GetName(); + if (grpName == groupName) + { + break; + } + else + groupDS = 0; + } + + const bool toFindVolumes = aMesh.NbVolumes() > 0; vector myPyramids; SMESH_MesherHelper helper(aMesh); @@ -930,6 +1026,9 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) if ( !myElemSearcher ) myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS ); SMESH_ElementSearcher* searcher = const_cast(myElemSearcher); + SMESHUtils::Deleter + volSearcher( SMESH_MeshAlgos::GetElementSearcher( *meshDS )); + vector< const SMDS_MeshElement* > suspectFaces, foundVolumes; TColgp_Array1OfPnt PN(1,5); TColgp_Array1OfVec VN(1,4); @@ -949,7 +1048,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) if ( what == NOT_QUAD ) continue; if ( volumes[0] && volumes[1] ) - continue; // face is shared by two volumes - no space for a pyramid + continue; // face is shared by two volumes - no room for a pyramid if ( what == DEGEN_QUAD ) { @@ -971,11 +1070,11 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) gp_Pnt Pres1,Pres2; gp_Ax1 line( PC, VNorm ); - vector< const SMDS_MeshElement* > suspectElems; - searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems); + vector< const SMDS_MeshElement* > suspectFaces; + searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces); - for ( size_t iF = 0; iF < suspectElems.size(); ++iF ) { - const SMDS_MeshElement* F = suspectElems[iF]; + for ( size_t iF = 0; iF < suspectFaces.size(); ++iF ) { + const SMDS_MeshElement* F = suspectFaces[iF]; if ( F == face ) continue; aContour.Clear(); for ( int i = 0; i < 4; ++i ) @@ -1026,7 +1125,9 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) continue; } + // ----------------------------------- // Case of non-degenerated quadrangle + // ----------------------------------- // Find pyramid peak @@ -1048,6 +1149,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) } // Restrict pyramid height by intersection with other faces + gp_Vec tmpDir(PC,PCbest); tmpDir.Normalize(); double tmp = PN(1).Distance(PN(3)) + PN(2).Distance(PN(4)); // far points: in (PC, PCbest) direction and vice-versa @@ -1055,16 +1157,33 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) PC.XYZ() - tmpDir.XYZ() * tmp * 1.e6 }; // check intersection for farPnt1 and farPnt2 bool intersected[2] = { false, false }; - double dist [2] = { RealLast(), RealLast() }; - gp_Pnt intPnt[2]; + double dist2int [2] = { RealLast(), RealLast() }; + gp_Pnt intPnt [2]; + int intFaceInd [2] = { 0, 0 }; + + if ( toFindVolumes && 0 ) // non-conformal mesh is not suitable for any mesher so far + { + // there are volumes in the mesh, in a non-conformal mesh an neighbor + // volume can be not found yet + for ( int isRev = 0; isRev < 2; ++isRev ) + { + if ( volumes[isRev] ) continue; + gp_Pnt testPnt = PC.XYZ() + tmpDir.XYZ() * height * ( isRev ? -0.1: 0.1 ); + foundVolumes.clear(); + if ( volSearcher->FindElementsByPoint( testPnt, SMDSAbs_Volume, foundVolumes )) + volumes[isRev] = foundVolumes[0]; + } + if ( volumes[0] && volumes[1] ) + continue; // no room for a pyramid + } gp_Ax1 line( PC, tmpDir ); - vector< const SMDS_MeshElement* > suspectElems; - searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems); + suspectFaces.clear(); + searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces); - for ( size_t iF = 0; iF < suspectElems.size(); ++iF ) + for ( size_t iF = 0; iF < suspectFaces.size(); ++iF ) { - const SMDS_MeshElement* F = suspectElems[iF]; + const SMDS_MeshElement* F = suspectFaces[iF]; if ( F == face ) continue; aContour.Clear(); int nbN = F->NbCornerNodes(); @@ -1073,31 +1192,33 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) gp_Pnt intP; for ( int isRev = 0; isRev < 2; ++isRev ) { - if( !volumes[isRev] && HasIntersection(farPnt[isRev], PC, intP, aContour) ) { - intersected[isRev] = true; + if( !volumes[isRev] && HasIntersection(farPnt[isRev], PC, intP, aContour) ) + { double d = PC.Distance( intP ); - if( d < dist[isRev] ) + if ( d < dist2int[isRev] ) { - intPnt[isRev] = intP; - dist [isRev] = d; + intersected[isRev] = true; + intPnt [isRev] = intP; + dist2int [isRev] = d; + intFaceInd [isRev] = iF; } } } } // if the face belong to the group of skinFaces, do not build a pyramid outside - if (groupDS && groupDS->Contains(face)) + if ( groupDS && groupDS->Contains(face) ) { intersected[0] = false; } else if ( intersected[0] && intersected[1] ) // check if one of pyramids is in a hole { - gp_Pnt P ( PC.XYZ() + tmpDir.XYZ() * 0.5 * PC.Distance( intPnt[0] )); + gp_Pnt P ( PC.XYZ() + tmpDir.XYZ() * 0.5 * dist2int[0] ); if ( searcher->GetPointState( P ) == TopAbs_OUT ) intersected[0] = false; else { - P = ( PC.XYZ() - tmpDir.XYZ() * 0.5 * PC.Distance( intPnt[1] )); + P = ( PC.XYZ() - tmpDir.XYZ() * 0.5 * dist2int[1] ); if ( searcher->GetPointState( P ) == TopAbs_OUT ) intersected[1] = false; } @@ -1107,15 +1228,20 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) for ( int isRev = 0; isRev < 2; ++isRev ) { - if( !intersected[isRev] ) continue; - double pyramidH = Min( height, PC.Distance(intPnt[isRev])/3.); - PCbest = PC.XYZ() + tmpDir.XYZ() * (isRev ? -pyramidH : pyramidH); + if ( !intersected[isRev] ) continue; + double pyramidH = Min( height, dist2int[isRev]/3. ); + gp_Pnt Papex = PC.XYZ() + tmpDir.XYZ() * (isRev ? -pyramidH : pyramidH); + if ( pyramidH < 1e-2 * height ) + return overlapError( aMesh, face, suspectFaces[ intFaceInd[isRev] ] ); + + if ( !LimitHeight( Papex, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/false )) + return false; - // create node for PCbest - SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() ); + // create node for Papex + SMDS_MeshNode* NewNode = helper.AddNode( Papex.X(), Papex.Y(), Papex.Z() ); // add triangles to result map - for(i=0; i<4; i++) { + for ( i = 0; i < 4; i++) { SMDS_MeshFace* NewFace; if(isRev) NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ); @@ -1146,16 +1272,19 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh, const vector& myPyramids) { - if(myPyramids.empty()) + if ( myPyramids.empty() ) return true; SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); size_t i, j, k; - int myShapeID = myPyramids[0]->GetNode(4)->getshapeId(); - - if ( myElemSearcher ) delete myElemSearcher; - myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS ); - SMESH_ElementSearcher* searcher = const_cast(myElemSearcher); + //int myShapeID = myPyramids[0]->GetNode(4)->getshapeId(); + { + SMDS_ElemIteratorPtr + pyramIt( new SMDS_ElementVectorIterator( myPyramids.begin(), myPyramids.end() )); + if ( myElemSearcher ) delete myElemSearcher; + myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, pyramIt ); + } + SMESH_ElementSearcher* searcher = const_cast( myElemSearcher ); set nodesToMove; @@ -1167,17 +1296,19 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& MergeAdjacent( PrmI, nodesToMove ); } - // iterate on all pyramids + // iterate on all new pyramids + vector< const SMDS_MeshElement* > suspectPyrams; for ( i = 0; i < myPyramids.size(); ++i ) { - const SMDS_MeshElement* PrmI = myPyramids[i]; + const SMDS_MeshElement* PrmI = myPyramids[i]; + const SMDS_MeshNode* apexI = PrmI->GetNode( PYRAM_APEX ); // compare PrmI with all the rest pyramids // collect adjacent pyramids and nodes coordinates of PrmI set checkedPyrams; - vector PsI(5); - for(k=0; k<5; k++) // loop on 4 base nodes of PrmI + gp_Pnt PsI[5]; + for ( k = 0; k < 5; k++ ) { const SMDS_MeshNode* n = PrmI->GetNode(k); PsI[k] = SMESH_TNodeXYZ( n ); @@ -1190,70 +1321,77 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& } } + // get pyramids to check + gp_XYZ PC = ( PsI[0].XYZ() + PsI[1].XYZ() + PsI[2].XYZ() + PsI[3].XYZ() ) / 4.; + gp_XYZ ray = PsI[4].XYZ() - PC; + gp_XYZ center = PC + 0.5 * ray; + double diameter = Max( PsI[0].Distance(PsI[2]), PsI[1].Distance(PsI[3])); + suspectPyrams.clear(); + searcher->GetElementsInSphere( center, diameter * 0.6, SMDSAbs_Volume, suspectPyrams); + // check intersection with distant pyramids - for(k=0; k<4; k++) // loop on 4 base nodes of PrmI + for ( j = 0; j < suspectPyrams.size(); ++j ) { - gp_Vec Vtmp(PsI[k],PsI[4]); - gp_Ax1 line( PsI[k], Vtmp ); - vector< const SMDS_MeshElement* > suspectPyrams; - searcher->GetElementsNearLine( line, SMDSAbs_Volume, suspectPyrams); + const SMDS_MeshElement* PrmJ = suspectPyrams[j]; + if ( PrmJ == PrmI ) + continue; + if ( apexI == PrmJ->GetNode( PYRAM_APEX )) + continue; // pyramids PrmI and PrmJ already merged + if ( !checkedPyrams.insert( PrmJ ).second ) + continue; // already checked - for ( j = 0; j < suspectPyrams.size(); ++j ) - { - const SMDS_MeshElement* PrmJ = suspectPyrams[j]; - if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 ) - continue; - if ( myShapeID != PrmJ->GetNode(4)->getshapeId()) - continue; // pyramid from other SOLID - if ( PrmI->GetNode(4) == PrmJ->GetNode(4) ) - continue; // pyramids PrmI and PrmJ already merged - if ( !checkedPyrams.insert( PrmJ ).second ) - continue; // already checked - - TXyzIterator xyzIt( PrmJ->nodesIterator() ); - vector PsJ( xyzIt, TXyzIterator() ); + gp_Pnt PsJ[5]; + for ( k = 0; k < 5; k++ ) + PsJ[k] = SMESH_TNodeXYZ( PrmJ->GetNode(k) ); + if ( ray * ( PsJ[4].XYZ() - PC ) < 0. ) + continue; // PrmJ is below PrmI + + for ( k = 0; k < 4; k++ ) // loop on 4 base nodes of PrmI + { gp_Pnt Pint; bool hasInt=false; - for(k=0; k<4 && !hasInt; k++) { - gp_Vec Vtmp(PsI[k],PsI[4]); + for ( k = 0; k < 4 && !hasInt; k++ ) + { + gp_Vec Vtmp( PsI[k], PsI[ PYRAM_APEX ]); gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01; // base node moved a bit to apex hasInt = - ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[4]) || - HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[4]) || - HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[4]) || - HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[4]) ); + ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[PYRAM_APEX]) || + HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[PYRAM_APEX]) || + HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[PYRAM_APEX]) || + HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[PYRAM_APEX]) ); } - for(k=0; k<4 && !hasInt; k++) { - gp_Vec Vtmp(PsJ[k],PsJ[4]); + for ( k = 0; k < 4 && !hasInt; k++ ) + { + gp_Vec Vtmp( PsJ[k], PsJ[ PYRAM_APEX ]); gp_Pnt Pshift = PsJ[k].XYZ() + Vtmp.XYZ() * 0.01; hasInt = - ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[4]) || - HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[4]) || - HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[4]) || - HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[4]) ); + ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[PYRAM_APEX]) || + HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[PYRAM_APEX]) || + HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[PYRAM_APEX]) || + HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[PYRAM_APEX]) ); } if ( hasInt ) { // count common nodes of base faces of two pyramids int nbc = 0; - for (k=0; k<4; k++) + for ( k = 0; k < 4; k++ ) nbc += int ( PrmI->GetNodeIndex( PrmJ->GetNode(k) ) >= 0 ); if ( nbc == 4 ) continue; // pyrams have a common base face - if(nbc>0) + if ( nbc > 0 ) { // Merge the two pyramids and others already merged with them MergePiramids( PrmI, PrmJ, nodesToMove ); } - else { // nbc==0 - + else // nbc==0 + { // decrease height of pyramids gp_XYZ PCi(0,0,0), PCj(0,0,0); - for(k=0; k<4; k++) { + for ( k = 0; k < 4; k++ ) { PCi += PsI[k].XYZ(); PCj += PsJ[k].XYZ(); } @@ -1272,9 +1410,9 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& VN1.Scale(coef1); VN2.Scale(coef2); - SMDS_MeshNode* aNode1 = const_cast(PrmI->GetNode(4)); + SMDS_MeshNode* aNode1 = const_cast( apexI ); aNode1->setXYZ( PCi.X()+VN1.X(), PCi.Y()+VN1.Y(), PCi.Z()+VN1.Z() ); - SMDS_MeshNode* aNode2 = const_cast(PrmJ->GetNode(4)); + SMDS_MeshNode* aNode2 = const_cast(PrmJ->GetNode( PYRAM_APEX )); aNode2->setXYZ( PCj.X()+VN2.X(), PCj.Y()+VN2.Y(), PCj.Z()+VN2.Z() ); nodesToMove.insert( aNode1 ); nodesToMove.insert( aNode2 );