ForkMaster/smesh
1
0
Fork 0
mirror of https://git.salome-platform.org/gitpub/modules/smesh.git synced 2025-01-26 13:20:34 +05:00
Code Issues Packages Projects Releases Wiki Activity

use again latest code from master but with fix to use the point with the maximum shift as pyramid apex

Browse Source
This commit is contained in:
mbs 2024-10-21 21:42:37 +01:00
parent 34a6e0e20f
commit 89ce13b860
1 changed files with 77 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

199
src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx
View File

@ -60,22 +60,22 @@ enum EQuadNature { NOT_QUAD, QUAD, DEGEN_QUAD, PYRAM_APEX = 4, TRIA_APEX = 0 };
// std-like iterator used to get coordinates of nodes of mesh element // std-like iterator used to get coordinates of nodes of mesh element
typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator; typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
// //================================================================================ //================================================================================
// /*! /*!
// * \brief Return ID of pyramid base face, for debug * \brief Return ID of pyramid base face, for debug
// */ */
// //================================================================================ //================================================================================
//
// int getFaceID(const SMDS_MeshElement* pyram) int getFaceID(const SMDS_MeshElement* pyram)
// { {
// if ( pyram ) if ( pyram )
// if ( const SMDS_MeshElement* f = SMDS_Mesh::FindFace( pyram->GetNode(0), if ( const SMDS_MeshElement* f = SMDS_Mesh::FindFace( pyram->GetNode(0),
// pyram->GetNode(1), pyram->GetNode(1),
// pyram->GetNode(2), pyram->GetNode(2),
// pyram->GetNode(3))) pyram->GetNode(3)))
// return f->GetID(); return f->GetID();
// return -1; return -1;
// } }
namespace namespace
{ {
@ -537,49 +537,49 @@ void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement* PrmI
*/ */
//================================================================================ //================================================================================
// bool StdMeshers_QuadToTriaAdaptor::DecreaseHeightDifference( const SMDS_MeshElement* thePyram, bool StdMeshers_QuadToTriaAdaptor::DecreaseHeightDifference( const SMDS_MeshElement* thePyram,
// const double theH2 ) const double theH2 )
// { {
// const double allowedFactor2 = 2. * 2.; const double allowedFactor2 = 2. * 2.;
//
// bool modif = false; bool modif = false;
// myNodes[0] = thePyram->GetNode( 3 ); myNodes[0] = thePyram->GetNode( 3 );
// for ( int i = 0; i < 4; ++i ) for ( int i = 0; i < 4; ++i )
// { {
// myNodes[1] = thePyram->GetNode( i ); myNodes[1] = thePyram->GetNode( i );
// SMDS_Mesh::GetElementsByNodes( myNodes, myAdjPyrams, SMDSAbs_Volume ); SMDS_Mesh::GetElementsByNodes( myNodes, myAdjPyrams, SMDSAbs_Volume );
// myNodes[0] = myNodes[1]; myNodes[0] = myNodes[1];
//
// for ( const SMDS_MeshElement* pyramAdj : myAdjPyrams ) for ( const SMDS_MeshElement* pyramAdj : myAdjPyrams )
// { {
// if ( pyramAdj == thePyram ) if ( pyramAdj == thePyram )
// continue; continue;
// if ( !myPyramHeight2.IsBound( pyramAdj )) if ( !myPyramHeight2.IsBound( pyramAdj ))
// continue; continue;
// double h2Adj = Abs( myPyramHeight2( pyramAdj )); double h2Adj = Abs( myPyramHeight2( pyramAdj ));
// double h2 = Abs( theH2 ); double h2 = Abs( theH2 );
// if ( h2Adj > h2 ) if ( h2Adj > h2 )
// { {
// if ( h2 * allowedFactor2 < h2Adj ) if ( h2 * allowedFactor2 < h2Adj )
// { {
// // bind negative value to allow finding pyramids whose height must change // bind negative value to allow finding pyramids whose height must change
// myPyramHeight2.Bind( pyramAdj, - h2 * allowedFactor2 ); myPyramHeight2.Bind( pyramAdj, - h2 * allowedFactor2 );
// modif = true; modif = true;
// } }
// } }
// else else
// { {
// if ( h2Adj * allowedFactor2 < h2 ) if ( h2Adj * allowedFactor2 < h2 )
// { {
// // bind negative value to allow finding pyramids whose height must change // bind negative value to allow finding pyramids whose height must change
// myPyramHeight2.Bind( thePyram, - h2Adj * allowedFactor2 ); myPyramHeight2.Bind( thePyram, - h2Adj * allowedFactor2 );
// modif = true; modif = true;
// } }
// } }
// } }
// } }
// return modif; return modif;
// } }
//================================================================================ //================================================================================
@ -614,11 +614,11 @@ StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor()
//======================================================================= //=======================================================================
static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2, static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2,
const gp_Pnt& PC, const gp_Vec& V/*, const gp_Pnt& PC, const gp_Vec& V,
double & shift*/) double & shift)
{ {
gp_Pnt Pbest = PC; gp_Pnt Pbest = PC;
//shift = 0; shift = 0;
const double a2 = P1.SquareDistance(P2); const double a2 = P1.SquareDistance(P2);
const double b2 = P1.SquareDistance(PC); const double b2 = P1.SquareDistance(PC);
const double c2 = P2.SquareDistance(PC); const double c2 = P2.SquareDistance(PC);
@ -629,7 +629,7 @@ static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2,
const double Vsize = V.Magnitude(); const double Vsize = V.Magnitude();
if ( fabs( Vsize ) > std::numeric_limits<double>::min() ) if ( fabs( Vsize ) > std::numeric_limits<double>::min() )
{ {
const double shift = sqrt( a2 + (b2-c2)*(b2-c2)/16/a2 - (b2+c2)/2 ); shift = sqrt( a2 + (b2-c2)*(b2-c2)/16/a2 - (b2+c2)/2 );
Pbest.ChangeCoord() += shift * V.XYZ() / Vsize; Pbest.ChangeCoord() += shift * V.XYZ() / Vsize;
} }
} }
@ -1103,9 +1103,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh,
case QUAD: case QUAD:
{ {
if(!isRev) VNorm.Reverse(); if(!isRev) VNorm.Reverse();
#if 0 // Find the point with the maximum shift
//double xc = 0., yc = 0., zc = 0.; double h, hMax = 0.;
double h, hMin = Precision::Infinite();
gp_Pnt PCbest = PC; gp_Pnt PCbest = PC;
int i = 1; int i = 1;
for(; i<=4; i++) { for(; i<=4; i++) {
@ -1114,28 +1113,12 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh,
Pbest = FindBestPoint(PN(i), PN(i+1), PC, VN(i).Reversed(), h); Pbest = FindBestPoint(PN(i), PN(i+1), PC, VN(i).Reversed(), h);
else else
Pbest = FindBestPoint(PN(i), PN(i+1), PC, VN(i), h); Pbest = FindBestPoint(PN(i), PN(i+1), PC, VN(i), h);
if ( 0 < h && h < hMin ) if ( h > hMax )
{ {
PCbest = Pbest; PCbest = Pbest;
hMin = h; hMax = h;
} }
} }
//gp_Pnt PCbest(xc/4., yc/4., zc/4.);
#else
double xc = 0., yc = 0., zc = 0.;
int i = 1;
for(; i<=4; i++) {
gp_Pnt Pbest;
if(!isRev)
Pbest = FindBestPoint(PN(i), PN(i+1), PC, VN(i).Reversed());
else
Pbest = FindBestPoint(PN(i), PN(i+1), PC, VN(i));
xc += Pbest.X();
yc += Pbest.Y();
zc += Pbest.Z();
}
gp_Pnt PCbest(xc/4., yc/4., zc/4.);
#endif
// check PCbest // check PCbest
double height = PCbest.Distance(PC); double height = PCbest.Distance(PC);
@ -1161,7 +1144,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh,
myPyramids.push_back(aPyram); myPyramids.push_back(aPyram);
//cout << "F" << face->GetID() << " - V" << aPyram->GetID() << endl; //cout << "F" << face->GetID() << " - V" << aPyram->GetID() << endl;
// myPyramHeight2.Bind( aPyram, PCbest.SquareDistance( PC )); myPyramHeight2.Bind( aPyram, PCbest.SquareDistance( PC ));
// add triangles to result map // add triangles to result map
helper->SetElementsOnShape( false ); helper->SetElementsOnShape( false );
@ -1370,29 +1353,19 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
// ----------------------------------- // -----------------------------------
// Find pyramid peak // Find pyramid peak
#if 0
gp_XYZ PCbest = PC.XYZ();//(0., 0., 0.); // pyramid peak // Find the point with the maximum shift
double h, hMin = Precision::Infinite(); double h, hMax = 0.;
gp_Pnt PCbest = PC;
int i = 1; int i = 1;
for ( ; i <= 4; i++ ) { for(; i<=4; i++) {
gp_Pnt Pbest = FindBestPoint(PN(i), PN(i+1), PC, VN(i), h); gp_Pnt Pbest = FindBestPoint(PN(i), PN(i+1), PC, VN(i), h);
if ( 0 < h && h < hMin ) if ( h > hMax )
{ {
PCbest = Pbest.XYZ(); PCbest = Pbest;
h = hMin; hMax = h;
} }
//PCbest += Pbest.XYZ();
} }
//PCbest /= 4;
#else
gp_XYZ PCbest(0., 0., 0.); // pyramid peak
int i = 1;
for ( ; i <= 4; i++ ) {
gp_Pnt Pbest = FindBestPoint(PN(i), PN(i+1), PC, VN(i));
PCbest += Pbest.XYZ();
}
PCbest /= 4;
#endif
double height = PC.Distance(PCbest); // pyramid height to precise double height = PC.Distance(PCbest); // pyramid height to precise
if ( height < 1.e-6 ) { if ( height < 1.e-6 ) {
@ -1550,13 +1523,6 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
// check adjacent pyramids // check adjacent pyramids
#if 1
for ( i = 0; i < myPyramids.size(); ++i )
{
const SMDS_MeshElement* PrmI = myPyramids[i];
MergeAdjacent( PrmI, nodesToMove );
}
#else
// Fix adjacent pyramids whose heights differ too much // Fix adjacent pyramids whose heights differ too much
myNodes.resize(2); myNodes.resize(2);
@ -1583,7 +1549,6 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
gp_Pnt newApex = gp_Pnt( PC ).Translated( h * V.Normalized() ); gp_Pnt newApex = gp_Pnt( PC ).Translated( h * V.Normalized() );
meshDS->MoveNode( Papex.Node(), newApex.X(), newApex.Y(), newApex.Z() ); meshDS->MoveNode( Papex.Node(), newApex.X(), newApex.Y(), newApex.Z() );
} }
#endif
// iterate on all new pyramids // iterate on all new pyramids
vector< const SMDS_MeshElement* > suspectPyrams; vector< const SMDS_MeshElement* > suspectPyrams;
@ -1671,7 +1636,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
if ( nbc == 4 ) if ( nbc == 4 )
continue; // pyrams have a common base face continue; // pyrams have a common base face
if (nbc > 0) //( nbc > 1 ) if ( nbc > 1 )
{ {
// Merge the two pyramids and others already merged with them // Merge the two pyramids and others already merged with them
MergePiramids( PrmI, PrmJ, nodesToMove ); MergePiramids( PrmI, PrmJ, nodesToMove );
@ -1693,9 +1658,6 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
double ang2 = fabs(VN2.Angle(VI2)); double ang2 = fabs(VN2.Angle(VI2));
double coef1 = 0.5 - (( ang1 < M_PI/3. ) ? cos(ang1)*0.25 : 0 ); double coef1 = 0.5 - (( ang1 < M_PI/3. ) ? cos(ang1)*0.25 : 0 );
double coef2 = 0.5 - (( ang2 < M_PI/3. ) ? cos(ang2)*0.25 : 0 ); // cos(ang2) ? double coef2 = 0.5 - (( ang2 < M_PI/3. ) ? cos(ang2)*0.25 : 0 ); // cos(ang2) ?
// double coef2 = 0.5;
// if(ang2<PI/3)
// coef2 -= cos(ang1)*0.25;
VN1.Scale(coef1); VN1.Scale(coef1);
VN2.Scale(coef2); VN2.Scale(coef2);
@ -1707,11 +1669,6 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
nodesToMove.insert( aNode2 ); nodesToMove.insert( aNode2 );
//cout << "Limit H F" << getFaceID( PrmI ) << " - F" << getFaceID( PrmJ ) << endl; //cout << "Limit H F" << getFaceID( PrmI ) << " - F" << getFaceID( PrmJ ) << endl;
} }
#if 1
// fix intersections that can appear after apex movement
MergeAdjacent( PrmI, nodesToMove );
MergeAdjacent( PrmJ, nodesToMove );
#endif
apexI = PrmI->GetNode( PYRAM_APEX ); // apexI can be removed by merge apexI = PrmI->GetNode( PYRAM_APEX ); // apexI can be removed by merge
@ -1721,14 +1678,12 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
} // loop on all pyramids } // loop on all pyramids
#if 0
//smIdType nbNodes = aMesh.NbNodes(); //smIdType nbNodes = aMesh.NbNodes();
for ( i = 0; i < myPyramids.size(); ++i ) for ( i = 0; i < myPyramids.size(); ++i )
{ {
const SMDS_MeshElement* PrmI = myPyramids[i]; const SMDS_MeshElement* PrmI = myPyramids[i];
MergeAdjacent( PrmI, nodesToMove ); MergeAdjacent( PrmI, nodesToMove );
} }
#endif
if( !nodesToMove.empty() && !meshDS->IsEmbeddedMode() ) if( !nodesToMove.empty() && !meshDS->IsEmbeddedMode() )
{ {