X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_QuadToTriaAdaptor.cxx;h=9afac27333671429736827228f0ba52d197f49b6;hp=f4c519ae7a63df1b8f2729542efeb14219df71b9;hb=cfa95f04778f656d7b5805b920f5e85d58631b8e;hpb=59627b07d70f4caa4c768be6805334d2610fa54c diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx index f4c519ae7..9afac2733 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-2019 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 @@ -169,9 +174,10 @@ namespace PrmJ->GetNodeIndex( otherFaceNode ) >= 0 )) 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 ) + // check projections of face direction (baOFN) to triangle normals (nI and nJ) + 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; @@ -252,8 +258,86 @@ namespace meshDS->RemoveFreeNode( nodesToRemove[i], sm, /*fromGroups=*/false); } } + return; } + //================================================================================ + /*! + * \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() ) + { + SMESH_BadInputElements* badElems = + new SMESH_BadInputElements( mesh.GetMeshDS(),COMPERR_BAD_INPUT_MESH, msg, sm->GetAlgo() ); + badElems->add( face1 ); + badElems->add( face2 ); + err.reset( badElems ); + } + } + + return false; // == "algo fails" + } + + //================================================================================ + /*! + * \brief Check if a face is in a SOLID + */ + //================================================================================ + + bool isInSolid( vector & faceNodes, + const int nbNodes, + const int solidID ) + { + if ( !faceNodes[0] ) + return true; // NOT_QUAD + for ( int i = 0; i < nbNodes; ++i ) + { + int shapeID = faceNodes[i]->GetShapeID(); + if ( shapeID == solidID ) + return true; + } + faceNodes.resize( nbNodes ); + std::vector vols; + SMDS_Mesh::GetElementsByNodes( faceNodes, vols, SMDSAbs_Volume ); + bool inSolid = false; + for ( size_t i = 0; i < vols.size() && !inSolid; ++i ) + { + int shapeID = vols[i]->GetShapeID(); + inSolid = ( shapeID == solidID ); + } + faceNodes.push_back( faceNodes[0] ); + return inSolid; + } } //================================================================================ @@ -266,6 +350,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)->GetID() << " " << PrmJ->GetNode(4)->GetID() << endl; const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove //int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume ); SMESH_TNodeXYZ Pj( Nrem ); @@ -284,11 +370,14 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator; TStdElemIterator itEnd; + typedef std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > TNNMap; + TNNMap mediumReplaceMap; + // find and remove coincided faces of merged pyramids 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; @@ -301,6 +390,11 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr } if ( FJEqual ) { + if ( FJEqual->NbNodes() == 6 ) // find medium nodes to replace + { + mediumReplaceMap.insert( std::make_pair( FJEqual->GetNode(3), FI->GetNode(5) )); + mediumReplaceMap.insert( std::make_pair( FJEqual->GetNode(5), FI->GetNode(3) )); + } removeTmpElement( FI ); removeTmpElement( FJEqual ); myRemovedTrias.insert( FI ); @@ -309,18 +403,34 @@ 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; + if ( !mediumReplaceMap.empty() ) + for ( size_t iN = elem->NbCornerNodes(); iN < nodes.size(); ++iN ) + { + TNNMap::iterator n2n = mediumReplaceMap.find( nodes[iN] ); + if ( n2n != mediumReplaceMap.end() ) + nodes[iN] = n2n->second; + } GetMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size()); } ASSERT( Nrem->NbInverseElements() == 0 ); GetMeshDS()->RemoveFreeNode( Nrem, GetMeshDS()->MeshElements( Nrem->getshapeId()), /*fromGroups=*/false); + if ( !mediumReplaceMap.empty() ) + for ( TNNMap::iterator n2n = mediumReplaceMap.begin(); n2n != mediumReplaceMap.end(); ++n2n ) + { + const SMDS_MeshNode* remNode = n2n->first; + if ( !remNode->IsNull() && remNode->NbInverseElements() == 0 ) + GetMeshDS()->RemoveFreeNode( remNode, 0, /*fromGroups=*/false); + } + return; } //================================================================================ @@ -330,34 +440,36 @@ 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 ); } + return; } //================================================================================ @@ -395,17 +507,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 +526,79 @@ 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, +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"< -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; + + bool hasInt = ( t > 0. && t < segLen ); + + if ( hasInt && det < EPSILON ) // t is inaccurate, additionally check + { + gp_XYZ triNorm = edge1 ^ edge2; + gp_XYZ int0vec = Pint.XYZ() - vert0; + gp_XYZ in = triNorm ^ edge1; // dir inside triangle from edge1 + double dot = int0vec * in; + if ( dot < 0 && dot / triNorm.Modulus() < -EPSILON ) + return false; + in = edge2 ^ triNorm; + dot = int0vec * in; + if ( dot < 0 && dot / triNorm.Modulus() < -EPSILON ) + return false; + gp_XYZ int1vec = Pint.XYZ() - vert1; + in = triNorm ^ ( vert2 - vert1 ); + dot = int1vec * in; + if ( dot < 0 && dot / triNorm.Modulus() < -EPSILON ) return false; - if( IAICQ.NbPoints() == 1 ) { - gp_Pnt PIn = IAICQ.Point(1); - const double preci = 1.e-10 * P.Distance(PC); - // check if this point is internal for segment [PC,P] - bool IsExternal = - ( (PC.X()-PIn.X())*(P.X()-PIn.X()) > 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() 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,59 +633,117 @@ 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 ) + { + 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 = suspectElems[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; } //================================================================================ /*! - * \brief Prepare data for the given face + * \brief Retrieve data of the given face * \param PN - coordinates of face nodes * \param VN - cross products of vectors (PC-PN(i)) ^ (PC-PN(i+1)) * \param FNodes - face nodes @@ -612,7 +782,8 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement* face int nbp = 4; int j = 0; - for(i=1; i<4; i++) { + for ( i = 1; i < 4; i++ ) + { j = i+1; for(; j<=4; j++) { if( PN(i).Distance(PN(j)) < 1.e-6 ) @@ -620,11 +791,10 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement* face } if(j<=4) break; } - //int deg_num = IsDegenarate(PN); - //if(deg_num>0) { + bool hasdeg = false; - if(i<4) { - //cout<<"find degeneration"<X(),N->Y(),N->Z()); if(Pdeg.Distance(Ptmp)<1.e-6) { DegNode = N; - //DegNode = const_cast(N); break; } } @@ -655,6 +824,7 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement* face PN.SetValue(nbp+1,PN(1)); FNodes[nbp] = FNodes[0]; + // find normal direction gp_Vec V1(PC,PN(nbp)); gp_Vec V2(PC,PN(1)); @@ -696,7 +866,6 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement* face } } - //cout<<" VNorm("< myPyramids; + const SMESHDS_SubMesh * aSubMeshDSFace; SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); - SMESH_MesherHelper helper(aMesh); - helper.IsQuadraticSubMesh(aShape); - helper.SetElementsOnShape( true ); + SMESH_MesherHelper helper1(aMesh); + helper1.IsQuadraticSubMesh(aShape); 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; + const int solidID = meshDS->ShapeToIndex( aShape ); - 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 ) aSubMeshDSFace = aProxyMesh->GetSubMesh( aShapeFace ); else aSubMeshDSFace = meshDS->MeshElements( aShapeFace ); + if ( !aSubMeshDSFace ) + continue; vector trias, quads; bool hasNewTrias = false; - if ( aSubMeshDSFace ) + const bool toCheckFaceInSolid = + aProxyMesh ? aProxyMesh->HasPrismsOnTwoSides( meshDS->MeshElements( aShapeFace )) : false; + if ( toCheckFaceInSolid && !dynamic_cast< const SMESH_ProxyMesh::SubMesh* >( aSubMeshDSFace )) + continue; // no room for pyramids as prisms are on both sides + { - bool isRev = false; - if ( helper.NbAncestors( aShapeFace, aMesh, aShape.ShapeType() ) > 1 ) - isRev = helper.IsReversedSubMesh( TopoDS::Face(aShapeFace) ); + bool isRevGlob = false; + SMESH_MesherHelper helper2( aMesh ); + PShapeIteratorPtr sIt = helper2.GetAncestors( aShapeFace, aMesh, aShape.ShapeType() ); + while ( const TopoDS_Shape * solid = sIt->next() ) + if ( !solid->IsSame( aShape )) + { + isRevGlob = helper2.IsReversedSubMesh( TopoDS::Face( aShapeFace )); + if ( toCheckFaceInSolid ) + helper2.IsQuadraticSubMesh( *solid ); + break; + } SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); while ( iteratorElem->more() ) // loop on elements on a geometrical face { const SMDS_MeshElement* face = iteratorElem->next(); + // preparation step to get face info - int stat = Preparation(face, PN, VN, FNodes, PC, VNorm); + int stat = Preparation( face, PN, VN, FNodes, PC, VNorm ); + + bool isRev = isRevGlob; + SMESH_MesherHelper* helper = &helper1; + if ( toCheckFaceInSolid && !isInSolid( FNodes, face->NbCornerNodes(), solidID )) + { + isRev = !isRevGlob; + helper = &helper2; + } + switch ( stat ) { case NOT_QUAD: @@ -773,10 +976,11 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, // degenerate face // add triangles to result map SMDS_MeshFace* NewFace; + helper->SetElementsOnShape( false ); if(!isRev) - NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] ); + NewFace = helper->AddFace( FNodes[0], FNodes[1], FNodes[2] ); else - NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] ); + NewFace = helper->AddFace( FNodes[0], FNodes[2], FNodes[1] ); storeTmpElement( NewFace ); trias.push_back ( NewFace ); quads.push_back( face ); @@ -809,32 +1013,28 @@ 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() ); + // create node at PCbest + helper->SetElementsOnShape( true ); + SMDS_MeshNode* NewNode = helper->AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() ); + + // create a pyramid + SMDS_MeshVolume* aPyram; + if ( isRev ) + aPyram = helper->AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode ); + else + aPyram = helper->AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); + myPyramids.push_back(aPyram); // add triangles to result map - for(i=0; i<4; i++) + helper->SetElementsOnShape( false ); + for ( i = 0; i < 4; i++ ) { - trias.push_back ( meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] )); + trias.push_back ( helper->AddFace( NewNode, FNodes[i], FNodes[i+1] )); storeTmpElement( trias.back() ); } - // create a pyramid - if ( isRev ) swap( FNodes[1], FNodes[3]); - SMDS_MeshVolume* aPyram = - helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); - myPyramids.push_back(aPyram); quads.push_back( face ); hasNewTrias = true; @@ -894,49 +1094,52 @@ 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); helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh()); - helper.SetElementsOnShape( true ); SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh(); if ( !myElemSearcher ) myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS ); - SMESH_ElementSearcher* searcher = const_cast(myElemSearcher); + 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); vector FNodes(5); TColgp_SequenceOfPnt aContour; - SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(/*idInceasingOrder=*/true); + SMDS_FaceIteratorPtr fIt = meshDS->facesIterator(); while( fIt->more()) { const SMDS_MeshElement* face = fIt->next(); @@ -945,11 +1148,11 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) gp_Pnt PC; gp_Vec VNorm; const SMDS_MeshElement* volumes[2]; - int what = Preparation(face, PN, VN, FNodes, PC, VNorm, volumes); + int what = Preparation( face, PN, VN, FNodes, PC, VNorm, volumes ); 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 +1174,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 ) @@ -1017,16 +1220,19 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) IsRev = true; } } + helper.SetElementsOnShape( false ); if(!IsRev) - NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] ); + NewFace = helper.AddFace( FNodes[0], FNodes[1], FNodes[2] ); else - NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] ); + NewFace = helper.AddFace( FNodes[0], FNodes[2], FNodes[1] ); storeTmpElement( NewFace ); prxSubMesh->AddElement( NewFace ); continue; } + // ----------------------------------- // Case of non-degenerated quadrangle + // ----------------------------------- // Find pyramid peak @@ -1048,6 +1254,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 +1262,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 a 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 +1297,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,23 +1333,19 @@ 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] ] ); - // create node for PCbest - SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() ); + if ( !LimitHeight( Papex, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/false )) + return false; + + // create node for Papex + helper.SetElementsOnShape( true ); + SMDS_MeshNode* NewNode = helper.AddNode( Papex.X(), Papex.Y(), Papex.Z() ); - // add triangles to result map - for(i=0; i<4; i++) { - SMDS_MeshFace* NewFace; - if(isRev) - NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ); - else - NewFace = meshDS->AddFace( NewNode, FNodes[i+1], FNodes[i] ); - storeTmpElement( NewFace ); - prxSubMesh->AddElement( NewFace ); - } // create a pyramid SMDS_MeshVolume* aPyram; if(isRev) @@ -1131,6 +1353,18 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) else aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode ); myPyramids.push_back(aPyram); + + // add triangles to result map + helper.SetElementsOnShape( false ); + for ( i = 0; i < 4; i++) { + SMDS_MeshFace* NewFace; + if(isRev) + NewFace = helper.AddFace( NewNode, FNodes[i], FNodes[i+1] ); + else + NewFace = helper.AddFace( NewNode, FNodes[i+1], FNodes[i] ); + storeTmpElement( NewFace ); + prxSubMesh->AddElement( NewFace ); + } } } // end loop on all faces @@ -1146,16 +1380,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 +1404,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 ); @@ -1185,75 +1424,82 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& while ( vIt->more() ) { const SMDS_MeshElement* PrmJ = vIt->next(); - if ( SMESH_MeshAlgos::GetCommonNodes( PrmI, PrmJ ).size() > 1 ) + if ( SMESH_MeshAlgos::NbCommonNodes( PrmI, PrmJ ) > 1 ) checkedPyrams.insert( PrmJ ); } } + // 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 +1518,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 );