From 86b8303fcf1c4a377b8bc2b99472a438249a4a9f Mon Sep 17 00:00:00 2001 From: eap Date: Wed, 30 Nov 2011 13:49:21 +0000 Subject: [PATCH] 0021422: EDF 1963 SMESH: Viscous layer algorithm fails in some cases (cylindre_partition.py) fix shrinking on cancave FACEs --- src/StdMeshers/StdMeshers_ViscousLayers.cxx | 330 ++++++++++++++++++-- 1 file changed, 299 insertions(+), 31 deletions(-) diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index ee32b9c47..4952afab6 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -32,16 +32,18 @@ #include "SMESHDS_Hypothesis.hxx" #include "SMESH_Algo.hxx" #include "SMESH_ComputeError.hxx" +#include "SMESH_ControlsDef.hxx" #include "SMESH_Gen.hxx" #include "SMESH_Group.hxx" #include "SMESH_Mesh.hxx" #include "SMESH_MesherHelper.hxx" +#include "SMESH_ProxyMesh.hxx" #include "SMESH_subMesh.hxx" #include "SMESH_subMeshEventListener.hxx" -#include "SMESH_ProxyMesh.hxx" #include "utilities.h" +#include #include #include #include @@ -56,6 +58,8 @@ #include #include #include +#include +#include #include #include #include @@ -67,7 +71,6 @@ #include #include #include -#include #include #include @@ -236,6 +239,10 @@ namespace VISCOUS double d = v1 ^ v2; return d*refSign > 1e-100; } + bool IsNeighbour(const _Simplex& other) const + { + return _nPrev == other._nNext || _nNext == other._nPrev; + } }; //-------------------------------------------------------------------------------- /*! @@ -412,6 +419,7 @@ namespace VISCOUS Handle(Geom_Surface)& surface, SMESH_MesherHelper& helper, const double refSign, + bool isCentroidal, bool set3D); }; //-------------------------------------------------------------------------------- @@ -445,7 +453,8 @@ namespace VISCOUS _SolidData& data); void getSimplices( const SMDS_MeshNode* node, vector<_Simplex>& simplices, const set& ingnoreShapes, - const _SolidData* dataToCheckOri = 0); + const _SolidData* dataToCheckOri = 0, + const bool toSort = false); bool sortEdges( _SolidData& data, vector< vector<_LayerEdge*> >& edgesByGeom); void limitStepSize( _SolidData& data, @@ -466,6 +475,7 @@ namespace VISCOUS bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F, SMESH_MesherHelper& helper, const SMESHDS_SubMesh* faceSubMesh ); + void fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper); bool addBoundaryElements(); bool error( const string& text, int solidID=-1 ); @@ -727,6 +737,67 @@ namespace } return dir; } + //================================================================================ + /*! + * \brief Returns true if a FACE is bound by a concave EDGE + */ + //================================================================================ + + bool isConcave( const TopoDS_Face& F, SMESH_MesherHelper& helper ) + { + gp_Vec2d drv1, drv2; + gp_Pnt2d p; + TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE ); + for ( ; eExp.More(); eExp.Next() ) + { + const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() ); + if ( BRep_Tool::Degenerated( E )) continue; + // check if 2D curve is concave + BRepAdaptor_Curve2d curve( E, F ); + const int nbIntervals = curve.NbIntervals( GeomAbs_C2 ); + TColStd_Array1OfReal intervals(1, nbIntervals + 1 ); + curve.Intervals( intervals, GeomAbs_C2 ); + bool isConvex = true; + for ( int i = 1; i <= nbIntervals && isConvex; ++i ) + { + double u1 = intervals( i ); + double u2 = intervals( i+1 ); + curve.D2( 0.5*( u1+u2 ), p, drv1, drv2 ); + double cross = drv2 ^ drv1; + if ( E.Orientation() == TopAbs_REVERSED ) + cross = -cross; + isConvex = ( cross < 1e-9 ); + } + // check if concavity is strong enough to care about it + //const double maxAngle = 5 * Standard_PI180; + if ( !isConvex ) + { + //cout << "Concave FACE " << helper.GetMeshDS()->ShapeToIndex( F ) << endl; + return true; + // map< double, const SMDS_MeshNode* > u2nodes; + // if ( !SMESH_Algo::GetSortedNodesOnEdge( helper.GetMeshDS(), E, + // /*ignoreMedium=*/true, u2nodes)) + // continue; + // map< double, const SMDS_MeshNode* >::iterator u2n = u2nodes.begin(); + // gp_Pnt2d uvPrev = helper.GetNodeUV( F, u2n->second ); + // double uPrev = u2n->first; + // for ( ++u2n; u2n != u2nodes.end(); ++u2n ) + // { + // gp_Pnt2d uv = helper.GetNodeUV( F, u2n->second ); + // gp_Vec2d segmentDir( uvPrev, uv ); + // curve.D1( uPrev, p, drv1 ); + // try { + // if ( fabs( segmentDir.Angle( drv1 )) > maxAngle ) + // return true; + // } + // catch ( ... ) {} + // uvPrev = uv; + // uPrev = u2n->first; + // } + } + } + return false; + } //-------------------------------------------------------------------------------- // DEBUG. Dump intermediate node positions into a python script #ifdef __myDEBUG @@ -738,31 +809,38 @@ namespace py = new ofstream(fname); *py << "from smesh import *" << endl << "meshSO = GetCurrentStudy().FindObjectID('0:1:2:3')" << endl - << "mesh = Mesh( meshSO.GetObject()._narrow( SMESH.SMESH_Mesh ))"< nodesToSmooth( smoothNodes.size() ); { @@ -3456,7 +3561,7 @@ bool _ViscousBuilder::shrink() const SMDS_MeshNode* n = smoothNodes[i]; nodesToSmooth[ i ]._node = n; // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices - getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes ); + getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, isConcaveFace ); // fix up incorrect uv of nodes on the FACE helper.GetNodeUV( F, n, 0, &isOkUV); dumpMove( n ); @@ -3521,7 +3626,8 @@ bool _ViscousBuilder::shrink() moved = false; for ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) { - moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,/*set3D=*/false ); + moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign, + /*isCentroidal=*/isConcaveFace,/*set3D=*/false ); } if ( badNb < oldBadNb ) nbNoImpSteps = 0; @@ -3532,9 +3638,6 @@ bool _ViscousBuilder::shrink() } if ( badNb > 0 ) return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first ); - - if ( !shrinked ) - break; } // No wrongly shaped faces remain; final smooth. Set node XYZ. // First, find out a needed quality of smoothing (high for quadrangles only) @@ -3554,11 +3657,14 @@ bool _ViscousBuilder::shrink() highQuality = ( *nbNodesSet.begin() == 4 ); } } + if ( !highQuality && isConcaveFace ) + fixBadFaces( F, helper ); // fix narrow faces by swaping diagonals for ( int st = highQuality ? 10 : 3; st; --st ) { dumpFunction(SMESH_Comment("shrinkFace")<first<<"_st"<<++smooStep); // debug for ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) - nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,/*set3D=*/st==1 ); + nodesToSmooth[i].Smooth( badNb,surface,helper,refSign, + /*isCentroidal=*/isConcaveFace,/*set3D=*/st==1 ); dumpFunctionEnd(); } // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh @@ -3750,6 +3856,124 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, // return true; } +//================================================================================ +/*! + * \brief Try to fix triangles with high aspect ratio by swaping diagonals + */ +//================================================================================ + +void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper) +{ + SMESH::Controls::AspectRatio qualifier; + SMESH::Controls::TSequenceOfXYZ points(3), points1(3), points2(3); + const double maxAspectRatio = 4.; + + // find bad triangles + + vector< const SMDS_MeshElement* > badTrias; + vector< double > badAspects; + SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( F ); + SMDS_ElemIteratorPtr fIt = sm->GetElements(); + while ( fIt->more() ) + { + const SMDS_MeshElement * f = fIt->next(); + if ( f->NbCornerNodes() != 3 ) continue; + for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = SMESH_TNodeXYZ( f->GetNode(iP)); + double aspect = qualifier.GetValue( points ); + if ( aspect > maxAspectRatio ) + { + badTrias.push_back( f ); + badAspects.push_back( aspect ); + } + } + if ( badTrias.empty() ) + return; + + // find couples of faces to swap diagonal + + typedef pair < const SMDS_MeshElement* , const SMDS_MeshElement* > T2Trias; + vector< T2Trias > triaCouples; + + TIDSortedElemSet involvedFaces, emptySet; + for ( size_t iTia = 0; iTia < badTrias.size(); ++iTia ) + { + T2Trias trias [3]; + double aspRatio [3]; + int i1, i2, i3; + + involvedFaces.insert( badTrias[iTia] ); + for ( int iP = 0; iP < 3; ++iP ) + points(iP+1) = SMESH_TNodeXYZ( badTrias[iTia]->GetNode(iP)); + + // find triangles adjacent to badTrias[iTia] with better aspect ratio after diag-swaping + int bestCouple = -1; + for ( int iSide = 0; iSide < 3; ++iSide ) + { + const SMDS_MeshNode* n1 = badTrias[iTia]->GetNode( iSide ); + const SMDS_MeshNode* n2 = badTrias[iTia]->GetNode(( iSide+1 ) % 3 ); + trias [iSide].first = badTrias[iTia]; + trias [iSide].second = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, involvedFaces, + & i1, & i2 ); + if ( ! trias[iSide].second || trias[iSide].second->NbCornerNodes() != 3 ) + continue; + + // aspect ratio of an adjacent tria + for ( int iP = 0; iP < 3; ++iP ) + points2(iP+1) = SMESH_TNodeXYZ( trias[iSide].second->GetNode(iP)); + double aspectInit = qualifier.GetValue( points2 ); + + // arrange nodes as after diag-swaping + if ( helper.WrapIndex( i1+1, 3 ) == i2 ) + i3 = helper.WrapIndex( i1-1, 3 ); + else + i3 = helper.WrapIndex( i1+1, 3 ); + points1 = points; + points1( 1+ iSide ) = points2( 1+ i3 ); + points2( 1+ i2 ) = points1( 1+ ( iSide+2 ) % 3 ); + + // aspect ratio after diag-swaping + aspRatio[ iSide ] = qualifier.GetValue( points1 ) + qualifier.GetValue( points2 ); + if ( aspRatio[ iSide ] > aspectInit + badAspects[ iTia ] ) + continue; + + if ( bestCouple < 0 || aspRatio[ bestCouple ] > aspRatio[ iSide ] ) + bestCouple = iSide; + } + + if ( bestCouple >= 0 ) + { + triaCouples.push_back( trias[bestCouple] ); + involvedFaces.insert ( trias[bestCouple].second ); + } + else + { + involvedFaces.erase( badTrias[iTia] ); + } + } + if ( triaCouples.empty() ) + return; + + // swap diagonals + + SMESH_MeshEditor editor( helper.GetMesh() ); + dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")< uv( _simplices.size() ); + for ( size_t i = 0; i < _simplices.size(); ++i ) + uv[i] = helper.GetNodeUV( face, _simplices[i]._nPrev, _node ); + // compute new UV for the node gp_XY newPos (0,0); - for ( unsigned i = 0; i < _simplices.size(); ++i ) - newPos += helper.GetNodeUV( face, _simplices[i]._nPrev, _node ); - newPos /= _simplices.size(); + if ( isCentroidal && _simplices.size() > 3 ) + { + // average centers of diagonals wieghted with their reciprocal lengths + if ( _simplices.size() == 4 ) + { + double w1 = 1. / ( uv[2]-uv[0] ).SquareModulus(); + double w2 = 1. / ( uv[3]-uv[1] ).SquareModulus(); + newPos = ( w1 * ( uv[2]+uv[0] ) + w2 * ( uv[3]+uv[1] )) / ( w1+w2 ) / 2; + } + else + { + double sumWeight = 0; + int nb = _simplices.size() == 4 ? 2 : _simplices.size(); + for ( int i = 0; i < nb; ++i ) + { + int iFrom = i + 2; + int iTo = i + _simplices.size() - 1; + for ( int j = iFrom; j < iTo; ++j ) + { + int i2 = SMESH_MesherHelper::WrapIndex( j, _simplices.size() ); + double w = 1. / ( uv[i]-uv[i2] ).SquareModulus(); + sumWeight += w; + newPos += w * ( uv[i]+uv[i2] ); + } + } + newPos /= 2 * sumWeight; + } + } + else + { + // Laplacian smooth + isCentroidal = false; + for ( size_t i = 0; i < _simplices.size(); ++i ) + newPos += uv[i]; + newPos /= _simplices.size(); + } // count quality metrics (orientation) of triangles around the node int nbOkBefore = 0; @@ -3873,7 +4136,12 @@ bool _SmoothNode::Smooth(int& badNb, nbOkAfter += _simplices[i].IsForward( newPos, _node, face, helper, refSign ); if ( nbOkAfter < nbOkBefore ) + { + // if ( isCentroidal ) + // return Smooth( badNb, surface, helper, refSign, !isCentroidal, set3D ); + badNb += _simplices.size() - nbOkBefore; return false; + } SMDS_FacePosition* pos = static_cast( _node->GetPosition() ); pos->SetUParameter( newPos.X() ); -- 2.30.2