c = new _Curvature;
c->_r = avgDist * avgDist / avgNormProj;
c->_k = avgDist * avgDist / c->_r / c->_r;
+ //c->_k = avgNormProj / c->_r;
c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive
c->_h2lenRatio = avgNormProj / ( avgDist + avgDist );
}
TopoDS_Shape _hypShape;
_MeshOfSolid* _proxyMesh;
set<TGeomID> _reversedFaceIds;
- set<TGeomID> _ignoreFaceIds;
+ set<TGeomID> _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDS
double _stepSize, _stepSizeCoeff;
const SMDS_MeshNode* _stepSizeNodes[2];
bool GetShapeEdges(const TGeomID shapeID, size_t& edgeEnd, int* iBeg=0, int* iEnd=0 ) const;
- void AddFacesToSmooth( const set< TGeomID >& faceIDs );
+ void AddShapesToSmooth( const set< TGeomID >& shapeIDs );
};
//--------------------------------------------------------------------------------
/*!
void limitStepSizeByCurvature( _SolidData& data );
void limitStepSize( _SolidData& data,
const SMDS_MeshElement* face,
- const double cosin);
+ const _LayerEdge* maxCosinEdge );
void limitStepSize( _SolidData& data, const double minSize);
bool inflate(_SolidData& data);
bool smoothAndCheck(_SolidData& data, const int nbSteps, double & distToIntersection);
if ( SMESH_Algo::isDegenerated( e )) continue;
TopExp::Vertices( e, VV[0], VV[1], /*CumOri=*/true );
if ( VV[1].IsSame( fromV )) {
+ nbEdges += edges[ 0 ].IsNull();
edges[ 0 ] = e;
- nbEdges++;
}
else if ( VV[0].IsSame( fromV )) {
+ nbEdges += edges[ 1 ].IsNull();
edges[ 1 ] = e;
- nbEdges++;
}
}
}
// get angle between the 2 edges
gp_Vec faceNormal;
- double angle = helper.GetAngle( edges[0], edges[1], faceFrw, &faceNormal );
+ double angle = helper.GetAngle( edges[0], edges[1], faceFrw, fromV, &faceNormal );
if ( Abs( angle ) < 5 * M_PI/180 )
{
dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] );
}
else if ( nbEdges == 1 )
{
- dir = getFaceDir( faceFrw, edges[0], node, helper, ok );
+ dir = getFaceDir( faceFrw, edges[ edges[0].IsNull() ], node, helper, ok );
if ( cosin ) *cosin = 1.;
}
else
while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 )))
iE2 = ( iE2 + 1 ) % nbEdges;
double angle = helper.GetAngle( wires[iW]->Edge( iE1 ),
- wires[iW]->Edge( iE2 ), F );
+ wires[iW]->Edge( iE2 ), F,
+ wires[iW]->FirstVertex( iE2 ));
if ( angle < -5. * M_PI / 180. )
return true;
}
}
//--------------------------------------------------------------------------------
// DEBUG. Dump intermediate node positions into a python script
+ // HOWTO use: run python commands written in a console to see
+ // construction steps of viscous layers
#ifdef __myDEBUG
ofstream* py;
+ int theNbFunc;
struct PyDump {
PyDump() {
const char* fname = "/tmp/viscous.py";
<< "smesh = smeshBuilder.New(salome.myStudy)" << endl
<< "meshSO = smesh.GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
<< "mesh = smesh.Mesh( meshSO.GetObject() )"<<endl;
+ theNbFunc = 0;
}
void Finish() {
- if (py)
- *py << "mesh.MakeGroup('Viscous Prisms',SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA)"<<endl;
+ if (py) {
+ *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Viscous Prisms',"
+ "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA))"<<endl;
+ *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Neg Volumes',"
+ "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_Volume3D,'<',0))"<<endl;
+ }
delete py; py=0;
}
- ~PyDump() { Finish(); }
+ ~PyDump() { Finish(); cout << "NB FUNCTIONS: " << theNbFunc << endl; }
};
#define dumpFunction(f) { _dumpFunction(f, __LINE__);}
#define dumpMove(n) { _dumpMove(n, __LINE__);}
#define dumpCmd(txt) { _dumpCmd(txt, __LINE__);}
void _dumpFunction(const string& fun, int ln)
- { if (py) *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl;}
+ { if (py) *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl; ++theNbFunc; }
void _dumpMove(const SMDS_MeshNode* n, int ln)
{ if (py) *py<< " mesh.MoveNode( "<<n->GetID()<< ", "<< n->X()
<< ", "<<n->Y()<<", "<< n->Z()<< ")\t\t # "<< ln <<endl; }
const SMDS_MeshElement* face = eIt->next();
newNodes.resize( face->NbCornerNodes() );
double faceMaxCosin = -1;
+ _LayerEdge* maxCosinEdge = 0;
for ( int i = 0 ; i < face->NbCornerNodes(); ++i )
{
const SMDS_MeshNode* n = face->GetNode(i);
return false;
}
dumpMove(edge->_nodes.back());
- if ( edge->_cosin > 0.01 )
+
+ if ( edge->_cosin > faceMaxCosin )
{
- if ( edge->_cosin > faceMaxCosin )
- faceMaxCosin = edge->_cosin;
+ faceMaxCosin = edge->_cosin;
+ maxCosinEdge = edge;
}
}
newNodes[ i ] = n2e->second->_nodes.back();
// compute inflation step size by min size of element on a convex surface
if ( faceMaxCosin > theMinSmoothCosin )
- limitStepSize( data, face, faceMaxCosin );
+ limitStepSize( data, face, maxCosinEdge );
} // loop on 2D elements on a FACE
} // loop on FACEs of a SOLID
void _ViscousBuilder::limitStepSize( _SolidData& data,
const SMDS_MeshElement* face,
- const double cosin)
+ const _LayerEdge* maxCosinEdge )
{
int iN = 0;
double minSize = 10 * data._stepSize;
for ( int i = 0; i < nbNodes; ++i )
{
const SMDS_MeshNode* nextN = face->GetNode( SMESH_MesherHelper::WrapIndex( i+1, nbNodes ));
- const SMDS_MeshNode* curN = face->GetNode( i );
+ const SMDS_MeshNode* curN = face->GetNode( i );
if ( nextN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ||
- curN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
+ curN-> GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
{
- double dist = SMESH_TNodeXYZ( face->GetNode(i)).Distance( nextN );
+ double dist = SMESH_TNodeXYZ( curN ).Distance( nextN );
if ( dist < minSize )
minSize = dist, iN = i;
}
}
- double newStep = 0.8 * minSize / cosin;
+ double newStep = 0.8 * minSize / maxCosinEdge->_lenFactor;
if ( newStep < data._stepSize )
{
data._stepSize = newStep;
- data._stepSizeCoeff = 0.8 / cosin;
+ data._stepSizeCoeff = 0.8 / maxCosinEdge->_lenFactor;
data._stepSizeNodes[0] = face->GetNode( iN );
data._stepSizeNodes[1] = face->GetNode( SMESH_MesherHelper::WrapIndex( iN+1, nbNodes ));
}
TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
vector<_LayerEdge*>& eV = edgesByGeom[ iV ];
if ( eV.empty() ) continue;
- double cosin = eV[0]->_cosin;
- bool badCosin =
- ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge));
- if ( badCosin )
- {
- gp_Vec dir1, dir2;
- if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE )
- dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ));
- else
- dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ),
- eV[0]->_nodes[0], helper, ok);
- dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
- double angle = dir1.Angle( dir2 );
- cosin = cos( angle );
- }
+ // double cosin = eV[0]->_cosin;
+ // bool badCosin =
+ // ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge));
+ // if ( badCosin )
+ // {
+ // gp_Vec dir1, dir2;
+ // if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE )
+ // dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ));
+ // else
+ // dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ),
+ // eV[0]->_nodes[0], helper, ok);
+ // dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
+ // double angle = dir1.Angle( dir2 );
+ // cosin = cos( angle );
+ // }
+ gp_Vec eDir = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
+ double angle = eDir.Angle( eV[0]->_normal );
+ double cosin = Cos( angle );
needSmooth = ( cosin > theMinSmoothCosin );
}
break;
const SMDS_MeshNode* node = edge._nodes[0]; // source node
SMDS_TypeOfPosition posType = node->GetPosition()->GetTypeOfPosition();
- edge._len = 0;
- edge._2neibors = 0;
+ edge._len = 0;
+ edge._2neibors = 0;
edge._curvature = 0;
// --------------------------
edge._normal.SetCoord(0,0,0);
int totalNbFaces = 0;
- gp_Pnt p;
- gp_Vec du, dv, geomNorm;
+ gp_Vec geomNorm;
bool normOK = true;
- TGeomID shapeInd = node->getshapeId();
+ const TGeomID shapeInd = node->getshapeId();
map< TGeomID, TopoDS_Shape >::const_iterator s2s = data._shrinkShape2Shape.find( shapeInd );
- bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
- TopoDS_Shape vertEdge;
+ const bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
if ( onShrinkShape ) // one of faces the node is on has no layers
{
- vertEdge = getMeshDS()->IndexToShape( s2s->first ); // vertex or edge
+ TopoDS_Shape vertEdge = getMeshDS()->IndexToShape( s2s->first ); // vertex or edge
if ( s2s->second.ShapeType() == TopAbs_EDGE )
{
// inflate from VERTEX along EDGE
isOK = false;
Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
- if ( GeomLib::NormEstim( surface, uv, 1e-10, normal ) < 3 )
+ int pointKind = GeomLib::NormEstim( surface, uv, 1e-5, normal );
+ enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE };
+ if ( pointKind < IMPOSSIBLE )
{
- normal;
+ if ( pointKind != REGULAR && !shiftInside )
+ {
+ gp_XYZ normShift = getFaceNormal( node, face, helper, isOK, /*shiftInside=*/true );
+ if ( normShift * normal.XYZ() < 0. )
+ normal = normShift;
+ }
isOK = true;
}
- else // hard singularity
+ else // hard singularity, to call with shiftInside=true ?
{
const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face );
}
else
{
- TopoDS_Vertex v10 = SMESH_MesherHelper::IthVertex( 1, ee[ 0 ]);
- TopoDS_Vertex v01 = SMESH_MesherHelper::IthVertex( 0, ee[ 1 ]);
- if ( !v10.IsSame( v01 ))
+ if ( !V.IsSame( SMESH_MesherHelper::IthVertex( 0, ee[ 1 ] )))
std::swap( ee[0], ee[1] );
}
- angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F );
+ angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F, TopoDS::Vertex( V ));
}
// compute a weighted normal
int oldBadNb = badNb;
badNb = 0;
moved = false;
- for ( int i = iBeg; i < iEnd; ++i )
- moved |= data._edges[i]->Smooth(badNb);
+ if ( step % 2 )
+ for ( int i = iBeg; i < iEnd; ++i )
+ moved |= data._edges[i]->Smooth(badNb);
+ else
+ for ( int i = iEnd-1; i >= iBeg; --i )
+ moved |= data._edges[i]->Smooth(badNb);
improved = ( badNb < oldBadNb );
- // issue 22576. no bad faces but still there are intersections to fix
+ // issue 22576 -- no bad faces but still there are intersections to fix
if ( improved && badNb == 0 )
stepLimit = step + 3;
int iLE = 0;
for ( size_t i = 0; i < data._edges.size(); ++i )
{
+ if ( !data._edges[i]->_sWOL.IsNull() )
+ continue;
if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace ))
return false;
if ( distToIntersection > dist )
*/
//================================================================================
-void _SolidData::AddFacesToSmooth( const set< TGeomID >& faceIDs )
+void _SolidData::AddShapesToSmooth( const set< TGeomID >& faceIDs )
{
// convert faceIDs to indices in _endEdgeOnShape
set< size_t > iEnds;
set< size_t >::iterator endsIt = iEnds.begin();
// "add" by move of _nbShapesToSmooth
- int nbFacesToAdd = faceIDs.size();
+ int nbFacesToAdd = iEnds.size();
while ( endsIt != iEnds.end() && *endsIt == _nbShapesToSmooth )
{
++endsIt;
for ( size_t i = 0; i < data._edges.size(); ++i )
{
_LayerEdge* edge = data._edges[i];
- if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
+ if (( !edge->IsOnEdge() ) &&
+ ( edge->_sWOL.IsNull() || edge->_sWOL.ShapeType() != TopAbs_FACE ))
+ continue;
if ( edge->FindIntersection( *searcher, dist, eps, &face ))
{
const _TmpMeshFaceOnEdge* f = (const _TmpMeshFaceOnEdge*) face;
{
dumpFunction(SMESH_Comment("updateNormals")<<data._index);
+ set< TGeomID > shapesToSmooth;
+
+ // vector to store new _normal and _cosin for each edge in edge2CloseEdge
+ vector< pair< _LayerEdge*, _LayerEdge > > edge2newEdge( edge2CloseEdge.size() );
+
TLEdge2LEdgeSet::iterator e2ee = edge2CloseEdge.begin();
- for ( ; e2ee != edge2CloseEdge.end(); ++e2ee )
+ for ( size_t iE = 0; e2ee != edge2CloseEdge.end(); ++e2ee, ++iE )
{
- _LayerEdge* edge1 = e2ee->first;
- _LayerEdge* edge2 = 0;
- set< _LayerEdge*, _LayerEdgeCmp >& ee = e2ee->second;
+ _LayerEdge* edge1 = e2ee->first;
+ _LayerEdge* edge2 = 0;
+ set< _LayerEdge*, _LayerEdgeCmp >& ee = e2ee->second;
+
+ edge2newEdge[ iE ].first = NULL;
// find EDGEs the edges reside
- TopoDS_Edge E1, E2;
- TopoDS_Shape S = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
- if ( S.ShapeType() != TopAbs_EDGE )
- continue; // TODO: find EDGE by VERTEX
- E1 = TopoDS::Edge( S );
+ // TopoDS_Edge E1, E2;
+ // TopoDS_Shape S = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
+ // if ( S.ShapeType() != TopAbs_EDGE )
+ // continue; // TODO: find EDGE by VERTEX
+ // E1 = TopoDS::Edge( S );
set< _LayerEdge*, _LayerEdgeCmp >::iterator eIt = ee.begin();
- while ( E2.IsNull() && eIt != ee.end())
+ for ( ; !edge2 && eIt != ee.end(); ++eIt )
{
- _LayerEdge* e2 = *eIt++;
- TopoDS_Shape S = helper.GetSubShapeByNode( e2->_nodes[0], getMeshDS() );
- if ( S.ShapeType() == TopAbs_EDGE )
- E2 = TopoDS::Edge( S ), edge2 = e2;
+ if ( edge1->_sWOL == (*eIt)->_sWOL )
+ edge2 = *eIt;
}
- if ( E2.IsNull() ) continue; // TODO: find EDGE by VERTEX
+ if ( !edge2 ) continue;
+
+ edge2newEdge[ iE ].first = edge1;
+ _LayerEdge& newEdge = edge2newEdge[ iE ].second;
+ // while ( E2.IsNull() && eIt != ee.end())
+ // {
+ // _LayerEdge* e2 = *eIt++;
+ // TopoDS_Shape S = helper.GetSubShapeByNode( e2->_nodes[0], getMeshDS() );
+ // if ( S.ShapeType() == TopAbs_EDGE )
+ // E2 = TopoDS::Edge( S ), edge2 = e2;
+ // }
+ // if ( E2.IsNull() ) continue; // TODO: find EDGE by VERTEX
// find 3 FACEs sharing 2 EDGEs
- TopoDS_Face FF1[2], FF2[2];
- PShapeIteratorPtr fIt = helper.GetAncestors(E1, *_mesh, TopAbs_FACE);
- while ( fIt->more() && FF1[1].IsNull())
+ // TopoDS_Face FF1[2], FF2[2];
+ // PShapeIteratorPtr fIt = helper.GetAncestors(E1, *_mesh, TopAbs_FACE);
+ // while ( fIt->more() && FF1[1].IsNull() )
+ // {
+ // const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
+ // if ( helper.IsSubShape( *F, data._solid))
+ // FF1[ FF1[0].IsNull() ? 0 : 1 ] = *F;
+ // }
+ // fIt = helper.GetAncestors(E2, *_mesh, TopAbs_FACE);
+ // while ( fIt->more() && FF2[1].IsNull())
+ // {
+ // const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
+ // if ( helper.IsSubShape( *F, data._solid))
+ // FF2[ FF2[0].IsNull() ? 0 : 1 ] = *F;
+ // }
+ // // exclude a FACE common to E1 and E2 (put it to FFn[1] )
+ // if ( FF1[0].IsSame( FF2[0]) || FF1[0].IsSame( FF2[1]))
+ // std::swap( FF1[0], FF1[1] );
+ // if ( FF2[0].IsSame( FF1[0]) )
+ // std::swap( FF2[0], FF2[1] );
+ // if ( FF1[0].IsNull() || FF2[0].IsNull() )
+ // continue;
+
+ // get a new normal for edge1
+ //bool ok;
+ gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal;
+ // if ( edge1->_cosin < 0 )
+ // dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized();
+ // if ( edge2->_cosin < 0 )
+ // dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized();
+
+ double cos1 = Abs( edge1->_cosin ), cos2 = Abs( edge2->_cosin );
+ double wgt1 = ( cos1 + 0.001 ) / ( cos1 + cos2 + 0.002 );
+ double wgt2 = ( cos2 + 0.001 ) / ( cos1 + cos2 + 0.002 );
+ newEdge._normal = ( wgt1 * dir1 + wgt2 * dir2 ).XYZ();
+ newEdge._normal.Normalize();
+
+ // cout << edge1->_nodes[0]->GetID() << " "
+ // << edge2->_nodes[0]->GetID() << " NORM: "
+ // << newEdge._normal.X() << ", " << newEdge._normal.Y() << ", " << newEdge._normal.Z() << endl;
+
+ // get new cosin
+ if ( cos1 < theMinSmoothCosin )
{
- const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
- if ( helper.IsSubShape( *F, data._solid))
- FF1[ FF1[0].IsNull() ? 0 : 1 ] = *F;
+ newEdge._cosin = edge2->_cosin;
}
- fIt = helper.GetAncestors(E2, *_mesh, TopAbs_FACE);
- while ( fIt->more() && FF2[1].IsNull())
+ else if ( cos2 > theMinSmoothCosin ) // both cos1 and cos2 > theMinSmoothCosin
{
- const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
- if ( helper.IsSubShape( *F, data._solid))
- FF2[ FF2[0].IsNull() ? 0 : 1 ] = *F;
+ // gp_Vec dirInFace;
+ // if ( edge1->_cosin < 0 )
+ // dirInFace = dir1;
+ // else
+ // dirInFace = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
+ // double angle = dirInFace.Angle( edge1->_normal ); // [0,PI]
+ // edge1->SetCosin( Cos( angle ));
+ //newEdge._cosin = 0; // ???????????
+ newEdge._cosin = ( wgt1 * cos1 + wgt2 * cos2 ) * edge1->_cosin / cos1;
}
- // exclude a FACE common to E1 and E2 (put it at [1] in FF* )
- if ( FF1[0].IsSame( FF2[0]) || FF1[0].IsSame( FF2[1]))
- std::swap( FF1[0], FF1[1] );
- if ( FF2[0].IsSame( FF1[0]) )
- std::swap( FF2[0], FF2[1] );
- if ( FF1[0].IsNull() || FF2[0].IsNull() )
- continue;
-
- // get a new normal for edge1
- bool ok;
- gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal;
- if ( edge1->_cosin < 0 )
- dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized();
- if ( edge2->_cosin < 0 )
- dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized();
- // gp_Vec dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
- // gp_Vec dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok2 );
- // double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
- // double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
- // gp_Vec newNorm = wgt1 * dir1 + wgt2 * dir2;
- // newNorm.Normalize();
-
- double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
- double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
- gp_Vec newNorm = wgt1 * dir1 + wgt2 * dir2;
- newNorm.Normalize();
-
- edge1->_normal = newNorm.XYZ();
-
- // update data of edge1 depending on _normal
- const SMDS_MeshNode *n1, *n2;
- n1 = edge1->_2neibors->_edges[0]->_nodes[0];
- n2 = edge1->_2neibors->_edges[1]->_nodes[0];
- edge1->SetDataByNeighbors( n1, n2, helper );
- gp_Vec dirInFace;
- if ( edge1->_cosin < 0 )
- dirInFace = dir1;
else
- dirInFace = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
- double angle = dir1.Angle( edge1->_normal ); // [0,PI]
- edge1->SetCosin( cos( angle ));
+ {
+ newEdge._cosin = edge1->_cosin;
+ }
- // limit data._stepSize
- if ( edge1->_cosin > theMinSmoothCosin )
+ // find shapes that need smoothing due to change of _normal
+ if ( edge1->_cosin < theMinSmoothCosin &&
+ newEdge._cosin > theMinSmoothCosin )
{
- SMDS_ElemIteratorPtr fIt = edge1->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
- while ( fIt->more() )
- limitStepSize( data, fIt->next(), edge1->_cosin );
+ if ( edge1->_sWOL.IsNull() )
+ {
+ SMDS_ElemIteratorPtr fIt = edge1->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() )
+ shapesToSmooth.insert( fIt->next()->getshapeId() );
+ //limitStepSize( data, fIt->next(), edge1->_cosin ); // too late
+ }
+ else // edge1 inflates along a FACE
+ {
+ TopoDS_Shape V = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
+ PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE );
+ while ( const TopoDS_Shape* E = eIt->next() )
+ {
+ if ( !helper.IsSubShape( *E, /*FACE=*/edge1->_sWOL ))
+ continue;
+ gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V ));
+ double angle = edgeDir.Angle( newEdge._normal ); // [0,PI]
+ if ( angle < M_PI / 2 )
+ shapesToSmooth.insert( getMeshDS()->ShapeToIndex( *E ));
+ }
+ }
}
- // set new XYZ of target node
+ }
+
+ data.AddShapesToSmooth( shapesToSmooth );
+
+ // Update data of edges depending on a new _normal
+
+ for ( size_t iE = 0; iE < edge2newEdge.size(); ++iE )
+ {
+ _LayerEdge* edge1 = edge2newEdge[ iE ].first;
+ _LayerEdge& newEdge = edge2newEdge[ iE ].second;
+ if ( !edge1 ) continue;
+
+ edge1->_normal = newEdge._normal;
+ edge1->SetCosin( newEdge._cosin );
edge1->InvalidateStep( 1 );
edge1->_len = 0;
edge1->SetNewLength( data._stepSize, helper );
- }
+ if ( edge1->IsOnEdge() )
+ {
+ const SMDS_MeshNode * n1 = edge1->_2neibors->_edges[0]->_nodes[0];
+ const SMDS_MeshNode * n2 = edge1->_2neibors->_edges[1]->_nodes[0];
+ edge1->SetDataByNeighbors( n1, n2, helper );
+ }
- // Update normals and other dependent data of not intersecting _LayerEdge's
- // neighboring the intersecting ones
+ // Update normals and other dependent data of not intersecting _LayerEdge's
+ // neighboring the intersecting ones
- for ( e2ee = edge2CloseEdge.begin(); e2ee != edge2CloseEdge.end(); ++e2ee )
- {
- _LayerEdge* edge1 = e2ee->first;
if ( !edge1->_2neibors )
continue;
for ( int j = 0; j < 2; ++j ) // loop on 2 neighbors
if ( edge2CloseEdge.count ( neighbor ))
continue; // j-th neighbor is also intersected
_LayerEdge* prevEdge = edge1;
- const int nbSteps = 6;
+ const int nbSteps = 10;
for ( int step = nbSteps; step; --step ) // step from edge1 in j-th direction
{
if ( !neighbor->_2neibors )
for ( ; iBeg < iEnd; ++iBeg )
avgNormal += data._edges[ iBeg ]->_normal;
}
- avgNormal.Normalize();
+ double normSize = avgNormal.SquareModulus();
+ if ( normSize < 1e-200 )
+ {
+ debugMsg( "updateNormalsOfConvexFaces(): zero avgNormal" );
+ return false;
+ }
+ avgNormal /= Sqrt( normSize );
// compute new _LayerEdge::_cosin on EDGEs
double avgCosin = 0;
}
}
}
- data.AddFacesToSmooth( adjFacesToSmooth );
+ data.AddShapesToSmooth( adjFacesToSmooth );
dumpFunctionEnd();
newNormal += norm;
double sz = newNormal.Modulus();
- if ( Abs ( sz ) < 1e-200 )
+ if ( sz < 1e-200 )
break;
newNormal /= sz;
return true;
bool segmentIntersected = false;
distance = Precision::Infinite();
int iFace = -1; // intersected face
- for ( size_t j = 0 ; j < suspectFaces.size() && !segmentIntersected; ++j )
+ for ( size_t j = 0 ; j < suspectFaces.size() /*&& !segmentIntersected*/; ++j )
{
const SMDS_MeshElement* face = suspectFaces[j];
if ( face->GetNodeIndex( _nodes.back() ) >= 0 ||
{
const SMDS_MeshNode* tria[3];
tria[0] = *nIt++;
- tria[1] = *nIt++;;
+ tria[1] = *nIt++;
for ( int n2 = 2; n2 < nbNodes && !intFound; ++n2 )
{
tria[2] = *nIt++;
}
if ( intFound )
{
- if ( dist < segLen*(1.01) && dist > -(_len-segLen) )
+ if ( dist < segLen*(1.01) && dist > -(_len*_lenFactor-segLen) )
segmentIntersected = true;
if ( distance > dist )
distance = dist, iFace = j;
/* calculate distance from vert0 to ray origin */
gp_XYZ tvec = orig - vert0;
- if ( tvec * dir > EPSILON )
+ //if ( tvec * dir > EPSILON )
// intersected face is at back side of the temporary face this _LayerEdge belongs to
- return false;
+ //return false;
gp_XYZ edge1 = vert1 - vert0;
gp_XYZ edge2 = vert2 - vert0;
double det = edge1 * pvec;
if (det > -EPSILON && det < EPSILON)
- return 0;
+ return false;
double inv_det = 1.0 / det;
/* calculate U parameter and test bounds */
double u = ( tvec * pvec ) * inv_det;
- if (u < 0.0 || u > 1.0)
- return 0;
+ //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) * inv_det;
- if ( v < 0.0 || u + v > 1.0 )
- return 0;
+ //if ( v < 0.0 || u + v > 1.0 )
+ if ( v < -EPSILON || u + v > 1.0 + EPSILON)
+ return false;
/* calculate t, ray intersects triangle */
t = (edge2 * qvec) * inv_det;
- return true;
+ //return true;
+ return t > 0.;
}
//================================================================================
newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev );
newPos /= _simplices.size();
+ const gp_XYZ& curPos ( _pos.back() );
+ const gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
if ( _curvature )
- newPos += _normal * _curvature->lenDelta( _len );
-
- gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
+ {
+ double delta = _curvature->lenDelta( _len );
+ if ( delta > 0 )
+ newPos += _normal * delta;
+ else
+ {
+ double segLen = _normal * ( newPos - prevPos.XYZ() );
+ if ( segLen + delta > 0 )
+ newPos += _normal * delta;
+ }
+ // double segLenChange = _normal * ( curPos - newPos );
+ // newPos += 0.5 * _normal * segLenChange;
+ }
// count quality metrics (orientation) of tetras around _tgtNode
int nbOkBefore = 0;
- SMESH_TNodeXYZ tgtXYZ( _nodes.back() );
for ( size_t i = 0; i < _simplices.size(); ++i )
- nbOkBefore += _simplices[i].IsForward( _nodes[0], &tgtXYZ );
+ nbOkBefore += _simplices[i].IsForward( _nodes[0], &curPos );
int nbOkAfter = 0;
for ( size_t i = 0; i < _simplices.size(); ++i )
dumpChangeNodes( f );
// Replace source nodes by target nodes in mesh faces to shrink
+ dumpFunction(SMESH_Comment("replNodesOnFace")<<f2sd->first); // debug
const SMDS_MeshNode* nodes[20];
for ( size_t i = 0; i < lEdges.size(); ++i )
{
nodes[iN] = ( n == srcNode ? tgtNode : n );
}
helper.GetMeshDS()->ChangeElementNodes( f, nodes, f->NbNodes() );
+ dumpChangeNodes( f );
}
}
// Create _SmoothNode's on face F
vector< _SmoothNode > nodesToSmooth( smoothNodes.size() );
{
+ dumpFunction(SMESH_Comment("fixUVOnFace")<<f2sd->first); // debug
const bool sortSimplices = isConcaveFace;
for ( size_t i = 0; i < smoothNodes.size(); ++i )
{
}
}
}
+ // TODO: check effect of this additional smooth
+ // additional laplacian smooth to increase allowed shrink step
+ // for ( int st = 1; st; --st )
+ // {
+ // dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
+ // for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
+ // {
+ // nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
+ // _SmoothNode::LAPLACIAN,/*set3D=*/false);
+ // }
+ // }
} // while ( shrinked )
// No wrongly shaped faces remain; final smooth. Set node XYZ.
gp_Vec2d uvDir( srcUV, tgtUV );
double uvLen = uvDir.Magnitude();
uvDir /= uvLen;
- edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0);
+ edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0 );
edge._len = uvLen;
edge._pos.resize(1);
const double kSafe = Max( 0.5, 1. - 0.1 * _simplices.size() );
// Select shrinking step such that not to make faces with wrong orientation.
- double stepSize = uvLen;
+ double stepSize = 1e100;
for ( size_t i = 0; i < _simplices.size(); ++i )
{
// find intersection of 2 lines: curUV-tgtUV and that connecting simplex nodes
stepSize = Min( step, stepSize );
}
gp_Pnt2d newUV;
- if ( uvLen - stepSize < _len / 200. )
+ if ( uvLen <= stepSize )
{
newUV = tgtUV;
_pos.clear();