X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_QuadToTriaAdaptor.cxx;h=f0d50dc0f321a45ae71992d52cff279c1dcab037;hp=a03ed6a82c19a97058955131261d1797f3db6daa;hb=dfcc0eb72b126380c8da3d725f90c9ee43f59827;hpb=b09212811c70415272e6b875949d71a1e2ae7025 diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx index a03ed6a82..f0d50dc0f 100644 --- a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx +++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2016 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 @@ -174,7 +174,7 @@ namespace PrmJ->GetNodeIndex( otherFaceNode ) >= 0 )) continue; // f is a base quadrangle - // check projections of face direction (baOFN) to triange normals (nI and nJ) + // 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 @@ -258,6 +258,7 @@ namespace meshDS->RemoveFreeNode( nodesToRemove[i], sm, /*fromGroups=*/false); } } + return; } //================================================================================ @@ -296,12 +297,13 @@ namespace 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 ); + SMESH_BadInputElements* badElems = + new SMESH_BadInputElements( mesh.GetMeshDS(),COMPERR_BAD_INPUT_MESH, msg, sm->GetAlgo() ); + badElems->add( face1 ); + badElems->add( face2 ); + err.reset( badElems ); } } - //throw SALOME_Exception( msg.c_str() ); return false; // == "algo fails" } @@ -318,7 +320,7 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr set & nodesToMove) { // cout << endl << "Merge " << PrmI->GetID() << " " << PrmJ->GetID() << " " - // << PrmI->GetNode(4) << PrmJ->GetNode(4) << endl; + // << 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 ); @@ -337,6 +339,9 @@ 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 @@ -354,6 +359,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 ); @@ -369,12 +379,27 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr const SMDS_MeshElement* elem = inverseElems[i]; 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; } //================================================================================ @@ -389,7 +414,7 @@ void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement* PrmI { 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); SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); @@ -413,6 +438,7 @@ void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement* PrmI for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm) MergeAdjacent( *prm, nodesToMove, true ); } + return; } //================================================================================ @@ -473,7 +499,7 @@ static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2, // 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) { const double EPSILON = 1e-6; @@ -485,9 +511,6 @@ static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, gp_XYZ vert1 = P2.XYZ(); gp_XYZ vert2 = P3.XYZ(); - /* calculate distance from vert0 to ray origin */ - gp_XYZ tvec = orig - vert0; - gp_XYZ edge1 = vert1 - vert0; gp_XYZ edge2 = vert2 - vert0; @@ -497,9 +520,13 @@ static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, /* if determinant is near zero, ray lies in plane of triangle */ double det = edge1 * pvec; - if (det > -EPSILON && det < EPSILON) + const double ANGL_EPSILON = 1e-12; + if ( det > -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) @@ -520,7 +547,27 @@ static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, Pint = orig + dir * t; - return ( t > 0. && t < segLen ); + 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; + } + return hasInt; } //======================================================================= @@ -665,7 +712,7 @@ bool StdMeshers_QuadToTriaAdaptor::LimitHeight (gp_Pnt& //================================================================================ /*! - * \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 @@ -704,7 +751,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 ) @@ -712,11 +760,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; } } @@ -747,6 +793,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)); @@ -788,7 +835,6 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement* face } } - //cout<<" VNorm("< itVec; @@ -876,10 +921,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 ); @@ -915,20 +961,25 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, if ( !LimitHeight( PCbest, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/true, aShape )) return false; } - // create node for PCbest + // 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; @@ -1012,24 +1063,28 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) 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(); @@ -1038,11 +1093,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 ) { @@ -1110,10 +1165,11 @@ 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; @@ -1143,6 +1199,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 @@ -1154,8 +1211,24 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) 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* > suspectFaces; + suspectFaces.clear(); searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces); for ( size_t iF = 0; iF < suspectFaces.size(); ++iF ) @@ -1169,7 +1242,7 @@ 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) ) + if( !volumes[isRev] && HasIntersection( farPnt[isRev], PC, intP, aContour )) { double d = PC.Distance( intP ); if ( d < dist2int[isRev] ) @@ -1215,25 +1288,28 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) return false; // create node for Papex + helper.SetElementsOnShape( true ); SMDS_MeshNode* NewNode = helper.AddNode( Papex.X(), Papex.Y(), Papex.Z() ); + // create a pyramid + SMDS_MeshVolume* aPyram; + if(isRev) + aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); + 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 = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ); + NewFace = helper.AddFace( NewNode, FNodes[i], FNodes[i+1] ); else - NewFace = meshDS->AddFace( NewNode, FNodes[i+1], FNodes[i] ); + NewFace = helper.AddFace( NewNode, FNodes[i+1], FNodes[i] ); storeTmpElement( NewFace ); prxSubMesh->AddElement( NewFace ); } - // create a pyramid - SMDS_MeshVolume* aPyram; - if(isRev) - aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode ); - else - aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode ); - myPyramids.push_back(aPyram); } } // end loop on all faces