#include "SMESH_subMeshEventListener.hxx"
#include "StdMeshers_FaceSide.hxx"
+#include <Adaptor3d_HSurface.hxx>
#include <BRepAdaptor_Curve2d.hxx>
#include <BRepAdaptor_Surface.hxx>
#include <BRepLProp_SLProps.hxx>
#include <TopoDS_Face.hxx>
#include <TopoDS_Vertex.hxx>
#include <gp_Ax1.hxx>
+#include <gp_Cone.hxx>
+#include <gp_Sphere.hxx>
#include <gp_Vec.hxx>
#include <gp_XY.hxx>
const double theSmoothThickToElemSizeRatio = 0.3;
// what part of thickness is allowed till intersection
- // defined by SALOME_TESTS/Grids/smesh/viscous_layers_00/A5
+ // (defined by SALOME_TESTS/Grids/smesh/viscous_layers_00/A5)
const double theThickToIntersection = 1.5;
bool needSmoothing( double cosin, double tgtThick, double elemSize )
void InvalidateStep( int curStep, bool restoreLength=false );
void ChooseSmooFunction(const set< TGeomID >& concaveVertices,
const TNode2Edge& n2eMap);
- bool Smooth(int& badNb, const int step, const bool isConcaveFace);
+ int Smooth(const int step, const bool isConcaveFace, const bool findBest);
bool SmoothOnEdge(Handle(Geom_Surface)& surface,
const TopoDS_Face& F,
SMESH_MesherHelper& helper);
SMESH_MesherHelper& helper,
bool& isOK,
bool shiftInside=false);
- gp_XYZ getWeigthedNormal( const SMDS_MeshNode* n,
- std::pair< TGeomID, gp_XYZ > fId2Normal[],
- int nbFaces );
+ bool getFaceNormalAtSingularity(const gp_XY& uv,
+ const TopoDS_Face& face,
+ SMESH_MesherHelper& helper,
+ gp_Dir& normal );
+ gp_XYZ getWeigthedNormal( const SMDS_MeshNode* n,
+ std::pair< TopoDS_Face, gp_XYZ > fId2Normal[],
+ int nbFaces );
bool findNeiborsOnEdge(const _LayerEdge* edge,
const SMDS_MeshNode*& n1,
const SMDS_MeshNode*& n2,
}
return done;
}
+ //================================================================================
+ /*!
+ * \brief Return direction of axis or revolution of a surface
+ */
+ //================================================================================
+
+ bool getRovolutionAxis( const Adaptor3d_Surface& surface,
+ gp_Dir & axis )
+ {
+ switch ( surface.GetType() ) {
+ case GeomAbs_Cone:
+ {
+ gp_Cone cone = surface.Cone();
+ axis = cone.Axis().Direction();
+ break;
+ }
+ case GeomAbs_Sphere:
+ {
+ gp_Sphere sphere = surface.Sphere();
+ axis = sphere.Position().Direction();
+ break;
+ }
+ case GeomAbs_SurfaceOfRevolution:
+ {
+ axis = surface.AxeOfRevolution().Direction();
+ break;
+ }
+ //case GeomAbs_SurfaceOfExtrusion:
+ case GeomAbs_OffsetSurface:
+ {
+ Handle(Adaptor3d_HSurface) base = surface.BasisSurface();
+ return getRovolutionAxis( base->Surface(), axis );
+ }
+ default: return false;
+ }
+ return true;
+ }
//--------------------------------------------------------------------------------
// DEBUG. Dump intermediate node positions into a python script
subIds = data._noShrinkShapes;
TopExp_Explorer exp( data._solid, TopAbs_FACE );
for ( ; exp.More(); exp.Next() )
+ {
+ SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
+ if ( ! data._ignoreFaceIds.count( fSubM->GetId() ))
{
- SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
- if ( ! data._ignoreFaceIds.count( fSubM->GetId() ))
- faceIds.insert( fSubM->GetId() );
+ faceIds.insert( fSubM->GetId() );
SMESH_subMeshIteratorPtr subIt = fSubM->getDependsOnIterator(/*includeSelf=*/true);
while ( subIt->more() )
subIds.insert( subIt->next()->GetId() );
}
+ }
// make a map to find new nodes on sub-shapes shared with other SOLID
map< TGeomID, TNode2Edge* >::iterator s2ne;
// Create temporary faces and _LayerEdge's
- dumpFunction(SMESH_Comment("makeLayers_")<<data._index);
+ dumpFunction(SMESH_Comment("makeLayers_")<<data._index);
data._stepSize = Precision::Infinite();
data._stepSizeNodes[0] = 0;
SMESH_MeshEditor editor(_mesh);
const SMDS_MeshNode* node = edge._nodes[0]; // source node
- SMDS_TypeOfPosition posType = node->GetPosition()->GetTypeOfPosition();
+ const SMDS_TypeOfPosition posType = node->GetPosition()->GetTypeOfPosition();
edge._len = 0;
edge._2neibors = 0;
edge._normal.SetCoord(0,0,0);
int totalNbFaces = 0;
+ TopoDS_Face F;
+ std::pair< TopoDS_Face, gp_XYZ > face2Norm[20];
gp_Vec geomNorm;
bool normOK = true;
+ // get geom FACEs the node lies on
+ {
+ set<TGeomID> faceIds;
+ if ( posType == SMDS_TOP_FACE )
+ {
+ faceIds.insert( node->getshapeId() );
+ }
+ else
+ {
+ SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() )
+ faceIds.insert( editor.FindShape(fIt->next()));
+ }
+ set<TGeomID>::iterator id = faceIds.begin();
+ for ( ; id != faceIds.end(); ++id )
+ {
+ const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id );
+ if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
+ continue;
+ F = TopoDS::Face( s );
+ face2Norm[ totalNbFaces ].first = F;
+ totalNbFaces++;
+ }
+ }
+
const TGeomID shapeInd = node->getshapeId();
map< TGeomID, TopoDS_Shape >::const_iterator s2s = data._shrinkShape2Shape.find( shapeInd );
const bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
+ // find _normal
if ( onShrinkShape ) // one of faces the node is on has no layers
{
TopoDS_Shape vertEdge = getMeshDS()->IndexToShape( s2s->first ); // vertex or edge
}
else // layers are on all faces of SOLID the node is on
{
- // find indices of geom faces the node lies on
- set<TGeomID> faceIds;
- if ( posType == SMDS_TOP_FACE )
+ int nbOkNorms = 0;
+ for ( int iF = 0; iF < totalNbFaces; ++iF )
{
- faceIds.insert( node->getshapeId() );
- }
- else
- {
- SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
- while ( fIt->more() )
- faceIds.insert( editor.FindShape(fIt->next()));
- }
-
- set<TGeomID>::iterator id = faceIds.begin();
- TopoDS_Face F;
- std::pair< TGeomID, gp_XYZ > id2Norm[20];
- for ( ; id != faceIds.end(); ++id )
- {
- const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id );
- if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
- continue;
- F = TopoDS::Face( s );
+ F = TopoDS::Face( face2Norm[ iF ].first );
geomNorm = getFaceNormal( node, F, helper, normOK );
if ( !normOK ) continue;
+ nbOkNorms++;
if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
geomNorm.Reverse();
- id2Norm[ totalNbFaces ].first = *id;
- id2Norm[ totalNbFaces ].second = geomNorm.XYZ();
- totalNbFaces++;
+ face2Norm[ iF ].second = geomNorm.XYZ();
edge._normal += geomNorm.XYZ();
}
- if ( totalNbFaces == 0 )
+ if ( nbOkNorms == 0 )
return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index);
- if ( normOK && edge._normal.Modulus() < 1e-3 && totalNbFaces > 1 )
+ if ( edge._normal.Modulus() < 1e-3 && nbOkNorms > 1 )
{
// opposite normals, re-get normals at shifted positions (IPAL 52426)
edge._normal.SetCoord( 0,0,0 );
- for ( int i = 0; i < totalNbFaces; ++i )
+ for ( int iF = 0; iF < totalNbFaces; ++iF )
{
- const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( id2Norm[i].first ));
+ const TopoDS_Face& F = face2Norm[iF].first;
geomNorm = getFaceNormal( node, F, helper, normOK, /*shiftInside=*/true );
if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
geomNorm.Reverse();
if ( normOK )
- id2Norm[ i ].second = geomNorm.XYZ();
- edge._normal += id2Norm[ i ].second;
+ face2Norm[ iF ].second = geomNorm.XYZ();
+ edge._normal += face2Norm[ iF ].second;
}
}
}
else
{
- edge._normal = getWeigthedNormal( node, id2Norm, totalNbFaces );
+ edge._normal = getWeigthedNormal( node, face2Norm, totalNbFaces );
}
+ }
- switch ( posType )
- {
- case SMDS_TOP_FACE:
- edge._cosin = 0; break;
-
- case SMDS_TOP_EDGE: {
- TopoDS_Edge E = TopoDS::Edge( helper.GetSubShapeByNode( node, getMeshDS()));
- gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK );
- double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
- edge._cosin = Cos( angle );
- //cout << "Cosin on EDGE " << edge._cosin << " node " << node->GetID() << endl;
- break;
- }
- case SMDS_TOP_VERTEX: {
- TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS()));
- gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK );
- double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
- edge._cosin = Cos( angle );
- if ( totalNbFaces > 2 || helper.IsSeamShape( node->getshapeId() ))
- for ( int iF = totalNbFaces-2; iF >=0; --iF )
- {
- F = TopoDS::Face( getMeshDS()->IndexToShape( id2Norm[ iF ].first ));
- inFaceDir = getFaceDir( F, V, node, helper, normOK );
- if ( normOK ) {
- double angle = inFaceDir.Angle( edge._normal );
- edge._cosin = Max( edge._cosin, Cos( angle ));
- }
+ // set _cosin
+ switch ( posType )
+ {
+ case SMDS_TOP_FACE: {
+ edge._cosin = 0;
+ break;
+ }
+ case SMDS_TOP_EDGE: {
+ TopoDS_Edge E = TopoDS::Edge( helper.GetSubShapeByNode( node, getMeshDS()));
+ gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK );
+ double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
+ edge._cosin = Cos( angle );
+ //cout << "Cosin on EDGE " << edge._cosin << " node " << node->GetID() << endl;
+ break;
+ }
+ case SMDS_TOP_VERTEX: {
+ TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS()));
+ gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK );
+ double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
+ edge._cosin = Cos( angle );
+ if ( totalNbFaces > 2 || helper.IsSeamShape( node->getshapeId() ))
+ for ( int iF = totalNbFaces-2; iF >=0; --iF )
+ {
+ F = face2Norm[ iF ].first;
+ inFaceDir = getFaceDir( F, V, node, helper, normOK );
+ if ( normOK ) {
+ double angle = inFaceDir.Angle( edge._normal );
+ edge._cosin = Max( edge._cosin, Cos( angle ));
}
- //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl;
- break;
- }
- default:
- return error(SMESH_Comment("Invalid shape position of node ")<<node, data._index);
- }
- } // case _sWOL.IsNull()
+ }
+ //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl;
+ break;
+ }
+ default:
+ return error(SMESH_Comment("Invalid shape position of node ")<<node, data._index);
+ }
double normSize = edge._normal.SquareModulus();
if ( normSize < numeric_limits<double>::min() )
isOK = false;
Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
+
+ if ( !shiftInside &&
+ helper.IsDegenShape( node->getshapeId() ) &&
+ getFaceNormalAtSingularity( uv, face, helper, normal ))
+ {
+ isOK = true;
+ return normal.XYZ();
+ }
+
int pointKind = GeomLib::NormEstim( surface, uv, 1e-5, normal );
enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE };
return normal.XYZ();
}
+//================================================================================
+/*!
+ * \brief Try to get normal at a singularity of a surface basing on it's nature
+ */
+//================================================================================
+
+bool _ViscousBuilder::getFaceNormalAtSingularity( const gp_XY& uv,
+ const TopoDS_Face& face,
+ SMESH_MesherHelper& helper,
+ gp_Dir& normal )
+{
+ BRepAdaptor_Surface surface( face );
+ gp_Dir axis;
+ if ( !getRovolutionAxis( surface, axis ))
+ return false;
+
+ double f,l, d, du, dv;
+ f = surface.FirstUParameter();
+ l = surface.LastUParameter();
+ d = ( uv.X() - f ) / ( l - f );
+ du = ( d < 0.5 ? +1. : -1 ) * 1e-5 * ( l - f );
+ f = surface.FirstVParameter();
+ l = surface.LastVParameter();
+ d = ( uv.Y() - f ) / ( l - f );
+ dv = ( d < 0.5 ? +1. : -1 ) * 1e-5 * ( l - f );
+
+ gp_Dir refDir;
+ gp_Pnt2d testUV = uv;
+ enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE };
+ double tol = 1e-5;
+ Handle(Geom_Surface) geomsurf = surface.Surface().Surface();
+ for ( int iLoop = 0; true ; ++iLoop )
+ {
+ testUV.SetCoord( testUV.X() + du, testUV.Y() + dv );
+ if ( GeomLib::NormEstim( geomsurf, testUV, tol, refDir ) == REGULAR )
+ break;
+ if ( iLoop > 20 )
+ return false;
+ tol /= 10.;
+ }
+
+ if ( axis * refDir < 0. )
+ axis.Reverse();
+
+ normal = axis;
+
+ return true;
+}
+
//================================================================================
/*!
* \brief Return a normal at a node weighted with angles taken by FACEs
*/
//================================================================================
-gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode* n,
- std::pair< TGeomID, gp_XYZ > fId2Normal[],
- int nbFaces )
+gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode* n,
+ std::pair< TopoDS_Face, gp_XYZ > fId2Normal[],
+ int nbFaces )
{
gp_XYZ resNorm(0,0,0);
TopoDS_Shape V = SMESH_MesherHelper::GetSubShapeByNode( n, getMeshDS() );
}
// exclude equal normals
- int nbUniqNorms = nbFaces;
+ //int nbUniqNorms = nbFaces;
for ( int i = 0; i < nbFaces; ++i )
for ( int j = i+1; j < nbFaces; ++j )
if ( fId2Normal[i].second.IsEqual( fId2Normal[j].second, 0.1 ))
{
fId2Normal[i].second.SetCoord( 0,0,0 );
- --nbUniqNorms;
+ //--nbUniqNorms;
break;
}
//if ( nbUniqNorms < 3 )
double angles[30];
for ( int i = 0; i < nbFaces; ++i )
{
- const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( fId2Normal[i].first ));
+ const TopoDS_Face& F = fId2Normal[i].first;
// look for two EDGEs shared by F and other FACEs within fId2Normal
TopoDS_Edge ee[2];
for ( int j = 0; j < nbFaces && !isSharedEdge; ++j )
{
if ( i == j ) continue;
- const TopoDS_Shape& otherF = getMeshDS()->IndexToShape( fId2Normal[j].first );
+ const TopoDS_Shape& otherF = fId2Normal[j].first;
isSharedEdge = SMESH_MesherHelper::IsSubShape( *E, otherF );
}
if ( !isSharedEdge )
return true; // no shapes needing smoothing
bool moved, improved;
+ vector< _LayerEdge* > badSmooEdges;
SMESH_MesherHelper helper(*_mesh);
Handle(Geom_Surface) surface;
const bool isConcaveFace = data._concaveFaces.count( sInd );
- int step = 0, stepLimit = 5, badNb = 0; moved = true;
- while (( ++step <= stepLimit && moved ) || improved )
+ int step = 0, stepLimit = 5, badNb = 0;
+ while (( ++step <= stepLimit ) || improved )
{
dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
<<"_InfStep"<<nbSteps<<"_"<<step); // debug
int oldBadNb = badNb;
- badNb = 0;
- moved = false;
- if ( step % 2 )
+ badSmooEdges.clear();
+
+ if ( step % 2 ) {
for ( int i = iBeg; i < iEnd; ++i ) // iterate forward
- moved |= data._edges[i]->Smooth( badNb, step, isConcaveFace );
- else
+ if ( data._edges[i]->Smooth( step, isConcaveFace, false ))
+ badSmooEdges.push_back( data._edges[i] );
+ }
+ else {
for ( int i = iEnd-1; i >= iBeg; --i ) // iterate backward
- moved |= data._edges[i]->Smooth( badNb, step, isConcaveFace );
+ if ( data._edges[i]->Smooth( step, isConcaveFace, false ))
+ badSmooEdges.push_back( data._edges[i] );
+ }
+ badNb = badSmooEdges.size();
improved = ( badNb < oldBadNb );
+ if ( !badSmooEdges.empty() && step >= stepLimit / 2 )
+ {
+ // look for the best smooth of _LayerEdge's neighboring badSmooEdges
+ vector<_Simplex> simplices;
+ for ( size_t i = 0; i < badSmooEdges.size(); ++i )
+ {
+ _LayerEdge* ledge = badSmooEdges[i];
+ _Simplex::GetSimplices( ledge->_nodes[0], simplices, data._ignoreFaceIds );
+ for ( size_t iS = 0; iS < simplices.size(); ++iS )
+ {
+ TNode2Edge::iterator n2e = data._n2eMap.find( simplices[iS]._nNext );
+ if ( n2e != data._n2eMap.end()) {
+ _LayerEdge* ledge2 = n2e->second;
+ if ( ledge2->_nodes[0]->getshapeId() == sInd )
+ ledge2->Smooth( step, isConcaveFace, /*findBest=*/true );
+ }
+ }
+ }
+ }
// issue 22576 -- no bad faces but still there are intersections to fix
// if ( improved && badNb == 0 )
// stepLimit = step + 3;
*/
//================================================================================
-bool _LayerEdge::Smooth(int& badNb, const int step, const bool isConcaveFace )
+int _LayerEdge::Smooth(const int step, const bool isConcaveFace, const bool findBest )
{
- bool moved = false;
if ( _simplices.size() < 2 )
- return moved; // _LayerEdge inflated along EDGE or FACE
+ return 0; // _LayerEdge inflated along EDGE or FACE
const gp_XYZ& curPos ( _pos.back() );
const gp_XYZ& prevPos( _pos[ _pos.size()-2 ]);
int nbBad = _simplices.size() - nbOkBefore;
// compute new position for the last _pos using different _funs
- gp_XYZ newPos, bestNewPos;
+ gp_XYZ newPos;
for ( int iFun = -1; iFun < theNbSmooFuns; ++iFun )
{
if ( iFun < 0 )
// get worse?
if ( nbOkAfter < nbOkBefore )
continue;
- if (( isConcaveFace ) &&
+ if (( isConcaveFace || findBest ) &&
( nbOkAfter == nbOkBefore ) &&
//( iFun > -1 || nbOkAfter < _simplices.size() ) &&
( minVolAfter <= minVolBefore ))
n->setXYZ( newPos.X(), newPos.Y(), newPos.Z());
_pos.back() = newPos;
- moved = true;
dumpMoveComm( n, _funNames[ iFun < 0 ? smooFunID() : iFun ]);
nbBad = _simplices.size() - nbOkAfter;
continue; // look for a better function
}
- break;
+ if ( !findBest )
+ break;
} // loop on smoothing functions
- badNb += nbBad;
- return moved;
+ return nbBad;
}
//================================================================================