1 // Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 // File : StdMeshers_ViscousLayers.cxx
21 // Created : Wed Dec 1 15:15:34 2010
22 // Author : Edward AGAPOV (eap)
24 #include "StdMeshers_ViscousLayers.hxx"
26 #include "SMDS_EdgePosition.hxx"
27 #include "SMDS_FaceOfNodes.hxx"
28 #include "SMDS_FacePosition.hxx"
29 #include "SMDS_MeshNode.hxx"
30 #include "SMDS_SetIterator.hxx"
31 #include "SMESHDS_Group.hxx"
32 #include "SMESHDS_Hypothesis.hxx"
33 #include "SMESH_Algo.hxx"
34 #include "SMESH_ComputeError.hxx"
35 #include "SMESH_ControlsDef.hxx"
36 #include "SMESH_Gen.hxx"
37 #include "SMESH_Group.hxx"
38 #include "SMESH_HypoFilter.hxx"
39 #include "SMESH_Mesh.hxx"
40 #include "SMESH_MeshAlgos.hxx"
41 #include "SMESH_MesherHelper.hxx"
42 #include "SMESH_ProxyMesh.hxx"
43 #include "SMESH_subMesh.hxx"
44 #include "SMESH_subMeshEventListener.hxx"
45 #include "StdMeshers_FaceSide.hxx"
47 #include <BRepAdaptor_Curve2d.hxx>
48 #include <BRepAdaptor_Surface.hxx>
49 #include <BRepLProp_SLProps.hxx>
50 #include <BRep_Tool.hxx>
51 #include <Bnd_B2d.hxx>
52 #include <Bnd_B3d.hxx>
54 #include <GCPnts_AbscissaPoint.hxx>
55 #include <Geom2d_Circle.hxx>
56 #include <Geom2d_Line.hxx>
57 #include <Geom2d_TrimmedCurve.hxx>
58 #include <GeomAdaptor_Curve.hxx>
59 #include <GeomLib.hxx>
60 #include <Geom_Circle.hxx>
61 #include <Geom_Curve.hxx>
62 #include <Geom_Line.hxx>
63 #include <Geom_TrimmedCurve.hxx>
64 #include <Precision.hxx>
65 #include <Standard_ErrorHandler.hxx>
66 #include <Standard_Failure.hxx>
67 #include <TColStd_Array1OfReal.hxx>
69 #include <TopExp_Explorer.hxx>
70 #include <TopTools_IndexedMapOfShape.hxx>
71 #include <TopTools_ListOfShape.hxx>
72 #include <TopTools_MapOfShape.hxx>
74 #include <TopoDS_Edge.hxx>
75 #include <TopoDS_Face.hxx>
76 #include <TopoDS_Vertex.hxx>
90 //================================================================================
95 enum UIndex { U_TGT = 1, U_SRC, LEN_TGT };
97 const double theMinSmoothCosin = 0.1;
100 * \brief SMESH_ProxyMesh computed by _ViscousBuilder for a SOLID.
101 * It is stored in a SMESH_subMesh of the SOLID as SMESH_subMeshEventListenerData
103 struct _MeshOfSolid : public SMESH_ProxyMesh,
104 public SMESH_subMeshEventListenerData
106 bool _n2nMapComputed;
108 _MeshOfSolid( SMESH_Mesh* mesh)
109 :SMESH_subMeshEventListenerData( /*isDeletable=*/true),_n2nMapComputed(false)
111 SMESH_ProxyMesh::setMesh( *mesh );
114 // returns submesh for a geom face
115 SMESH_ProxyMesh::SubMesh* getFaceSubM(const TopoDS_Face& F, bool create=false)
117 TGeomID i = SMESH_ProxyMesh::shapeIndex(F);
118 return create ? SMESH_ProxyMesh::getProxySubMesh(i) : findProxySubMesh(i);
120 void setNode2Node(const SMDS_MeshNode* srcNode,
121 const SMDS_MeshNode* proxyNode,
122 const SMESH_ProxyMesh::SubMesh* subMesh)
124 SMESH_ProxyMesh::setNode2Node( srcNode,proxyNode,subMesh);
127 //--------------------------------------------------------------------------------
129 * \brief Listener of events of 3D sub-meshes computed with viscous layers.
130 * It is used to clear an inferior dim sub-meshes modified by viscous layers
132 class _ShrinkShapeListener : SMESH_subMeshEventListener
134 _ShrinkShapeListener()
135 : SMESH_subMeshEventListener(/*isDeletable=*/false,
136 "StdMeshers_ViscousLayers::_ShrinkShapeListener") {}
138 static SMESH_subMeshEventListener* Get() { static _ShrinkShapeListener l; return &l; }
139 virtual void ProcessEvent(const int event,
141 SMESH_subMesh* solidSM,
142 SMESH_subMeshEventListenerData* data,
143 const SMESH_Hypothesis* hyp)
145 if ( SMESH_subMesh::COMPUTE_EVENT == eventType && solidSM->IsEmpty() && data )
147 SMESH_subMeshEventListener::ProcessEvent(event,eventType,solidSM,data,hyp);
151 //--------------------------------------------------------------------------------
153 * \brief Listener of events of 3D sub-meshes computed with viscous layers.
154 * It is used to store data computed by _ViscousBuilder for a sub-mesh and to
155 * delete the data as soon as it has been used
157 class _ViscousListener : SMESH_subMeshEventListener
160 SMESH_subMeshEventListener(/*isDeletable=*/false,
161 "StdMeshers_ViscousLayers::_ViscousListener") {}
162 static SMESH_subMeshEventListener* Get() { static _ViscousListener l; return &l; }
164 virtual void ProcessEvent(const int event,
166 SMESH_subMesh* subMesh,
167 SMESH_subMeshEventListenerData* data,
168 const SMESH_Hypothesis* hyp)
170 if ( SMESH_subMesh::COMPUTE_EVENT == eventType )
172 // delete SMESH_ProxyMesh containing temporary faces
173 subMesh->DeleteEventListener( this );
176 // Finds or creates proxy mesh of the solid
177 static _MeshOfSolid* GetSolidMesh(SMESH_Mesh* mesh,
178 const TopoDS_Shape& solid,
181 if ( !mesh ) return 0;
182 SMESH_subMesh* sm = mesh->GetSubMesh(solid);
183 _MeshOfSolid* data = (_MeshOfSolid*) sm->GetEventListenerData( Get() );
184 if ( !data && toCreate )
186 data = new _MeshOfSolid(mesh);
187 data->mySubMeshes.push_back( sm ); // to find SOLID by _MeshOfSolid
188 sm->SetEventListener( Get(), data, sm );
192 // Removes proxy mesh of the solid
193 static void RemoveSolidMesh(SMESH_Mesh* mesh, const TopoDS_Shape& solid)
195 mesh->GetSubMesh(solid)->DeleteEventListener( _ViscousListener::Get() );
199 //================================================================================
201 * \brief sets a sub-mesh event listener to clear sub-meshes of sub-shapes of
202 * the main shape when sub-mesh of the main shape is cleared,
203 * for example to clear sub-meshes of FACEs when sub-mesh of a SOLID
206 //================================================================================
208 void ToClearSubWithMain( SMESH_subMesh* sub, const TopoDS_Shape& main)
210 SMESH_subMesh* mainSM = sub->GetFather()->GetSubMesh( main );
211 SMESH_subMeshEventListenerData* data =
212 mainSM->GetEventListenerData( _ShrinkShapeListener::Get());
215 if ( find( data->mySubMeshes.begin(), data->mySubMeshes.end(), sub ) ==
216 data->mySubMeshes.end())
217 data->mySubMeshes.push_back( sub );
221 data = SMESH_subMeshEventListenerData::MakeData( /*dependent=*/sub );
222 sub->SetEventListener( _ShrinkShapeListener::Get(), data, /*whereToListenTo=*/mainSM );
225 //--------------------------------------------------------------------------------
227 * \brief Simplex (triangle or tetrahedron) based on 1 (tria) or 2 (tet) nodes of
228 * _LayerEdge and 2 nodes of the mesh surface beening smoothed.
229 * The class is used to check validity of face or volumes around a smoothed node;
230 * it stores only 2 nodes as the other nodes are stored by _LayerEdge.
234 const SMDS_MeshNode *_nPrev, *_nNext; // nodes on a smoothed mesh surface
235 const SMDS_MeshNode *_nOpp; // in 2D case, a node opposite to a smoothed node in QUAD
236 _Simplex(const SMDS_MeshNode* nPrev=0,
237 const SMDS_MeshNode* nNext=0,
238 const SMDS_MeshNode* nOpp=0)
239 : _nPrev(nPrev), _nNext(nNext), _nOpp(nOpp) {}
240 bool IsForward(const SMDS_MeshNode* nSrc, const gp_XYZ* pntTgt) const
242 const double M[3][3] =
243 {{ _nNext->X() - nSrc->X(), _nNext->Y() - nSrc->Y(), _nNext->Z() - nSrc->Z() },
244 { pntTgt->X() - nSrc->X(), pntTgt->Y() - nSrc->Y(), pntTgt->Z() - nSrc->Z() },
245 { _nPrev->X() - nSrc->X(), _nPrev->Y() - nSrc->Y(), _nPrev->Z() - nSrc->Z() }};
246 double determinant = ( + M[0][0]*M[1][1]*M[2][2]
247 + M[0][1]*M[1][2]*M[2][0]
248 + M[0][2]*M[1][0]*M[2][1]
249 - M[0][0]*M[1][2]*M[2][1]
250 - M[0][1]*M[1][0]*M[2][2]
251 - M[0][2]*M[1][1]*M[2][0]);
252 return determinant > 1e-100;
254 bool IsForward(const gp_XY& tgtUV,
255 const SMDS_MeshNode* smoothedNode,
256 const TopoDS_Face& face,
257 SMESH_MesherHelper& helper,
258 const double refSign) const
260 gp_XY prevUV = helper.GetNodeUV( face, _nPrev, smoothedNode );
261 gp_XY nextUV = helper.GetNodeUV( face, _nNext, smoothedNode );
262 gp_Vec2d v1( tgtUV, prevUV ), v2( tgtUV, nextUV );
264 return d*refSign > 1e-100;
266 bool IsNeighbour(const _Simplex& other) const
268 return _nPrev == other._nNext || _nNext == other._nPrev;
271 //--------------------------------------------------------------------------------
273 * Structure used to take into account surface curvature while smoothing
278 double _k; // factor to correct node smoothed position
279 double _h2lenRatio; // avgNormProj / (2*avgDist)
281 static _Curvature* New( double avgNormProj, double avgDist )
284 if ( fabs( avgNormProj / avgDist ) > 1./200 )
287 c->_r = avgDist * avgDist / avgNormProj;
288 c->_k = avgDist * avgDist / c->_r / c->_r;
289 //c->_k = avgNormProj / c->_r;
290 c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive
291 c->_h2lenRatio = avgNormProj / ( avgDist + avgDist );
295 double lenDelta(double len) const { return _k * ( _r + len ); }
296 double lenDeltaByDist(double dist) const { return dist * _h2lenRatio; }
299 //--------------------------------------------------------------------------------
301 * \brief Edge normal to surface, connecting a node on solid surface (_nodes[0])
302 * and a node of the most internal layer (_nodes.back())
306 vector< const SMDS_MeshNode*> _nodes;
308 gp_XYZ _normal; // to solid surface
309 vector<gp_XYZ> _pos; // points computed during inflation
310 double _len; // length achived with the last inflation step
311 double _cosin; // of angle (_normal ^ surface)
312 double _lenFactor; // to compute _len taking _cosin into account
314 // face or edge w/o layer along or near which _LayerEdge is inflated
316 // simplices connected to the source node (_nodes[0]);
317 // used for smoothing and quality check of _LayerEdge's based on the FACE
318 vector<_Simplex> _simplices;
319 // data for smoothing of _LayerEdge's based on the EDGE
320 _2NearEdges* _2neibors;
322 _Curvature* _curvature;
323 // TODO:: detele _Curvature, _plnNorm
325 void SetNewLength( double len, SMESH_MesherHelper& helper );
326 bool SetNewLength2d( Handle(Geom_Surface)& surface,
327 const TopoDS_Face& F,
328 SMESH_MesherHelper& helper );
329 void SetDataByNeighbors( const SMDS_MeshNode* n1,
330 const SMDS_MeshNode* n2,
331 SMESH_MesherHelper& helper);
332 void InvalidateStep( int curStep, bool restoreLength=false );
333 bool Smooth(int& badNb);
334 bool SmoothOnEdge(Handle(Geom_Surface)& surface,
335 const TopoDS_Face& F,
336 SMESH_MesherHelper& helper);
337 bool FindIntersection( SMESH_ElementSearcher& searcher,
339 const double& epsilon,
340 const SMDS_MeshElement** face = 0);
341 bool SegTriaInter( const gp_Ax1& lastSegment,
342 const SMDS_MeshNode* n0,
343 const SMDS_MeshNode* n1,
344 const SMDS_MeshNode* n2,
346 const double& epsilon) const;
347 gp_Ax1 LastSegment(double& segLen) const;
348 gp_XY LastUV( const TopoDS_Face& F ) const;
349 bool IsOnEdge() const { return _2neibors; }
350 gp_XYZ Copy( _LayerEdge& other, SMESH_MesherHelper& helper );
351 void SetCosin( double cosin );
355 bool operator () (const _LayerEdge* e1, const _LayerEdge* e2) const
357 const bool cmpNodes = ( e1 && e2 && e1->_nodes.size() && e2->_nodes.size() );
358 return cmpNodes ? ( e1->_nodes[0]->GetID() < e2->_nodes[0]->GetID()) : ( e1 < e2 );
362 //--------------------------------------------------------------------------------
364 * Structure used to smooth a _LayerEdge based on an EDGE.
368 double _wgt [2]; // weights of _nodes
369 _LayerEdge* _edges[2];
371 // normal to plane passing through _LayerEdge._normal and tangent of EDGE
374 _2NearEdges() { _edges[0]=_edges[1]=0; _plnNorm = 0; }
375 const SMDS_MeshNode* tgtNode(bool is2nd) {
376 return _edges[is2nd] ? _edges[is2nd]->_nodes.back() : 0;
378 const SMDS_MeshNode* srcNode(bool is2nd) {
379 return _edges[is2nd] ? _edges[is2nd]->_nodes[0] : 0;
382 std::swap( _wgt [0], _wgt [1] );
383 std::swap( _edges[0], _edges[1] );
386 //--------------------------------------------------------------------------------
388 * \brief Convex FACE whose radius of curvature is less than the thickness of
389 * layers. It is used to detect distortion of prisms based on a convex
390 * FACE and to update normals to enable further increasing the thickness
396 // edges whose _simplices are used to detect prism destorsion
397 vector< _LayerEdge* > _simplexTestEdges;
399 // map a sub-shape to it's index in _SolidData::_endEdgeOnShape vector
400 map< TGeomID, int > _subIdToEdgeEnd;
404 bool GetCenterOfCurvature( _LayerEdge* ledge,
405 BRepLProp_SLProps& surfProp,
406 SMESH_MesherHelper& helper,
407 gp_Pnt & center ) const;
408 bool CheckPrisms() const;
411 //--------------------------------------------------------------------------------
413 typedef map< const SMDS_MeshNode*, _LayerEdge*, TIDCompare > TNode2Edge;
415 //--------------------------------------------------------------------------------
417 * \brief Data of a SOLID
422 const StdMeshers_ViscousLayers* _hyp;
423 TopoDS_Shape _hypShape;
424 _MeshOfSolid* _proxyMesh;
425 set<TGeomID> _reversedFaceIds;
426 set<TGeomID> _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDS
428 double _stepSize, _stepSizeCoeff;
429 const SMDS_MeshNode* _stepSizeNodes[2];
431 TNode2Edge _n2eMap; // nodes and _LayerEdge's based on them
433 // map to find _n2eMap of another _SolidData by a shrink shape shared by two _SolidData's
434 map< TGeomID, TNode2Edge* > _s2neMap;
435 // edges of _n2eMap. We keep same data in two containers because
436 // iteration over the map is 5 time longer than over the vector
437 vector< _LayerEdge* > _edges;
439 // key: an id of shape (EDGE or VERTEX) shared by a FACE with
440 // layers and a FACE w/o layers
441 // value: the shape (FACE or EDGE) to shrink mesh on.
442 // _LayerEdge's basing on nodes on key shape are inflated along the value shape
443 map< TGeomID, TopoDS_Shape > _shrinkShape2Shape;
445 // Convex FACEs whose radius of curvature is less than the thickness of layers
446 map< TGeomID, _ConvexFace > _convexFaces;
448 // shapes (EDGEs and VERTEXes) srink from which is forbiden due to collisions with
449 // the adjacent SOLID
450 set< TGeomID > _noShrinkShapes;
452 // <EDGE to smooth on> to <it's curve> -- for analytic smooth
453 map< TGeomID,Handle(Geom_Curve)> _edge2curve;
455 // end indices in _edges of _LayerEdge on each shape, first go shapes to smooth
456 vector< int > _endEdgeOnShape;
457 int _nbShapesToSmooth;
459 double _epsilon; // precision for SegTriaInter()
461 TGeomID _index; // SOLID id, for debug
463 _SolidData(const TopoDS_Shape& s=TopoDS_Shape(),
464 const StdMeshers_ViscousLayers* h=0,
465 const TopoDS_Shape& hs=TopoDS_Shape(),
467 :_solid(s), _hyp(h), _hypShape(hs), _proxyMesh(m) {}
470 Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge& E,
473 Handle(Geom_Surface)& surface,
474 const TopoDS_Face& F,
475 SMESH_MesherHelper& helper);
477 void SortOnEdge( const TopoDS_Edge& E,
480 SMESH_MesherHelper& helper);
482 _ConvexFace* GetConvexFace( const TGeomID faceID )
484 map< TGeomID, _ConvexFace >::iterator id2face = _convexFaces.find( faceID );
485 return id2face == _convexFaces.end() ? 0 : & id2face->second;
487 void GetEdgesOnShape( size_t end, int & iBeg, int & iEnd )
489 iBeg = end > 0 ? _endEdgeOnShape[ end-1 ] : 0;
490 iEnd = _endEdgeOnShape[ end ];
493 bool GetShapeEdges(const TGeomID shapeID, size_t& edgeEnd, int* iBeg=0, int* iEnd=0 ) const;
495 void AddShapesToSmooth( const set< TGeomID >& shapeIDs );
497 //--------------------------------------------------------------------------------
499 * \brief Container of centers of curvature at nodes on an EDGE bounding _ConvexFace
501 struct _CentralCurveOnEdge
504 vector< gp_Pnt > _curvaCenters;
505 vector< _LayerEdge* > _ledges;
506 vector< gp_XYZ > _normals; // new normal for each of _ledges
507 vector< double > _segLength2;
510 TopoDS_Face _adjFace;
511 bool _adjFaceToSmooth;
513 void Append( const gp_Pnt& center, _LayerEdge* ledge )
515 if ( _curvaCenters.size() > 0 )
516 _segLength2.push_back( center.SquareDistance( _curvaCenters.back() ));
517 _curvaCenters.push_back( center );
518 _ledges.push_back( ledge );
519 _normals.push_back( ledge->_normal );
521 bool FindNewNormal( const gp_Pnt& center, gp_XYZ& newNormal );
522 void SetShapes( const TopoDS_Edge& edge,
523 const _ConvexFace& convFace,
524 const _SolidData& data,
525 SMESH_MesherHelper& helper);
527 //--------------------------------------------------------------------------------
529 * \brief Data of node on a shrinked FACE
533 const SMDS_MeshNode* _node;
534 vector<_Simplex> _simplices; // for quality check
536 enum SmoothType { LAPLACIAN, CENTROIDAL, ANGULAR, TFI };
538 bool Smooth(int& badNb,
539 Handle(Geom_Surface)& surface,
540 SMESH_MesherHelper& helper,
541 const double refSign,
545 gp_XY computeAngularPos(vector<gp_XY>& uv,
546 const gp_XY& uvToFix,
547 const double refSign );
549 //--------------------------------------------------------------------------------
551 * \brief Builder of viscous layers
553 class _ViscousBuilder
558 SMESH_ComputeErrorPtr Compute(SMESH_Mesh& mesh,
559 const TopoDS_Shape& shape);
561 // restore event listeners used to clear an inferior dim sub-mesh modified by viscous layers
562 void RestoreListeners();
564 // computes SMESH_ProxyMesh::SubMesh::_n2n;
565 bool MakeN2NMap( _MeshOfSolid* pm );
569 bool findSolidsWithLayers();
570 bool findFacesWithLayers();
571 bool makeLayer(_SolidData& data);
572 bool setEdgeData(_LayerEdge& edge, const set<TGeomID>& subIds,
573 SMESH_MesherHelper& helper, _SolidData& data);
574 gp_XYZ getFaceNormal(const SMDS_MeshNode* n,
575 const TopoDS_Face& face,
576 SMESH_MesherHelper& helper,
578 bool shiftInside=false);
579 gp_XYZ getWeigthedNormal( const SMDS_MeshNode* n,
580 std::pair< TGeomID, gp_XYZ > fId2Normal[],
582 bool findNeiborsOnEdge(const _LayerEdge* edge,
583 const SMDS_MeshNode*& n1,
584 const SMDS_MeshNode*& n2,
586 void getSimplices( const SMDS_MeshNode* node, vector<_Simplex>& simplices,
587 const set<TGeomID>& ingnoreShapes,
588 const _SolidData* dataToCheckOri = 0,
589 const bool toSort = false);
590 void findSimplexTestEdges( _SolidData& data,
591 vector< vector<_LayerEdge*> >& edgesByGeom);
592 bool sortEdges( _SolidData& data,
593 vector< vector<_LayerEdge*> >& edgesByGeom);
594 void limitStepSizeByCurvature( _SolidData& data );
595 void limitStepSize( _SolidData& data,
596 const SMDS_MeshElement* face,
597 const _LayerEdge* maxCosinEdge );
598 void limitStepSize( _SolidData& data, const double minSize);
599 bool inflate(_SolidData& data);
600 bool smoothAndCheck(_SolidData& data, const int nbSteps, double & distToIntersection);
601 bool smoothAnalyticEdge( _SolidData& data,
604 Handle(Geom_Surface)& surface,
605 const TopoDS_Face& F,
606 SMESH_MesherHelper& helper);
607 bool updateNormals( _SolidData& data, SMESH_MesherHelper& helper, int stepNb );
608 bool updateNormalsOfConvexFaces( _SolidData& data,
609 SMESH_MesherHelper& helper,
611 bool refine(_SolidData& data);
613 bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F,
614 SMESH_MesherHelper& helper,
615 const SMESHDS_SubMesh* faceSubMesh );
616 void restoreNoShrink( _LayerEdge& edge ) const;
617 void fixBadFaces(const TopoDS_Face& F,
618 SMESH_MesherHelper& helper,
621 set<const SMDS_MeshNode*> * involvedNodes=NULL);
622 bool addBoundaryElements();
624 bool error( const string& text, int solidID=-1 );
625 SMESHDS_Mesh* getMeshDS() const { return _mesh->GetMeshDS(); }
628 void makeGroupOfLE();
631 SMESH_ComputeErrorPtr _error;
633 vector< _SolidData > _sdVec;
636 //--------------------------------------------------------------------------------
638 * \brief Shrinker of nodes on the EDGE
642 vector<double> _initU;
643 vector<double> _normPar;
644 vector<const SMDS_MeshNode*> _nodes;
645 const _LayerEdge* _edges[2];
648 void AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper );
649 void Compute(bool set3D, SMESH_MesherHelper& helper);
650 void RestoreParams();
651 void SwapSrcTgtNodes(SMESHDS_Mesh* mesh);
653 //--------------------------------------------------------------------------------
655 * \brief Class of temporary mesh face.
656 * We can't use SMDS_FaceOfNodes since it's impossible to set it's ID which is
657 * needed because SMESH_ElementSearcher internaly uses set of elements sorted by ID
659 struct _TmpMeshFace : public SMDS_MeshElement
661 vector<const SMDS_MeshNode* > _nn;
662 _TmpMeshFace( const vector<const SMDS_MeshNode*>& nodes, int id, int faceID=-1):
663 SMDS_MeshElement(id), _nn(nodes) { setShapeId(faceID); }
664 virtual const SMDS_MeshNode* GetNode(const int ind) const { return _nn[ind]; }
665 virtual SMDSAbs_ElementType GetType() const { return SMDSAbs_Face; }
666 virtual vtkIdType GetVtkType() const { return -1; }
667 virtual SMDSAbs_EntityType GetEntityType() const { return SMDSEntity_Last; }
668 virtual SMDSAbs_GeometryType GetGeomType() const
669 { return _nn.size() == 3 ? SMDSGeom_TRIANGLE : SMDSGeom_QUADRANGLE; }
670 virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType) const
671 { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));}
673 //--------------------------------------------------------------------------------
675 * \brief Class of temporary mesh face storing _LayerEdge it's based on
677 struct _TmpMeshFaceOnEdge : public _TmpMeshFace
679 _LayerEdge *_le1, *_le2;
680 _TmpMeshFaceOnEdge( _LayerEdge* le1, _LayerEdge* le2, int ID ):
681 _TmpMeshFace( vector<const SMDS_MeshNode*>(4), ID ), _le1(le1), _le2(le2)
683 _nn[0]=_le1->_nodes[0];
684 _nn[1]=_le1->_nodes.back();
685 _nn[2]=_le2->_nodes.back();
686 _nn[3]=_le2->_nodes[0];
689 //--------------------------------------------------------------------------------
691 * \brief Retriever of node coordinates either directly of from a surface by node UV.
692 * \warning Location of a surface is ignored
694 struct _NodeCoordHelper
696 SMESH_MesherHelper& _helper;
697 const TopoDS_Face& _face;
698 Handle(Geom_Surface) _surface;
699 gp_XYZ (_NodeCoordHelper::* _fun)(const SMDS_MeshNode* n) const;
701 _NodeCoordHelper(const TopoDS_Face& F, SMESH_MesherHelper& helper, bool is2D)
702 : _helper( helper ), _face( F )
707 _surface = BRep_Tool::Surface( _face, loc );
709 if ( _surface.IsNull() )
710 _fun = & _NodeCoordHelper::direct;
712 _fun = & _NodeCoordHelper::byUV;
714 gp_XYZ operator()(const SMDS_MeshNode* n) const { return (this->*_fun)( n ); }
717 gp_XYZ direct(const SMDS_MeshNode* n) const
719 return SMESH_TNodeXYZ( n );
721 gp_XYZ byUV (const SMDS_MeshNode* n) const
723 gp_XY uv = _helper.GetNodeUV( _face, n );
724 return _surface->Value( uv.X(), uv.Y() ).XYZ();
727 } // namespace VISCOUS_3D
731 //================================================================================
732 // StdMeshers_ViscousLayers hypothesis
734 StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen)
735 :SMESH_Hypothesis(hypId, studyId, gen),
736 _isToIgnoreShapes(1), _nbLayers(1), _thickness(1), _stretchFactor(1)
738 _name = StdMeshers_ViscousLayers::GetHypType();
739 _param_algo_dim = -3; // auxiliary hyp used by 3D algos
740 } // --------------------------------------------------------------------------------
741 void StdMeshers_ViscousLayers::SetBndShapes(const std::vector<int>& faceIds, bool toIgnore)
743 if ( faceIds != _shapeIds )
744 _shapeIds = faceIds, NotifySubMeshesHypothesisModification();
745 if ( _isToIgnoreShapes != toIgnore )
746 _isToIgnoreShapes = toIgnore, NotifySubMeshesHypothesisModification();
747 } // --------------------------------------------------------------------------------
748 void StdMeshers_ViscousLayers::SetTotalThickness(double thickness)
750 if ( thickness != _thickness )
751 _thickness = thickness, NotifySubMeshesHypothesisModification();
752 } // --------------------------------------------------------------------------------
753 void StdMeshers_ViscousLayers::SetNumberLayers(int nb)
755 if ( _nbLayers != nb )
756 _nbLayers = nb, NotifySubMeshesHypothesisModification();
757 } // --------------------------------------------------------------------------------
758 void StdMeshers_ViscousLayers::SetStretchFactor(double factor)
760 if ( _stretchFactor != factor )
761 _stretchFactor = factor, NotifySubMeshesHypothesisModification();
762 } // --------------------------------------------------------------------------------
764 StdMeshers_ViscousLayers::Compute(SMESH_Mesh& theMesh,
765 const TopoDS_Shape& theShape,
766 const bool toMakeN2NMap) const
768 using namespace VISCOUS_3D;
769 _ViscousBuilder bulder;
770 SMESH_ComputeErrorPtr err = bulder.Compute( theMesh, theShape );
771 if ( err && !err->IsOK() )
772 return SMESH_ProxyMesh::Ptr();
774 vector<SMESH_ProxyMesh::Ptr> components;
775 TopExp_Explorer exp( theShape, TopAbs_SOLID );
776 for ( ; exp.More(); exp.Next() )
778 if ( _MeshOfSolid* pm =
779 _ViscousListener::GetSolidMesh( &theMesh, exp.Current(), /*toCreate=*/false))
781 if ( toMakeN2NMap && !pm->_n2nMapComputed )
782 if ( !bulder.MakeN2NMap( pm ))
783 return SMESH_ProxyMesh::Ptr();
784 components.push_back( SMESH_ProxyMesh::Ptr( pm ));
785 pm->myIsDeletable = false; // it will de deleted by boost::shared_ptr
787 _ViscousListener::RemoveSolidMesh ( &theMesh, exp.Current() );
789 switch ( components.size() )
793 case 1: return components[0];
795 default: return SMESH_ProxyMesh::Ptr( new SMESH_ProxyMesh( components ));
797 return SMESH_ProxyMesh::Ptr();
798 } // --------------------------------------------------------------------------------
799 std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save)
801 save << " " << _nbLayers
803 << " " << _stretchFactor
804 << " " << _shapeIds.size();
805 for ( size_t i = 0; i < _shapeIds.size(); ++i )
806 save << " " << _shapeIds[i];
807 save << " " << !_isToIgnoreShapes; // negate to keep the behavior in old studies.
809 } // --------------------------------------------------------------------------------
810 std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load)
812 int nbFaces, faceID, shapeToTreat;
813 load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces;
814 while ( _shapeIds.size() < nbFaces && load >> faceID )
815 _shapeIds.push_back( faceID );
816 if ( load >> shapeToTreat )
817 _isToIgnoreShapes = !shapeToTreat;
819 _isToIgnoreShapes = true; // old behavior
821 } // --------------------------------------------------------------------------------
822 bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh* theMesh,
823 const TopoDS_Shape& theShape)
828 // END StdMeshers_ViscousLayers hypothesis
829 //================================================================================
833 gp_XYZ getEdgeDir( const TopoDS_Edge& E, const TopoDS_Vertex& fromV )
837 Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
838 gp_Pnt p = BRep_Tool::Pnt( fromV );
839 double distF = p.SquareDistance( c->Value( f ));
840 double distL = p.SquareDistance( c->Value( l ));
841 c->D1(( distF < distL ? f : l), p, dir );
842 if ( distL < distF ) dir.Reverse();
845 //--------------------------------------------------------------------------------
846 gp_XYZ getEdgeDir( const TopoDS_Edge& E, const SMDS_MeshNode* atNode,
847 SMESH_MesherHelper& helper)
850 double f,l; gp_Pnt p;
851 Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
852 if ( c.IsNull() ) return gp_XYZ( 1e100, 1e100, 1e100 );
853 double u = helper.GetNodeU( E, atNode );
857 //--------------------------------------------------------------------------------
858 gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
859 const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok,
861 //--------------------------------------------------------------------------------
862 gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE,
863 const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok)
866 Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l );
869 TopoDS_Vertex v = helper.IthVertex( 0, fromE );
870 return getFaceDir( F, v, node, helper, ok );
872 gp_XY uv = helper.GetNodeUV( F, node, 0, &ok );
873 Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
874 gp_Pnt p; gp_Vec du, dv, norm;
875 surface->D1( uv.X(),uv.Y(), p, du,dv );
878 double u = helper.GetNodeU( fromE, node, 0, &ok );
880 TopAbs_Orientation o = helper.GetSubShapeOri( F.Oriented(TopAbs_FORWARD), fromE);
881 if ( o == TopAbs_REVERSED )
884 gp_Vec dir = norm ^ du;
886 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX &&
887 helper.IsClosedEdge( fromE ))
889 if ( fabs(u-f) < fabs(u-l)) c->D1( l, p, dv );
890 else c->D1( f, p, dv );
891 if ( o == TopAbs_REVERSED )
893 gp_Vec dir2 = norm ^ dv;
894 dir = dir.Normalized() + dir2.Normalized();
898 //--------------------------------------------------------------------------------
899 gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
900 const SMDS_MeshNode* node, SMESH_MesherHelper& helper,
901 bool& ok, double* cosin)
903 TopoDS_Face faceFrw = F;
904 faceFrw.Orientation( TopAbs_FORWARD );
905 double f,l; TopLoc_Location loc;
906 TopoDS_Edge edges[2]; // sharing a vertex
910 TopExp_Explorer exp( faceFrw, TopAbs_EDGE );
911 for ( ; exp.More() && nbEdges < 2; exp.Next() )
913 const TopoDS_Edge& e = TopoDS::Edge( exp.Current() );
914 if ( SMESH_Algo::isDegenerated( e )) continue;
915 TopExp::Vertices( e, VV[0], VV[1], /*CumOri=*/true );
916 if ( VV[1].IsSame( fromV )) {
917 nbEdges += edges[ 0 ].IsNull();
920 else if ( VV[0].IsSame( fromV )) {
921 nbEdges += edges[ 1 ].IsNull();
926 gp_XYZ dir(0,0,0), edgeDir[2];
929 // get dirs of edges going fromV
931 for ( size_t i = 0; i < nbEdges && ok; ++i )
933 edgeDir[i] = getEdgeDir( edges[i], fromV );
934 double size2 = edgeDir[i].SquareModulus();
935 if (( ok = size2 > numeric_limits<double>::min() ))
936 edgeDir[i] /= sqrt( size2 );
938 if ( !ok ) return dir;
940 // get angle between the 2 edges
942 double angle = helper.GetAngle( edges[0], edges[1], faceFrw, fromV, &faceNormal );
943 if ( Abs( angle ) < 5 * M_PI/180 )
945 dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] );
949 dir = edgeDir[0] + edgeDir[1];
954 double angle = gp_Vec( edgeDir[0] ).Angle( dir );
955 *cosin = Cos( angle );
958 else if ( nbEdges == 1 )
960 dir = getFaceDir( faceFrw, edges[ edges[0].IsNull() ], node, helper, ok );
961 if ( cosin ) *cosin = 1.;
970 //================================================================================
972 * \brief Returns true if a FACE is bound by a concave EDGE
974 //================================================================================
976 bool isConcave( const TopoDS_Face& F, SMESH_MesherHelper& helper )
978 // if ( helper.Count( F, TopAbs_WIRE, /*useMap=*/false) > 1 )
982 TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE );
983 for ( ; eExp.More(); eExp.Next() )
985 const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
986 if ( SMESH_Algo::isDegenerated( E )) continue;
987 // check if 2D curve is concave
988 BRepAdaptor_Curve2d curve( E, F );
989 const int nbIntervals = curve.NbIntervals( GeomAbs_C2 );
990 TColStd_Array1OfReal intervals(1, nbIntervals + 1 );
991 curve.Intervals( intervals, GeomAbs_C2 );
992 bool isConvex = true;
993 for ( int i = 1; i <= nbIntervals && isConvex; ++i )
995 double u1 = intervals( i );
996 double u2 = intervals( i+1 );
997 curve.D2( 0.5*( u1+u2 ), p, drv1, drv2 );
998 double cross = drv2 ^ drv1;
999 if ( E.Orientation() == TopAbs_REVERSED )
1001 isConvex = ( cross > 0.1 ); //-1e-9 );
1005 //cout << "Concave FACE " << helper.GetMeshDS()->ShapeToIndex( F ) << endl;
1009 // check angles at VERTEXes
1011 TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *helper.GetMesh(), 0, error );
1012 for ( size_t iW = 0; iW < wires.size(); ++iW )
1014 const int nbEdges = wires[iW]->NbEdges();
1015 if ( nbEdges < 2 && SMESH_Algo::isDegenerated( wires[iW]->Edge(0)))
1017 for ( int iE1 = 0; iE1 < nbEdges; ++iE1 )
1019 if ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE1 ))) continue;
1020 int iE2 = ( iE1 + 1 ) % nbEdges;
1021 while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 )))
1022 iE2 = ( iE2 + 1 ) % nbEdges;
1023 double angle = helper.GetAngle( wires[iW]->Edge( iE1 ),
1024 wires[iW]->Edge( iE2 ), F,
1025 wires[iW]->FirstVertex( iE2 ));
1026 if ( angle < -5. * M_PI / 180. )
1032 //--------------------------------------------------------------------------------
1033 // DEBUG. Dump intermediate node positions into a python script
1034 // HOWTO use: run python commands written in a console to see
1035 // construction steps of viscous layers
1041 const char* fname = "/tmp/viscous.py";
1042 cout << "execfile('"<<fname<<"')"<<endl;
1043 py = new ofstream(fname);
1044 *py << "import SMESH" << endl
1045 << "from salome.smesh import smeshBuilder" << endl
1046 << "smesh = smeshBuilder.New(salome.myStudy)" << endl
1047 << "meshSO = smesh.GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
1048 << "mesh = smesh.Mesh( meshSO.GetObject() )"<<endl;
1053 *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Viscous Prisms',"
1054 "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA))"<<endl;
1055 *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Neg Volumes',"
1056 "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_Volume3D,'<',0))"<<endl;
1060 ~PyDump() { Finish(); cout << "NB FUNCTIONS: " << theNbFunc << endl; }
1062 #define dumpFunction(f) { _dumpFunction(f, __LINE__);}
1063 #define dumpMove(n) { _dumpMove(n, __LINE__);}
1064 #define dumpCmd(txt) { _dumpCmd(txt, __LINE__);}
1065 void _dumpFunction(const string& fun, int ln)
1066 { if (py) *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl; ++theNbFunc; }
1067 void _dumpMove(const SMDS_MeshNode* n, int ln)
1068 { if (py) *py<< " mesh.MoveNode( "<<n->GetID()<< ", "<< n->X()
1069 << ", "<<n->Y()<<", "<< n->Z()<< ")\t\t # "<< ln <<endl; }
1070 void _dumpCmd(const string& txt, int ln)
1071 { if (py) *py<< " "<<txt<<" # "<< ln <<endl; }
1072 void dumpFunctionEnd()
1073 { if (py) *py<< " return"<< endl; }
1074 void dumpChangeNodes( const SMDS_MeshElement* f )
1075 { if (py) { *py<< " mesh.ChangeElemNodes( " << f->GetID()<<", [";
1076 for ( int i=1; i < f->NbNodes(); ++i ) *py << f->GetNode(i-1)->GetID()<<", ";
1077 *py << f->GetNode( f->NbNodes()-1 )->GetID() << " ])"<< endl; }}
1078 #define debugMsg( txt ) { cout << txt << " (line: " << __LINE__ << ")" << endl; }
1080 struct PyDump { void Finish() {} };
1081 #define dumpFunction(f) f
1083 #define dumpCmd(txt)
1084 #define dumpFunctionEnd()
1085 #define dumpChangeNodes(f)
1086 #define debugMsg( txt ) {}
1090 using namespace VISCOUS_3D;
1092 //================================================================================
1094 * \brief Constructor of _ViscousBuilder
1096 //================================================================================
1098 _ViscousBuilder::_ViscousBuilder()
1100 _error = SMESH_ComputeError::New(COMPERR_OK);
1104 //================================================================================
1106 * \brief Stores error description and returns false
1108 //================================================================================
1110 bool _ViscousBuilder::error(const string& text, int solidId )
1112 const string prefix = string("Viscous layers builder: ");
1113 _error->myName = COMPERR_ALGO_FAILED;
1114 _error->myComment = prefix + text;
1117 SMESH_subMesh* sm = _mesh->GetSubMeshContaining( solidId );
1118 if ( !sm && !_sdVec.empty() )
1119 sm = _mesh->GetSubMeshContaining( solidId = _sdVec[0]._index );
1120 if ( sm && sm->GetSubShape().ShapeType() == TopAbs_SOLID )
1122 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1123 if ( smError && smError->myAlgo )
1124 _error->myAlgo = smError->myAlgo;
1126 sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1128 // set KO to all solids
1129 for ( size_t i = 0; i < _sdVec.size(); ++i )
1131 if ( _sdVec[i]._index == solidId )
1133 sm = _mesh->GetSubMesh( _sdVec[i]._solid );
1134 if ( !sm->IsEmpty() )
1136 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1137 if ( !smError || smError->IsOK() )
1139 smError = SMESH_ComputeError::New( COMPERR_ALGO_FAILED, prefix + "failed");
1140 sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1144 makeGroupOfLE(); // debug
1149 //================================================================================
1151 * \brief At study restoration, restore event listeners used to clear an inferior
1152 * dim sub-mesh modified by viscous layers
1154 //================================================================================
1156 void _ViscousBuilder::RestoreListeners()
1161 //================================================================================
1163 * \brief computes SMESH_ProxyMesh::SubMesh::_n2n
1165 //================================================================================
1167 bool _ViscousBuilder::MakeN2NMap( _MeshOfSolid* pm )
1169 SMESH_subMesh* solidSM = pm->mySubMeshes.front();
1170 TopExp_Explorer fExp( solidSM->GetSubShape(), TopAbs_FACE );
1171 for ( ; fExp.More(); fExp.Next() )
1173 SMESHDS_SubMesh* srcSmDS = pm->GetMeshDS()->MeshElements( fExp.Current() );
1174 const SMESH_ProxyMesh::SubMesh* prxSmDS = pm->GetProxySubMesh( fExp.Current() );
1176 if ( !srcSmDS || !prxSmDS || !srcSmDS->NbElements() || !prxSmDS->NbElements() )
1178 if ( srcSmDS->GetElements()->next() == prxSmDS->GetElements()->next())
1181 if ( srcSmDS->NbElements() != prxSmDS->NbElements() )
1182 return error( "Different nb elements in a source and a proxy sub-mesh", solidSM->GetId());
1184 SMDS_ElemIteratorPtr srcIt = srcSmDS->GetElements();
1185 SMDS_ElemIteratorPtr prxIt = prxSmDS->GetElements();
1186 while( prxIt->more() )
1188 const SMDS_MeshElement* fSrc = srcIt->next();
1189 const SMDS_MeshElement* fPrx = prxIt->next();
1190 if ( fSrc->NbNodes() != fPrx->NbNodes())
1191 return error( "Different elements in a source and a proxy sub-mesh", solidSM->GetId());
1192 for ( int i = 0 ; i < fPrx->NbNodes(); ++i )
1193 pm->setNode2Node( fSrc->GetNode(i), fPrx->GetNode(i), prxSmDS );
1196 pm->_n2nMapComputed = true;
1200 //================================================================================
1202 * \brief Does its job
1204 //================================================================================
1206 SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh,
1207 const TopoDS_Shape& theShape)
1209 // TODO: set priority of solids during Gen::Compute()
1213 // check if proxy mesh already computed
1214 TopExp_Explorer exp( theShape, TopAbs_SOLID );
1216 return error("No SOLID's in theShape"), _error;
1218 if ( _ViscousListener::GetSolidMesh( _mesh, exp.Current(), /*toCreate=*/false))
1219 return SMESH_ComputeErrorPtr(); // everything already computed
1223 // TODO: ignore already computed SOLIDs
1224 if ( !findSolidsWithLayers())
1227 if ( !findFacesWithLayers() )
1230 for ( size_t i = 0; i < _sdVec.size(); ++i )
1232 if ( ! makeLayer(_sdVec[i]) )
1235 if ( _sdVec[i]._edges.size() == 0 )
1238 if ( ! inflate(_sdVec[i]) )
1241 if ( ! refine(_sdVec[i]) )
1247 addBoundaryElements();
1249 makeGroupOfLE(); // debug
1255 //================================================================================
1257 * \brief Finds SOLIDs to compute using viscous layers. Fills _sdVec
1259 //================================================================================
1261 bool _ViscousBuilder::findSolidsWithLayers()
1264 TopTools_IndexedMapOfShape allSolids;
1265 TopExp::MapShapes( _mesh->GetShapeToMesh(), TopAbs_SOLID, allSolids );
1266 _sdVec.reserve( allSolids.Extent());
1268 SMESH_Gen* gen = _mesh->GetGen();
1269 SMESH_HypoFilter filter;
1270 for ( int i = 1; i <= allSolids.Extent(); ++i )
1272 // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
1273 SMESH_Algo* algo = gen->GetAlgo( *_mesh, allSolids(i) );
1274 if ( !algo ) continue;
1275 // TODO: check if algo is hidden
1276 const list <const SMESHDS_Hypothesis *> & allHyps =
1277 algo->GetUsedHypothesis(*_mesh, allSolids(i), /*ignoreAuxiliary=*/false);
1278 list< const SMESHDS_Hypothesis *>::const_iterator hyp = allHyps.begin();
1279 const StdMeshers_ViscousLayers* viscHyp = 0;
1280 for ( ; hyp != allHyps.end() && !viscHyp; ++hyp )
1281 viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp );
1284 TopoDS_Shape hypShape;
1285 filter.Init( filter.Is( viscHyp ));
1286 _mesh->GetHypothesis( allSolids(i), filter, true, &hypShape );
1288 _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
1291 _sdVec.push_back( _SolidData( allSolids(i), viscHyp, hypShape, proxyMesh ));
1292 _sdVec.back()._index = getMeshDS()->ShapeToIndex( allSolids(i));
1295 if ( _sdVec.empty() )
1297 ( SMESH_Comment(StdMeshers_ViscousLayers::GetHypType()) << " hypothesis not found",0);
1302 //================================================================================
1306 //================================================================================
1308 bool _ViscousBuilder::findFacesWithLayers()
1310 SMESH_MesherHelper helper( *_mesh );
1311 TopExp_Explorer exp;
1312 TopTools_IndexedMapOfShape solids;
1314 // collect all faces to ignore defined by hyp
1315 for ( size_t i = 0; i < _sdVec.size(); ++i )
1317 solids.Add( _sdVec[i]._solid );
1319 vector<TGeomID> ids = _sdVec[i]._hyp->GetBndShapes();
1320 if ( _sdVec[i]._hyp->IsToIgnoreShapes() ) // FACEs to ignore are given
1322 for ( size_t ii = 0; ii < ids.size(); ++ii )
1324 const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[ii] );
1325 if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
1326 _sdVec[i]._ignoreFaceIds.insert( ids[ii] );
1329 else // FACEs with layers are given
1331 exp.Init( _sdVec[i]._solid, TopAbs_FACE );
1332 for ( ; exp.More(); exp.Next() )
1334 TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
1335 if ( find( ids.begin(), ids.end(), faceInd ) == ids.end() )
1336 _sdVec[i]._ignoreFaceIds.insert( faceInd );
1340 // ignore internal FACEs if inlets and outlets are specified
1342 TopTools_IndexedDataMapOfShapeListOfShape solidsOfFace;
1343 if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
1344 TopExp::MapShapesAndAncestors( _sdVec[i]._hypShape,
1345 TopAbs_FACE, TopAbs_SOLID, solidsOfFace);
1347 exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
1348 for ( ; exp.More(); exp.Next() )
1350 const TopoDS_Face& face = TopoDS::Face( exp.Current() );
1351 if ( helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) < 2 )
1354 const TGeomID faceInd = getMeshDS()->ShapeToIndex( face );
1355 if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
1357 int nbSolids = solidsOfFace.FindFromKey( face ).Extent();
1359 _sdVec[i]._ignoreFaceIds.insert( faceInd );
1362 if ( helper.IsReversedSubMesh( face ))
1364 _sdVec[i]._reversedFaceIds.insert( faceInd );
1370 // Find faces to shrink mesh on (solution 2 in issue 0020832);
1371 TopTools_IndexedMapOfShape shapes;
1372 for ( size_t i = 0; i < _sdVec.size(); ++i )
1375 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_EDGE, shapes);
1376 for ( int iE = 1; iE <= shapes.Extent(); ++iE )
1378 const TopoDS_Shape& edge = shapes(iE);
1379 // find 2 faces sharing an edge
1381 PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE);
1382 while ( fIt->more())
1384 const TopoDS_Shape* f = fIt->next();
1385 if ( helper.IsSubShape( *f, _sdVec[i]._solid))
1386 FF[ int( !FF[0].IsNull()) ] = *f;
1388 if( FF[1].IsNull() ) continue; // seam edge can be shared by 1 FACE only
1389 // check presence of layers on them
1391 for ( int j = 0; j < 2; ++j )
1392 ignore[j] = _sdVec[i]._ignoreFaceIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
1393 if ( ignore[0] == ignore[1] )
1394 continue; // nothing interesting
1395 TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ];
1396 // check presence of layers on fWOL within an adjacent SOLID
1397 bool collision = false;
1398 PShapeIteratorPtr sIt = helper.GetAncestors( fWOL, *_mesh, TopAbs_SOLID );
1399 while ( const TopoDS_Shape* solid = sIt->next() )
1400 if ( !solid->IsSame( _sdVec[i]._solid ))
1402 int iSolid = solids.FindIndex( *solid );
1403 int iFace = getMeshDS()->ShapeToIndex( fWOL );
1404 if ( iSolid > 0 && !_sdVec[ iSolid-1 ]._ignoreFaceIds.count( iFace ))
1406 //_sdVec[i]._noShrinkShapes.insert( iFace );
1412 if ( !fWOL.IsNull())
1414 TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
1415 _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
1418 // _shrinkShape2Shape will be used to temporary inflate _LayerEdge's based
1419 // on the edge but shrink won't be performed
1420 _sdVec[i]._noShrinkShapes.insert( edgeInd );
1425 // Exclude from _shrinkShape2Shape FACE's that can't be shrinked since
1426 // the algo of the SOLID sharing the FACE does not support it
1427 set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D");
1428 for ( size_t i = 0; i < _sdVec.size(); ++i )
1430 map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
1431 for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
1433 const TopoDS_Shape& fWOL = e2f->second;
1434 const TGeomID edgeID = e2f->first;
1435 bool notShrinkFace = false;
1436 PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
1437 while ( soIt->more() )
1439 const TopoDS_Shape* solid = soIt->next();
1440 if ( _sdVec[i]._solid.IsSame( *solid )) continue;
1441 SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid );
1442 if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue;
1443 notShrinkFace = true;
1445 for ( ; iSolid < _sdVec.size(); ++iSolid )
1447 if ( _sdVec[iSolid]._solid.IsSame( *solid ) ) {
1448 if ( _sdVec[iSolid]._shrinkShape2Shape.count( edgeID ))
1449 notShrinkFace = false;
1453 if ( notShrinkFace )
1455 _sdVec[i]._noShrinkShapes.insert( edgeID );
1457 // add VERTEXes of the edge in _noShrinkShapes
1458 TopoDS_Shape edge = getMeshDS()->IndexToShape( edgeID );
1459 for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() )
1460 _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( vIt.Value() ));
1462 // check if there is a collision with to-shrink-from EDGEs in iSolid
1463 if ( iSolid == _sdVec.size() )
1464 continue; // no VL in the solid
1466 TopExp::MapShapes( fWOL, TopAbs_EDGE, shapes);
1467 for ( int iE = 1; iE <= shapes.Extent(); ++iE )
1469 const TopoDS_Edge& E = TopoDS::Edge( shapes( iE ));
1470 const TGeomID eID = getMeshDS()->ShapeToIndex( E );
1471 if ( eID == edgeID ||
1472 !_sdVec[iSolid]._shrinkShape2Shape.count( eID ) ||
1473 _sdVec[i]._noShrinkShapes.count( eID ))
1475 for ( int is1st = 0; is1st < 2; ++is1st )
1477 TopoDS_Vertex V = helper.IthVertex( is1st, E );
1478 if ( _sdVec[i]._noShrinkShapes.count( getMeshDS()->ShapeToIndex( V ) ))
1480 // _sdVec[i]._noShrinkShapes.insert( eID );
1481 // V = helper.IthVertex( !is1st, E );
1482 // _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( V ));
1483 //iE = 0; // re-start the loop on EDGEs of fWOL
1484 return error("No way to make a conformal mesh with "
1485 "the given set of faces with layers", _sdVec[i]._index);
1491 } // while ( soIt->more() )
1492 } // loop on _sdVec[i]._shrinkShape2Shape
1493 } // loop on _sdVec to fill in _SolidData::_noShrinkShapes
1495 // Find the SHAPE along which to inflate _LayerEdge based on VERTEX
1497 for ( size_t i = 0; i < _sdVec.size(); ++i )
1500 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_VERTEX, shapes);
1501 for ( int iV = 1; iV <= shapes.Extent(); ++iV )
1503 const TopoDS_Shape& vertex = shapes(iV);
1504 // find faces WOL sharing the vertex
1505 vector< TopoDS_Shape > facesWOL;
1506 int totalNbFaces = 0;
1507 PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE);
1508 while ( fIt->more())
1510 const TopoDS_Shape* f = fIt->next();
1511 if ( helper.IsSubShape( *f, _sdVec[i]._solid ) )
1514 const int fID = getMeshDS()->ShapeToIndex( *f );
1515 if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&&
1516 !_sdVec[i]._noShrinkShapes.count( fID )*/)
1517 facesWOL.push_back( *f );
1520 if ( facesWOL.size() == totalNbFaces || facesWOL.empty() )
1521 continue; // no layers at this vertex or no WOL
1522 TGeomID vInd = getMeshDS()->ShapeToIndex( vertex );
1523 switch ( facesWOL.size() )
1527 helper.SetSubShape( facesWOL[0] );
1528 if ( helper.IsRealSeam( vInd )) // inflate along a seam edge?
1530 TopoDS_Shape seamEdge;
1531 PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1532 while ( eIt->more() && seamEdge.IsNull() )
1534 const TopoDS_Shape* e = eIt->next();
1535 if ( helper.IsRealSeam( *e ) )
1538 if ( !seamEdge.IsNull() )
1540 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge ));
1544 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] ));
1549 // find an edge shared by 2 faces
1550 PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1551 while ( eIt->more())
1553 const TopoDS_Shape* e = eIt->next();
1554 if ( helper.IsSubShape( *e, facesWOL[0]) &&
1555 helper.IsSubShape( *e, facesWOL[1]))
1557 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break;
1563 return error("Not yet supported case", _sdVec[i]._index);
1568 // add FACEs of other SOLIDs to _ignoreFaceIds
1569 for ( size_t i = 0; i < _sdVec.size(); ++i )
1572 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_FACE, shapes);
1574 for ( exp.Init( _mesh->GetShapeToMesh(), TopAbs_FACE ); exp.More(); exp.Next() )
1576 if ( !shapes.Contains( exp.Current() ))
1577 _sdVec[i]._ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( exp.Current() ));
1584 //================================================================================
1586 * \brief Create the inner surface of the viscous layer and prepare data for infation
1588 //================================================================================
1590 bool _ViscousBuilder::makeLayer(_SolidData& data)
1592 // get all sub-shapes to make layers on
1593 set<TGeomID> subIds, faceIds;
1594 subIds = data._noShrinkShapes;
1595 TopExp_Explorer exp( data._solid, TopAbs_FACE );
1596 for ( ; exp.More(); exp.Next() )
1598 SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
1599 if ( ! data._ignoreFaceIds.count( fSubM->GetId() ))
1600 faceIds.insert( fSubM->GetId() );
1601 SMESH_subMeshIteratorPtr subIt = fSubM->getDependsOnIterator(/*includeSelf=*/true);
1602 while ( subIt->more() )
1603 subIds.insert( subIt->next()->GetId() );
1606 // make a map to find new nodes on sub-shapes shared with other SOLID
1607 map< TGeomID, TNode2Edge* >::iterator s2ne;
1608 map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
1609 for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
1611 TGeomID shapeInd = s2s->first;
1612 for ( size_t i = 0; i < _sdVec.size(); ++i )
1614 if ( _sdVec[i]._index == data._index ) continue;
1615 map< TGeomID, TopoDS_Shape >::iterator s2s2 = _sdVec[i]._shrinkShape2Shape.find( shapeInd );
1616 if ( s2s2 != _sdVec[i]._shrinkShape2Shape.end() &&
1617 *s2s == *s2s2 && !_sdVec[i]._n2eMap.empty() )
1619 data._s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap ));
1625 // Create temporary faces and _LayerEdge's
1627 dumpFunction(SMESH_Comment("makeLayers_")<<data._index);
1629 data._stepSize = Precision::Infinite();
1630 data._stepSizeNodes[0] = 0;
1632 SMESH_MesherHelper helper( *_mesh );
1633 helper.SetSubShape( data._solid );
1634 helper.SetElementsOnShape( true );
1636 vector< const SMDS_MeshNode*> newNodes; // of a mesh face
1637 TNode2Edge::iterator n2e2;
1639 // collect _LayerEdge's of shapes they are based on
1640 const int nbShapes = getMeshDS()->MaxShapeIndex();
1641 vector< vector<_LayerEdge*> > edgesByGeom( nbShapes+1 );
1643 for ( set<TGeomID>::iterator id = faceIds.begin(); id != faceIds.end(); ++id )
1645 SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( *id );
1646 if ( !smDS ) return error(SMESH_Comment("Not meshed face ") << *id, data._index );
1648 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( *id ));
1649 SMESH_ProxyMesh::SubMesh* proxySub =
1650 data._proxyMesh->getFaceSubM( F, /*create=*/true);
1652 SMDS_ElemIteratorPtr eIt = smDS->GetElements();
1653 while ( eIt->more() )
1655 const SMDS_MeshElement* face = eIt->next();
1656 double faceMaxCosin = -1;
1657 _LayerEdge* maxCosinEdge = 0;
1658 int nbDegenNodes = 0;
1660 newNodes.resize( face->NbCornerNodes() );
1661 for ( size_t i = 0 ; i < newNodes.size(); ++i )
1663 const SMDS_MeshNode* n = face->GetNode( i );
1664 const int shapeID = n->getshapeId();
1665 const bool onDegenShap = helper.IsDegenShape( shapeID );
1666 const bool onDegenEdge = ( onDegenShap && n->GetPosition()->GetDim() == 1 );
1671 // substitute n on a degenerated EDGE with a node on a corresponding VERTEX
1672 const TopoDS_Shape& E = getMeshDS()->IndexToShape( shapeID );
1673 TopoDS_Vertex V = helper.IthVertex( 0, TopoDS::Edge( E ));
1674 if ( const SMDS_MeshNode* vN = SMESH_Algo::VertexNode( V, getMeshDS() )) {
1684 TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first;
1685 if ( !(*n2e).second )
1688 _LayerEdge* edge = new _LayerEdge();
1689 edge->_nodes.push_back( n );
1691 edgesByGeom[ shapeID ].push_back( edge );
1692 const bool noShrink = data._noShrinkShapes.count( shapeID );
1694 SMESH_TNodeXYZ xyz( n );
1696 // set edge data or find already refined _LayerEdge and get data from it
1697 if (( !noShrink ) &&
1698 ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE ) &&
1699 (( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() ) &&
1700 (( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end() ))
1702 _LayerEdge* foundEdge = (*n2e2).second;
1703 gp_XYZ lastPos = edge->Copy( *foundEdge, helper );
1704 foundEdge->_pos.push_back( lastPos );
1705 // location of the last node is modified and we restore it by foundEdge->_pos.back()
1706 const_cast< SMDS_MeshNode* >
1707 ( edge->_nodes.back() )->setXYZ( xyz.X(), xyz.Y(), xyz.Z() );
1713 edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ));
1715 if ( !setEdgeData( *edge, subIds, helper, data ))
1718 dumpMove(edge->_nodes.back());
1720 if ( edge->_cosin > faceMaxCosin )
1722 faceMaxCosin = edge->_cosin;
1723 maxCosinEdge = edge;
1726 newNodes[ i ] = n2e->second->_nodes.back();
1729 data._n2eMap.insert( make_pair( face->GetNode( i ), n2e->second ));
1731 if ( newNodes.size() - nbDegenNodes < 2 )
1734 // create a temporary face
1735 const SMDS_MeshElement* newFace =
1736 new _TmpMeshFace( newNodes, --_tmpFaceID, face->getshapeId() );
1737 proxySub->AddElement( newFace );
1739 // compute inflation step size by min size of element on a convex surface
1740 if ( faceMaxCosin > theMinSmoothCosin )
1741 limitStepSize( data, face, maxCosinEdge );
1743 } // loop on 2D elements on a FACE
1744 } // loop on FACEs of a SOLID
1746 data._epsilon = 1e-7;
1747 if ( data._stepSize < 1. )
1748 data._epsilon *= data._stepSize;
1750 // Put _LayerEdge's into the vector data._edges
1751 if ( !sortEdges( data, edgesByGeom ))
1754 // limit data._stepSize depending on surface curvature and fill data._convexFaces
1755 limitStepSizeByCurvature( data ); // !!! it must be before node substitution in _Simplex
1757 // Set target nodes into _Simplex and _LayerEdge's to _2NearEdges
1758 TNode2Edge::iterator n2e;
1759 const SMDS_MeshNode* nn[2];
1760 for ( size_t i = 0; i < data._edges.size(); ++i )
1762 _LayerEdge* edge = data._edges[i];
1763 if ( edge->IsOnEdge() )
1765 // get neighbor nodes
1766 bool hasData = ( edge->_2neibors->_edges[0] );
1767 if ( hasData ) // _LayerEdge is a copy of another one
1769 nn[0] = edge->_2neibors->srcNode(0);
1770 nn[1] = edge->_2neibors->srcNode(1);
1772 else if ( !findNeiborsOnEdge( edge, nn[0],nn[1], data ))
1776 // set neighbor _LayerEdge's
1777 for ( int j = 0; j < 2; ++j )
1779 if (( n2e = data._n2eMap.find( nn[j] )) == data._n2eMap.end() )
1780 return error("_LayerEdge not found by src node", data._index);
1781 edge->_2neibors->_edges[j] = n2e->second;
1784 edge->SetDataByNeighbors( nn[0], nn[1], helper);
1787 for ( size_t j = 0; j < edge->_simplices.size(); ++j )
1789 _Simplex& s = edge->_simplices[j];
1790 s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
1791 s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
1794 // For an _LayerEdge on a degenerated EDGE, copy some data from
1795 // a corresponding _LayerEdge on a VERTEX
1796 // (issue 52453, pb on a downloaded SampleCase2-Tet-netgen-mephisto.hdf)
1797 if ( helper.IsDegenShape( edge->_nodes[0]->getshapeId() ))
1799 // Generally we should not get here
1800 const TopoDS_Shape& E = getMeshDS()->IndexToShape( edge->_nodes[0]->getshapeId() );
1801 if ( E.ShapeType() != TopAbs_EDGE )
1803 TopoDS_Vertex V = helper.IthVertex( 0, TopoDS::Edge( E ));
1804 const SMDS_MeshNode* vN = SMESH_Algo::VertexNode( V, getMeshDS() );
1805 if (( n2e = data._n2eMap.find( vN )) == data._n2eMap.end() )
1807 const _LayerEdge* vEdge = n2e->second;
1808 edge->_normal = vEdge->_normal;
1809 edge->_lenFactor = vEdge->_lenFactor;
1810 edge->_cosin = vEdge->_cosin;
1818 //================================================================================
1820 * \brief Compute inflation step size by min size of element on a convex surface
1822 //================================================================================
1824 void _ViscousBuilder::limitStepSize( _SolidData& data,
1825 const SMDS_MeshElement* face,
1826 const _LayerEdge* maxCosinEdge )
1829 double minSize = 10 * data._stepSize;
1830 const int nbNodes = face->NbCornerNodes();
1831 for ( int i = 0; i < nbNodes; ++i )
1833 const SMDS_MeshNode* nextN = face->GetNode( SMESH_MesherHelper::WrapIndex( i+1, nbNodes ));
1834 const SMDS_MeshNode* curN = face->GetNode( i );
1835 if ( nextN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ||
1836 curN-> GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
1838 double dist = SMESH_TNodeXYZ( curN ).Distance( nextN );
1839 if ( dist < minSize )
1840 minSize = dist, iN = i;
1843 double newStep = 0.8 * minSize / maxCosinEdge->_lenFactor;
1844 if ( newStep < data._stepSize )
1846 data._stepSize = newStep;
1847 data._stepSizeCoeff = 0.8 / maxCosinEdge->_lenFactor;
1848 data._stepSizeNodes[0] = face->GetNode( iN );
1849 data._stepSizeNodes[1] = face->GetNode( SMESH_MesherHelper::WrapIndex( iN+1, nbNodes ));
1853 //================================================================================
1855 * \brief Compute inflation step size by min size of element on a convex surface
1857 //================================================================================
1859 void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize )
1861 if ( minSize < data._stepSize )
1863 data._stepSize = minSize;
1864 if ( data._stepSizeNodes[0] )
1867 SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
1868 data._stepSizeCoeff = data._stepSize / dist;
1873 //================================================================================
1875 * \brief Limit data._stepSize by evaluating curvature of shapes and fill data._convexFaces
1877 //================================================================================
1879 void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
1881 const int nbTestPnt = 5; // on a FACE sub-shape
1882 const double minCurvature = 0.9 / data._hyp->GetTotalThickness();
1884 BRepLProp_SLProps surfProp( 2, 1e-6 );
1885 SMESH_MesherHelper helper( *_mesh );
1887 data._convexFaces.clear();
1889 TopExp_Explorer face( data._solid, TopAbs_FACE );
1890 for ( ; face.More(); face.Next() )
1892 const TopoDS_Face& F = TopoDS::Face( face.Current() );
1893 SMESH_subMesh * sm = _mesh->GetSubMesh( F );
1894 const TGeomID faceID = sm->GetId();
1895 if ( data._ignoreFaceIds.count( faceID )) continue;
1897 BRepAdaptor_Surface surface( F, false );
1898 surfProp.SetSurface( surface );
1900 bool isTooCurved = false;
1903 _ConvexFace cnvFace;
1904 const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. );
1905 SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
1906 while ( smIt->more() )
1909 const TGeomID subID = sm->GetId();
1910 // find _LayerEdge's of a sub-shape
1912 if ( data.GetShapeEdges( subID, edgesEnd, &iBeg, &iEnd ))
1913 cnvFace._subIdToEdgeEnd.insert( make_pair( subID, edgesEnd ));
1916 // check concavity and curvature and limit data._stepSize
1917 int nbLEdges = iEnd - iBeg;
1918 int iStep = Max( 1, nbLEdges / nbTestPnt );
1919 for ( ; iBeg < iEnd; iBeg += iStep )
1921 gp_XY uv = helper.GetNodeUV( F, data._edges[ iBeg ]->_nodes[0] );
1922 surfProp.SetParameters( uv.X(), uv.Y() );
1923 if ( !surfProp.IsCurvatureDefined() )
1925 if ( surfProp.MaxCurvature() * oriFactor > minCurvature )
1927 limitStepSize( data, 0.9 / surfProp.MaxCurvature() * oriFactor );
1930 if ( surfProp.MinCurvature() * oriFactor > minCurvature )
1932 limitStepSize( data, 0.9 / surfProp.MinCurvature() * oriFactor );
1936 } // loop on sub-shapes of the FACE
1938 if ( !isTooCurved ) continue;
1940 _ConvexFace & convFace =
1941 data._convexFaces.insert( make_pair( faceID, cnvFace )).first->second;
1944 convFace._normalsFixed = false;
1946 // Fill _ConvexFace::_simplexTestEdges. These _LayerEdge's are used to detect
1947 // prism distortion.
1948 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.find( faceID );
1949 if ( id2end != convFace._subIdToEdgeEnd.end() )
1951 // there are _LayerEdge's on the FACE it-self;
1952 // select _LayerEdge's near EDGEs
1953 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
1954 for ( ; iBeg < iEnd; ++iBeg )
1956 _LayerEdge* ledge = data._edges[ iBeg ];
1957 for ( size_t j = 0; j < ledge->_simplices.size(); ++j )
1958 if ( ledge->_simplices[j]._nNext->GetPosition()->GetDim() < 2 )
1960 convFace._simplexTestEdges.push_back( ledge );
1967 // where there are no _LayerEdge's on a _ConvexFace,
1968 // as e.g. on a fillet surface with no internal nodes - issue 22580,
1969 // so that collision of viscous internal faces is not detected by check of
1970 // intersection of _LayerEdge's with the viscous internal faces.
1972 set< const SMDS_MeshNode* > usedNodes;
1974 // look for _LayerEdge's with null _sWOL
1975 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.begin();
1976 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
1978 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
1979 if ( iBeg >= iEnd || !data._edges[ iBeg ]->_sWOL.IsNull() )
1981 for ( ; iBeg < iEnd; ++iBeg )
1983 _LayerEdge* ledge = data._edges[ iBeg ];
1984 const SMDS_MeshNode* srcNode = ledge->_nodes[0];
1985 if ( !usedNodes.insert( srcNode ).second ) continue;
1987 getSimplices( srcNode, ledge->_simplices, data._ignoreFaceIds, &data );
1988 for ( size_t i = 0; i < ledge->_simplices.size(); ++i )
1990 usedNodes.insert( ledge->_simplices[i]._nPrev );
1991 usedNodes.insert( ledge->_simplices[i]._nNext );
1993 convFace._simplexTestEdges.push_back( ledge );
1997 } // loop on FACEs of data._solid
2000 //================================================================================
2002 * \brief Separate shapes (and _LayerEdge's on them) to smooth from the rest ones
2004 //================================================================================
2006 bool _ViscousBuilder::sortEdges( _SolidData& data,
2007 vector< vector<_LayerEdge*> >& edgesByGeom)
2009 // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's
2010 // boundry inclined at a sharp angle to the shape
2012 list< TGeomID > shapesToSmooth;
2014 SMESH_MesherHelper helper( *_mesh );
2017 for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
2019 vector<_LayerEdge*>& eS = edgesByGeom[iS];
2020 if ( eS.empty() ) continue;
2021 const TopoDS_Shape& S = getMeshDS()->IndexToShape( iS );
2022 bool needSmooth = false;
2023 switch ( S.ShapeType() )
2027 if ( SMESH_Algo::isDegenerated( TopoDS::Edge( S )))
2029 //bool isShrinkEdge = !eS[0]->_sWOL.IsNull();
2030 for ( TopoDS_Iterator vIt( S ); vIt.More() && !needSmooth; vIt.Next() )
2032 TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
2033 vector<_LayerEdge*>& eV = edgesByGeom[ iV ];
2034 if ( eV.empty() ) continue;
2035 // double cosin = eV[0]->_cosin;
2037 // ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge));
2040 // gp_Vec dir1, dir2;
2041 // if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE )
2042 // dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ));
2044 // dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ),
2045 // eV[0]->_nodes[0], helper, ok);
2046 // dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
2047 // double angle = dir1.Angle( dir2 );
2048 // cosin = cos( angle );
2050 gp_Vec eDir = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
2051 double angle = eDir.Angle( eV[0]->_normal );
2052 double cosin = Cos( angle );
2053 needSmooth = ( cosin > theMinSmoothCosin );
2059 for ( TopExp_Explorer eExp( S, TopAbs_EDGE ); eExp.More() && !needSmooth; eExp.Next() )
2061 TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
2062 vector<_LayerEdge*>& eE = edgesByGeom[ iE ];
2063 if ( eE.empty() ) continue;
2064 if ( eE[0]->_sWOL.IsNull() )
2066 for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
2067 needSmooth = ( eE[i]->_cosin > theMinSmoothCosin );
2071 const TopoDS_Face& F1 = TopoDS::Face( S );
2072 const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
2073 const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
2074 for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
2076 gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
2077 gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
2078 double angle = dir1.Angle( dir2 );
2079 double cosin = cos( angle );
2080 needSmooth = ( cosin > theMinSmoothCosin );
2092 if ( S.ShapeType() == TopAbs_EDGE ) shapesToSmooth.push_front( iS );
2093 else shapesToSmooth.push_back ( iS );
2096 } // loop on edgesByGeom
2098 data._edges.reserve( data._n2eMap.size() );
2099 data._endEdgeOnShape.clear();
2101 // first we put _LayerEdge's on shapes to smooth
2102 data._nbShapesToSmooth = 0;
2103 list< TGeomID >::iterator gIt = shapesToSmooth.begin();
2104 for ( ; gIt != shapesToSmooth.end(); ++gIt )
2106 vector<_LayerEdge*>& eVec = edgesByGeom[ *gIt ];
2107 if ( eVec.empty() ) continue;
2108 data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
2109 data._endEdgeOnShape.push_back( data._edges.size() );
2110 data._nbShapesToSmooth++;
2114 // then the rest _LayerEdge's
2115 for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
2117 vector<_LayerEdge*>& eVec = edgesByGeom[iS];
2118 if ( eVec.empty() ) continue;
2119 data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
2120 data._endEdgeOnShape.push_back( data._edges.size() );
2127 //================================================================================
2129 * \brief Set data of _LayerEdge needed for smoothing
2130 * \param subIds - ids of sub-shapes of a SOLID to take into account faces from
2132 //================================================================================
2134 bool _ViscousBuilder::setEdgeData(_LayerEdge& edge,
2135 const set<TGeomID>& subIds,
2136 SMESH_MesherHelper& helper,
2139 SMESH_MeshEditor editor(_mesh);
2141 const SMDS_MeshNode* node = edge._nodes[0]; // source node
2142 SMDS_TypeOfPosition posType = node->GetPosition()->GetTypeOfPosition();
2146 edge._curvature = 0;
2148 // --------------------------
2149 // Compute _normal and _cosin
2150 // --------------------------
2153 edge._normal.SetCoord(0,0,0);
2155 int totalNbFaces = 0;
2159 const TGeomID shapeInd = node->getshapeId();
2160 map< TGeomID, TopoDS_Shape >::const_iterator s2s = data._shrinkShape2Shape.find( shapeInd );
2161 const bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
2163 if ( onShrinkShape ) // one of faces the node is on has no layers
2165 TopoDS_Shape vertEdge = getMeshDS()->IndexToShape( s2s->first ); // vertex or edge
2166 if ( s2s->second.ShapeType() == TopAbs_EDGE )
2168 // inflate from VERTEX along EDGE
2169 edge._normal = getEdgeDir( TopoDS::Edge( s2s->second ), TopoDS::Vertex( vertEdge ));
2171 else if ( vertEdge.ShapeType() == TopAbs_VERTEX )
2173 // inflate from VERTEX along FACE
2174 edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Vertex( vertEdge ),
2175 node, helper, normOK, &edge._cosin);
2179 // inflate from EDGE along FACE
2180 edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Edge( vertEdge ),
2181 node, helper, normOK);
2184 else // layers are on all faces of SOLID the node is on
2186 // find indices of geom faces the node lies on
2187 set<TGeomID> faceIds;
2188 if ( posType == SMDS_TOP_FACE )
2190 faceIds.insert( node->getshapeId() );
2194 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2195 while ( fIt->more() )
2196 faceIds.insert( editor.FindShape(fIt->next()));
2199 set<TGeomID>::iterator id = faceIds.begin();
2201 std::pair< TGeomID, gp_XYZ > id2Norm[20];
2202 for ( ; id != faceIds.end(); ++id )
2204 const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id );
2205 if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
2207 F = TopoDS::Face( s );
2208 geomNorm = getFaceNormal( node, F, helper, normOK );
2209 if ( !normOK ) continue;
2211 if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
2213 id2Norm[ totalNbFaces ].first = *id;
2214 id2Norm[ totalNbFaces ].second = geomNorm.XYZ();
2216 edge._normal += geomNorm.XYZ();
2218 if ( totalNbFaces == 0 )
2219 return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index);
2221 if ( normOK && edge._normal.Modulus() < 1e-3 && totalNbFaces > 1 )
2223 // opposite normals, re-get normals at shifted positions (IPAL 52426)
2224 edge._normal.SetCoord( 0,0,0 );
2225 for ( int i = 0; i < totalNbFaces; ++i )
2227 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( id2Norm[i].first ));
2228 geomNorm = getFaceNormal( node, F, helper, normOK, /*shiftInside=*/true );
2229 if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
2232 id2Norm[ i ].second = geomNorm.XYZ();
2233 edge._normal += id2Norm[ i ].second;
2237 if ( totalNbFaces < 3 )
2239 //edge._normal /= totalNbFaces;
2243 edge._normal = getWeigthedNormal( node, id2Norm, totalNbFaces );
2249 edge._cosin = 0; break;
2251 case SMDS_TOP_EDGE: {
2252 TopoDS_Edge E = TopoDS::Edge( helper.GetSubShapeByNode( node, getMeshDS()));
2253 gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK );
2254 double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
2255 edge._cosin = cos( angle );
2256 //cout << "Cosin on EDGE " << edge._cosin << " node " << node->GetID() << endl;
2259 case SMDS_TOP_VERTEX: {
2260 TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS()));
2261 gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK );
2262 double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
2263 edge._cosin = cos( angle );
2264 //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl;
2268 return error(SMESH_Comment("Invalid shape position of node ")<<node, data._index);
2270 } // case _sWOL.IsNull()
2272 double normSize = edge._normal.SquareModulus();
2273 if ( normSize < numeric_limits<double>::min() )
2274 return error(SMESH_Comment("Bad normal at node ")<< node->GetID(), data._index );
2276 edge._normal /= sqrt( normSize );
2278 // TODO: if ( !normOK ) then get normal by mesh faces
2280 // Set the rest data
2281 // --------------------
2282 if ( onShrinkShape )
2284 edge._sWOL = (*s2s).second;
2286 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( edge._nodes.back() );
2287 if ( SMESHDS_SubMesh* sm = getMeshDS()->MeshElements( data._solid ))
2288 sm->RemoveNode( tgtNode , /*isNodeDeleted=*/false );
2290 // set initial position which is parameters on _sWOL in this case
2291 if ( edge._sWOL.ShapeType() == TopAbs_EDGE )
2293 double u = helper.GetNodeU( TopoDS::Edge( edge._sWOL ), node, 0, &normOK );
2294 edge._pos.push_back( gp_XYZ( u, 0, 0 ));
2295 if ( edge._nodes.size() > 1 )
2296 getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( edge._sWOL ), u );
2300 gp_XY uv = helper.GetNodeUV( TopoDS::Face( edge._sWOL ), node, 0, &normOK );
2301 edge._pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
2302 if ( edge._nodes.size() > 1 )
2303 getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( edge._sWOL ), uv.X(), uv.Y() );
2308 edge._pos.push_back( SMESH_TNodeXYZ( node ));
2310 if ( posType == SMDS_TOP_FACE )
2312 getSimplices( node, edge._simplices, data._ignoreFaceIds, &data );
2313 double avgNormProj = 0, avgLen = 0;
2314 for ( size_t i = 0; i < edge._simplices.size(); ++i )
2316 gp_XYZ vec = edge._pos.back() - SMESH_TNodeXYZ( edge._simplices[i]._nPrev );
2317 avgNormProj += edge._normal * vec;
2318 avgLen += vec.Modulus();
2320 avgNormProj /= edge._simplices.size();
2321 avgLen /= edge._simplices.size();
2322 edge._curvature = _Curvature::New( avgNormProj, avgLen );
2326 // Set neighbour nodes for a _LayerEdge based on EDGE
2328 if ( posType == SMDS_TOP_EDGE /*||
2329 ( onShrinkShape && posType == SMDS_TOP_VERTEX && fabs( edge._cosin ) < 1e-10 )*/)
2331 edge._2neibors = new _2NearEdges;
2332 // target node instead of source ones will be set later
2333 // if ( ! findNeiborsOnEdge( &edge,
2334 // edge._2neibors->_nodes[0],
2335 // edge._2neibors->_nodes[1],
2338 // edge.SetDataByNeighbors( edge._2neibors->_nodes[0],
2339 // edge._2neibors->_nodes[1],
2343 edge.SetCosin( edge._cosin ); // to update edge._lenFactor
2348 //================================================================================
2350 * \brief Return normal to a FACE at a node
2351 * \param [in] n - node
2352 * \param [in] face - FACE
2353 * \param [in] helper - helper
2354 * \param [out] isOK - true or false
2355 * \param [in] shiftInside - to find normal at a position shifted inside the face
2356 * \return gp_XYZ - normal
2358 //================================================================================
2360 gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
2361 const TopoDS_Face& face,
2362 SMESH_MesherHelper& helper,
2369 // get a shifted position
2370 gp_Pnt p = SMESH_TNodeXYZ( node );
2371 gp_XYZ shift( 0,0,0 );
2372 TopoDS_Shape S = helper.GetSubShapeByNode( node, helper.GetMeshDS() );
2373 switch ( S.ShapeType() ) {
2376 shift = getFaceDir( face, TopoDS::Vertex( S ), node, helper, isOK );
2381 shift = getFaceDir( face, TopoDS::Edge( S ), node, helper, isOK );
2389 p.Translate( shift * 1e-5 );
2391 TopLoc_Location loc;
2392 GeomAPI_ProjectPointOnSurf& projector = helper.GetProjector( face, loc, 1e-7 );
2394 if ( !loc.IsIdentity() ) p.Transform( loc.Transformation().Inverted() );
2396 projector.Perform( p );
2397 if ( !projector.IsDone() || projector.NbPoints() < 1 )
2402 Quantity_Parameter U,V;
2403 projector.LowerDistanceParameters(U,V);
2408 uv = helper.GetNodeUV( face, node, 0, &isOK );
2414 Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
2415 int pointKind = GeomLib::NormEstim( surface, uv, 1e-5, normal );
2416 enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE };
2418 if ( pointKind == IMPOSSIBLE &&
2419 node->GetPosition()->GetDim() == 2 ) // node inside the FACE
2421 // probably NormEstim() failed due to a too high tolerance
2422 pointKind = GeomLib::NormEstim( surface, uv, 1e-20, normal );
2423 isOK = ( pointKind < IMPOSSIBLE );
2425 if ( pointKind < IMPOSSIBLE )
2427 if ( pointKind != REGULAR &&
2429 node->GetPosition()->GetDim() < 2 ) // FACE boundary
2431 gp_XYZ normShift = getFaceNormal( node, face, helper, isOK, /*shiftInside=*/true );
2432 if ( normShift * normal.XYZ() < 0. )
2438 if ( !isOK ) // hard singularity, to call with shiftInside=true ?
2440 const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face );
2442 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2443 while ( fIt->more() )
2445 const SMDS_MeshElement* f = fIt->next();
2446 if ( f->getshapeId() == faceID )
2448 isOK = SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) normal.XYZ(), /*normalized=*/true );
2451 TopoDS_Face ff = face;
2452 ff.Orientation( TopAbs_FORWARD );
2453 if ( helper.IsReversedSubMesh( ff ))
2460 return normal.XYZ();
2463 //================================================================================
2465 * \brief Return a normal at a node weighted with angles taken by FACEs
2466 * \param [in] n - the node
2467 * \param [in] fId2Normal - FACE ids and normals
2468 * \param [in] nbFaces - nb of FACEs meeting at the node
2469 * \return gp_XYZ - computed normal
2471 //================================================================================
2473 gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode* n,
2474 std::pair< TGeomID, gp_XYZ > fId2Normal[],
2477 gp_XYZ resNorm(0,0,0);
2478 TopoDS_Shape V = SMESH_MesherHelper::GetSubShapeByNode( n, getMeshDS() );
2479 if ( V.ShapeType() != TopAbs_VERTEX )
2481 for ( int i = 0; i < nbFaces; ++i )
2482 resNorm += fId2Normal[i].second / nbFaces ;
2487 for ( int i = 0; i < nbFaces; ++i )
2489 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( fId2Normal[i].first ));
2491 // look for two EDGEs shared by F and other FACEs within fId2Normal
2494 PShapeIteratorPtr eIt = SMESH_MesherHelper::GetAncestors( V, *_mesh, TopAbs_EDGE );
2495 while ( const TopoDS_Shape* E = eIt->next() )
2497 if ( !SMESH_MesherHelper::IsSubShape( *E, F ))
2499 bool isSharedEdge = false;
2500 for ( int j = 0; j < nbFaces && !isSharedEdge; ++j )
2502 if ( i == j ) continue;
2503 const TopoDS_Shape& otherF = getMeshDS()->IndexToShape( fId2Normal[j].first );
2504 isSharedEdge = SMESH_MesherHelper::IsSubShape( *E, otherF );
2506 if ( !isSharedEdge )
2508 ee[ nbE ] = TopoDS::Edge( *E );
2509 ee[ nbE ].Orientation( SMESH_MesherHelper::GetSubShapeOri( F, *E ));
2514 // get an angle between the two EDGEs
2516 if ( nbE < 1 ) continue;
2523 if ( !V.IsSame( SMESH_MesherHelper::IthVertex( 0, ee[ 1 ] )))
2524 std::swap( ee[0], ee[1] );
2526 angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F, TopoDS::Vertex( V ));
2529 // compute a weighted normal
2530 double sumAngle = 0;
2531 for ( int i = 0; i < nbFaces; ++i )
2533 angles[i] = ( angles[i] > 2*M_PI ) ? 0 : M_PI - angles[i];
2534 sumAngle += angles[i];
2536 for ( int i = 0; i < nbFaces; ++i )
2537 resNorm += angles[i] / sumAngle * fId2Normal[i].second;
2542 //================================================================================
2544 * \brief Find 2 neigbor nodes of a node on EDGE
2546 //================================================================================
2548 bool _ViscousBuilder::findNeiborsOnEdge(const _LayerEdge* edge,
2549 const SMDS_MeshNode*& n1,
2550 const SMDS_MeshNode*& n2,
2553 const SMDS_MeshNode* node = edge->_nodes[0];
2554 const int shapeInd = node->getshapeId();
2555 SMESHDS_SubMesh* edgeSM = 0;
2556 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE )
2558 edgeSM = getMeshDS()->MeshElements( shapeInd );
2559 if ( !edgeSM || edgeSM->NbElements() == 0 )
2560 return error(SMESH_Comment("Not meshed EDGE ") << shapeInd, data._index);
2564 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Edge);
2565 while ( eIt->more() && !n2 )
2567 const SMDS_MeshElement* e = eIt->next();
2568 const SMDS_MeshNode* nNeibor = e->GetNode( 0 );
2569 if ( nNeibor == node ) nNeibor = e->GetNode( 1 );
2572 if (!edgeSM->Contains(e)) continue;
2576 TopoDS_Shape s = SMESH_MesherHelper::GetSubShapeByNode(nNeibor, getMeshDS() );
2577 if ( !SMESH_MesherHelper::IsSubShape( s, edge->_sWOL )) continue;
2579 ( iN++ ? n2 : n1 ) = nNeibor;
2582 return error(SMESH_Comment("Wrongly meshed EDGE ") << shapeInd, data._index);
2586 //================================================================================
2588 * \brief Set _curvature and _2neibors->_plnNorm by 2 neigbor nodes residing the same EDGE
2590 //================================================================================
2592 void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
2593 const SMDS_MeshNode* n2,
2594 SMESH_MesherHelper& helper)
2596 if ( _nodes[0]->GetPosition()->GetTypeOfPosition() != SMDS_TOP_EDGE )
2599 gp_XYZ pos = SMESH_TNodeXYZ( _nodes[0] );
2600 gp_XYZ vec1 = pos - SMESH_TNodeXYZ( n1 );
2601 gp_XYZ vec2 = pos - SMESH_TNodeXYZ( n2 );
2605 double sumLen = vec1.Modulus() + vec2.Modulus();
2606 _2neibors->_wgt[0] = 1 - vec1.Modulus() / sumLen;
2607 _2neibors->_wgt[1] = 1 - vec2.Modulus() / sumLen;
2608 double avgNormProj = 0.5 * ( _normal * vec1 + _normal * vec2 );
2609 double avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() );
2610 if ( _curvature ) delete _curvature;
2611 _curvature = _Curvature::New( avgNormProj, avgLen );
2612 // if ( _curvature )
2613 // debugMsg( _nodes[0]->GetID()
2614 // << " CURV r,k: " << _curvature->_r<<","<<_curvature->_k
2615 // << " proj = "<<avgNormProj<< " len = " << avgLen << "| lenDelta(0) = "
2616 // << _curvature->lenDelta(0) );
2620 if ( _sWOL.IsNull() )
2622 TopoDS_Shape S = helper.GetSubShapeByNode( _nodes[0], helper.GetMeshDS() );
2623 TopoDS_Edge E = TopoDS::Edge( S );
2624 // if ( SMESH_Algo::isDegenerated( E ))
2626 gp_XYZ dirE = getEdgeDir( E, _nodes[0], helper );
2627 gp_XYZ plnNorm = dirE ^ _normal;
2628 double proj0 = plnNorm * vec1;
2629 double proj1 = plnNorm * vec2;
2630 if ( fabs( proj0 ) > 1e-10 || fabs( proj1 ) > 1e-10 )
2632 if ( _2neibors->_plnNorm ) delete _2neibors->_plnNorm;
2633 _2neibors->_plnNorm = new gp_XYZ( plnNorm.Normalized() );
2638 //================================================================================
2640 * \brief Copy data from a _LayerEdge of other SOLID and based on the same node;
2641 * this and other _LayerEdge's are inflated along a FACE or an EDGE
2643 //================================================================================
2645 gp_XYZ _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
2647 _nodes = other._nodes;
2648 _normal = other._normal;
2650 _lenFactor = other._lenFactor;
2651 _cosin = other._cosin;
2652 _sWOL = other._sWOL;
2653 _2neibors = other._2neibors;
2654 _curvature = 0; std::swap( _curvature, other._curvature );
2655 _2neibors = 0; std::swap( _2neibors, other._2neibors );
2657 gp_XYZ lastPos( 0,0,0 );
2658 if ( _sWOL.ShapeType() == TopAbs_EDGE )
2660 double u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes[0] );
2661 _pos.push_back( gp_XYZ( u, 0, 0));
2663 u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes.back() );
2668 gp_XY uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes[0]);
2669 _pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
2671 uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes.back() );
2672 lastPos.SetX( uv.X() );
2673 lastPos.SetY( uv.Y() );
2678 //================================================================================
2680 * \brief Set _cosin and _lenFactor
2682 //================================================================================
2684 void _LayerEdge::SetCosin( double cosin )
2687 cosin = Abs( _cosin );
2688 _lenFactor = ( /*0.1 < cosin &&*/ cosin < 1-1e-12 ) ? 1./sqrt(1-cosin*cosin) : 1.0;
2691 //================================================================================
2693 * \brief Fills a vector<_Simplex >
2695 //================================================================================
2697 void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node,
2698 vector<_Simplex>& simplices,
2699 const set<TGeomID>& ingnoreShapes,
2700 const _SolidData* dataToCheckOri,
2704 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2705 while ( fIt->more() )
2707 const SMDS_MeshElement* f = fIt->next();
2708 const TGeomID shapeInd = f->getshapeId();
2709 if ( ingnoreShapes.count( shapeInd )) continue;
2710 const int nbNodes = f->NbCornerNodes();
2711 const int srcInd = f->GetNodeIndex( node );
2712 const SMDS_MeshNode* nPrev = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd-1, nbNodes ));
2713 const SMDS_MeshNode* nNext = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+1, nbNodes ));
2714 const SMDS_MeshNode* nOpp = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+2, nbNodes ));
2715 if ( dataToCheckOri && dataToCheckOri->_reversedFaceIds.count( shapeInd ))
2716 std::swap( nPrev, nNext );
2717 simplices.push_back( _Simplex( nPrev, nNext, nOpp ));
2722 vector<_Simplex> sortedSimplices( simplices.size() );
2723 sortedSimplices[0] = simplices[0];
2725 for ( size_t i = 1; i < simplices.size(); ++i )
2727 for ( size_t j = 1; j < simplices.size(); ++j )
2728 if ( sortedSimplices[i-1]._nNext == simplices[j]._nPrev )
2730 sortedSimplices[i] = simplices[j];
2735 if ( nbFound == simplices.size() - 1 )
2736 simplices.swap( sortedSimplices );
2740 //================================================================================
2742 * \brief DEBUG. Create groups contating temorary data of _LayerEdge's
2744 //================================================================================
2746 void _ViscousBuilder::makeGroupOfLE()
2749 for ( size_t i = 0 ; i < _sdVec.size(); ++i )
2751 if ( _sdVec[i]._edges.empty() ) continue;
2753 dumpFunction( SMESH_Comment("make_LayerEdge_") << i );
2754 for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
2756 _LayerEdge* le = _sdVec[i]._edges[j];
2757 for ( size_t iN = 1; iN < le->_nodes.size(); ++iN )
2758 dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<le->_nodes[iN-1]->GetID()
2759 << ", " << le->_nodes[iN]->GetID() <<"])");
2763 dumpFunction( SMESH_Comment("makeNormals") << i );
2764 for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
2766 _LayerEdge& edge = *_sdVec[i]._edges[j];
2767 SMESH_TNodeXYZ nXYZ( edge._nodes[0] );
2768 nXYZ += edge._normal * _sdVec[i]._stepSize;
2769 dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<edge._nodes[0]->GetID()
2770 << ", mesh.AddNode( " << nXYZ.X()<<","<< nXYZ.Y()<<","<< nXYZ.Z()<<")])");
2774 dumpFunction( SMESH_Comment("makeTmpFaces_") << i );
2775 TopExp_Explorer fExp( _sdVec[i]._solid, TopAbs_FACE );
2776 for ( ; fExp.More(); fExp.Next() )
2778 if (const SMESHDS_SubMesh* sm = _sdVec[i]._proxyMesh->GetProxySubMesh( fExp.Current()))
2780 SMDS_ElemIteratorPtr fIt = sm->GetElements();
2781 while ( fIt->more())
2783 const SMDS_MeshElement* e = fIt->next();
2784 SMESH_Comment cmd("mesh.AddFace([");
2785 for ( int j=0; j < e->NbCornerNodes(); ++j )
2786 cmd << e->GetNode(j)->GetID() << (j+1<e->NbCornerNodes() ? ",": "])");
2796 //================================================================================
2798 * \brief Increase length of _LayerEdge's to reach the required thickness of layers
2800 //================================================================================
2802 bool _ViscousBuilder::inflate(_SolidData& data)
2804 SMESH_MesherHelper helper( *_mesh );
2806 // Limit inflation step size by geometry size found by itersecting
2807 // normals of _LayerEdge's with mesh faces
2808 double geomSize = Precision::Infinite(), intersecDist;
2809 auto_ptr<SMESH_ElementSearcher> searcher
2810 ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
2811 data._proxyMesh->GetFaces( data._solid )) );
2812 for ( size_t i = 0; i < data._edges.size(); ++i )
2814 if ( data._edges[i]->IsOnEdge() ) continue;
2815 data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon );
2816 if ( geomSize > intersecDist && intersecDist > 0 )
2817 geomSize = intersecDist;
2819 if ( data._stepSize > 0.3 * geomSize )
2820 limitStepSize( data, 0.3 * geomSize );
2822 const double tgtThick = data._hyp->GetTotalThickness();
2823 if ( data._stepSize > tgtThick )
2824 limitStepSize( data, tgtThick );
2826 if ( data._stepSize < 1. )
2827 data._epsilon = data._stepSize * 1e-7;
2829 debugMsg( "-- geomSize = " << geomSize << ", stepSize = " << data._stepSize );
2831 double avgThick = 0, curThick = 0, distToIntersection = Precision::Infinite();
2832 int nbSteps = 0, nbRepeats = 0;
2833 while ( 1.01 * avgThick < tgtThick )
2835 // new target length
2836 curThick += data._stepSize;
2837 if ( curThick > tgtThick )
2839 curThick = tgtThick + ( tgtThick-avgThick ) * nbRepeats;
2843 // Elongate _LayerEdge's
2844 dumpFunction(SMESH_Comment("inflate")<<data._index<<"_step"<<nbSteps); // debug
2845 for ( size_t i = 0; i < data._edges.size(); ++i )
2847 data._edges[i]->SetNewLength( curThick, helper );
2851 if ( !updateNormals( data, helper, nbSteps ))
2854 // Improve and check quality
2855 if ( !smoothAndCheck( data, nbSteps, distToIntersection ))
2859 dumpFunction(SMESH_Comment("invalidate")<<data._index<<"_step"<<nbSteps); // debug
2860 for ( size_t i = 0; i < data._edges.size(); ++i )
2862 data._edges[i]->InvalidateStep( nbSteps+1 );
2866 break; // no more inflating possible
2870 // Evaluate achieved thickness
2872 for ( size_t i = 0; i < data._edges.size(); ++i )
2873 avgThick += data._edges[i]->_len;
2874 avgThick /= data._edges.size();
2875 debugMsg( "-- Thickness " << avgThick << " reached" );
2877 if ( distToIntersection < avgThick*1.5 )
2879 debugMsg( "-- Stop inflation since "
2880 << " distToIntersection( "<<distToIntersection<<" ) < avgThick( "
2881 << avgThick << " ) * 1.5" );
2885 limitStepSize( data, 0.25 * distToIntersection );
2886 if ( data._stepSizeNodes[0] )
2887 data._stepSize = data._stepSizeCoeff *
2888 SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
2890 } // while ( 1.01 * avgThick < tgtThick )
2893 return error("failed at the very first inflation step", data._index);
2895 if ( 1.01 * avgThick < tgtThick )
2896 if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( data._index ))
2898 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
2899 if ( !smError || smError->IsOK() )
2901 ( new SMESH_ComputeError (COMPERR_WARNING,
2902 SMESH_Comment("Thickness ") << tgtThick <<
2903 " of viscous layers not reached,"
2904 " average reached thickness is " << avgThick ));
2908 // Restore position of src nodes moved by infaltion on _noShrinkShapes
2909 dumpFunction(SMESH_Comment("restoNoShrink_So")<<data._index); // debug
2911 for ( int iS = 0; iS < data._endEdgeOnShape.size(); ++iS )
2914 iEnd = data._endEdgeOnShape[ iS ];
2915 if ( data._edges[ iBeg ]->_nodes.size() == 1 )
2916 for ( ; iBeg < iEnd; ++iBeg )
2918 restoreNoShrink( *data._edges[ iBeg ] );
2926 //================================================================================
2928 * \brief Improve quality of layer inner surface and check intersection
2930 //================================================================================
2932 bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
2934 double & distToIntersection)
2936 if ( data._nbShapesToSmooth == 0 )
2937 return true; // no shapes needing smoothing
2939 bool moved, improved;
2941 SMESH_MesherHelper helper(*_mesh);
2942 Handle(Geom_Surface) surface;
2946 for ( int iS = 0; iS < data._nbShapesToSmooth; ++iS )
2949 iEnd = data._endEdgeOnShape[ iS ];
2951 if ( !data._edges[ iBeg ]->_sWOL.IsNull() &&
2952 data._edges[ iBeg ]->_sWOL.ShapeType() == TopAbs_FACE )
2954 if ( !F.IsSame( data._edges[ iBeg ]->_sWOL )) {
2955 F = TopoDS::Face( data._edges[ iBeg ]->_sWOL );
2956 helper.SetSubShape( F );
2957 surface = BRep_Tool::Surface( F );
2962 F.Nullify(); surface.Nullify();
2964 TGeomID sInd = data._edges[ iBeg ]->_nodes[0]->getshapeId();
2966 if ( data._edges[ iBeg ]->IsOnEdge() )
2968 dumpFunction(SMESH_Comment("smooth")<<data._index << "_Ed"<<sInd <<"_InfStep"<<nbSteps);
2970 // try a simple solution on an analytic EDGE
2971 if ( !smoothAnalyticEdge( data, iBeg, iEnd, surface, F, helper ))
2977 for ( int i = iBeg; i < iEnd; ++i )
2979 moved |= data._edges[i]->SmoothOnEdge(surface, F, helper);
2981 dumpCmd( SMESH_Comment("# end step ")<<step);
2983 while ( moved && step++ < 5 );
2990 int step = 0, stepLimit = 5, badNb = 0; moved = true;
2991 while (( ++step <= stepLimit && moved ) || improved )
2993 dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
2994 <<"_InfStep"<<nbSteps<<"_"<<step); // debug
2995 int oldBadNb = badNb;
2999 for ( int i = iBeg; i < iEnd; ++i )
3000 moved |= data._edges[i]->Smooth(badNb);
3002 for ( int i = iEnd-1; i >= iBeg; --i )
3003 moved |= data._edges[i]->Smooth(badNb);
3004 improved = ( badNb < oldBadNb );
3006 // issue 22576 -- no bad faces but still there are intersections to fix
3007 if ( improved && badNb == 0 )
3008 stepLimit = step + 3;
3015 for ( int i = iBeg; i < iEnd; ++i )
3017 _LayerEdge* edge = data._edges[i];
3018 SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
3019 for ( size_t j = 0; j < edge->_simplices.size(); ++j )
3020 if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
3022 cout << "Bad simplex ( " << edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
3023 << " "<< edge->_simplices[j]._nPrev->GetID()
3024 << " "<< edge->_simplices[j]._nNext->GetID() << " )" << endl;
3032 } // loop on shapes to smooth
3034 // Check orientation of simplices of _ConvexFace::_simplexTestEdges
3035 map< TGeomID, _ConvexFace >::iterator id2face = data._convexFaces.begin();
3036 for ( ; id2face != data._convexFaces.end(); ++id2face )
3038 _ConvexFace & convFace = (*id2face).second;
3039 if ( !convFace._simplexTestEdges.empty() &&
3040 convFace._simplexTestEdges[0]->_nodes[0]->GetPosition()->GetDim() == 2 )
3041 continue; // _simplexTestEdges are based on FACE -- already checked while smoothing
3043 if ( !convFace.CheckPrisms() )
3047 // Check if the last segments of _LayerEdge intersects 2D elements;
3048 // checked elements are either temporary faces or faces on surfaces w/o the layers
3050 auto_ptr<SMESH_ElementSearcher> searcher
3051 ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
3052 data._proxyMesh->GetFaces( data._solid )) );
3054 distToIntersection = Precision::Infinite();
3056 const SMDS_MeshElement* intFace = 0;
3057 const SMDS_MeshElement* closestFace = 0;
3059 for ( size_t i = 0; i < data._edges.size(); ++i )
3061 if ( !data._edges[i]->_sWOL.IsNull() )
3063 if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace ))
3065 if ( distToIntersection > dist )
3067 // ignore intersection of a _LayerEdge based on a _ConvexFace with a face
3068 // lying on this _ConvexFace
3069 if ( _ConvexFace* convFace = data.GetConvexFace( intFace->getshapeId() ))
3070 if ( convFace->_subIdToEdgeEnd.count ( data._edges[i]->_nodes[0]->getshapeId() ))
3073 distToIntersection = dist;
3075 closestFace = intFace;
3081 SMDS_MeshElement::iterator nIt = closestFace->begin_nodes();
3082 cout << "Shortest distance: _LayerEdge nodes: tgt " << data._edges[iLE]->_nodes.back()->GetID()
3083 << " src " << data._edges[iLE]->_nodes[0]->GetID()<< ", intersection with face ("
3084 << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
3085 << ") distance = " << distToIntersection<< endl;
3092 //================================================================================
3094 * \brief Return a curve of the EDGE to be used for smoothing and arrange
3095 * _LayerEdge's to be in a consequent order
3097 //================================================================================
3099 Handle(Geom_Curve) _SolidData::CurveForSmooth( const TopoDS_Edge& E,
3102 Handle(Geom_Surface)& surface,
3103 const TopoDS_Face& F,
3104 SMESH_MesherHelper& helper)
3106 TGeomID eIndex = helper.GetMeshDS()->ShapeToIndex( E );
3108 map< TGeomID, Handle(Geom_Curve)>::iterator i2curve = _edge2curve.find( eIndex );
3110 if ( i2curve == _edge2curve.end() )
3112 // sort _LayerEdge's by position on the EDGE
3113 SortOnEdge( E, iFrom, iTo, helper );
3115 SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( eIndex );
3117 TopLoc_Location loc; double f,l;
3119 Handle(Geom_Line) line;
3120 Handle(Geom_Circle) circle;
3121 bool isLine, isCirc;
3122 if ( F.IsNull() ) // 3D case
3124 // check if the EDGE is a line
3125 Handle(Geom_Curve) curve = BRep_Tool::Curve( E, loc, f, l);
3126 if ( curve->IsKind( STANDARD_TYPE( Geom_TrimmedCurve )))
3127 curve = Handle(Geom_TrimmedCurve)::DownCast( curve )->BasisCurve();
3129 line = Handle(Geom_Line)::DownCast( curve );
3130 circle = Handle(Geom_Circle)::DownCast( curve );
3131 isLine = (!line.IsNull());
3132 isCirc = (!circle.IsNull());
3134 if ( !isLine && !isCirc ) // Check if the EDGE is close to a line
3137 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
3138 while ( nIt->more() )
3139 bndBox.Add( SMESH_TNodeXYZ( nIt->next() ));
3140 gp_XYZ size = bndBox.CornerMax() - bndBox.CornerMin();
3142 SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->tgtNode(0) );
3143 SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->tgtNode(1) );
3144 const double lineTol = 1e-2 * ( p0 - p1 ).Modulus();
3145 for ( int i = 0; i < 3 && !isLine; ++i )
3146 isLine = ( size.Coord( i+1 ) <= lineTol );
3148 if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle
3155 // check if the EDGE is a line
3156 Handle(Geom2d_Curve) curve = BRep_Tool::CurveOnSurface( E, F, f, l);
3157 if ( curve->IsKind( STANDARD_TYPE( Geom2d_TrimmedCurve )))
3158 curve = Handle(Geom2d_TrimmedCurve)::DownCast( curve )->BasisCurve();
3160 Handle(Geom2d_Line) line2d = Handle(Geom2d_Line)::DownCast( curve );
3161 Handle(Geom2d_Circle) circle2d = Handle(Geom2d_Circle)::DownCast( curve );
3162 isLine = (!line2d.IsNull());
3163 isCirc = (!circle2d.IsNull());
3165 if ( !isLine && !isCirc) // Check if the EDGE is close to a line
3168 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
3169 while ( nIt->more() )
3170 bndBox.Add( helper.GetNodeUV( F, nIt->next() ));
3171 gp_XY size = bndBox.CornerMax() - bndBox.CornerMin();
3173 const double lineTol = 1e-2 * sqrt( bndBox.SquareExtent() );
3174 for ( int i = 0; i < 2 && !isLine; ++i )
3175 isLine = ( size.Coord( i+1 ) <= lineTol );
3177 if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle
3183 line = new Geom_Line( gp::OX() ); // only type does matter
3187 gp_Pnt2d p = circle2d->Location();
3188 gp_Ax2 ax( gp_Pnt( p.X(), p.Y(), 0), gp::DX());
3189 circle = new Geom_Circle( ax, 1.); // only center position does matter
3193 Handle(Geom_Curve)& res = _edge2curve[ eIndex ];
3201 return i2curve->second;
3204 //================================================================================
3206 * \brief Sort _LayerEdge's by a parameter on a given EDGE
3208 //================================================================================
3210 void _SolidData::SortOnEdge( const TopoDS_Edge& E,
3213 SMESH_MesherHelper& helper)
3215 map< double, _LayerEdge* > u2edge;
3216 for ( int i = iFrom; i < iTo; ++i )
3217 u2edge.insert( make_pair( helper.GetNodeU( E, _edges[i]->_nodes[0] ), _edges[i] ));
3219 ASSERT( u2edge.size() == iTo - iFrom );
3220 map< double, _LayerEdge* >::iterator u2e = u2edge.begin();
3221 for ( int i = iFrom; i < iTo; ++i, ++u2e )
3222 _edges[i] = u2e->second;
3224 // set _2neibors according to the new order
3225 for ( int i = iFrom; i < iTo-1; ++i )
3226 if ( _edges[i]->_2neibors->tgtNode(1) != _edges[i+1]->_nodes.back() )
3227 _edges[i]->_2neibors->reverse();
3228 if ( u2edge.size() > 1 &&
3229 _edges[iTo-1]->_2neibors->tgtNode(0) != _edges[iTo-2]->_nodes.back() )
3230 _edges[iTo-1]->_2neibors->reverse();
3233 //================================================================================
3235 * \brief Return index corresponding to the shape in _endEdgeOnShape
3237 //================================================================================
3239 bool _SolidData::GetShapeEdges(const TGeomID shapeID,
3244 int beg = 0, end = 0;
3245 for ( edgesEnd = 0; edgesEnd < _endEdgeOnShape.size(); ++edgesEnd )
3247 end = _endEdgeOnShape[ edgesEnd ];
3248 TGeomID sID = _edges[ beg ]->_nodes[0]->getshapeId();
3249 if ( sID == shapeID )
3251 if ( iBeg ) *iBeg = beg;
3252 if ( iEnd ) *iEnd = end;
3260 //================================================================================
3262 * \brief Add faces for smoothing
3264 //================================================================================
3266 void _SolidData::AddShapesToSmooth( const set< TGeomID >& faceIDs )
3268 // convert faceIDs to indices in _endEdgeOnShape
3269 set< size_t > iEnds;
3271 set< TGeomID >::const_iterator fId = faceIDs.begin();
3272 for ( ; fId != faceIDs.end(); ++fId )
3273 if ( GetShapeEdges( *fId, end ) && end >= _nbShapesToSmooth )
3274 iEnds.insert( end );
3276 set< size_t >::iterator endsIt = iEnds.begin();
3278 // "add" by move of _nbShapesToSmooth
3279 int nbFacesToAdd = iEnds.size();
3280 while ( endsIt != iEnds.end() && *endsIt == _nbShapesToSmooth )
3283 ++_nbShapesToSmooth;
3286 if ( endsIt == iEnds.end() )
3289 // Move _LayerEdge's on FACEs just after _nbShapesToSmooth
3291 vector< _LayerEdge* > nonSmoothLE, smoothLE;
3292 size_t lastSmooth = *iEnds.rbegin();
3294 for ( size_t i = _nbShapesToSmooth; i <= lastSmooth; ++i )
3296 vector< _LayerEdge* > & edgesVec = iEnds.count(i) ? smoothLE : nonSmoothLE;
3297 iBeg = i ? _endEdgeOnShape[ i-1 ] : 0;
3298 iEnd = _endEdgeOnShape[ i ];
3299 edgesVec.insert( edgesVec.end(), _edges.begin() + iBeg, _edges.begin() + iEnd );
3302 iBeg = _nbShapesToSmooth ? _endEdgeOnShape[ _nbShapesToSmooth-1 ] : 0;
3303 std::copy( smoothLE.begin(), smoothLE.end(), &_edges[ iBeg ] );
3304 std::copy( nonSmoothLE.begin(), nonSmoothLE.end(), &_edges[ iBeg + smoothLE.size()]);
3306 // update _endEdgeOnShape
3307 for ( size_t i = _nbShapesToSmooth; i < _endEdgeOnShape.size(); ++i )
3309 TGeomID curShape = _edges[ iBeg ]->_nodes[0]->getshapeId();
3310 while ( ++iBeg < _edges.size() &&
3311 curShape == _edges[ iBeg ]->_nodes[0]->getshapeId() );
3313 _endEdgeOnShape[ i ] = iBeg;
3316 _nbShapesToSmooth += nbFacesToAdd;
3319 //================================================================================
3321 * \brief smooth _LayerEdge's on a staight EDGE or circular EDGE
3323 //================================================================================
3325 bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data,
3328 Handle(Geom_Surface)& surface,
3329 const TopoDS_Face& F,
3330 SMESH_MesherHelper& helper)
3332 TopoDS_Shape S = helper.GetSubShapeByNode( data._edges[ iFrom ]->_nodes[0],
3333 helper.GetMeshDS());
3334 TopoDS_Edge E = TopoDS::Edge( S );
3336 Handle(Geom_Curve) curve = data.CurveForSmooth( E, iFrom, iTo, surface, F, helper );
3337 if ( curve.IsNull() ) return false;
3339 // compute a relative length of segments
3340 vector< double > len( iTo-iFrom+1 );
3342 double curLen, prevLen = len[0] = 1.0;
3343 for ( int i = iFrom; i < iTo; ++i )
3345 curLen = prevLen * data._edges[i]->_2neibors->_wgt[0] / data._edges[i]->_2neibors->_wgt[1];
3346 len[i-iFrom+1] = len[i-iFrom] + curLen;
3351 if ( curve->IsKind( STANDARD_TYPE( Geom_Line )))
3353 if ( F.IsNull() ) // 3D
3355 SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->tgtNode(0));
3356 SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->tgtNode(1));
3357 for ( int i = iFrom; i < iTo; ++i )
3359 double r = len[i-iFrom] / len.back();
3360 gp_XYZ newPos = p0 * ( 1. - r ) + p1 * r;
3361 data._edges[i]->_pos.back() = newPos;
3362 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3363 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3364 dumpMove( tgtNode );
3369 // gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->tgtNode(0));
3370 // gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->tgtNode(1));
3371 gp_XY uv0 = data._edges[iFrom]->_2neibors->_edges[0]->LastUV( F );
3372 gp_XY uv1 = data._edges[iTo-1]->_2neibors->_edges[1]->LastUV( F );
3373 if ( data._edges[iFrom]->_2neibors->tgtNode(0) ==
3374 data._edges[iTo-1]->_2neibors->tgtNode(1) ) // closed edge
3376 int iPeriodic = helper.GetPeriodicIndex();
3377 if ( iPeriodic == 1 || iPeriodic == 2 )
3379 uv1.SetCoord( iPeriodic, helper.GetOtherParam( uv1.Coord( iPeriodic )));
3380 if ( uv0.Coord( iPeriodic ) > uv1.Coord( iPeriodic ))
3381 std::swap( uv0, uv1 );
3384 const gp_XY rangeUV = uv1 - uv0;
3385 for ( int i = iFrom; i < iTo; ++i )
3387 double r = len[i-iFrom] / len.back();
3388 gp_XY newUV = uv0 + r * rangeUV;
3389 data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
3391 gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
3392 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3393 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3394 dumpMove( tgtNode );
3396 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
3397 pos->SetUParameter( newUV.X() );
3398 pos->SetVParameter( newUV.Y() );
3404 if ( curve->IsKind( STANDARD_TYPE( Geom_Circle )))
3406 Handle(Geom_Circle) circle = Handle(Geom_Circle)::DownCast( curve );
3407 gp_Pnt center3D = circle->Location();
3409 if ( F.IsNull() ) // 3D
3411 if ( data._edges[iFrom]->_2neibors->tgtNode(0) ==
3412 data._edges[iTo-1]->_2neibors->tgtNode(1) )
3413 return true; // closed EDGE - nothing to do
3415 return false; // TODO ???
3419 const gp_XY center( center3D.X(), center3D.Y() );
3421 gp_XY uv0 = data._edges[iFrom]->_2neibors->_edges[0]->LastUV( F );
3422 gp_XY uvM = data._edges[iFrom]->LastUV( F );
3423 gp_XY uv1 = data._edges[iTo-1]->_2neibors->_edges[1]->LastUV( F );
3424 // gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->tgtNode(0));
3425 // gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back());
3426 // gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->tgtNode(1));
3427 gp_Vec2d vec0( center, uv0 );
3428 gp_Vec2d vecM( center, uvM );
3429 gp_Vec2d vec1( center, uv1 );
3430 double uLast = vec0.Angle( vec1 ); // -PI - +PI
3431 double uMidl = vec0.Angle( vecM );
3432 if ( uLast * uMidl <= 0. )
3433 uLast += ( uMidl > 0 ? +2. : -2. ) * M_PI;
3434 const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() );
3436 gp_Ax2d axis( center, vec0 );
3437 gp_Circ2d circ( axis, radius );
3438 for ( int i = iFrom; i < iTo; ++i )
3440 double newU = uLast * len[i-iFrom] / len.back();
3441 gp_Pnt2d newUV = ElCLib::Value( newU, circ );
3442 data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
3444 gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
3445 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3446 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3447 dumpMove( tgtNode );
3449 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
3450 pos->SetUParameter( newUV.X() );
3451 pos->SetVParameter( newUV.Y() );
3460 //================================================================================
3462 * \brief Modify normals of _LayerEdge's on EDGE's to avoid intersection with
3463 * _LayerEdge's on neighbor EDGE's
3465 //================================================================================
3467 bool _ViscousBuilder::updateNormals( _SolidData& data,
3468 SMESH_MesherHelper& helper,
3472 return updateNormalsOfConvexFaces( data, helper, stepNb );
3474 // make temporary quadrangles got by extrusion of
3475 // mesh edges along _LayerEdge._normal's
3477 vector< const SMDS_MeshElement* > tmpFaces;
3479 set< SMESH_TLink > extrudedLinks; // contains target nodes
3480 vector< const SMDS_MeshNode*> nodes(4); // of a tmp mesh face
3482 dumpFunction(SMESH_Comment("makeTmpFacesOnEdges")<<data._index);
3483 for ( size_t i = 0; i < data._edges.size(); ++i )
3485 _LayerEdge* edge = data._edges[i];
3486 if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
3487 const SMDS_MeshNode* tgt1 = edge->_nodes.back();
3488 for ( int j = 0; j < 2; ++j ) // loop on _2NearEdges
3490 const SMDS_MeshNode* tgt2 = edge->_2neibors->tgtNode(j);
3491 pair< set< SMESH_TLink >::iterator, bool > link_isnew =
3492 extrudedLinks.insert( SMESH_TLink( tgt1, tgt2 ));
3493 if ( !link_isnew.second )
3495 extrudedLinks.erase( link_isnew.first );
3496 continue; // already extruded and will no more encounter
3498 // a _LayerEdge containg tgt2
3499 _LayerEdge* neiborEdge = edge->_2neibors->_edges[j];
3501 _TmpMeshFaceOnEdge* f = new _TmpMeshFaceOnEdge( edge, neiborEdge, --_tmpFaceID );
3502 tmpFaces.push_back( f );
3504 dumpCmd(SMESH_Comment("mesh.AddFace([ ")
3505 <<f->_nn[0]->GetID()<<", "<<f->_nn[1]->GetID()<<", "
3506 <<f->_nn[2]->GetID()<<", "<<f->_nn[3]->GetID()<<" ])");
3511 // Check if _LayerEdge's based on EDGE's intersects tmpFaces.
3512 // Perform two loops on _LayerEdge on EDGE's:
3513 // 1) to find and fix intersection
3514 // 2) to check that no new intersection appears as result of 1)
3516 SMDS_ElemIteratorPtr fIt( new SMDS_ElementVectorIterator( tmpFaces.begin(),
3518 auto_ptr<SMESH_ElementSearcher> searcher
3519 ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), fIt ));
3521 // 1) Find intersections
3523 const SMDS_MeshElement* face;
3524 typedef map< _LayerEdge*, set< _LayerEdge*, _LayerEdgeCmp >, _LayerEdgeCmp > TLEdge2LEdgeSet;
3525 TLEdge2LEdgeSet edge2CloseEdge;
3527 const double eps = data._epsilon * data._epsilon;
3528 for ( size_t i = 0; i < data._edges.size(); ++i )
3530 _LayerEdge* edge = data._edges[i];
3531 if (( !edge->IsOnEdge() ) &&
3532 ( edge->_sWOL.IsNull() || edge->_sWOL.ShapeType() != TopAbs_FACE ))
3534 if ( edge->FindIntersection( *searcher, dist, eps, &face ))
3536 const _TmpMeshFaceOnEdge* f = (const _TmpMeshFaceOnEdge*) face;
3537 set< _LayerEdge*, _LayerEdgeCmp > & ee = edge2CloseEdge[ edge ];
3538 ee.insert( f->_le1 );
3539 ee.insert( f->_le2 );
3540 if ( f->_le1->IsOnEdge() && f->_le1->_sWOL.IsNull() )
3541 edge2CloseEdge[ f->_le1 ].insert( edge );
3542 if ( f->_le2->IsOnEdge() && f->_le2->_sWOL.IsNull() )
3543 edge2CloseEdge[ f->_le2 ].insert( edge );
3547 // Set _LayerEdge._normal
3549 if ( !edge2CloseEdge.empty() )
3551 dumpFunction(SMESH_Comment("updateNormals")<<data._index);
3553 set< TGeomID > shapesToSmooth;
3555 // vector to store new _normal and _cosin for each edge in edge2CloseEdge
3556 vector< pair< _LayerEdge*, _LayerEdge > > edge2newEdge( edge2CloseEdge.size() );
3558 TLEdge2LEdgeSet::iterator e2ee = edge2CloseEdge.begin();
3559 for ( size_t iE = 0; e2ee != edge2CloseEdge.end(); ++e2ee, ++iE )
3561 _LayerEdge* edge1 = e2ee->first;
3562 _LayerEdge* edge2 = 0;
3563 set< _LayerEdge*, _LayerEdgeCmp >& ee = e2ee->second;
3565 edge2newEdge[ iE ].first = NULL;
3567 // find EDGEs the edges reside
3568 // TopoDS_Edge E1, E2;
3569 // TopoDS_Shape S = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
3570 // if ( S.ShapeType() != TopAbs_EDGE )
3571 // continue; // TODO: find EDGE by VERTEX
3572 // E1 = TopoDS::Edge( S );
3573 set< _LayerEdge*, _LayerEdgeCmp >::iterator eIt = ee.begin();
3574 for ( ; !edge2 && eIt != ee.end(); ++eIt )
3576 if ( edge1->_sWOL == (*eIt)->_sWOL )
3579 if ( !edge2 ) continue;
3581 edge2newEdge[ iE ].first = edge1;
3582 _LayerEdge& newEdge = edge2newEdge[ iE ].second;
3583 // while ( E2.IsNull() && eIt != ee.end())
3585 // _LayerEdge* e2 = *eIt++;
3586 // TopoDS_Shape S = helper.GetSubShapeByNode( e2->_nodes[0], getMeshDS() );
3587 // if ( S.ShapeType() == TopAbs_EDGE )
3588 // E2 = TopoDS::Edge( S ), edge2 = e2;
3590 // if ( E2.IsNull() ) continue; // TODO: find EDGE by VERTEX
3592 // find 3 FACEs sharing 2 EDGEs
3594 // TopoDS_Face FF1[2], FF2[2];
3595 // PShapeIteratorPtr fIt = helper.GetAncestors(E1, *_mesh, TopAbs_FACE);
3596 // while ( fIt->more() && FF1[1].IsNull() )
3598 // const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
3599 // if ( helper.IsSubShape( *F, data._solid))
3600 // FF1[ FF1[0].IsNull() ? 0 : 1 ] = *F;
3602 // fIt = helper.GetAncestors(E2, *_mesh, TopAbs_FACE);
3603 // while ( fIt->more() && FF2[1].IsNull())
3605 // const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
3606 // if ( helper.IsSubShape( *F, data._solid))
3607 // FF2[ FF2[0].IsNull() ? 0 : 1 ] = *F;
3609 // // exclude a FACE common to E1 and E2 (put it to FFn[1] )
3610 // if ( FF1[0].IsSame( FF2[0]) || FF1[0].IsSame( FF2[1]))
3611 // std::swap( FF1[0], FF1[1] );
3612 // if ( FF2[0].IsSame( FF1[0]) )
3613 // std::swap( FF2[0], FF2[1] );
3614 // if ( FF1[0].IsNull() || FF2[0].IsNull() )
3617 // get a new normal for edge1
3619 gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal;
3620 // if ( edge1->_cosin < 0 )
3621 // dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized();
3622 // if ( edge2->_cosin < 0 )
3623 // dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized();
3625 double cos1 = Abs( edge1->_cosin ), cos2 = Abs( edge2->_cosin );
3626 double wgt1 = ( cos1 + 0.001 ) / ( cos1 + cos2 + 0.002 );
3627 double wgt2 = ( cos2 + 0.001 ) / ( cos1 + cos2 + 0.002 );
3628 newEdge._normal = ( wgt1 * dir1 + wgt2 * dir2 ).XYZ();
3629 newEdge._normal.Normalize();
3631 // cout << edge1->_nodes[0]->GetID() << " "
3632 // << edge2->_nodes[0]->GetID() << " NORM: "
3633 // << newEdge._normal.X() << ", " << newEdge._normal.Y() << ", " << newEdge._normal.Z() << endl;
3636 if ( cos1 < theMinSmoothCosin )
3638 newEdge._cosin = edge2->_cosin;
3640 else if ( cos2 > theMinSmoothCosin ) // both cos1 and cos2 > theMinSmoothCosin
3642 // gp_Vec dirInFace;
3643 // if ( edge1->_cosin < 0 )
3644 // dirInFace = dir1;
3646 // dirInFace = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
3647 // double angle = dirInFace.Angle( edge1->_normal ); // [0,PI]
3648 // edge1->SetCosin( Cos( angle ));
3649 //newEdge._cosin = 0; // ???????????
3650 newEdge._cosin = ( wgt1 * cos1 + wgt2 * cos2 ) * edge1->_cosin / cos1;
3654 newEdge._cosin = edge1->_cosin;
3657 // find shapes that need smoothing due to change of _normal
3658 if ( edge1->_cosin < theMinSmoothCosin &&
3659 newEdge._cosin > theMinSmoothCosin )
3661 if ( edge1->_sWOL.IsNull() )
3663 SMDS_ElemIteratorPtr fIt = edge1->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
3664 while ( fIt->more() )
3665 shapesToSmooth.insert( fIt->next()->getshapeId() );
3666 //limitStepSize( data, fIt->next(), edge1->_cosin ); // too late
3668 else // edge1 inflates along a FACE
3670 TopoDS_Shape V = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
3671 PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE );
3672 while ( const TopoDS_Shape* E = eIt->next() )
3674 if ( !helper.IsSubShape( *E, /*FACE=*/edge1->_sWOL ))
3676 gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V ));
3677 double angle = edgeDir.Angle( newEdge._normal ); // [0,PI]
3678 if ( angle < M_PI / 2 )
3679 shapesToSmooth.insert( getMeshDS()->ShapeToIndex( *E ));
3685 data.AddShapesToSmooth( shapesToSmooth );
3687 // Update data of edges depending on a new _normal
3689 for ( size_t iE = 0; iE < edge2newEdge.size(); ++iE )
3691 _LayerEdge* edge1 = edge2newEdge[ iE ].first;
3692 _LayerEdge& newEdge = edge2newEdge[ iE ].second;
3693 if ( !edge1 ) continue;
3695 edge1->_normal = newEdge._normal;
3696 edge1->SetCosin( newEdge._cosin );
3697 edge1->InvalidateStep( 1 );
3699 edge1->SetNewLength( data._stepSize, helper );
3700 if ( edge1->IsOnEdge() )
3702 const SMDS_MeshNode * n1 = edge1->_2neibors->srcNode(0);
3703 const SMDS_MeshNode * n2 = edge1->_2neibors->srcNode(1);
3704 edge1->SetDataByNeighbors( n1, n2, helper );
3707 // Update normals and other dependent data of not intersecting _LayerEdge's
3708 // neighboring the intersecting ones
3710 if ( !edge1->_2neibors )
3712 for ( int j = 0; j < 2; ++j ) // loop on 2 neighbors
3714 _LayerEdge* neighbor = edge1->_2neibors->_edges[j];
3715 if ( edge2CloseEdge.count ( neighbor ))
3716 continue; // j-th neighbor is also intersected
3717 _LayerEdge* prevEdge = edge1;
3718 const int nbSteps = 10;
3719 for ( int step = nbSteps; step; --step ) // step from edge1 in j-th direction
3721 if ( !neighbor->_2neibors )
3722 break; // neighbor is on VERTEX
3724 _LayerEdge* nextEdge = neighbor->_2neibors->_edges[iNext];
3725 if ( nextEdge == prevEdge )
3726 nextEdge = neighbor->_2neibors->_edges[ ++iNext ];
3727 double r = double(step-1)/nbSteps;
3728 if ( !nextEdge->_2neibors )
3731 gp_XYZ newNorm = prevEdge->_normal * r + nextEdge->_normal * (1-r);
3732 newNorm.Normalize();
3734 neighbor->_normal = newNorm;
3735 neighbor->SetCosin( prevEdge->_cosin * r + nextEdge->_cosin * (1-r) );
3736 neighbor->SetDataByNeighbors( prevEdge->_nodes[0], nextEdge->_nodes[0], helper );
3738 neighbor->InvalidateStep( 1 );
3740 neighbor->SetNewLength( data._stepSize, helper );
3742 // goto the next neighbor
3743 prevEdge = neighbor;
3744 neighbor = nextEdge;
3750 // 2) Check absence of intersections
3753 for ( size_t i = 0 ; i < tmpFaces.size(); ++i )
3759 //================================================================================
3761 * \brief Modify normals of _LayerEdge's on _ConvexFace's
3763 //================================================================================
3765 bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData& data,
3766 SMESH_MesherHelper& helper,
3769 SMESHDS_Mesh* meshDS = helper.GetMeshDS();
3772 map< TGeomID, _ConvexFace >::iterator id2face = data._convexFaces.begin();
3773 for ( ; id2face != data._convexFaces.end(); ++id2face )
3775 _ConvexFace & convFace = (*id2face).second;
3776 if ( convFace._normalsFixed )
3777 continue; // already fixed
3778 if ( convFace.CheckPrisms() )
3779 continue; // nothing to fix
3781 convFace._normalsFixed = true;
3783 BRepAdaptor_Surface surface ( convFace._face, false );
3784 BRepLProp_SLProps surfProp( surface, 2, 1e-6 );
3786 // check if the convex FACE is of spherical shape
3788 Bnd_B3d centersBox; // bbox of centers of curvature of _LayerEdge's on VERTEXes
3793 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.begin();
3794 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3796 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3798 if ( meshDS->IndexToShape( id2end->first ).ShapeType() == TopAbs_VERTEX )
3800 _LayerEdge* ledge = data._edges[ iBeg ];
3801 if ( convFace.GetCenterOfCurvature( ledge, surfProp, helper, center ))
3802 centersBox.Add( center );
3804 for ( ; iBeg < iEnd; ++iBeg )
3805 nodesBox.Add( SMESH_TNodeXYZ( data._edges[ iBeg ]->_nodes[0] ));
3807 if ( centersBox.IsVoid() )
3809 debugMsg( "Error: centersBox.IsVoid()" );
3812 const bool isSpherical =
3813 ( centersBox.SquareExtent() < 1e-6 * nodesBox.SquareExtent() );
3815 int nbEdges = helper.Count( convFace._face, TopAbs_EDGE, /*ignoreSame=*/false );
3816 vector < _CentralCurveOnEdge > centerCurves( nbEdges );
3820 // set _LayerEdge::_normal as average of all normals
3822 // WARNING: different density of nodes on EDGEs is not taken into account that
3823 // can lead to an improper new normal
3825 gp_XYZ avgNormal( 0,0,0 );
3827 id2end = convFace._subIdToEdgeEnd.begin();
3828 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3830 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3831 // set data of _CentralCurveOnEdge
3832 const TopoDS_Shape& S = meshDS->IndexToShape( id2end->first );
3833 if ( S.ShapeType() == TopAbs_EDGE )
3835 _CentralCurveOnEdge& ceCurve = centerCurves[ nbEdges++ ];
3836 ceCurve.SetShapes( TopoDS::Edge(S), convFace, data, helper );
3837 if ( !data._edges[ iBeg ]->_sWOL.IsNull() )
3838 ceCurve._adjFace.Nullify();
3840 ceCurve._ledges.insert( ceCurve._ledges.end(),
3841 &data._edges[ iBeg ], &data._edges[ iEnd ]);
3843 // summarize normals
3844 for ( ; iBeg < iEnd; ++iBeg )
3845 avgNormal += data._edges[ iBeg ]->_normal;
3847 double normSize = avgNormal.SquareModulus();
3848 if ( normSize < 1e-200 )
3850 debugMsg( "updateNormalsOfConvexFaces(): zero avgNormal" );
3853 avgNormal /= Sqrt( normSize );
3855 // compute new _LayerEdge::_cosin on EDGEs
3856 double avgCosin = 0;
3859 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3861 _CentralCurveOnEdge& ceCurve = centerCurves[ iE ];
3862 if ( ceCurve._adjFace.IsNull() )
3864 for ( size_t iLE = 0; iLE < ceCurve._ledges.size(); ++iLE )
3866 const SMDS_MeshNode* node = ceCurve._ledges[ iLE ]->_nodes[0];
3867 inFaceDir = getFaceDir( ceCurve._adjFace, ceCurve._edge, node, helper, isOK );
3870 double angle = inFaceDir.Angle( avgNormal ); // [0,PI]
3871 ceCurve._ledges[ iLE ]->_cosin = Cos( angle );
3872 avgCosin += ceCurve._ledges[ iLE ]->_cosin;
3878 avgCosin /= nbCosin;
3880 // set _LayerEdge::_normal = avgNormal
3881 id2end = convFace._subIdToEdgeEnd.begin();
3882 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3884 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3885 const TopoDS_Shape& S = meshDS->IndexToShape( id2end->first );
3886 if ( S.ShapeType() != TopAbs_EDGE )
3887 for ( int i = iBeg; i < iEnd; ++i )
3888 data._edges[ i ]->_cosin = avgCosin;
3890 for ( ; iBeg < iEnd; ++iBeg )
3891 data._edges[ iBeg ]->_normal = avgNormal;
3894 else // if ( isSpherical )
3896 // We suppose that centers of curvature at all points of the FACE
3897 // lie on some curve, let's call it "central curve". For all _LayerEdge's
3898 // having a common center of curvature we define the same new normal
3899 // as a sum of normals of _LayerEdge's on EDGEs among them.
3901 // get all centers of curvature for each EDGE
3903 helper.SetSubShape( convFace._face );
3904 _LayerEdge* vertexLEdges[2], **edgeLEdge, **edgeLEdgeEnd;
3906 TopExp_Explorer edgeExp( convFace._face, TopAbs_EDGE );
3907 for ( int iE = 0; edgeExp.More(); edgeExp.Next(), ++iE )
3909 const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
3911 // set adjacent FACE
3912 centerCurves[ iE ].SetShapes( edge, convFace, data, helper );
3914 // get _LayerEdge's of the EDGE
3915 TGeomID edgeID = meshDS->ShapeToIndex( edge );
3916 id2end = convFace._subIdToEdgeEnd.find( edgeID );
3917 if ( id2end == convFace._subIdToEdgeEnd.end() )
3919 // no _LayerEdge's on EDGE, use _LayerEdge's on VERTEXes
3920 for ( int iV = 0; iV < 2; ++iV )
3922 TopoDS_Vertex v = helper.IthVertex( iV, edge );
3923 TGeomID vID = meshDS->ShapeToIndex( v );
3924 int end = convFace._subIdToEdgeEnd[ vID ];
3925 int iBeg = end > 0 ? data._endEdgeOnShape[ end-1 ] : 0;
3926 vertexLEdges[ iV ] = data._edges[ iBeg ];
3928 edgeLEdge = &vertexLEdges[0];
3929 edgeLEdgeEnd = edgeLEdge + 2;
3931 centerCurves[ iE ]._adjFace.Nullify();
3935 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3936 if ( id2end->second >= data._nbShapesToSmooth )
3937 data.SortOnEdge( edge, iBeg, iEnd, helper );
3938 edgeLEdge = &data._edges[ iBeg ];
3939 edgeLEdgeEnd = edgeLEdge + iEnd - iBeg;
3940 vertexLEdges[0] = data._edges[ iBeg ]->_2neibors->_edges[0];
3941 vertexLEdges[1] = data._edges[ iEnd-1 ]->_2neibors->_edges[1];
3943 if ( ! data._edges[ iBeg ]->_sWOL.IsNull() )
3944 centerCurves[ iE ]._adjFace.Nullify();
3947 // Get curvature centers
3951 if ( edgeLEdge[0]->IsOnEdge() &&
3952 convFace.GetCenterOfCurvature( vertexLEdges[0], surfProp, helper, center ))
3954 centerCurves[ iE ].Append( center, vertexLEdges[0] );
3955 centersBox.Add( center );
3957 for ( ; edgeLEdge < edgeLEdgeEnd; ++edgeLEdge )
3958 if ( convFace.GetCenterOfCurvature( *edgeLEdge, surfProp, helper, center ))
3959 { // EDGE or VERTEXes
3960 centerCurves[ iE ].Append( center, *edgeLEdge );
3961 centersBox.Add( center );
3963 if ( edgeLEdge[-1]->IsOnEdge() &&
3964 convFace.GetCenterOfCurvature( vertexLEdges[1], surfProp, helper, center ))
3966 centerCurves[ iE ].Append( center, vertexLEdges[1] );
3967 centersBox.Add( center );
3969 centerCurves[ iE ]._isDegenerated =
3970 ( centersBox.IsVoid() || centersBox.SquareExtent() < 1e-6 * nodesBox.SquareExtent() );
3972 } // loop on EDGES of convFace._face to set up data of centerCurves
3974 // Compute new normals for _LayerEdge's on EDGEs
3976 double avgCosin = 0;
3979 for ( size_t iE1 = 0; iE1 < centerCurves.size(); ++iE1 )
3981 _CentralCurveOnEdge& ceCurve = centerCurves[ iE1 ];
3982 if ( ceCurve._isDegenerated )
3984 const vector< gp_Pnt >& centers = ceCurve._curvaCenters;
3985 vector< gp_XYZ > & newNormals = ceCurve._normals;
3986 for ( size_t iC1 = 0; iC1 < centers.size(); ++iC1 )
3989 for ( size_t iE2 = 0; iE2 < centerCurves.size() && !isOK; ++iE2 )
3992 isOK = centerCurves[ iE2 ].FindNewNormal( centers[ iC1 ], newNormals[ iC1 ]);
3994 if ( isOK && !ceCurve._adjFace.IsNull() )
3996 // compute new _LayerEdge::_cosin
3997 const SMDS_MeshNode* node = ceCurve._ledges[ iC1 ]->_nodes[0];
3998 inFaceDir = getFaceDir( ceCurve._adjFace, ceCurve._edge, node, helper, isOK );
4001 double angle = inFaceDir.Angle( newNormals[ iC1 ] ); // [0,PI]
4002 ceCurve._ledges[ iC1 ]->_cosin = Cos( angle );
4003 avgCosin += ceCurve._ledges[ iC1 ]->_cosin;
4009 // set new normals to _LayerEdge's of NOT degenerated central curves
4010 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
4012 if ( centerCurves[ iE ]._isDegenerated )
4014 for ( size_t iLE = 0; iLE < centerCurves[ iE ]._ledges.size(); ++iLE )
4015 centerCurves[ iE ]._ledges[ iLE ]->_normal = centerCurves[ iE ]._normals[ iLE ];
4017 // set new normals to _LayerEdge's of degenerated central curves
4018 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
4020 if ( !centerCurves[ iE ]._isDegenerated ||
4021 centerCurves[ iE ]._ledges.size() < 3 )
4023 // new normal is an average of new normals at VERTEXes that
4024 // was computed on non-degenerated _CentralCurveOnEdge's
4025 gp_XYZ newNorm = ( centerCurves[ iE ]._ledges.front()->_normal +
4026 centerCurves[ iE ]._ledges.back ()->_normal );
4027 double sz = newNorm.Modulus();
4031 double newCosin = ( 0.5 * centerCurves[ iE ]._ledges.front()->_cosin +
4032 0.5 * centerCurves[ iE ]._ledges.back ()->_cosin );
4033 for ( size_t iLE = 1, nb = centerCurves[ iE ]._ledges.size() - 1; iLE < nb; ++iLE )
4035 centerCurves[ iE ]._ledges[ iLE ]->_normal = newNorm;
4036 centerCurves[ iE ]._ledges[ iLE ]->_cosin = newCosin;
4040 // Find new normals for _LayerEdge's based on FACE
4043 avgCosin /= nbCosin;
4044 const TGeomID faceID = meshDS->ShapeToIndex( convFace._face );
4045 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.find( faceID );
4046 if ( id2end != convFace._subIdToEdgeEnd.end() )
4050 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
4051 for ( ; iBeg < iEnd; ++iBeg )
4053 _LayerEdge* ledge = data._edges[ iBeg ];
4054 if ( !convFace.GetCenterOfCurvature( ledge, surfProp, helper, center ))
4056 for ( size_t i = 0; i < centerCurves.size(); ++i, ++iE )
4058 iE = iE % centerCurves.size();
4059 if ( centerCurves[ iE ]._isDegenerated )
4061 newNorm.SetCoord( 0,0,0 );
4062 if ( centerCurves[ iE ].FindNewNormal( center, newNorm ))
4064 ledge->_normal = newNorm;
4065 ledge->_cosin = avgCosin;
4072 } // not a quasi-spherical FACE
4074 // Update _LayerEdge's data according to a new normal
4076 dumpFunction(SMESH_Comment("updateNormalsOfConvexFaces")<<data._index
4077 <<"_F"<<meshDS->ShapeToIndex( convFace._face ));
4079 id2end = convFace._subIdToEdgeEnd.begin();
4080 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
4082 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
4083 for ( ; iBeg < iEnd; ++iBeg )
4085 _LayerEdge* & ledge = data._edges[ iBeg ];
4086 double len = ledge->_len;
4087 ledge->InvalidateStep( stepNb + 1, /*restoreLength=*/true );
4088 ledge->SetCosin( ledge->_cosin );
4089 ledge->SetNewLength( len, helper );
4092 } // loop on sub-shapes of convFace._face
4094 // Find FACEs adjacent to convFace._face that got necessity to smooth
4095 // as a result of normals modification
4097 set< TGeomID > adjFacesToSmooth;
4098 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
4100 if ( centerCurves[ iE ]._adjFace.IsNull() ||
4101 centerCurves[ iE ]._adjFaceToSmooth )
4103 for ( size_t iLE = 0; iLE < centerCurves[ iE ]._ledges.size(); ++iLE )
4105 if ( centerCurves[ iE ]._ledges[ iLE ]->_cosin > theMinSmoothCosin )
4107 adjFacesToSmooth.insert( meshDS->ShapeToIndex( centerCurves[ iE ]._adjFace ));
4112 data.AddShapesToSmooth( adjFacesToSmooth );
4117 } // loop on data._convexFaces
4122 //================================================================================
4124 * \brief Finds a center of curvature of a surface at a _LayerEdge
4126 //================================================================================
4128 bool _ConvexFace::GetCenterOfCurvature( _LayerEdge* ledge,
4129 BRepLProp_SLProps& surfProp,
4130 SMESH_MesherHelper& helper,
4131 gp_Pnt & center ) const
4133 gp_XY uv = helper.GetNodeUV( _face, ledge->_nodes[0] );
4134 surfProp.SetParameters( uv.X(), uv.Y() );
4135 if ( !surfProp.IsCurvatureDefined() )
4138 const double oriFactor = ( _face.Orientation() == TopAbs_REVERSED ? +1. : -1. );
4139 double surfCurvatureMax = surfProp.MaxCurvature() * oriFactor;
4140 double surfCurvatureMin = surfProp.MinCurvature() * oriFactor;
4141 if ( surfCurvatureMin > surfCurvatureMax )
4142 center = surfProp.Value().Translated( surfProp.Normal().XYZ() / surfCurvatureMin * oriFactor );
4144 center = surfProp.Value().Translated( surfProp.Normal().XYZ() / surfCurvatureMax * oriFactor );
4149 //================================================================================
4151 * \brief Check that prisms are not distorted
4153 //================================================================================
4155 bool _ConvexFace::CheckPrisms() const
4157 for ( size_t i = 0; i < _simplexTestEdges.size(); ++i )
4159 const _LayerEdge* edge = _simplexTestEdges[i];
4160 SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
4161 for ( size_t j = 0; j < edge->_simplices.size(); ++j )
4162 if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
4164 debugMsg( "Bad simplex of _simplexTestEdges ("
4165 << " "<< edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
4166 << " "<< edge->_simplices[j]._nPrev->GetID()
4167 << " "<< edge->_simplices[j]._nNext->GetID() << " )" );
4174 //================================================================================
4176 * \brief Try to compute a new normal by interpolating normals of _LayerEdge's
4177 * stored in this _CentralCurveOnEdge.
4178 * \param [in] center - curvature center of a point of another _CentralCurveOnEdge.
4179 * \param [in,out] newNormal - current normal at this point, to be redefined
4180 * \return bool - true if succeeded.
4182 //================================================================================
4184 bool _CentralCurveOnEdge::FindNewNormal( const gp_Pnt& center, gp_XYZ& newNormal )
4186 if ( this->_isDegenerated )
4189 // find two centers the given one lies between
4191 for ( size_t i = 0, nb = _curvaCenters.size()-1; i < nb; ++i )
4193 double sl2 = 1.001 * _segLength2[ i ];
4195 double d1 = center.SquareDistance( _curvaCenters[ i ]);
4199 double d2 = center.SquareDistance( _curvaCenters[ i+1 ]);
4200 if ( d2 > sl2 || d2 + d1 < 1e-100 )
4205 double r = d1 / ( d1 + d2 );
4206 gp_XYZ norm = (( 1. - r ) * _ledges[ i ]->_normal +
4207 ( r ) * _ledges[ i+1 ]->_normal );
4211 double sz = newNormal.Modulus();
4220 //================================================================================
4222 * \brief Set shape members
4224 //================================================================================
4226 void _CentralCurveOnEdge::SetShapes( const TopoDS_Edge& edge,
4227 const _ConvexFace& convFace,
4228 const _SolidData& data,
4229 SMESH_MesherHelper& helper)
4233 PShapeIteratorPtr fIt = helper.GetAncestors( edge, *helper.GetMesh(), TopAbs_FACE );
4234 while ( const TopoDS_Shape* F = fIt->next())
4235 if ( !convFace._face.IsSame( *F ))
4237 _adjFace = TopoDS::Face( *F );
4238 _adjFaceToSmooth = false;
4239 // _adjFace already in a smoothing queue ?
4241 TGeomID adjFaceID = helper.GetMeshDS()->ShapeToIndex( *F );
4242 if ( data.GetShapeEdges( adjFaceID, end ))
4243 _adjFaceToSmooth = ( end < data._nbShapesToSmooth );
4248 //================================================================================
4250 * \brief Looks for intersection of it's last segment with faces
4251 * \param distance - returns shortest distance from the last node to intersection
4253 //================================================================================
4255 bool _LayerEdge::FindIntersection( SMESH_ElementSearcher& searcher,
4257 const double& epsilon,
4258 const SMDS_MeshElement** face)
4260 vector< const SMDS_MeshElement* > suspectFaces;
4262 gp_Ax1 lastSegment = LastSegment(segLen);
4263 searcher.GetElementsNearLine( lastSegment, SMDSAbs_Face, suspectFaces );
4265 bool segmentIntersected = false;
4266 distance = Precision::Infinite();
4267 int iFace = -1; // intersected face
4268 for ( size_t j = 0 ; j < suspectFaces.size() /*&& !segmentIntersected*/; ++j )
4270 const SMDS_MeshElement* face = suspectFaces[j];
4271 if ( face->GetNodeIndex( _nodes.back() ) >= 0 ||
4272 face->GetNodeIndex( _nodes[0] ) >= 0 )
4273 continue; // face sharing _LayerEdge node
4274 const int nbNodes = face->NbCornerNodes();
4275 bool intFound = false;
4277 SMDS_MeshElement::iterator nIt = face->begin_nodes();
4280 intFound = SegTriaInter( lastSegment, *nIt++, *nIt++, *nIt++, dist, epsilon );
4284 const SMDS_MeshNode* tria[3];
4287 for ( int n2 = 2; n2 < nbNodes && !intFound; ++n2 )
4290 intFound = SegTriaInter(lastSegment, tria[0], tria[1], tria[2], dist, epsilon );
4296 if ( dist < segLen*(1.01) && dist > -(_len*_lenFactor-segLen) )
4297 segmentIntersected = true;
4298 if ( distance > dist )
4299 distance = dist, iFace = j;
4302 if ( iFace != -1 && face ) *face = suspectFaces[iFace];
4304 if ( segmentIntersected )
4307 SMDS_MeshElement::iterator nIt = suspectFaces[iFace]->begin_nodes();
4308 gp_XYZ intP( lastSegment.Location().XYZ() + lastSegment.Direction().XYZ() * distance );
4309 cout << "nodes: tgt " << _nodes.back()->GetID() << " src " << _nodes[0]->GetID()
4310 << ", intersection with face ("
4311 << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
4312 << ") at point (" << intP.X() << ", " << intP.Y() << ", " << intP.Z()
4313 << ") distance = " << distance - segLen<< endl;
4319 return segmentIntersected;
4322 //================================================================================
4324 * \brief Returns size and direction of the last segment
4326 //================================================================================
4328 gp_Ax1 _LayerEdge::LastSegment(double& segLen) const
4330 // find two non-coincident positions
4331 gp_XYZ orig = _pos.back();
4333 int iPrev = _pos.size() - 2;
4334 while ( iPrev >= 0 )
4336 dir = orig - _pos[iPrev];
4337 if ( dir.SquareModulus() > 1e-100 )
4347 segDir.SetLocation( SMESH_TNodeXYZ( _nodes[0] ));
4348 segDir.SetDirection( _normal );
4353 gp_Pnt pPrev = _pos[ iPrev ];
4354 if ( !_sWOL.IsNull() )
4356 TopLoc_Location loc;
4357 if ( _sWOL.ShapeType() == TopAbs_EDGE )
4360 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
4361 pPrev = curve->Value( pPrev.X() ).Transformed( loc );
4365 Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
4366 pPrev = surface->Value( pPrev.X(), pPrev.Y() ).Transformed( loc );
4368 dir = SMESH_TNodeXYZ( _nodes.back() ) - pPrev.XYZ();
4370 segDir.SetLocation( pPrev );
4371 segDir.SetDirection( dir );
4372 segLen = dir.Modulus();
4378 //================================================================================
4380 * \brief Return the last position of the target node on a FACE.
4381 * \param [in] F - the FACE this _LayerEdge is inflated along
4382 * \return gp_XY - result UV
4384 //================================================================================
4386 gp_XY _LayerEdge::LastUV( const TopoDS_Face& F ) const
4388 if ( F.IsSame( _sWOL )) // F is my FACE
4389 return gp_XY( _pos.back().X(), _pos.back().Y() );
4391 if ( _sWOL.IsNull() || _sWOL.ShapeType() != TopAbs_EDGE ) // wrong call
4392 return gp_XY( 1e100, 1e100 );
4394 // _sWOL is EDGE of F; _pos.back().X() is the last U on the EDGE
4395 double f, l, u = _pos.back().X();
4396 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( TopoDS::Edge(_sWOL), F, f,l);
4397 if ( !C2d.IsNull() && f <= u && u <= l )
4398 return C2d->Value( u ).XY();
4400 return gp_XY( 1e100, 1e100 );
4403 //================================================================================
4405 * \brief Test intersection of the last segment with a given triangle
4406 * using Moller-Trumbore algorithm
4407 * Intersection is detected if distance to intersection is less than _LayerEdge._len
4409 //================================================================================
4411 bool _LayerEdge::SegTriaInter( const gp_Ax1& lastSegment,
4412 const SMDS_MeshNode* n0,
4413 const SMDS_MeshNode* n1,
4414 const SMDS_MeshNode* n2,
4416 const double& EPSILON) const
4418 //const double EPSILON = 1e-6;
4420 gp_XYZ orig = lastSegment.Location().XYZ();
4421 gp_XYZ dir = lastSegment.Direction().XYZ();
4423 SMESH_TNodeXYZ vert0( n0 );
4424 SMESH_TNodeXYZ vert1( n1 );
4425 SMESH_TNodeXYZ vert2( n2 );
4427 /* calculate distance from vert0 to ray origin */
4428 gp_XYZ tvec = orig - vert0;
4430 //if ( tvec * dir > EPSILON )
4431 // intersected face is at back side of the temporary face this _LayerEdge belongs to
4434 gp_XYZ edge1 = vert1 - vert0;
4435 gp_XYZ edge2 = vert2 - vert0;
4437 /* begin calculating determinant - also used to calculate U parameter */
4438 gp_XYZ pvec = dir ^ edge2;
4440 /* if determinant is near zero, ray lies in plane of triangle */
4441 double det = edge1 * pvec;
4443 if (det > -EPSILON && det < EPSILON)
4445 double inv_det = 1.0 / det;
4447 /* calculate U parameter and test bounds */
4448 double u = ( tvec * pvec ) * inv_det;
4449 //if (u < 0.0 || u > 1.0)
4450 if (u < -EPSILON || u > 1.0 + EPSILON)
4453 /* prepare to test V parameter */
4454 gp_XYZ qvec = tvec ^ edge1;
4456 /* calculate V parameter and test bounds */
4457 double v = (dir * qvec) * inv_det;
4458 //if ( v < 0.0 || u + v > 1.0 )
4459 if ( v < -EPSILON || u + v > 1.0 + EPSILON)
4462 /* calculate t, ray intersects triangle */
4463 t = (edge2 * qvec) * inv_det;
4469 //================================================================================
4471 * \brief Perform smooth of _LayerEdge's based on EDGE's
4472 * \retval bool - true if node has been moved
4474 //================================================================================
4476 bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface,
4477 const TopoDS_Face& F,
4478 SMESH_MesherHelper& helper)
4480 ASSERT( IsOnEdge() );
4482 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( _nodes.back() );
4483 SMESH_TNodeXYZ oldPos( tgtNode );
4484 double dist01, distNewOld;
4486 SMESH_TNodeXYZ p0( _2neibors->tgtNode(0));
4487 SMESH_TNodeXYZ p1( _2neibors->tgtNode(1));
4488 dist01 = p0.Distance( _2neibors->tgtNode(1) );
4490 gp_Pnt newPos = p0 * _2neibors->_wgt[0] + p1 * _2neibors->_wgt[1];
4491 double lenDelta = 0;
4494 //lenDelta = _curvature->lenDelta( _len );
4495 lenDelta = _curvature->lenDeltaByDist( dist01 );
4496 newPos.ChangeCoord() += _normal * lenDelta;
4499 distNewOld = newPos.Distance( oldPos );
4503 if ( _2neibors->_plnNorm )
4505 // put newPos on the plane defined by source node and _plnNorm
4506 gp_XYZ new2src = SMESH_TNodeXYZ( _nodes[0] ) - newPos.XYZ();
4507 double new2srcProj = (*_2neibors->_plnNorm) * new2src;
4508 newPos.ChangeCoord() += (*_2neibors->_plnNorm) * new2srcProj;
4510 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4511 _pos.back() = newPos.XYZ();
4515 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4516 gp_XY uv( Precision::Infinite(), 0 );
4517 helper.CheckNodeUV( F, tgtNode, uv, 1e-10, /*force=*/true );
4518 _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
4520 newPos = surface->Value( uv.X(), uv.Y() );
4521 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4524 if ( _curvature && lenDelta < 0 )
4526 gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
4527 _len -= prevPos.Distance( oldPos );
4528 _len += prevPos.Distance( newPos );
4530 bool moved = distNewOld > dist01/50;
4532 dumpMove( tgtNode ); // debug
4537 //================================================================================
4539 * \brief Perform laplacian smooth in 3D of nodes inflated from FACE
4540 * \retval bool - true if _tgtNode has been moved
4542 //================================================================================
4544 bool _LayerEdge::Smooth(int& badNb)
4546 if ( _simplices.size() < 2 )
4547 return false; // _LayerEdge inflated along EDGE or FACE
4549 // compute new position for the last _pos
4550 gp_XYZ newPos (0,0,0);
4551 for ( size_t i = 0; i < _simplices.size(); ++i )
4552 newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev );
4553 newPos /= _simplices.size();
4555 const gp_XYZ& curPos ( _pos.back() );
4556 const gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
4559 double delta = _curvature->lenDelta( _len );
4561 newPos += _normal * delta;
4564 double segLen = _normal * ( newPos - prevPos.XYZ() );
4565 if ( segLen + delta > 0 )
4566 newPos += _normal * delta;
4568 // double segLenChange = _normal * ( curPos - newPos );
4569 // newPos += 0.5 * _normal * segLenChange;
4572 // count quality metrics (orientation) of tetras around _tgtNode
4574 for ( size_t i = 0; i < _simplices.size(); ++i )
4575 nbOkBefore += _simplices[i].IsForward( _nodes[0], &curPos );
4578 for ( size_t i = 0; i < _simplices.size(); ++i )
4579 nbOkAfter += _simplices[i].IsForward( _nodes[0], &newPos );
4581 if ( nbOkAfter < nbOkBefore )
4584 SMDS_MeshNode* n = const_cast< SMDS_MeshNode* >( _nodes.back() );
4586 _len -= prevPos.Distance(SMESH_TNodeXYZ( n ));
4587 _len += prevPos.Distance(newPos);
4589 n->setXYZ( newPos.X(), newPos.Y(), newPos.Z());
4590 _pos.back() = newPos;
4592 badNb += _simplices.size() - nbOkAfter;
4599 //================================================================================
4601 * \brief Add a new segment to _LayerEdge during inflation
4603 //================================================================================
4605 void _LayerEdge::SetNewLength( double len, SMESH_MesherHelper& helper )
4607 if ( _len - len > -1e-6 )
4609 _pos.push_back( _pos.back() );
4613 SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
4614 SMESH_TNodeXYZ oldXYZ( n );
4615 gp_XYZ nXYZ = oldXYZ + _normal * ( len - _len ) * _lenFactor;
4616 n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
4618 _pos.push_back( nXYZ );
4620 if ( !_sWOL.IsNull() )
4623 if ( _sWOL.ShapeType() == TopAbs_EDGE )
4625 double u = Precision::Infinite(); // to force projection w/o distance check
4626 helper.CheckNodeU( TopoDS::Edge( _sWOL ), n, u, 1e-10, /*force=*/true, distXYZ );
4627 _pos.back().SetCoord( u, 0, 0 );
4628 if ( _nodes.size() > 1 )
4630 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
4631 pos->SetUParameter( u );
4636 gp_XY uv( Precision::Infinite(), 0 );
4637 helper.CheckNodeUV( TopoDS::Face( _sWOL ), n, uv, 1e-10, /*force=*/true, distXYZ );
4638 _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
4639 if ( _nodes.size() > 1 )
4641 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
4642 pos->SetUParameter( uv.X() );
4643 pos->SetVParameter( uv.Y() );
4646 n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]);
4648 dumpMove( n ); //debug
4651 //================================================================================
4653 * \brief Remove last inflation step
4655 //================================================================================
4657 void _LayerEdge::InvalidateStep( int curStep, bool restoreLength )
4659 if ( _pos.size() > curStep )
4661 if ( restoreLength )
4662 _len -= ( _pos[ curStep-1 ] - _pos.back() ).Modulus();
4664 _pos.resize( curStep );
4665 gp_Pnt nXYZ = _pos.back();
4666 SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
4667 if ( !_sWOL.IsNull() )
4669 TopLoc_Location loc;
4670 if ( _sWOL.ShapeType() == TopAbs_EDGE )
4672 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
4673 pos->SetUParameter( nXYZ.X() );
4675 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
4676 nXYZ = curve->Value( nXYZ.X() ).Transformed( loc );
4680 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
4681 pos->SetUParameter( nXYZ.X() );
4682 pos->SetVParameter( nXYZ.Y() );
4683 Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
4684 nXYZ = surface->Value( nXYZ.X(), nXYZ.Y() ).Transformed( loc );
4687 n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
4692 //================================================================================
4694 * \brief Create layers of prisms
4696 //================================================================================
4698 bool _ViscousBuilder::refine(_SolidData& data)
4700 SMESH_MesherHelper helper( *_mesh );
4701 helper.SetSubShape( data._solid );
4702 helper.SetElementsOnShape(false);
4704 Handle(Geom_Curve) curve;
4705 Handle(Geom_Surface) surface;
4706 TopoDS_Edge geomEdge;
4707 TopoDS_Face geomFace;
4708 TopoDS_Shape prevSWOL;
4709 TopLoc_Location loc;
4713 TGeomID prevBaseId = -1;
4714 TNode2Edge* n2eMap = 0;
4715 TNode2Edge::iterator n2e;
4717 // Create intermediate nodes on each _LayerEdge
4719 for ( size_t i = 0; i < data._edges.size(); ++i )
4721 _LayerEdge& edge = *data._edges[i];
4723 if ( edge._nodes.size() < 2 )
4724 continue; // on _noShrinkShapes
4726 // get accumulated length of segments
4727 vector< double > segLen( edge._pos.size() );
4729 for ( size_t j = 1; j < edge._pos.size(); ++j )
4730 segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus();
4732 // allocate memory for new nodes if it is not yet refined
4733 const SMDS_MeshNode* tgtNode = edge._nodes.back();
4734 if ( edge._nodes.size() == 2 )
4736 edge._nodes.resize( data._hyp->GetNumberLayers() + 1, 0 );
4738 edge._nodes.back() = tgtNode;
4740 // get data of a shrink shape
4741 if ( !edge._sWOL.IsNull() && edge._sWOL != prevSWOL )
4743 isOnEdge = ( edge._sWOL.ShapeType() == TopAbs_EDGE );
4746 geomEdge = TopoDS::Edge( edge._sWOL );
4747 curve = BRep_Tool::Curve( geomEdge, loc, f,l);
4751 geomFace = TopoDS::Face( edge._sWOL );
4752 surface = BRep_Tool::Surface( geomFace, loc );
4754 prevSWOL = edge._sWOL;
4756 // restore shapePos of the last node by already treated _LayerEdge of another _SolidData
4757 const TGeomID baseShapeId = edge._nodes[0]->getshapeId();
4758 if ( baseShapeId != prevBaseId )
4760 map< TGeomID, TNode2Edge* >::iterator s2ne = data._s2neMap.find( baseShapeId );
4761 n2eMap = ( s2ne == data._s2neMap.end() ) ? 0 : n2eMap = s2ne->second;
4762 prevBaseId = baseShapeId;
4764 _LayerEdge* edgeOnSameNode = 0;
4765 if ( n2eMap && (( n2e = n2eMap->find( edge._nodes[0] )) != n2eMap->end() ))
4767 edgeOnSameNode = n2e->second;
4768 const gp_XYZ& otherTgtPos = edgeOnSameNode->_pos.back();
4769 SMDS_PositionPtr lastPos = tgtNode->GetPosition();
4772 SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( lastPos );
4773 epos->SetUParameter( otherTgtPos.X() );
4777 SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( lastPos );
4778 fpos->SetUParameter( otherTgtPos.X() );
4779 fpos->SetVParameter( otherTgtPos.Y() );
4782 // calculate height of the first layer
4784 const double T = segLen.back(); //data._hyp.GetTotalThickness();
4785 const double f = data._hyp->GetStretchFactor();
4786 const int N = data._hyp->GetNumberLayers();
4787 const double fPowN = pow( f, N );
4788 if ( fPowN - 1 <= numeric_limits<double>::min() )
4791 h0 = T * ( f - 1 )/( fPowN - 1 );
4793 const double zeroLen = std::numeric_limits<double>::min();
4795 // create intermediate nodes
4796 double hSum = 0, hi = h0/f;
4798 for ( size_t iStep = 1; iStep < edge._nodes.size(); ++iStep )
4800 // compute an intermediate position
4803 while ( hSum > segLen[iSeg] && iSeg < segLen.size()-1)
4805 int iPrevSeg = iSeg-1;
4806 while ( fabs( segLen[iPrevSeg] - segLen[iSeg]) <= zeroLen && iPrevSeg > 0 )
4808 double r = ( segLen[iSeg] - hSum ) / ( segLen[iSeg] - segLen[iPrevSeg] );
4809 gp_Pnt pos = r * edge._pos[iPrevSeg] + (1-r) * edge._pos[iSeg];
4811 SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >( edge._nodes[ iStep ]);
4812 if ( !edge._sWOL.IsNull() )
4814 // compute XYZ by parameters <pos>
4819 pos = curve->Value( u ).Transformed(loc);
4823 uv.SetCoord( pos.X(), pos.Y() );
4825 pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc);
4828 // create or update the node
4831 node = helper.AddNode( pos.X(), pos.Y(), pos.Z());
4832 if ( !edge._sWOL.IsNull() )
4835 getMeshDS()->SetNodeOnEdge( node, geomEdge, u );
4837 getMeshDS()->SetNodeOnFace( node, geomFace, uv.X(), uv.Y() );
4841 getMeshDS()->SetNodeInVolume( node, helper.GetSubShapeID() );
4846 if ( !edge._sWOL.IsNull() )
4848 // make average pos from new and current parameters
4851 u = 0.5 * ( u + helper.GetNodeU( geomEdge, node ));
4852 pos = curve->Value( u ).Transformed(loc);
4854 SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( node->GetPosition() );
4855 epos->SetUParameter( u );
4859 uv = 0.5 * ( uv + helper.GetNodeUV( geomFace, node ));
4860 pos = surface->Value( uv.X(), uv.Y()).Transformed(loc);
4862 SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( node->GetPosition() );
4863 fpos->SetUParameter( uv.X() );
4864 fpos->SetVParameter( uv.Y() );
4867 node->setXYZ( pos.X(), pos.Y(), pos.Z() );
4869 } // loop on edge._nodes
4871 if ( !edge._sWOL.IsNull() ) // prepare for shrink()
4874 edge._pos.back().SetCoord( u, 0,0);
4876 edge._pos.back().SetCoord( uv.X(), uv.Y() ,0);
4878 if ( edgeOnSameNode )
4879 edgeOnSameNode->_pos.back() = edge._pos.back();
4882 } // loop on data._edges to create nodes
4884 if ( !getMeshDS()->IsEmbeddedMode() )
4885 // Log node movement
4886 for ( size_t i = 0; i < data._edges.size(); ++i )
4888 _LayerEdge& edge = *data._edges[i];
4889 SMESH_TNodeXYZ p ( edge._nodes.back() );
4890 getMeshDS()->MoveNode( p._node, p.X(), p.Y(), p.Z() );
4895 helper.SetElementsOnShape(true);
4897 vector< vector<const SMDS_MeshNode*>* > nnVec;
4898 set< vector<const SMDS_MeshNode*>* > nnSet;
4899 set< int > degenEdgeInd;
4900 vector<const SMDS_MeshElement*> degenVols;
4902 TopExp_Explorer exp( data._solid, TopAbs_FACE );
4903 for ( ; exp.More(); exp.Next() )
4905 if ( data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
4907 SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
4908 SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
4909 while ( fIt->more() )
4911 const SMDS_MeshElement* face = fIt->next();
4912 const int nbNodes = face->NbCornerNodes();
4913 nnVec.resize( nbNodes );
4915 degenEdgeInd.clear();
4917 SMDS_NodeIteratorPtr nIt = face->nodeIterator();
4918 for ( int iN = 0; iN < nbNodes; ++iN )
4920 const SMDS_MeshNode* n = nIt->next();
4921 nnVec[ iN ] = & data._n2eMap[ n ]->_nodes;
4922 if ( nnVec[ iN ]->size() < 2 )
4923 degenEdgeInd.insert( iN );
4925 nbZ = nnVec[ iN ]->size();
4927 if ( helper.HasDegeneratedEdges() )
4928 nnSet.insert( nnVec[ iN ]);
4932 if ( 0 < nnSet.size() && nnSet.size() < 3 )
4938 switch ( degenEdgeInd.size() )
4942 for ( int iZ = 1; iZ < nbZ; ++iZ )
4943 helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
4944 (*nnVec[0])[iZ], (*nnVec[1])[iZ], (*nnVec[2])[iZ]);
4949 int i2 = *degenEdgeInd.begin();
4950 int i0 = helper.WrapIndex( i2 - 1, nbNodes );
4951 int i1 = helper.WrapIndex( i2 + 1, nbNodes );
4952 for ( int iZ = 1; iZ < nbZ; ++iZ )
4953 helper.AddVolume( (*nnVec[i0])[iZ-1], (*nnVec[i1])[iZ-1],
4954 (*nnVec[i1])[iZ], (*nnVec[i0])[iZ], (*nnVec[i2])[0]);
4959 int i3 = !degenEdgeInd.count(0) ? 0 : !degenEdgeInd.count(1) ? 1 : 2;
4960 for ( int iZ = 1; iZ < nbZ; ++iZ )
4961 helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
4969 switch ( degenEdgeInd.size() )
4973 for ( int iZ = 1; iZ < nbZ; ++iZ )
4974 helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1],
4975 (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1],
4976 (*nnVec[0])[iZ], (*nnVec[1])[iZ],
4977 (*nnVec[2])[iZ], (*nnVec[3])[iZ]);
4982 int i2 = *degenEdgeInd.begin();
4983 int i3 = *degenEdgeInd.rbegin();
4984 bool ok = ( i3 - i2 == 1 );
4985 if ( i2 == 0 && i3 == 3 ) { i2 = 3; i3 = 0; ok = true; }
4986 int i0 = helper.WrapIndex( i3 + 1, nbNodes );
4987 int i1 = helper.WrapIndex( i0 + 1, nbNodes );
4988 for ( int iZ = 1; iZ < nbZ; ++iZ )
4990 const SMDS_MeshElement* vol =
4991 helper.AddVolume( (*nnVec[i3])[0], (*nnVec[i0])[iZ], (*nnVec[i0])[iZ-1],
4992 (*nnVec[i2])[0], (*nnVec[i1])[iZ], (*nnVec[i1])[iZ-1]);
4994 degenVols.push_back( vol );
4998 case 3: // degen HEX
5000 const SMDS_MeshNode* nn[8];
5001 for ( int iZ = 1; iZ < nbZ; ++iZ )
5003 const SMDS_MeshElement* vol =
5004 helper.AddVolume( nnVec[0]->size() > 1 ? (*nnVec[0])[iZ-1] : (*nnVec[0])[0],
5005 nnVec[1]->size() > 1 ? (*nnVec[1])[iZ-1] : (*nnVec[1])[0],
5006 nnVec[2]->size() > 1 ? (*nnVec[2])[iZ-1] : (*nnVec[2])[0],
5007 nnVec[3]->size() > 1 ? (*nnVec[3])[iZ-1] : (*nnVec[3])[0],
5008 nnVec[0]->size() > 1 ? (*nnVec[0])[iZ] : (*nnVec[0])[0],
5009 nnVec[1]->size() > 1 ? (*nnVec[1])[iZ] : (*nnVec[1])[0],
5010 nnVec[2]->size() > 1 ? (*nnVec[2])[iZ] : (*nnVec[2])[0],
5011 nnVec[3]->size() > 1 ? (*nnVec[3])[iZ] : (*nnVec[3])[0]);
5012 degenVols.push_back( vol );
5020 return error("Not supported type of element", data._index);
5022 } // switch ( nbNodes )
5023 } // while ( fIt->more() )
5026 if ( !degenVols.empty() )
5028 SMESH_ComputeErrorPtr& err = _mesh->GetSubMesh( data._solid )->GetComputeError();
5029 if ( !err || err->IsOK() )
5031 err.reset( new SMESH_ComputeError( COMPERR_WARNING,
5032 "Degenerated volumes created" ));
5033 err->myBadElements.insert( err->myBadElements.end(),
5034 degenVols.begin(),degenVols.end() );
5041 //================================================================================
5043 * \brief Shrink 2D mesh on faces to let space for inflated layers
5045 //================================================================================
5047 bool _ViscousBuilder::shrink()
5049 // make map of (ids of FACEs to shrink mesh on) to (_SolidData containing _LayerEdge's
5050 // inflated along FACE or EDGE)
5051 map< TGeomID, _SolidData* > f2sdMap;
5052 for ( size_t i = 0 ; i < _sdVec.size(); ++i )
5054 _SolidData& data = _sdVec[i];
5055 TopTools_MapOfShape FFMap;
5056 map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
5057 for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
5058 if ( s2s->second.ShapeType() == TopAbs_FACE )
5060 f2sdMap.insert( make_pair( getMeshDS()->ShapeToIndex( s2s->second ), &data ));
5062 if ( FFMap.Add( (*s2s).second ))
5063 // Put mesh faces on the shrinked FACE to the proxy sub-mesh to avoid
5064 // usage of mesh faces made in addBoundaryElements() by the 3D algo or
5065 // by StdMeshers_QuadToTriaAdaptor
5066 if ( SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( s2s->second ))
5068 SMESH_ProxyMesh::SubMesh* proxySub =
5069 data._proxyMesh->getFaceSubM( TopoDS::Face( s2s->second ), /*create=*/true);
5070 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
5071 while ( fIt->more() )
5072 proxySub->AddElement( fIt->next() );
5073 // as a result 3D algo will use elements from proxySub and not from smDS
5078 SMESH_MesherHelper helper( *_mesh );
5079 helper.ToFixNodeParameters( true );
5082 map< TGeomID, _Shrinker1D > e2shrMap;
5083 vector< _LayerEdge* > lEdges;
5085 // loop on FACES to srink mesh on
5086 map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin();
5087 for ( ; f2sd != f2sdMap.end(); ++f2sd )
5089 _SolidData& data = *f2sd->second;
5090 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first ));
5091 SMESH_subMesh* sm = _mesh->GetSubMesh( F );
5092 SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
5094 Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
5096 helper.SetSubShape(F);
5098 // ===========================
5099 // Prepare data for shrinking
5100 // ===========================
5102 // Collect nodes to smooth, as src nodes are not yet replaced by tgt ones
5103 // and thus all nodes on a FACE connected to 2d elements are to be smoothed
5104 vector < const SMDS_MeshNode* > smoothNodes;
5106 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
5107 while ( nIt->more() )
5109 const SMDS_MeshNode* n = nIt->next();
5110 if ( n->NbInverseElements( SMDSAbs_Face ) > 0 )
5111 smoothNodes.push_back( n );
5114 // Find out face orientation
5116 const set<TGeomID> ignoreShapes;
5118 if ( !smoothNodes.empty() )
5120 vector<_Simplex> simplices;
5121 getSimplices( smoothNodes[0], simplices, ignoreShapes );
5122 helper.GetNodeUV( F, simplices[0]._nPrev, 0, &isOkUV ); // fix UV of silpmex nodes
5123 helper.GetNodeUV( F, simplices[0]._nNext, 0, &isOkUV );
5124 gp_XY uv = helper.GetNodeUV( F, smoothNodes[0], 0, &isOkUV );
5125 if ( !simplices[0].IsForward(uv, smoothNodes[0], F, helper,refSign) )
5129 // Find _LayerEdge's inflated along F
5132 set< TGeomID > subIDs;
5133 SMESH_subMeshIteratorPtr subIt = sm->getDependsOnIterator(/*includeSelf=*/false);
5134 while ( subIt->more() )
5135 subIDs.insert( subIt->next()->GetId() );
5138 for ( int iS = 0; iS < data._endEdgeOnShape.size() && !subIDs.empty(); ++iS )
5141 iEnd = data._endEdgeOnShape[ iS ];
5142 TGeomID shapeID = data._edges[ iBeg ]->_nodes[0]->getshapeId();
5143 set< TGeomID >::iterator idIt = subIDs.find( shapeID );
5144 if ( idIt == subIDs.end() ||
5145 data._edges[ iBeg ]->_sWOL.IsNull() ) continue;
5146 subIDs.erase( idIt );
5148 if ( !data._noShrinkShapes.count( shapeID ))
5149 for ( ; iBeg < iEnd; ++iBeg )
5151 lEdges.push_back( data._edges[ iBeg ] );
5152 prepareEdgeToShrink( *data._edges[ iBeg ], F, helper, smDS );
5157 dumpFunction(SMESH_Comment("beforeShrinkFace")<<f2sd->first); // debug
5158 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
5159 while ( fIt->more() )
5160 if ( const SMDS_MeshElement* f = fIt->next() )
5161 dumpChangeNodes( f );
5163 // Replace source nodes by target nodes in mesh faces to shrink
5164 dumpFunction(SMESH_Comment("replNodesOnFace")<<f2sd->first); // debug
5165 const SMDS_MeshNode* nodes[20];
5166 for ( size_t i = 0; i < lEdges.size(); ++i )
5168 _LayerEdge& edge = *lEdges[i];
5169 const SMDS_MeshNode* srcNode = edge._nodes[0];
5170 const SMDS_MeshNode* tgtNode = edge._nodes.back();
5171 SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
5172 while ( fIt->more() )
5174 const SMDS_MeshElement* f = fIt->next();
5175 if ( !smDS->Contains( f ))
5177 SMDS_NodeIteratorPtr nIt = f->nodeIterator();
5178 for ( int iN = 0; nIt->more(); ++iN )
5180 const SMDS_MeshNode* n = nIt->next();
5181 nodes[iN] = ( n == srcNode ? tgtNode : n );
5183 helper.GetMeshDS()->ChangeElementNodes( f, nodes, f->NbNodes() );
5184 dumpChangeNodes( f );
5188 // find out if a FACE is concave
5189 const bool isConcaveFace = isConcave( F, helper );
5191 // Create _SmoothNode's on face F
5192 vector< _SmoothNode > nodesToSmooth( smoothNodes.size() );
5194 dumpFunction(SMESH_Comment("fixUVOnFace")<<f2sd->first); // debug
5195 const bool sortSimplices = isConcaveFace;
5196 for ( size_t i = 0; i < smoothNodes.size(); ++i )
5198 const SMDS_MeshNode* n = smoothNodes[i];
5199 nodesToSmooth[ i ]._node = n;
5200 // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices
5201 getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, sortSimplices );
5202 // fix up incorrect uv of nodes on the FACE
5203 helper.GetNodeUV( F, n, 0, &isOkUV);
5207 //if ( nodesToSmooth.empty() ) continue;
5209 // Find EDGE's to shrink and set simpices to LayerEdge's
5210 set< _Shrinker1D* > eShri1D;
5212 for ( size_t i = 0; i < lEdges.size(); ++i )
5214 _LayerEdge* edge = lEdges[i];
5215 if ( edge->_sWOL.ShapeType() == TopAbs_EDGE )
5217 TGeomID edgeIndex = getMeshDS()->ShapeToIndex( edge->_sWOL );
5218 _Shrinker1D& srinker = e2shrMap[ edgeIndex ];
5219 eShri1D.insert( & srinker );
5220 srinker.AddEdge( edge, helper );
5221 VISCOUS_3D::ToClearSubWithMain( _mesh->GetSubMesh( edge->_sWOL ), data._solid );
5222 // restore params of nodes on EGDE if the EDGE has been already
5223 // srinked while srinking another FACE
5224 srinker.RestoreParams();
5226 getSimplices( /*tgtNode=*/edge->_nodes.back(), edge->_simplices, ignoreShapes );
5230 bool toFixTria = false; // to improve quality of trias by diagonal swap
5231 if ( isConcaveFace )
5233 const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles();
5234 if ( hasTria != hasQuad ) {
5235 toFixTria = hasTria;
5238 set<int> nbNodesSet;
5239 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
5240 while ( fIt->more() && nbNodesSet.size() < 2 )
5241 nbNodesSet.insert( fIt->next()->NbCornerNodes() );
5242 toFixTria = ( *nbNodesSet.begin() == 3 );
5246 // ==================
5247 // Perform shrinking
5248 // ==================
5250 bool shrinked = true;
5251 int badNb, shriStep=0, smooStep=0;
5252 _SmoothNode::SmoothType smoothType
5253 = isConcaveFace ? _SmoothNode::ANGULAR : _SmoothNode::LAPLACIAN;
5257 // Move boundary nodes (actually just set new UV)
5258 // -----------------------------------------------
5259 dumpFunction(SMESH_Comment("moveBoundaryOnF")<<f2sd->first<<"_st"<<shriStep ); // debug
5261 for ( size_t i = 0; i < lEdges.size(); ++i )
5263 shrinked |= lEdges[i]->SetNewLength2d( surface,F,helper );
5267 // Move nodes on EDGE's
5268 // (XYZ is set as soon as a needed length reached in SetNewLength2d())
5269 set< _Shrinker1D* >::iterator shr = eShri1D.begin();
5270 for ( ; shr != eShri1D.end(); ++shr )
5271 (*shr)->Compute( /*set3D=*/false, helper );
5274 // -----------------
5275 int nbNoImpSteps = 0;
5278 while (( nbNoImpSteps < 5 && badNb > 0) && moved)
5280 dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
5282 int oldBadNb = badNb;
5285 for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5287 moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
5288 smoothType, /*set3D=*/isConcaveFace);
5290 if ( badNb < oldBadNb )
5298 return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first );
5299 if ( shriStep > 200 )
5300 return error(SMESH_Comment("Infinite loop at shrinking 2D mesh on face ") << f2sd->first );
5302 // Fix narrow triangles by swapping diagonals
5303 // ---------------------------------------
5306 set<const SMDS_MeshNode*> usedNodes;
5307 fixBadFaces( F, helper, /*is2D=*/true, shriStep, & usedNodes); // swap diagonals
5309 // update working data
5310 set<const SMDS_MeshNode*>::iterator n;
5311 for ( size_t i = 0; i < nodesToSmooth.size() && !usedNodes.empty(); ++i )
5313 n = usedNodes.find( nodesToSmooth[ i ]._node );
5314 if ( n != usedNodes.end())
5316 getSimplices( nodesToSmooth[ i ]._node,
5317 nodesToSmooth[ i ]._simplices,
5319 /*sortSimplices=*/ smoothType == _SmoothNode::ANGULAR );
5320 usedNodes.erase( n );
5323 for ( size_t i = 0; i < lEdges.size() && !usedNodes.empty(); ++i )
5325 n = usedNodes.find( /*tgtNode=*/ lEdges[i]->_nodes.back() );
5326 if ( n != usedNodes.end())
5328 getSimplices( lEdges[i]->_nodes.back(),
5329 lEdges[i]->_simplices,
5331 usedNodes.erase( n );
5335 // TODO: check effect of this additional smooth
5336 // additional laplacian smooth to increase allowed shrink step
5337 // for ( int st = 1; st; --st )
5339 // dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
5340 // for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5342 // nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
5343 // _SmoothNode::LAPLACIAN,/*set3D=*/false);
5346 } // while ( shrinked )
5348 // No wrongly shaped faces remain; final smooth. Set node XYZ.
5349 bool isStructuredFixed = false;
5350 if ( SMESH_2D_Algo* algo = dynamic_cast<SMESH_2D_Algo*>( sm->GetAlgo() ))
5351 isStructuredFixed = algo->FixInternalNodes( *data._proxyMesh, F );
5352 if ( !isStructuredFixed )
5354 if ( isConcaveFace ) // fix narrow faces by swapping diagonals
5355 fixBadFaces( F, helper, /*is2D=*/false, ++shriStep );
5357 for ( int st = 3; st; --st )
5360 case 1: smoothType = _SmoothNode::LAPLACIAN; break;
5361 case 2: smoothType = _SmoothNode::LAPLACIAN; break;
5362 case 3: smoothType = _SmoothNode::ANGULAR; break;
5364 dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
5365 for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5367 nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
5368 smoothType,/*set3D=*/st==1 );
5373 // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh
5374 VISCOUS_3D::ToClearSubWithMain( sm, data._solid );
5376 if ( !getMeshDS()->IsEmbeddedMode() )
5377 // Log node movement
5378 for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5380 SMESH_TNodeXYZ p ( nodesToSmooth[i]._node );
5381 getMeshDS()->MoveNode( nodesToSmooth[i]._node, p.X(), p.Y(), p.Z() );
5384 } // loop on FACES to srink mesh on
5387 // Replace source nodes by target nodes in shrinked mesh edges
5389 map< int, _Shrinker1D >::iterator e2shr = e2shrMap.begin();
5390 for ( ; e2shr != e2shrMap.end(); ++e2shr )
5391 e2shr->second.SwapSrcTgtNodes( getMeshDS() );
5396 //================================================================================
5398 * \brief Computes 2d shrink direction and finds nodes limiting shrinking
5400 //================================================================================
5402 bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge,
5403 const TopoDS_Face& F,
5404 SMESH_MesherHelper& helper,
5405 const SMESHDS_SubMesh* faceSubMesh)
5407 const SMDS_MeshNode* srcNode = edge._nodes[0];
5408 const SMDS_MeshNode* tgtNode = edge._nodes.back();
5410 if ( edge._sWOL.ShapeType() == TopAbs_FACE )
5412 gp_XY srcUV( edge._pos[0].X(), edge._pos[0].Y() );//helper.GetNodeUV( F, srcNode );
5413 gp_XY tgtUV = edge.LastUV( F ); //helper.GetNodeUV( F, tgtNode );
5414 gp_Vec2d uvDir( srcUV, tgtUV );
5415 double uvLen = uvDir.Magnitude();
5417 edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0 );
5420 edge._pos.resize(1);
5421 edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 );
5423 // set UV of source node to target node
5424 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
5425 pos->SetUParameter( srcUV.X() );
5426 pos->SetVParameter( srcUV.Y() );
5428 else // _sWOL is TopAbs_EDGE
5430 const TopoDS_Edge& E = TopoDS::Edge( edge._sWOL );
5431 SMESHDS_SubMesh* edgeSM = getMeshDS()->MeshElements( E );
5432 if ( !edgeSM || edgeSM->NbElements() == 0 )
5433 return error(SMESH_Comment("Not meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
5435 const SMDS_MeshNode* n2 = 0;
5436 SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
5437 while ( eIt->more() && !n2 )
5439 const SMDS_MeshElement* e = eIt->next();
5440 if ( !edgeSM->Contains(e)) continue;
5441 n2 = e->GetNode( 0 );
5442 if ( n2 == srcNode ) n2 = e->GetNode( 1 );
5445 return error(SMESH_Comment("Wrongly meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
5447 double uSrc = helper.GetNodeU( E, srcNode, n2 );
5448 double uTgt = helper.GetNodeU( E, tgtNode, srcNode );
5449 double u2 = helper.GetNodeU( E, n2, srcNode );
5453 if ( fabs( uSrc-uTgt ) < 0.99 * fabs( uSrc-u2 ))
5455 // tgtNode is located so that it does not make faces with wrong orientation
5458 edge._pos.resize(1);
5459 edge._pos[0].SetCoord( U_TGT, uTgt );
5460 edge._pos[0].SetCoord( U_SRC, uSrc );
5461 edge._pos[0].SetCoord( LEN_TGT, fabs( uSrc-uTgt ));
5463 edge._simplices.resize( 1 );
5464 edge._simplices[0]._nPrev = n2;
5466 // set U of source node to the target node
5467 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
5468 pos->SetUParameter( uSrc );
5473 //================================================================================
5475 * \brief Restore position of a sole node of a _LayerEdge based on _noShrinkShapes
5477 //================================================================================
5479 void _ViscousBuilder::restoreNoShrink( _LayerEdge& edge ) const
5481 if ( edge._nodes.size() == 1 )
5486 const SMDS_MeshNode* srcNode = edge._nodes[0];
5487 TopoDS_Shape S = SMESH_MesherHelper::GetSubShapeByNode( srcNode, getMeshDS() );
5488 if ( S.IsNull() ) return;
5492 switch ( S.ShapeType() )
5497 TopLoc_Location loc;
5498 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( S ), loc, f, l );
5499 if ( curve.IsNull() ) return;
5500 SMDS_EdgePosition* ePos = static_cast<SMDS_EdgePosition*>( srcNode->GetPosition() );
5501 p = curve->Value( ePos->GetUParameter() );
5506 p = BRep_Tool::Pnt( TopoDS::Vertex( S ));
5511 getMeshDS()->MoveNode( srcNode, p.X(), p.Y(), p.Z() );
5512 dumpMove( srcNode );
5516 //================================================================================
5518 * \brief Try to fix triangles with high aspect ratio by swaping diagonals
5520 //================================================================================
5522 void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F,
5523 SMESH_MesherHelper& helper,
5526 set<const SMDS_MeshNode*> * involvedNodes)
5528 SMESH::Controls::AspectRatio qualifier;
5529 SMESH::Controls::TSequenceOfXYZ points(3), points1(3), points2(3);
5530 const double maxAspectRatio = is2D ? 4. : 2;
5531 _NodeCoordHelper xyz( F, helper, is2D );
5533 // find bad triangles
5535 vector< const SMDS_MeshElement* > badTrias;
5536 vector< double > badAspects;
5537 SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( F );
5538 SMDS_ElemIteratorPtr fIt = sm->GetElements();
5539 while ( fIt->more() )
5541 const SMDS_MeshElement * f = fIt->next();
5542 if ( f->NbCornerNodes() != 3 ) continue;
5543 for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = xyz( f->GetNode(iP));
5544 double aspect = qualifier.GetValue( points );
5545 if ( aspect > maxAspectRatio )
5547 badTrias.push_back( f );
5548 badAspects.push_back( aspect );
5553 dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID());
5554 SMDS_ElemIteratorPtr fIt = sm->GetElements();
5555 while ( fIt->more() )
5557 const SMDS_MeshElement * f = fIt->next();
5558 if ( f->NbCornerNodes() == 3 )
5559 dumpChangeNodes( f );
5563 if ( badTrias.empty() )
5566 // find couples of faces to swap diagonal
5568 typedef pair < const SMDS_MeshElement* , const SMDS_MeshElement* > T2Trias;
5569 vector< T2Trias > triaCouples;
5571 TIDSortedElemSet involvedFaces, emptySet;
5572 for ( size_t iTia = 0; iTia < badTrias.size(); ++iTia )
5575 double aspRatio [3];
5578 if ( !involvedFaces.insert( badTrias[iTia] ).second )
5580 for ( int iP = 0; iP < 3; ++iP )
5581 points(iP+1) = xyz( badTrias[iTia]->GetNode(iP));
5583 // find triangles adjacent to badTrias[iTia] with better aspect ratio after diag-swaping
5584 int bestCouple = -1;
5585 for ( int iSide = 0; iSide < 3; ++iSide )
5587 const SMDS_MeshNode* n1 = badTrias[iTia]->GetNode( iSide );
5588 const SMDS_MeshNode* n2 = badTrias[iTia]->GetNode(( iSide+1 ) % 3 );
5589 trias [iSide].first = badTrias[iTia];
5590 trias [iSide].second = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, involvedFaces,
5592 if (( ! trias[iSide].second ) ||
5593 ( trias[iSide].second->NbCornerNodes() != 3 ) ||
5594 ( ! sm->Contains( trias[iSide].second )))
5597 // aspect ratio of an adjacent tria
5598 for ( int iP = 0; iP < 3; ++iP )
5599 points2(iP+1) = xyz( trias[iSide].second->GetNode(iP));
5600 double aspectInit = qualifier.GetValue( points2 );
5602 // arrange nodes as after diag-swaping
5603 if ( helper.WrapIndex( i1+1, 3 ) == i2 )
5604 i3 = helper.WrapIndex( i1-1, 3 );
5606 i3 = helper.WrapIndex( i1+1, 3 );
5608 points1( 1+ iSide ) = points2( 1+ i3 );
5609 points2( 1+ i2 ) = points1( 1+ ( iSide+2 ) % 3 );
5611 // aspect ratio after diag-swaping
5612 aspRatio[ iSide ] = qualifier.GetValue( points1 ) + qualifier.GetValue( points2 );
5613 if ( aspRatio[ iSide ] > aspectInit + badAspects[ iTia ] )
5616 // prevent inversion of a triangle
5617 gp_Vec norm1 = gp_Vec( points1(1), points1(3) ) ^ gp_Vec( points1(1), points1(2) );
5618 gp_Vec norm2 = gp_Vec( points2(1), points2(3) ) ^ gp_Vec( points2(1), points2(2) );
5619 if ( norm1 * norm2 < 0. && norm1.Angle( norm2 ) > 70./180.*M_PI )
5622 if ( bestCouple < 0 || aspRatio[ bestCouple ] > aspRatio[ iSide ] )
5626 if ( bestCouple >= 0 )
5628 triaCouples.push_back( trias[bestCouple] );
5629 involvedFaces.insert ( trias[bestCouple].second );
5633 involvedFaces.erase( badTrias[iTia] );
5636 if ( triaCouples.empty() )
5641 SMESH_MeshEditor editor( helper.GetMesh() );
5642 dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
5643 for ( size_t i = 0; i < triaCouples.size(); ++i )
5645 dumpChangeNodes( triaCouples[i].first );
5646 dumpChangeNodes( triaCouples[i].second );
5647 editor.InverseDiag( triaCouples[i].first, triaCouples[i].second );
5650 if ( involvedNodes )
5651 for ( size_t i = 0; i < triaCouples.size(); ++i )
5653 involvedNodes->insert( triaCouples[i].first->begin_nodes(),
5654 triaCouples[i].first->end_nodes() );
5655 involvedNodes->insert( triaCouples[i].second->begin_nodes(),
5656 triaCouples[i].second->end_nodes() );
5659 // just for debug dump resulting triangles
5660 dumpFunction(SMESH_Comment("swapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
5661 for ( size_t i = 0; i < triaCouples.size(); ++i )
5663 dumpChangeNodes( triaCouples[i].first );
5664 dumpChangeNodes( triaCouples[i].second );
5668 //================================================================================
5670 * \brief Move target node to it's final position on the FACE during shrinking
5672 //================================================================================
5674 bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
5675 const TopoDS_Face& F,
5676 SMESH_MesherHelper& helper )
5679 return false; // already at the target position
5681 SMDS_MeshNode* tgtNode = const_cast< SMDS_MeshNode*& >( _nodes.back() );
5683 if ( _sWOL.ShapeType() == TopAbs_FACE )
5685 gp_XY curUV = helper.GetNodeUV( F, tgtNode );
5686 gp_Pnt2d tgtUV( _pos[0].X(), _pos[0].Y() );
5687 gp_Vec2d uvDir( _normal.X(), _normal.Y() );
5688 const double uvLen = tgtUV.Distance( curUV );
5689 const double kSafe = Max( 0.5, 1. - 0.1 * _simplices.size() );
5691 // Select shrinking step such that not to make faces with wrong orientation.
5692 double stepSize = 1e100;
5693 for ( size_t i = 0; i < _simplices.size(); ++i )
5695 // find intersection of 2 lines: curUV-tgtUV and that connecting simplex nodes
5696 gp_XY uvN1 = helper.GetNodeUV( F, _simplices[i]._nPrev );
5697 gp_XY uvN2 = helper.GetNodeUV( F, _simplices[i]._nNext );
5698 gp_XY dirN = uvN2 - uvN1;
5699 double det = uvDir.Crossed( dirN );
5700 if ( Abs( det ) < std::numeric_limits<double>::min() ) continue;
5701 gp_XY dirN2Cur = curUV - uvN1;
5702 double step = dirN.Crossed( dirN2Cur ) / det;
5704 stepSize = Min( step, stepSize );
5707 if ( uvLen <= stepSize )
5712 else if ( stepSize > 0 )
5714 newUV = curUV + uvDir.XY() * stepSize * kSafe;
5720 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
5721 pos->SetUParameter( newUV.X() );
5722 pos->SetVParameter( newUV.Y() );
5725 gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
5726 tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
5727 dumpMove( tgtNode );
5730 else // _sWOL is TopAbs_EDGE
5732 const TopoDS_Edge& E = TopoDS::Edge( _sWOL );
5733 const SMDS_MeshNode* n2 = _simplices[0]._nPrev;
5734 SMDS_EdgePosition* tgtPos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
5736 const double u2 = helper.GetNodeU( E, n2, tgtNode );
5737 const double uSrc = _pos[0].Coord( U_SRC );
5738 const double lenTgt = _pos[0].Coord( LEN_TGT );
5740 double newU = _pos[0].Coord( U_TGT );
5741 if ( lenTgt < 0.99 * fabs( uSrc-u2 )) // n2 got out of src-tgt range
5747 newU = 0.1 * tgtPos->GetUParameter() + 0.9 * u2;
5749 tgtPos->SetUParameter( newU );
5751 gp_XY newUV = helper.GetNodeUV( F, tgtNode, _nodes[0]);
5752 gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
5753 tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
5754 dumpMove( tgtNode );
5760 //================================================================================
5762 * \brief Perform smooth on the FACE
5763 * \retval bool - true if the node has been moved
5765 //================================================================================
5767 bool _SmoothNode::Smooth(int& badNb,
5768 Handle(Geom_Surface)& surface,
5769 SMESH_MesherHelper& helper,
5770 const double refSign,
5774 const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() );
5776 // get uv of surrounding nodes
5777 vector<gp_XY> uv( _simplices.size() );
5778 for ( size_t i = 0; i < _simplices.size(); ++i )
5779 uv[i] = helper.GetNodeUV( face, _simplices[i]._nPrev, _node );
5781 // compute new UV for the node
5783 if ( how == TFI && _simplices.size() == 4 )
5786 for ( size_t i = 0; i < _simplices.size(); ++i )
5787 if ( _simplices[i]._nOpp )
5788 corners[i] = helper.GetNodeUV( face, _simplices[i]._nOpp, _node );
5790 throw SALOME_Exception(LOCALIZED("TFI smoothing: _Simplex::_nOpp not set!"));
5792 newPos = helper.calcTFI ( 0.5, 0.5,
5793 corners[0], corners[1], corners[2], corners[3],
5794 uv[1], uv[2], uv[3], uv[0] );
5796 else if ( how == ANGULAR )
5798 newPos = computeAngularPos( uv, helper.GetNodeUV( face, _node ), refSign );
5800 else if ( how == CENTROIDAL && _simplices.size() > 3 )
5802 // average centers of diagonals wieghted with their reciprocal lengths
5803 if ( _simplices.size() == 4 )
5805 double w1 = 1. / ( uv[2]-uv[0] ).SquareModulus();
5806 double w2 = 1. / ( uv[3]-uv[1] ).SquareModulus();
5807 newPos = ( w1 * ( uv[2]+uv[0] ) + w2 * ( uv[3]+uv[1] )) / ( w1+w2 ) / 2;
5811 double sumWeight = 0;
5812 int nb = _simplices.size() == 4 ? 2 : _simplices.size();
5813 for ( int i = 0; i < nb; ++i )
5816 int iTo = i + _simplices.size() - 1;
5817 for ( int j = iFrom; j < iTo; ++j )
5819 int i2 = SMESH_MesherHelper::WrapIndex( j, _simplices.size() );
5820 double w = 1. / ( uv[i]-uv[i2] ).SquareModulus();
5822 newPos += w * ( uv[i]+uv[i2] );
5825 newPos /= 2 * sumWeight; // 2 is to get a middle between uv's
5831 for ( size_t i = 0; i < _simplices.size(); ++i )
5833 newPos /= _simplices.size();
5836 // count quality metrics (orientation) of triangles around the node
5838 gp_XY tgtUV = helper.GetNodeUV( face, _node );
5839 for ( size_t i = 0; i < _simplices.size(); ++i )
5840 nbOkBefore += _simplices[i].IsForward( tgtUV, _node, face, helper, refSign );
5843 for ( size_t i = 0; i < _simplices.size(); ++i )
5844 nbOkAfter += _simplices[i].IsForward( newPos, _node, face, helper, refSign );
5846 if ( nbOkAfter < nbOkBefore )
5848 badNb += _simplices.size() - nbOkBefore;
5852 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( _node->GetPosition() );
5853 pos->SetUParameter( newPos.X() );
5854 pos->SetVParameter( newPos.Y() );
5861 gp_Pnt p = surface->Value( newPos.X(), newPos.Y() );
5862 const_cast< SMDS_MeshNode* >( _node )->setXYZ( p.X(), p.Y(), p.Z() );
5866 badNb += _simplices.size() - nbOkAfter;
5867 return ( (tgtUV-newPos).SquareModulus() > 1e-10 );
5870 //================================================================================
5872 * \brief Computes new UV using angle based smoothing technic
5874 //================================================================================
5876 gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
5877 const gp_XY& uvToFix,
5878 const double refSign)
5880 uv.push_back( uv.front() );
5882 vector< gp_XY > edgeDir ( uv.size() );
5883 vector< double > edgeSize( uv.size() );
5884 for ( size_t i = 1; i < edgeDir.size(); ++i )
5886 edgeDir [i-1] = uv[i] - uv[i-1];
5887 edgeSize[i-1] = edgeDir[i-1].Modulus();
5888 if ( edgeSize[i-1] < numeric_limits<double>::min() )
5889 edgeDir[i-1].SetX( 100 );
5891 edgeDir[i-1] /= edgeSize[i-1] * refSign;
5893 edgeDir.back() = edgeDir.front();
5894 edgeSize.back() = edgeSize.front();
5899 for ( size_t i = 1; i < edgeDir.size(); ++i )
5901 if ( edgeDir[i-1].X() > 1. ) continue;
5903 while ( edgeDir[i].X() > 1. && ++i < edgeDir.size() );
5904 if ( i == edgeDir.size() ) break;
5906 gp_XY norm1( -edgeDir[i1].Y(), edgeDir[i1].X() );
5907 gp_XY norm2( -edgeDir[i].Y(), edgeDir[i].X() );
5908 gp_XY bisec = norm1 + norm2;
5909 double bisecSize = bisec.Modulus();
5910 if ( bisecSize < numeric_limits<double>::min() )
5912 bisec = -edgeDir[i1] + edgeDir[i];
5913 bisecSize = bisec.Modulus();
5917 gp_XY dirToN = uvToFix - p;
5918 double distToN = dirToN.Modulus();
5919 if ( bisec * dirToN < 0 )
5922 newPos += ( p + bisec * distToN ) * ( edgeSize[i1] + edgeSize[i] );
5924 sumSize += edgeSize[i1] + edgeSize[i];
5926 newPos /= /*nbEdges * */sumSize;
5930 //================================================================================
5932 * \brief Delete _SolidData
5934 //================================================================================
5936 _SolidData::~_SolidData()
5938 for ( size_t i = 0; i < _edges.size(); ++i )
5940 if ( _edges[i] && _edges[i]->_2neibors )
5941 delete _edges[i]->_2neibors;
5946 //================================================================================
5948 * \brief Add a _LayerEdge inflated along the EDGE
5950 //================================================================================
5952 void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper )
5955 if ( _nodes.empty() )
5957 _edges[0] = _edges[1] = 0;
5961 if ( e == _edges[0] || e == _edges[1] )
5963 if ( e->_sWOL.IsNull() || e->_sWOL.ShapeType() != TopAbs_EDGE )
5964 throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
5965 if ( _edges[0] && _edges[0]->_sWOL != e->_sWOL )
5966 throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
5969 const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
5971 BRep_Tool::Range( E, f,l );
5972 double u = helper.GetNodeU( E, e->_nodes[0], e->_nodes.back());
5973 _edges[ u < 0.5*(f+l) ? 0 : 1 ] = e;
5977 const SMDS_MeshNode* tgtNode0 = _edges[0] ? _edges[0]->_nodes.back() : 0;
5978 const SMDS_MeshNode* tgtNode1 = _edges[1] ? _edges[1]->_nodes.back() : 0;
5980 if ( _nodes.empty() )
5982 SMESHDS_SubMesh * eSubMesh = helper.GetMeshDS()->MeshElements( E );
5983 if ( !eSubMesh || eSubMesh->NbNodes() < 1 )
5985 TopLoc_Location loc;
5986 Handle(Geom_Curve) C = BRep_Tool::Curve(E, loc, f,l);
5987 GeomAdaptor_Curve aCurve(C, f,l);
5988 const double totLen = GCPnts_AbscissaPoint::Length(aCurve, f, l);
5990 int nbExpectNodes = eSubMesh->NbNodes();
5991 _initU .reserve( nbExpectNodes );
5992 _normPar.reserve( nbExpectNodes );
5993 _nodes .reserve( nbExpectNodes );
5994 SMDS_NodeIteratorPtr nIt = eSubMesh->GetNodes();
5995 while ( nIt->more() )
5997 const SMDS_MeshNode* node = nIt->next();
5998 if ( node->NbInverseElements(SMDSAbs_Edge) == 0 ||
5999 node == tgtNode0 || node == tgtNode1 )
6000 continue; // refinement nodes
6001 _nodes.push_back( node );
6002 _initU.push_back( helper.GetNodeU( E, node ));
6003 double len = GCPnts_AbscissaPoint::Length(aCurve, f, _initU.back());
6004 _normPar.push_back( len / totLen );
6009 // remove target node of the _LayerEdge from _nodes
6011 for ( size_t i = 0; i < _nodes.size(); ++i )
6012 if ( !_nodes[i] || _nodes[i] == tgtNode0 || _nodes[i] == tgtNode1 )
6013 _nodes[i] = 0, nbFound++;
6014 if ( nbFound == _nodes.size() )
6019 //================================================================================
6021 * \brief Move nodes on EDGE from ends where _LayerEdge's are inflated
6023 //================================================================================
6025 void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
6027 if ( _done || _nodes.empty())
6029 const _LayerEdge* e = _edges[0];
6030 if ( !e ) e = _edges[1];
6033 _done = (( !_edges[0] || _edges[0]->_pos.empty() ) &&
6034 ( !_edges[1] || _edges[1]->_pos.empty() ));
6036 const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
6038 if ( set3D || _done )
6040 Handle(Geom_Curve) C = BRep_Tool::Curve(E, f,l);
6041 GeomAdaptor_Curve aCurve(C, f,l);
6044 f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
6046 l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
6047 double totLen = GCPnts_AbscissaPoint::Length( aCurve, f, l );
6049 for ( size_t i = 0; i < _nodes.size(); ++i )
6051 if ( !_nodes[i] ) continue;
6052 double len = totLen * _normPar[i];
6053 GCPnts_AbscissaPoint discret( aCurve, len, f );
6054 if ( !discret.IsDone() )
6055 return throw SALOME_Exception(LOCALIZED("GCPnts_AbscissaPoint failed"));
6056 double u = discret.Parameter();
6057 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
6058 pos->SetUParameter( u );
6059 gp_Pnt p = C->Value( u );
6060 const_cast< SMDS_MeshNode*>( _nodes[i] )->setXYZ( p.X(), p.Y(), p.Z() );
6065 BRep_Tool::Range( E, f,l );
6067 f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
6069 l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
6071 for ( size_t i = 0; i < _nodes.size(); ++i )
6073 if ( !_nodes[i] ) continue;
6074 double u = f * ( 1-_normPar[i] ) + l * _normPar[i];
6075 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
6076 pos->SetUParameter( u );
6081 //================================================================================
6083 * \brief Restore initial parameters of nodes on EDGE
6085 //================================================================================
6087 void _Shrinker1D::RestoreParams()
6090 for ( size_t i = 0; i < _nodes.size(); ++i )
6092 if ( !_nodes[i] ) continue;
6093 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
6094 pos->SetUParameter( _initU[i] );
6099 //================================================================================
6101 * \brief Replace source nodes by target nodes in shrinked mesh edges
6103 //================================================================================
6105 void _Shrinker1D::SwapSrcTgtNodes( SMESHDS_Mesh* mesh )
6107 const SMDS_MeshNode* nodes[3];
6108 for ( int i = 0; i < 2; ++i )
6110 if ( !_edges[i] ) continue;
6112 SMESHDS_SubMesh * eSubMesh = mesh->MeshElements( _edges[i]->_sWOL );
6113 if ( !eSubMesh ) return;
6114 const SMDS_MeshNode* srcNode = _edges[i]->_nodes[0];
6115 const SMDS_MeshNode* tgtNode = _edges[i]->_nodes.back();
6116 SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
6117 while ( eIt->more() )
6119 const SMDS_MeshElement* e = eIt->next();
6120 if ( !eSubMesh->Contains( e ))
6122 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
6123 for ( int iN = 0; iN < e->NbNodes(); ++iN )
6125 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
6126 nodes[iN] = ( n == srcNode ? tgtNode : n );
6128 mesh->ChangeElementNodes( e, nodes, e->NbNodes() );
6133 //================================================================================
6135 * \brief Creates 2D and 1D elements on boundaries of new prisms
6137 //================================================================================
6139 bool _ViscousBuilder::addBoundaryElements()
6141 SMESH_MesherHelper helper( *_mesh );
6143 for ( size_t i = 0; i < _sdVec.size(); ++i )
6145 _SolidData& data = _sdVec[i];
6146 TopTools_IndexedMapOfShape geomEdges;
6147 TopExp::MapShapes( data._solid, TopAbs_EDGE, geomEdges );
6148 for ( int iE = 1; iE <= geomEdges.Extent(); ++iE )
6150 const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE));
6151 if ( data._noShrinkShapes.count( getMeshDS()->ShapeToIndex( E )))
6154 // Get _LayerEdge's based on E
6156 map< double, const SMDS_MeshNode* > u2nodes;
6157 if ( !SMESH_Algo::GetSortedNodesOnEdge( getMeshDS(), E, /*ignoreMedium=*/false, u2nodes))
6160 vector< _LayerEdge* > ledges; ledges.reserve( u2nodes.size() );
6161 TNode2Edge & n2eMap = data._n2eMap;
6162 map< double, const SMDS_MeshNode* >::iterator u2n = u2nodes.begin();
6164 //check if 2D elements are needed on E
6165 TNode2Edge::iterator n2e = n2eMap.find( u2n->second );
6166 if ( n2e == n2eMap.end() ) continue; // no layers on vertex
6167 ledges.push_back( n2e->second );
6169 if (( n2e = n2eMap.find( u2n->second )) == n2eMap.end() )
6170 continue; // no layers on E
6171 ledges.push_back( n2eMap[ u2n->second ]);
6173 const SMDS_MeshNode* tgtN0 = ledges[0]->_nodes.back();
6174 const SMDS_MeshNode* tgtN1 = ledges[1]->_nodes.back();
6175 int nbSharedPyram = 0;
6176 SMDS_ElemIteratorPtr vIt = tgtN0->GetInverseElementIterator(SMDSAbs_Volume);
6177 while ( vIt->more() )
6179 const SMDS_MeshElement* v = vIt->next();
6180 nbSharedPyram += int( v->GetNodeIndex( tgtN1 ) >= 0 );
6182 if ( nbSharedPyram > 1 )
6183 continue; // not free border of the pyramid
6185 if ( getMeshDS()->FindFace( ledges[0]->_nodes[0], ledges[0]->_nodes[1],
6186 ledges[1]->_nodes[0], ledges[1]->_nodes[1]))
6187 continue; // faces already created
6189 for ( ++u2n; u2n != u2nodes.end(); ++u2n )
6190 ledges.push_back( n2eMap[ u2n->second ]);
6192 // Find out orientation and type of face to create
6194 bool reverse = false, isOnFace;
6196 map< TGeomID, TopoDS_Shape >::iterator e2f =
6197 data._shrinkShape2Shape.find( getMeshDS()->ShapeToIndex( E ));
6199 if (( isOnFace = ( e2f != data._shrinkShape2Shape.end() )))
6201 F = e2f->second.Oriented( TopAbs_FORWARD );
6202 reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED );
6203 if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED )
6204 reverse = !reverse, F.Reverse();
6205 if ( helper.IsReversedSubMesh( TopoDS::Face(F) ))
6210 // find FACE with layers sharing E
6211 PShapeIteratorPtr fIt = helper.GetAncestors( E, *_mesh, TopAbs_FACE );
6212 while ( fIt->more() && F.IsNull() )
6214 const TopoDS_Shape* pF = fIt->next();
6215 if ( helper.IsSubShape( *pF, data._solid) &&
6216 !data._ignoreFaceIds.count( e2f->first ))
6220 // Find the sub-mesh to add new faces
6221 SMESHDS_SubMesh* sm = 0;
6223 sm = getMeshDS()->MeshElements( F );
6225 sm = data._proxyMesh->getFaceSubM( TopoDS::Face(F), /*create=*/true );
6227 return error("error in addBoundaryElements()", data._index);
6230 const int dj1 = reverse ? 0 : 1;
6231 const int dj2 = reverse ? 1 : 0;
6232 for ( size_t j = 1; j < ledges.size(); ++j )
6234 vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes;
6235 vector< const SMDS_MeshNode*>& nn2 = ledges[j-dj2]->_nodes;
6236 if ( nn1.size() == nn2.size() )
6239 for ( size_t z = 1; z < nn1.size(); ++z )
6240 sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
6242 for ( size_t z = 1; z < nn1.size(); ++z )
6243 sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
6245 else if ( nn1.size() == 1 )
6248 for ( size_t z = 1; z < nn2.size(); ++z )
6249 sm->AddElement( getMeshDS()->AddFace( nn1[0], nn2[z-1], nn2[z] ));
6251 for ( size_t z = 1; z < nn2.size(); ++z )
6252 sm->AddElement( new SMDS_FaceOfNodes( nn1[0], nn2[z-1], nn2[z] ));
6257 for ( size_t z = 1; z < nn1.size(); ++z )
6258 sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[0], nn1[z] ));
6260 for ( size_t z = 1; z < nn1.size(); ++z )
6261 sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[0], nn2[z] ));
6266 for ( int isFirst = 0; isFirst < 2; ++isFirst )
6268 _LayerEdge* edge = isFirst ? ledges.front() : ledges.back();
6269 if ( !edge->_sWOL.IsNull() && edge->_sWOL.ShapeType() == TopAbs_EDGE )
6271 vector< const SMDS_MeshNode*>& nn = edge->_nodes;
6272 if ( nn.size() < 2 || nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() )
6274 helper.SetSubShape( edge->_sWOL );
6275 helper.SetElementsOnShape( true );
6276 for ( size_t z = 1; z < nn.size(); ++z )
6277 helper.AddEdge( nn[z-1], nn[z] );