00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include "StdMeshers_ViscousLayers.hxx"
00025
00026 #include "SMDS_EdgePosition.hxx"
00027 #include "SMDS_FaceOfNodes.hxx"
00028 #include "SMDS_FacePosition.hxx"
00029 #include "SMDS_MeshNode.hxx"
00030 #include "SMDS_SetIterator.hxx"
00031 #include "SMESHDS_Group.hxx"
00032 #include "SMESHDS_Hypothesis.hxx"
00033 #include "SMESH_Algo.hxx"
00034 #include "SMESH_ComputeError.hxx"
00035 #include "SMESH_Gen.hxx"
00036 #include "SMESH_Group.hxx"
00037 #include "SMESH_Mesh.hxx"
00038 #include "SMESH_MesherHelper.hxx"
00039 #include "SMESH_subMesh.hxx"
00040 #include "SMESH_subMeshEventListener.hxx"
00041 #include "SMESH_ProxyMesh.hxx"
00042
00043 #include "utilities.h"
00044
00045 #include <BRep_Tool.hxx>
00046 #include <Bnd_B2d.hxx>
00047 #include <Bnd_B3d.hxx>
00048 #include <ElCLib.hxx>
00049 #include <GCPnts_AbscissaPoint.hxx>
00050 #include <Geom2d_Circle.hxx>
00051 #include <Geom2d_Line.hxx>
00052 #include <Geom2d_TrimmedCurve.hxx>
00053 #include <GeomAdaptor_Curve.hxx>
00054 #include <Geom_Circle.hxx>
00055 #include <Geom_Curve.hxx>
00056 #include <Geom_Line.hxx>
00057 #include <Geom_TrimmedCurve.hxx>
00058 #include <Precision.hxx>
00059 #include <TopExp.hxx>
00060 #include <TopExp_Explorer.hxx>
00061 #include <TopTools_IndexedMapOfShape.hxx>
00062 #include <TopTools_MapOfShape.hxx>
00063 #include <TopoDS.hxx>
00064 #include <TopoDS_Edge.hxx>
00065 #include <TopoDS_Face.hxx>
00066 #include <TopoDS_Vertex.hxx>
00067 #include <gp_Ax1.hxx>
00068 #include <gp_Vec.hxx>
00069 #include <gp_XY.hxx>
00070 #include <gp_XYZ.hxx>
00071
00072 #include <list>
00073 #include <string>
00074 #include <math.h>
00075 #include <limits>
00076
00077
00078
00079 using namespace std;
00080
00081
00082 namespace VISCOUS
00083 {
00084 typedef int TGeomID;
00085
00086 enum UIndex { U_TGT = 1, U_SRC, LEN_TGT };
00087
00092 struct _MeshOfSolid : public SMESH_ProxyMesh,
00093 public SMESH_subMeshEventListenerData
00094 {
00095 bool _n2nMapComputed;
00096
00097 _MeshOfSolid( SMESH_Mesh* mesh)
00098 :SMESH_subMeshEventListenerData( true),_n2nMapComputed(false)
00099 {
00100 SMESH_ProxyMesh::setMesh( *mesh );
00101 }
00102
00103
00104 SMESH_ProxyMesh::SubMesh* getFaceSubM(const TopoDS_Face& F, bool create=false)
00105 {
00106 TGeomID i = SMESH_ProxyMesh::shapeIndex(F);
00107 return create ? SMESH_ProxyMesh::getProxySubMesh(i) : findProxySubMesh(i);
00108 }
00109 void setNode2Node(const SMDS_MeshNode* srcNode,
00110 const SMDS_MeshNode* proxyNode,
00111 const SMESH_ProxyMesh::SubMesh* subMesh)
00112 {
00113 SMESH_ProxyMesh::setNode2Node( srcNode,proxyNode,subMesh);
00114 }
00115 };
00116
00121 class _SrinkShapeListener : SMESH_subMeshEventListener
00122 {
00123 _SrinkShapeListener(): SMESH_subMeshEventListener(false) {}
00124 static SMESH_subMeshEventListener* Get() { static _SrinkShapeListener l; return &l; }
00125 public:
00126 virtual void ProcessEvent(const int event,
00127 const int eventType,
00128 SMESH_subMesh* solidSM,
00129 SMESH_subMeshEventListenerData* data,
00130 const SMESH_Hypothesis* hyp)
00131 {
00132 if ( SMESH_subMesh::COMPUTE_EVENT == eventType && solidSM->IsEmpty() && data )
00133 {
00134 SMESH_subMeshEventListener::ProcessEvent(event,eventType,solidSM,data,hyp);
00135 }
00136 }
00137 static void ToClearSubMeshWithSolid( SMESH_subMesh* sm,
00138 const TopoDS_Shape& solid)
00139 {
00140 SMESH_subMesh* solidSM = sm->GetFather()->GetSubMesh( solid );
00141 SMESH_subMeshEventListenerData* data = solidSM->GetEventListenerData( Get());
00142 if ( data )
00143 {
00144 if ( find( data->mySubMeshes.begin(), data->mySubMeshes.end(), sm ) ==
00145 data->mySubMeshes.end())
00146 data->mySubMeshes.push_back( sm );
00147 }
00148 else
00149 {
00150 data = SMESH_subMeshEventListenerData::MakeData( sm );
00151 sm->SetEventListener( Get(), data, solidSM );
00152 }
00153 }
00154 };
00155
00161 class _ViscousListener : SMESH_subMeshEventListener
00162 {
00163 _ViscousListener(): SMESH_subMeshEventListener(false) {}
00164 static SMESH_subMeshEventListener* Get() { static _ViscousListener l; return &l; }
00165 public:
00166 virtual void ProcessEvent(const int event,
00167 const int eventType,
00168 SMESH_subMesh* subMesh,
00169 SMESH_subMeshEventListenerData* data,
00170 const SMESH_Hypothesis* hyp)
00171 {
00172 if ( SMESH_subMesh::COMPUTE_EVENT == eventType )
00173 {
00174
00175 subMesh->DeleteEventListener( this );
00176 }
00177 }
00178
00179 static _MeshOfSolid* GetSolidMesh(SMESH_Mesh* mesh,
00180 const TopoDS_Shape& solid,
00181 bool toCreate=false)
00182 {
00183 if ( !mesh ) return 0;
00184 SMESH_subMesh* sm = mesh->GetSubMesh(solid);
00185 _MeshOfSolid* data = (_MeshOfSolid*) sm->GetEventListenerData( Get() );
00186 if ( !data && toCreate )
00187 {
00188 data = new _MeshOfSolid(mesh);
00189 data->mySubMeshes.push_back( sm );
00190 sm->SetEventListener( Get(), data, sm );
00191 }
00192 return data;
00193 }
00194
00195 static void RemoveSolidMesh(SMESH_Mesh* mesh, const TopoDS_Shape& solid)
00196 {
00197 mesh->GetSubMesh(solid)->DeleteEventListener( _ViscousListener::Get() );
00198 }
00199 };
00200
00201
00208 struct _Simplex
00209 {
00210 const SMDS_MeshNode *_nPrev, *_nNext;
00211 _Simplex(const SMDS_MeshNode* nPrev=0, const SMDS_MeshNode* nNext=0)
00212 : _nPrev(nPrev), _nNext(nNext) {}
00213 bool IsForward(const SMDS_MeshNode* nSrc, const gp_XYZ* pntTgt) const
00214 {
00215 const double M[3][3] =
00216 {{ _nNext->X() - nSrc->X(), _nNext->Y() - nSrc->Y(), _nNext->Z() - nSrc->Z() },
00217 { pntTgt->X() - nSrc->X(), pntTgt->Y() - nSrc->Y(), pntTgt->Z() - nSrc->Z() },
00218 { _nPrev->X() - nSrc->X(), _nPrev->Y() - nSrc->Y(), _nPrev->Z() - nSrc->Z() }};
00219 double determinant = ( + M[0][0]*M[1][1]*M[2][2]
00220 + M[0][1]*M[1][2]*M[2][0]
00221 + M[0][2]*M[1][0]*M[2][1]
00222 - M[0][0]*M[1][2]*M[2][1]
00223 - M[0][1]*M[1][0]*M[2][2]
00224 - M[0][2]*M[1][1]*M[2][0]);
00225 return determinant > 1e-100;
00226 }
00227 bool IsForward(const gp_XY& tgtUV,
00228 const TopoDS_Face& face,
00229 SMESH_MesherHelper& helper,
00230 const double refSign) const
00231 {
00232 gp_XY prevUV = helper.GetNodeUV( face, _nPrev );
00233 gp_XY nextUV = helper.GetNodeUV( face, _nNext );
00234 gp_Vec2d v1( tgtUV, prevUV ), v2( tgtUV, nextUV );
00235 double d = v1 ^ v2;
00236 return d*refSign > 1e-100;
00237 }
00238 };
00239
00243 struct _Curvature
00244 {
00245 double _r;
00246 double _k;
00247 public:
00248 static _Curvature* New( double avgNormProj, double avgDist )
00249 {
00250 _Curvature* c = 0;
00251 if ( fabs( avgNormProj / avgDist ) > 1./200 )
00252 {
00253 c = new _Curvature;
00254 c->_r = avgDist * avgDist / avgNormProj;
00255 c->_k = avgDist * avgDist / c->_r / c->_r;
00256 c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 );
00257 }
00258 return c;
00259 }
00260 double lenDelta(double len) const { return _k * ( _r + len ); }
00261 };
00262 struct _LayerEdge;
00263
00267 struct _2NearEdges
00268 {
00269
00270 const SMDS_MeshNode* _nodes[2];
00271
00272
00273 double _wgt[2];
00274 _LayerEdge* _edges[2];
00275
00276
00277 gp_XYZ* _plnNorm;
00278
00279 _2NearEdges() { _nodes[0]=_nodes[1]=0; _plnNorm = 0; }
00280 void reverse() {
00281 std::swap( _nodes[0], _nodes[1] );
00282 std::swap( _wgt[0], _wgt[1] );
00283 }
00284 };
00285
00290 struct _LayerEdge
00291 {
00292 vector< const SMDS_MeshNode*> _nodes;
00293
00294 gp_XYZ _normal;
00295 vector<gp_XYZ> _pos;
00296 double _len;
00297 double _cosin;
00298 double _lenFactor;
00299
00300
00301 TopoDS_Shape _sWOL;
00302
00303
00304 vector<_Simplex> _simplices;
00305
00306 _2NearEdges* _2neibors;
00307
00308 _Curvature* _curvature;
00309
00310
00311 void SetNewLength( double len, SMESH_MesherHelper& helper );
00312 bool SetNewLength2d( Handle(Geom_Surface)& surface,
00313 const TopoDS_Face& F,
00314 SMESH_MesherHelper& helper );
00315 void SetDataByNeighbors( const SMDS_MeshNode* n1,
00316 const SMDS_MeshNode* n2,
00317 SMESH_MesherHelper& helper);
00318 void InvalidateStep( int curStep );
00319 bool Smooth(int& badNb);
00320 bool SmoothOnEdge(Handle(Geom_Surface)& surface,
00321 const TopoDS_Face& F,
00322 SMESH_MesherHelper& helper);
00323 bool FindIntersection( SMESH_ElementSearcher& searcher,
00324 double & distance,
00325 const double& epsilon,
00326 const SMDS_MeshElement** face = 0);
00327 bool SegTriaInter( const gp_Ax1& lastSegment,
00328 const SMDS_MeshNode* n0,
00329 const SMDS_MeshNode* n1,
00330 const SMDS_MeshNode* n2,
00331 double& dist,
00332 const double& epsilon) const;
00333 gp_Ax1 LastSegment(double& segLen) const;
00334 bool IsOnEdge() const { return _2neibors; }
00335 void Copy( _LayerEdge& other, SMESH_MesherHelper& helper );
00336 void SetCosin( double cosin );
00337 };
00338 struct _LayerEdgeCmp
00339 {
00340 bool operator () (const _LayerEdge* e1, const _LayerEdge* e2) const
00341 {
00342 const bool cmpNodes = ( e1 && e2 && e1->_nodes.size() && e2->_nodes.size() );
00343 return cmpNodes ? ( e1->_nodes[0]->GetID() < e2->_nodes[0]->GetID()) : ( e1 < e2 );
00344 }
00345 };
00346
00347
00348 typedef map< const SMDS_MeshNode*, _LayerEdge*, TIDCompare > TNode2Edge;
00349
00350
00354 struct _SolidData
00355 {
00356 TopoDS_Shape _solid;
00357 const StdMeshers_ViscousLayers* _hyp;
00358 _MeshOfSolid* _proxyMesh;
00359 set<TGeomID> _reversedFaceIds;
00360
00361 double _stepSize, _stepSizeCoeff;
00362 const SMDS_MeshNode* _stepSizeNodes[2];
00363
00364 TNode2Edge _n2eMap;
00365
00366
00367 vector< _LayerEdge* > _edges;
00368
00369
00370
00371
00372
00373 map< TGeomID, TopoDS_Shape > _shrinkShape2Shape;
00374
00375
00376 set< TGeomID > _noShrinkFaces;
00377
00378
00379 map< TGeomID,Handle(Geom_Curve)> _edge2curve;
00380
00381
00382 vector< int > _endEdgeToSmooth;
00383
00384 double _epsilon;
00385
00386 int _index;
00387
00388 _SolidData(const TopoDS_Shape& s=TopoDS_Shape(),
00389 const StdMeshers_ViscousLayers* h=0,
00390 _MeshOfSolid* m=0) :_solid(s), _hyp(h), _proxyMesh(m) {}
00391 ~_SolidData();
00392
00393 Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge& E,
00394 const int iFrom,
00395 const int iTo,
00396 Handle(Geom_Surface)& surface,
00397 const TopoDS_Face& F,
00398 SMESH_MesherHelper& helper);
00399 };
00400
00404 struct _SmoothNode
00405 {
00406 const SMDS_MeshNode* _node;
00407
00408 vector<_Simplex> _simplices;
00409
00410 bool Smooth(int& badNb,
00411 Handle(Geom_Surface)& surface,
00412 SMESH_MesherHelper& helper,
00413 const double refSign,
00414 bool set3D);
00415 };
00416
00420 class _ViscousBuilder
00421 {
00422 public:
00423 _ViscousBuilder();
00424
00425 SMESH_ComputeErrorPtr Compute(SMESH_Mesh& mesh,
00426 const TopoDS_Shape& shape);
00427
00428
00429 void RestoreListeners();
00430
00431
00432 bool MakeN2NMap( _MeshOfSolid* pm );
00433
00434 private:
00435
00436 bool findSolidsWithLayers();
00437 bool findFacesWithLayers();
00438 bool makeLayer(_SolidData& data);
00439 bool setEdgeData(_LayerEdge& edge, const set<TGeomID>& subIds,
00440 SMESH_MesherHelper& helper, _SolidData& data);
00441 bool findNeiborsOnEdge(const _LayerEdge* edge,
00442 const SMDS_MeshNode*& n1,
00443 const SMDS_MeshNode*& n2,
00444 _SolidData& data);
00445 void getSimplices( const SMDS_MeshNode* node, vector<_Simplex>& simplices,
00446 const set<TGeomID>& ingnoreShapes,
00447 const _SolidData* dataToCheckOri = 0);
00448 bool sortEdges( _SolidData& data,
00449 vector< vector<_LayerEdge*> >& edgesByGeom);
00450 void limitStepSize( _SolidData& data,
00451 const SMDS_MeshElement* face,
00452 const double cosin);
00453 void limitStepSize( _SolidData& data, const double minSize);
00454 bool inflate(_SolidData& data);
00455 bool smoothAndCheck(_SolidData& data, const int nbSteps, double & distToIntersection);
00456 bool smoothAnalyticEdge( _SolidData& data,
00457 const int iFrom,
00458 const int iTo,
00459 Handle(Geom_Surface)& surface,
00460 const TopoDS_Face& F,
00461 SMESH_MesherHelper& helper);
00462 bool updateNormals( _SolidData& data, SMESH_MesherHelper& helper );
00463 bool refine(_SolidData& data);
00464 bool shrink();
00465 bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F,
00466 SMESH_MesherHelper& helper,
00467 const SMESHDS_SubMesh* faceSubMesh );
00468 bool addBoundaryElements();
00469
00470 bool error( const string& text, int solidID=-1 );
00471 SMESHDS_Mesh* getMeshDS() { return _mesh->GetMeshDS(); }
00472
00473
00474 void makeGroupOfLE();
00475
00476 SMESH_Mesh* _mesh;
00477 SMESH_ComputeErrorPtr _error;
00478
00479 vector< _SolidData > _sdVec;
00480 set<TGeomID> _ignoreShapeIds;
00481 int _tmpFaceID;
00482 };
00483
00487 class _Shrinker1D
00488 {
00489 vector<double> _initU;
00490 vector<double> _normPar;
00491 vector<const SMDS_MeshNode*> _nodes;
00492 const _LayerEdge* _edges[2];
00493 bool _done;
00494 public:
00495 void AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper );
00496 void Compute(bool set3D, SMESH_MesherHelper& helper);
00497 void RestoreParams();
00498 void SwapSrcTgtNodes(SMESHDS_Mesh* mesh);
00499 };
00500
00506 struct TmpMeshFace : public SMDS_MeshElement
00507 {
00508 vector<const SMDS_MeshNode* > _nn;
00509 TmpMeshFace( const vector<const SMDS_MeshNode*>& nodes, int id):
00510 SMDS_MeshElement(id), _nn(nodes) {}
00511 virtual const SMDS_MeshNode* GetNode(const int ind) const { return _nn[ind]; }
00512 virtual SMDSAbs_ElementType GetType() const { return SMDSAbs_Face; }
00513 virtual vtkIdType GetVtkType() const { return -1; }
00514 virtual SMDSAbs_EntityType GetEntityType() const { return SMDSEntity_Last; }
00515 virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const
00516 { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));}
00517 };
00518
00522 struct TmpMeshFaceOnEdge : public TmpMeshFace
00523 {
00524 _LayerEdge *_le1, *_le2;
00525 TmpMeshFaceOnEdge( _LayerEdge* le1, _LayerEdge* le2, int ID ):
00526 TmpMeshFace( vector<const SMDS_MeshNode*>(4), ID ), _le1(le1), _le2(le2)
00527 {
00528 _nn[0]=_le1->_nodes[0];
00529 _nn[1]=_le1->_nodes.back();
00530 _nn[2]=_le2->_nodes.back();
00531 _nn[3]=_le2->_nodes[0];
00532 }
00533 };
00534 }
00535
00536
00537
00538
00539 StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen)
00540 :SMESH_Hypothesis(hypId, studyId, gen),
00541 _nbLayers(1), _thickness(1), _stretchFactor(1)
00542 {
00543 _name = StdMeshers_ViscousLayers::GetHypType();
00544 _param_algo_dim = -3;
00545 }
00546 void StdMeshers_ViscousLayers::SetIgnoreFaces(const std::vector<int>& faceIds)
00547 {
00548 if ( faceIds != _ignoreFaceIds )
00549 _ignoreFaceIds = faceIds, NotifySubMeshesHypothesisModification();
00550 }
00551 void StdMeshers_ViscousLayers::SetTotalThickness(double thickness)
00552 {
00553 if ( thickness != _thickness )
00554 _thickness = thickness, NotifySubMeshesHypothesisModification();
00555 }
00556 void StdMeshers_ViscousLayers::SetNumberLayers(int nb)
00557 {
00558 if ( _nbLayers != nb )
00559 _nbLayers = nb, NotifySubMeshesHypothesisModification();
00560 }
00561 void StdMeshers_ViscousLayers::SetStretchFactor(double factor)
00562 {
00563 if ( _stretchFactor != factor )
00564 _stretchFactor = factor, NotifySubMeshesHypothesisModification();
00565 }
00566 SMESH_ProxyMesh::Ptr
00567 StdMeshers_ViscousLayers::Compute(SMESH_Mesh& theMesh,
00568 const TopoDS_Shape& theShape,
00569 const bool toMakeN2NMap) const
00570 {
00571 using namespace VISCOUS;
00572 _ViscousBuilder bulder;
00573 SMESH_ComputeErrorPtr err = bulder.Compute( theMesh, theShape );
00574 if ( err && !err->IsOK() )
00575 return SMESH_ProxyMesh::Ptr();
00576
00577 vector<SMESH_ProxyMesh::Ptr> components;
00578 TopExp_Explorer exp( theShape, TopAbs_SOLID );
00579 for ( ; exp.More(); exp.Next() )
00580 {
00581 if ( _MeshOfSolid* pm =
00582 _ViscousListener::GetSolidMesh( &theMesh, exp.Current(), false))
00583 {
00584 if ( toMakeN2NMap && !pm->_n2nMapComputed )
00585 if ( !bulder.MakeN2NMap( pm ))
00586 return SMESH_ProxyMesh::Ptr();
00587 components.push_back( SMESH_ProxyMesh::Ptr( pm ));
00588 pm->myIsDeletable = false;
00589 }
00590 _ViscousListener::RemoveSolidMesh ( &theMesh, exp.Current() );
00591 }
00592 switch ( components.size() )
00593 {
00594 case 0: break;
00595
00596 case 1: return components[0];
00597
00598 default: return SMESH_ProxyMesh::Ptr( new SMESH_ProxyMesh( components ));
00599 }
00600 return SMESH_ProxyMesh::Ptr();
00601 }
00602 std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save)
00603 {
00604 save << " " << _nbLayers
00605 << " " << _thickness
00606 << " " << _stretchFactor
00607 << " " << _ignoreFaceIds.size();
00608 for ( unsigned i = 0; i < _ignoreFaceIds.size(); ++i )
00609 save << " " << _ignoreFaceIds[i];
00610 return save;
00611 }
00612 std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load)
00613 {
00614 int nbFaces, faceID;
00615 load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces;
00616 while ( _ignoreFaceIds.size() < nbFaces && load >> faceID )
00617 _ignoreFaceIds.push_back( faceID );
00618 return load;
00619 }
00620 bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh* theMesh,
00621 const TopoDS_Shape& theShape)
00622 {
00623
00624 return false;
00625 }
00626
00627
00628
00629 namespace
00630 {
00631 gp_XYZ getEdgeDir( const TopoDS_Edge& E, const TopoDS_Vertex& fromV )
00632 {
00633 gp_Vec dir;
00634 double f,l;
00635 Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
00636 gp_Pnt p = BRep_Tool::Pnt( fromV );
00637 double distF = p.SquareDistance( c->Value( f ));
00638 double distL = p.SquareDistance( c->Value( l ));
00639 c->D1(( distF < distL ? f : l), p, dir );
00640 if ( distL < distF ) dir.Reverse();
00641 return dir.XYZ();
00642 }
00643
00644 gp_XYZ getEdgeDir( const TopoDS_Edge& E, const SMDS_MeshNode* atNode,
00645 SMESH_MesherHelper& helper)
00646 {
00647 gp_Vec dir;
00648 double f,l; gp_Pnt p;
00649 Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
00650 double u = helper.GetNodeU( E, atNode );
00651 c->D1( u, p, dir );
00652 return dir.XYZ();
00653 }
00654
00655 gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE,
00656 const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok)
00657 {
00658 gp_XY uv = helper.GetNodeUV( F, node, 0, &ok );
00659 Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
00660 gp_Pnt p; gp_Vec du, dv, norm;
00661 surface->D1( uv.X(),uv.Y(), p, du,dv );
00662 norm = du ^ dv;
00663
00664 double f,l;
00665 Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l );
00666 double u = helper.GetNodeU( fromE, node, 0, &ok );
00667 c->D1( u, p, du );
00668 TopAbs_Orientation o = helper.GetSubShapeOri( F.Oriented(TopAbs_FORWARD), fromE);
00669 if ( o == TopAbs_REVERSED )
00670 du.Reverse();
00671
00672 gp_Vec dir = norm ^ du;
00673
00674 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX &&
00675 helper.IsClosedEdge( fromE ))
00676 {
00677 if ( fabs(u-f) < fabs(u-l )) c->D1( l, p, dv );
00678 else c->D1( f, p, dv );
00679 if ( o == TopAbs_REVERSED )
00680 dv.Reverse();
00681 gp_Vec dir2 = norm ^ dv;
00682 dir = dir.Normalized() + dir2.Normalized();
00683 }
00684 return dir.XYZ();
00685 }
00686
00687 gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
00688 const SMDS_MeshNode* node, SMESH_MesherHelper& helper,
00689 bool& ok, double* cosin=0)
00690 {
00691 double f,l; TopLoc_Location loc;
00692 vector< TopoDS_Edge > edges;
00693 PShapeIteratorPtr eIt = helper.GetAncestors( fromV, *helper.GetMesh(), TopAbs_EDGE);
00694 while ( eIt->more())
00695 {
00696 const TopoDS_Edge* e = static_cast<const TopoDS_Edge*>( eIt->next() );
00697 if ( helper.IsSubShape( *e, F ) && BRep_Tool::Curve( *e, loc,f,l))
00698 edges.push_back( *e );
00699 }
00700 gp_XYZ dir(0,0,0);
00701 if ( !( ok = ( edges.size() > 0 ))) return dir;
00702
00703 gp_Vec edgeDir;
00704 for ( unsigned i = 0; i < edges.size(); ++i )
00705 {
00706 edgeDir = getEdgeDir( edges[i], fromV );
00707 double size2 = edgeDir.SquareMagnitude();
00708 if ( size2 > numeric_limits<double>::min() )
00709 edgeDir /= sqrt( size2 );
00710 else
00711 ok = false;
00712 dir += edgeDir.XYZ();
00713 }
00714 gp_XYZ fromEdgeDir = getFaceDir( F, edges[0], node, helper, ok );
00715 if ( edges.size() == 1 || dir.SquareModulus() < 1e-10)
00716 dir = fromEdgeDir;
00717 else if ( dir * fromEdgeDir < 0 )
00718 dir *= -1;
00719 if ( ok )
00720 {
00721
00722 if ( cosin ) {
00723 double angle = edgeDir.Angle( dir );
00724 *cosin = cos( angle );
00725 }
00726 }
00727 return dir;
00728 }
00729
00730
00731 #ifdef __myDEBUG
00732 ofstream* py;
00733 struct PyDump {
00734 PyDump() {
00735 const char* fname = "/tmp/viscous.py";
00736 cout << "execfile('"<<fname<<"')"<<endl;
00737 py = new ofstream(fname);
00738 *py << "from smesh import *" << endl
00739 << "meshSO = GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
00740 << "mesh = Mesh( meshSO.GetObject()._narrow( SMESH.SMESH_Mesh ))"<<endl;
00741 }
00742 ~PyDump() {
00743 *py << "mesh.MakeGroup('Prisms of viscous layers',VOLUME,FT_ElemGeomType,'=',Geom_PENTA)"
00744 <<endl; delete py; py=0;
00745 }
00746 };
00747 #define dumpFunction(f) { _dumpFunction(f, __LINE__);}
00748 #define dumpMove(n) { _dumpMove(n, __LINE__);}
00749 #define dumpCmd(txt) { _dumpCmd(txt, __LINE__);}
00750 void _dumpFunction(const string& fun, int ln)
00751 { *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl;}
00752 void _dumpMove(const SMDS_MeshNode* n, int ln)
00753 { *py<< " mesh.MoveNode( "<<n->GetID()<< ", "<< n->X()
00754 << ", "<<n->Y()<<", "<< n->Z()<< ")\t\t # "<< ln <<endl; }
00755 void _dumpCmd(const string& txt, int ln)
00756 { *py<< " "<<txt<<" # "<< ln <<endl; }
00757 void dumpFunctionEnd()
00758 { *py<< " return"<< endl; }
00759 #else
00760 struct PyDump { PyDump() {} };
00761 void dumpFunction(const string& fun ){}
00762 void dumpFunctionEnd() {}
00763 void dumpMove(const SMDS_MeshNode* n ){}
00764 void dumpCmd(const string& txt){}
00765 #endif
00766 }
00767
00768 using namespace VISCOUS;
00769
00770
00774
00775
00776 _ViscousBuilder::_ViscousBuilder()
00777 {
00778 _error = SMESH_ComputeError::New(COMPERR_OK);
00779 _tmpFaceID = 0;
00780 }
00781
00782
00786
00787
00788 bool _ViscousBuilder::error(const string& text, int solidId )
00789 {
00790 _error->myName = COMPERR_ALGO_FAILED;
00791 _error->myComment = string("Viscous layers builder: ") + text;
00792 if ( _mesh )
00793 {
00794 SMESH_subMesh* sm = _mesh->GetSubMeshContaining( solidId );
00795 if ( !sm && !_sdVec.empty() )
00796 sm = _mesh->GetSubMeshContaining( _sdVec[0]._index );
00797 if ( sm && sm->GetSubShape().ShapeType() == TopAbs_SOLID )
00798 {
00799 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
00800 if ( smError && smError->myAlgo )
00801 _error->myAlgo = smError->myAlgo;
00802 smError = _error;
00803 }
00804 }
00805 makeGroupOfLE();
00806
00807 return false;
00808 }
00809
00810
00815
00816
00817 void _ViscousBuilder::RestoreListeners()
00818 {
00819
00820 }
00821
00822
00826
00827
00828 bool _ViscousBuilder::MakeN2NMap( _MeshOfSolid* pm )
00829 {
00830 SMESH_subMesh* solidSM = pm->mySubMeshes.front();
00831 TopExp_Explorer fExp( solidSM->GetSubShape(), TopAbs_FACE );
00832 for ( ; fExp.More(); fExp.Next() )
00833 {
00834 SMESHDS_SubMesh* srcSmDS = pm->GetMeshDS()->MeshElements( fExp.Current() );
00835 const SMESH_ProxyMesh::SubMesh* prxSmDS = pm->GetProxySubMesh( fExp.Current() );
00836
00837 if ( !srcSmDS || !prxSmDS || !srcSmDS->NbElements() || !prxSmDS->NbElements() )
00838 continue;
00839 if ( srcSmDS->GetElements()->next() == prxSmDS->GetElements()->next())
00840 continue;
00841
00842 if ( srcSmDS->NbElements() != prxSmDS->NbElements() )
00843 return error( "Different nb elements in a source and a proxy sub-mesh", solidSM->GetId());
00844
00845 SMDS_ElemIteratorPtr srcIt = srcSmDS->GetElements();
00846 SMDS_ElemIteratorPtr prxIt = prxSmDS->GetElements();
00847 while( prxIt->more() )
00848 {
00849 const SMDS_MeshElement* fSrc = srcIt->next();
00850 const SMDS_MeshElement* fPrx = prxIt->next();
00851 if ( fSrc->NbNodes() != fPrx->NbNodes())
00852 return error( "Different elements in a source and a proxy sub-mesh", solidSM->GetId());
00853 for ( int i = 0 ; i < fPrx->NbNodes(); ++i )
00854 pm->setNode2Node( fSrc->GetNode(i), fPrx->GetNode(i), prxSmDS );
00855 }
00856 }
00857 pm->_n2nMapComputed = true;
00858 return true;
00859 }
00860
00861
00865
00866
00867 SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh,
00868 const TopoDS_Shape& theShape)
00869 {
00870
00871
00872 _mesh = & theMesh;
00873
00874
00875 TopExp_Explorer exp( theShape, TopAbs_SOLID );
00876 if ( !exp.More() )
00877 return error("No SOLID's in theShape"), _error;
00878
00879 if ( _ViscousListener::GetSolidMesh( _mesh, exp.Current(), false))
00880 return SMESH_ComputeErrorPtr();
00881
00882 PyDump debugDump;
00883
00884
00885 if ( !findSolidsWithLayers())
00886 return _error;
00887
00888 if ( !findFacesWithLayers() )
00889 return _error;
00890
00891 for ( unsigned i = 0; i < _sdVec.size(); ++i )
00892 {
00893 if ( ! makeLayer(_sdVec[i]) )
00894 return _error;
00895
00896 if ( ! inflate(_sdVec[i]) )
00897 return _error;
00898
00899 if ( ! refine(_sdVec[i]) )
00900 return _error;
00901 }
00902 if ( !shrink() )
00903 return _error;
00904
00905 addBoundaryElements();
00906
00907 makeGroupOfLE();
00908
00909 return _error;
00910 }
00911
00912
00916
00917
00918 bool _ViscousBuilder::findSolidsWithLayers()
00919 {
00920
00921 TopTools_IndexedMapOfShape allSolids;
00922 TopExp::MapShapes( _mesh->GetShapeToMesh(), TopAbs_SOLID, allSolids );
00923 _sdVec.reserve( allSolids.Extent());
00924
00925 SMESH_Gen* gen = _mesh->GetGen();
00926 for ( int i = 1; i <= allSolids.Extent(); ++i )
00927 {
00928
00929 SMESH_Algo* algo = gen->GetAlgo( *_mesh, allSolids(i) );
00930 if ( !algo ) continue;
00931
00932 const list <const SMESHDS_Hypothesis *> & allHyps =
00933 algo->GetUsedHypothesis(*_mesh, allSolids(i), false);
00934 list< const SMESHDS_Hypothesis *>::const_iterator hyp = allHyps.begin();
00935 const StdMeshers_ViscousLayers* viscHyp = 0;
00936 for ( ; hyp != allHyps.end() && !viscHyp; ++hyp )
00937 viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp );
00938 if ( viscHyp )
00939 {
00940 _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
00941 allSolids(i),
00942 true);
00943 _sdVec.push_back( _SolidData( allSolids(i), viscHyp, proxyMesh ));
00944 _sdVec.back()._index = getMeshDS()->ShapeToIndex( allSolids(i));
00945 }
00946 }
00947 if ( _sdVec.empty() )
00948 return error
00949 ( SMESH_Comment(StdMeshers_ViscousLayers::GetHypType()) << " hypothesis not found",0);
00950
00951 return true;
00952 }
00953
00954
00958
00959
00960 bool _ViscousBuilder::findFacesWithLayers()
00961 {
00962
00963 vector<TopoDS_Shape> ignoreFaces;
00964 for ( unsigned i = 0; i < _sdVec.size(); ++i )
00965 {
00966 vector<TGeomID> ids = _sdVec[i]._hyp->GetIgnoreFaces();
00967 for ( unsigned i = 0; i < ids.size(); ++i )
00968 {
00969 const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[i] );
00970 if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
00971 {
00972 _ignoreShapeIds.insert( ids[i] );
00973 ignoreFaces.push_back( s );
00974 }
00975 }
00976 }
00977
00978
00979 SMESH_MesherHelper helper( *_mesh );
00980 TopExp_Explorer exp;
00981 for ( unsigned i = 0; i < _sdVec.size(); ++i )
00982 {
00983 exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
00984 for ( ; exp.More(); exp.Next() )
00985 {
00986 TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
00987 if ( helper.NbAncestors( exp.Current(), *_mesh, TopAbs_SOLID ) > 1 )
00988 {
00989 _ignoreShapeIds.insert( faceInd );
00990 ignoreFaces.push_back( exp.Current() );
00991 if ( SMESH_Algo::IsReversedSubMesh( TopoDS::Face( exp.Current() ), getMeshDS()))
00992 _sdVec[i]._reversedFaceIds.insert( faceInd );
00993 }
00994 }
00995 }
00996
00997
00998 TopTools_IndexedMapOfShape shapes;
00999 for ( unsigned i = 0; i < _sdVec.size(); ++i )
01000 {
01001 shapes.Clear();
01002 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_EDGE, shapes);
01003 for ( int iE = 1; iE <= shapes.Extent(); ++iE )
01004 {
01005 const TopoDS_Shape& edge = shapes(iE);
01006
01007 TopoDS_Shape FF[2];
01008 PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE);
01009 while ( fIt->more())
01010 {
01011 const TopoDS_Shape* f = fIt->next();
01012 if ( helper.IsSubShape( *f, _sdVec[i]._solid))
01013 FF[ int( !FF[0].IsNull()) ] = *f;
01014 }
01015 if( FF[1].IsNull() ) continue;
01016
01017 int ignore[2];
01018 for ( int j = 0; j < 2; ++j )
01019 ignore[j] = _ignoreShapeIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
01020 if ( ignore[0] == ignore[1] ) continue;
01021 TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ];
01022
01023 TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
01024 _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
01025 }
01026 }
01027
01028
01029 set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D");
01030 for ( unsigned i = 0; i < _sdVec.size(); ++i )
01031 {
01032 TopTools_MapOfShape noShrinkVertices;
01033 map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
01034 for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
01035 {
01036 const TopoDS_Shape& fWOL = e2f->second;
01037 TGeomID edgeID = e2f->first;
01038 bool notShrinkFace = false;
01039 PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
01040 while ( soIt->more())
01041 {
01042 const TopoDS_Shape* solid = soIt->next();
01043 if ( _sdVec[i]._solid.IsSame( *solid )) continue;
01044 SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid );
01045 if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue;
01046 notShrinkFace = true;
01047 for ( unsigned j = 0; j < _sdVec.size(); ++j )
01048 {
01049 if ( _sdVec[j]._solid.IsSame( *solid ) )
01050 if ( _sdVec[j]._shrinkShape2Shape.count( edgeID ))
01051 notShrinkFace = false;
01052 }
01053 }
01054 if ( notShrinkFace )
01055 {
01056 _sdVec[i]._noShrinkFaces.insert( getMeshDS()->ShapeToIndex( fWOL ));
01057 for ( TopExp_Explorer vExp( fWOL, TopAbs_VERTEX ); vExp.More(); vExp.Next() )
01058 noShrinkVertices.Add( vExp.Current() );
01059 }
01060 }
01061
01062
01063 e2f = _sdVec[i]._shrinkShape2Shape.begin();
01064 for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); )
01065 {
01066 TGeomID edgeID = e2f->first;
01067 TopoDS_Vertex VV[2];
01068 TopExp::Vertices( TopoDS::Edge( getMeshDS()->IndexToShape( edgeID )),VV[0],VV[1]);
01069 if ( noShrinkVertices.Contains( VV[0] ) || noShrinkVertices.Contains( VV[1] ))
01070 _sdVec[i]._shrinkShape2Shape.erase( e2f++ );
01071 else
01072 e2f++;
01073 }
01074 }
01075
01076
01077
01078 for ( unsigned i = 0; i < _sdVec.size(); ++i )
01079 {
01080 shapes.Clear();
01081 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_VERTEX, shapes);
01082 for ( int iV = 1; iV <= shapes.Extent(); ++iV )
01083 {
01084 const TopoDS_Shape& vertex = shapes(iV);
01085
01086 vector< TopoDS_Shape > facesWOL;
01087 int totalNbFaces = 0;
01088 PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE);
01089 while ( fIt->more())
01090 {
01091 const TopoDS_Shape* f = fIt->next();
01092 const int fID = getMeshDS()->ShapeToIndex( *f );
01093 if ( helper.IsSubShape( *f, _sdVec[i]._solid ) )
01094 {
01095 totalNbFaces++;
01096 if ( _ignoreShapeIds.count ( fID ) && ! _sdVec[i]._noShrinkFaces.count( fID ))
01097 facesWOL.push_back( *f );
01098 }
01099 }
01100 if ( facesWOL.size() == totalNbFaces || facesWOL.empty() )
01101 continue;
01102 TGeomID vInd = getMeshDS()->ShapeToIndex( vertex );
01103 switch ( facesWOL.size() )
01104 {
01105 case 1:
01106 {
01107 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] )); break;
01108 }
01109 case 2:
01110 {
01111
01112 PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
01113 while ( eIt->more())
01114 {
01115 const TopoDS_Shape* e = eIt->next();
01116 if ( helper.IsSubShape( *e, facesWOL[0]) &&
01117 helper.IsSubShape( *e, facesWOL[1]))
01118 {
01119 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break;
01120 }
01121 }
01122 break;
01123 }
01124 default:
01125 return error("Not yet supported case", _sdVec[i]._index);
01126 }
01127 }
01128 }
01129
01130 return true;
01131 }
01132
01133
01137
01138
01139 bool _ViscousBuilder::makeLayer(_SolidData& data)
01140 {
01141
01142 set<TGeomID> subIds, faceIds;
01143 subIds = data._noShrinkFaces;
01144 TopExp_Explorer exp( data._solid, TopAbs_FACE );
01145 for ( ; exp.More(); exp.Next() )
01146 if ( ! _ignoreShapeIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
01147 {
01148 SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
01149 faceIds.insert( fSubM->GetId() );
01150 SMESH_subMeshIteratorPtr subIt =
01151 fSubM->getDependsOnIterator(true, false);
01152 while ( subIt->more() )
01153 subIds.insert( subIt->next()->GetId() );
01154 }
01155
01156
01157 map< TGeomID, TNode2Edge* > s2neMap;
01158 map< TGeomID, TNode2Edge* >::iterator s2ne;
01159 map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
01160 for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
01161 {
01162 TGeomID shapeInd = s2s->first;
01163 for ( unsigned i = 0; i < _sdVec.size(); ++i )
01164 {
01165 if ( _sdVec[i]._index == data._index ) continue;
01166 map< TGeomID, TopoDS_Shape >::iterator s2s2 = _sdVec[i]._shrinkShape2Shape.find( shapeInd );
01167 if ( s2s2 != _sdVec[i]._shrinkShape2Shape.end() &&
01168 *s2s == *s2s2 && !_sdVec[i]._n2eMap.empty() )
01169 {
01170 s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap ));
01171 break;
01172 }
01173 }
01174 }
01175
01176
01177
01178 dumpFunction(SMESH_Comment("makeLayers_")<<data._index);
01179
01180 data._stepSize = Precision::Infinite();
01181 data._stepSizeNodes[0] = 0;
01182
01183 SMESH_MesherHelper helper( *_mesh );
01184 helper.SetSubShape( data._solid );
01185 helper.SetElementsOnShape(true);
01186
01187 vector< const SMDS_MeshNode*> newNodes;
01188 TNode2Edge::iterator n2e2;
01189
01190
01191 const int nbShapes = getMeshDS()->MaxShapeIndex();
01192 vector< vector<_LayerEdge*> > edgesByGeom( nbShapes+1 );
01193
01194 for ( set<TGeomID>::iterator id = faceIds.begin(); id != faceIds.end(); ++id )
01195 {
01196 SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( *id );
01197 if ( !smDS ) return error(SMESH_Comment("Not meshed face ") << *id, data._index );
01198
01199 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( *id ));
01200 SMESH_ProxyMesh::SubMesh* proxySub =
01201 data._proxyMesh->getFaceSubM( F, true);
01202
01203 SMDS_ElemIteratorPtr eIt = smDS->GetElements();
01204 while ( eIt->more() )
01205 {
01206 const SMDS_MeshElement* face = eIt->next();
01207 newNodes.resize( face->NbCornerNodes() );
01208 double faceMaxCosin = -1;
01209 for ( int i = 0 ; i < face->NbCornerNodes(); ++i )
01210 {
01211 const SMDS_MeshNode* n = face->GetNode(i);
01212 TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first;
01213 if ( !(*n2e).second )
01214 {
01215
01216 _LayerEdge* edge = new _LayerEdge();
01217 n2e->second = edge;
01218 edge->_nodes.push_back( n );
01219 const int shapeID = n->getshapeId();
01220 edgesByGeom[ shapeID ].push_back( edge );
01221
01222
01223 if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
01224 ( s2ne = s2neMap.find( shapeID )) != s2neMap.end() &&
01225 ( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end())
01226 {
01227 _LayerEdge* foundEdge = (*n2e2).second;
01228 edge->Copy( *foundEdge, helper );
01229
01230
01231 const_cast< SMDS_MeshNode* >
01232 ( edge->_nodes.back() )->setXYZ( n->X(), n->Y(), n->Z() );
01233 }
01234 else
01235 {
01236 edge->_nodes.push_back( helper.AddNode( n->X(), n->Y(), n->Z() ));
01237 if ( !setEdgeData( *edge, subIds, helper, data ))
01238 return false;
01239 }
01240 dumpMove(edge->_nodes.back());
01241 if ( edge->_cosin > 0.01 )
01242 {
01243 if ( edge->_cosin > faceMaxCosin )
01244 faceMaxCosin = edge->_cosin;
01245 }
01246 }
01247 newNodes[ i ] = n2e->second->_nodes.back();
01248 }
01249
01250 const SMDS_MeshElement* newFace = new TmpMeshFace( newNodes, --_tmpFaceID );
01251 proxySub->AddElement( newFace );
01252
01253
01254 if ( faceMaxCosin > 0.1 )
01255 limitStepSize( data, face, faceMaxCosin );
01256 }
01257 }
01258
01259 data._epsilon = 1e-7;
01260 if ( data._stepSize < 1. )
01261 data._epsilon *= data._stepSize;
01262
01263
01264
01265 if ( !sortEdges( data, edgesByGeom ))
01266 return false;
01267
01268
01269 TNode2Edge::iterator n2e;
01270 for ( unsigned i = 0; i < data._edges.size(); ++i )
01271 {
01272 if ( data._edges[i]->IsOnEdge())
01273 for ( int j = 0; j < 2; ++j )
01274 {
01275 if ( data._edges[i]->_nodes.back()->NbInverseElements(SMDSAbs_Volume) > 0 )
01276 break;
01277 const SMDS_MeshNode* & n = data._edges[i]->_2neibors->_nodes[j];
01278 if (( n2e = data._n2eMap.find( n )) == data._n2eMap.end() )
01279 return error("_LayerEdge not found by src node", data._index);
01280 n = (*n2e).second->_nodes.back();
01281 data._edges[i]->_2neibors->_edges[j] = n2e->second;
01282 }
01283 else
01284 for ( unsigned j = 0; j < data._edges[i]->_simplices.size(); ++j )
01285 {
01286 _Simplex& s = data._edges[i]->_simplices[j];
01287 s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
01288 s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
01289 }
01290 }
01291
01292 dumpFunctionEnd();
01293 return true;
01294 }
01295
01296
01300
01301
01302 void _ViscousBuilder::limitStepSize( _SolidData& data,
01303 const SMDS_MeshElement* face,
01304 const double cosin)
01305 {
01306 int iN = 0;
01307 double minSize = 10 * data._stepSize;
01308 const int nbNodes = face->NbCornerNodes();
01309 for ( int i = 0; i < nbNodes; ++i )
01310 {
01311 const SMDS_MeshNode* nextN = face->GetNode( SMESH_MesherHelper::WrapIndex( i+1, nbNodes ));
01312 const SMDS_MeshNode* curN = face->GetNode( i );
01313 if ( nextN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ||
01314 curN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
01315 {
01316 double dist = SMESH_TNodeXYZ( face->GetNode(i)).Distance( nextN );
01317 if ( dist < minSize )
01318 minSize = dist, iN = i;
01319 }
01320 }
01321 double newStep = 0.8 * minSize / cosin;
01322 if ( newStep < data._stepSize )
01323 {
01324 data._stepSize = newStep;
01325 data._stepSizeCoeff = 0.8 / cosin;
01326 data._stepSizeNodes[0] = face->GetNode( iN );
01327 data._stepSizeNodes[1] = face->GetNode( SMESH_MesherHelper::WrapIndex( iN+1, nbNodes ));
01328 }
01329 }
01330
01331
01335
01336
01337 void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize)
01338 {
01339 if ( minSize < data._stepSize )
01340 {
01341 data._stepSize = minSize;
01342 if ( data._stepSizeNodes[0] )
01343 {
01344 double dist =
01345 SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
01346 data._stepSizeCoeff = data._stepSize / dist;
01347 }
01348 }
01349 }
01350
01351
01355
01356
01357 bool _ViscousBuilder::sortEdges( _SolidData& data,
01358 vector< vector<_LayerEdge*> >& edgesByGeom)
01359 {
01360
01361
01362
01363 list< TGeomID > shapesToSmooth;
01364
01365 SMESH_MesherHelper helper( *_mesh );
01366 bool ok;
01367
01368 for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS )
01369 {
01370 vector<_LayerEdge*>& eS = edgesByGeom[iS];
01371 if ( eS.empty() ) continue;
01372 TopoDS_Shape S = getMeshDS()->IndexToShape( iS );
01373 bool needSmooth = false;
01374 switch ( S.ShapeType() )
01375 {
01376 case TopAbs_EDGE: {
01377
01378 bool isShrinkEdge = !eS[0]->_sWOL.IsNull();
01379 for ( TopoDS_Iterator vIt( S ); vIt.More() && !needSmooth; vIt.Next() )
01380 {
01381 TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
01382 vector<_LayerEdge*>& eV = edgesByGeom[ iV ];
01383 if ( eV.empty() ) continue;
01384 double cosin = eV[0]->_cosin;
01385 bool badCosin =
01386 ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge));
01387 if ( badCosin )
01388 {
01389 gp_Vec dir1, dir2;
01390 if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE )
01391 dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ));
01392 else
01393 dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ),
01394 eV[0]->_nodes[0], helper, ok);
01395 dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
01396 double angle = dir1.Angle( dir2 );
01397 cosin = cos( angle );
01398 }
01399 needSmooth = ( cosin > 0.1 );
01400 }
01401 break;
01402 }
01403 case TopAbs_FACE: {
01404
01405 for ( TopExp_Explorer eExp( S, TopAbs_EDGE ); eExp.More() && !needSmooth; eExp.Next() )
01406 {
01407 TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
01408 vector<_LayerEdge*>& eE = edgesByGeom[ iE ];
01409 if ( eE.empty() ) continue;
01410 if ( eE[0]->_sWOL.IsNull() )
01411 {
01412 for ( unsigned i = 0; i < eE.size() && !needSmooth; ++i )
01413 needSmooth = ( eE[i]->_cosin > 0.1 );
01414 }
01415 else
01416 {
01417 const TopoDS_Face& F1 = TopoDS::Face( S );
01418 const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
01419 const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
01420 for ( unsigned i = 0; i < eE.size() && !needSmooth; ++i )
01421 {
01422 gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
01423 gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
01424 double angle = dir1.Angle( dir2 );
01425 double cosin = cos( angle );
01426 needSmooth = ( cosin > 0.1 );
01427 }
01428 }
01429 }
01430 break;
01431 }
01432 case TopAbs_VERTEX:
01433 continue;
01434 default:;
01435 }
01436 if ( needSmooth )
01437 {
01438 if ( S.ShapeType() == TopAbs_EDGE ) shapesToSmooth.push_front( iS );
01439 else shapesToSmooth.push_back ( iS );
01440 }
01441
01442 }
01443
01444 data._edges.reserve( data._n2eMap.size() );
01445 data._endEdgeToSmooth.clear();
01446
01447
01448 list< TGeomID >::iterator gIt = shapesToSmooth.begin();
01449 for ( ; gIt != shapesToSmooth.end(); ++gIt )
01450 {
01451 vector<_LayerEdge*>& eVec = edgesByGeom[ *gIt ];
01452 if ( eVec.empty() ) continue;
01453 data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
01454 data._endEdgeToSmooth.push_back( data._edges.size() );
01455 eVec.clear();
01456 }
01457
01458
01459 for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS )
01460 {
01461 vector<_LayerEdge*>& eVec = edgesByGeom[iS];
01462 data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
01463 eVec.clear();
01464 }
01465
01466 return ok;
01467 }
01468
01469
01474
01475
01476 bool _ViscousBuilder::setEdgeData(_LayerEdge& edge,
01477 const set<TGeomID>& subIds,
01478 SMESH_MesherHelper& helper,
01479 _SolidData& data)
01480 {
01481 SMESH_MeshEditor editor(_mesh);
01482
01483 const SMDS_MeshNode* node = edge._nodes[0];
01484 SMDS_TypeOfPosition posType = node->GetPosition()->GetTypeOfPosition();
01485
01486 edge._len = 0;
01487 edge._2neibors = 0;
01488 edge._curvature = 0;
01489
01490
01491
01492
01493
01494 edge._cosin = 0;
01495 edge._normal.SetCoord(0,0,0);
01496
01497 int totalNbFaces = 0;
01498 gp_Pnt p;
01499 gp_Vec du, dv, geomNorm;
01500 bool normOK = true;
01501
01502 TGeomID shapeInd = node->getshapeId();
01503 map< TGeomID, TopoDS_Shape >::const_iterator s2s = data._shrinkShape2Shape.find( shapeInd );
01504 bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
01505 TopoDS_Shape vertEdge;
01506
01507 if ( onShrinkShape )
01508 {
01509 vertEdge = getMeshDS()->IndexToShape( s2s->first );
01510 if ( s2s->second.ShapeType() == TopAbs_EDGE )
01511 {
01512
01513 edge._normal = getEdgeDir( TopoDS::Edge( s2s->second ), TopoDS::Vertex( vertEdge ));
01514 }
01515 else if ( vertEdge.ShapeType() == TopAbs_VERTEX )
01516 {
01517
01518 edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Vertex( vertEdge ),
01519 node, helper, normOK, &edge._cosin);
01520 }
01521 else
01522 {
01523
01524 edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Edge( vertEdge ),
01525 node, helper, normOK);
01526 }
01527 }
01528 else
01529 {
01530
01531 set<TGeomID> faceIds;
01532 if ( posType == SMDS_TOP_FACE )
01533 {
01534 faceIds.insert( node->getshapeId() );
01535 }
01536 else
01537 {
01538 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
01539 while ( fIt->more() )
01540 faceIds.insert( editor.FindShape(fIt->next()));
01541 }
01542
01543 set<TGeomID>::iterator id = faceIds.begin();
01544 TopoDS_Face F;
01545 for ( ; id != faceIds.end(); ++id )
01546 {
01547 const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id );
01548 if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
01549 continue;
01550 totalNbFaces++;
01551
01552 F = TopoDS::Face( s );
01553
01554 gp_XY uv = helper.GetNodeUV( F, node, 0, &normOK );
01555 Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
01556 surface->D1( uv.X(),uv.Y(), p, du,dv );
01557 geomNorm = du ^ dv;
01558 double size2 = geomNorm.SquareMagnitude();
01559 if ( size2 > numeric_limits<double>::min() )
01560 geomNorm /= sqrt( size2 );
01561 else
01562 normOK = false;
01563 if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
01564 geomNorm.Reverse();
01565 edge._normal += geomNorm.XYZ();
01566 }
01567 if ( totalNbFaces == 0 )
01568 return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index);
01569
01570 edge._normal /= totalNbFaces;
01571
01572 switch ( posType )
01573 {
01574 case SMDS_TOP_FACE:
01575 edge._cosin = 0; break;
01576
01577 case SMDS_TOP_EDGE: {
01578 TopoDS_Edge E = TopoDS::Edge( helper.GetSubShapeByNode( node, getMeshDS()));
01579 gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK);
01580 double angle = inFaceDir.Angle( edge._normal );
01581 edge._cosin = cos( angle );
01582
01583 break;
01584 }
01585 case SMDS_TOP_VERTEX: {
01586 TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS()));
01587 gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK);
01588 double angle = inFaceDir.Angle( edge._normal );
01589 edge._cosin = cos( angle );
01590
01591 break;
01592 }
01593 default:
01594 return error(SMESH_Comment("Invalid shape position of node ")<<node, data._index);
01595 }
01596 }
01597
01598 double normSize = edge._normal.SquareModulus();
01599 if ( normSize < numeric_limits<double>::min() )
01600 return error(SMESH_Comment("Bad normal at node ")<< node->GetID(), data._index );
01601
01602 edge._normal /= sqrt( normSize );
01603
01604
01605
01606
01607
01608 if ( onShrinkShape )
01609 {
01610 edge._sWOL = (*s2s).second;
01611
01612 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( edge._nodes.back() );
01613 if ( SMESHDS_SubMesh* sm = getMeshDS()->MeshElements( data._solid ))
01614 sm->RemoveNode( tgtNode , false );
01615
01616
01617 if ( edge._sWOL.ShapeType() == TopAbs_EDGE )
01618 {
01619 double u = helper.GetNodeU( TopoDS::Edge( edge._sWOL ), node, 0, &normOK );
01620 edge._pos.push_back( gp_XYZ( u, 0, 0));
01621 getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( edge._sWOL ), u );
01622 }
01623 else
01624 {
01625 gp_XY uv = helper.GetNodeUV( TopoDS::Face( edge._sWOL ), node, 0, &normOK );
01626 edge._pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
01627 getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( edge._sWOL ), uv.X(), uv.Y() );
01628 }
01629 }
01630 else
01631 {
01632 edge._pos.push_back( SMESH_TNodeXYZ( node ));
01633
01634 if ( posType == SMDS_TOP_FACE )
01635 {
01636 getSimplices( node, edge._simplices, _ignoreShapeIds, &data );
01637 double avgNormProj = 0, avgLen = 0;
01638 for ( unsigned i = 0; i < edge._simplices.size(); ++i )
01639 {
01640 gp_XYZ vec = edge._pos.back() - SMESH_TNodeXYZ( edge._simplices[i]._nPrev );
01641 avgNormProj += edge._normal * vec;
01642 avgLen += vec.Modulus();
01643 }
01644 avgNormProj /= edge._simplices.size();
01645 avgLen /= edge._simplices.size();
01646 edge._curvature = _Curvature::New( avgNormProj, avgLen );
01647 }
01648 }
01649
01650
01651
01652 if ( posType == SMDS_TOP_EDGE
01653 )
01654 {
01655 edge._2neibors = new _2NearEdges;
01656
01657 if ( ! findNeiborsOnEdge( &edge,
01658 edge._2neibors->_nodes[0],
01659 edge._2neibors->_nodes[1],
01660 data))
01661 return false;
01662 edge.SetDataByNeighbors( edge._2neibors->_nodes[0],
01663 edge._2neibors->_nodes[1],
01664 helper);
01665 }
01666
01667 edge.SetCosin( edge._cosin );
01668
01669 return true;
01670 }
01671
01672
01676
01677
01678 bool _ViscousBuilder::findNeiborsOnEdge(const _LayerEdge* edge,
01679 const SMDS_MeshNode*& n1,
01680 const SMDS_MeshNode*& n2,
01681 _SolidData& data)
01682 {
01683 const SMDS_MeshNode* node = edge->_nodes[0];
01684 const int shapeInd = node->getshapeId();
01685 SMESHDS_SubMesh* edgeSM = 0;
01686 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE )
01687 {
01688
01689 edgeSM = getMeshDS()->MeshElements( shapeInd );
01690 if ( !edgeSM || edgeSM->NbElements() == 0 )
01691 return error(SMESH_Comment("Not meshed EDGE ") << shapeInd, data._index);
01692 }
01693 int iN = 0;
01694 n2 = 0;
01695 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Edge);
01696 while ( eIt->more() && !n2 )
01697 {
01698 const SMDS_MeshElement* e = eIt->next();
01699 const SMDS_MeshNode* nNeibor = e->GetNode( 0 );
01700 if ( nNeibor == node ) nNeibor = e->GetNode( 1 );
01701 if ( edgeSM )
01702 {
01703 if (!edgeSM->Contains(e)) continue;
01704 }
01705 else
01706 {
01707 TopoDS_Shape s = SMESH_MesherHelper::GetSubShapeByNode(nNeibor, getMeshDS() );
01708 if ( !SMESH_MesherHelper::IsSubShape( s, edge->_sWOL )) continue;
01709 }
01710 ( iN++ ? n2 : n1 ) = nNeibor;
01711 }
01712 if ( !n2 )
01713 return error(SMESH_Comment("Wrongly meshed EDGE ") << shapeInd, data._index);
01714 return true;
01715 }
01716
01717
01721
01722
01723 void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
01724 const SMDS_MeshNode* n2,
01725 SMESH_MesherHelper& helper)
01726 {
01727 if ( _nodes[0]->GetPosition()->GetTypeOfPosition() != SMDS_TOP_EDGE )
01728 return;
01729
01730 gp_XYZ pos = SMESH_TNodeXYZ( _nodes[0] );
01731 gp_XYZ vec1 = pos - SMESH_TNodeXYZ( n1 );
01732 gp_XYZ vec2 = pos - SMESH_TNodeXYZ( n2 );
01733
01734
01735
01736 double sumLen = vec1.Modulus() + vec2.Modulus();
01737 _2neibors->_wgt[0] = 1 - vec1.Modulus() / sumLen;
01738 _2neibors->_wgt[1] = 1 - vec2.Modulus() / sumLen;
01739 double avgNormProj = 0.5 * ( _normal * vec1 + _normal * vec2 );
01740 double avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() );
01741 if ( _curvature ) delete _curvature;
01742 _curvature = _Curvature::New( avgNormProj, avgLen );
01743 #ifdef __myDEBUG
01744
01745
01746
01747
01748
01749 #endif
01750
01751
01752
01753 if ( _sWOL.IsNull() )
01754 {
01755 TopoDS_Shape S = helper.GetSubShapeByNode( _nodes[0], helper.GetMeshDS() );
01756 gp_XYZ dirE = getEdgeDir( TopoDS::Edge( S ), _nodes[0], helper );
01757 gp_XYZ plnNorm = dirE ^ _normal;
01758 double proj0 = plnNorm * vec1;
01759 double proj1 = plnNorm * vec2;
01760 if ( fabs( proj0 ) > 1e-10 || fabs( proj1 ) > 1e-10 )
01761 {
01762 if ( _2neibors->_plnNorm ) delete _2neibors->_plnNorm;
01763 _2neibors->_plnNorm = new gp_XYZ( plnNorm.Normalized() );
01764 }
01765 }
01766 }
01767
01768
01773
01774
01775 void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
01776 {
01777 _nodes = other._nodes;
01778 _normal = other._normal;
01779 _len = 0;
01780 _lenFactor = other._lenFactor;
01781 _cosin = other._cosin;
01782 _sWOL = other._sWOL;
01783 _2neibors = other._2neibors;
01784 _curvature = 0; std::swap( _curvature, other._curvature );
01785 _2neibors = 0; std::swap( _2neibors, other._2neibors );
01786
01787 if ( _sWOL.ShapeType() == TopAbs_EDGE )
01788 {
01789 double u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes[0] );
01790 _pos.push_back( gp_XYZ( u, 0, 0));
01791 }
01792 else
01793 {
01794 gp_XY uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes[0]);
01795 _pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
01796 }
01797 }
01798
01799
01803
01804
01805 void _LayerEdge::SetCosin( double cosin )
01806 {
01807 _cosin = cosin;
01808 _lenFactor = ( _cosin > 0.1 ) ? 1./sqrt(1-_cosin*_cosin) : 1.0;
01809 }
01810
01811
01815
01816
01817 void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node,
01818 vector<_Simplex>& simplices,
01819 const set<TGeomID>& ingnoreShapes,
01820 const _SolidData* dataToCheckOri)
01821 {
01822 SMESH_MeshEditor editor( _mesh );
01823 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
01824 while ( fIt->more() )
01825 {
01826 const SMDS_MeshElement* f = fIt->next();
01827 const TGeomID shapeInd = editor.FindShape( f );
01828 if ( ingnoreShapes.count( shapeInd )) continue;
01829 const int nbNodes = f->NbCornerNodes();
01830 int srcInd = f->GetNodeIndex( node );
01831 const SMDS_MeshNode* nPrev = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd-1, nbNodes ));
01832 const SMDS_MeshNode* nNext = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+1, nbNodes ));
01833 if ( dataToCheckOri && dataToCheckOri->_reversedFaceIds.count( shapeInd ))
01834 std::swap( nPrev, nNext );
01835 simplices.push_back( _Simplex( nPrev, nNext ));
01836 }
01837 simplices.resize( simplices.size() );
01838 }
01839
01840
01844
01845
01846 void _ViscousBuilder::makeGroupOfLE()
01847 {
01848 #ifdef _DEBUG_
01849 for ( unsigned i = 0 ; i < _sdVec.size(); ++i )
01850 {
01851 if ( _sdVec[i]._edges.empty() ) continue;
01852
01853
01854
01855
01856
01857
01858 dumpFunction( SMESH_Comment("make_LayerEdge_") << i );
01859 for ( unsigned j = 0 ; j < _sdVec[i]._edges.size(); ++j )
01860 {
01861 _LayerEdge* le = _sdVec[i]._edges[j];
01862 for ( unsigned iN = 1; iN < le->_nodes.size(); ++iN )
01863 dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<le->_nodes[iN-1]->GetID()
01864 << ", " << le->_nodes[iN]->GetID() <<"])");
01865
01866 }
01867 dumpFunctionEnd();
01868
01869 dumpFunction( SMESH_Comment("makeNormals") << i );
01870 for ( unsigned j = 0 ; j < _sdVec[i]._edges.size(); ++j )
01871 {
01872 _LayerEdge& edge = *_sdVec[i]._edges[j];
01873 SMESH_TNodeXYZ nXYZ( edge._nodes[0] );
01874 nXYZ += edge._normal * _sdVec[i]._stepSize;
01875 dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<edge._nodes[0]->GetID()
01876 << ", mesh.AddNode( " << nXYZ.X()<<","<< nXYZ.Y()<<","<< nXYZ.Z()<<")])");
01877 }
01878 dumpFunctionEnd();
01879
01880
01881
01882
01883
01884 dumpFunction( SMESH_Comment("makeTmpFaces_") << i );
01885 TopExp_Explorer fExp( _sdVec[i]._solid, TopAbs_FACE );
01886 for ( ; fExp.More(); fExp.Next() )
01887 {
01888 if (const SMESHDS_SubMesh* sm = _sdVec[i]._proxyMesh->GetProxySubMesh( fExp.Current()))
01889 {
01890 SMDS_ElemIteratorPtr fIt = sm->GetElements();
01891 while ( fIt->more())
01892 {
01893 const SMDS_MeshElement* e = fIt->next();
01894 SMESH_Comment cmd("mesh.AddFace([");
01895 for ( int j=0; j < e->NbCornerNodes(); ++j )
01896 cmd << e->GetNode(j)->GetID() << (j+1<e->NbCornerNodes() ? ",": "])");
01897 dumpCmd( cmd );
01898
01899
01900 }
01901 }
01902 }
01903 dumpFunctionEnd();
01904 }
01905 #endif
01906 }
01907
01908
01912
01913
01914 bool _ViscousBuilder::inflate(_SolidData& data)
01915 {
01916 SMESH_MesherHelper helper( *_mesh );
01917
01918
01919
01920 double geomSize = Precision::Infinite(), intersecDist;
01921 SMESH_MeshEditor editor( _mesh );
01922 auto_ptr<SMESH_ElementSearcher> searcher
01923 ( editor.GetElementSearcher( data._proxyMesh->GetFaces( data._solid )) );
01924 for ( unsigned i = 0; i < data._edges.size(); ++i )
01925 {
01926 if ( data._edges[i]->IsOnEdge() ) continue;
01927 data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon );
01928 if ( geomSize > intersecDist )
01929 geomSize = intersecDist;
01930 }
01931 if ( data._stepSize > 0.3 * geomSize )
01932 limitStepSize( data, 0.3 * geomSize );
01933
01934 const double tgtThick = data._hyp->GetTotalThickness();
01935 if ( data._stepSize > tgtThick )
01936 limitStepSize( data, tgtThick );
01937
01938 if ( data._stepSize < 1. )
01939 data._epsilon = data._stepSize * 1e-7;
01940
01941 #ifdef __myDEBUG
01942 cout << "-- geomSize = " << geomSize << ", stepSize = " << data._stepSize << endl;
01943 #endif
01944
01945 double avgThick = 0, curThick = 0, distToIntersection = Precision::Infinite();
01946 int nbSteps = 0, nbRepeats = 0;
01947 while ( 1.01 * avgThick < tgtThick )
01948 {
01949
01950 curThick += data._stepSize;
01951 if ( curThick > tgtThick )
01952 {
01953 curThick = tgtThick + ( tgtThick-avgThick ) * nbRepeats;
01954 nbRepeats++;
01955 }
01956
01957
01958 dumpFunction(SMESH_Comment("inflate")<<data._index<<"_step"<<nbSteps);
01959 for ( unsigned i = 0; i < data._edges.size(); ++i )
01960 {
01961 data._edges[i]->SetNewLength( curThick, helper );
01962 }
01963 dumpFunctionEnd();
01964
01965 if ( !nbSteps )
01966 if ( !updateNormals( data, helper ) )
01967 return false;
01968
01969
01970 if ( !smoothAndCheck( data, nbSteps, distToIntersection ))
01971 {
01972 if ( nbSteps > 0 )
01973 {
01974 dumpFunction(SMESH_Comment("invalidate")<<data._index<<"_step"<<nbSteps);
01975 for ( unsigned i = 0; i < data._edges.size(); ++i )
01976 {
01977 data._edges[i]->InvalidateStep( nbSteps+1 );
01978 }
01979 dumpFunctionEnd();
01980 }
01981 break;
01982 }
01983 nbSteps++;
01984
01985
01986 avgThick = 0;
01987 for ( unsigned i = 0; i < data._edges.size(); ++i )
01988 avgThick += data._edges[i]->_len;
01989 avgThick /= data._edges.size();
01990 #ifdef __myDEBUG
01991 cout << "-- Thickness " << avgThick << " reached" << endl;
01992 #endif
01993
01994 if ( distToIntersection < avgThick*1.5 )
01995 {
01996 #ifdef __myDEBUG
01997 cout << "-- Stop inflation since distToIntersection( "<<distToIntersection<<" ) < avgThick( "
01998 << avgThick << " ) * 1.5" << endl;
01999 #endif
02000 break;
02001 }
02002
02003 limitStepSize( data, 0.25 * distToIntersection );
02004 if ( data._stepSizeNodes[0] )
02005 data._stepSize = data._stepSizeCoeff *
02006 SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
02007 }
02008
02009 if (nbSteps == 0 )
02010 return error("failed at the very first inflation step", data._index);
02011
02012 return true;
02013 }
02014
02015
02019
02020
02021 bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
02022 const int nbSteps,
02023 double & distToIntersection)
02024 {
02025 if ( data._endEdgeToSmooth.empty() )
02026 return true;
02027
02028 bool moved, improved;
02029
02030 SMESH_MesherHelper helper(*_mesh);
02031 Handle(Geom_Surface) surface;
02032 TopoDS_Face F;
02033
02034 int iBeg, iEnd = 0;
02035 for ( unsigned iS = 0; iS < data._endEdgeToSmooth.size(); ++iS )
02036 {
02037 iBeg = iEnd;
02038 iEnd = data._endEdgeToSmooth[ iS ];
02039
02040 if ( !data._edges[ iBeg ]->_sWOL.IsNull() &&
02041 data._edges[ iBeg ]->_sWOL.ShapeType() == TopAbs_FACE )
02042 {
02043 if ( !F.IsSame( data._edges[ iBeg ]->_sWOL )) {
02044 F = TopoDS::Face( data._edges[ iBeg ]->_sWOL );
02045 helper.SetSubShape( F );
02046 surface = BRep_Tool::Surface( F );
02047 }
02048 }
02049 else
02050 {
02051 F.Nullify(); surface.Nullify();
02052 }
02053 TGeomID sInd = data._edges[ iBeg ]->_nodes[0]->getshapeId();
02054
02055 if ( data._edges[ iBeg ]->IsOnEdge() )
02056 {
02057 if ( !smoothAnalyticEdge( data, iBeg, iEnd, surface, F, helper ))
02058 {
02059 dumpFunction(SMESH_Comment("smooth")<<data._index << "_Ed"<<sInd <<"_InfStep"<<nbSteps);
02060
02061 int step = 0;
02062 do {
02063 moved = false;
02064 for ( int i = iBeg; i < iEnd; ++i )
02065 {
02066 moved |= data._edges[i]->SmoothOnEdge(surface, F, helper);
02067 }
02068 dumpCmd( SMESH_Comment("# end step ")<<step);
02069 }
02070 while ( moved && step++ < 5 );
02071
02072
02073 dumpFunctionEnd();
02074 }
02075 }
02076 else
02077 {
02078
02079 int step = 0, badNb = 0; moved = true;
02080 while (( ++step <= 5 && moved ) || improved )
02081 {
02082 dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
02083 <<"_InfStep"<<nbSteps<<"_"<<step);
02084 int oldBadNb = badNb;
02085 badNb = 0;
02086 moved = false;
02087 for ( int i = iBeg; i < iEnd; ++i )
02088 moved |= data._edges[i]->Smooth(badNb);
02089 improved = ( badNb < oldBadNb );
02090
02091 dumpFunctionEnd();
02092 }
02093 if ( badNb > 0 )
02094 {
02095 #ifdef __myDEBUG
02096 for ( int i = iBeg; i < iEnd; ++i )
02097 {
02098 _LayerEdge* edge = data._edges[i];
02099 SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
02100 for ( unsigned j = 0; j < edge->_simplices.size(); ++j )
02101 if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
02102 {
02103 cout << "Bad simplex ( " << edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
02104 << " "<< edge->_simplices[j]._nPrev->GetID()
02105 << " "<< edge->_simplices[j]._nNext->GetID() << " )" << endl;
02106 return false;
02107 }
02108 }
02109 #endif
02110 return false;
02111 }
02112 }
02113 }
02114
02115
02116
02117
02118 SMESH_MeshEditor editor( _mesh );
02119 auto_ptr<SMESH_ElementSearcher> searcher
02120 ( editor.GetElementSearcher( data._proxyMesh->GetFaces( data._solid )) );
02121
02122 distToIntersection = Precision::Infinite();
02123 double dist;
02124 const SMDS_MeshElement* intFace = 0;
02125 #ifdef __myDEBUG
02126 const SMDS_MeshElement* closestFace = 0;
02127 int iLE = 0;
02128 #endif
02129 for ( unsigned i = 0; i < data._edges.size(); ++i )
02130 {
02131 if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace ))
02132 return false;
02133 if ( distToIntersection > dist )
02134 {
02135 distToIntersection = dist;
02136 #ifdef __myDEBUG
02137 iLE = i;
02138 closestFace = intFace;
02139 #endif
02140 }
02141 }
02142 #ifdef __myDEBUG
02143 if ( closestFace )
02144 {
02145 SMDS_MeshElement::iterator nIt = closestFace->begin_nodes();
02146 cout << "Shortest distance: _LayerEdge nodes: tgt " << data._edges[iLE]->_nodes.back()->GetID()
02147 << " src " << data._edges[iLE]->_nodes[0]->GetID()<< ", intersection with face ("
02148 << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
02149 << ") distance = " << distToIntersection<< endl;
02150 }
02151 #endif
02152
02153 return true;
02154 }
02155
02156
02161
02162
02163 Handle(Geom_Curve) _SolidData::CurveForSmooth( const TopoDS_Edge& E,
02164 const int iFrom,
02165 const int iTo,
02166 Handle(Geom_Surface)& surface,
02167 const TopoDS_Face& F,
02168 SMESH_MesherHelper& helper)
02169 {
02170 TGeomID eIndex = helper.GetMeshDS()->ShapeToIndex( E );
02171
02172 map< TGeomID, Handle(Geom_Curve)>::iterator i2curve = _edge2curve.find( eIndex );
02173
02174 if ( i2curve == _edge2curve.end() )
02175 {
02176
02177 {
02178 map< double, _LayerEdge* > u2edge;
02179 for ( int i = iFrom; i < iTo; ++i )
02180 u2edge.insert( make_pair( helper.GetNodeU( E, _edges[i]->_nodes[0] ), _edges[i] ));
02181
02182 ASSERT( u2edge.size() == iTo - iFrom );
02183 map< double, _LayerEdge* >::iterator u2e = u2edge.begin();
02184 for ( int i = iFrom; i < iTo; ++i, ++u2e )
02185 _edges[i] = u2e->second;
02186
02187
02188 for ( int i = iFrom; i < iTo-1; ++i )
02189 if ( _edges[i]->_2neibors->_nodes[1] != _edges[i+1]->_nodes.back() )
02190 _edges[i]->_2neibors->reverse();
02191 if ( u2edge.size() > 1 &&
02192 _edges[iTo-1]->_2neibors->_nodes[0] != _edges[iTo-2]->_nodes.back() )
02193 _edges[iTo-1]->_2neibors->reverse();
02194 }
02195
02196 SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( eIndex );
02197
02198 TopLoc_Location loc; double f,l;
02199
02200 Handle(Geom_Line) line;
02201 Handle(Geom_Circle) circle;
02202 bool isLine, isCirc;
02203 if ( F.IsNull() )
02204 {
02205
02206 Handle(Geom_Curve) curve = BRep_Tool::Curve( E, loc, f, l);
02207 if ( curve->IsKind( STANDARD_TYPE( Geom_TrimmedCurve )))
02208 curve = Handle(Geom_TrimmedCurve)::DownCast( curve )->BasisCurve();
02209
02210 line = Handle(Geom_Line)::DownCast( curve );
02211 circle = Handle(Geom_Circle)::DownCast( curve );
02212 isLine = (!line.IsNull());
02213 isCirc = (!circle.IsNull());
02214
02215 if ( !isLine && !isCirc )
02216 {
02217 Bnd_B3d bndBox;
02218 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
02219 while ( nIt->more() )
02220 bndBox.Add( SMESH_TNodeXYZ( nIt->next() ));
02221 gp_XYZ size = bndBox.CornerMax() - bndBox.CornerMin();
02222
02223 SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->_nodes[0] );
02224 SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->_nodes[1] );
02225 const double lineTol = 1e-2 * ( p0 - p1 ).Modulus();
02226 for ( int i = 0; i < 3 && !isLine; ++i )
02227 isLine = ( size.Coord( i+1 ) <= lineTol );
02228 }
02229 if ( !isLine && !isCirc && iTo-iFrom > 2)
02230 {
02231
02232 }
02233 }
02234 else
02235 {
02236
02237 Handle(Geom2d_Curve) curve = BRep_Tool::CurveOnSurface( E, F, f, l);
02238 if ( curve->IsKind( STANDARD_TYPE( Geom2d_TrimmedCurve )))
02239 curve = Handle(Geom2d_TrimmedCurve)::DownCast( curve )->BasisCurve();
02240
02241 Handle(Geom2d_Line) line2d = Handle(Geom2d_Line)::DownCast( curve );
02242 Handle(Geom2d_Circle) circle2d = Handle(Geom2d_Circle)::DownCast( curve );
02243 isLine = (!line2d.IsNull());
02244 isCirc = (!circle2d.IsNull());
02245
02246 if ( !isLine && !isCirc)
02247 {
02248 Bnd_B2d bndBox;
02249 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
02250 while ( nIt->more() )
02251 bndBox.Add( helper.GetNodeUV( F, nIt->next() ));
02252 gp_XY size = bndBox.CornerMax() - bndBox.CornerMin();
02253
02254 const double lineTol = 1e-2 * sqrt( bndBox.SquareExtent() );
02255 for ( int i = 0; i < 2 && !isLine; ++i )
02256 isLine = ( size.Coord( i+1 ) <= lineTol );
02257 }
02258 if ( !isLine && !isCirc && iTo-iFrom > 2)
02259 {
02260
02261 }
02262 if ( isLine )
02263 {
02264 line = new Geom_Line( gp::OX() );
02265 }
02266 else if ( isCirc )
02267 {
02268 gp_Pnt2d p = circle2d->Location();
02269 gp_Ax2 ax( gp_Pnt( p.X(), p.Y(), 0), gp::DX());
02270 circle = new Geom_Circle( ax, 1.);
02271 }
02272 }
02273
02274 Handle(Geom_Curve)& res = _edge2curve[ eIndex ];
02275 if ( isLine )
02276 res = line;
02277 else if ( isCirc )
02278 res = circle;
02279
02280 return res;
02281 }
02282 return i2curve->second;
02283 }
02284
02285
02289
02290
02291 bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data,
02292 const int iFrom,
02293 const int iTo,
02294 Handle(Geom_Surface)& surface,
02295 const TopoDS_Face& F,
02296 SMESH_MesherHelper& helper)
02297 {
02298 TopoDS_Shape S = helper.GetSubShapeByNode( data._edges[ iFrom ]->_nodes[0],
02299 helper.GetMeshDS());
02300 TopoDS_Edge E = TopoDS::Edge( S );
02301
02302 Handle(Geom_Curve) curve = data.CurveForSmooth( E, iFrom, iTo, surface, F, helper );
02303 if ( curve.IsNull() ) return false;
02304
02305
02306 vector< double > len( iTo-iFrom+1 );
02307 {
02308 double curLen, prevLen = len[0] = 1.0;
02309 for ( int i = iFrom; i < iTo; ++i )
02310 {
02311 curLen = prevLen * data._edges[i]->_2neibors->_wgt[0] / data._edges[i]->_2neibors->_wgt[1];
02312 len[i-iFrom+1] = len[i-iFrom] + curLen;
02313 prevLen = curLen;
02314 }
02315 }
02316
02317 if ( curve->IsKind( STANDARD_TYPE( Geom_Line )))
02318 {
02319 if ( F.IsNull() )
02320 {
02321 SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->_nodes[0]);
02322 SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->_nodes[1]);
02323 for ( int i = iFrom; i < iTo; ++i )
02324 {
02325 double r = len[i-iFrom] / len.back();
02326 gp_XYZ newPos = p0 * ( 1. - r ) + p1 * r;
02327 data._edges[i]->_pos.back() = newPos;
02328 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
02329 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
02330 }
02331 }
02332 else
02333 {
02334 gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]);
02335 gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]);
02336 for ( int i = iFrom; i < iTo; ++i )
02337 {
02338 double r = len[i-iFrom] / len.back();
02339 gp_XY newUV = uv0 * ( 1. - r ) + uv1 * r;
02340 data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
02341
02342 gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
02343 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
02344 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
02345
02346 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
02347 pos->SetUParameter( newUV.X() );
02348 pos->SetVParameter( newUV.Y() );
02349 }
02350 }
02351 return true;
02352 }
02353
02354 if ( curve->IsKind( STANDARD_TYPE( Geom_Circle )))
02355 {
02356 Handle(Geom_Circle) circle = Handle(Geom_Circle)::DownCast( curve );
02357 gp_Pnt center3D = circle->Location();
02358
02359 if ( F.IsNull() )
02360 {
02361 return false;
02362 }
02363 else
02364 {
02365 const gp_XY center( center3D.X(), center3D.Y() );
02366
02367 gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]);
02368 gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back());
02369 gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]);
02370 gp_Vec2d vec0( center, uv0 );
02371 gp_Vec2d vecM( center, uvM);
02372 gp_Vec2d vec1( center, uv1 );
02373 double uLast = vec0.Angle( vec1 );
02374 double uMidl = vec0.Angle( vecM );
02375 if ( uLast < 0 ) uLast += 2*PI;
02376 if ( uMidl < 0 ) uMidl += 2*PI;
02377 const bool sense = ( uMidl < uLast );
02378 const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() );
02379
02380 gp_Ax2d axis( center, vec0 );
02381 gp_Circ2d circ ( axis, radius, sense );
02382 for ( int i = iFrom; i < iTo; ++i )
02383 {
02384 double newU = uLast * len[i-iFrom] / len.back();
02385 gp_Pnt2d newUV = ElCLib::Value( newU, circ );
02386 data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
02387
02388 gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
02389 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
02390 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
02391
02392 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
02393 pos->SetUParameter( newUV.X() );
02394 pos->SetVParameter( newUV.Y() );
02395 }
02396 }
02397 return true;
02398 }
02399
02400 return false;
02401 }
02402
02403
02408
02409
02410 bool _ViscousBuilder::updateNormals( _SolidData& data,
02411 SMESH_MesherHelper& helper )
02412 {
02413
02414
02415
02416 vector< const SMDS_MeshElement* > tmpFaces;
02417 {
02418 set< SMESH_TLink > extrudedLinks;
02419 vector< const SMDS_MeshNode*> nodes(4);
02420
02421 dumpFunction(SMESH_Comment("makeTmpFacesOnEdges")<<data._index);
02422 for ( unsigned i = 0; i < data._edges.size(); ++i )
02423 {
02424 _LayerEdge* edge = data._edges[i];
02425 if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
02426 const SMDS_MeshNode* tgt1 = edge->_nodes.back();
02427 for ( int j = 0; j < 2; ++j )
02428 {
02429 const SMDS_MeshNode* tgt2 = edge->_2neibors->_nodes[j];
02430 pair< set< SMESH_TLink >::iterator, bool > link_isnew =
02431 extrudedLinks.insert( SMESH_TLink( tgt1, tgt2 ));
02432 if ( !link_isnew.second )
02433 {
02434 extrudedLinks.erase( link_isnew.first );
02435 continue;
02436 }
02437
02438
02439
02440
02441
02442
02443
02444
02445
02446
02447
02448
02449 _LayerEdge* neiborEdge = edge->_2neibors->_edges[j];
02450
02451 TmpMeshFaceOnEdge* f = new TmpMeshFaceOnEdge( edge, neiborEdge, --_tmpFaceID );
02452 tmpFaces.push_back( f );
02453
02454 dumpCmd(SMESH_Comment("mesh.AddFace([ ")
02455 <<f->_nn[0]->GetID()<<", "<<f->_nn[1]->GetID()<<", "
02456 <<f->_nn[2]->GetID()<<", "<<f->_nn[3]->GetID()<<" ])");
02457 }
02458 }
02459 dumpFunctionEnd();
02460 }
02461
02462
02463
02464
02465
02466 SMESH_MeshEditor editor( _mesh );
02467 SMDS_ElemIteratorPtr fIt( new SMDS_ElementVectorIterator( tmpFaces.begin(),
02468 tmpFaces.end()));
02469 auto_ptr<SMESH_ElementSearcher> searcher ( editor.GetElementSearcher( fIt ));
02470
02471
02472 double dist;
02473 const SMDS_MeshElement* face;
02474 typedef map< _LayerEdge*, set< _LayerEdge*, _LayerEdgeCmp >, _LayerEdgeCmp > TLEdge2LEdgeSet;
02475 TLEdge2LEdgeSet edge2CloseEdge;
02476
02477 const double eps = data._epsilon * data._epsilon;
02478 for ( unsigned i = 0; i < data._edges.size(); ++i )
02479 {
02480 _LayerEdge* edge = data._edges[i];
02481 if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
02482 if ( edge->FindIntersection( *searcher, dist, eps, &face ))
02483 {
02484 const TmpMeshFaceOnEdge* f = (const TmpMeshFaceOnEdge*) face;
02485 set< _LayerEdge*, _LayerEdgeCmp > & ee = edge2CloseEdge[ edge ];
02486 ee.insert( f->_le1 );
02487 ee.insert( f->_le2 );
02488 if ( f->_le1->IsOnEdge() && f->_le1->_sWOL.IsNull() )
02489 edge2CloseEdge[ f->_le1 ].insert( edge );
02490 if ( f->_le2->IsOnEdge() && f->_le2->_sWOL.IsNull() )
02491 edge2CloseEdge[ f->_le2 ].insert( edge );
02492 }
02493 }
02494
02495
02496
02497 if ( !edge2CloseEdge.empty() )
02498 {
02499 dumpFunction(SMESH_Comment("updateNormals")<<data._index);
02500
02501 TLEdge2LEdgeSet::iterator e2ee = edge2CloseEdge.begin();
02502 for ( ; e2ee != edge2CloseEdge.end(); ++e2ee )
02503 {
02504 _LayerEdge* edge1 = e2ee->first;
02505 _LayerEdge* edge2 = 0;
02506 set< _LayerEdge*, _LayerEdgeCmp >& ee = e2ee->second;
02507
02508
02509 TopoDS_Edge E1, E2;
02510 TopoDS_Shape S = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
02511 if ( S.ShapeType() != TopAbs_EDGE )
02512 continue;
02513 E1 = TopoDS::Edge( S );
02514 set< _LayerEdge* >::iterator eIt = ee.begin();
02515 while ( E2.IsNull() && eIt != ee.end())
02516 {
02517 _LayerEdge* e2 = *eIt++;
02518 TopoDS_Shape S = helper.GetSubShapeByNode( e2->_nodes[0], getMeshDS() );
02519 if ( S.ShapeType() == TopAbs_EDGE )
02520 E2 = TopoDS::Edge( S ), edge2 = e2;
02521 }
02522 if ( E2.IsNull() ) continue;
02523
02524
02525
02526 TopoDS_Face FF1[2], FF2[2];
02527 PShapeIteratorPtr fIt = helper.GetAncestors(E1, *_mesh, TopAbs_FACE);
02528 while ( fIt->more() && FF1[1].IsNull())
02529 {
02530 const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
02531 if ( helper.IsSubShape( *F, data._solid))
02532 FF1[ FF1[0].IsNull() ? 0 : 1 ] = *F;
02533 }
02534 fIt = helper.GetAncestors(E2, *_mesh, TopAbs_FACE);
02535 while ( fIt->more() && FF2[1].IsNull())
02536 {
02537 const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
02538 if ( helper.IsSubShape( *F, data._solid))
02539 FF2[ FF2[0].IsNull() ? 0 : 1 ] = *F;
02540 }
02541
02542 if ( FF1[0].IsSame( FF2[0]) || FF1[0].IsSame( FF2[1]))
02543 std::swap( FF1[0], FF1[1] );
02544 if ( FF2[0].IsSame( FF1[0]) )
02545 std::swap( FF2[0], FF2[1] );
02546 if ( FF1[0].IsNull() || FF2[0].IsNull() )
02547 continue;
02548
02549
02550 bool ok;
02551 gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal;
02552 if ( edge1->_cosin < 0 )
02553 dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized();
02554 if ( edge2->_cosin < 0 )
02555 dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized();
02556
02557
02558
02559
02560
02561
02562
02563 double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
02564 double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 );
02565 gp_Vec newNorm = wgt1 * dir1 + wgt2 * dir2;
02566 newNorm.Normalize();
02567
02568 edge1->_normal = newNorm.XYZ();
02569
02570
02571 const SMDS_MeshNode *n1, *n2;
02572 n1 = edge1->_2neibors->_edges[0]->_nodes[0];
02573 n2 = edge1->_2neibors->_edges[1]->_nodes[0];
02574
02575
02576 edge1->SetDataByNeighbors( n1, n2, helper );
02577 gp_Vec dirInFace;
02578 if ( edge1->_cosin < 0 )
02579 dirInFace = dir1;
02580 else
02581 getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
02582 double angle = dir1.Angle( edge1->_normal );
02583 edge1->SetCosin( cos( angle ));
02584
02585
02586 if ( edge1->_cosin > 0.1 )
02587 {
02588 SMDS_ElemIteratorPtr fIt = edge1->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
02589 while ( fIt->more() )
02590 limitStepSize( data, fIt->next(), edge1->_cosin );
02591 }
02592
02593 edge1->InvalidateStep( 1 );
02594 edge1->_len = 0;
02595 edge1->SetNewLength( data._stepSize, helper );
02596 }
02597
02598
02599
02600
02601 for ( e2ee = edge2CloseEdge.begin(); e2ee != edge2CloseEdge.end(); ++e2ee )
02602 {
02603 _LayerEdge* edge1 = e2ee->first;
02604 if ( !edge1->_2neibors )
02605 continue;
02606 for ( int j = 0; j < 2; ++j )
02607 {
02608 _LayerEdge* neighbor = edge1->_2neibors->_edges[j];
02609 if ( edge2CloseEdge.count ( neighbor ))
02610 continue;
02611 _LayerEdge* prevEdge = edge1;
02612 const int nbSteps = 6;
02613 for ( int step = nbSteps; step; --step )
02614 {
02615 if ( !neighbor->_2neibors )
02616 break;
02617 int iNext = 0;
02618 _LayerEdge* nextEdge = neighbor->_2neibors->_edges[iNext];
02619 if ( nextEdge == prevEdge )
02620 nextEdge = neighbor->_2neibors->_edges[ ++iNext ];
02621
02622
02623 double r = double(step-1)/nbSteps;
02624 if ( !nextEdge->_2neibors )
02625 r = 0.5;
02626
02627 gp_XYZ newNorm = prevEdge->_normal * r + nextEdge->_normal * (1-r);
02628 newNorm.Normalize();
02629
02630 neighbor->_normal = newNorm;
02631 neighbor->SetCosin( prevEdge->_cosin * r + nextEdge->_cosin * (1-r) );
02632 neighbor->SetDataByNeighbors( prevEdge->_nodes[0], nextEdge->_nodes[0], helper );
02633
02634 neighbor->InvalidateStep( 1 );
02635 neighbor->_len = 0;
02636 neighbor->SetNewLength( data._stepSize, helper );
02637
02638
02639 prevEdge = neighbor;
02640 neighbor = nextEdge;
02641 }
02642 }
02643 }
02644 dumpFunctionEnd();
02645 }
02646
02647
02648
02649 for ( unsigned i = 0 ; i < tmpFaces.size(); ++i )
02650 delete tmpFaces[i];
02651
02652 return true;
02653 }
02654
02655
02660
02661
02662 bool _LayerEdge::FindIntersection( SMESH_ElementSearcher& searcher,
02663 double & distance,
02664 const double& epsilon,
02665 const SMDS_MeshElement** face)
02666 {
02667 vector< const SMDS_MeshElement* > suspectFaces;
02668 double segLen;
02669 gp_Ax1 lastSegment = LastSegment(segLen);
02670 searcher.GetElementsNearLine( lastSegment, SMDSAbs_Face, suspectFaces );
02671
02672 bool segmentIntersected = false;
02673 distance = Precision::Infinite();
02674 int iFace = -1;
02675 for ( unsigned j = 0 ; j < suspectFaces.size() && !segmentIntersected; ++j )
02676 {
02677 const SMDS_MeshElement* face = suspectFaces[j];
02678 if ( face->GetNodeIndex( _nodes.back() ) >= 0 ||
02679 face->GetNodeIndex( _nodes[0] ) >= 0 )
02680 continue;
02681 const int nbNodes = face->NbCornerNodes();
02682 bool intFound = false;
02683 double dist;
02684 SMDS_MeshElement::iterator nIt = face->begin_nodes();
02685 if ( nbNodes == 3 )
02686 {
02687 intFound = SegTriaInter( lastSegment, *nIt++, *nIt++, *nIt++, dist, epsilon );
02688 }
02689 else
02690 {
02691 const SMDS_MeshNode* tria[3];
02692 tria[0] = *nIt++;
02693 tria[1] = *nIt++;;
02694 for ( int n2 = 2; n2 < nbNodes && !intFound; ++n2 )
02695 {
02696 tria[2] = *nIt++;
02697 intFound = SegTriaInter(lastSegment, tria[0], tria[1], tria[2], dist, epsilon );
02698 tria[1] = tria[2];
02699 }
02700 }
02701 if ( intFound )
02702 {
02703 if ( dist < segLen*(1.01))
02704 segmentIntersected = true;
02705 if ( distance > dist )
02706 distance = dist, iFace = j;
02707 }
02708 }
02709 if ( iFace != -1 && face ) *face = suspectFaces[iFace];
02710
02711
02712
02713
02714
02715
02716
02717
02718 if ( segmentIntersected )
02719 {
02720 #ifdef __myDEBUG
02721 SMDS_MeshElement::iterator nIt = suspectFaces[iFace]->begin_nodes();
02722 gp_XYZ intP( lastSegment.Location().XYZ() + lastSegment.Direction().XYZ() * distance );
02723 cout << "nodes: tgt " << _nodes.back()->GetID() << " src " << _nodes[0]->GetID()
02724 << ", intersection with face ("
02725 << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
02726 << ") at point (" << intP.X() << ", " << intP.Y() << ", " << intP.Z()
02727 << ") distance = " << distance - segLen<< endl;
02728 #endif
02729 }
02730
02731 distance -= segLen;
02732
02733 return segmentIntersected;
02734 }
02735
02736
02740
02741
02742 gp_Ax1 _LayerEdge::LastSegment(double& segLen) const
02743 {
02744
02745 gp_XYZ orig = _pos.back();
02746 gp_XYZ dir;
02747 int iPrev = _pos.size() - 2;
02748 while ( iPrev >= 0 )
02749 {
02750 dir = orig - _pos[iPrev];
02751 if ( dir.SquareModulus() > 1e-100 )
02752 break;
02753 else
02754 iPrev--;
02755 }
02756
02757
02758 gp_Ax1 segDir;
02759 if ( iPrev < 0 )
02760 {
02761 segDir.SetLocation( SMESH_TNodeXYZ( _nodes[0] ));
02762 segDir.SetDirection( _normal );
02763 segLen = 0;
02764 }
02765 else
02766 {
02767 gp_Pnt pPrev = _pos[ iPrev ];
02768 if ( !_sWOL.IsNull() )
02769 {
02770 TopLoc_Location loc;
02771 if ( _sWOL.ShapeType() == TopAbs_EDGE )
02772 {
02773 double f,l;
02774 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
02775 pPrev = curve->Value( pPrev.X() ).Transformed( loc );
02776 }
02777 else
02778 {
02779 Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
02780 pPrev = surface->Value( pPrev.X(), pPrev.Y() ).Transformed( loc );
02781 }
02782 dir = SMESH_TNodeXYZ( _nodes.back() ) - pPrev.XYZ();
02783 }
02784 segDir.SetLocation( pPrev );
02785 segDir.SetDirection( dir );
02786 segLen = dir.Modulus();
02787 }
02788
02789 return segDir;
02790 }
02791
02792
02798
02799
02800 bool _LayerEdge::SegTriaInter( const gp_Ax1& lastSegment,
02801 const SMDS_MeshNode* n0,
02802 const SMDS_MeshNode* n1,
02803 const SMDS_MeshNode* n2,
02804 double& t,
02805 const double& EPSILON) const
02806 {
02807
02808
02809 gp_XYZ orig = lastSegment.Location().XYZ();
02810 gp_XYZ dir = lastSegment.Direction().XYZ();
02811
02812 SMESH_TNodeXYZ vert0( n0 );
02813 SMESH_TNodeXYZ vert1( n1 );
02814 SMESH_TNodeXYZ vert2( n2 );
02815
02816
02817 gp_XYZ tvec = orig - vert0;
02818
02819 if ( tvec * dir > EPSILON )
02820
02821 return false;
02822
02823 gp_XYZ edge1 = vert1 - vert0;
02824 gp_XYZ edge2 = vert2 - vert0;
02825
02826
02827 gp_XYZ pvec = dir ^ edge2;
02828
02829
02830 double det = edge1 * pvec;
02831
02832 if (det > -EPSILON && det < EPSILON)
02833 return 0;
02834 double inv_det = 1.0 / det;
02835
02836
02837 double u = ( tvec * pvec ) * inv_det;
02838 if (u < 0.0 || u > 1.0)
02839 return 0;
02840
02841
02842 gp_XYZ qvec = tvec ^ edge1;
02843
02844
02845 double v = (dir * qvec) * inv_det;
02846 if ( v < 0.0 || u + v > 1.0 )
02847 return 0;
02848
02849
02850 t = (edge2 * qvec) * inv_det;
02851
02852
02853
02854
02855
02856
02857
02858
02859
02860
02861
02862
02863
02864
02865
02866
02867
02868
02869
02870
02871
02872
02873
02874
02875
02876
02877
02878 return true;
02879 }
02880
02881
02886
02887
02888 bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface,
02889 const TopoDS_Face& F,
02890 SMESH_MesherHelper& helper)
02891 {
02892 ASSERT( IsOnEdge() );
02893
02894 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( _nodes.back() );
02895 SMESH_TNodeXYZ oldPos( tgtNode );
02896 double dist01, distNewOld;
02897
02898 SMESH_TNodeXYZ p0( _2neibors->_nodes[0]);
02899 SMESH_TNodeXYZ p1( _2neibors->_nodes[1]);
02900 dist01 = p0.Distance( _2neibors->_nodes[1] );
02901
02902 gp_Pnt newPos = p0 * _2neibors->_wgt[0] + p1 * _2neibors->_wgt[1];
02903 double lenDelta = 0;
02904 if ( _curvature )
02905 {
02906 lenDelta = _curvature->lenDelta( _len );
02907 newPos.ChangeCoord() += _normal * lenDelta;
02908 }
02909
02910 distNewOld = newPos.Distance( oldPos );
02911
02912 if ( F.IsNull() )
02913 {
02914 if ( _2neibors->_plnNorm )
02915 {
02916
02917 gp_XYZ new2src = SMESH_TNodeXYZ( _nodes[0] ) - newPos.XYZ();
02918 double new2srcProj = (*_2neibors->_plnNorm) * new2src;
02919 newPos.ChangeCoord() += (*_2neibors->_plnNorm) * new2srcProj;
02920 }
02921 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
02922 _pos.back() = newPos.XYZ();
02923 }
02924 else
02925 {
02926 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
02927 gp_XY uv( Precision::Infinite(), 0 );
02928 helper.CheckNodeUV( F, tgtNode, uv, 1e-10, true );
02929 _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
02930
02931 newPos = surface->Value( uv.X(), uv.Y() );
02932 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
02933 }
02934
02935 if ( _curvature && lenDelta < 0 )
02936 {
02937 gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
02938 _len -= prevPos.Distance( oldPos );
02939 _len += prevPos.Distance( newPos );
02940 }
02941 bool moved = distNewOld > dist01/50;
02942
02943 dumpMove( tgtNode );
02944
02945 return moved;
02946 }
02947
02948
02953
02954
02955 bool _LayerEdge::Smooth(int& badNb)
02956 {
02957 if ( _simplices.size() < 2 )
02958 return false;
02959
02960
02961 gp_XYZ newPos (0,0,0);
02962 for ( unsigned i = 0; i < _simplices.size(); ++i )
02963 newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev );
02964 newPos /= _simplices.size();
02965
02966 if ( _curvature )
02967 newPos += _normal * _curvature->lenDelta( _len );
02968
02969 gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
02970
02971
02972
02973
02974
02975
02976
02977
02978
02979
02980
02981
02982
02983 int nbOkBefore = 0;
02984 SMESH_TNodeXYZ tgtXYZ( _nodes.back() );
02985 for ( unsigned i = 0; i < _simplices.size(); ++i )
02986 nbOkBefore += _simplices[i].IsForward( _nodes[0], &tgtXYZ );
02987
02988 int nbOkAfter = 0;
02989 for ( unsigned i = 0; i < _simplices.size(); ++i )
02990 nbOkAfter += _simplices[i].IsForward( _nodes[0], &newPos );
02991
02992 if ( nbOkAfter < nbOkBefore )
02993 return false;
02994
02995 SMDS_MeshNode* n = const_cast< SMDS_MeshNode* >( _nodes.back() );
02996
02997 _len -= prevPos.Distance(SMESH_TNodeXYZ( n ));
02998 _len += prevPos.Distance(newPos);
02999
03000 n->setXYZ( newPos.X(), newPos.Y(), newPos.Z());
03001 _pos.back() = newPos;
03002
03003 badNb += _simplices.size() - nbOkAfter;
03004
03005 dumpMove( n );
03006
03007 return true;
03008 }
03009
03010
03014
03015
03016 void _LayerEdge::SetNewLength( double len, SMESH_MesherHelper& helper )
03017 {
03018 if ( _len - len > -1e-6 )
03019 {
03020 _pos.push_back( _pos.back() );
03021 return;
03022 }
03023
03024 SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
03025 SMESH_TNodeXYZ oldXYZ( n );
03026 gp_XYZ nXYZ = oldXYZ + _normal * ( len - _len ) * _lenFactor;
03027 n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
03028
03029 _pos.push_back( nXYZ );
03030 _len = len;
03031 if ( !_sWOL.IsNull() )
03032 {
03033 double distXYZ[4];
03034 if ( _sWOL.ShapeType() == TopAbs_EDGE )
03035 {
03036 double u = Precision::Infinite();
03037 helper.CheckNodeU( TopoDS::Edge( _sWOL ), n, u, 1e-10, true, distXYZ );
03038 _pos.back().SetCoord( u, 0, 0 );
03039 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
03040 pos->SetUParameter( u );
03041 }
03042 else
03043 {
03044 gp_XY uv( Precision::Infinite(), 0 );
03045 helper.CheckNodeUV( TopoDS::Face( _sWOL ), n, uv, 1e-10, true, distXYZ );
03046 _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
03047 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
03048 pos->SetUParameter( uv.X() );
03049 pos->SetVParameter( uv.Y() );
03050 }
03051 n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]);
03052 }
03053 dumpMove( n );
03054 }
03055
03056
03060
03061
03062 void _LayerEdge::InvalidateStep( int curStep )
03063 {
03064 if ( _pos.size() > curStep )
03065 {
03066 _pos.resize( curStep );
03067 gp_Pnt nXYZ = _pos.back();
03068 SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
03069 if ( !_sWOL.IsNull() )
03070 {
03071 TopLoc_Location loc;
03072 if ( _sWOL.ShapeType() == TopAbs_EDGE )
03073 {
03074 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
03075 pos->SetUParameter( nXYZ.X() );
03076 double f,l;
03077 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
03078 nXYZ = curve->Value( nXYZ.X() ).Transformed( loc );
03079 }
03080 else
03081 {
03082 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
03083 pos->SetUParameter( nXYZ.X() );
03084 pos->SetVParameter( nXYZ.Y() );
03085 Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
03086 nXYZ = surface->Value( nXYZ.X(), nXYZ.Y() ).Transformed( loc );
03087 }
03088 }
03089 n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
03090 dumpMove( n );
03091 }
03092 }
03093
03094
03098
03099
03100 bool _ViscousBuilder::refine(_SolidData& data)
03101 {
03102 SMESH_MesherHelper helper( *_mesh );
03103 helper.SetSubShape( data._solid );
03104 helper.SetElementsOnShape(false);
03105
03106 Handle(Geom_Curve) curve;
03107 Handle(Geom_Surface) surface;
03108 TopoDS_Edge geomEdge;
03109 TopoDS_Face geomFace;
03110 TopLoc_Location loc;
03111 double f,l, u;
03112 gp_XY uv;
03113 bool isOnEdge;
03114
03115 for ( unsigned i = 0; i < data._edges.size(); ++i )
03116 {
03117 _LayerEdge& edge = *data._edges[i];
03118
03119
03120 vector< double > segLen( edge._pos.size() );
03121 segLen[0] = 0.0;
03122 for ( unsigned j = 1; j < edge._pos.size(); ++j )
03123 segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus();
03124
03125
03126 const SMDS_MeshNode* tgtNode = edge._nodes.back();
03127 if ( edge._nodes.size() == 2 )
03128 {
03129 edge._nodes.resize( data._hyp->GetNumberLayers() + 1, 0 );
03130 edge._nodes[1] = 0;
03131 edge._nodes.back() = tgtNode;
03132 }
03133 if ( !edge._sWOL.IsNull() )
03134 {
03135 isOnEdge = ( edge._sWOL.ShapeType() == TopAbs_EDGE );
03136
03137
03138 if ( isOnEdge )
03139 {
03140 geomEdge = TopoDS::Edge( edge._sWOL );
03141 curve = BRep_Tool::Curve( geomEdge, loc, f,l);
03142
03143
03144 }
03145 else
03146 {
03147 geomFace = TopoDS::Face( edge._sWOL );
03148 surface = BRep_Tool::Surface( geomFace, loc );
03149
03150
03151 }
03152
03153
03154 }
03155
03156 double h0;
03157 const double T = segLen.back();
03158 const double f = data._hyp->GetStretchFactor();
03159 const int N = data._hyp->GetNumberLayers();
03160 const double fPowN = pow( f, N );
03161 if ( fPowN - 1 <= numeric_limits<double>::min() )
03162 h0 = T / N;
03163 else
03164 h0 = T * ( f - 1 )/( fPowN - 1 );
03165
03166 const double zeroLen = std::numeric_limits<double>::min();
03167
03168
03169 double hSum = 0, hi = h0/f;
03170 unsigned iSeg = 1;
03171 for ( unsigned iStep = 1; iStep < edge._nodes.size(); ++iStep )
03172 {
03173
03174 hi *= f;
03175 hSum += hi;
03176 while ( hSum > segLen[iSeg] && iSeg < segLen.size()-1)
03177 ++iSeg;
03178 int iPrevSeg = iSeg-1;
03179 while ( fabs( segLen[iPrevSeg] - segLen[iSeg]) <= zeroLen && iPrevSeg > 0 )
03180 --iPrevSeg;
03181 double r = ( segLen[iSeg] - hSum ) / ( segLen[iSeg] - segLen[iPrevSeg] );
03182 gp_Pnt pos = r * edge._pos[iPrevSeg] + (1-r) * edge._pos[iSeg];
03183
03184 SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >(edge._nodes[ iStep ]);
03185 if ( !edge._sWOL.IsNull() )
03186 {
03187
03188 if ( isOnEdge )
03189 {
03190 u = pos.X();
03191 pos = curve->Value( u ).Transformed(loc);
03192 }
03193 else
03194 {
03195 uv.SetCoord( pos.X(), pos.Y() );
03196 pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc);
03197 }
03198 }
03199
03200 if ( !node )
03201 {
03202 node = helper.AddNode( pos.X(), pos.Y(), pos.Z());
03203 if ( !edge._sWOL.IsNull() )
03204 {
03205 if ( isOnEdge )
03206 getMeshDS()->SetNodeOnEdge( node, geomEdge, u );
03207 else
03208 getMeshDS()->SetNodeOnFace( node, geomFace, uv.X(), uv.Y() );
03209 }
03210 else
03211 {
03212 getMeshDS()->SetNodeInVolume( node, helper.GetSubShapeID() );
03213 }
03214 }
03215 else
03216 {
03217 if ( !edge._sWOL.IsNull() )
03218 {
03219
03220 if ( isOnEdge )
03221 {
03222 u = 0.5 * ( u + helper.GetNodeU( geomEdge, node ));
03223 pos = curve->Value( u ).Transformed(loc);
03224 }
03225 else
03226 {
03227 uv = 0.5 * ( uv + helper.GetNodeUV( geomFace, node ));
03228 pos = surface->Value( uv.X(), uv.Y()).Transformed(loc);
03229 }
03230 }
03231 node->setXYZ( pos.X(), pos.Y(), pos.Z() );
03232 }
03233 }
03234 }
03235
03236
03237
03238 helper.SetElementsOnShape(true);
03239
03240 TopExp_Explorer exp( data._solid, TopAbs_FACE );
03241 for ( ; exp.More(); exp.Next() )
03242 {
03243 if ( _ignoreShapeIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
03244 continue;
03245 SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
03246 SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
03247 vector< vector<const SMDS_MeshNode*>* > nnVec;
03248 while ( fIt->more() )
03249 {
03250 const SMDS_MeshElement* face = fIt->next();
03251 int nbNodes = face->NbCornerNodes();
03252 nnVec.resize( nbNodes );
03253 SMDS_ElemIteratorPtr nIt = face->nodesIterator();
03254 for ( int iN = 0; iN < nbNodes; ++iN )
03255 {
03256 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
03257 nnVec[ iN ] = & data._n2eMap[ n ]->_nodes;
03258 }
03259
03260 int nbZ = nnVec[0]->size();
03261 switch ( nbNodes )
03262 {
03263 case 3:
03264 for ( int iZ = 1; iZ < nbZ; ++iZ )
03265 helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
03266 (*nnVec[0])[iZ], (*nnVec[1])[iZ], (*nnVec[2])[iZ]);
03267 break;
03268 case 4:
03269 for ( int iZ = 1; iZ < nbZ; ++iZ )
03270 helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1],
03271 (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1],
03272 (*nnVec[0])[iZ], (*nnVec[1])[iZ],
03273 (*nnVec[2])[iZ], (*nnVec[3])[iZ]);
03274 break;
03275 default:
03276 return error("Not supported type of element", data._index);
03277 }
03278 }
03279 }
03280 return true;
03281 }
03282
03283
03287
03288
03289 bool _ViscousBuilder::shrink()
03290 {
03291
03292
03293 map< TGeomID, _SolidData* > f2sdMap;
03294 for ( unsigned i = 0 ; i < _sdVec.size(); ++i )
03295 {
03296 _SolidData& data = _sdVec[i];
03297 TopTools_MapOfShape FFMap;
03298 map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
03299 for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
03300 if ( s2s->second.ShapeType() == TopAbs_FACE )
03301 {
03302 f2sdMap.insert( make_pair( getMeshDS()->ShapeToIndex( s2s->second ), &data ));
03303
03304 if ( FFMap.Add( (*s2s).second ))
03305
03306
03307
03308 if ( SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( s2s->second ))
03309 {
03310 SMESH_ProxyMesh::SubMesh* proxySub =
03311 data._proxyMesh->getFaceSubM( TopoDS::Face( s2s->second ), true);
03312 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
03313 while ( fIt->more() )
03314 proxySub->AddElement( fIt->next() );
03315 }
03316 }
03317 }
03318
03319 SMESH_MesherHelper helper( *_mesh );
03320
03321
03322 map< int, _Shrinker1D > e2shrMap;
03323
03324
03325 map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin();
03326 for ( ; f2sd != f2sdMap.end(); ++f2sd )
03327 {
03328 _SolidData& data = *f2sd->second;
03329 TNode2Edge& n2eMap = data._n2eMap;
03330 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first ));
03331 const bool reverse = ( data._reversedFaceIds.count( f2sd->first ));
03332
03333 Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
03334
03335 SMESH_subMesh* sm = _mesh->GetSubMesh( F );
03336 SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
03337
03338 helper.SetSubShape(F);
03339
03340
03341
03342
03343
03344
03345
03346 vector < const SMDS_MeshNode* > smoothNodes;
03347 {
03348 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
03349 while ( nIt->more() )
03350 {
03351 const SMDS_MeshNode* n = nIt->next();
03352 if ( n->NbInverseElements( SMDSAbs_Face ) > 0 )
03353 smoothNodes.push_back( n );
03354 }
03355 }
03356
03357 double refSign = 1;
03358 const set<TGeomID> ignoreShapes;
03359 if ( !smoothNodes.empty() )
03360 {
03361 gp_XY uv = helper.GetNodeUV( F, smoothNodes[0] );
03362 vector<_Simplex> simplices;
03363 getSimplices( smoothNodes[0], simplices, ignoreShapes );
03364 if ( simplices[0].IsForward(uv, F, helper,refSign) != (!reverse))
03365 refSign = -1;
03366 }
03367
03368
03369 vector< _LayerEdge* > lEdges;
03370 {
03371 SMESH_subMeshIteratorPtr subIt =
03372 sm->getDependsOnIterator(false, false);
03373 while ( subIt->more() )
03374 {
03375 SMESH_subMesh* sub = subIt->next();
03376 SMESHDS_SubMesh* subDS = sub->GetSubMeshDS();
03377 if ( subDS->NbNodes() == 0 || !n2eMap.count( subDS->GetNodes()->next() ))
03378 continue;
03379 SMDS_NodeIteratorPtr nIt = subDS->GetNodes();
03380 while ( nIt->more() )
03381 {
03382 _LayerEdge* edge = n2eMap[ nIt->next() ];
03383 lEdges.push_back( edge );
03384 prepareEdgeToShrink( *edge, F, helper, smDS );
03385 }
03386 }
03387 }
03388
03389
03390 const SMDS_MeshNode* nodes[20];
03391 for ( unsigned i = 0; i < lEdges.size(); ++i )
03392 {
03393 _LayerEdge& edge = *lEdges[i];
03394 const SMDS_MeshNode* srcNode = edge._nodes[0];
03395 const SMDS_MeshNode* tgtNode = edge._nodes.back();
03396 SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
03397 while ( fIt->more() )
03398 {
03399 const SMDS_MeshElement* f = fIt->next();
03400 if ( !smDS->Contains( f ))
03401 continue;
03402 SMDS_ElemIteratorPtr nIt = f->nodesIterator();
03403 for ( int iN = 0; iN < f->NbNodes(); ++iN )
03404 {
03405 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
03406 nodes[iN] = ( n == srcNode ? tgtNode : n );
03407 }
03408 helper.GetMeshDS()->ChangeElementNodes( f, nodes, f->NbNodes() );
03409 }
03410 }
03411
03412
03413 vector< _SmoothNode > nodesToSmooth( smoothNodes.size() );
03414 {
03415 dumpFunction(SMESH_Comment("beforeShrinkFace")<<f2sd->first);
03416 for ( unsigned i = 0; i < smoothNodes.size(); ++i )
03417 {
03418 const SMDS_MeshNode* n = smoothNodes[i];
03419 nodesToSmooth[ i ]._node = n;
03420
03421 getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes );
03422 dumpMove( n );
03423 }
03424 dumpFunctionEnd();
03425 }
03426
03427
03428
03429 set< _Shrinker1D* > eShri1D;
03430 {
03431 for ( unsigned i = 0; i < lEdges.size(); ++i )
03432 {
03433 _LayerEdge* edge = lEdges[i];
03434 if ( edge->_sWOL.ShapeType() == TopAbs_EDGE )
03435 {
03436 TGeomID edgeIndex = getMeshDS()->ShapeToIndex( edge->_sWOL );
03437 _Shrinker1D& srinker = e2shrMap[ edgeIndex ];
03438 eShri1D.insert( & srinker );
03439 srinker.AddEdge( edge, helper );
03440
03441
03442 srinker.RestoreParams();
03443 }
03444 }
03445 }
03446
03447
03448
03449
03450
03451 bool shrinked = true;
03452 int badNb, shriStep=0, smooStep=0;
03453 while ( shrinked )
03454 {
03455
03456
03457 dumpFunction(SMESH_Comment("moveBoundaryOnF")<<f2sd->first<<"_st"<<shriStep++ );
03458 shrinked = false;
03459 for ( unsigned i = 0; i < lEdges.size(); ++i )
03460 {
03461 shrinked |= lEdges[i]->SetNewLength2d( surface,F,helper );
03462 }
03463 dumpFunctionEnd();
03464 if ( !shrinked )
03465 break;
03466
03467
03468 set< _Shrinker1D* >::iterator shr = eShri1D.begin();
03469 for ( ; shr != eShri1D.end(); ++shr )
03470 (*shr)->Compute( false, helper );
03471
03472
03473
03474 int nbNoImpSteps = 0;
03475 bool moved = true;
03476 badNb = 1;
03477 while (( nbNoImpSteps < 5 && badNb > 0) && moved)
03478 {
03479 dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep);
03480
03481 int oldBadNb = badNb;
03482 badNb = 0;
03483 moved = false;
03484 for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
03485 {
03486 moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,false );
03487 }
03488 if ( badNb < oldBadNb )
03489 nbNoImpSteps = 0;
03490 else
03491 nbNoImpSteps++;
03492
03493 dumpFunctionEnd();
03494 }
03495 if ( badNb > 0 )
03496 return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first );
03497 }
03498
03499
03500 bool highQuality;
03501 {
03502 const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles();
03503 if ( hasTria != hasQuad )
03504 {
03505 highQuality = hasQuad;
03506 }
03507 else
03508 {
03509 set<int> nbNodesSet;
03510 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
03511 while ( fIt->more() && nbNodesSet.size() < 2 )
03512 nbNodesSet.insert( fIt->next()->NbCornerNodes() );
03513 highQuality = ( *nbNodesSet.begin() == 4 );
03514 }
03515 }
03516 for ( int st = highQuality ? 8 : 3; st; --st )
03517 {
03518 dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep);
03519 for ( unsigned i = 0; i < nodesToSmooth.size(); ++i )
03520 nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,st==1 );
03521 dumpFunctionEnd();
03522 }
03523
03524 _SrinkShapeListener::ToClearSubMeshWithSolid( sm, data._solid );
03525
03526 }
03527
03528
03529
03530
03531 map< int, _Shrinker1D >::iterator e2shr = e2shrMap.begin();
03532 for ( ; e2shr != e2shrMap.end(); ++e2shr )
03533 e2shr->second.SwapSrcTgtNodes( getMeshDS() );
03534
03535 return true;
03536 }
03537
03538
03542
03543
03544 bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge,
03545 const TopoDS_Face& F,
03546 SMESH_MesherHelper& helper,
03547 const SMESHDS_SubMesh* faceSubMesh)
03548 {
03549 const SMDS_MeshNode* srcNode = edge._nodes[0];
03550 const SMDS_MeshNode* tgtNode = edge._nodes.back();
03551
03552 edge._pos.clear();
03553
03554 if ( edge._sWOL.ShapeType() == TopAbs_FACE )
03555 {
03556 gp_XY srcUV = helper.GetNodeUV( F, srcNode );
03557 gp_XY tgtUV = helper.GetNodeUV( F, tgtNode );
03558 gp_Vec2d uvDir( srcUV, tgtUV );
03559 double uvLen = uvDir.Magnitude();
03560 uvDir /= uvLen;
03561 edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0);
03562
03563
03564 vector<const SMDS_MeshElement*> faces;
03565 multimap< double, const SMDS_MeshNode* > proj2node;
03566 SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
03567 while ( fIt->more() )
03568 {
03569 const SMDS_MeshElement* f = fIt->next();
03570 if ( faceSubMesh->Contains( f ))
03571 faces.push_back( f );
03572 }
03573 for ( unsigned i = 0; i < faces.size(); ++i )
03574 {
03575 const int nbNodes = faces[i]->NbCornerNodes();
03576 for ( int j = 0; j < nbNodes; ++j )
03577 {
03578 const SMDS_MeshNode* n = faces[i]->GetNode(j);
03579 if ( n == srcNode ) continue;
03580 if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
03581 ( faces.size() > 1 || nbNodes > 3 ))
03582 continue;
03583 gp_Pnt2d uv = helper.GetNodeUV( F, n );
03584 gp_Vec2d uvDirN( srcUV, uv );
03585 double proj = uvDirN * uvDir;
03586 proj2node.insert( make_pair( proj, n ));
03587 }
03588 }
03589
03590 multimap< double, const SMDS_MeshNode* >::iterator p2n = proj2node.begin(), p2nEnd;
03591 const double minProj = p2n->first;
03592 const double projThreshold = 1.1 * uvLen;
03593 if ( minProj > projThreshold )
03594 {
03595
03596 return true;
03597 }
03598 edge._pos.resize(1);
03599 edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 );
03600
03601
03602 p2nEnd = proj2node.lower_bound( projThreshold );
03603 int nbSimpl = ( std::distance( p2n, p2nEnd ) + 1) / 2;
03604 edge._simplices.resize( nbSimpl );
03605 for ( int i = 0; i < nbSimpl; ++i )
03606 {
03607 edge._simplices[i]._nPrev = p2n->second;
03608 if ( ++p2n != p2nEnd )
03609 edge._simplices[i]._nNext = p2n->second;
03610 }
03611
03612 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
03613 pos->SetUParameter( srcUV.X() );
03614 pos->SetVParameter( srcUV.Y() );
03615 }
03616 else
03617 {
03618 TopoDS_Edge E = TopoDS::Edge( edge._sWOL);
03619 SMESHDS_SubMesh* edgeSM = getMeshDS()->MeshElements( E );
03620 if ( !edgeSM || edgeSM->NbElements() == 0 )
03621 return error(SMESH_Comment("Not meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
03622
03623 const SMDS_MeshNode* n2 = 0;
03624 SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
03625 while ( eIt->more() && !n2 )
03626 {
03627 const SMDS_MeshElement* e = eIt->next();
03628 if ( !edgeSM->Contains(e)) continue;
03629 n2 = e->GetNode( 0 );
03630 if ( n2 == srcNode ) n2 = e->GetNode( 1 );
03631 }
03632 if ( !n2 )
03633 return error(SMESH_Comment("Wrongly meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
03634
03635 double uSrc = helper.GetNodeU( E, srcNode, n2 );
03636 double uTgt = helper.GetNodeU( E, tgtNode, srcNode );
03637 double u2 = helper.GetNodeU( E, n2, srcNode );
03638
03639 if ( fabs( uSrc-uTgt ) < 0.99 * fabs( uSrc-u2 ))
03640 {
03641
03642 return true;
03643 }
03644 edge._pos.resize(1);
03645 edge._pos[0].SetCoord( U_TGT, uTgt );
03646 edge._pos[0].SetCoord( U_SRC, uSrc );
03647 edge._pos[0].SetCoord( LEN_TGT, fabs( uSrc-uTgt ));
03648
03649 edge._simplices.resize( 1 );
03650 edge._simplices[0]._nPrev = n2;
03651
03652
03653 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
03654 pos->SetUParameter( uSrc );
03655 }
03656 return true;
03657
03658
03662
03663
03664
03665
03666
03667
03668
03669
03670
03671
03672
03673
03674
03675
03676
03677
03678
03679
03680
03681
03682
03683
03684
03685
03686
03687
03688
03689
03690
03691
03692
03693
03694
03695
03696
03697
03698
03699
03700
03701
03702
03703
03704
03705
03706
03707
03708
03709
03710 }
03711
03712
03716
03717
03718 bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
03719 const TopoDS_Face& F,
03720 SMESH_MesherHelper& helper )
03721 {
03722 if ( _pos.empty() )
03723 return false;
03724
03725 SMDS_MeshNode* tgtNode = const_cast< SMDS_MeshNode*& >( _nodes.back() );
03726
03727 if ( _sWOL.ShapeType() == TopAbs_FACE )
03728 {
03729 gp_XY curUV = helper.GetNodeUV( F, tgtNode );
03730 gp_Pnt2d tgtUV( _pos[0].X(), _pos[0].Y());
03731 gp_Vec2d uvDir( _normal.X(), _normal.Y() );
03732 const double uvLen = tgtUV.Distance( curUV );
03733
03734
03735 const double kSafe = 0.8;
03736 const double minStepSize = uvLen / 10;
03737 double stepSize = uvLen;
03738 for ( unsigned i = 0; i < _simplices.size(); ++i )
03739 {
03740 const SMDS_MeshNode* nn[2] = { _simplices[i]._nPrev, _simplices[i]._nNext };
03741 for ( int j = 0; j < 2; ++j )
03742 if ( const SMDS_MeshNode* n = nn[j] )
03743 {
03744 gp_XY uv = helper.GetNodeUV( F, n );
03745 gp_Vec2d uvDirN( curUV, uv );
03746 double proj = uvDirN * uvDir * kSafe;
03747 if ( proj < stepSize && proj > minStepSize )
03748 stepSize = proj;
03749 }
03750 }
03751
03752 gp_Pnt2d newUV;
03753 if ( stepSize == uvLen )
03754 {
03755 newUV = tgtUV;
03756 _pos.clear();
03757 }
03758 else
03759 {
03760 newUV = curUV + uvDir.XY() * stepSize;
03761 }
03762
03763 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
03764 pos->SetUParameter( newUV.X() );
03765 pos->SetVParameter( newUV.Y() );
03766
03767 #ifdef __myDEBUG
03768 gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
03769 tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
03770 dumpMove( tgtNode );
03771 #endif
03772 }
03773 else
03774 {
03775 TopoDS_Edge E = TopoDS::Edge( _sWOL );
03776 const SMDS_MeshNode* n2 = _simplices[0]._nPrev;
03777
03778 const double u2 = helper.GetNodeU( E, n2, tgtNode );
03779 const double uSrc = _pos[0].Coord( U_SRC );
03780 const double lenTgt = _pos[0].Coord( LEN_TGT );
03781
03782 double newU = _pos[0].Coord( U_TGT );
03783 if ( lenTgt < 0.99 * fabs( uSrc-u2 ))
03784 {
03785 _pos.clear();
03786 }
03787 else
03788 {
03789 newU = 0.1 * uSrc + 0.9 * u2;
03790 }
03791 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
03792 pos->SetUParameter( newU );
03793 #ifdef __myDEBUG
03794 gp_XY newUV = helper.GetNodeUV( F, tgtNode, _nodes[0]);
03795 gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
03796 tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
03797 dumpMove( tgtNode );
03798 #endif
03799 }
03800 return true;
03801 }
03802
03803
03808
03809
03810 bool _SmoothNode::Smooth(int& badNb,
03811 Handle(Geom_Surface)& surface,
03812 SMESH_MesherHelper& helper,
03813 const double refSign,
03814 bool set3D)
03815 {
03816 const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() );
03817
03818
03819 gp_XY newPos (0,0);
03820 for ( unsigned i = 0; i < _simplices.size(); ++i )
03821 newPos += helper.GetNodeUV( face, _simplices[i]._nPrev );
03822 newPos /= _simplices.size();
03823
03824
03825 int nbOkBefore = 0;
03826 gp_XY tgtUV = helper.GetNodeUV( face, _node );
03827 for ( unsigned i = 0; i < _simplices.size(); ++i )
03828 nbOkBefore += _simplices[i].IsForward( tgtUV, face, helper, refSign );
03829
03830 int nbOkAfter = 0;
03831 for ( unsigned i = 0; i < _simplices.size(); ++i )
03832 nbOkAfter += _simplices[i].IsForward( newPos, face, helper, refSign );
03833
03834 if ( nbOkAfter < nbOkBefore )
03835 return false;
03836
03837 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( _node->GetPosition() );
03838 pos->SetUParameter( newPos.X() );
03839 pos->SetVParameter( newPos.Y() );
03840
03841 #ifdef __myDEBUG
03842 set3D = true;
03843 #endif
03844 if ( set3D )
03845 {
03846 gp_Pnt p = surface->Value( newPos.X(), newPos.Y() );
03847 const_cast< SMDS_MeshNode* >( _node )->setXYZ( p.X(), p.Y(), p.Z() );
03848 dumpMove( _node );
03849 }
03850
03851 badNb += _simplices.size() - nbOkAfter;
03852 return ( (tgtUV-newPos).SquareModulus() > 1e-10 );
03853 }
03854
03855
03859
03860
03861 _SolidData::~_SolidData()
03862 {
03863 for ( unsigned i = 0; i < _edges.size(); ++i )
03864 {
03865 if ( _edges[i] && _edges[i]->_2neibors )
03866 delete _edges[i]->_2neibors;
03867 delete _edges[i];
03868 }
03869 _edges.clear();
03870 }
03871
03875
03876
03877 void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper )
03878 {
03879
03880 if ( _nodes.empty() )
03881 {
03882 _edges[0] = _edges[1] = 0;
03883 _done = false;
03884 }
03885
03886 if ( e == _edges[0] || e == _edges[1] )
03887 return;
03888 if ( e->_sWOL.IsNull() || e->_sWOL.ShapeType() != TopAbs_EDGE )
03889 throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
03890 if ( _edges[0] && _edges[0]->_sWOL != e->_sWOL )
03891 throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
03892
03893
03894 const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
03895 double f,l;
03896 BRep_Tool::Range( E, f,l );
03897 double u = helper.GetNodeU( E, e->_nodes[0], e->_nodes.back());
03898 _edges[ u < 0.5*(f+l) ? 0 : 1 ] = e;
03899
03900
03901
03902 const SMDS_MeshNode* tgtNode0 = _edges[0] ? _edges[0]->_nodes.back() : 0;
03903 const SMDS_MeshNode* tgtNode1 = _edges[1] ? _edges[1]->_nodes.back() : 0;
03904
03905 if ( _nodes.empty() )
03906 {
03907 SMESHDS_SubMesh * eSubMesh = helper.GetMeshDS()->MeshElements( E );
03908 if ( !eSubMesh || eSubMesh->NbNodes() < 1 )
03909 return;
03910 TopLoc_Location loc;
03911 Handle(Geom_Curve) C = BRep_Tool::Curve(E, loc, f,l);
03912 GeomAdaptor_Curve aCurve(C);
03913 const double totLen = GCPnts_AbscissaPoint::Length(aCurve, f, l);
03914
03915 int nbExpectNodes = eSubMesh->NbNodes() - e->_nodes.size();
03916 _initU .reserve( nbExpectNodes );
03917 _normPar.reserve( nbExpectNodes );
03918 _nodes .reserve( nbExpectNodes );
03919 SMDS_NodeIteratorPtr nIt = eSubMesh->GetNodes();
03920 while ( nIt->more() )
03921 {
03922 const SMDS_MeshNode* node = nIt->next();
03923 if ( node->NbInverseElements(SMDSAbs_Edge) == 0 ||
03924 node == tgtNode0 || node == tgtNode1 )
03925 continue;
03926 _nodes.push_back( node );
03927 _initU.push_back( helper.GetNodeU( E, node ));
03928 double len = GCPnts_AbscissaPoint::Length(aCurve, f, _initU.back());
03929 _normPar.push_back( len / totLen );
03930 }
03931 }
03932 else
03933 {
03934
03935 int nbFound = 0;
03936 for ( unsigned i = 0; i < _nodes.size(); ++i )
03937 if ( !_nodes[i] || _nodes[i] == tgtNode0 || _nodes[i] == tgtNode1 )
03938 _nodes[i] = 0, nbFound++;
03939 if ( nbFound == _nodes.size() )
03940 _nodes.clear();
03941 }
03942 }
03943
03944
03948
03949
03950 void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
03951 {
03952 if ( _done || _nodes.empty())
03953 return;
03954 const _LayerEdge* e = _edges[0];
03955 if ( !e ) e = _edges[1];
03956 if ( !e ) return;
03957
03958 _done = (( !_edges[0] || _edges[0]->_pos.empty() ) &&
03959 ( !_edges[1] || _edges[1]->_pos.empty() ));
03960
03961 const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
03962 double f,l;
03963 if ( set3D || _done )
03964 {
03965 Handle(Geom_Curve) C = BRep_Tool::Curve(E, f,l);
03966 GeomAdaptor_Curve aCurve(C);
03967
03968 if ( _edges[0] )
03969 f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
03970 if ( _edges[1] )
03971 l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
03972 double totLen = GCPnts_AbscissaPoint::Length( aCurve, f, l );
03973
03974 for ( unsigned i = 0; i < _nodes.size(); ++i )
03975 {
03976 if ( !_nodes[i] ) continue;
03977 double len = totLen * _normPar[i];
03978 GCPnts_AbscissaPoint discret( aCurve, len, f );
03979 if ( !discret.IsDone() )
03980 return throw SALOME_Exception(LOCALIZED("GCPnts_AbscissaPoint failed"));
03981 double u = discret.Parameter();
03982 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
03983 pos->SetUParameter( u );
03984 gp_Pnt p = C->Value( u );
03985 const_cast< SMDS_MeshNode*>( _nodes[i] )->setXYZ( p.X(), p.Y(), p.Z() );
03986 }
03987 }
03988 else
03989 {
03990 BRep_Tool::Range( E, f,l );
03991 if ( _edges[0] )
03992 f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
03993 if ( _edges[1] )
03994 l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
03995
03996 for ( unsigned i = 0; i < _nodes.size(); ++i )
03997 {
03998 if ( !_nodes[i] ) continue;
03999 double u = f * ( 1-_normPar[i] ) + l * _normPar[i];
04000 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
04001 pos->SetUParameter( u );
04002 }
04003 }
04004 }
04005
04006
04010
04011
04012 void _Shrinker1D::RestoreParams()
04013 {
04014 if ( _done )
04015 for ( unsigned i = 0; i < _nodes.size(); ++i )
04016 {
04017 if ( !_nodes[i] ) continue;
04018 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
04019 pos->SetUParameter( _initU[i] );
04020 }
04021 _done = false;
04022 }
04023
04027
04028
04029 void _Shrinker1D::SwapSrcTgtNodes( SMESHDS_Mesh* mesh )
04030 {
04031 const SMDS_MeshNode* nodes[3];
04032 for ( int i = 0; i < 2; ++i )
04033 {
04034 if ( !_edges[i] ) continue;
04035
04036 SMESHDS_SubMesh * eSubMesh = mesh->MeshElements( _edges[i]->_sWOL );
04037 if ( !eSubMesh ) return;
04038 const SMDS_MeshNode* srcNode = _edges[i]->_nodes[0];
04039 const SMDS_MeshNode* tgtNode = _edges[i]->_nodes.back();
04040 SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
04041 while ( eIt->more() )
04042 {
04043 const SMDS_MeshElement* e = eIt->next();
04044 if ( !eSubMesh->Contains( e ))
04045 continue;
04046 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
04047 for ( int iN = 0; iN < e->NbNodes(); ++iN )
04048 {
04049 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
04050 nodes[iN] = ( n == srcNode ? tgtNode : n );
04051 }
04052 mesh->ChangeElementNodes( e, nodes, e->NbNodes() );
04053 }
04054 }
04055 }
04056
04057
04061
04062
04063 bool _ViscousBuilder::addBoundaryElements()
04064 {
04065 SMESH_MesherHelper helper( *_mesh );
04066
04067 for ( unsigned i = 0; i < _sdVec.size(); ++i )
04068 {
04069 _SolidData& data = _sdVec[i];
04070 TopTools_IndexedMapOfShape geomEdges;
04071 TopExp::MapShapes( data._solid, TopAbs_EDGE, geomEdges );
04072 for ( int iE = 1; iE <= geomEdges.Extent(); ++iE )
04073 {
04074 const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE));
04075
04076
04077
04078 map< double, const SMDS_MeshNode* > u2nodes;
04079 if ( !SMESH_Algo::GetSortedNodesOnEdge( getMeshDS(), E, false, u2nodes))
04080 continue;
04081
04082 vector< _LayerEdge* > ledges; ledges.reserve( u2nodes.size() );
04083 TNode2Edge & n2eMap = data._n2eMap;
04084 map< double, const SMDS_MeshNode* >::iterator u2n = u2nodes.begin();
04085 {
04086
04087 TNode2Edge::iterator n2e = n2eMap.find( u2n->second );
04088 if ( n2e == n2eMap.end() ) continue;
04089 ledges.push_back( n2e->second );
04090 u2n++;
04091 if (( n2e = n2eMap.find( u2n->second )) == n2eMap.end() )
04092 continue;
04093 ledges.push_back( n2eMap[ u2n->second ]);
04094
04095 const SMDS_MeshNode* tgtN0 = ledges[0]->_nodes.back();
04096 const SMDS_MeshNode* tgtN1 = ledges[1]->_nodes.back();
04097 int nbSharedPyram = 0;
04098 SMDS_ElemIteratorPtr vIt = tgtN0->GetInverseElementIterator(SMDSAbs_Volume);
04099 while ( vIt->more() )
04100 {
04101 const SMDS_MeshElement* v = vIt->next();
04102 nbSharedPyram += int( v->GetNodeIndex( tgtN1 ) >= 0 );
04103 }
04104 if ( nbSharedPyram > 1 )
04105 continue;
04106
04107 if ( getMeshDS()->FindFace( ledges[0]->_nodes[0], ledges[0]->_nodes[1],
04108 ledges[1]->_nodes[0], ledges[1]->_nodes[1]))
04109 continue;
04110 }
04111 for ( ++u2n; u2n != u2nodes.end(); ++u2n )
04112 ledges.push_back( n2eMap[ u2n->second ]);
04113
04114
04115
04116 bool reverse = false, isOnFace;
04117
04118 map< TGeomID, TopoDS_Shape >::iterator e2f =
04119 data._shrinkShape2Shape.find( getMeshDS()->ShapeToIndex( E ));
04120 TopoDS_Shape F;
04121 if (( isOnFace = ( e2f != data._shrinkShape2Shape.end() )))
04122 {
04123 F = e2f->second.Oriented( TopAbs_FORWARD );
04124 reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED );
04125 if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED )
04126 reverse = !reverse;
04127 }
04128 else
04129 {
04130
04131 PShapeIteratorPtr fIt = helper.GetAncestors( E, *_mesh, TopAbs_FACE );
04132 while ( fIt->more() && F.IsNull() )
04133 {
04134 const TopoDS_Shape* pF = fIt->next();
04135 if ( helper.IsSubShape( *pF, data._solid) &&
04136 !_ignoreShapeIds.count( e2f->first ))
04137 F = *pF;
04138 }
04139 }
04140
04141 SMESHDS_SubMesh* sm = 0;
04142 if ( isOnFace )
04143 sm = getMeshDS()->MeshElements( F );
04144 else
04145 sm = data._proxyMesh->getFaceSubM( TopoDS::Face(F), true );
04146 if ( !sm )
04147 return error("error in addBoundaryElements()", data._index);
04148
04149
04150 const int dj1 = reverse ? 0 : 1;
04151 const int dj2 = reverse ? 1 : 0;
04152 for ( unsigned j = 1; j < ledges.size(); ++j )
04153 {
04154 vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes;
04155 vector< const SMDS_MeshNode*>& nn2 = ledges[j-dj2]->_nodes;
04156 if ( isOnFace )
04157 for ( unsigned z = 1; z < nn1.size(); ++z )
04158 sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
04159 else
04160 for ( unsigned z = 1; z < nn1.size(); ++z )
04161 sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z]));
04162 }
04163 }
04164 }
04165
04166 return true;
04167 }