-// Copyright (C) 2007-2011 CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
#include "SMESHDS_Hypothesis.hxx"
#include "SMESH_Algo.hxx"
#include "SMESH_ComputeError.hxx"
+#include "SMESH_ControlsDef.hxx"
#include "SMESH_Gen.hxx"
#include "SMESH_Group.hxx"
#include "SMESH_Mesh.hxx"
#include "SMESH_MesherHelper.hxx"
+#include "SMESH_ProxyMesh.hxx"
#include "SMESH_subMesh.hxx"
#include "SMESH_subMeshEventListener.hxx"
-#include "SMESH_ProxyMesh.hxx"
#include "utilities.h"
+#include <BRepAdaptor_Curve2d.hxx>
#include <BRep_Tool.hxx>
#include <Bnd_B2d.hxx>
#include <Bnd_B3d.hxx>
#include <Geom_Line.hxx>
#include <Geom_TrimmedCurve.hxx>
#include <Precision.hxx>
+#include <Standard_ErrorHandler.hxx>
+#include <TColStd_Array1OfReal.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>
#include <TopTools_IndexedMapOfShape.hxx>
#include <gp_Ax1.hxx>
#include <gp_Vec.hxx>
#include <gp_XY.hxx>
-#include <gp_XYZ.hxx>
#include <list>
#include <string>
-#include <math.h>
+#include <cmath>
#include <limits>
//#define __myDEBUG
using namespace std;
//================================================================================
-namespace VISCOUS
+namespace VISCOUS_3D
{
typedef int TGeomID;
//--------------------------------------------------------------------------------
/*!
* \brief Listener of events of 3D sub-meshes computed with viscous layers.
- * It is used to clear an inferior dim sub-mesh modified by viscous layers
+ * It is used to clear an inferior dim sub-meshes modified by viscous layers
*/
- class _SrinkShapeListener : SMESH_subMeshEventListener
+ class _ShrinkShapeListener : SMESH_subMeshEventListener
{
- _SrinkShapeListener(): SMESH_subMeshEventListener(/*isDeletable=*/false) {}
- static SMESH_subMeshEventListener* Get() { static _SrinkShapeListener l; return &l; }
+ _ShrinkShapeListener()
+ : SMESH_subMeshEventListener(/*isDeletable=*/false,
+ "StdMeshers_ViscousLayers::_ShrinkShapeListener") {}
public:
+ static SMESH_subMeshEventListener* Get() { static _ShrinkShapeListener l; return &l; }
virtual void ProcessEvent(const int event,
const int eventType,
SMESH_subMesh* solidSM,
SMESH_subMeshEventListener::ProcessEvent(event,eventType,solidSM,data,hyp);
}
}
- static void ToClearSubMeshWithSolid( SMESH_subMesh* sm,
- const TopoDS_Shape& solid)
- {
- SMESH_subMesh* solidSM = sm->GetFather()->GetSubMesh( solid );
- SMESH_subMeshEventListenerData* data = solidSM->GetEventListenerData( Get());
- if ( data )
- {
- if ( find( data->mySubMeshes.begin(), data->mySubMeshes.end(), sm ) ==
- data->mySubMeshes.end())
- data->mySubMeshes.push_back( sm );
- }
- else
- {
- data = SMESH_subMeshEventListenerData::MakeData( /*dependent=*/sm );
- sm->SetEventListener( Get(), data, /*whereToListenTo=*/solidSM );
- }
- }
};
//--------------------------------------------------------------------------------
/*!
*/
class _ViscousListener : SMESH_subMeshEventListener
{
- _ViscousListener(): SMESH_subMeshEventListener(/*isDeletable=*/false) {}
+ _ViscousListener():
+ SMESH_subMeshEventListener(/*isDeletable=*/false,
+ "StdMeshers_ViscousLayers::_ViscousListener") {}
static SMESH_subMeshEventListener* Get() { static _ViscousListener l; return &l; }
public:
virtual void ProcessEvent(const int event,
}
};
+ //================================================================================
+ /*!
+ * \brief sets a sub-mesh event listener to clear sub-meshes of sub-shapes of
+ * the main shape when sub-mesh of the main shape is cleared,
+ * for example to clear sub-meshes of FACEs when sub-mesh of a SOLID
+ * is cleared
+ */
+ //================================================================================
+
+ void ToClearSubWithMain( SMESH_subMesh* sub, const TopoDS_Shape& main)
+ {
+ SMESH_subMesh* mainSM = sub->GetFather()->GetSubMesh( main );
+ SMESH_subMeshEventListenerData* data =
+ mainSM->GetEventListenerData( _ShrinkShapeListener::Get());
+ if ( data )
+ {
+ if ( find( data->mySubMeshes.begin(), data->mySubMeshes.end(), sub ) ==
+ data->mySubMeshes.end())
+ data->mySubMeshes.push_back( sub );
+ }
+ else
+ {
+ data = SMESH_subMeshEventListenerData::MakeData( /*dependent=*/sub );
+ sub->SetEventListener( _ShrinkShapeListener::Get(), data, /*whereToListenTo=*/mainSM );
+ }
+ }
//--------------------------------------------------------------------------------
/*!
* \brief Simplex (triangle or tetrahedron) based on 1 (tria) or 2 (tet) nodes of
- M[0][2]*M[1][1]*M[2][0]);
return determinant > 1e-100;
}
- bool IsForward(const gp_XY& tgtUV,
- const TopoDS_Face& face,
- SMESH_MesherHelper& helper,
- const double refSign) const
+ bool IsForward(const gp_XY& tgtUV,
+ const SMDS_MeshNode* smoothedNode,
+ const TopoDS_Face& face,
+ SMESH_MesherHelper& helper,
+ const double refSign) const
{
- gp_XY prevUV = helper.GetNodeUV( face, _nPrev );
- gp_XY nextUV = helper.GetNodeUV( face, _nNext );
+ gp_XY prevUV = helper.GetNodeUV( face, _nPrev, smoothedNode );
+ gp_XY nextUV = helper.GetNodeUV( face, _nNext, smoothedNode );
gp_Vec2d v1( tgtUV, prevUV ), v2( tgtUV, nextUV );
double d = v1 ^ v2;
return d*refSign > 1e-100;
}
+ bool IsNeighbour(const _Simplex& other) const
+ {
+ return _nPrev == other._nNext || _nNext == other._nPrev;
+ }
};
//--------------------------------------------------------------------------------
/*!
{
double _r; // radius
double _k; // factor to correct node smoothed position
+ double _h2lenRatio; // avgNormProj / (2*avgDist)
public:
static _Curvature* New( double avgNormProj, double avgDist )
{
c->_r = avgDist * avgDist / avgNormProj;
c->_k = avgDist * avgDist / c->_r / c->_r;
c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive
+ c->_h2lenRatio = avgNormProj / ( avgDist + avgDist );
}
return c;
}
double lenDelta(double len) const { return _k * ( _r + len ); }
+ double lenDeltaByDist(double dist) const { return dist * _h2lenRatio; }
};
struct _LayerEdge;
//--------------------------------------------------------------------------------
Handle(Geom_Surface)& surface,
SMESH_MesherHelper& helper,
const double refSign,
+ bool isCentroidal,
bool set3D);
};
//--------------------------------------------------------------------------------
_SolidData& data);
void getSimplices( const SMDS_MeshNode* node, vector<_Simplex>& simplices,
const set<TGeomID>& ingnoreShapes,
- const _SolidData* dataToCheckOri = 0);
+ const _SolidData* dataToCheckOri = 0,
+ const bool toSort = false);
bool sortEdges( _SolidData& data,
vector< vector<_LayerEdge*> >& edgesByGeom);
void limitStepSize( _SolidData& data,
bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F,
SMESH_MesherHelper& helper,
const SMESHDS_SubMesh* faceSubMesh );
+ void fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper);
bool addBoundaryElements();
bool error( const string& text, int solidID=-1 );
virtual SMDSAbs_ElementType GetType() const { return SMDSAbs_Face; }
virtual vtkIdType GetVtkType() const { return -1; }
virtual SMDSAbs_EntityType GetEntityType() const { return SMDSEntity_Last; }
- virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const
+ virtual SMDSAbs_GeometryType GetGeomType() const { return SMDSGeom_TRIANGLE; }
+virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const
{ return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));}
};
//--------------------------------------------------------------------------------
_nn[3]=_le2->_nodes[0];
}
};
-} // namespace VISCOUS
+} // namespace VISCOUS_3D
//================================================================================
// StdMeshers_ViscousLayers hypothesis
_name = StdMeshers_ViscousLayers::GetHypType();
_param_algo_dim = -3; // auxiliary hyp used by 3D algos
} // --------------------------------------------------------------------------------
-void StdMeshers_ViscousLayers::SetIgnoreFaces(const std::vector<int>& faceIds)
+void StdMeshers_ViscousLayers::SetBndShapesToIgnore(const std::vector<int>& faceIds)
+{
+ if ( faceIds != _ignoreBndShapeIds )
+ _ignoreBndShapeIds = faceIds, NotifySubMeshesHypothesisModification();
+} // --------------------------------------------------------------------------------
+bool StdMeshers_ViscousLayers::IsIgnoredShape(const int shapeID) const
{
- if ( faceIds != _ignoreFaceIds )
- _ignoreFaceIds = faceIds, NotifySubMeshesHypothesisModification();
+ return ( find( _ignoreBndShapeIds.begin(), _ignoreBndShapeIds.end(), shapeID )
+ != _ignoreBndShapeIds.end() );
} // --------------------------------------------------------------------------------
void StdMeshers_ViscousLayers::SetTotalThickness(double thickness)
{
const TopoDS_Shape& theShape,
const bool toMakeN2NMap) const
{
- using namespace VISCOUS;
+ using namespace VISCOUS_3D;
_ViscousBuilder bulder;
SMESH_ComputeErrorPtr err = bulder.Compute( theMesh, theShape );
if ( err && !err->IsOK() )
save << " " << _nbLayers
<< " " << _thickness
<< " " << _stretchFactor
- << " " << _ignoreFaceIds.size();
- for ( unsigned i = 0; i < _ignoreFaceIds.size(); ++i )
- save << " " << _ignoreFaceIds[i];
+ << " " << _ignoreBndShapeIds.size();
+ for ( unsigned i = 0; i < _ignoreBndShapeIds.size(); ++i )
+ save << " " << _ignoreBndShapeIds[i];
return save;
} // --------------------------------------------------------------------------------
std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load)
{
int nbFaces, faceID;
load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces;
- while ( _ignoreFaceIds.size() < nbFaces && load >> faceID )
- _ignoreFaceIds.push_back( faceID );
+ while ( _ignoreBndShapeIds.size() < nbFaces && load >> faceID )
+ _ignoreBndShapeIds.push_back( faceID );
return load;
} // --------------------------------------------------------------------------------
bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh* theMesh,
while ( eIt->more())
{
const TopoDS_Edge* e = static_cast<const TopoDS_Edge*>( eIt->next() );
- if ( helper.IsSubShape( *e, F ) && BRep_Tool::Curve( *e, loc,f,l))
+ if ( helper.IsSubShape( *e, F ) && !BRep_Tool::Curve( *e, loc,f,l).IsNull() )
edges.push_back( *e );
}
gp_XYZ dir(0,0,0);
if ( !( ok = ( edges.size() > 0 ))) return dir;
// get average dir of edges going fromV
- gp_Vec edgeDir;
- for ( unsigned i = 0; i < edges.size(); ++i )
- {
- edgeDir = getEdgeDir( edges[i], fromV );
- double size2 = edgeDir.SquareMagnitude();
- if ( size2 > numeric_limits<double>::min() )
- edgeDir /= sqrt( size2 );
- else
- ok = false;
- dir += edgeDir.XYZ();
- }
+ gp_XYZ edgeDir;
+ //if ( edges.size() > 1 )
+ for ( unsigned i = 0; i < edges.size(); ++i )
+ {
+ edgeDir = getEdgeDir( edges[i], fromV );
+ double size2 = edgeDir.SquareModulus();
+ if ( size2 > numeric_limits<double>::min() )
+ edgeDir /= sqrt( size2 );
+ else
+ ok = false;
+ dir += edgeDir;
+ }
gp_XYZ fromEdgeDir = getFaceDir( F, edges[0], node, helper, ok );
- if ( edges.size() == 1 || dir.SquareModulus() < 1e-10)
+ if ( edges.size() == 1 )
dir = fromEdgeDir;
+ else if ( dir.SquareModulus() < 0.1 ) // ~< 20 degrees
+ dir = fromEdgeDir + getFaceDir( F, edges[1], node, helper, ok );
else if ( dir * fromEdgeDir < 0 )
dir *= -1;
if ( ok )
{
//dir /= edges.size();
if ( cosin ) {
- double angle = edgeDir.Angle( dir );
+ double angle = gp_Vec( edgeDir ).Angle( dir );
*cosin = cos( angle );
}
}
return dir;
}
+ //================================================================================
+ /*!
+ * \brief Returns true if a FACE is bound by a concave EDGE
+ */
+ //================================================================================
+
+ bool isConcave( const TopoDS_Face& F, SMESH_MesherHelper& helper )
+ {
+ gp_Vec2d drv1, drv2;
+ gp_Pnt2d p;
+ TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE );
+ for ( ; eExp.More(); eExp.Next() )
+ {
+ const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
+ if ( BRep_Tool::Degenerated( E )) continue;
+ // check if 2D curve is concave
+ BRepAdaptor_Curve2d curve( E, F );
+ const int nbIntervals = curve.NbIntervals( GeomAbs_C2 );
+ TColStd_Array1OfReal intervals(1, nbIntervals + 1 );
+ curve.Intervals( intervals, GeomAbs_C2 );
+ bool isConvex = true;
+ for ( int i = 1; i <= nbIntervals && isConvex; ++i )
+ {
+ double u1 = intervals( i );
+ double u2 = intervals( i+1 );
+ curve.D2( 0.5*( u1+u2 ), p, drv1, drv2 );
+ double cross = drv2 ^ drv1;
+ if ( E.Orientation() == TopAbs_REVERSED )
+ cross = -cross;
+ isConvex = ( cross < 1e-9 );
+ }
+ // check if concavity is strong enough to care about it
+ //const double maxAngle = 5 * Standard_PI180;
+ if ( !isConvex )
+ {
+ //cout << "Concave FACE " << helper.GetMeshDS()->ShapeToIndex( F ) << endl;
+ return true;
+ // map< double, const SMDS_MeshNode* > u2nodes;
+ // if ( !SMESH_Algo::GetSortedNodesOnEdge( helper.GetMeshDS(), E,
+ // /*ignoreMedium=*/true, u2nodes))
+ // continue;
+ // map< double, const SMDS_MeshNode* >::iterator u2n = u2nodes.begin();
+ // gp_Pnt2d uvPrev = helper.GetNodeUV( F, u2n->second );
+ // double uPrev = u2n->first;
+ // for ( ++u2n; u2n != u2nodes.end(); ++u2n )
+ // {
+ // gp_Pnt2d uv = helper.GetNodeUV( F, u2n->second );
+ // gp_Vec2d segmentDir( uvPrev, uv );
+ // curve.D1( uPrev, p, drv1 );
+ // try {
+ // if ( fabs( segmentDir.Angle( drv1 )) > maxAngle )
+ // return true;
+ // }
+ // catch ( ... ) {}
+ // uvPrev = uv;
+ // uPrev = u2n->first;
+ // }
+ }
+ }
+ return false;
+ }
//--------------------------------------------------------------------------------
// DEBUG. Dump intermediate node positions into a python script
#ifdef __myDEBUG
py = new ofstream(fname);
*py << "from smesh import *" << endl
<< "meshSO = GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
- << "mesh = Mesh( meshSO.GetObject()._narrow( SMESH.SMESH_Mesh ))"<<endl;
+ << "mesh = Mesh( meshSO.GetObject() )"<<endl;
}
- ~PyDump() {
- *py << "mesh.MakeGroup('Prisms of viscous layers',VOLUME,FT_ElemGeomType,'=',Geom_PENTA)"
- <<endl; delete py; py=0;
+ void Finish() {
+ if (py)
+ *py << "mesh.MakeGroup('Viscous Prisms',VOLUME,FT_ElemGeomType,'=',Geom_PENTA)"<<endl;
+ delete py; py=0;
}
+ ~PyDump() { Finish(); }
};
#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)
- { *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl;}
+ { if (py) *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl;}
void _dumpMove(const SMDS_MeshNode* n, int ln)
- { *py<< " mesh.MoveNode( "<<n->GetID()<< ", "<< n->X()
- << ", "<<n->Y()<<", "<< n->Z()<< ")\t\t # "<< ln <<endl; }
+ { if (py) *py<< " mesh.MoveNode( "<<n->GetID()<< ", "<< n->X()
+ << ", "<<n->Y()<<", "<< n->Z()<< ")\t\t # "<< ln <<endl; }
void _dumpCmd(const string& txt, int ln)
- { *py<< " "<<txt<<" # "<< ln <<endl; }
+ { if (py) *py<< " "<<txt<<" # "<< ln <<endl; }
void dumpFunctionEnd()
- { *py<< " return"<< endl; }
+ { if (py) *py<< " return"<< endl; }
+ void dumpChangeNodes( const SMDS_MeshElement* f )
+ { if (py) { *py<< " mesh.ChangeElemNodes( " << f->GetID()<<", [";
+ for ( int i=1; i < f->NbNodes(); ++i ) *py << f->GetNode(i-1)->GetID()<<", ";
+ *py << f->GetNode( f->NbNodes()-1 )->GetID() << " ])"<< endl; }}
#else
- struct PyDump { PyDump() {} };
- void dumpFunction(const string& fun ){}
- void dumpFunctionEnd() {}
- void dumpMove(const SMDS_MeshNode* n ){}
- void dumpCmd(const string& txt){}
+ struct PyDump { void Finish() {} };
+#define dumpFunction(f) f
+#define dumpMove(n)
+#define dumpCmd(txt)
+#define dumpFunctionEnd()
+#define dumpChangeNodes(f)
#endif
}
-using namespace VISCOUS;
+using namespace VISCOUS_3D;
//================================================================================
/*!
{
if ( ! makeLayer(_sdVec[i]) )
return _error;
+
+ if ( _sdVec[i]._edges.size() == 0 )
+ continue;
if ( ! inflate(_sdVec[i]) )
return _error;
return _error;
addBoundaryElements();
-
+
makeGroupOfLE(); // debug
+ debugDump.Finish();
return _error;
}
//================================================================================
/*!
- * \brief Finds SOLIDs to compute using viscous layers. Fill _sdVec
+ * \brief Finds SOLIDs to compute using viscous layers. Fills _sdVec
*/
//================================================================================
vector<TopoDS_Shape> ignoreFaces;
for ( unsigned i = 0; i < _sdVec.size(); ++i )
{
- vector<TGeomID> ids = _sdVec[i]._hyp->GetIgnoreFaces();
+ vector<TGeomID> ids = _sdVec[i]._hyp->GetBndShapesToIgnore();
for ( unsigned i = 0; i < ids.size(); ++i )
{
const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[i] );
{
_ignoreShapeIds.insert( faceInd );
ignoreFaces.push_back( exp.Current() );
- if ( SMESH_Algo::IsReversedSubMesh( TopoDS::Face( exp.Current() ), getMeshDS()))
+ if ( helper.IsReversedSubMesh( TopoDS::Face( exp.Current() )))
_sdVec[i]._reversedFaceIds.insert( faceInd );
}
}
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++;
+ }
}
}
{
case 1:
{
- _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] )); break;
+ helper.SetSubShape( facesWOL[0] );
+ if ( helper.IsRealSeam( vInd )) // inflate along a seam edge?
+ {
+ TopoDS_Shape seamEdge;
+ PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
+ while ( eIt->more() && seamEdge.IsNull() )
+ {
+ const TopoDS_Shape* e = eIt->next();
+ if ( helper.IsRealSeam( *e ) )
+ seamEdge = *e;
+ }
+ if ( !seamEdge.IsNull() )
+ {
+ _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge ));
+ break;
+ }
+ }
+ _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] ));
+ break;
}
case 2:
{
if ( data._stepSize < 1. )
data._epsilon *= data._stepSize;
- // Put _LayerEdge's into a vector
+ // Put _LayerEdge's into the vector data._edges
if ( !sortEdges( data, edgesByGeom ))
return false;
list< TGeomID > shapesToSmooth;
SMESH_MesherHelper helper( *_mesh );
- bool ok;
+ bool ok = true;
for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS )
{
}
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]
+ gp_XYZ inFaceDir = getFaceDir( F, V, node, helper, normOK);
+ double angle = gp_Vec( inFaceDir).Angle( edge._normal ); // [0,PI]
edge._cosin = cos( angle );
//cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl;
break;
void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node,
vector<_Simplex>& simplices,
const set<TGeomID>& ingnoreShapes,
- const _SolidData* dataToCheckOri)
+ const _SolidData* dataToCheckOri,
+ const bool toSort)
{
- SMESH_MeshEditor editor( _mesh );
SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
while ( fIt->more() )
{
const SMDS_MeshElement* f = fIt->next();
- const TGeomID shapeInd = editor.FindShape( f );
+ const TGeomID shapeInd = f->getshapeId();
if ( ingnoreShapes.count( shapeInd )) continue;
const int nbNodes = f->NbCornerNodes();
int srcInd = f->GetNodeIndex( node );
std::swap( nPrev, nNext );
simplices.push_back( _Simplex( nPrev, nNext ));
}
- simplices.resize( simplices.size() );
+
+ if ( toSort )
+ {
+ vector<_Simplex> sortedSimplices( simplices.size() );
+ sortedSimplices[0] = simplices[0];
+ int nbFound = 0;
+ for ( size_t i = 1; i < simplices.size(); ++i )
+ {
+ for ( size_t j = 1; j < simplices.size(); ++j )
+ if ( sortedSimplices[i-1]._nNext == simplices[j]._nPrev )
+ {
+ sortedSimplices[i] = simplices[j];
+ nbFound++;
+ break;
+ }
+ }
+ if ( nbFound == simplices.size() - 1 )
+ simplices.swap( sortedSimplices );
+ }
}
//================================================================================
{
if ( data._edges[i]->IsOnEdge() ) continue;
data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon );
- if ( geomSize > intersecDist )
+ if ( geomSize > intersecDist && intersecDist > 0 )
geomSize = intersecDist;
}
if ( data._stepSize > 0.3 * geomSize )
TGeomID sInd = data._edges[ iBeg ]->_nodes[0]->getshapeId();
if ( data._edges[ iBeg ]->IsOnEdge() )
- { // try a simple solution on an analytic EDGE
+ {
+ dumpFunction(SMESH_Comment("smooth")<<data._index << "_Ed"<<sInd <<"_InfStep"<<nbSteps);
+
+ // try a simple solution on an analytic EDGE
if ( !smoothAnalyticEdge( data, iBeg, iEnd, surface, F, helper ))
{
- dumpFunction(SMESH_Comment("smooth")<<data._index << "_Ed"<<sInd <<"_InfStep"<<nbSteps);
// smooth on EDGE's
int step = 0;
do {
}
while ( moved && step++ < 5 );
//cout << " NB STEPS: " << step << endl;
-
- dumpFunctionEnd();
}
+ dumpFunctionEnd();
}
else
{
data._edges[i]->_pos.back() = newPos;
SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
+ dumpMove( tgtNode );
}
}
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
+ {
+ int iPeriodic = helper.GetPeriodicIndex();
+ if ( iPeriodic == 1 || iPeriodic == 2 )
+ {
+ uv1.SetCoord( iPeriodic, helper.GetOtherParam( uv1.Coord( iPeriodic )));
+ if ( uv0.Coord( iPeriodic ) > uv1.Coord( iPeriodic ))
+ std::swap( uv0, uv1 );
+ }
+ }
+ const gp_XY rangeUV = uv1 - uv0;
for ( int i = iFrom; i < iTo; ++i )
{
double r = len[i-iFrom] / len.back();
- gp_XY newUV = uv0 * ( 1. - r ) + uv1 * r;
+ gp_XY newUV = uv0 + r * rangeUV;
data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
+ dumpMove( tgtNode );
SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
pos->SetUParameter( newUV.X() );
else // 2D
{
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_Vec2d vec0( center, uv0 );
- gp_Vec2d vecM( center, uvM);
+ gp_Vec2d vecM( center, uvM );
gp_Vec2d vec1( center, uv1 );
double uLast = vec0.Angle( vec1 ); // -PI - +PI
double uMidl = vec0.Angle( vecM );
- if ( uLast < 0 ) uLast += 2*PI; // 0.0 - 2*PI
- if ( uMidl < 0 ) uMidl += 2*PI;
- const bool sense = ( uMidl < uLast );
+ if ( uLast * uMidl < 0. )
+ uLast += ( uMidl > 0 ? +2. : -2. ) * M_PI;
const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() );
- gp_Ax2d axis( center, vec0 );
- gp_Circ2d circ ( axis, radius, sense );
+ gp_Ax2d axis( center, vec0 );
+ gp_Circ2d circ( axis, radius );
for ( int i = iFrom; i < iTo; ++i )
{
double newU = uLast * len[i-iFrom] / len.back();
gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
+ dumpMove( tgtNode );
SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
pos->SetUParameter( newUV.X() );
if ( S.ShapeType() != TopAbs_EDGE )
continue; // TODO: find EDGE by VERTEX
E1 = TopoDS::Edge( S );
- set< _LayerEdge* >::iterator eIt = ee.begin();
+ set< _LayerEdge*, _LayerEdgeCmp >::iterator eIt = ee.begin();
while ( E2.IsNull() && eIt != ee.end())
{
_LayerEdge* e2 = *eIt++;
}
if ( intFound )
{
- if ( dist < segLen*(1.01))
+ if ( dist < segLen*(1.01) && dist > -(_len-segLen) )
segmentIntersected = true;
if ( distance > dist )
distance = dist, iFace = j;
double lenDelta = 0;
if ( _curvature )
{
- lenDelta = _curvature->lenDelta( _len );
+ //lenDelta = _curvature->lenDelta( _len );
+ lenDelta = _curvature->lenDeltaByDist( dist01 );
newPos.ChangeCoord() += _normal * lenDelta;
}
}
}
+ if ( !getMeshDS()->IsEmbeddedMode() )
+ // Log node movement
+ for ( unsigned i = 0; i < data._edges.size(); ++i )
+ {
+ _LayerEdge& edge = *data._edges[i];
+ SMESH_TNodeXYZ p ( edge._nodes.back() );
+ getMeshDS()->MoveNode( p._node, p.X(), p.Y(), p.Z() );
+ }
+
// TODO: make quadratic prisms and polyhedrons(?)
helper.SetElementsOnShape(true);
SMDS_ElemIteratorPtr fIt = smDS->GetElements();
while ( fIt->more() )
proxySub->AddElement( fIt->next() );
+ // as a result 3D algo will use elements from proxySub and not from smDS
}
}
}
SMESH_MesherHelper helper( *_mesh );
+ helper.ToFixNodeParameters( true );
// EDGE's to shrink
- map< int, _Shrinker1D > e2shrMap;
+ map< TGeomID, _Shrinker1D > e2shrMap;
// loop on FACES to srink mesh on
map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin();
_SolidData& data = *f2sd->second;
TNode2Edge& n2eMap = data._n2eMap;
const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first ));
- const bool reverse = ( data._reversedFaceIds.count( f2sd->first ));
Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
// Prepare data for shrinking
// ===========================
- // Collect nodes to smooth as src nodes are not yet replaced by tgt ones
- // and thus all nodes on FACE connected to 2d elements are to be smoothed
+ // Collect nodes to smooth, as src nodes are not yet replaced by tgt ones
+ // and thus all nodes on a FACE connected to 2d elements are to be smoothed
vector < const SMDS_MeshNode* > smoothNodes;
{
SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
// Find out face orientation
double refSign = 1;
const set<TGeomID> ignoreShapes;
+ bool isOkUV;
if ( !smoothNodes.empty() )
{
- gp_XY uv = helper.GetNodeUV( F, smoothNodes[0] );
vector<_Simplex> simplices;
getSimplices( smoothNodes[0], simplices, ignoreShapes );
- if ( simplices[0].IsForward(uv, F, helper,refSign) != (!reverse))
+ helper.GetNodeUV( F, simplices[0]._nPrev, 0, &isOkUV ); // fix UV of silpmex nodes
+ helper.GetNodeUV( F, simplices[0]._nNext, 0, &isOkUV );
+ gp_XY uv = helper.GetNodeUV( F, smoothNodes[0], 0, &isOkUV );
+ if ( !simplices[0].IsForward(uv, smoothNodes[0], F, helper,refSign) )
refSign = -1;
}
}
}
+ // find out if a FACE is concave
+ const bool isConcaveFace = isConcave( F, helper );
+
// Create _SmoothNode's on face F
vector< _SmoothNode > nodesToSmooth( smoothNodes.size() );
{
const SMDS_MeshNode* n = smoothNodes[i];
nodesToSmooth[ i ]._node = n;
// src nodes must be replaced by tgt nodes to have tgt nodes in _simplices
- getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes );
+ getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, isConcaveFace );
+ // fix up incorrect uv of nodes on the FACE
+ helper.GetNodeUV( F, n, 0, &isOkUV);
dumpMove( n );
}
dumpFunctionEnd();
_Shrinker1D& srinker = e2shrMap[ edgeIndex ];
eShri1D.insert( & srinker );
srinker.AddEdge( edge, helper );
+ VISCOUS_3D::ToClearSubWithMain( _mesh->GetSubMesh( edge->_sWOL ), data._solid );
// restore params of nodes on EGDE if the EDGE has been already
// srinked while srinking another FACE
srinker.RestoreParams();
shrinked |= lEdges[i]->SetNewLength2d( surface,F,helper );
}
dumpFunctionEnd();
- if ( !shrinked )
- break;
// Move nodes on EDGE's
set< _Shrinker1D* >::iterator shr = eShri1D.begin();
moved = false;
for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
{
- moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,/*set3D=*/false );
+ moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
+ /*isCentroidal=*/isConcaveFace,
+ /*set3D=*/isConcaveFace);
}
if ( badNb < oldBadNb )
nbNoImpSteps = 0;
bool highQuality;
{
const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles();
- if ( hasTria != hasQuad )
- {
+ if ( hasTria != hasQuad ) {
highQuality = hasQuad;
}
- else
- {
+ else {
set<int> nbNodesSet;
SMDS_ElemIteratorPtr fIt = smDS->GetElements();
while ( fIt->more() && nbNodesSet.size() < 2 )
highQuality = ( *nbNodesSet.begin() == 4 );
}
}
- for ( int st = highQuality ? 8 : 3; st; --st )
+ if ( !highQuality && isConcaveFace )
+ fixBadFaces( F, helper ); // fix narrow faces by swaping diagonals
+ for ( int st = highQuality ? 10 : 3; st; --st )
{
dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
- nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,/*set3D=*/st==1 );
+ nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
+ /*isCentroidal=*/isConcaveFace,/*set3D=*/st==1 );
dumpFunctionEnd();
}
// Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh
- _SrinkShapeListener::ToClearSubMeshWithSolid( sm, data._solid );
+ VISCOUS_3D::ToClearSubWithMain( sm, data._solid );
- }// loop on FACES to srink mesh on
+ if ( !getMeshDS()->IsEmbeddedMode() )
+ // Log node movement
+ for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
+ {
+ SMESH_TNodeXYZ p ( nodesToSmooth[i]._node );
+ getMeshDS()->MoveNode( nodesToSmooth[i]._node, p.X(), p.Y(), p.Z() );
+ }
+
+ } // loop on FACES to srink mesh on
// Replace source nodes by target nodes in shrinked mesh edges
// return true;
}
+//================================================================================
+/*!
+ * \brief Try to fix triangles with high aspect ratio by swaping diagonals
+ */
+//================================================================================
+
+void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper)
+{
+ SMESH::Controls::AspectRatio qualifier;
+ SMESH::Controls::TSequenceOfXYZ points(3), points1(3), points2(3);
+ const double maxAspectRatio = 4.;
+
+ // find bad triangles
+
+ vector< const SMDS_MeshElement* > badTrias;
+ vector< double > badAspects;
+ SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( F );
+ SMDS_ElemIteratorPtr fIt = sm->GetElements();
+ while ( fIt->more() )
+ {
+ const SMDS_MeshElement * f = fIt->next();
+ if ( f->NbCornerNodes() != 3 ) continue;
+ for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = SMESH_TNodeXYZ( f->GetNode(iP));
+ double aspect = qualifier.GetValue( points );
+ if ( aspect > maxAspectRatio )
+ {
+ badTrias.push_back( f );
+ badAspects.push_back( aspect );
+ }
+ }
+ if ( badTrias.empty() )
+ return;
+
+ // find couples of faces to swap diagonal
+
+ typedef pair < const SMDS_MeshElement* , const SMDS_MeshElement* > T2Trias;
+ vector< T2Trias > triaCouples;
+
+ TIDSortedElemSet involvedFaces, emptySet;
+ for ( size_t iTia = 0; iTia < badTrias.size(); ++iTia )
+ {
+ T2Trias trias [3];
+ double aspRatio [3];
+ int i1, i2, i3;
+
+ involvedFaces.insert( badTrias[iTia] );
+ for ( int iP = 0; iP < 3; ++iP )
+ points(iP+1) = SMESH_TNodeXYZ( badTrias[iTia]->GetNode(iP));
+
+ // find triangles adjacent to badTrias[iTia] with better aspect ratio after diag-swaping
+ int bestCouple = -1;
+ for ( int iSide = 0; iSide < 3; ++iSide )
+ {
+ const SMDS_MeshNode* n1 = badTrias[iTia]->GetNode( iSide );
+ const SMDS_MeshNode* n2 = badTrias[iTia]->GetNode(( iSide+1 ) % 3 );
+ trias [iSide].first = badTrias[iTia];
+ trias [iSide].second = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, involvedFaces,
+ & i1, & i2 );
+ if ( ! trias[iSide].second || trias[iSide].second->NbCornerNodes() != 3 )
+ continue;
+
+ // aspect ratio of an adjacent tria
+ for ( int iP = 0; iP < 3; ++iP )
+ points2(iP+1) = SMESH_TNodeXYZ( trias[iSide].second->GetNode(iP));
+ double aspectInit = qualifier.GetValue( points2 );
+
+ // arrange nodes as after diag-swaping
+ if ( helper.WrapIndex( i1+1, 3 ) == i2 )
+ i3 = helper.WrapIndex( i1-1, 3 );
+ else
+ i3 = helper.WrapIndex( i1+1, 3 );
+ points1 = points;
+ points1( 1+ iSide ) = points2( 1+ i3 );
+ points2( 1+ i2 ) = points1( 1+ ( iSide+2 ) % 3 );
+
+ // aspect ratio after diag-swaping
+ aspRatio[ iSide ] = qualifier.GetValue( points1 ) + qualifier.GetValue( points2 );
+ if ( aspRatio[ iSide ] > aspectInit + badAspects[ iTia ] )
+ continue;
+
+ if ( bestCouple < 0 || aspRatio[ bestCouple ] > aspRatio[ iSide ] )
+ bestCouple = iSide;
+ }
+
+ if ( bestCouple >= 0 )
+ {
+ triaCouples.push_back( trias[bestCouple] );
+ involvedFaces.insert ( trias[bestCouple].second );
+ }
+ else
+ {
+ involvedFaces.erase( badTrias[iTia] );
+ }
+ }
+ if ( triaCouples.empty() )
+ return;
+
+ // swap diagonals
+
+ SMESH_MeshEditor editor( helper.GetMesh() );
+ dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID());
+ for ( size_t i = 0; i < triaCouples.size(); ++i )
+ {
+ dumpChangeNodes( triaCouples[i].first );
+ dumpChangeNodes( triaCouples[i].second );
+ editor.InverseDiag( triaCouples[i].first, triaCouples[i].second );
+ }
+ dumpFunctionEnd();
+
+ // just for debug dump resulting triangles
+ dumpFunction(SMESH_Comment("swapDiagonals_F")<<helper.GetSubShapeID());
+ for ( size_t i = 0; i < triaCouples.size(); ++i )
+ {
+ dumpChangeNodes( triaCouples[i].first );
+ dumpChangeNodes( triaCouples[i].second );
+ }
+}
+
//================================================================================
/*!
* \brief Move target node to it's final position on the FACE during shrinking
//================================================================================
/*!
- * \brief Perform laplacian smooth on the FACE
+ * \brief Perform smooth on the FACE
* \retval bool - true if the node has been moved
*/
//================================================================================
Handle(Geom_Surface)& surface,
SMESH_MesherHelper& helper,
const double refSign,
+ bool isCentroidal,
bool set3D)
{
const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() );
+ // get uv of surrounding nodes
+ vector<gp_XY> uv( _simplices.size() );
+ for ( size_t i = 0; i < _simplices.size(); ++i )
+ uv[i] = helper.GetNodeUV( face, _simplices[i]._nPrev, _node );
+
// compute new UV for the node
gp_XY newPos (0,0);
- for ( unsigned i = 0; i < _simplices.size(); ++i )
- newPos += helper.GetNodeUV( face, _simplices[i]._nPrev );
- newPos /= _simplices.size();
+ if ( isCentroidal && _simplices.size() > 3 )
+ {
+ // average centers of diagonals wieghted with their reciprocal lengths
+ if ( _simplices.size() == 4 )
+ {
+ double w1 = 1. / ( uv[2]-uv[0] ).SquareModulus();
+ double w2 = 1. / ( uv[3]-uv[1] ).SquareModulus();
+ newPos = ( w1 * ( uv[2]+uv[0] ) + w2 * ( uv[3]+uv[1] )) / ( w1+w2 ) / 2;
+ }
+ else
+ {
+ double sumWeight = 0;
+ int nb = _simplices.size() == 4 ? 2 : _simplices.size();
+ for ( int i = 0; i < nb; ++i )
+ {
+ int iFrom = i + 2;
+ int iTo = i + _simplices.size() - 1;
+ for ( int j = iFrom; j < iTo; ++j )
+ {
+ int i2 = SMESH_MesherHelper::WrapIndex( j, _simplices.size() );
+ double w = 1. / ( uv[i]-uv[i2] ).SquareModulus();
+ sumWeight += w;
+ newPos += w * ( uv[i]+uv[i2] );
+ }
+ }
+ newPos /= 2 * sumWeight;
+ }
+ }
+ else
+ {
+ // Laplacian smooth
+ isCentroidal = false;
+ for ( size_t i = 0; i < _simplices.size(); ++i )
+ newPos += uv[i];
+ newPos /= _simplices.size();
+ }
// count quality metrics (orientation) of triangles around the node
int nbOkBefore = 0;
gp_XY tgtUV = helper.GetNodeUV( face, _node );
for ( unsigned i = 0; i < _simplices.size(); ++i )
- nbOkBefore += _simplices[i].IsForward( tgtUV, face, helper, refSign );
+ nbOkBefore += _simplices[i].IsForward( tgtUV, _node, face, helper, refSign );
int nbOkAfter = 0;
for ( unsigned i = 0; i < _simplices.size(); ++i )
- nbOkAfter += _simplices[i].IsForward( newPos, face, helper, refSign );
+ nbOkAfter += _simplices[i].IsForward( newPos, _node, face, helper, refSign );
if ( nbOkAfter < nbOkBefore )
+ {
+ // if ( isCentroidal )
+ // return Smooth( badNb, surface, helper, refSign, !isCentroidal, set3D );
+ badNb += _simplices.size() - nbOkBefore;
return false;
+ }
SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( _node->GetPosition() );
pos->SetUParameter( newPos.X() );
GeomAdaptor_Curve aCurve(C, f,l);
const double totLen = GCPnts_AbscissaPoint::Length(aCurve, f, l);
- int nbExpectNodes = eSubMesh->NbNodes() - e->_nodes.size();
+ int nbExpectNodes = eSubMesh->NbNodes();
_initU .reserve( nbExpectNodes );
_normPar.reserve( nbExpectNodes );
_nodes .reserve( nbExpectNodes );
}
_done = false;
}
+
//================================================================================
/*!
* \brief Replace source nodes by target nodes in shrinked mesh edges
F = e2f->second.Oriented( TopAbs_FORWARD );
reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED );
if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED )
+ reverse = !reverse, F.Reverse();
+ if ( helper.IsReversedSubMesh( TopoDS::Face(F) ))
reverse = !reverse;
}
else
vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes;
vector< const SMDS_MeshNode*>& nn2 = ledges[j-dj2]->_nodes;
if ( isOnFace )
- for ( unsigned z = 1; z < nn1.size(); ++z )
+ 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 ( unsigned z = 1; z < nn1.size(); ++z )
+ for ( size_t z = 1; z < nn1.size(); ++z )
sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z]));
}
+
+ // Make edges
+ for ( int isFirst = 0; isFirst < 2; ++isFirst )
+ {
+ _LayerEdge* edge = isFirst ? ledges.front() : ledges.back();
+ if ( !edge->_sWOL.IsNull() && edge->_sWOL.ShapeType() == TopAbs_EDGE )
+ {
+ vector< const SMDS_MeshNode*>& nn = edge->_nodes;
+ if ( nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() )
+ continue;
+ helper.SetSubShape( edge->_sWOL );
+ helper.SetElementsOnShape( true );
+ for ( size_t z = 1; z < nn.size(); ++z )
+ helper.AddEdge( nn[z-1], nn[z] );
+ }
+ }
}
}
