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 );
}
double lenDelta(double len) const { return _k * ( _r + len ); }
double lenDeltaByDist(double dist) const { return dist * _h2lenRatio; }
};
- struct _LayerEdge;
- //--------------------------------------------------------------------------------
- /*!
- * Structure used to smooth a _LayerEdge (master) based on an EDGE.
- */
- struct _2NearEdges
- {
- // target nodes of 2 neighbour _LayerEdge's based on the same EDGE
- const SMDS_MeshNode* _nodes[2];
- // vectors from source nodes of 2 _LayerEdge's to the source node of master _LayerEdge
- //gp_XYZ _vec[2];
- double _wgt[2]; // weights of _nodes
- _LayerEdge* _edges[2];
-
- // normal to plane passing through _LayerEdge._normal and tangent of EDGE
- gp_XYZ* _plnNorm;
-
- _2NearEdges() { _nodes[0]=_nodes[1]=0; _plnNorm = 0; }
- void reverse() {
- std::swap( _nodes[0], _nodes[1] );
- std::swap( _wgt [0], _wgt [1] );
- std::swap( _edges[0], _edges[1] );
- }
- };
+ struct _2NearEdges;
//--------------------------------------------------------------------------------
/*!
* \brief Edge normal to surface, connecting a node on solid surface (_nodes[0])
double& dist,
const double& epsilon) const;
gp_Ax1 LastSegment(double& segLen) const;
+ gp_XY LastUV( const TopoDS_Face& F ) const;
bool IsOnEdge() const { return _2neibors; }
gp_XYZ Copy( _LayerEdge& other, SMESH_MesherHelper& helper );
void SetCosin( double cosin );
return cmpNodes ? ( e1->_nodes[0]->GetID() < e2->_nodes[0]->GetID()) : ( e1 < e2 );
}
};
+ struct _LayerEdge;
+ //--------------------------------------------------------------------------------
+ /*!
+ * Structure used to smooth a _LayerEdge based on an EDGE.
+ */
+ struct _2NearEdges
+ {
+ double _wgt [2]; // weights of _nodes
+ _LayerEdge* _edges[2];
+
+ // normal to plane passing through _LayerEdge._normal and tangent of EDGE
+ gp_XYZ* _plnNorm;
+
+ _2NearEdges() { _edges[0]=_edges[1]=0; _plnNorm = 0; }
+ const SMDS_MeshNode* tgtNode(bool is2nd) {
+ return _edges[is2nd] ? _edges[is2nd]->_nodes.back() : 0;
+ }
+ const SMDS_MeshNode* srcNode(bool is2nd) {
+ return _edges[is2nd] ? _edges[is2nd]->_nodes[0] : 0;
+ }
+ void reverse() {
+ std::swap( _wgt [0], _wgt [1] );
+ std::swap( _edges[0], _edges[1] );
+ }
+ };
//--------------------------------------------------------------------------------
/*!
* \brief Convex FACE whose radius of curvature is less than the thickness of
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];
- TNode2Edge _n2eMap;
+ TNode2Edge _n2eMap; // nodes and _LayerEdge's based on them
+
// map to find _n2eMap of another _SolidData by a shrink shape shared by two _SolidData's
map< TGeomID, TNode2Edge* > _s2neMap;
// edges of _n2eMap. We keep same data in two containers because
// Convex FACEs whose radius of curvature is less than the thickness of layers
map< TGeomID, _ConvexFace > _convexFaces;
- // FACE's WOL, srink on which is forbiden due to algo on the adjacent SOLID
- set< TGeomID > _noShrinkFaces;
+ // shapes (EDGEs and VERTEXes) srink from which is forbiden due to collisions with
+ // the adjacent SOLID
+ set< TGeomID > _noShrinkShapes;
// <EDGE to smooth on> to <it's curve> -- for analytic smooth
map< TGeomID,Handle(Geom_Curve)> _edge2curve;
double _epsilon; // precision for SegTriaInter()
- int _index; // for debug
+ TGeomID _index; // SOLID id, for debug
_SolidData(const TopoDS_Shape& s=TopoDS_Shape(),
const StdMeshers_ViscousLayers* h=0,
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);
bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F,
SMESH_MesherHelper& helper,
const SMESHDS_SubMesh* faceSubMesh );
+ void restoreNoShrink( _LayerEdge& edge ) const;
void fixBadFaces(const TopoDS_Face& F,
SMESH_MesherHelper& helper,
const bool is2D,
bool addBoundaryElements();
bool error( const string& text, int solidID=-1 );
- SMESHDS_Mesh* getMeshDS() { return _mesh->GetMeshDS(); }
+ SMESHDS_Mesh* getMeshDS() const { return _mesh->GetMeshDS(); }
// debug
void makeGroupOfLE();
virtual SMDSAbs_ElementType GetType() const { return SMDSAbs_Face; }
virtual vtkIdType GetVtkType() const { return -1; }
virtual SMDSAbs_EntityType GetEntityType() const { return SMDSEntity_Last; }
- virtual SMDSAbs_GeometryType GetGeomType() const { return SMDSGeom_TRIANGLE; }
+ virtual SMDSAbs_GeometryType GetGeomType() const
+ { return _nn.size() == 3 ? SMDSGeom_TRIANGLE : SMDSGeom_QUADRANGLE; }
virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType) const
{ return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));}
};
return dir.XYZ();
}
//--------------------------------------------------------------------------------
+ gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
+ const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok,
+ double* cosin=0);
+ //--------------------------------------------------------------------------------
gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE,
const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok)
{
+ double f,l;
+ Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l );
+ if ( c.IsNull() )
+ {
+ TopoDS_Vertex v = helper.IthVertex( 0, fromE );
+ return getFaceDir( F, v, node, helper, ok );
+ }
gp_XY uv = helper.GetNodeUV( F, node, 0, &ok );
Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
gp_Pnt p; gp_Vec du, dv, norm;
surface->D1( uv.X(),uv.Y(), p, du,dv );
norm = du ^ dv;
- double f,l;
- Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l );
double u = helper.GetNodeU( fromE, node, 0, &ok );
c->D1( u, p, du );
TopAbs_Orientation o = helper.GetSubShapeOri( F.Oriented(TopAbs_FORWARD), fromE);
//--------------------------------------------------------------------------------
gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
const SMDS_MeshNode* node, SMESH_MesherHelper& helper,
- bool& ok, double* cosin=0)
+ bool& ok, double* cosin)
{
TopoDS_Face faceFrw = F;
faceFrw.Orientation( TopAbs_FORWARD );
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; }
bool _ViscousBuilder::error(const string& text, int solidId )
{
+ const string prefix = string("Viscous layers builder: ");
_error->myName = COMPERR_ALGO_FAILED;
- _error->myComment = string("Viscous layers builder: ") + text;
+ _error->myComment = prefix + text;
if ( _mesh )
{
SMESH_subMesh* sm = _mesh->GetSubMeshContaining( solidId );
if ( !sm && !_sdVec.empty() )
- sm = _mesh->GetSubMeshContaining( _sdVec[0]._index );
+ sm = _mesh->GetSubMeshContaining( solidId = _sdVec[0]._index );
if ( sm && sm->GetSubShape().ShapeType() == TopAbs_SOLID )
{
SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
if ( smError && smError->myAlgo )
_error->myAlgo = smError->myAlgo;
smError = _error;
+ sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
+ }
+ // set KO to all solids
+ for ( size_t i = 0; i < _sdVec.size(); ++i )
+ {
+ if ( _sdVec[i]._index == solidId )
+ continue;
+ sm = _mesh->GetSubMesh( _sdVec[i]._solid );
+ if ( !sm->IsEmpty() )
+ continue;
+ SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+ if ( !smError || smError->IsOK() )
+ {
+ smError = SMESH_ComputeError::New( COMPERR_ALGO_FAILED, prefix + "failed");
+ sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
+ }
}
}
makeGroupOfLE(); // debug
continue; // nothing interesting
TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ];
// check presence of layers on fWOL within an adjacent SOLID
+ bool collision = false;
PShapeIteratorPtr sIt = helper.GetAncestors( fWOL, *_mesh, TopAbs_SOLID );
while ( const TopoDS_Shape* solid = sIt->next() )
if ( !solid->IsSame( _sdVec[i]._solid ))
int iFace = getMeshDS()->ShapeToIndex( fWOL );
if ( iSolid > 0 && !_sdVec[ iSolid-1 ]._ignoreFaceIds.count( iFace ))
{
- _sdVec[i]._noShrinkFaces.insert( iFace );
- fWOL.Nullify();
+ //_sdVec[i]._noShrinkShapes.insert( iFace );
+ //fWOL.Nullify();
+ collision = true;
}
}
// add edge to maps
{
TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
_sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
+ if ( collision )
+ {
+ // _shrinkShape2Shape will be used to temporary inflate _LayerEdge's based
+ // on the edge but shrink won't be performed
+ _sdVec[i]._noShrinkShapes.insert( edgeInd );
+ }
}
}
}
set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D");
for ( size_t i = 0; i < _sdVec.size(); ++i )
{
- TopTools_MapOfShape noShrinkVertices;
map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
{
const TopoDS_Shape& fWOL = e2f->second;
- TGeomID edgeID = e2f->first;
+ const TGeomID edgeID = e2f->first;
bool notShrinkFace = false;
PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
- while ( soIt->more())
+ while ( soIt->more() )
{
const TopoDS_Shape* solid = soIt->next();
if ( _sdVec[i]._solid.IsSame( *solid )) continue;
SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid );
if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue;
notShrinkFace = true;
- for ( size_t j = 0; j < _sdVec.size(); ++j )
+ size_t iSolid = 0;
+ for ( ; iSolid < _sdVec.size(); ++iSolid )
{
- if ( _sdVec[j]._solid.IsSame( *solid ) )
- if ( _sdVec[j]._shrinkShape2Shape.count( edgeID ))
+ if ( _sdVec[iSolid]._solid.IsSame( *solid ) ) {
+ if ( _sdVec[iSolid]._shrinkShape2Shape.count( edgeID ))
notShrinkFace = false;
+ break;
+ }
}
- }
- if ( notShrinkFace )
- {
- _sdVec[i]._noShrinkFaces.insert( getMeshDS()->ShapeToIndex( fWOL ));
- for ( TopExp_Explorer vExp( fWOL, TopAbs_VERTEX ); vExp.More(); vExp.Next() )
- noShrinkVertices.Add( vExp.Current() );
- }
- }
- // erase from _shrinkShape2Shape all srink EDGE's of a SOLID connected
- // to the found not shrinked fWOL's
- e2f = _sdVec[i]._shrinkShape2Shape.begin();
- for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); )
- {
- TGeomID edgeID = e2f->first;
- TopoDS_Vertex VV[2];
- TopExp::Vertices( TopoDS::Edge( getMeshDS()->IndexToShape( edgeID )),VV[0],VV[1]);
- if ( noShrinkVertices.Contains( VV[0] ) || noShrinkVertices.Contains( VV[1] ))
- {
- _sdVec[i]._noShrinkFaces.insert( getMeshDS()->ShapeToIndex( e2f->second ));
- _sdVec[i]._shrinkShape2Shape.erase( e2f++ );
- }
- else
- {
- e2f++;
- }
- }
- }
+ if ( notShrinkFace )
+ {
+ _sdVec[i]._noShrinkShapes.insert( edgeID );
+
+ // add VERTEXes of the edge in _noShrinkShapes
+ TopoDS_Shape edge = getMeshDS()->IndexToShape( edgeID );
+ for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() )
+ _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( vIt.Value() ));
+
+ // check if there is a collision with to-shrink-from EDGEs in iSolid
+ if ( iSolid == _sdVec.size() )
+ continue; // no VL in the solid
+ shapes.Clear();
+ TopExp::MapShapes( fWOL, TopAbs_EDGE, shapes);
+ for ( int iE = 1; iE <= shapes.Extent(); ++iE )
+ {
+ const TopoDS_Edge& E = TopoDS::Edge( shapes( iE ));
+ const TGeomID eID = getMeshDS()->ShapeToIndex( E );
+ if ( eID == edgeID ||
+ !_sdVec[iSolid]._shrinkShape2Shape.count( eID ) ||
+ _sdVec[i]._noShrinkShapes.count( eID ))
+ continue;
+ for ( int is1st = 0; is1st < 2; ++is1st )
+ {
+ TopoDS_Vertex V = helper.IthVertex( is1st, E );
+ if ( _sdVec[i]._noShrinkShapes.count( getMeshDS()->ShapeToIndex( V ) ))
+ {
+ // _sdVec[i]._noShrinkShapes.insert( eID );
+ // V = helper.IthVertex( !is1st, E );
+ // _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( V ));
+ //iE = 0; // re-start the loop on EDGEs of fWOL
+ return error("No way to make a conformal mesh with "
+ "the given set of faces with layers", _sdVec[i]._index);
+ }
+ }
+ }
+ }
+
+ } // while ( soIt->more() )
+ } // loop on _sdVec[i]._shrinkShape2Shape
+ } // loop on _sdVec to fill in _SolidData::_noShrinkShapes
// Find the SHAPE along which to inflate _LayerEdge based on VERTEX
{
totalNbFaces++;
const int fID = getMeshDS()->ShapeToIndex( *f );
- if ( _sdVec[i]._ignoreFaceIds.count ( fID ) &&
- !_sdVec[i]._noShrinkFaces.count( fID ))
+ if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&&
+ !_sdVec[i]._noShrinkShapes.count( fID )*/)
facesWOL.push_back( *f );
}
}
{
// get all sub-shapes to make layers on
set<TGeomID> subIds, faceIds;
- subIds = data._noShrinkFaces;
+ 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( getMeshDS()->ShapeToIndex( exp.Current() )))
+ if ( ! data._ignoreFaceIds.count( fSubM->GetId() ))
faceIds.insert( fSubM->GetId() );
- SMESH_subMeshIteratorPtr subIt =
- fSubM->getDependsOnIterator(/*includeSelf=*/true, /*complexShapeFirst=*/false);
+ SMESH_subMeshIteratorPtr subIt = fSubM->getDependsOnIterator(/*includeSelf=*/true);
while ( subIt->more() )
subIds.insert( subIt->next()->GetId() );
}
SMESH_MesherHelper helper( *_mesh );
helper.SetSubShape( data._solid );
- helper.SetElementsOnShape(true);
+ helper.SetElementsOnShape( true );
vector< const SMDS_MeshNode*> newNodes; // of a mesh face
TNode2Edge::iterator n2e2;
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);
_LayerEdge* edge = new _LayerEdge();
n2e->second = edge;
edge->_nodes.push_back( n );
- const int shapeID = n->getshapeId();
+ const int shapeID = n->getshapeId();
+ const bool noShrink = data._noShrinkShapes.count( shapeID );
edgesByGeom[ shapeID ].push_back( edge );
SMESH_TNodeXYZ xyz( n );
// set edge data or find already refined _LayerEdge and get data from it
- if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
- ( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() &&
- ( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end())
+ if (( !noShrink ) &&
+ ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE ) &&
+ (( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() ) &&
+ (( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end() ))
{
_LayerEdge* foundEdge = (*n2e2).second;
gp_XYZ lastPos = edge->Copy( *foundEdge, helper );
}
else
{
- edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ));
+ if ( !noShrink )
+ {
+ edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ));
+ }
if ( !setEdgeData( *edge, subIds, helper, data ))
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
// limit data._stepSize depending on surface curvature and fill data._convexFaces
limitStepSizeByCurvature( data ); // !!! it must be before node substitution in _Simplex
- // Set target nodes into _Simplex and _2NearEdges of _LayerEdge's
+ // Set target nodes into _Simplex and _LayerEdge's to _2NearEdges
TNode2Edge::iterator n2e;
+ const SMDS_MeshNode* nn[2];
for ( size_t i = 0; i < data._edges.size(); ++i )
{
- if ( data._edges[i]->IsOnEdge())
+ if ( data._edges[i]->IsOnEdge() )
+ {
+ // get neighbor nodes
+ bool hasData = ( data._edges[i]->_2neibors->_edges[0] );
+ if ( hasData ) // _LayerEdge is a copy of another one
+ {
+ nn[0] = data._edges[i]->_2neibors->srcNode(0);
+ nn[1] = data._edges[i]->_2neibors->srcNode(1);
+ }
+ else if ( !findNeiborsOnEdge( data._edges[i], nn[0],nn[1], data ))
+ {
+ return false;
+ }
+ // set neighbor _LayerEdge's
for ( int j = 0; j < 2; ++j )
{
- if ( data._edges[i]->_nodes.back()->NbInverseElements(SMDSAbs_Volume) > 0 )
- break; // _LayerEdge is shared by two _SolidData's
- const SMDS_MeshNode* & n = data._edges[i]->_2neibors->_nodes[j];
- if (( n2e = data._n2eMap.find( n )) == data._n2eMap.end() )
+ if (( n2e = data._n2eMap.find( nn[j] )) == data._n2eMap.end() )
return error("_LayerEdge not found by src node", data._index);
- n = (*n2e).second->_nodes.back();
data._edges[i]->_2neibors->_edges[j] = n2e->second;
}
- //else
+ if ( !hasData )
+ data._edges[i]->SetDataByNeighbors( nn[0], nn[1], helper);
+ }
+
for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j )
{
_Simplex& s = data._edges[i]->_simplices[j];
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 ));
}
for ( ; face.More(); face.Next() )
{
const TopoDS_Face& F = TopoDS::Face( face.Current() );
+ SMESH_subMesh * sm = _mesh->GetSubMesh( F );
+ const TGeomID faceID = sm->GetId();
+ if ( data._ignoreFaceIds.count( faceID )) continue;
+
BRepAdaptor_Surface surface( F, false );
surfProp.SetSurface( surface );
int iBeg, iEnd;
_ConvexFace cnvFace;
- SMESH_subMesh * sm = _mesh->GetSubMesh( F );
- const TGeomID faceID = sm->GetId();
- if ( data._ignoreFaceIds.count( faceID )) continue;
const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. );
SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
while ( smIt->more() )
continue;
// check concavity and curvature and limit data._stepSize
int nbLEdges = iEnd - iBeg;
- int step = Max( 1, nbLEdges / nbTestPnt );
- for ( ; iBeg < iEnd; iBeg += step )
+ int iStep = Max( 1, nbLEdges / nbTestPnt );
+ for ( ; iBeg < iEnd; iBeg += iStep )
{
gp_XY uv = helper.GetNodeUV( F, data._edges[ iBeg ]->_nodes[0] );
surfProp.SetParameters( uv.X(), uv.Y() );
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
if ( edge._sWOL.ShapeType() == TopAbs_EDGE )
{
double u = helper.GetNodeU( TopoDS::Edge( edge._sWOL ), node, 0, &normOK );
- edge._pos.push_back( gp_XYZ( u, 0, 0));
- getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( edge._sWOL ), u );
+ edge._pos.push_back( gp_XYZ( u, 0, 0 ));
+ if ( edge._nodes.size() > 1 )
+ getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( edge._sWOL ), u );
}
else // TopAbs_FACE
{
gp_XY uv = helper.GetNodeUV( TopoDS::Face( edge._sWOL ), node, 0, &normOK );
edge._pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
- getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( edge._sWOL ), uv.X(), uv.Y() );
+ if ( edge._nodes.size() > 1 )
+ getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( edge._sWOL ), uv.X(), uv.Y() );
}
}
else
{
edge._2neibors = new _2NearEdges;
// target node instead of source ones will be set later
- if ( ! findNeiborsOnEdge( &edge,
- edge._2neibors->_nodes[0],
- edge._2neibors->_nodes[1],
- data))
- return false;
- edge.SetDataByNeighbors( edge._2neibors->_nodes[0],
- edge._2neibors->_nodes[1],
- helper);
+ // if ( ! findNeiborsOnEdge( &edge,
+ // edge._2neibors->_nodes[0],
+ // edge._2neibors->_nodes[1],
+ // data))
+ // return false;
+ // edge.SetDataByNeighbors( edge._2neibors->_nodes[0],
+ // edge._2neibors->_nodes[1],
+ // helper);
}
edge.SetCosin( edge._cosin ); // to update edge._lenFactor
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 &&
+ node->GetPosition()->GetDim() == 2 ) // node inside the FACE
{
- normal;
+ // probably NormEstim() failed due to a too high tolerance
+ pointKind = GeomLib::NormEstim( surface, uv, 1e-20, normal );
+ isOK = ( pointKind < IMPOSSIBLE );
+ }
+ if ( pointKind < IMPOSSIBLE )
+ {
+ if ( pointKind != REGULAR &&
+ !shiftInside &&
+ node->GetPosition()->GetDim() < 2 ) // FACE boundary
+ {
+ gp_XYZ normShift = getFaceNormal( node, face, helper, isOK, /*shiftInside=*/true );
+ if ( normShift * normal.XYZ() < 0. )
+ normal = normShift;
+ }
isOK = true;
}
- else // hard singularity
+
+ if ( !isOK ) // hard singularity, to call with shiftInside=true ?
{
const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face );
isOK = SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) normal.XYZ(), /*normalized=*/true );
if ( isOK )
{
- if ( helper.IsReversedSubMesh( face ))
+ TopoDS_Face ff = face;
+ ff.Orientation( TopAbs_FORWARD );
+ if ( helper.IsReversedSubMesh( ff ))
normal.Reverse();
break;
}
}
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
_SolidData& data)
{
const SMDS_MeshNode* node = edge->_nodes[0];
- const int shapeInd = node->getshapeId();
- SMESHDS_SubMesh* edgeSM = 0;
+ const int shapeInd = node->getshapeId();
+ SMESHDS_SubMesh* edgeSM = 0;
if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE )
{
-
edgeSM = getMeshDS()->MeshElements( shapeInd );
if ( !edgeSM || edgeSM->NbElements() == 0 )
return error(SMESH_Comment("Not meshed EDGE ") << shapeInd, data._index);
" average reached thickness is " << avgThick ));
}
+
+ // Restore position of src nodes moved by infaltion on _noShrinkShapes
+ dumpFunction(SMESH_Comment("restoNoShrink_So")<<data._index); // debug
+ int iBeg, iEnd = 0;
+ for ( int iS = 0; iS < data._endEdgeOnShape.size(); ++iS )
+ {
+ iBeg = iEnd;
+ iEnd = data._endEdgeOnShape[ iS ];
+ if ( data._edges[ iBeg ]->_nodes.size() == 1 )
+ for ( ; iBeg < iEnd; ++iBeg )
+ {
+ restoreNoShrink( *data._edges[ iBeg ] );
+ }
+ }
+ dumpFunctionEnd();
+
return true;
}
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 )
bndBox.Add( SMESH_TNodeXYZ( nIt->next() ));
gp_XYZ size = bndBox.CornerMax() - bndBox.CornerMin();
- SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->_nodes[0] );
- SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->_nodes[1] );
+ SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->tgtNode(0) );
+ SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->tgtNode(1) );
const double lineTol = 1e-2 * ( p0 - p1 ).Modulus();
for ( int i = 0; i < 3 && !isLine; ++i )
isLine = ( size.Coord( i+1 ) <= lineTol );
// set _2neibors according to the new order
for ( int i = iFrom; i < iTo-1; ++i )
- if ( _edges[i]->_2neibors->_nodes[1] != _edges[i+1]->_nodes.back() )
+ if ( _edges[i]->_2neibors->tgtNode(1) != _edges[i+1]->_nodes.back() )
_edges[i]->_2neibors->reverse();
if ( u2edge.size() > 1 &&
- _edges[iTo-1]->_2neibors->_nodes[0] != _edges[iTo-2]->_nodes.back() )
+ _edges[iTo-1]->_2neibors->tgtNode(0) != _edges[iTo-2]->_nodes.back() )
_edges[iTo-1]->_2neibors->reverse();
}
*/
//================================================================================
-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;
{
if ( F.IsNull() ) // 3D
{
- SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->_nodes[0]);
- SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->_nodes[1]);
+ SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->tgtNode(0));
+ SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->tgtNode(1));
for ( int i = iFrom; i < iTo; ++i )
{
double r = len[i-iFrom] / len.back();
}
else
{
- gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]);
- gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]);
- if ( data._edges[iFrom]->_2neibors->_nodes[0] ==
- data._edges[iTo-1]->_2neibors->_nodes[1] ) // closed edge
+ // gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->tgtNode(0));
+ // gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->tgtNode(1));
+ gp_XY uv0 = data._edges[iFrom]->_2neibors->_edges[0]->LastUV( F );
+ gp_XY uv1 = data._edges[iTo-1]->_2neibors->_edges[1]->LastUV( F );
+ if ( data._edges[iFrom]->_2neibors->tgtNode(0) ==
+ data._edges[iTo-1]->_2neibors->tgtNode(1) ) // closed edge
{
int iPeriodic = helper.GetPeriodicIndex();
if ( iPeriodic == 1 || iPeriodic == 2 )
if ( F.IsNull() ) // 3D
{
- if ( data._edges[iFrom]->_2neibors->_nodes[0] ==
- data._edges[iTo-1]->_2neibors->_nodes[1] )
+ if ( data._edges[iFrom]->_2neibors->tgtNode(0) ==
+ data._edges[iTo-1]->_2neibors->tgtNode(1) )
return true; // closed EDGE - nothing to do
return false; // TODO ???
{
const gp_XY center( center3D.X(), center3D.Y() );
- gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]);
- gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back());
- gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]);
+ gp_XY uv0 = data._edges[iFrom]->_2neibors->_edges[0]->LastUV( F );
+ gp_XY uvM = data._edges[iFrom]->LastUV( F );
+ gp_XY uv1 = data._edges[iTo-1]->_2neibors->_edges[1]->LastUV( F );
+ // gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->tgtNode(0));
+ // gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back());
+ // gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->tgtNode(1));
gp_Vec2d vec0( center, uv0 );
gp_Vec2d vecM( center, uvM );
gp_Vec2d vec1( center, uv1 );
const SMDS_MeshNode* tgt1 = edge->_nodes.back();
for ( int j = 0; j < 2; ++j ) // loop on _2NearEdges
{
- const SMDS_MeshNode* tgt2 = edge->_2neibors->_nodes[j];
+ const SMDS_MeshNode* tgt2 = edge->_2neibors->tgtNode(j);
pair< set< SMESH_TLink >::iterator, bool > link_isnew =
extrudedLinks.insert( SMESH_TLink( tgt1, tgt2 ));
if ( !link_isnew.second )
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->srcNode(0);
+ const SMDS_MeshNode * n2 = edge1->_2neibors->srcNode(1);
+ 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;
return segDir;
}
+//================================================================================
+/*!
+ * \brief Return the last position of the target node on a FACE.
+ * \param [in] F - the FACE this _LayerEdge is inflated along
+ * \return gp_XY - result UV
+ */
+//================================================================================
+
+gp_XY _LayerEdge::LastUV( const TopoDS_Face& F ) const
+{
+ if ( F.IsSame( _sWOL )) // F is my FACE
+ return gp_XY( _pos.back().X(), _pos.back().Y() );
+
+ if ( _sWOL.IsNull() || _sWOL.ShapeType() != TopAbs_EDGE ) // wrong call
+ return gp_XY( 1e100, 1e100 );
+
+ // _sWOL is EDGE of F; _pos.back().X() is the last U on the EDGE
+ double f, l, u = _pos.back().X();
+ Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( TopoDS::Edge(_sWOL), F, f,l);
+ if ( !C2d.IsNull() && f <= u && u <= l )
+ return C2d->Value( u ).XY();
+
+ return gp_XY( 1e100, 1e100 );
+}
+
//================================================================================
/*!
* \brief Test intersection of the last segment with a given triangle
/* 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.;
}
//================================================================================
SMESH_TNodeXYZ oldPos( tgtNode );
double dist01, distNewOld;
- SMESH_TNodeXYZ p0( _2neibors->_nodes[0]);
- SMESH_TNodeXYZ p1( _2neibors->_nodes[1]);
- dist01 = p0.Distance( _2neibors->_nodes[1] );
+ SMESH_TNodeXYZ p0( _2neibors->tgtNode(0));
+ SMESH_TNodeXYZ p1( _2neibors->tgtNode(1));
+ dist01 = p0.Distance( _2neibors->tgtNode(1) );
gp_Pnt newPos = p0 * _2neibors->_wgt[0] + p1 * _2neibors->_wgt[1];
double lenDelta = 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 )
double u = Precision::Infinite(); // to force projection w/o distance check
helper.CheckNodeU( TopoDS::Edge( _sWOL ), n, u, 1e-10, /*force=*/true, distXYZ );
_pos.back().SetCoord( u, 0, 0 );
- SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
- pos->SetUParameter( u );
+ if ( _nodes.size() > 1 )
+ {
+ SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
+ pos->SetUParameter( u );
+ }
}
else // TopAbs_FACE
{
gp_XY uv( Precision::Infinite(), 0 );
helper.CheckNodeUV( TopoDS::Face( _sWOL ), n, uv, 1e-10, /*force=*/true, distXYZ );
_pos.back().SetCoord( uv.X(), uv.Y(), 0 );
- SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
- pos->SetUParameter( uv.X() );
- pos->SetVParameter( uv.Y() );
+ if ( _nodes.size() > 1 )
+ {
+ SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
+ pos->SetUParameter( uv.X() );
+ pos->SetVParameter( uv.Y() );
+ }
}
n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]);
}
TNode2Edge* n2eMap = 0;
TNode2Edge::iterator n2e;
+ // Create intermediate nodes on each _LayerEdge
+
for ( size_t i = 0; i < data._edges.size(); ++i )
{
_LayerEdge& edge = *data._edges[i];
+ if ( edge._nodes.size() < 2 )
+ continue; // on _noShrinkShapes
+
// get accumulated length of segments
vector< double > segLen( edge._pos.size() );
segLen[0] = 0.0;
n2eMap = ( s2ne == data._s2neMap.end() ) ? 0 : n2eMap = s2ne->second;
prevBaseId = baseShapeId;
}
+ _LayerEdge* edgeOnSameNode = 0;
if ( n2eMap && (( n2e = n2eMap->find( edge._nodes[0] )) != n2eMap->end() ))
{
- _LayerEdge* foundEdge = n2e->second;
- const gp_XYZ& foundPos = foundEdge->_pos.back();
- SMDS_PositionPtr lastPos = tgtNode->GetPosition();
+ edgeOnSameNode = n2e->second;
+ const gp_XYZ& otherTgtPos = edgeOnSameNode->_pos.back();
+ SMDS_PositionPtr lastPos = tgtNode->GetPosition();
if ( isOnEdge )
{
SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( lastPos );
- epos->SetUParameter( foundPos.X() );
+ epos->SetUParameter( otherTgtPos.X() );
}
else
{
SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( lastPos );
- fpos->SetUParameter( foundPos.X() );
- fpos->SetVParameter( foundPos.Y() );
+ fpos->SetUParameter( otherTgtPos.X() );
+ fpos->SetVParameter( otherTgtPos.Y() );
}
}
// calculate height of the first layer
double r = ( segLen[iSeg] - hSum ) / ( segLen[iSeg] - segLen[iPrevSeg] );
gp_Pnt pos = r * edge._pos[iPrevSeg] + (1-r) * edge._pos[iSeg];
- SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >(edge._nodes[ iStep ]);
+ SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >( edge._nodes[ iStep ]);
if ( !edge._sWOL.IsNull() )
{
// compute XYZ by parameters <pos>
}
node->setXYZ( pos.X(), pos.Y(), pos.Z() );
}
+ } // loop on edge._nodes
+
+ if ( !edge._sWOL.IsNull() ) // prepare for shrink()
+ {
+ if ( isOnEdge )
+ edge._pos.back().SetCoord( u, 0,0);
+ else
+ edge._pos.back().SetCoord( uv.X(), uv.Y() ,0);
+
+ if ( edgeOnSameNode )
+ edgeOnSameNode->_pos.back() = edge._pos.back();
}
- }
+
+ } // loop on data._edges to create nodes
if ( !getMeshDS()->IsEmbeddedMode() )
// Log node movement
getMeshDS()->MoveNode( p._node, p.X(), p.Y(), p.Z() );
}
- // TODO: make quadratic prisms and polyhedrons(?)
+ // Create volumes
helper.SetElementsOnShape(true);
+ vector< vector<const SMDS_MeshNode*>* > nnVec;
+ set< int > degenEdgeInd;
+ vector<const SMDS_MeshElement*> degenVols;
+
TopExp_Explorer exp( data._solid, TopAbs_FACE );
for ( ; exp.More(); exp.Next() )
{
if ( data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
continue;
- SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
+ SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
- vector< vector<const SMDS_MeshNode*>* > nnVec;
while ( fIt->more() )
{
const SMDS_MeshElement* face = fIt->next();
- int nbNodes = face->NbCornerNodes();
+ const int nbNodes = face->NbCornerNodes();
nnVec.resize( nbNodes );
+ degenEdgeInd.clear();
+ int nbZ = 0;
SMDS_ElemIteratorPtr nIt = face->nodesIterator();
for ( int iN = 0; iN < nbNodes; ++iN )
{
const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
nnVec[ iN ] = & data._n2eMap[ n ]->_nodes;
+ if ( nnVec[ iN ]->size() < 2 )
+ degenEdgeInd.insert( iN );
+ else
+ nbZ = nnVec[ iN ]->size();
}
+ if ( nbZ == 0 )
+ continue;
- int nbZ = nnVec[0]->size();
switch ( nbNodes )
{
case 3:
- for ( int iZ = 1; iZ < nbZ; ++iZ )
- helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
- (*nnVec[0])[iZ], (*nnVec[1])[iZ], (*nnVec[2])[iZ]);
+ switch ( degenEdgeInd.size() )
+ {
+ case 0: // PENTA
+ {
+ for ( int iZ = 1; iZ < nbZ; ++iZ )
+ helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
+ (*nnVec[0])[iZ], (*nnVec[1])[iZ], (*nnVec[2])[iZ]);
+ break;
+ }
+ case 1: // PYRAM
+ {
+ int i2 = *degenEdgeInd.begin();
+ int i0 = helper.WrapIndex( i2 - 1, nbNodes );
+ int i1 = helper.WrapIndex( i2 + 1, nbNodes );
+ for ( int iZ = 1; iZ < nbZ; ++iZ )
+ helper.AddVolume( (*nnVec[i0])[iZ-1], (*nnVec[i1])[iZ-1],
+ (*nnVec[i1])[iZ], (*nnVec[i0])[iZ], (*nnVec[i2])[0]);
+ break;
+ }
+ case 2: // TETRA
+ {
+ int i3 = !degenEdgeInd.count(0) ? 0 : !degenEdgeInd.count(1) ? 1 : 2;
+ for ( int iZ = 1; iZ < nbZ; ++iZ )
+ helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
+ (*nnVec[i3])[iZ]);
+ break;
+ }
+ }
break;
+
case 4:
- for ( int iZ = 1; iZ < nbZ; ++iZ )
- helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1],
- (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1],
- (*nnVec[0])[iZ], (*nnVec[1])[iZ],
- (*nnVec[2])[iZ], (*nnVec[3])[iZ]);
+ switch ( degenEdgeInd.size() )
+ {
+ case 0: // HEX
+ {
+ for ( int iZ = 1; iZ < nbZ; ++iZ )
+ helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1],
+ (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1],
+ (*nnVec[0])[iZ], (*nnVec[1])[iZ],
+ (*nnVec[2])[iZ], (*nnVec[3])[iZ]);
+ break;
+ }
+ case 2: // PENTA?
+ {
+ int i2 = *degenEdgeInd.begin();
+ int i3 = *degenEdgeInd.rbegin();
+ bool ok = ( i3 - i2 == 1 );
+ if ( i2 == 0 && i3 == 3 ) { i2 = 3; i3 = 0; ok = true; }
+ int i0 = helper.WrapIndex( i3 + 1, nbNodes );
+ int i1 = helper.WrapIndex( i0 + 1, nbNodes );
+ for ( int iZ = 1; iZ < nbZ; ++iZ )
+ {
+ const SMDS_MeshElement* vol =
+ helper.AddVolume( (*nnVec[i3])[0], (*nnVec[i0])[iZ], (*nnVec[i0])[iZ-1],
+ (*nnVec[i2])[0], (*nnVec[i1])[iZ], (*nnVec[i1])[iZ-1]);
+ if ( !ok && vol )
+ degenVols.push_back( vol );
+ }
+ break;
+ }
+ case 3: // degen HEX
+ {
+ const SMDS_MeshNode* nn[8];
+ for ( int iZ = 1; iZ < nbZ; ++iZ )
+ {
+ const SMDS_MeshElement* vol =
+ helper.AddVolume( nnVec[0]->size() > 1 ? (*nnVec[0])[iZ-1] : (*nnVec[0])[0],
+ nnVec[1]->size() > 1 ? (*nnVec[1])[iZ-1] : (*nnVec[1])[0],
+ nnVec[2]->size() > 1 ? (*nnVec[2])[iZ-1] : (*nnVec[2])[0],
+ nnVec[3]->size() > 1 ? (*nnVec[3])[iZ-1] : (*nnVec[3])[0],
+ nnVec[0]->size() > 1 ? (*nnVec[0])[iZ] : (*nnVec[0])[0],
+ nnVec[1]->size() > 1 ? (*nnVec[1])[iZ] : (*nnVec[1])[0],
+ nnVec[2]->size() > 1 ? (*nnVec[2])[iZ] : (*nnVec[2])[0],
+ nnVec[3]->size() > 1 ? (*nnVec[3])[iZ] : (*nnVec[3])[0]);
+ degenVols.push_back( vol );
+ }
+ }
+ break;
+ }
break;
+
default:
return error("Not supported type of element", data._index);
- }
+
+ } // switch ( nbNodes )
+ } // while ( fIt->more() )
+ } // loop on FACEs
+
+ if ( !degenVols.empty() )
+ {
+ SMESH_ComputeErrorPtr& err = _mesh->GetSubMesh( data._solid )->GetComputeError();
+ if ( !err || err->IsOK() )
+ {
+ err.reset( new SMESH_ComputeError( COMPERR_WARNING,
+ "Degenerated volumes created" ));
+ err->myBadElements.insert( err->myBadElements.end(),
+ degenVols.begin(),degenVols.end() );
}
}
+
return true;
}
// EDGE's to shrink
map< TGeomID, _Shrinker1D > e2shrMap;
+ vector< _LayerEdge* > lEdges;
// loop on FACES to srink mesh on
map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin();
for ( ; f2sd != f2sdMap.end(); ++f2sd )
{
- _SolidData& data = *f2sd->second;
- TNode2Edge& n2eMap = data._n2eMap;
- const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first ));
-
- Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
-
+ _SolidData& data = *f2sd->second;
+ const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first ));
SMESH_subMesh* sm = _mesh->GetSubMesh( F );
SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
+ Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
+
helper.SetSubShape(F);
// ===========================
}
// Find _LayerEdge's inflated along F
- vector< _LayerEdge* > lEdges;
+ lEdges.clear();
{
- SMESH_subMeshIteratorPtr subIt =
- sm->getDependsOnIterator(/*includeSelf=*/false, /*complexShapeFirst=*/false);
+ set< TGeomID > subIDs;
+ SMESH_subMeshIteratorPtr subIt = sm->getDependsOnIterator(/*includeSelf=*/false);
while ( subIt->more() )
+ subIDs.insert( subIt->next()->GetId() );
+
+ int iBeg, iEnd = 0;
+ for ( int iS = 0; iS < data._endEdgeOnShape.size() && !subIDs.empty(); ++iS )
{
- SMESH_subMesh* sub = subIt->next();
- SMESHDS_SubMesh* subDS = sub->GetSubMeshDS();
- if ( subDS->NbNodes() == 0 || !n2eMap.count( subDS->GetNodes()->next() ))
- continue;
- SMDS_NodeIteratorPtr nIt = subDS->GetNodes();
- while ( nIt->more() )
- {
- _LayerEdge* edge = n2eMap[ nIt->next() ];
- lEdges.push_back( edge );
- prepareEdgeToShrink( *edge, F, helper, smDS );
- }
+ iBeg = iEnd;
+ iEnd = data._endEdgeOnShape[ iS ];
+ TGeomID shapeID = data._edges[ iBeg ]->_nodes[0]->getshapeId();
+ set< TGeomID >::iterator idIt = subIDs.find( shapeID );
+ if ( idIt == subIDs.end() ||
+ data._edges[ iBeg ]->_sWOL.IsNull() ) continue;
+ subIDs.erase( idIt );
+
+ if ( !data._noShrinkShapes.count( shapeID ))
+ for ( ; iBeg < iEnd; ++iBeg )
+ {
+ lEdges.push_back( data._edges[ iBeg ] );
+ prepareEdgeToShrink( *data._edges[ iBeg ], F, helper, smDS );
+ }
}
}
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.
const SMDS_MeshNode* srcNode = edge._nodes[0];
const SMDS_MeshNode* tgtNode = edge._nodes.back();
- edge._pos.clear();
-
if ( edge._sWOL.ShapeType() == TopAbs_FACE )
{
- gp_XY srcUV = helper.GetNodeUV( F, srcNode );
- gp_XY tgtUV = helper.GetNodeUV( F, tgtNode );
+ gp_XY srcUV( edge._pos[0].X(), edge._pos[0].Y() );//helper.GetNodeUV( F, srcNode );
+ gp_XY tgtUV = edge.LastUV( F ); //helper.GetNodeUV( F, tgtNode );
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);
}
else // _sWOL is TopAbs_EDGE
{
- TopoDS_Edge E = TopoDS::Edge( edge._sWOL);
+ const TopoDS_Edge& E = TopoDS::Edge( edge._sWOL );
SMESHDS_SubMesh* edgeSM = getMeshDS()->MeshElements( E );
if ( !edgeSM || edgeSM->NbElements() == 0 )
return error(SMESH_Comment("Not meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
double uSrc = helper.GetNodeU( E, srcNode, n2 );
double uTgt = helper.GetNodeU( E, tgtNode, srcNode );
- double u2 = helper.GetNodeU( E, n2, srcNode );
+ double u2 = helper.GetNodeU( E, n2, srcNode );
+
+ edge._pos.clear();
if ( fabs( uSrc-uTgt ) < 0.99 * fabs( uSrc-u2 ))
{
edge._simplices.resize( 1 );
edge._simplices[0]._nPrev = n2;
- // set UV of source node to target node
+ // set U of source node to the target node
SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
pos->SetUParameter( uSrc );
}
return true;
}
+//================================================================================
+/*!
+ * \brief Restore position of a sole node of a _LayerEdge based on _noShrinkShapes
+ */
+//================================================================================
+
+void _ViscousBuilder::restoreNoShrink( _LayerEdge& edge ) const
+{
+ if ( edge._nodes.size() == 1 )
+ {
+ edge._pos.clear();
+ edge._len = 0;
+
+ const SMDS_MeshNode* srcNode = edge._nodes[0];
+ TopoDS_Shape S = SMESH_MesherHelper::GetSubShapeByNode( srcNode, getMeshDS() );
+ if ( S.IsNull() ) return;
+
+ gp_Pnt p;
+
+ switch ( S.ShapeType() )
+ {
+ case TopAbs_EDGE:
+ {
+ double f,l;
+ TopLoc_Location loc;
+ Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( S ), loc, f, l );
+ if ( curve.IsNull() ) return;
+ SMDS_EdgePosition* ePos = static_cast<SMDS_EdgePosition*>( srcNode->GetPosition() );
+ p = curve->Value( ePos->GetUParameter() );
+ break;
+ }
+ case TopAbs_VERTEX:
+ {
+ p = BRep_Tool::Pnt( TopoDS::Vertex( S ));
+ break;
+ }
+ default: return;
+ }
+ getMeshDS()->MoveNode( srcNode, p.X(), p.Y(), p.Z() );
+ dumpMove( srcNode );
+ }
+}
+
//================================================================================
/*!
* \brief Try to fix triangles with high aspect ratio by swaping diagonals
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();
}
else // _sWOL is TopAbs_EDGE
{
- TopoDS_Edge E = TopoDS::Edge( _sWOL );
- const SMDS_MeshNode* n2 = _simplices[0]._nPrev;
+ const TopoDS_Edge& E = TopoDS::Edge( _sWOL );
+ const SMDS_MeshNode* n2 = _simplices[0]._nPrev;
SMDS_EdgePosition* tgtPos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
- const double u2 = helper.GetNodeU( E, n2, tgtNode );
+ const double u2 = helper.GetNodeU( E, n2, tgtNode );
const double uSrc = _pos[0].Coord( U_SRC );
const double lenTgt = _pos[0].Coord( LEN_TGT );
for ( int iE = 1; iE <= geomEdges.Extent(); ++iE )
{
const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE));
+ if ( data._noShrinkShapes.count( getMeshDS()->ShapeToIndex( E )))
+ continue;
// Get _LayerEdge's based on E
{
vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes;
vector< const SMDS_MeshNode*>& nn2 = ledges[j-dj2]->_nodes;
- if ( isOnFace )
- for ( size_t z = 1; z < nn1.size(); ++z )
- sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
+ if ( nn1.size() == nn2.size() )
+ {
+ if ( isOnFace )
+ for ( size_t z = 1; z < nn1.size(); ++z )
+ sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
+ else
+ for ( size_t z = 1; z < nn1.size(); ++z )
+ sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
+ }
+ else if ( nn1.size() == 1 )
+ {
+ if ( isOnFace )
+ for ( size_t z = 1; z < nn2.size(); ++z )
+ sm->AddElement( getMeshDS()->AddFace( nn1[0], nn2[z-1], nn2[z] ));
+ else
+ for ( size_t z = 1; z < nn2.size(); ++z )
+ sm->AddElement( new SMDS_FaceOfNodes( nn1[0], nn2[z-1], nn2[z] ));
+ }
else
- for ( size_t z = 1; z < nn1.size(); ++z )
- sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z]));
+ {
+ if ( isOnFace )
+ for ( size_t z = 1; z < nn1.size(); ++z )
+ sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[0], nn1[z] ));
+ else
+ for ( size_t z = 1; z < nn1.size(); ++z )
+ sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[0], nn2[z] ));
+ }
}
// Make edges
if ( !edge->_sWOL.IsNull() && edge->_sWOL.ShapeType() == TopAbs_EDGE )
{
vector< const SMDS_MeshNode*>& nn = edge->_nodes;
- if ( nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() )
+ if ( nn.size() < 2 || nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() )
continue;
helper.SetSubShape( edge->_sWOL );
helper.SetElementsOnShape( true );