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 double u = helper.GetNodeU( E, atNode );
856 //--------------------------------------------------------------------------------
857 gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE,
858 const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok)
860 gp_XY uv = helper.GetNodeUV( F, node, 0, &ok );
861 Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
862 gp_Pnt p; gp_Vec du, dv, norm;
863 surface->D1( uv.X(),uv.Y(), p, du,dv );
867 Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l );
868 double u = helper.GetNodeU( fromE, node, 0, &ok );
870 TopAbs_Orientation o = helper.GetSubShapeOri( F.Oriented(TopAbs_FORWARD), fromE);
871 if ( o == TopAbs_REVERSED )
874 gp_Vec dir = norm ^ du;
876 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX &&
877 helper.IsClosedEdge( fromE ))
879 if ( fabs(u-f) < fabs(u-l)) c->D1( l, p, dv );
880 else c->D1( f, p, dv );
881 if ( o == TopAbs_REVERSED )
883 gp_Vec dir2 = norm ^ dv;
884 dir = dir.Normalized() + dir2.Normalized();
888 //--------------------------------------------------------------------------------
889 gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
890 const SMDS_MeshNode* node, SMESH_MesherHelper& helper,
891 bool& ok, double* cosin=0)
893 TopoDS_Face faceFrw = F;
894 faceFrw.Orientation( TopAbs_FORWARD );
895 double f,l; TopLoc_Location loc;
896 TopoDS_Edge edges[2]; // sharing a vertex
900 TopExp_Explorer exp( faceFrw, TopAbs_EDGE );
901 for ( ; exp.More() && nbEdges < 2; exp.Next() )
903 const TopoDS_Edge& e = TopoDS::Edge( exp.Current() );
904 if ( SMESH_Algo::isDegenerated( e )) continue;
905 TopExp::Vertices( e, VV[0], VV[1], /*CumOri=*/true );
906 if ( VV[1].IsSame( fromV )) {
907 nbEdges += edges[ 0 ].IsNull();
910 else if ( VV[0].IsSame( fromV )) {
911 nbEdges += edges[ 1 ].IsNull();
916 gp_XYZ dir(0,0,0), edgeDir[2];
919 // get dirs of edges going fromV
921 for ( size_t i = 0; i < nbEdges && ok; ++i )
923 edgeDir[i] = getEdgeDir( edges[i], fromV );
924 double size2 = edgeDir[i].SquareModulus();
925 if (( ok = size2 > numeric_limits<double>::min() ))
926 edgeDir[i] /= sqrt( size2 );
928 if ( !ok ) return dir;
930 // get angle between the 2 edges
932 double angle = helper.GetAngle( edges[0], edges[1], faceFrw, fromV, &faceNormal );
933 if ( Abs( angle ) < 5 * M_PI/180 )
935 dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] );
939 dir = edgeDir[0] + edgeDir[1];
944 double angle = gp_Vec( edgeDir[0] ).Angle( dir );
945 *cosin = Cos( angle );
948 else if ( nbEdges == 1 )
950 dir = getFaceDir( faceFrw, edges[ edges[0].IsNull() ], node, helper, ok );
951 if ( cosin ) *cosin = 1.;
960 //================================================================================
962 * \brief Returns true if a FACE is bound by a concave EDGE
964 //================================================================================
966 bool isConcave( const TopoDS_Face& F, SMESH_MesherHelper& helper )
968 // if ( helper.Count( F, TopAbs_WIRE, /*useMap=*/false) > 1 )
972 TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE );
973 for ( ; eExp.More(); eExp.Next() )
975 const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
976 if ( SMESH_Algo::isDegenerated( E )) continue;
977 // check if 2D curve is concave
978 BRepAdaptor_Curve2d curve( E, F );
979 const int nbIntervals = curve.NbIntervals( GeomAbs_C2 );
980 TColStd_Array1OfReal intervals(1, nbIntervals + 1 );
981 curve.Intervals( intervals, GeomAbs_C2 );
982 bool isConvex = true;
983 for ( int i = 1; i <= nbIntervals && isConvex; ++i )
985 double u1 = intervals( i );
986 double u2 = intervals( i+1 );
987 curve.D2( 0.5*( u1+u2 ), p, drv1, drv2 );
988 double cross = drv2 ^ drv1;
989 if ( E.Orientation() == TopAbs_REVERSED )
991 isConvex = ( cross > 0.1 ); //-1e-9 );
995 //cout << "Concave FACE " << helper.GetMeshDS()->ShapeToIndex( F ) << endl;
999 // check angles at VERTEXes
1001 TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *helper.GetMesh(), 0, error );
1002 for ( size_t iW = 0; iW < wires.size(); ++iW )
1004 const int nbEdges = wires[iW]->NbEdges();
1005 if ( nbEdges < 2 && SMESH_Algo::isDegenerated( wires[iW]->Edge(0)))
1007 for ( int iE1 = 0; iE1 < nbEdges; ++iE1 )
1009 if ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE1 ))) continue;
1010 int iE2 = ( iE1 + 1 ) % nbEdges;
1011 while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 )))
1012 iE2 = ( iE2 + 1 ) % nbEdges;
1013 double angle = helper.GetAngle( wires[iW]->Edge( iE1 ),
1014 wires[iW]->Edge( iE2 ), F,
1015 wires[iW]->FirstVertex( iE2 ));
1016 if ( angle < -5. * M_PI / 180. )
1022 //--------------------------------------------------------------------------------
1023 // DEBUG. Dump intermediate node positions into a python script
1024 // HOWTO use: run python commands written in a console to see
1025 // construction steps of viscous layers
1031 const char* fname = "/tmp/viscous.py";
1032 cout << "execfile('"<<fname<<"')"<<endl;
1033 py = new ofstream(fname);
1034 *py << "import SMESH" << endl
1035 << "from salome.smesh import smeshBuilder" << endl
1036 << "smesh = smeshBuilder.New(salome.myStudy)" << endl
1037 << "meshSO = smesh.GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
1038 << "mesh = smesh.Mesh( meshSO.GetObject() )"<<endl;
1043 *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Viscous Prisms',"
1044 "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA))"<<endl;
1045 *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Neg Volumes',"
1046 "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_Volume3D,'<',0))"<<endl;
1050 ~PyDump() { Finish(); cout << "NB FUNCTIONS: " << theNbFunc << endl; }
1052 #define dumpFunction(f) { _dumpFunction(f, __LINE__);}
1053 #define dumpMove(n) { _dumpMove(n, __LINE__);}
1054 #define dumpCmd(txt) { _dumpCmd(txt, __LINE__);}
1055 void _dumpFunction(const string& fun, int ln)
1056 { if (py) *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl; ++theNbFunc; }
1057 void _dumpMove(const SMDS_MeshNode* n, int ln)
1058 { if (py) *py<< " mesh.MoveNode( "<<n->GetID()<< ", "<< n->X()
1059 << ", "<<n->Y()<<", "<< n->Z()<< ")\t\t # "<< ln <<endl; }
1060 void _dumpCmd(const string& txt, int ln)
1061 { if (py) *py<< " "<<txt<<" # "<< ln <<endl; }
1062 void dumpFunctionEnd()
1063 { if (py) *py<< " return"<< endl; }
1064 void dumpChangeNodes( const SMDS_MeshElement* f )
1065 { if (py) { *py<< " mesh.ChangeElemNodes( " << f->GetID()<<", [";
1066 for ( int i=1; i < f->NbNodes(); ++i ) *py << f->GetNode(i-1)->GetID()<<", ";
1067 *py << f->GetNode( f->NbNodes()-1 )->GetID() << " ])"<< endl; }}
1068 #define debugMsg( txt ) { cout << txt << " (line: " << __LINE__ << ")" << endl; }
1070 struct PyDump { void Finish() {} };
1071 #define dumpFunction(f) f
1073 #define dumpCmd(txt)
1074 #define dumpFunctionEnd()
1075 #define dumpChangeNodes(f)
1076 #define debugMsg( txt ) {}
1080 using namespace VISCOUS_3D;
1082 //================================================================================
1084 * \brief Constructor of _ViscousBuilder
1086 //================================================================================
1088 _ViscousBuilder::_ViscousBuilder()
1090 _error = SMESH_ComputeError::New(COMPERR_OK);
1094 //================================================================================
1096 * \brief Stores error description and returns false
1098 //================================================================================
1100 bool _ViscousBuilder::error(const string& text, int solidId )
1102 const string prefix = string("Viscous layers builder: ");
1103 _error->myName = COMPERR_ALGO_FAILED;
1104 _error->myComment = prefix + text;
1107 SMESH_subMesh* sm = _mesh->GetSubMeshContaining( solidId );
1108 if ( !sm && !_sdVec.empty() )
1109 sm = _mesh->GetSubMeshContaining( solidId = _sdVec[0]._index );
1110 if ( sm && sm->GetSubShape().ShapeType() == TopAbs_SOLID )
1112 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1113 if ( smError && smError->myAlgo )
1114 _error->myAlgo = smError->myAlgo;
1116 sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1118 // set KO to all solids
1119 for ( size_t i = 0; i < _sdVec.size(); ++i )
1121 if ( _sdVec[i]._index == solidId )
1123 sm = _mesh->GetSubMesh( _sdVec[i]._solid );
1124 if ( !sm->IsEmpty() )
1126 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1127 if ( !smError || smError->IsOK() )
1129 smError = SMESH_ComputeError::New( COMPERR_ALGO_FAILED, prefix + "failed");
1130 sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1134 makeGroupOfLE(); // debug
1139 //================================================================================
1141 * \brief At study restoration, restore event listeners used to clear an inferior
1142 * dim sub-mesh modified by viscous layers
1144 //================================================================================
1146 void _ViscousBuilder::RestoreListeners()
1151 //================================================================================
1153 * \brief computes SMESH_ProxyMesh::SubMesh::_n2n
1155 //================================================================================
1157 bool _ViscousBuilder::MakeN2NMap( _MeshOfSolid* pm )
1159 SMESH_subMesh* solidSM = pm->mySubMeshes.front();
1160 TopExp_Explorer fExp( solidSM->GetSubShape(), TopAbs_FACE );
1161 for ( ; fExp.More(); fExp.Next() )
1163 SMESHDS_SubMesh* srcSmDS = pm->GetMeshDS()->MeshElements( fExp.Current() );
1164 const SMESH_ProxyMesh::SubMesh* prxSmDS = pm->GetProxySubMesh( fExp.Current() );
1166 if ( !srcSmDS || !prxSmDS || !srcSmDS->NbElements() || !prxSmDS->NbElements() )
1168 if ( srcSmDS->GetElements()->next() == prxSmDS->GetElements()->next())
1171 if ( srcSmDS->NbElements() != prxSmDS->NbElements() )
1172 return error( "Different nb elements in a source and a proxy sub-mesh", solidSM->GetId());
1174 SMDS_ElemIteratorPtr srcIt = srcSmDS->GetElements();
1175 SMDS_ElemIteratorPtr prxIt = prxSmDS->GetElements();
1176 while( prxIt->more() )
1178 const SMDS_MeshElement* fSrc = srcIt->next();
1179 const SMDS_MeshElement* fPrx = prxIt->next();
1180 if ( fSrc->NbNodes() != fPrx->NbNodes())
1181 return error( "Different elements in a source and a proxy sub-mesh", solidSM->GetId());
1182 for ( int i = 0 ; i < fPrx->NbNodes(); ++i )
1183 pm->setNode2Node( fSrc->GetNode(i), fPrx->GetNode(i), prxSmDS );
1186 pm->_n2nMapComputed = true;
1190 //================================================================================
1192 * \brief Does its job
1194 //================================================================================
1196 SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh,
1197 const TopoDS_Shape& theShape)
1199 // TODO: set priority of solids during Gen::Compute()
1203 // check if proxy mesh already computed
1204 TopExp_Explorer exp( theShape, TopAbs_SOLID );
1206 return error("No SOLID's in theShape"), _error;
1208 if ( _ViscousListener::GetSolidMesh( _mesh, exp.Current(), /*toCreate=*/false))
1209 return SMESH_ComputeErrorPtr(); // everything already computed
1213 // TODO: ignore already computed SOLIDs
1214 if ( !findSolidsWithLayers())
1217 if ( !findFacesWithLayers() )
1220 for ( size_t i = 0; i < _sdVec.size(); ++i )
1222 if ( ! makeLayer(_sdVec[i]) )
1225 if ( _sdVec[i]._edges.size() == 0 )
1228 if ( ! inflate(_sdVec[i]) )
1231 if ( ! refine(_sdVec[i]) )
1237 addBoundaryElements();
1239 makeGroupOfLE(); // debug
1245 //================================================================================
1247 * \brief Finds SOLIDs to compute using viscous layers. Fills _sdVec
1249 //================================================================================
1251 bool _ViscousBuilder::findSolidsWithLayers()
1254 TopTools_IndexedMapOfShape allSolids;
1255 TopExp::MapShapes( _mesh->GetShapeToMesh(), TopAbs_SOLID, allSolids );
1256 _sdVec.reserve( allSolids.Extent());
1258 SMESH_Gen* gen = _mesh->GetGen();
1259 SMESH_HypoFilter filter;
1260 for ( int i = 1; i <= allSolids.Extent(); ++i )
1262 // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
1263 SMESH_Algo* algo = gen->GetAlgo( *_mesh, allSolids(i) );
1264 if ( !algo ) continue;
1265 // TODO: check if algo is hidden
1266 const list <const SMESHDS_Hypothesis *> & allHyps =
1267 algo->GetUsedHypothesis(*_mesh, allSolids(i), /*ignoreAuxiliary=*/false);
1268 list< const SMESHDS_Hypothesis *>::const_iterator hyp = allHyps.begin();
1269 const StdMeshers_ViscousLayers* viscHyp = 0;
1270 for ( ; hyp != allHyps.end() && !viscHyp; ++hyp )
1271 viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp );
1274 TopoDS_Shape hypShape;
1275 filter.Init( filter.Is( viscHyp ));
1276 _mesh->GetHypothesis( allSolids(i), filter, true, &hypShape );
1278 _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
1281 _sdVec.push_back( _SolidData( allSolids(i), viscHyp, hypShape, proxyMesh ));
1282 _sdVec.back()._index = getMeshDS()->ShapeToIndex( allSolids(i));
1285 if ( _sdVec.empty() )
1287 ( SMESH_Comment(StdMeshers_ViscousLayers::GetHypType()) << " hypothesis not found",0);
1292 //================================================================================
1296 //================================================================================
1298 bool _ViscousBuilder::findFacesWithLayers()
1300 SMESH_MesherHelper helper( *_mesh );
1301 TopExp_Explorer exp;
1302 TopTools_IndexedMapOfShape solids;
1304 // collect all faces to ignore defined by hyp
1305 for ( size_t i = 0; i < _sdVec.size(); ++i )
1307 solids.Add( _sdVec[i]._solid );
1309 vector<TGeomID> ids = _sdVec[i]._hyp->GetBndShapes();
1310 if ( _sdVec[i]._hyp->IsToIgnoreShapes() ) // FACEs to ignore are given
1312 for ( size_t ii = 0; ii < ids.size(); ++ii )
1314 const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[ii] );
1315 if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
1316 _sdVec[i]._ignoreFaceIds.insert( ids[ii] );
1319 else // FACEs with layers are given
1321 exp.Init( _sdVec[i]._solid, TopAbs_FACE );
1322 for ( ; exp.More(); exp.Next() )
1324 TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
1325 if ( find( ids.begin(), ids.end(), faceInd ) == ids.end() )
1326 _sdVec[i]._ignoreFaceIds.insert( faceInd );
1330 // ignore internal FACEs if inlets and outlets are specified
1332 TopTools_IndexedDataMapOfShapeListOfShape solidsOfFace;
1333 if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
1334 TopExp::MapShapesAndAncestors( _sdVec[i]._hypShape,
1335 TopAbs_FACE, TopAbs_SOLID, solidsOfFace);
1337 exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
1338 for ( ; exp.More(); exp.Next() )
1340 const TopoDS_Face& face = TopoDS::Face( exp.Current() );
1341 if ( helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) < 2 )
1344 const TGeomID faceInd = getMeshDS()->ShapeToIndex( face );
1345 if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
1347 int nbSolids = solidsOfFace.FindFromKey( face ).Extent();
1349 _sdVec[i]._ignoreFaceIds.insert( faceInd );
1352 if ( helper.IsReversedSubMesh( face ))
1354 _sdVec[i]._reversedFaceIds.insert( faceInd );
1360 // Find faces to shrink mesh on (solution 2 in issue 0020832);
1361 TopTools_IndexedMapOfShape shapes;
1362 for ( size_t i = 0; i < _sdVec.size(); ++i )
1365 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_EDGE, shapes);
1366 for ( int iE = 1; iE <= shapes.Extent(); ++iE )
1368 const TopoDS_Shape& edge = shapes(iE);
1369 // find 2 faces sharing an edge
1371 PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE);
1372 while ( fIt->more())
1374 const TopoDS_Shape* f = fIt->next();
1375 if ( helper.IsSubShape( *f, _sdVec[i]._solid))
1376 FF[ int( !FF[0].IsNull()) ] = *f;
1378 if( FF[1].IsNull() ) continue; // seam edge can be shared by 1 FACE only
1379 // check presence of layers on them
1381 for ( int j = 0; j < 2; ++j )
1382 ignore[j] = _sdVec[i]._ignoreFaceIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
1383 if ( ignore[0] == ignore[1] )
1384 continue; // nothing interesting
1385 TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ];
1386 // check presence of layers on fWOL within an adjacent SOLID
1387 bool collision = false;
1388 PShapeIteratorPtr sIt = helper.GetAncestors( fWOL, *_mesh, TopAbs_SOLID );
1389 while ( const TopoDS_Shape* solid = sIt->next() )
1390 if ( !solid->IsSame( _sdVec[i]._solid ))
1392 int iSolid = solids.FindIndex( *solid );
1393 int iFace = getMeshDS()->ShapeToIndex( fWOL );
1394 if ( iSolid > 0 && !_sdVec[ iSolid-1 ]._ignoreFaceIds.count( iFace ))
1396 //_sdVec[i]._noShrinkShapes.insert( iFace );
1402 if ( !fWOL.IsNull())
1404 TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
1405 _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
1408 // _shrinkShape2Shape will be used to temporary inflate _LayerEdge's based
1409 // on the edge but shrink won't be performed
1410 _sdVec[i]._noShrinkShapes.insert( edgeInd );
1415 // Exclude from _shrinkShape2Shape FACE's that can't be shrinked since
1416 // the algo of the SOLID sharing the FACE does not support it
1417 set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D");
1418 for ( size_t i = 0; i < _sdVec.size(); ++i )
1420 map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
1421 for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
1423 const TopoDS_Shape& fWOL = e2f->second;
1424 const TGeomID edgeID = e2f->first;
1425 bool notShrinkFace = false;
1426 PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
1427 while ( soIt->more() )
1429 const TopoDS_Shape* solid = soIt->next();
1430 if ( _sdVec[i]._solid.IsSame( *solid )) continue;
1431 SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid );
1432 if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue;
1433 notShrinkFace = true;
1435 for ( ; iSolid < _sdVec.size(); ++iSolid )
1437 if ( _sdVec[iSolid]._solid.IsSame( *solid ) ) {
1438 if ( _sdVec[iSolid]._shrinkShape2Shape.count( edgeID ))
1439 notShrinkFace = false;
1443 if ( notShrinkFace )
1445 _sdVec[i]._noShrinkShapes.insert( edgeID );
1447 // add VERTEXes of the edge in _noShrinkShapes
1448 TopoDS_Shape edge = getMeshDS()->IndexToShape( edgeID );
1449 for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() )
1450 _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( vIt.Value() ));
1452 // check if there is a collision with to-shrink-from EDGEs in iSolid
1453 if ( iSolid == _sdVec.size() )
1454 continue; // no VL in the solid
1456 TopExp::MapShapes( fWOL, TopAbs_EDGE, shapes);
1457 for ( int iE = 1; iE <= shapes.Extent(); ++iE )
1459 const TopoDS_Edge& E = TopoDS::Edge( shapes( iE ));
1460 const TGeomID eID = getMeshDS()->ShapeToIndex( E );
1461 if ( eID == edgeID ||
1462 !_sdVec[iSolid]._shrinkShape2Shape.count( eID ) ||
1463 _sdVec[i]._noShrinkShapes.count( eID ))
1465 for ( int is1st = 0; is1st < 2; ++is1st )
1467 TopoDS_Vertex V = helper.IthVertex( is1st, E );
1468 if ( _sdVec[i]._noShrinkShapes.count( getMeshDS()->ShapeToIndex( V ) ))
1470 // _sdVec[i]._noShrinkShapes.insert( eID );
1471 // V = helper.IthVertex( !is1st, E );
1472 // _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( V ));
1473 //iE = 0; // re-start the loop on EDGEs of fWOL
1474 return error("No way to make a conformal mesh with "
1475 "the given set of faces with layers", _sdVec[i]._index);
1481 } // while ( soIt->more() )
1482 } // loop on _sdVec[i]._shrinkShape2Shape
1483 } // loop on _sdVec to fill in _SolidData::_noShrinkShapes
1485 // Find the SHAPE along which to inflate _LayerEdge based on VERTEX
1487 for ( size_t i = 0; i < _sdVec.size(); ++i )
1490 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_VERTEX, shapes);
1491 for ( int iV = 1; iV <= shapes.Extent(); ++iV )
1493 const TopoDS_Shape& vertex = shapes(iV);
1494 // find faces WOL sharing the vertex
1495 vector< TopoDS_Shape > facesWOL;
1496 int totalNbFaces = 0;
1497 PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE);
1498 while ( fIt->more())
1500 const TopoDS_Shape* f = fIt->next();
1501 if ( helper.IsSubShape( *f, _sdVec[i]._solid ) )
1504 const int fID = getMeshDS()->ShapeToIndex( *f );
1505 if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&&
1506 !_sdVec[i]._noShrinkShapes.count( fID )*/)
1507 facesWOL.push_back( *f );
1510 if ( facesWOL.size() == totalNbFaces || facesWOL.empty() )
1511 continue; // no layers at this vertex or no WOL
1512 TGeomID vInd = getMeshDS()->ShapeToIndex( vertex );
1513 switch ( facesWOL.size() )
1517 helper.SetSubShape( facesWOL[0] );
1518 if ( helper.IsRealSeam( vInd )) // inflate along a seam edge?
1520 TopoDS_Shape seamEdge;
1521 PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1522 while ( eIt->more() && seamEdge.IsNull() )
1524 const TopoDS_Shape* e = eIt->next();
1525 if ( helper.IsRealSeam( *e ) )
1528 if ( !seamEdge.IsNull() )
1530 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge ));
1534 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] ));
1539 // find an edge shared by 2 faces
1540 PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1541 while ( eIt->more())
1543 const TopoDS_Shape* e = eIt->next();
1544 if ( helper.IsSubShape( *e, facesWOL[0]) &&
1545 helper.IsSubShape( *e, facesWOL[1]))
1547 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break;
1553 return error("Not yet supported case", _sdVec[i]._index);
1558 // add FACEs of other SOLIDs to _ignoreFaceIds
1559 for ( size_t i = 0; i < _sdVec.size(); ++i )
1562 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_FACE, shapes);
1564 for ( exp.Init( _mesh->GetShapeToMesh(), TopAbs_FACE ); exp.More(); exp.Next() )
1566 if ( !shapes.Contains( exp.Current() ))
1567 _sdVec[i]._ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( exp.Current() ));
1574 //================================================================================
1576 * \brief Create the inner surface of the viscous layer and prepare data for infation
1578 //================================================================================
1580 bool _ViscousBuilder::makeLayer(_SolidData& data)
1582 // get all sub-shapes to make layers on
1583 set<TGeomID> subIds, faceIds;
1584 subIds = data._noShrinkShapes;
1585 TopExp_Explorer exp( data._solid, TopAbs_FACE );
1586 for ( ; exp.More(); exp.Next() )
1588 SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
1589 if ( ! data._ignoreFaceIds.count( fSubM->GetId() ))
1590 faceIds.insert( fSubM->GetId() );
1591 SMESH_subMeshIteratorPtr subIt = fSubM->getDependsOnIterator(/*includeSelf=*/true);
1592 while ( subIt->more() )
1593 subIds.insert( subIt->next()->GetId() );
1596 // make a map to find new nodes on sub-shapes shared with other SOLID
1597 map< TGeomID, TNode2Edge* >::iterator s2ne;
1598 map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
1599 for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
1601 TGeomID shapeInd = s2s->first;
1602 for ( size_t i = 0; i < _sdVec.size(); ++i )
1604 if ( _sdVec[i]._index == data._index ) continue;
1605 map< TGeomID, TopoDS_Shape >::iterator s2s2 = _sdVec[i]._shrinkShape2Shape.find( shapeInd );
1606 if ( s2s2 != _sdVec[i]._shrinkShape2Shape.end() &&
1607 *s2s == *s2s2 && !_sdVec[i]._n2eMap.empty() )
1609 data._s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap ));
1615 // Create temporary faces and _LayerEdge's
1617 dumpFunction(SMESH_Comment("makeLayers_")<<data._index);
1619 data._stepSize = Precision::Infinite();
1620 data._stepSizeNodes[0] = 0;
1622 SMESH_MesherHelper helper( *_mesh );
1623 helper.SetSubShape( data._solid );
1624 helper.SetElementsOnShape( true );
1626 vector< const SMDS_MeshNode*> newNodes; // of a mesh face
1627 TNode2Edge::iterator n2e2;
1629 // collect _LayerEdge's of shapes they are based on
1630 const int nbShapes = getMeshDS()->MaxShapeIndex();
1631 vector< vector<_LayerEdge*> > edgesByGeom( nbShapes+1 );
1633 for ( set<TGeomID>::iterator id = faceIds.begin(); id != faceIds.end(); ++id )
1635 SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( *id );
1636 if ( !smDS ) return error(SMESH_Comment("Not meshed face ") << *id, data._index );
1638 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( *id ));
1639 SMESH_ProxyMesh::SubMesh* proxySub =
1640 data._proxyMesh->getFaceSubM( F, /*create=*/true);
1642 SMDS_ElemIteratorPtr eIt = smDS->GetElements();
1643 while ( eIt->more() )
1645 const SMDS_MeshElement* face = eIt->next();
1646 newNodes.resize( face->NbCornerNodes() );
1647 double faceMaxCosin = -1;
1648 _LayerEdge* maxCosinEdge = 0;
1649 for ( int i = 0 ; i < face->NbCornerNodes(); ++i )
1651 const SMDS_MeshNode* n = face->GetNode(i);
1652 TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first;
1653 if ( !(*n2e).second )
1656 _LayerEdge* edge = new _LayerEdge();
1658 edge->_nodes.push_back( n );
1659 const int shapeID = n->getshapeId();
1660 const bool noShrink = data._noShrinkShapes.count( shapeID );
1661 edgesByGeom[ shapeID ].push_back( edge );
1663 SMESH_TNodeXYZ xyz( n );
1665 // set edge data or find already refined _LayerEdge and get data from it
1666 if (( !noShrink ) &&
1667 ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE ) &&
1668 (( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() ) &&
1669 (( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end() ))
1671 _LayerEdge* foundEdge = (*n2e2).second;
1672 gp_XYZ lastPos = edge->Copy( *foundEdge, helper );
1673 foundEdge->_pos.push_back( lastPos );
1674 // location of the last node is modified and we restore it by foundEdge->_pos.back()
1675 const_cast< SMDS_MeshNode* >
1676 ( edge->_nodes.back() )->setXYZ( xyz.X(), xyz.Y(), xyz.Z() );
1682 edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ));
1684 if ( !setEdgeData( *edge, subIds, helper, data ))
1687 dumpMove(edge->_nodes.back());
1689 if ( edge->_cosin > faceMaxCosin )
1691 faceMaxCosin = edge->_cosin;
1692 maxCosinEdge = edge;
1695 newNodes[ i ] = n2e->second->_nodes.back();
1697 // create a temporary face
1698 const SMDS_MeshElement* newFace =
1699 new _TmpMeshFace( newNodes, --_tmpFaceID, face->getshapeId() );
1700 proxySub->AddElement( newFace );
1702 // compute inflation step size by min size of element on a convex surface
1703 if ( faceMaxCosin > theMinSmoothCosin )
1704 limitStepSize( data, face, maxCosinEdge );
1705 } // loop on 2D elements on a FACE
1706 } // loop on FACEs of a SOLID
1708 data._epsilon = 1e-7;
1709 if ( data._stepSize < 1. )
1710 data._epsilon *= data._stepSize;
1712 // Put _LayerEdge's into the vector data._edges
1713 if ( !sortEdges( data, edgesByGeom ))
1716 // limit data._stepSize depending on surface curvature and fill data._convexFaces
1717 limitStepSizeByCurvature( data ); // !!! it must be before node substitution in _Simplex
1719 // Set target nodes into _Simplex and _LayerEdge's to _2NearEdges
1720 TNode2Edge::iterator n2e;
1721 const SMDS_MeshNode* nn[2];
1722 for ( size_t i = 0; i < data._edges.size(); ++i )
1724 if ( data._edges[i]->IsOnEdge() )
1726 // get neighbor nodes
1727 bool hasData = ( data._edges[i]->_2neibors->_edges[0] );
1728 if ( hasData ) // _LayerEdge is a copy of another one
1730 nn[0] = data._edges[i]->_2neibors->srcNode(0);
1731 nn[1] = data._edges[i]->_2neibors->srcNode(1);
1733 else if ( !findNeiborsOnEdge( data._edges[i], nn[0],nn[1], data ))
1737 // set neighbor _LayerEdge's
1738 for ( int j = 0; j < 2; ++j )
1740 if (( n2e = data._n2eMap.find( nn[j] )) == data._n2eMap.end() )
1741 return error("_LayerEdge not found by src node", data._index);
1742 data._edges[i]->_2neibors->_edges[j] = n2e->second;
1745 data._edges[i]->SetDataByNeighbors( nn[0], nn[1], helper);
1748 for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j )
1750 _Simplex& s = data._edges[i]->_simplices[j];
1751 s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
1752 s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
1760 //================================================================================
1762 * \brief Compute inflation step size by min size of element on a convex surface
1764 //================================================================================
1766 void _ViscousBuilder::limitStepSize( _SolidData& data,
1767 const SMDS_MeshElement* face,
1768 const _LayerEdge* maxCosinEdge )
1771 double minSize = 10 * data._stepSize;
1772 const int nbNodes = face->NbCornerNodes();
1773 for ( int i = 0; i < nbNodes; ++i )
1775 const SMDS_MeshNode* nextN = face->GetNode( SMESH_MesherHelper::WrapIndex( i+1, nbNodes ));
1776 const SMDS_MeshNode* curN = face->GetNode( i );
1777 if ( nextN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ||
1778 curN-> GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
1780 double dist = SMESH_TNodeXYZ( curN ).Distance( nextN );
1781 if ( dist < minSize )
1782 minSize = dist, iN = i;
1785 double newStep = 0.8 * minSize / maxCosinEdge->_lenFactor;
1786 if ( newStep < data._stepSize )
1788 data._stepSize = newStep;
1789 data._stepSizeCoeff = 0.8 / maxCosinEdge->_lenFactor;
1790 data._stepSizeNodes[0] = face->GetNode( iN );
1791 data._stepSizeNodes[1] = face->GetNode( SMESH_MesherHelper::WrapIndex( iN+1, nbNodes ));
1795 //================================================================================
1797 * \brief Compute inflation step size by min size of element on a convex surface
1799 //================================================================================
1801 void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize )
1803 if ( minSize < data._stepSize )
1805 data._stepSize = minSize;
1806 if ( data._stepSizeNodes[0] )
1809 SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
1810 data._stepSizeCoeff = data._stepSize / dist;
1815 //================================================================================
1817 * \brief Limit data._stepSize by evaluating curvature of shapes and fill data._convexFaces
1819 //================================================================================
1821 void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
1823 const int nbTestPnt = 5; // on a FACE sub-shape
1824 const double minCurvature = 0.9 / data._hyp->GetTotalThickness();
1826 BRepLProp_SLProps surfProp( 2, 1e-6 );
1827 SMESH_MesherHelper helper( *_mesh );
1829 data._convexFaces.clear();
1831 TopExp_Explorer face( data._solid, TopAbs_FACE );
1832 for ( ; face.More(); face.Next() )
1834 const TopoDS_Face& F = TopoDS::Face( face.Current() );
1835 SMESH_subMesh * sm = _mesh->GetSubMesh( F );
1836 const TGeomID faceID = sm->GetId();
1837 if ( data._ignoreFaceIds.count( faceID )) continue;
1839 BRepAdaptor_Surface surface( F, false );
1840 surfProp.SetSurface( surface );
1842 bool isTooCurved = false;
1845 _ConvexFace cnvFace;
1846 const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. );
1847 SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
1848 while ( smIt->more() )
1851 const TGeomID subID = sm->GetId();
1852 // find _LayerEdge's of a sub-shape
1854 if ( data.GetShapeEdges( subID, edgesEnd, &iBeg, &iEnd ))
1855 cnvFace._subIdToEdgeEnd.insert( make_pair( subID, edgesEnd ));
1858 // check concavity and curvature and limit data._stepSize
1859 int nbLEdges = iEnd - iBeg;
1860 int iStep = Max( 1, nbLEdges / nbTestPnt );
1861 for ( ; iBeg < iEnd; iBeg += iStep )
1863 gp_XY uv = helper.GetNodeUV( F, data._edges[ iBeg ]->_nodes[0] );
1864 surfProp.SetParameters( uv.X(), uv.Y() );
1865 if ( !surfProp.IsCurvatureDefined() )
1867 if ( surfProp.MaxCurvature() * oriFactor > minCurvature )
1869 limitStepSize( data, 0.9 / surfProp.MaxCurvature() * oriFactor );
1872 if ( surfProp.MinCurvature() * oriFactor > minCurvature )
1874 limitStepSize( data, 0.9 / surfProp.MinCurvature() * oriFactor );
1878 } // loop on sub-shapes of the FACE
1880 if ( !isTooCurved ) continue;
1882 _ConvexFace & convFace =
1883 data._convexFaces.insert( make_pair( faceID, cnvFace )).first->second;
1886 convFace._normalsFixed = false;
1888 // Fill _ConvexFace::_simplexTestEdges. These _LayerEdge's are used to detect
1889 // prism distortion.
1890 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.find( faceID );
1891 if ( id2end != convFace._subIdToEdgeEnd.end() )
1893 // there are _LayerEdge's on the FACE it-self;
1894 // select _LayerEdge's near EDGEs
1895 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
1896 for ( ; iBeg < iEnd; ++iBeg )
1898 _LayerEdge* ledge = data._edges[ iBeg ];
1899 for ( size_t j = 0; j < ledge->_simplices.size(); ++j )
1900 if ( ledge->_simplices[j]._nNext->GetPosition()->GetDim() < 2 )
1902 convFace._simplexTestEdges.push_back( ledge );
1909 // where there are no _LayerEdge's on a _ConvexFace,
1910 // as e.g. on a fillet surface with no internal nodes - issue 22580,
1911 // so that collision of viscous internal faces is not detected by check of
1912 // intersection of _LayerEdge's with the viscous internal faces.
1914 set< const SMDS_MeshNode* > usedNodes;
1916 // look for _LayerEdge's with null _sWOL
1917 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.begin();
1918 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
1920 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
1921 if ( iBeg >= iEnd || !data._edges[ iBeg ]->_sWOL.IsNull() )
1923 for ( ; iBeg < iEnd; ++iBeg )
1925 _LayerEdge* ledge = data._edges[ iBeg ];
1926 const SMDS_MeshNode* srcNode = ledge->_nodes[0];
1927 if ( !usedNodes.insert( srcNode ).second ) continue;
1929 getSimplices( srcNode, ledge->_simplices, data._ignoreFaceIds, &data );
1930 for ( size_t i = 0; i < ledge->_simplices.size(); ++i )
1932 usedNodes.insert( ledge->_simplices[i]._nPrev );
1933 usedNodes.insert( ledge->_simplices[i]._nNext );
1935 convFace._simplexTestEdges.push_back( ledge );
1939 } // loop on FACEs of data._solid
1942 //================================================================================
1944 * \brief Separate shapes (and _LayerEdge's on them) to smooth from the rest ones
1946 //================================================================================
1948 bool _ViscousBuilder::sortEdges( _SolidData& data,
1949 vector< vector<_LayerEdge*> >& edgesByGeom)
1951 // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's
1952 // boundry inclined at a sharp angle to the shape
1954 list< TGeomID > shapesToSmooth;
1956 SMESH_MesherHelper helper( *_mesh );
1959 for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
1961 vector<_LayerEdge*>& eS = edgesByGeom[iS];
1962 if ( eS.empty() ) continue;
1963 const TopoDS_Shape& S = getMeshDS()->IndexToShape( iS );
1964 bool needSmooth = false;
1965 switch ( S.ShapeType() )
1969 bool isShrinkEdge = !eS[0]->_sWOL.IsNull();
1970 for ( TopoDS_Iterator vIt( S ); vIt.More() && !needSmooth; vIt.Next() )
1972 TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
1973 vector<_LayerEdge*>& eV = edgesByGeom[ iV ];
1974 if ( eV.empty() ) continue;
1975 // double cosin = eV[0]->_cosin;
1977 // ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge));
1980 // gp_Vec dir1, dir2;
1981 // if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE )
1982 // dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ));
1984 // dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ),
1985 // eV[0]->_nodes[0], helper, ok);
1986 // dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
1987 // double angle = dir1.Angle( dir2 );
1988 // cosin = cos( angle );
1990 gp_Vec eDir = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
1991 double angle = eDir.Angle( eV[0]->_normal );
1992 double cosin = Cos( angle );
1993 needSmooth = ( cosin > theMinSmoothCosin );
1999 for ( TopExp_Explorer eExp( S, TopAbs_EDGE ); eExp.More() && !needSmooth; eExp.Next() )
2001 TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
2002 vector<_LayerEdge*>& eE = edgesByGeom[ iE ];
2003 if ( eE.empty() ) continue;
2004 if ( eE[0]->_sWOL.IsNull() )
2006 for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
2007 needSmooth = ( eE[i]->_cosin > theMinSmoothCosin );
2011 const TopoDS_Face& F1 = TopoDS::Face( S );
2012 const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
2013 const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
2014 for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
2016 gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
2017 gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
2018 double angle = dir1.Angle( dir2 );
2019 double cosin = cos( angle );
2020 needSmooth = ( cosin > theMinSmoothCosin );
2032 if ( S.ShapeType() == TopAbs_EDGE ) shapesToSmooth.push_front( iS );
2033 else shapesToSmooth.push_back ( iS );
2036 } // loop on edgesByGeom
2038 data._edges.reserve( data._n2eMap.size() );
2039 data._endEdgeOnShape.clear();
2041 // first we put _LayerEdge's on shapes to smooth
2042 data._nbShapesToSmooth = 0;
2043 list< TGeomID >::iterator gIt = shapesToSmooth.begin();
2044 for ( ; gIt != shapesToSmooth.end(); ++gIt )
2046 vector<_LayerEdge*>& eVec = edgesByGeom[ *gIt ];
2047 if ( eVec.empty() ) continue;
2048 data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
2049 data._endEdgeOnShape.push_back( data._edges.size() );
2050 data._nbShapesToSmooth++;
2054 // then the rest _LayerEdge's
2055 for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
2057 vector<_LayerEdge*>& eVec = edgesByGeom[iS];
2058 if ( eVec.empty() ) continue;
2059 data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
2060 data._endEdgeOnShape.push_back( data._edges.size() );
2067 //================================================================================
2069 * \brief Set data of _LayerEdge needed for smoothing
2070 * \param subIds - ids of sub-shapes of a SOLID to take into account faces from
2072 //================================================================================
2074 bool _ViscousBuilder::setEdgeData(_LayerEdge& edge,
2075 const set<TGeomID>& subIds,
2076 SMESH_MesherHelper& helper,
2079 SMESH_MeshEditor editor(_mesh);
2081 const SMDS_MeshNode* node = edge._nodes[0]; // source node
2082 SMDS_TypeOfPosition posType = node->GetPosition()->GetTypeOfPosition();
2086 edge._curvature = 0;
2088 // --------------------------
2089 // Compute _normal and _cosin
2090 // --------------------------
2093 edge._normal.SetCoord(0,0,0);
2095 int totalNbFaces = 0;
2099 const TGeomID shapeInd = node->getshapeId();
2100 map< TGeomID, TopoDS_Shape >::const_iterator s2s = data._shrinkShape2Shape.find( shapeInd );
2101 const bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
2103 if ( onShrinkShape ) // one of faces the node is on has no layers
2105 TopoDS_Shape vertEdge = getMeshDS()->IndexToShape( s2s->first ); // vertex or edge
2106 if ( s2s->second.ShapeType() == TopAbs_EDGE )
2108 // inflate from VERTEX along EDGE
2109 edge._normal = getEdgeDir( TopoDS::Edge( s2s->second ), TopoDS::Vertex( vertEdge ));
2111 else if ( vertEdge.ShapeType() == TopAbs_VERTEX )
2113 // inflate from VERTEX along FACE
2114 edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Vertex( vertEdge ),
2115 node, helper, normOK, &edge._cosin);
2119 // inflate from EDGE along FACE
2120 edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Edge( vertEdge ),
2121 node, helper, normOK);
2124 else // layers are on all faces of SOLID the node is on
2126 // find indices of geom faces the node lies on
2127 set<TGeomID> faceIds;
2128 if ( posType == SMDS_TOP_FACE )
2130 faceIds.insert( node->getshapeId() );
2134 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2135 while ( fIt->more() )
2136 faceIds.insert( editor.FindShape(fIt->next()));
2139 set<TGeomID>::iterator id = faceIds.begin();
2141 std::pair< TGeomID, gp_XYZ > id2Norm[20];
2142 for ( ; id != faceIds.end(); ++id )
2144 const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id );
2145 if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
2147 F = TopoDS::Face( s );
2148 geomNorm = getFaceNormal( node, F, helper, normOK );
2149 if ( !normOK ) continue;
2151 if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
2153 id2Norm[ totalNbFaces ].first = *id;
2154 id2Norm[ totalNbFaces ].second = geomNorm.XYZ();
2156 edge._normal += geomNorm.XYZ();
2158 if ( totalNbFaces == 0 )
2159 return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index);
2161 if ( normOK && edge._normal.Modulus() < 1e-3 && totalNbFaces > 1 )
2163 // opposite normals, re-get normals at shifted positions (IPAL 52426)
2164 edge._normal.SetCoord( 0,0,0 );
2165 for ( int i = 0; i < totalNbFaces; ++i )
2167 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( id2Norm[i].first ));
2168 geomNorm = getFaceNormal( node, F, helper, normOK, /*shiftInside=*/true );
2169 if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
2172 id2Norm[ i ].second = geomNorm.XYZ();
2173 edge._normal += id2Norm[ i ].second;
2177 if ( totalNbFaces < 3 )
2179 //edge._normal /= totalNbFaces;
2183 edge._normal = getWeigthedNormal( node, id2Norm, totalNbFaces );
2189 edge._cosin = 0; break;
2191 case SMDS_TOP_EDGE: {
2192 TopoDS_Edge E = TopoDS::Edge( helper.GetSubShapeByNode( node, getMeshDS()));
2193 gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK );
2194 double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
2195 edge._cosin = cos( angle );
2196 //cout << "Cosin on EDGE " << edge._cosin << " node " << node->GetID() << endl;
2199 case SMDS_TOP_VERTEX: {
2200 TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS()));
2201 gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK );
2202 double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
2203 edge._cosin = cos( angle );
2204 //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl;
2208 return error(SMESH_Comment("Invalid shape position of node ")<<node, data._index);
2210 } // case _sWOL.IsNull()
2212 double normSize = edge._normal.SquareModulus();
2213 if ( normSize < numeric_limits<double>::min() )
2214 return error(SMESH_Comment("Bad normal at node ")<< node->GetID(), data._index );
2216 edge._normal /= sqrt( normSize );
2218 // TODO: if ( !normOK ) then get normal by mesh faces
2220 // Set the rest data
2221 // --------------------
2222 if ( onShrinkShape )
2224 edge._sWOL = (*s2s).second;
2226 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( edge._nodes.back() );
2227 if ( SMESHDS_SubMesh* sm = getMeshDS()->MeshElements( data._solid ))
2228 sm->RemoveNode( tgtNode , /*isNodeDeleted=*/false );
2230 // set initial position which is parameters on _sWOL in this case
2231 if ( edge._sWOL.ShapeType() == TopAbs_EDGE )
2233 double u = helper.GetNodeU( TopoDS::Edge( edge._sWOL ), node, 0, &normOK );
2234 edge._pos.push_back( gp_XYZ( u, 0, 0 ));
2235 if ( edge._nodes.size() > 1 )
2236 getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( edge._sWOL ), u );
2240 gp_XY uv = helper.GetNodeUV( TopoDS::Face( edge._sWOL ), node, 0, &normOK );
2241 edge._pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
2242 if ( edge._nodes.size() > 1 )
2243 getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( edge._sWOL ), uv.X(), uv.Y() );
2248 edge._pos.push_back( SMESH_TNodeXYZ( node ));
2250 if ( posType == SMDS_TOP_FACE )
2252 getSimplices( node, edge._simplices, data._ignoreFaceIds, &data );
2253 double avgNormProj = 0, avgLen = 0;
2254 for ( size_t i = 0; i < edge._simplices.size(); ++i )
2256 gp_XYZ vec = edge._pos.back() - SMESH_TNodeXYZ( edge._simplices[i]._nPrev );
2257 avgNormProj += edge._normal * vec;
2258 avgLen += vec.Modulus();
2260 avgNormProj /= edge._simplices.size();
2261 avgLen /= edge._simplices.size();
2262 edge._curvature = _Curvature::New( avgNormProj, avgLen );
2266 // Set neighbour nodes for a _LayerEdge based on EDGE
2268 if ( posType == SMDS_TOP_EDGE /*||
2269 ( onShrinkShape && posType == SMDS_TOP_VERTEX && fabs( edge._cosin ) < 1e-10 )*/)
2271 edge._2neibors = new _2NearEdges;
2272 // target node instead of source ones will be set later
2273 // if ( ! findNeiborsOnEdge( &edge,
2274 // edge._2neibors->_nodes[0],
2275 // edge._2neibors->_nodes[1],
2278 // edge.SetDataByNeighbors( edge._2neibors->_nodes[0],
2279 // edge._2neibors->_nodes[1],
2283 edge.SetCosin( edge._cosin ); // to update edge._lenFactor
2288 //================================================================================
2290 * \brief Return normal to a FACE at a node
2291 * \param [in] n - node
2292 * \param [in] face - FACE
2293 * \param [in] helper - helper
2294 * \param [out] isOK - true or false
2295 * \param [in] shiftInside - to find normal at a position shifted inside the face
2296 * \return gp_XYZ - normal
2298 //================================================================================
2300 gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
2301 const TopoDS_Face& face,
2302 SMESH_MesherHelper& helper,
2309 // get a shifted position
2310 gp_Pnt p = SMESH_TNodeXYZ( node );
2311 gp_XYZ shift( 0,0,0 );
2312 TopoDS_Shape S = helper.GetSubShapeByNode( node, helper.GetMeshDS() );
2313 switch ( S.ShapeType() ) {
2316 shift = getFaceDir( face, TopoDS::Vertex( S ), node, helper, isOK );
2321 shift = getFaceDir( face, TopoDS::Edge( S ), node, helper, isOK );
2329 p.Translate( shift * 1e-5 );
2331 TopLoc_Location loc;
2332 GeomAPI_ProjectPointOnSurf& projector = helper.GetProjector( face, loc, 1e-7 );
2334 if ( !loc.IsIdentity() ) p.Transform( loc.Transformation().Inverted() );
2336 projector.Perform( p );
2337 if ( !projector.IsDone() || projector.NbPoints() < 1 )
2342 Quantity_Parameter U,V;
2343 projector.LowerDistanceParameters(U,V);
2348 uv = helper.GetNodeUV( face, node, 0, &isOK );
2354 Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
2355 int pointKind = GeomLib::NormEstim( surface, uv, 1e-5, normal );
2356 enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE };
2357 if ( pointKind < IMPOSSIBLE )
2359 if ( pointKind != REGULAR && !shiftInside )
2361 gp_XYZ normShift = getFaceNormal( node, face, helper, isOK, /*shiftInside=*/true );
2362 if ( normShift * normal.XYZ() < 0. )
2367 else // hard singularity, to call with shiftInside=true ?
2369 const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face );
2371 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2372 while ( fIt->more() )
2374 const SMDS_MeshElement* f = fIt->next();
2375 if ( f->getshapeId() == faceID )
2377 isOK = SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) normal.XYZ(), /*normalized=*/true );
2380 if ( helper.IsReversedSubMesh( face ))
2387 return normal.XYZ();
2390 //================================================================================
2392 * \brief Return a normal at a node weighted with angles taken by FACEs
2393 * \param [in] n - the node
2394 * \param [in] fId2Normal - FACE ids and normals
2395 * \param [in] nbFaces - nb of FACEs meeting at the node
2396 * \return gp_XYZ - computed normal
2398 //================================================================================
2400 gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode* n,
2401 std::pair< TGeomID, gp_XYZ > fId2Normal[],
2404 gp_XYZ resNorm(0,0,0);
2405 TopoDS_Shape V = SMESH_MesherHelper::GetSubShapeByNode( n, getMeshDS() );
2406 if ( V.ShapeType() != TopAbs_VERTEX )
2408 for ( int i = 0; i < nbFaces; ++i )
2409 resNorm += fId2Normal[i].second / nbFaces ;
2414 for ( int i = 0; i < nbFaces; ++i )
2416 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( fId2Normal[i].first ));
2418 // look for two EDGEs shared by F and other FACEs within fId2Normal
2421 PShapeIteratorPtr eIt = SMESH_MesherHelper::GetAncestors( V, *_mesh, TopAbs_EDGE );
2422 while ( const TopoDS_Shape* E = eIt->next() )
2424 if ( !SMESH_MesherHelper::IsSubShape( *E, F ))
2426 bool isSharedEdge = false;
2427 for ( int j = 0; j < nbFaces && !isSharedEdge; ++j )
2429 if ( i == j ) continue;
2430 const TopoDS_Shape& otherF = getMeshDS()->IndexToShape( fId2Normal[j].first );
2431 isSharedEdge = SMESH_MesherHelper::IsSubShape( *E, otherF );
2433 if ( !isSharedEdge )
2435 ee[ nbE ] = TopoDS::Edge( *E );
2436 ee[ nbE ].Orientation( SMESH_MesherHelper::GetSubShapeOri( F, *E ));
2441 // get an angle between the two EDGEs
2443 if ( nbE < 1 ) continue;
2450 if ( !V.IsSame( SMESH_MesherHelper::IthVertex( 0, ee[ 1 ] )))
2451 std::swap( ee[0], ee[1] );
2453 angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F, TopoDS::Vertex( V ));
2456 // compute a weighted normal
2457 double sumAngle = 0;
2458 for ( int i = 0; i < nbFaces; ++i )
2460 angles[i] = ( angles[i] > 2*M_PI ) ? 0 : M_PI - angles[i];
2461 sumAngle += angles[i];
2463 for ( int i = 0; i < nbFaces; ++i )
2464 resNorm += angles[i] / sumAngle * fId2Normal[i].second;
2469 //================================================================================
2471 * \brief Find 2 neigbor nodes of a node on EDGE
2473 //================================================================================
2475 bool _ViscousBuilder::findNeiborsOnEdge(const _LayerEdge* edge,
2476 const SMDS_MeshNode*& n1,
2477 const SMDS_MeshNode*& n2,
2480 const SMDS_MeshNode* node = edge->_nodes[0];
2481 const int shapeInd = node->getshapeId();
2482 SMESHDS_SubMesh* edgeSM = 0;
2483 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE )
2485 edgeSM = getMeshDS()->MeshElements( shapeInd );
2486 if ( !edgeSM || edgeSM->NbElements() == 0 )
2487 return error(SMESH_Comment("Not meshed EDGE ") << shapeInd, data._index);
2491 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Edge);
2492 while ( eIt->more() && !n2 )
2494 const SMDS_MeshElement* e = eIt->next();
2495 const SMDS_MeshNode* nNeibor = e->GetNode( 0 );
2496 if ( nNeibor == node ) nNeibor = e->GetNode( 1 );
2499 if (!edgeSM->Contains(e)) continue;
2503 TopoDS_Shape s = SMESH_MesherHelper::GetSubShapeByNode(nNeibor, getMeshDS() );
2504 if ( !SMESH_MesherHelper::IsSubShape( s, edge->_sWOL )) continue;
2506 ( iN++ ? n2 : n1 ) = nNeibor;
2509 return error(SMESH_Comment("Wrongly meshed EDGE ") << shapeInd, data._index);
2513 //================================================================================
2515 * \brief Set _curvature and _2neibors->_plnNorm by 2 neigbor nodes residing the same EDGE
2517 //================================================================================
2519 void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
2520 const SMDS_MeshNode* n2,
2521 SMESH_MesherHelper& helper)
2523 if ( _nodes[0]->GetPosition()->GetTypeOfPosition() != SMDS_TOP_EDGE )
2526 gp_XYZ pos = SMESH_TNodeXYZ( _nodes[0] );
2527 gp_XYZ vec1 = pos - SMESH_TNodeXYZ( n1 );
2528 gp_XYZ vec2 = pos - SMESH_TNodeXYZ( n2 );
2532 double sumLen = vec1.Modulus() + vec2.Modulus();
2533 _2neibors->_wgt[0] = 1 - vec1.Modulus() / sumLen;
2534 _2neibors->_wgt[1] = 1 - vec2.Modulus() / sumLen;
2535 double avgNormProj = 0.5 * ( _normal * vec1 + _normal * vec2 );
2536 double avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() );
2537 if ( _curvature ) delete _curvature;
2538 _curvature = _Curvature::New( avgNormProj, avgLen );
2539 // if ( _curvature )
2540 // debugMsg( _nodes[0]->GetID()
2541 // << " CURV r,k: " << _curvature->_r<<","<<_curvature->_k
2542 // << " proj = "<<avgNormProj<< " len = " << avgLen << "| lenDelta(0) = "
2543 // << _curvature->lenDelta(0) );
2547 if ( _sWOL.IsNull() )
2549 TopoDS_Shape S = helper.GetSubShapeByNode( _nodes[0], helper.GetMeshDS() );
2550 gp_XYZ dirE = getEdgeDir( TopoDS::Edge( S ), _nodes[0], helper );
2551 gp_XYZ plnNorm = dirE ^ _normal;
2552 double proj0 = plnNorm * vec1;
2553 double proj1 = plnNorm * vec2;
2554 if ( fabs( proj0 ) > 1e-10 || fabs( proj1 ) > 1e-10 )
2556 if ( _2neibors->_plnNorm ) delete _2neibors->_plnNorm;
2557 _2neibors->_plnNorm = new gp_XYZ( plnNorm.Normalized() );
2562 //================================================================================
2564 * \brief Copy data from a _LayerEdge of other SOLID and based on the same node;
2565 * this and other _LayerEdge's are inflated along a FACE or an EDGE
2567 //================================================================================
2569 gp_XYZ _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
2571 _nodes = other._nodes;
2572 _normal = other._normal;
2574 _lenFactor = other._lenFactor;
2575 _cosin = other._cosin;
2576 _sWOL = other._sWOL;
2577 _2neibors = other._2neibors;
2578 _curvature = 0; std::swap( _curvature, other._curvature );
2579 _2neibors = 0; std::swap( _2neibors, other._2neibors );
2581 gp_XYZ lastPos( 0,0,0 );
2582 if ( _sWOL.ShapeType() == TopAbs_EDGE )
2584 double u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes[0] );
2585 _pos.push_back( gp_XYZ( u, 0, 0));
2587 u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes.back() );
2592 gp_XY uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes[0]);
2593 _pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
2595 uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes.back() );
2596 lastPos.SetX( uv.X() );
2597 lastPos.SetY( uv.Y() );
2602 //================================================================================
2604 * \brief Set _cosin and _lenFactor
2606 //================================================================================
2608 void _LayerEdge::SetCosin( double cosin )
2611 cosin = Abs( _cosin );
2612 _lenFactor = ( /*0.1 < cosin &&*/ cosin < 1-1e-12 ) ? 1./sqrt(1-cosin*cosin) : 1.0;
2615 //================================================================================
2617 * \brief Fills a vector<_Simplex >
2619 //================================================================================
2621 void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node,
2622 vector<_Simplex>& simplices,
2623 const set<TGeomID>& ingnoreShapes,
2624 const _SolidData* dataToCheckOri,
2628 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2629 while ( fIt->more() )
2631 const SMDS_MeshElement* f = fIt->next();
2632 const TGeomID shapeInd = f->getshapeId();
2633 if ( ingnoreShapes.count( shapeInd )) continue;
2634 const int nbNodes = f->NbCornerNodes();
2635 const int srcInd = f->GetNodeIndex( node );
2636 const SMDS_MeshNode* nPrev = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd-1, nbNodes ));
2637 const SMDS_MeshNode* nNext = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+1, nbNodes ));
2638 const SMDS_MeshNode* nOpp = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+2, nbNodes ));
2639 if ( dataToCheckOri && dataToCheckOri->_reversedFaceIds.count( shapeInd ))
2640 std::swap( nPrev, nNext );
2641 simplices.push_back( _Simplex( nPrev, nNext, nOpp ));
2646 vector<_Simplex> sortedSimplices( simplices.size() );
2647 sortedSimplices[0] = simplices[0];
2649 for ( size_t i = 1; i < simplices.size(); ++i )
2651 for ( size_t j = 1; j < simplices.size(); ++j )
2652 if ( sortedSimplices[i-1]._nNext == simplices[j]._nPrev )
2654 sortedSimplices[i] = simplices[j];
2659 if ( nbFound == simplices.size() - 1 )
2660 simplices.swap( sortedSimplices );
2664 //================================================================================
2666 * \brief DEBUG. Create groups contating temorary data of _LayerEdge's
2668 //================================================================================
2670 void _ViscousBuilder::makeGroupOfLE()
2673 for ( size_t i = 0 ; i < _sdVec.size(); ++i )
2675 if ( _sdVec[i]._edges.empty() ) continue;
2677 dumpFunction( SMESH_Comment("make_LayerEdge_") << i );
2678 for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
2680 _LayerEdge* le = _sdVec[i]._edges[j];
2681 for ( size_t iN = 1; iN < le->_nodes.size(); ++iN )
2682 dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<le->_nodes[iN-1]->GetID()
2683 << ", " << le->_nodes[iN]->GetID() <<"])");
2687 dumpFunction( SMESH_Comment("makeNormals") << i );
2688 for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
2690 _LayerEdge& edge = *_sdVec[i]._edges[j];
2691 SMESH_TNodeXYZ nXYZ( edge._nodes[0] );
2692 nXYZ += edge._normal * _sdVec[i]._stepSize;
2693 dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<edge._nodes[0]->GetID()
2694 << ", mesh.AddNode( " << nXYZ.X()<<","<< nXYZ.Y()<<","<< nXYZ.Z()<<")])");
2698 dumpFunction( SMESH_Comment("makeTmpFaces_") << i );
2699 TopExp_Explorer fExp( _sdVec[i]._solid, TopAbs_FACE );
2700 for ( ; fExp.More(); fExp.Next() )
2702 if (const SMESHDS_SubMesh* sm = _sdVec[i]._proxyMesh->GetProxySubMesh( fExp.Current()))
2704 SMDS_ElemIteratorPtr fIt = sm->GetElements();
2705 while ( fIt->more())
2707 const SMDS_MeshElement* e = fIt->next();
2708 SMESH_Comment cmd("mesh.AddFace([");
2709 for ( int j=0; j < e->NbCornerNodes(); ++j )
2710 cmd << e->GetNode(j)->GetID() << (j+1<e->NbCornerNodes() ? ",": "])");
2720 //================================================================================
2722 * \brief Increase length of _LayerEdge's to reach the required thickness of layers
2724 //================================================================================
2726 bool _ViscousBuilder::inflate(_SolidData& data)
2728 SMESH_MesherHelper helper( *_mesh );
2730 // Limit inflation step size by geometry size found by itersecting
2731 // normals of _LayerEdge's with mesh faces
2732 double geomSize = Precision::Infinite(), intersecDist;
2733 auto_ptr<SMESH_ElementSearcher> searcher
2734 ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
2735 data._proxyMesh->GetFaces( data._solid )) );
2736 for ( size_t i = 0; i < data._edges.size(); ++i )
2738 if ( data._edges[i]->IsOnEdge() ) continue;
2739 data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon );
2740 if ( geomSize > intersecDist && intersecDist > 0 )
2741 geomSize = intersecDist;
2743 if ( data._stepSize > 0.3 * geomSize )
2744 limitStepSize( data, 0.3 * geomSize );
2746 const double tgtThick = data._hyp->GetTotalThickness();
2747 if ( data._stepSize > tgtThick )
2748 limitStepSize( data, tgtThick );
2750 if ( data._stepSize < 1. )
2751 data._epsilon = data._stepSize * 1e-7;
2753 debugMsg( "-- geomSize = " << geomSize << ", stepSize = " << data._stepSize );
2755 double avgThick = 0, curThick = 0, distToIntersection = Precision::Infinite();
2756 int nbSteps = 0, nbRepeats = 0;
2757 while ( 1.01 * avgThick < tgtThick )
2759 // new target length
2760 curThick += data._stepSize;
2761 if ( curThick > tgtThick )
2763 curThick = tgtThick + ( tgtThick-avgThick ) * nbRepeats;
2767 // Elongate _LayerEdge's
2768 dumpFunction(SMESH_Comment("inflate")<<data._index<<"_step"<<nbSteps); // debug
2769 for ( size_t i = 0; i < data._edges.size(); ++i )
2771 data._edges[i]->SetNewLength( curThick, helper );
2775 if ( !updateNormals( data, helper, nbSteps ))
2778 // Improve and check quality
2779 if ( !smoothAndCheck( data, nbSteps, distToIntersection ))
2783 dumpFunction(SMESH_Comment("invalidate")<<data._index<<"_step"<<nbSteps); // debug
2784 for ( size_t i = 0; i < data._edges.size(); ++i )
2786 data._edges[i]->InvalidateStep( nbSteps+1 );
2790 break; // no more inflating possible
2794 // Evaluate achieved thickness
2796 for ( size_t i = 0; i < data._edges.size(); ++i )
2797 avgThick += data._edges[i]->_len;
2798 avgThick /= data._edges.size();
2799 debugMsg( "-- Thickness " << avgThick << " reached" );
2801 if ( distToIntersection < avgThick*1.5 )
2803 debugMsg( "-- Stop inflation since "
2804 << " distToIntersection( "<<distToIntersection<<" ) < avgThick( "
2805 << avgThick << " ) * 1.5" );
2809 limitStepSize( data, 0.25 * distToIntersection );
2810 if ( data._stepSizeNodes[0] )
2811 data._stepSize = data._stepSizeCoeff *
2812 SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
2814 } // while ( 1.01 * avgThick < tgtThick )
2817 return error("failed at the very first inflation step", data._index);
2819 if ( 1.01 * avgThick < tgtThick )
2820 if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( data._index ))
2822 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
2823 if ( !smError || smError->IsOK() )
2825 ( new SMESH_ComputeError (COMPERR_WARNING,
2826 SMESH_Comment("Thickness ") << tgtThick <<
2827 " of viscous layers not reached,"
2828 " average reached thickness is " << avgThick ));
2832 // Restore position of src nodes moved by infaltion on _noShrinkShapes
2833 dumpFunction(SMESH_Comment("restoNoShrink_So")<<data._index); // debug
2835 for ( int iS = 0; iS < data._endEdgeOnShape.size(); ++iS )
2838 iEnd = data._endEdgeOnShape[ iS ];
2839 if ( data._edges[ iBeg ]->_nodes.size() == 1 )
2840 for ( ; iBeg < iEnd; ++iBeg )
2842 restoreNoShrink( *data._edges[ iBeg ] );
2850 //================================================================================
2852 * \brief Improve quality of layer inner surface and check intersection
2854 //================================================================================
2856 bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
2858 double & distToIntersection)
2860 if ( data._nbShapesToSmooth == 0 )
2861 return true; // no shapes needing smoothing
2863 bool moved, improved;
2865 SMESH_MesherHelper helper(*_mesh);
2866 Handle(Geom_Surface) surface;
2870 for ( int iS = 0; iS < data._nbShapesToSmooth; ++iS )
2873 iEnd = data._endEdgeOnShape[ iS ];
2875 if ( !data._edges[ iBeg ]->_sWOL.IsNull() &&
2876 data._edges[ iBeg ]->_sWOL.ShapeType() == TopAbs_FACE )
2878 if ( !F.IsSame( data._edges[ iBeg ]->_sWOL )) {
2879 F = TopoDS::Face( data._edges[ iBeg ]->_sWOL );
2880 helper.SetSubShape( F );
2881 surface = BRep_Tool::Surface( F );
2886 F.Nullify(); surface.Nullify();
2888 TGeomID sInd = data._edges[ iBeg ]->_nodes[0]->getshapeId();
2890 if ( data._edges[ iBeg ]->IsOnEdge() )
2892 dumpFunction(SMESH_Comment("smooth")<<data._index << "_Ed"<<sInd <<"_InfStep"<<nbSteps);
2894 // try a simple solution on an analytic EDGE
2895 if ( !smoothAnalyticEdge( data, iBeg, iEnd, surface, F, helper ))
2901 for ( int i = iBeg; i < iEnd; ++i )
2903 moved |= data._edges[i]->SmoothOnEdge(surface, F, helper);
2905 dumpCmd( SMESH_Comment("# end step ")<<step);
2907 while ( moved && step++ < 5 );
2914 int step = 0, stepLimit = 5, badNb = 0; moved = true;
2915 while (( ++step <= stepLimit && moved ) || improved )
2917 dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
2918 <<"_InfStep"<<nbSteps<<"_"<<step); // debug
2919 int oldBadNb = badNb;
2923 for ( int i = iBeg; i < iEnd; ++i )
2924 moved |= data._edges[i]->Smooth(badNb);
2926 for ( int i = iEnd-1; i >= iBeg; --i )
2927 moved |= data._edges[i]->Smooth(badNb);
2928 improved = ( badNb < oldBadNb );
2930 // issue 22576 -- no bad faces but still there are intersections to fix
2931 if ( improved && badNb == 0 )
2932 stepLimit = step + 3;
2939 for ( int i = iBeg; i < iEnd; ++i )
2941 _LayerEdge* edge = data._edges[i];
2942 SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
2943 for ( size_t j = 0; j < edge->_simplices.size(); ++j )
2944 if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
2946 cout << "Bad simplex ( " << edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
2947 << " "<< edge->_simplices[j]._nPrev->GetID()
2948 << " "<< edge->_simplices[j]._nNext->GetID() << " )" << endl;
2956 } // loop on shapes to smooth
2958 // Check orientation of simplices of _ConvexFace::_simplexTestEdges
2959 map< TGeomID, _ConvexFace >::iterator id2face = data._convexFaces.begin();
2960 for ( ; id2face != data._convexFaces.end(); ++id2face )
2962 _ConvexFace & convFace = (*id2face).second;
2963 if ( !convFace._simplexTestEdges.empty() &&
2964 convFace._simplexTestEdges[0]->_nodes[0]->GetPosition()->GetDim() == 2 )
2965 continue; // _simplexTestEdges are based on FACE -- already checked while smoothing
2967 if ( !convFace.CheckPrisms() )
2971 // Check if the last segments of _LayerEdge intersects 2D elements;
2972 // checked elements are either temporary faces or faces on surfaces w/o the layers
2974 auto_ptr<SMESH_ElementSearcher> searcher
2975 ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
2976 data._proxyMesh->GetFaces( data._solid )) );
2978 distToIntersection = Precision::Infinite();
2980 const SMDS_MeshElement* intFace = 0;
2981 const SMDS_MeshElement* closestFace = 0;
2983 for ( size_t i = 0; i < data._edges.size(); ++i )
2985 if ( !data._edges[i]->_sWOL.IsNull() )
2987 if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace ))
2989 if ( distToIntersection > dist )
2991 // ignore intersection of a _LayerEdge based on a _ConvexFace with a face
2992 // lying on this _ConvexFace
2993 if ( _ConvexFace* convFace = data.GetConvexFace( intFace->getshapeId() ))
2994 if ( convFace->_subIdToEdgeEnd.count ( data._edges[i]->_nodes[0]->getshapeId() ))
2997 distToIntersection = dist;
2999 closestFace = intFace;
3005 SMDS_MeshElement::iterator nIt = closestFace->begin_nodes();
3006 cout << "Shortest distance: _LayerEdge nodes: tgt " << data._edges[iLE]->_nodes.back()->GetID()
3007 << " src " << data._edges[iLE]->_nodes[0]->GetID()<< ", intersection with face ("
3008 << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
3009 << ") distance = " << distToIntersection<< endl;
3016 //================================================================================
3018 * \brief Return a curve of the EDGE to be used for smoothing and arrange
3019 * _LayerEdge's to be in a consequent order
3021 //================================================================================
3023 Handle(Geom_Curve) _SolidData::CurveForSmooth( const TopoDS_Edge& E,
3026 Handle(Geom_Surface)& surface,
3027 const TopoDS_Face& F,
3028 SMESH_MesherHelper& helper)
3030 TGeomID eIndex = helper.GetMeshDS()->ShapeToIndex( E );
3032 map< TGeomID, Handle(Geom_Curve)>::iterator i2curve = _edge2curve.find( eIndex );
3034 if ( i2curve == _edge2curve.end() )
3036 // sort _LayerEdge's by position on the EDGE
3037 SortOnEdge( E, iFrom, iTo, helper );
3039 SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( eIndex );
3041 TopLoc_Location loc; double f,l;
3043 Handle(Geom_Line) line;
3044 Handle(Geom_Circle) circle;
3045 bool isLine, isCirc;
3046 if ( F.IsNull() ) // 3D case
3048 // check if the EDGE is a line
3049 Handle(Geom_Curve) curve = BRep_Tool::Curve( E, loc, f, l);
3050 if ( curve->IsKind( STANDARD_TYPE( Geom_TrimmedCurve )))
3051 curve = Handle(Geom_TrimmedCurve)::DownCast( curve )->BasisCurve();
3053 line = Handle(Geom_Line)::DownCast( curve );
3054 circle = Handle(Geom_Circle)::DownCast( curve );
3055 isLine = (!line.IsNull());
3056 isCirc = (!circle.IsNull());
3058 if ( !isLine && !isCirc ) // Check if the EDGE is close to a line
3061 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
3062 while ( nIt->more() )
3063 bndBox.Add( SMESH_TNodeXYZ( nIt->next() ));
3064 gp_XYZ size = bndBox.CornerMax() - bndBox.CornerMin();
3066 SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->tgtNode(0) );
3067 SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->tgtNode(1) );
3068 const double lineTol = 1e-2 * ( p0 - p1 ).Modulus();
3069 for ( int i = 0; i < 3 && !isLine; ++i )
3070 isLine = ( size.Coord( i+1 ) <= lineTol );
3072 if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle
3079 // check if the EDGE is a line
3080 Handle(Geom2d_Curve) curve = BRep_Tool::CurveOnSurface( E, F, f, l);
3081 if ( curve->IsKind( STANDARD_TYPE( Geom2d_TrimmedCurve )))
3082 curve = Handle(Geom2d_TrimmedCurve)::DownCast( curve )->BasisCurve();
3084 Handle(Geom2d_Line) line2d = Handle(Geom2d_Line)::DownCast( curve );
3085 Handle(Geom2d_Circle) circle2d = Handle(Geom2d_Circle)::DownCast( curve );
3086 isLine = (!line2d.IsNull());
3087 isCirc = (!circle2d.IsNull());
3089 if ( !isLine && !isCirc) // Check if the EDGE is close to a line
3092 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
3093 while ( nIt->more() )
3094 bndBox.Add( helper.GetNodeUV( F, nIt->next() ));
3095 gp_XY size = bndBox.CornerMax() - bndBox.CornerMin();
3097 const double lineTol = 1e-2 * sqrt( bndBox.SquareExtent() );
3098 for ( int i = 0; i < 2 && !isLine; ++i )
3099 isLine = ( size.Coord( i+1 ) <= lineTol );
3101 if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle
3107 line = new Geom_Line( gp::OX() ); // only type does matter
3111 gp_Pnt2d p = circle2d->Location();
3112 gp_Ax2 ax( gp_Pnt( p.X(), p.Y(), 0), gp::DX());
3113 circle = new Geom_Circle( ax, 1.); // only center position does matter
3117 Handle(Geom_Curve)& res = _edge2curve[ eIndex ];
3125 return i2curve->second;
3128 //================================================================================
3130 * \brief Sort _LayerEdge's by a parameter on a given EDGE
3132 //================================================================================
3134 void _SolidData::SortOnEdge( const TopoDS_Edge& E,
3137 SMESH_MesherHelper& helper)
3139 map< double, _LayerEdge* > u2edge;
3140 for ( int i = iFrom; i < iTo; ++i )
3141 u2edge.insert( make_pair( helper.GetNodeU( E, _edges[i]->_nodes[0] ), _edges[i] ));
3143 ASSERT( u2edge.size() == iTo - iFrom );
3144 map< double, _LayerEdge* >::iterator u2e = u2edge.begin();
3145 for ( int i = iFrom; i < iTo; ++i, ++u2e )
3146 _edges[i] = u2e->second;
3148 // set _2neibors according to the new order
3149 for ( int i = iFrom; i < iTo-1; ++i )
3150 if ( _edges[i]->_2neibors->tgtNode(1) != _edges[i+1]->_nodes.back() )
3151 _edges[i]->_2neibors->reverse();
3152 if ( u2edge.size() > 1 &&
3153 _edges[iTo-1]->_2neibors->tgtNode(0) != _edges[iTo-2]->_nodes.back() )
3154 _edges[iTo-1]->_2neibors->reverse();
3157 //================================================================================
3159 * \brief Return index corresponding to the shape in _endEdgeOnShape
3161 //================================================================================
3163 bool _SolidData::GetShapeEdges(const TGeomID shapeID,
3168 int beg = 0, end = 0;
3169 for ( edgesEnd = 0; edgesEnd < _endEdgeOnShape.size(); ++edgesEnd )
3171 end = _endEdgeOnShape[ edgesEnd ];
3172 TGeomID sID = _edges[ beg ]->_nodes[0]->getshapeId();
3173 if ( sID == shapeID )
3175 if ( iBeg ) *iBeg = beg;
3176 if ( iEnd ) *iEnd = end;
3184 //================================================================================
3186 * \brief Add faces for smoothing
3188 //================================================================================
3190 void _SolidData::AddShapesToSmooth( const set< TGeomID >& faceIDs )
3192 // convert faceIDs to indices in _endEdgeOnShape
3193 set< size_t > iEnds;
3195 set< TGeomID >::const_iterator fId = faceIDs.begin();
3196 for ( ; fId != faceIDs.end(); ++fId )
3197 if ( GetShapeEdges( *fId, end ) && end >= _nbShapesToSmooth )
3198 iEnds.insert( end );
3200 set< size_t >::iterator endsIt = iEnds.begin();
3202 // "add" by move of _nbShapesToSmooth
3203 int nbFacesToAdd = iEnds.size();
3204 while ( endsIt != iEnds.end() && *endsIt == _nbShapesToSmooth )
3207 ++_nbShapesToSmooth;
3210 if ( endsIt == iEnds.end() )
3213 // Move _LayerEdge's on FACEs just after _nbShapesToSmooth
3215 vector< _LayerEdge* > nonSmoothLE, smoothLE;
3216 size_t lastSmooth = *iEnds.rbegin();
3218 for ( size_t i = _nbShapesToSmooth; i <= lastSmooth; ++i )
3220 vector< _LayerEdge* > & edgesVec = iEnds.count(i) ? smoothLE : nonSmoothLE;
3221 iBeg = i ? _endEdgeOnShape[ i-1 ] : 0;
3222 iEnd = _endEdgeOnShape[ i ];
3223 edgesVec.insert( edgesVec.end(), _edges.begin() + iBeg, _edges.begin() + iEnd );
3226 iBeg = _nbShapesToSmooth ? _endEdgeOnShape[ _nbShapesToSmooth-1 ] : 0;
3227 std::copy( smoothLE.begin(), smoothLE.end(), &_edges[ iBeg ] );
3228 std::copy( nonSmoothLE.begin(), nonSmoothLE.end(), &_edges[ iBeg + smoothLE.size()]);
3230 // update _endEdgeOnShape
3231 for ( size_t i = _nbShapesToSmooth; i < _endEdgeOnShape.size(); ++i )
3233 TGeomID curShape = _edges[ iBeg ]->_nodes[0]->getshapeId();
3234 while ( ++iBeg < _edges.size() &&
3235 curShape == _edges[ iBeg ]->_nodes[0]->getshapeId() );
3237 _endEdgeOnShape[ i ] = iBeg;
3240 _nbShapesToSmooth += nbFacesToAdd;
3243 //================================================================================
3245 * \brief smooth _LayerEdge's on a staight EDGE or circular EDGE
3247 //================================================================================
3249 bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data,
3252 Handle(Geom_Surface)& surface,
3253 const TopoDS_Face& F,
3254 SMESH_MesherHelper& helper)
3256 TopoDS_Shape S = helper.GetSubShapeByNode( data._edges[ iFrom ]->_nodes[0],
3257 helper.GetMeshDS());
3258 TopoDS_Edge E = TopoDS::Edge( S );
3260 Handle(Geom_Curve) curve = data.CurveForSmooth( E, iFrom, iTo, surface, F, helper );
3261 if ( curve.IsNull() ) return false;
3263 // compute a relative length of segments
3264 vector< double > len( iTo-iFrom+1 );
3266 double curLen, prevLen = len[0] = 1.0;
3267 for ( int i = iFrom; i < iTo; ++i )
3269 curLen = prevLen * data._edges[i]->_2neibors->_wgt[0] / data._edges[i]->_2neibors->_wgt[1];
3270 len[i-iFrom+1] = len[i-iFrom] + curLen;
3275 if ( curve->IsKind( STANDARD_TYPE( Geom_Line )))
3277 if ( F.IsNull() ) // 3D
3279 SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->tgtNode(0));
3280 SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->tgtNode(1));
3281 for ( int i = iFrom; i < iTo; ++i )
3283 double r = len[i-iFrom] / len.back();
3284 gp_XYZ newPos = p0 * ( 1. - r ) + p1 * r;
3285 data._edges[i]->_pos.back() = newPos;
3286 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3287 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3288 dumpMove( tgtNode );
3293 // gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->tgtNode(0));
3294 // gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->tgtNode(1));
3295 gp_XY uv0 = data._edges[iFrom]->_2neibors->_edges[0]->LastUV( F );
3296 gp_XY uv1 = data._edges[iTo-1]->_2neibors->_edges[1]->LastUV( F );
3297 if ( data._edges[iFrom]->_2neibors->tgtNode(0) ==
3298 data._edges[iTo-1]->_2neibors->tgtNode(1) ) // closed edge
3300 int iPeriodic = helper.GetPeriodicIndex();
3301 if ( iPeriodic == 1 || iPeriodic == 2 )
3303 uv1.SetCoord( iPeriodic, helper.GetOtherParam( uv1.Coord( iPeriodic )));
3304 if ( uv0.Coord( iPeriodic ) > uv1.Coord( iPeriodic ))
3305 std::swap( uv0, uv1 );
3308 const gp_XY rangeUV = uv1 - uv0;
3309 for ( int i = iFrom; i < iTo; ++i )
3311 double r = len[i-iFrom] / len.back();
3312 gp_XY newUV = uv0 + r * rangeUV;
3313 data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
3315 gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
3316 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3317 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3318 dumpMove( tgtNode );
3320 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
3321 pos->SetUParameter( newUV.X() );
3322 pos->SetVParameter( newUV.Y() );
3328 if ( curve->IsKind( STANDARD_TYPE( Geom_Circle )))
3330 Handle(Geom_Circle) circle = Handle(Geom_Circle)::DownCast( curve );
3331 gp_Pnt center3D = circle->Location();
3333 if ( F.IsNull() ) // 3D
3335 if ( data._edges[iFrom]->_2neibors->tgtNode(0) ==
3336 data._edges[iTo-1]->_2neibors->tgtNode(1) )
3337 return true; // closed EDGE - nothing to do
3339 return false; // TODO ???
3343 const gp_XY center( center3D.X(), center3D.Y() );
3345 gp_XY uv0 = data._edges[iFrom]->_2neibors->_edges[0]->LastUV( F );
3346 gp_XY uvM = data._edges[iFrom]->LastUV( F );
3347 gp_XY uv1 = data._edges[iTo-1]->_2neibors->_edges[1]->LastUV( F );
3348 // gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->tgtNode(0));
3349 // gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back());
3350 // gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->tgtNode(1));
3351 gp_Vec2d vec0( center, uv0 );
3352 gp_Vec2d vecM( center, uvM );
3353 gp_Vec2d vec1( center, uv1 );
3354 double uLast = vec0.Angle( vec1 ); // -PI - +PI
3355 double uMidl = vec0.Angle( vecM );
3356 if ( uLast * uMidl <= 0. )
3357 uLast += ( uMidl > 0 ? +2. : -2. ) * M_PI;
3358 const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() );
3360 gp_Ax2d axis( center, vec0 );
3361 gp_Circ2d circ( axis, radius );
3362 for ( int i = iFrom; i < iTo; ++i )
3364 double newU = uLast * len[i-iFrom] / len.back();
3365 gp_Pnt2d newUV = ElCLib::Value( newU, circ );
3366 data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
3368 gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
3369 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3370 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3371 dumpMove( tgtNode );
3373 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
3374 pos->SetUParameter( newUV.X() );
3375 pos->SetVParameter( newUV.Y() );
3384 //================================================================================
3386 * \brief Modify normals of _LayerEdge's on EDGE's to avoid intersection with
3387 * _LayerEdge's on neighbor EDGE's
3389 //================================================================================
3391 bool _ViscousBuilder::updateNormals( _SolidData& data,
3392 SMESH_MesherHelper& helper,
3396 return updateNormalsOfConvexFaces( data, helper, stepNb );
3398 // make temporary quadrangles got by extrusion of
3399 // mesh edges along _LayerEdge._normal's
3401 vector< const SMDS_MeshElement* > tmpFaces;
3403 set< SMESH_TLink > extrudedLinks; // contains target nodes
3404 vector< const SMDS_MeshNode*> nodes(4); // of a tmp mesh face
3406 dumpFunction(SMESH_Comment("makeTmpFacesOnEdges")<<data._index);
3407 for ( size_t i = 0; i < data._edges.size(); ++i )
3409 _LayerEdge* edge = data._edges[i];
3410 if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
3411 const SMDS_MeshNode* tgt1 = edge->_nodes.back();
3412 for ( int j = 0; j < 2; ++j ) // loop on _2NearEdges
3414 const SMDS_MeshNode* tgt2 = edge->_2neibors->tgtNode(j);
3415 pair< set< SMESH_TLink >::iterator, bool > link_isnew =
3416 extrudedLinks.insert( SMESH_TLink( tgt1, tgt2 ));
3417 if ( !link_isnew.second )
3419 extrudedLinks.erase( link_isnew.first );
3420 continue; // already extruded and will no more encounter
3422 // a _LayerEdge containg tgt2
3423 _LayerEdge* neiborEdge = edge->_2neibors->_edges[j];
3425 _TmpMeshFaceOnEdge* f = new _TmpMeshFaceOnEdge( edge, neiborEdge, --_tmpFaceID );
3426 tmpFaces.push_back( f );
3428 dumpCmd(SMESH_Comment("mesh.AddFace([ ")
3429 <<f->_nn[0]->GetID()<<", "<<f->_nn[1]->GetID()<<", "
3430 <<f->_nn[2]->GetID()<<", "<<f->_nn[3]->GetID()<<" ])");
3435 // Check if _LayerEdge's based on EDGE's intersects tmpFaces.
3436 // Perform two loops on _LayerEdge on EDGE's:
3437 // 1) to find and fix intersection
3438 // 2) to check that no new intersection appears as result of 1)
3440 SMDS_ElemIteratorPtr fIt( new SMDS_ElementVectorIterator( tmpFaces.begin(),
3442 auto_ptr<SMESH_ElementSearcher> searcher
3443 ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), fIt ));
3445 // 1) Find intersections
3447 const SMDS_MeshElement* face;
3448 typedef map< _LayerEdge*, set< _LayerEdge*, _LayerEdgeCmp >, _LayerEdgeCmp > TLEdge2LEdgeSet;
3449 TLEdge2LEdgeSet edge2CloseEdge;
3451 const double eps = data._epsilon * data._epsilon;
3452 for ( size_t i = 0; i < data._edges.size(); ++i )
3454 _LayerEdge* edge = data._edges[i];
3455 if (( !edge->IsOnEdge() ) &&
3456 ( edge->_sWOL.IsNull() || edge->_sWOL.ShapeType() != TopAbs_FACE ))
3458 if ( edge->FindIntersection( *searcher, dist, eps, &face ))
3460 const _TmpMeshFaceOnEdge* f = (const _TmpMeshFaceOnEdge*) face;
3461 set< _LayerEdge*, _LayerEdgeCmp > & ee = edge2CloseEdge[ edge ];
3462 ee.insert( f->_le1 );
3463 ee.insert( f->_le2 );
3464 if ( f->_le1->IsOnEdge() && f->_le1->_sWOL.IsNull() )
3465 edge2CloseEdge[ f->_le1 ].insert( edge );
3466 if ( f->_le2->IsOnEdge() && f->_le2->_sWOL.IsNull() )
3467 edge2CloseEdge[ f->_le2 ].insert( edge );
3471 // Set _LayerEdge._normal
3473 if ( !edge2CloseEdge.empty() )
3475 dumpFunction(SMESH_Comment("updateNormals")<<data._index);
3477 set< TGeomID > shapesToSmooth;
3479 // vector to store new _normal and _cosin for each edge in edge2CloseEdge
3480 vector< pair< _LayerEdge*, _LayerEdge > > edge2newEdge( edge2CloseEdge.size() );
3482 TLEdge2LEdgeSet::iterator e2ee = edge2CloseEdge.begin();
3483 for ( size_t iE = 0; e2ee != edge2CloseEdge.end(); ++e2ee, ++iE )
3485 _LayerEdge* edge1 = e2ee->first;
3486 _LayerEdge* edge2 = 0;
3487 set< _LayerEdge*, _LayerEdgeCmp >& ee = e2ee->second;
3489 edge2newEdge[ iE ].first = NULL;
3491 // find EDGEs the edges reside
3492 // TopoDS_Edge E1, E2;
3493 // TopoDS_Shape S = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
3494 // if ( S.ShapeType() != TopAbs_EDGE )
3495 // continue; // TODO: find EDGE by VERTEX
3496 // E1 = TopoDS::Edge( S );
3497 set< _LayerEdge*, _LayerEdgeCmp >::iterator eIt = ee.begin();
3498 for ( ; !edge2 && eIt != ee.end(); ++eIt )
3500 if ( edge1->_sWOL == (*eIt)->_sWOL )
3503 if ( !edge2 ) continue;
3505 edge2newEdge[ iE ].first = edge1;
3506 _LayerEdge& newEdge = edge2newEdge[ iE ].second;
3507 // while ( E2.IsNull() && eIt != ee.end())
3509 // _LayerEdge* e2 = *eIt++;
3510 // TopoDS_Shape S = helper.GetSubShapeByNode( e2->_nodes[0], getMeshDS() );
3511 // if ( S.ShapeType() == TopAbs_EDGE )
3512 // E2 = TopoDS::Edge( S ), edge2 = e2;
3514 // if ( E2.IsNull() ) continue; // TODO: find EDGE by VERTEX
3516 // find 3 FACEs sharing 2 EDGEs
3518 // TopoDS_Face FF1[2], FF2[2];
3519 // PShapeIteratorPtr fIt = helper.GetAncestors(E1, *_mesh, TopAbs_FACE);
3520 // while ( fIt->more() && FF1[1].IsNull() )
3522 // const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
3523 // if ( helper.IsSubShape( *F, data._solid))
3524 // FF1[ FF1[0].IsNull() ? 0 : 1 ] = *F;
3526 // fIt = helper.GetAncestors(E2, *_mesh, TopAbs_FACE);
3527 // while ( fIt->more() && FF2[1].IsNull())
3529 // const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
3530 // if ( helper.IsSubShape( *F, data._solid))
3531 // FF2[ FF2[0].IsNull() ? 0 : 1 ] = *F;
3533 // // exclude a FACE common to E1 and E2 (put it to FFn[1] )
3534 // if ( FF1[0].IsSame( FF2[0]) || FF1[0].IsSame( FF2[1]))
3535 // std::swap( FF1[0], FF1[1] );
3536 // if ( FF2[0].IsSame( FF1[0]) )
3537 // std::swap( FF2[0], FF2[1] );
3538 // if ( FF1[0].IsNull() || FF2[0].IsNull() )
3541 // get a new normal for edge1
3543 gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal;
3544 // if ( edge1->_cosin < 0 )
3545 // dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized();
3546 // if ( edge2->_cosin < 0 )
3547 // dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized();
3549 double cos1 = Abs( edge1->_cosin ), cos2 = Abs( edge2->_cosin );
3550 double wgt1 = ( cos1 + 0.001 ) / ( cos1 + cos2 + 0.002 );
3551 double wgt2 = ( cos2 + 0.001 ) / ( cos1 + cos2 + 0.002 );
3552 newEdge._normal = ( wgt1 * dir1 + wgt2 * dir2 ).XYZ();
3553 newEdge._normal.Normalize();
3555 // cout << edge1->_nodes[0]->GetID() << " "
3556 // << edge2->_nodes[0]->GetID() << " NORM: "
3557 // << newEdge._normal.X() << ", " << newEdge._normal.Y() << ", " << newEdge._normal.Z() << endl;
3560 if ( cos1 < theMinSmoothCosin )
3562 newEdge._cosin = edge2->_cosin;
3564 else if ( cos2 > theMinSmoothCosin ) // both cos1 and cos2 > theMinSmoothCosin
3566 // gp_Vec dirInFace;
3567 // if ( edge1->_cosin < 0 )
3568 // dirInFace = dir1;
3570 // dirInFace = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
3571 // double angle = dirInFace.Angle( edge1->_normal ); // [0,PI]
3572 // edge1->SetCosin( Cos( angle ));
3573 //newEdge._cosin = 0; // ???????????
3574 newEdge._cosin = ( wgt1 * cos1 + wgt2 * cos2 ) * edge1->_cosin / cos1;
3578 newEdge._cosin = edge1->_cosin;
3581 // find shapes that need smoothing due to change of _normal
3582 if ( edge1->_cosin < theMinSmoothCosin &&
3583 newEdge._cosin > theMinSmoothCosin )
3585 if ( edge1->_sWOL.IsNull() )
3587 SMDS_ElemIteratorPtr fIt = edge1->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
3588 while ( fIt->more() )
3589 shapesToSmooth.insert( fIt->next()->getshapeId() );
3590 //limitStepSize( data, fIt->next(), edge1->_cosin ); // too late
3592 else // edge1 inflates along a FACE
3594 TopoDS_Shape V = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
3595 PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE );
3596 while ( const TopoDS_Shape* E = eIt->next() )
3598 if ( !helper.IsSubShape( *E, /*FACE=*/edge1->_sWOL ))
3600 gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V ));
3601 double angle = edgeDir.Angle( newEdge._normal ); // [0,PI]
3602 if ( angle < M_PI / 2 )
3603 shapesToSmooth.insert( getMeshDS()->ShapeToIndex( *E ));
3609 data.AddShapesToSmooth( shapesToSmooth );
3611 // Update data of edges depending on a new _normal
3613 for ( size_t iE = 0; iE < edge2newEdge.size(); ++iE )
3615 _LayerEdge* edge1 = edge2newEdge[ iE ].first;
3616 _LayerEdge& newEdge = edge2newEdge[ iE ].second;
3617 if ( !edge1 ) continue;
3619 edge1->_normal = newEdge._normal;
3620 edge1->SetCosin( newEdge._cosin );
3621 edge1->InvalidateStep( 1 );
3623 edge1->SetNewLength( data._stepSize, helper );
3624 if ( edge1->IsOnEdge() )
3626 const SMDS_MeshNode * n1 = edge1->_2neibors->srcNode(0);
3627 const SMDS_MeshNode * n2 = edge1->_2neibors->srcNode(1);
3628 edge1->SetDataByNeighbors( n1, n2, helper );
3631 // Update normals and other dependent data of not intersecting _LayerEdge's
3632 // neighboring the intersecting ones
3634 if ( !edge1->_2neibors )
3636 for ( int j = 0; j < 2; ++j ) // loop on 2 neighbors
3638 _LayerEdge* neighbor = edge1->_2neibors->_edges[j];
3639 if ( edge2CloseEdge.count ( neighbor ))
3640 continue; // j-th neighbor is also intersected
3641 _LayerEdge* prevEdge = edge1;
3642 const int nbSteps = 10;
3643 for ( int step = nbSteps; step; --step ) // step from edge1 in j-th direction
3645 if ( !neighbor->_2neibors )
3646 break; // neighbor is on VERTEX
3648 _LayerEdge* nextEdge = neighbor->_2neibors->_edges[iNext];
3649 if ( nextEdge == prevEdge )
3650 nextEdge = neighbor->_2neibors->_edges[ ++iNext ];
3651 double r = double(step-1)/nbSteps;
3652 if ( !nextEdge->_2neibors )
3655 gp_XYZ newNorm = prevEdge->_normal * r + nextEdge->_normal * (1-r);
3656 newNorm.Normalize();
3658 neighbor->_normal = newNorm;
3659 neighbor->SetCosin( prevEdge->_cosin * r + nextEdge->_cosin * (1-r) );
3660 neighbor->SetDataByNeighbors( prevEdge->_nodes[0], nextEdge->_nodes[0], helper );
3662 neighbor->InvalidateStep( 1 );
3664 neighbor->SetNewLength( data._stepSize, helper );
3666 // goto the next neighbor
3667 prevEdge = neighbor;
3668 neighbor = nextEdge;
3674 // 2) Check absence of intersections
3677 for ( size_t i = 0 ; i < tmpFaces.size(); ++i )
3683 //================================================================================
3685 * \brief Modify normals of _LayerEdge's on _ConvexFace's
3687 //================================================================================
3689 bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData& data,
3690 SMESH_MesherHelper& helper,
3693 SMESHDS_Mesh* meshDS = helper.GetMeshDS();
3696 map< TGeomID, _ConvexFace >::iterator id2face = data._convexFaces.begin();
3697 for ( ; id2face != data._convexFaces.end(); ++id2face )
3699 _ConvexFace & convFace = (*id2face).second;
3700 if ( convFace._normalsFixed )
3701 continue; // already fixed
3702 if ( convFace.CheckPrisms() )
3703 continue; // nothing to fix
3705 convFace._normalsFixed = true;
3707 BRepAdaptor_Surface surface ( convFace._face, false );
3708 BRepLProp_SLProps surfProp( surface, 2, 1e-6 );
3710 // check if the convex FACE is of spherical shape
3712 Bnd_B3d centersBox; // bbox of centers of curvature of _LayerEdge's on VERTEXes
3717 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.begin();
3718 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3720 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3722 if ( meshDS->IndexToShape( id2end->first ).ShapeType() == TopAbs_VERTEX )
3724 _LayerEdge* ledge = data._edges[ iBeg ];
3725 if ( convFace.GetCenterOfCurvature( ledge, surfProp, helper, center ))
3726 centersBox.Add( center );
3728 for ( ; iBeg < iEnd; ++iBeg )
3729 nodesBox.Add( SMESH_TNodeXYZ( data._edges[ iBeg ]->_nodes[0] ));
3731 if ( centersBox.IsVoid() )
3733 debugMsg( "Error: centersBox.IsVoid()" );
3736 const bool isSpherical =
3737 ( centersBox.SquareExtent() < 1e-6 * nodesBox.SquareExtent() );
3739 int nbEdges = helper.Count( convFace._face, TopAbs_EDGE, /*ignoreSame=*/false );
3740 vector < _CentralCurveOnEdge > centerCurves( nbEdges );
3744 // set _LayerEdge::_normal as average of all normals
3746 // WARNING: different density of nodes on EDGEs is not taken into account that
3747 // can lead to an improper new normal
3749 gp_XYZ avgNormal( 0,0,0 );
3751 id2end = convFace._subIdToEdgeEnd.begin();
3752 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3754 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3755 // set data of _CentralCurveOnEdge
3756 const TopoDS_Shape& S = meshDS->IndexToShape( id2end->first );
3757 if ( S.ShapeType() == TopAbs_EDGE )
3759 _CentralCurveOnEdge& ceCurve = centerCurves[ nbEdges++ ];
3760 ceCurve.SetShapes( TopoDS::Edge(S), convFace, data, helper );
3761 if ( !data._edges[ iBeg ]->_sWOL.IsNull() )
3762 ceCurve._adjFace.Nullify();
3764 ceCurve._ledges.insert( ceCurve._ledges.end(),
3765 &data._edges[ iBeg ], &data._edges[ iEnd ]);
3767 // summarize normals
3768 for ( ; iBeg < iEnd; ++iBeg )
3769 avgNormal += data._edges[ iBeg ]->_normal;
3771 double normSize = avgNormal.SquareModulus();
3772 if ( normSize < 1e-200 )
3774 debugMsg( "updateNormalsOfConvexFaces(): zero avgNormal" );
3777 avgNormal /= Sqrt( normSize );
3779 // compute new _LayerEdge::_cosin on EDGEs
3780 double avgCosin = 0;
3783 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3785 _CentralCurveOnEdge& ceCurve = centerCurves[ iE ];
3786 if ( ceCurve._adjFace.IsNull() )
3788 for ( size_t iLE = 0; iLE < ceCurve._ledges.size(); ++iLE )
3790 const SMDS_MeshNode* node = ceCurve._ledges[ iLE ]->_nodes[0];
3791 inFaceDir = getFaceDir( ceCurve._adjFace, ceCurve._edge, node, helper, isOK );
3794 double angle = inFaceDir.Angle( avgNormal ); // [0,PI]
3795 ceCurve._ledges[ iLE ]->_cosin = Cos( angle );
3796 avgCosin += ceCurve._ledges[ iLE ]->_cosin;
3802 avgCosin /= nbCosin;
3804 // set _LayerEdge::_normal = avgNormal
3805 id2end = convFace._subIdToEdgeEnd.begin();
3806 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3808 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3809 const TopoDS_Shape& S = meshDS->IndexToShape( id2end->first );
3810 if ( S.ShapeType() != TopAbs_EDGE )
3811 for ( int i = iBeg; i < iEnd; ++i )
3812 data._edges[ i ]->_cosin = avgCosin;
3814 for ( ; iBeg < iEnd; ++iBeg )
3815 data._edges[ iBeg ]->_normal = avgNormal;
3818 else // if ( isSpherical )
3820 // We suppose that centers of curvature at all points of the FACE
3821 // lie on some curve, let's call it "central curve". For all _LayerEdge's
3822 // having a common center of curvature we define the same new normal
3823 // as a sum of normals of _LayerEdge's on EDGEs among them.
3825 // get all centers of curvature for each EDGE
3827 helper.SetSubShape( convFace._face );
3828 _LayerEdge* vertexLEdges[2], **edgeLEdge, **edgeLEdgeEnd;
3830 TopExp_Explorer edgeExp( convFace._face, TopAbs_EDGE );
3831 for ( int iE = 0; edgeExp.More(); edgeExp.Next(), ++iE )
3833 const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
3835 // set adjacent FACE
3836 centerCurves[ iE ].SetShapes( edge, convFace, data, helper );
3838 // get _LayerEdge's of the EDGE
3839 TGeomID edgeID = meshDS->ShapeToIndex( edge );
3840 id2end = convFace._subIdToEdgeEnd.find( edgeID );
3841 if ( id2end == convFace._subIdToEdgeEnd.end() )
3843 // no _LayerEdge's on EDGE, use _LayerEdge's on VERTEXes
3844 for ( int iV = 0; iV < 2; ++iV )
3846 TopoDS_Vertex v = helper.IthVertex( iV, edge );
3847 TGeomID vID = meshDS->ShapeToIndex( v );
3848 int end = convFace._subIdToEdgeEnd[ vID ];
3849 int iBeg = end > 0 ? data._endEdgeOnShape[ end-1 ] : 0;
3850 vertexLEdges[ iV ] = data._edges[ iBeg ];
3852 edgeLEdge = &vertexLEdges[0];
3853 edgeLEdgeEnd = edgeLEdge + 2;
3855 centerCurves[ iE ]._adjFace.Nullify();
3859 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3860 if ( id2end->second >= data._nbShapesToSmooth )
3861 data.SortOnEdge( edge, iBeg, iEnd, helper );
3862 edgeLEdge = &data._edges[ iBeg ];
3863 edgeLEdgeEnd = edgeLEdge + iEnd - iBeg;
3864 vertexLEdges[0] = data._edges[ iBeg ]->_2neibors->_edges[0];
3865 vertexLEdges[1] = data._edges[ iEnd-1 ]->_2neibors->_edges[1];
3867 if ( ! data._edges[ iBeg ]->_sWOL.IsNull() )
3868 centerCurves[ iE ]._adjFace.Nullify();
3871 // Get curvature centers
3875 if ( edgeLEdge[0]->IsOnEdge() &&
3876 convFace.GetCenterOfCurvature( vertexLEdges[0], surfProp, helper, center ))
3878 centerCurves[ iE ].Append( center, vertexLEdges[0] );
3879 centersBox.Add( center );
3881 for ( ; edgeLEdge < edgeLEdgeEnd; ++edgeLEdge )
3882 if ( convFace.GetCenterOfCurvature( *edgeLEdge, surfProp, helper, center ))
3883 { // EDGE or VERTEXes
3884 centerCurves[ iE ].Append( center, *edgeLEdge );
3885 centersBox.Add( center );
3887 if ( edgeLEdge[-1]->IsOnEdge() &&
3888 convFace.GetCenterOfCurvature( vertexLEdges[1], surfProp, helper, center ))
3890 centerCurves[ iE ].Append( center, vertexLEdges[1] );
3891 centersBox.Add( center );
3893 centerCurves[ iE ]._isDegenerated =
3894 ( centersBox.IsVoid() || centersBox.SquareExtent() < 1e-6 * nodesBox.SquareExtent() );
3896 } // loop on EDGES of convFace._face to set up data of centerCurves
3898 // Compute new normals for _LayerEdge's on EDGEs
3900 double avgCosin = 0;
3903 for ( size_t iE1 = 0; iE1 < centerCurves.size(); ++iE1 )
3905 _CentralCurveOnEdge& ceCurve = centerCurves[ iE1 ];
3906 if ( ceCurve._isDegenerated )
3908 const vector< gp_Pnt >& centers = ceCurve._curvaCenters;
3909 vector< gp_XYZ > & newNormals = ceCurve._normals;
3910 for ( size_t iC1 = 0; iC1 < centers.size(); ++iC1 )
3913 for ( size_t iE2 = 0; iE2 < centerCurves.size() && !isOK; ++iE2 )
3916 isOK = centerCurves[ iE2 ].FindNewNormal( centers[ iC1 ], newNormals[ iC1 ]);
3918 if ( isOK && !ceCurve._adjFace.IsNull() )
3920 // compute new _LayerEdge::_cosin
3921 const SMDS_MeshNode* node = ceCurve._ledges[ iC1 ]->_nodes[0];
3922 inFaceDir = getFaceDir( ceCurve._adjFace, ceCurve._edge, node, helper, isOK );
3925 double angle = inFaceDir.Angle( newNormals[ iC1 ] ); // [0,PI]
3926 ceCurve._ledges[ iC1 ]->_cosin = Cos( angle );
3927 avgCosin += ceCurve._ledges[ iC1 ]->_cosin;
3933 // set new normals to _LayerEdge's of NOT degenerated central curves
3934 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3936 if ( centerCurves[ iE ]._isDegenerated )
3938 for ( size_t iLE = 0; iLE < centerCurves[ iE ]._ledges.size(); ++iLE )
3939 centerCurves[ iE ]._ledges[ iLE ]->_normal = centerCurves[ iE ]._normals[ iLE ];
3941 // set new normals to _LayerEdge's of degenerated central curves
3942 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3944 if ( !centerCurves[ iE ]._isDegenerated ||
3945 centerCurves[ iE ]._ledges.size() < 3 )
3947 // new normal is an average of new normals at VERTEXes that
3948 // was computed on non-degenerated _CentralCurveOnEdge's
3949 gp_XYZ newNorm = ( centerCurves[ iE ]._ledges.front()->_normal +
3950 centerCurves[ iE ]._ledges.back ()->_normal );
3951 double sz = newNorm.Modulus();
3955 double newCosin = ( 0.5 * centerCurves[ iE ]._ledges.front()->_cosin +
3956 0.5 * centerCurves[ iE ]._ledges.back ()->_cosin );
3957 for ( size_t iLE = 1, nb = centerCurves[ iE ]._ledges.size() - 1; iLE < nb; ++iLE )
3959 centerCurves[ iE ]._ledges[ iLE ]->_normal = newNorm;
3960 centerCurves[ iE ]._ledges[ iLE ]->_cosin = newCosin;
3964 // Find new normals for _LayerEdge's based on FACE
3967 avgCosin /= nbCosin;
3968 const TGeomID faceID = meshDS->ShapeToIndex( convFace._face );
3969 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.find( faceID );
3970 if ( id2end != convFace._subIdToEdgeEnd.end() )
3974 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3975 for ( ; iBeg < iEnd; ++iBeg )
3977 _LayerEdge* ledge = data._edges[ iBeg ];
3978 if ( !convFace.GetCenterOfCurvature( ledge, surfProp, helper, center ))
3980 for ( size_t i = 0; i < centerCurves.size(); ++i, ++iE )
3982 iE = iE % centerCurves.size();
3983 if ( centerCurves[ iE ]._isDegenerated )
3985 newNorm.SetCoord( 0,0,0 );
3986 if ( centerCurves[ iE ].FindNewNormal( center, newNorm ))
3988 ledge->_normal = newNorm;
3989 ledge->_cosin = avgCosin;
3996 } // not a quasi-spherical FACE
3998 // Update _LayerEdge's data according to a new normal
4000 dumpFunction(SMESH_Comment("updateNormalsOfConvexFaces")<<data._index
4001 <<"_F"<<meshDS->ShapeToIndex( convFace._face ));
4003 id2end = convFace._subIdToEdgeEnd.begin();
4004 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
4006 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
4007 for ( ; iBeg < iEnd; ++iBeg )
4009 _LayerEdge* & ledge = data._edges[ iBeg ];
4010 double len = ledge->_len;
4011 ledge->InvalidateStep( stepNb + 1, /*restoreLength=*/true );
4012 ledge->SetCosin( ledge->_cosin );
4013 ledge->SetNewLength( len, helper );
4016 } // loop on sub-shapes of convFace._face
4018 // Find FACEs adjacent to convFace._face that got necessity to smooth
4019 // as a result of normals modification
4021 set< TGeomID > adjFacesToSmooth;
4022 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
4024 if ( centerCurves[ iE ]._adjFace.IsNull() ||
4025 centerCurves[ iE ]._adjFaceToSmooth )
4027 for ( size_t iLE = 0; iLE < centerCurves[ iE ]._ledges.size(); ++iLE )
4029 if ( centerCurves[ iE ]._ledges[ iLE ]->_cosin > theMinSmoothCosin )
4031 adjFacesToSmooth.insert( meshDS->ShapeToIndex( centerCurves[ iE ]._adjFace ));
4036 data.AddShapesToSmooth( adjFacesToSmooth );
4041 } // loop on data._convexFaces
4046 //================================================================================
4048 * \brief Finds a center of curvature of a surface at a _LayerEdge
4050 //================================================================================
4052 bool _ConvexFace::GetCenterOfCurvature( _LayerEdge* ledge,
4053 BRepLProp_SLProps& surfProp,
4054 SMESH_MesherHelper& helper,
4055 gp_Pnt & center ) const
4057 gp_XY uv = helper.GetNodeUV( _face, ledge->_nodes[0] );
4058 surfProp.SetParameters( uv.X(), uv.Y() );
4059 if ( !surfProp.IsCurvatureDefined() )
4062 const double oriFactor = ( _face.Orientation() == TopAbs_REVERSED ? +1. : -1. );
4063 double surfCurvatureMax = surfProp.MaxCurvature() * oriFactor;
4064 double surfCurvatureMin = surfProp.MinCurvature() * oriFactor;
4065 if ( surfCurvatureMin > surfCurvatureMax )
4066 center = surfProp.Value().Translated( surfProp.Normal().XYZ() / surfCurvatureMin * oriFactor );
4068 center = surfProp.Value().Translated( surfProp.Normal().XYZ() / surfCurvatureMax * oriFactor );
4073 //================================================================================
4075 * \brief Check that prisms are not distorted
4077 //================================================================================
4079 bool _ConvexFace::CheckPrisms() const
4081 for ( size_t i = 0; i < _simplexTestEdges.size(); ++i )
4083 const _LayerEdge* edge = _simplexTestEdges[i];
4084 SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
4085 for ( size_t j = 0; j < edge->_simplices.size(); ++j )
4086 if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
4088 debugMsg( "Bad simplex of _simplexTestEdges ("
4089 << " "<< edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
4090 << " "<< edge->_simplices[j]._nPrev->GetID()
4091 << " "<< edge->_simplices[j]._nNext->GetID() << " )" );
4098 //================================================================================
4100 * \brief Try to compute a new normal by interpolating normals of _LayerEdge's
4101 * stored in this _CentralCurveOnEdge.
4102 * \param [in] center - curvature center of a point of another _CentralCurveOnEdge.
4103 * \param [in,out] newNormal - current normal at this point, to be redefined
4104 * \return bool - true if succeeded.
4106 //================================================================================
4108 bool _CentralCurveOnEdge::FindNewNormal( const gp_Pnt& center, gp_XYZ& newNormal )
4110 if ( this->_isDegenerated )
4113 // find two centers the given one lies between
4115 for ( size_t i = 0, nb = _curvaCenters.size()-1; i < nb; ++i )
4117 double sl2 = 1.001 * _segLength2[ i ];
4119 double d1 = center.SquareDistance( _curvaCenters[ i ]);
4123 double d2 = center.SquareDistance( _curvaCenters[ i+1 ]);
4124 if ( d2 > sl2 || d2 + d1 < 1e-100 )
4129 double r = d1 / ( d1 + d2 );
4130 gp_XYZ norm = (( 1. - r ) * _ledges[ i ]->_normal +
4131 ( r ) * _ledges[ i+1 ]->_normal );
4135 double sz = newNormal.Modulus();
4144 //================================================================================
4146 * \brief Set shape members
4148 //================================================================================
4150 void _CentralCurveOnEdge::SetShapes( const TopoDS_Edge& edge,
4151 const _ConvexFace& convFace,
4152 const _SolidData& data,
4153 SMESH_MesherHelper& helper)
4157 PShapeIteratorPtr fIt = helper.GetAncestors( edge, *helper.GetMesh(), TopAbs_FACE );
4158 while ( const TopoDS_Shape* F = fIt->next())
4159 if ( !convFace._face.IsSame( *F ))
4161 _adjFace = TopoDS::Face( *F );
4162 _adjFaceToSmooth = false;
4163 // _adjFace already in a smoothing queue ?
4165 TGeomID adjFaceID = helper.GetMeshDS()->ShapeToIndex( *F );
4166 if ( data.GetShapeEdges( adjFaceID, end ))
4167 _adjFaceToSmooth = ( end < data._nbShapesToSmooth );
4172 //================================================================================
4174 * \brief Looks for intersection of it's last segment with faces
4175 * \param distance - returns shortest distance from the last node to intersection
4177 //================================================================================
4179 bool _LayerEdge::FindIntersection( SMESH_ElementSearcher& searcher,
4181 const double& epsilon,
4182 const SMDS_MeshElement** face)
4184 vector< const SMDS_MeshElement* > suspectFaces;
4186 gp_Ax1 lastSegment = LastSegment(segLen);
4187 searcher.GetElementsNearLine( lastSegment, SMDSAbs_Face, suspectFaces );
4189 bool segmentIntersected = false;
4190 distance = Precision::Infinite();
4191 int iFace = -1; // intersected face
4192 for ( size_t j = 0 ; j < suspectFaces.size() /*&& !segmentIntersected*/; ++j )
4194 const SMDS_MeshElement* face = suspectFaces[j];
4195 if ( face->GetNodeIndex( _nodes.back() ) >= 0 ||
4196 face->GetNodeIndex( _nodes[0] ) >= 0 )
4197 continue; // face sharing _LayerEdge node
4198 const int nbNodes = face->NbCornerNodes();
4199 bool intFound = false;
4201 SMDS_MeshElement::iterator nIt = face->begin_nodes();
4204 intFound = SegTriaInter( lastSegment, *nIt++, *nIt++, *nIt++, dist, epsilon );
4208 const SMDS_MeshNode* tria[3];
4211 for ( int n2 = 2; n2 < nbNodes && !intFound; ++n2 )
4214 intFound = SegTriaInter(lastSegment, tria[0], tria[1], tria[2], dist, epsilon );
4220 if ( dist < segLen*(1.01) && dist > -(_len*_lenFactor-segLen) )
4221 segmentIntersected = true;
4222 if ( distance > dist )
4223 distance = dist, iFace = j;
4226 if ( iFace != -1 && face ) *face = suspectFaces[iFace];
4228 if ( segmentIntersected )
4231 SMDS_MeshElement::iterator nIt = suspectFaces[iFace]->begin_nodes();
4232 gp_XYZ intP( lastSegment.Location().XYZ() + lastSegment.Direction().XYZ() * distance );
4233 cout << "nodes: tgt " << _nodes.back()->GetID() << " src " << _nodes[0]->GetID()
4234 << ", intersection with face ("
4235 << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
4236 << ") at point (" << intP.X() << ", " << intP.Y() << ", " << intP.Z()
4237 << ") distance = " << distance - segLen<< endl;
4243 return segmentIntersected;
4246 //================================================================================
4248 * \brief Returns size and direction of the last segment
4250 //================================================================================
4252 gp_Ax1 _LayerEdge::LastSegment(double& segLen) const
4254 // find two non-coincident positions
4255 gp_XYZ orig = _pos.back();
4257 int iPrev = _pos.size() - 2;
4258 while ( iPrev >= 0 )
4260 dir = orig - _pos[iPrev];
4261 if ( dir.SquareModulus() > 1e-100 )
4271 segDir.SetLocation( SMESH_TNodeXYZ( _nodes[0] ));
4272 segDir.SetDirection( _normal );
4277 gp_Pnt pPrev = _pos[ iPrev ];
4278 if ( !_sWOL.IsNull() )
4280 TopLoc_Location loc;
4281 if ( _sWOL.ShapeType() == TopAbs_EDGE )
4284 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
4285 pPrev = curve->Value( pPrev.X() ).Transformed( loc );
4289 Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
4290 pPrev = surface->Value( pPrev.X(), pPrev.Y() ).Transformed( loc );
4292 dir = SMESH_TNodeXYZ( _nodes.back() ) - pPrev.XYZ();
4294 segDir.SetLocation( pPrev );
4295 segDir.SetDirection( dir );
4296 segLen = dir.Modulus();
4302 //================================================================================
4304 * \brief Return the last position of the target node on a FACE.
4305 * \param [in] F - the FACE this _LayerEdge is inflated along
4306 * \return gp_XY - result UV
4308 //================================================================================
4310 gp_XY _LayerEdge::LastUV( const TopoDS_Face& F ) const
4312 if ( F.IsSame( _sWOL )) // F is my FACE
4313 return gp_XY( _pos.back().X(), _pos.back().Y() );
4315 if ( _sWOL.IsNull() || _sWOL.ShapeType() != TopAbs_EDGE ) // wrong call
4316 return gp_XY( 1e100, 1e100 );
4318 // _sWOL is EDGE of F; _pos.back().X() is the last U on the EDGE
4319 double f, l, u = _pos.back().X();
4320 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( TopoDS::Edge(_sWOL), F, f,l);
4321 if ( !C2d.IsNull() && f <= u && u <= l )
4322 return C2d->Value( u ).XY();
4324 return gp_XY( 1e100, 1e100 );
4327 //================================================================================
4329 * \brief Test intersection of the last segment with a given triangle
4330 * using Moller-Trumbore algorithm
4331 * Intersection is detected if distance to intersection is less than _LayerEdge._len
4333 //================================================================================
4335 bool _LayerEdge::SegTriaInter( const gp_Ax1& lastSegment,
4336 const SMDS_MeshNode* n0,
4337 const SMDS_MeshNode* n1,
4338 const SMDS_MeshNode* n2,
4340 const double& EPSILON) const
4342 //const double EPSILON = 1e-6;
4344 gp_XYZ orig = lastSegment.Location().XYZ();
4345 gp_XYZ dir = lastSegment.Direction().XYZ();
4347 SMESH_TNodeXYZ vert0( n0 );
4348 SMESH_TNodeXYZ vert1( n1 );
4349 SMESH_TNodeXYZ vert2( n2 );
4351 /* calculate distance from vert0 to ray origin */
4352 gp_XYZ tvec = orig - vert0;
4354 //if ( tvec * dir > EPSILON )
4355 // intersected face is at back side of the temporary face this _LayerEdge belongs to
4358 gp_XYZ edge1 = vert1 - vert0;
4359 gp_XYZ edge2 = vert2 - vert0;
4361 /* begin calculating determinant - also used to calculate U parameter */
4362 gp_XYZ pvec = dir ^ edge2;
4364 /* if determinant is near zero, ray lies in plane of triangle */
4365 double det = edge1 * pvec;
4367 if (det > -EPSILON && det < EPSILON)
4369 double inv_det = 1.0 / det;
4371 /* calculate U parameter and test bounds */
4372 double u = ( tvec * pvec ) * inv_det;
4373 //if (u < 0.0 || u > 1.0)
4374 if (u < -EPSILON || u > 1.0 + EPSILON)
4377 /* prepare to test V parameter */
4378 gp_XYZ qvec = tvec ^ edge1;
4380 /* calculate V parameter and test bounds */
4381 double v = (dir * qvec) * inv_det;
4382 //if ( v < 0.0 || u + v > 1.0 )
4383 if ( v < -EPSILON || u + v > 1.0 + EPSILON)
4386 /* calculate t, ray intersects triangle */
4387 t = (edge2 * qvec) * inv_det;
4393 //================================================================================
4395 * \brief Perform smooth of _LayerEdge's based on EDGE's
4396 * \retval bool - true if node has been moved
4398 //================================================================================
4400 bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface,
4401 const TopoDS_Face& F,
4402 SMESH_MesherHelper& helper)
4404 ASSERT( IsOnEdge() );
4406 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( _nodes.back() );
4407 SMESH_TNodeXYZ oldPos( tgtNode );
4408 double dist01, distNewOld;
4410 SMESH_TNodeXYZ p0( _2neibors->tgtNode(0));
4411 SMESH_TNodeXYZ p1( _2neibors->tgtNode(1));
4412 dist01 = p0.Distance( _2neibors->tgtNode(1) );
4414 gp_Pnt newPos = p0 * _2neibors->_wgt[0] + p1 * _2neibors->_wgt[1];
4415 double lenDelta = 0;
4418 //lenDelta = _curvature->lenDelta( _len );
4419 lenDelta = _curvature->lenDeltaByDist( dist01 );
4420 newPos.ChangeCoord() += _normal * lenDelta;
4423 distNewOld = newPos.Distance( oldPos );
4427 if ( _2neibors->_plnNorm )
4429 // put newPos on the plane defined by source node and _plnNorm
4430 gp_XYZ new2src = SMESH_TNodeXYZ( _nodes[0] ) - newPos.XYZ();
4431 double new2srcProj = (*_2neibors->_plnNorm) * new2src;
4432 newPos.ChangeCoord() += (*_2neibors->_plnNorm) * new2srcProj;
4434 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4435 _pos.back() = newPos.XYZ();
4439 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4440 gp_XY uv( Precision::Infinite(), 0 );
4441 helper.CheckNodeUV( F, tgtNode, uv, 1e-10, /*force=*/true );
4442 _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
4444 newPos = surface->Value( uv.X(), uv.Y() );
4445 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4448 if ( _curvature && lenDelta < 0 )
4450 gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
4451 _len -= prevPos.Distance( oldPos );
4452 _len += prevPos.Distance( newPos );
4454 bool moved = distNewOld > dist01/50;
4456 dumpMove( tgtNode ); // debug
4461 //================================================================================
4463 * \brief Perform laplacian smooth in 3D of nodes inflated from FACE
4464 * \retval bool - true if _tgtNode has been moved
4466 //================================================================================
4468 bool _LayerEdge::Smooth(int& badNb)
4470 if ( _simplices.size() < 2 )
4471 return false; // _LayerEdge inflated along EDGE or FACE
4473 // compute new position for the last _pos
4474 gp_XYZ newPos (0,0,0);
4475 for ( size_t i = 0; i < _simplices.size(); ++i )
4476 newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev );
4477 newPos /= _simplices.size();
4479 const gp_XYZ& curPos ( _pos.back() );
4480 const gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
4483 double delta = _curvature->lenDelta( _len );
4485 newPos += _normal * delta;
4488 double segLen = _normal * ( newPos - prevPos.XYZ() );
4489 if ( segLen + delta > 0 )
4490 newPos += _normal * delta;
4492 // double segLenChange = _normal * ( curPos - newPos );
4493 // newPos += 0.5 * _normal * segLenChange;
4496 // count quality metrics (orientation) of tetras around _tgtNode
4498 for ( size_t i = 0; i < _simplices.size(); ++i )
4499 nbOkBefore += _simplices[i].IsForward( _nodes[0], &curPos );
4502 for ( size_t i = 0; i < _simplices.size(); ++i )
4503 nbOkAfter += _simplices[i].IsForward( _nodes[0], &newPos );
4505 if ( nbOkAfter < nbOkBefore )
4508 SMDS_MeshNode* n = const_cast< SMDS_MeshNode* >( _nodes.back() );
4510 _len -= prevPos.Distance(SMESH_TNodeXYZ( n ));
4511 _len += prevPos.Distance(newPos);
4513 n->setXYZ( newPos.X(), newPos.Y(), newPos.Z());
4514 _pos.back() = newPos;
4516 badNb += _simplices.size() - nbOkAfter;
4523 //================================================================================
4525 * \brief Add a new segment to _LayerEdge during inflation
4527 //================================================================================
4529 void _LayerEdge::SetNewLength( double len, SMESH_MesherHelper& helper )
4531 if ( _len - len > -1e-6 )
4533 _pos.push_back( _pos.back() );
4537 SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
4538 SMESH_TNodeXYZ oldXYZ( n );
4539 gp_XYZ nXYZ = oldXYZ + _normal * ( len - _len ) * _lenFactor;
4540 n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
4542 _pos.push_back( nXYZ );
4544 if ( !_sWOL.IsNull() )
4547 if ( _sWOL.ShapeType() == TopAbs_EDGE )
4549 double u = Precision::Infinite(); // to force projection w/o distance check
4550 helper.CheckNodeU( TopoDS::Edge( _sWOL ), n, u, 1e-10, /*force=*/true, distXYZ );
4551 _pos.back().SetCoord( u, 0, 0 );
4552 if ( _nodes.size() > 1 )
4554 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
4555 pos->SetUParameter( u );
4560 gp_XY uv( Precision::Infinite(), 0 );
4561 helper.CheckNodeUV( TopoDS::Face( _sWOL ), n, uv, 1e-10, /*force=*/true, distXYZ );
4562 _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
4563 if ( _nodes.size() > 1 )
4565 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
4566 pos->SetUParameter( uv.X() );
4567 pos->SetVParameter( uv.Y() );
4570 n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]);
4572 dumpMove( n ); //debug
4575 //================================================================================
4577 * \brief Remove last inflation step
4579 //================================================================================
4581 void _LayerEdge::InvalidateStep( int curStep, bool restoreLength )
4583 if ( _pos.size() > curStep )
4585 if ( restoreLength )
4586 _len -= ( _pos[ curStep-1 ] - _pos.back() ).Modulus();
4588 _pos.resize( curStep );
4589 gp_Pnt nXYZ = _pos.back();
4590 SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
4591 if ( !_sWOL.IsNull() )
4593 TopLoc_Location loc;
4594 if ( _sWOL.ShapeType() == TopAbs_EDGE )
4596 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
4597 pos->SetUParameter( nXYZ.X() );
4599 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
4600 nXYZ = curve->Value( nXYZ.X() ).Transformed( loc );
4604 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
4605 pos->SetUParameter( nXYZ.X() );
4606 pos->SetVParameter( nXYZ.Y() );
4607 Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
4608 nXYZ = surface->Value( nXYZ.X(), nXYZ.Y() ).Transformed( loc );
4611 n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
4616 //================================================================================
4618 * \brief Create layers of prisms
4620 //================================================================================
4622 bool _ViscousBuilder::refine(_SolidData& data)
4624 SMESH_MesherHelper helper( *_mesh );
4625 helper.SetSubShape( data._solid );
4626 helper.SetElementsOnShape(false);
4628 Handle(Geom_Curve) curve;
4629 Handle(Geom_Surface) surface;
4630 TopoDS_Edge geomEdge;
4631 TopoDS_Face geomFace;
4632 TopoDS_Shape prevSWOL;
4633 TopLoc_Location loc;
4637 TGeomID prevBaseId = -1;
4638 TNode2Edge* n2eMap = 0;
4639 TNode2Edge::iterator n2e;
4641 // Create intermediate nodes on each _LayerEdge
4643 for ( size_t i = 0; i < data._edges.size(); ++i )
4645 _LayerEdge& edge = *data._edges[i];
4647 if ( edge._nodes.size() < 2 )
4648 continue; // on _noShrinkShapes
4650 // get accumulated length of segments
4651 vector< double > segLen( edge._pos.size() );
4653 for ( size_t j = 1; j < edge._pos.size(); ++j )
4654 segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus();
4656 // allocate memory for new nodes if it is not yet refined
4657 const SMDS_MeshNode* tgtNode = edge._nodes.back();
4658 if ( edge._nodes.size() == 2 )
4660 edge._nodes.resize( data._hyp->GetNumberLayers() + 1, 0 );
4662 edge._nodes.back() = tgtNode;
4664 // get data of a shrink shape
4665 if ( !edge._sWOL.IsNull() && edge._sWOL != prevSWOL )
4667 isOnEdge = ( edge._sWOL.ShapeType() == TopAbs_EDGE );
4670 geomEdge = TopoDS::Edge( edge._sWOL );
4671 curve = BRep_Tool::Curve( geomEdge, loc, f,l);
4675 geomFace = TopoDS::Face( edge._sWOL );
4676 surface = BRep_Tool::Surface( geomFace, loc );
4678 prevSWOL = edge._sWOL;
4680 // restore shapePos of the last node by already treated _LayerEdge of another _SolidData
4681 const TGeomID baseShapeId = edge._nodes[0]->getshapeId();
4682 if ( baseShapeId != prevBaseId )
4684 map< TGeomID, TNode2Edge* >::iterator s2ne = data._s2neMap.find( baseShapeId );
4685 n2eMap = ( s2ne == data._s2neMap.end() ) ? 0 : n2eMap = s2ne->second;
4686 prevBaseId = baseShapeId;
4688 _LayerEdge* edgeOnSameNode = 0;
4689 if ( n2eMap && (( n2e = n2eMap->find( edge._nodes[0] )) != n2eMap->end() ))
4691 edgeOnSameNode = n2e->second;
4692 const gp_XYZ& otherTgtPos = edgeOnSameNode->_pos.back();
4693 SMDS_PositionPtr lastPos = tgtNode->GetPosition();
4696 SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( lastPos );
4697 epos->SetUParameter( otherTgtPos.X() );
4701 SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( lastPos );
4702 fpos->SetUParameter( otherTgtPos.X() );
4703 fpos->SetVParameter( otherTgtPos.Y() );
4706 // calculate height of the first layer
4708 const double T = segLen.back(); //data._hyp.GetTotalThickness();
4709 const double f = data._hyp->GetStretchFactor();
4710 const int N = data._hyp->GetNumberLayers();
4711 const double fPowN = pow( f, N );
4712 if ( fPowN - 1 <= numeric_limits<double>::min() )
4715 h0 = T * ( f - 1 )/( fPowN - 1 );
4717 const double zeroLen = std::numeric_limits<double>::min();
4719 // create intermediate nodes
4720 double hSum = 0, hi = h0/f;
4722 for ( size_t iStep = 1; iStep < edge._nodes.size(); ++iStep )
4724 // compute an intermediate position
4727 while ( hSum > segLen[iSeg] && iSeg < segLen.size()-1)
4729 int iPrevSeg = iSeg-1;
4730 while ( fabs( segLen[iPrevSeg] - segLen[iSeg]) <= zeroLen && iPrevSeg > 0 )
4732 double r = ( segLen[iSeg] - hSum ) / ( segLen[iSeg] - segLen[iPrevSeg] );
4733 gp_Pnt pos = r * edge._pos[iPrevSeg] + (1-r) * edge._pos[iSeg];
4735 SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >( edge._nodes[ iStep ]);
4736 if ( !edge._sWOL.IsNull() )
4738 // compute XYZ by parameters <pos>
4743 pos = curve->Value( u ).Transformed(loc);
4747 uv.SetCoord( pos.X(), pos.Y() );
4749 pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc);
4752 // create or update the node
4755 node = helper.AddNode( pos.X(), pos.Y(), pos.Z());
4756 if ( !edge._sWOL.IsNull() )
4759 getMeshDS()->SetNodeOnEdge( node, geomEdge, u );
4761 getMeshDS()->SetNodeOnFace( node, geomFace, uv.X(), uv.Y() );
4765 getMeshDS()->SetNodeInVolume( node, helper.GetSubShapeID() );
4770 if ( !edge._sWOL.IsNull() )
4772 // make average pos from new and current parameters
4775 u = 0.5 * ( u + helper.GetNodeU( geomEdge, node ));
4776 pos = curve->Value( u ).Transformed(loc);
4778 SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( node->GetPosition() );
4779 epos->SetUParameter( u );
4783 uv = 0.5 * ( uv + helper.GetNodeUV( geomFace, node ));
4784 pos = surface->Value( uv.X(), uv.Y()).Transformed(loc);
4786 SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( node->GetPosition() );
4787 fpos->SetUParameter( uv.X() );
4788 fpos->SetVParameter( uv.Y() );
4791 node->setXYZ( pos.X(), pos.Y(), pos.Z() );
4793 } // loop on edge._nodes
4795 if ( !edge._sWOL.IsNull() ) // prepare for shrink()
4798 edge._pos.back().SetCoord( u, 0,0);
4800 edge._pos.back().SetCoord( uv.X(), uv.Y() ,0);
4802 if ( edgeOnSameNode )
4803 edgeOnSameNode->_pos.back() = edge._pos.back();
4806 } // loop on data._edges to create nodes
4808 if ( !getMeshDS()->IsEmbeddedMode() )
4809 // Log node movement
4810 for ( size_t i = 0; i < data._edges.size(); ++i )
4812 _LayerEdge& edge = *data._edges[i];
4813 SMESH_TNodeXYZ p ( edge._nodes.back() );
4814 getMeshDS()->MoveNode( p._node, p.X(), p.Y(), p.Z() );
4819 helper.SetElementsOnShape(true);
4821 vector< vector<const SMDS_MeshNode*>* > nnVec;
4822 set< int > degenEdgeInd;
4823 vector<const SMDS_MeshElement*> degenVols;
4825 TopExp_Explorer exp( data._solid, TopAbs_FACE );
4826 for ( ; exp.More(); exp.Next() )
4828 if ( data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
4830 SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
4831 SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
4832 while ( fIt->more() )
4834 const SMDS_MeshElement* face = fIt->next();
4835 const int nbNodes = face->NbCornerNodes();
4836 nnVec.resize( nbNodes );
4837 degenEdgeInd.clear();
4839 SMDS_ElemIteratorPtr nIt = face->nodesIterator();
4840 for ( int iN = 0; iN < nbNodes; ++iN )
4842 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
4843 nnVec[ iN ] = & data._n2eMap[ n ]->_nodes;
4844 if ( nnVec[ iN ]->size() < 2 )
4845 degenEdgeInd.insert( iN );
4847 nbZ = nnVec[ iN ]->size();
4855 switch ( degenEdgeInd.size() )
4859 for ( int iZ = 1; iZ < nbZ; ++iZ )
4860 helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
4861 (*nnVec[0])[iZ], (*nnVec[1])[iZ], (*nnVec[2])[iZ]);
4866 int i2 = *degenEdgeInd.begin();
4867 int i0 = helper.WrapIndex( i2 - 1, nbNodes );
4868 int i1 = helper.WrapIndex( i2 + 1, nbNodes );
4869 for ( int iZ = 1; iZ < nbZ; ++iZ )
4870 helper.AddVolume( (*nnVec[i0])[iZ-1], (*nnVec[i1])[iZ-1],
4871 (*nnVec[i1])[iZ], (*nnVec[i0])[iZ], (*nnVec[i2])[0]);
4876 int i3 = !degenEdgeInd.count(0) ? 0 : !degenEdgeInd.count(1) ? 1 : 2;
4877 for ( int iZ = 1; iZ < nbZ; ++iZ )
4878 helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
4886 switch ( degenEdgeInd.size() )
4890 for ( int iZ = 1; iZ < nbZ; ++iZ )
4891 helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1],
4892 (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1],
4893 (*nnVec[0])[iZ], (*nnVec[1])[iZ],
4894 (*nnVec[2])[iZ], (*nnVec[3])[iZ]);
4899 int i2 = *degenEdgeInd.begin();
4900 int i3 = *degenEdgeInd.rbegin();
4901 bool ok = ( i3 - i2 == 1 );
4902 if ( i2 == 0 && i3 == 3 ) { i2 = 3; i3 = 0; ok = true; }
4903 int i0 = helper.WrapIndex( i3 + 1, nbNodes );
4904 int i1 = helper.WrapIndex( i0 + 1, nbNodes );
4905 for ( int iZ = 1; iZ < nbZ; ++iZ )
4907 const SMDS_MeshElement* vol =
4908 helper.AddVolume( (*nnVec[i3])[0], (*nnVec[i0])[iZ], (*nnVec[i0])[iZ-1],
4909 (*nnVec[i2])[0], (*nnVec[i1])[iZ], (*nnVec[i1])[iZ-1]);
4911 degenVols.push_back( vol );
4915 case 3: // degen HEX
4917 const SMDS_MeshNode* nn[8];
4918 for ( int iZ = 1; iZ < nbZ; ++iZ )
4920 const SMDS_MeshElement* vol =
4921 helper.AddVolume( nnVec[0]->size() > 1 ? (*nnVec[0])[iZ-1] : (*nnVec[0])[0],
4922 nnVec[1]->size() > 1 ? (*nnVec[1])[iZ-1] : (*nnVec[1])[0],
4923 nnVec[2]->size() > 1 ? (*nnVec[2])[iZ-1] : (*nnVec[2])[0],
4924 nnVec[3]->size() > 1 ? (*nnVec[3])[iZ-1] : (*nnVec[3])[0],
4925 nnVec[0]->size() > 1 ? (*nnVec[0])[iZ] : (*nnVec[0])[0],
4926 nnVec[1]->size() > 1 ? (*nnVec[1])[iZ] : (*nnVec[1])[0],
4927 nnVec[2]->size() > 1 ? (*nnVec[2])[iZ] : (*nnVec[2])[0],
4928 nnVec[3]->size() > 1 ? (*nnVec[3])[iZ] : (*nnVec[3])[0]);
4929 degenVols.push_back( vol );
4937 return error("Not supported type of element", data._index);
4939 } // switch ( nbNodes )
4940 } // while ( fIt->more() )
4943 if ( !degenVols.empty() )
4945 SMESH_ComputeErrorPtr& err = _mesh->GetSubMesh( data._solid )->GetComputeError();
4946 if ( !err || err->IsOK() )
4948 err.reset( new SMESH_ComputeError( COMPERR_WARNING,
4949 "Degenerated volumes created" ));
4950 err->myBadElements.insert( err->myBadElements.end(),
4951 degenVols.begin(),degenVols.end() );
4958 //================================================================================
4960 * \brief Shrink 2D mesh on faces to let space for inflated layers
4962 //================================================================================
4964 bool _ViscousBuilder::shrink()
4966 // make map of (ids of FACEs to shrink mesh on) to (_SolidData containing _LayerEdge's
4967 // inflated along FACE or EDGE)
4968 map< TGeomID, _SolidData* > f2sdMap;
4969 for ( size_t i = 0 ; i < _sdVec.size(); ++i )
4971 _SolidData& data = _sdVec[i];
4972 TopTools_MapOfShape FFMap;
4973 map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
4974 for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
4975 if ( s2s->second.ShapeType() == TopAbs_FACE )
4977 f2sdMap.insert( make_pair( getMeshDS()->ShapeToIndex( s2s->second ), &data ));
4979 if ( FFMap.Add( (*s2s).second ))
4980 // Put mesh faces on the shrinked FACE to the proxy sub-mesh to avoid
4981 // usage of mesh faces made in addBoundaryElements() by the 3D algo or
4982 // by StdMeshers_QuadToTriaAdaptor
4983 if ( SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( s2s->second ))
4985 SMESH_ProxyMesh::SubMesh* proxySub =
4986 data._proxyMesh->getFaceSubM( TopoDS::Face( s2s->second ), /*create=*/true);
4987 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
4988 while ( fIt->more() )
4989 proxySub->AddElement( fIt->next() );
4990 // as a result 3D algo will use elements from proxySub and not from smDS
4995 SMESH_MesherHelper helper( *_mesh );
4996 helper.ToFixNodeParameters( true );
4999 map< TGeomID, _Shrinker1D > e2shrMap;
5000 vector< _LayerEdge* > lEdges;
5002 // loop on FACES to srink mesh on
5003 map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin();
5004 for ( ; f2sd != f2sdMap.end(); ++f2sd )
5006 _SolidData& data = *f2sd->second;
5007 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first ));
5008 SMESH_subMesh* sm = _mesh->GetSubMesh( F );
5009 SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
5011 Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
5013 helper.SetSubShape(F);
5015 // ===========================
5016 // Prepare data for shrinking
5017 // ===========================
5019 // Collect nodes to smooth, as src nodes are not yet replaced by tgt ones
5020 // and thus all nodes on a FACE connected to 2d elements are to be smoothed
5021 vector < const SMDS_MeshNode* > smoothNodes;
5023 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
5024 while ( nIt->more() )
5026 const SMDS_MeshNode* n = nIt->next();
5027 if ( n->NbInverseElements( SMDSAbs_Face ) > 0 )
5028 smoothNodes.push_back( n );
5031 // Find out face orientation
5033 const set<TGeomID> ignoreShapes;
5035 if ( !smoothNodes.empty() )
5037 vector<_Simplex> simplices;
5038 getSimplices( smoothNodes[0], simplices, ignoreShapes );
5039 helper.GetNodeUV( F, simplices[0]._nPrev, 0, &isOkUV ); // fix UV of silpmex nodes
5040 helper.GetNodeUV( F, simplices[0]._nNext, 0, &isOkUV );
5041 gp_XY uv = helper.GetNodeUV( F, smoothNodes[0], 0, &isOkUV );
5042 if ( !simplices[0].IsForward(uv, smoothNodes[0], F, helper,refSign) )
5046 // Find _LayerEdge's inflated along F
5049 set< TGeomID > subIDs;
5050 SMESH_subMeshIteratorPtr subIt = sm->getDependsOnIterator(/*includeSelf=*/false);
5051 while ( subIt->more() )
5052 subIDs.insert( subIt->next()->GetId() );
5055 for ( int iS = 0; iS < data._endEdgeOnShape.size() && !subIDs.empty(); ++iS )
5058 iEnd = data._endEdgeOnShape[ iS ];
5059 TGeomID shapeID = data._edges[ iBeg ]->_nodes[0]->getshapeId();
5060 set< TGeomID >::iterator idIt = subIDs.find( shapeID );
5061 if ( idIt == subIDs.end() ||
5062 data._edges[ iBeg ]->_sWOL.IsNull() ) continue;
5063 subIDs.erase( idIt );
5065 if ( !data._noShrinkShapes.count( shapeID ))
5066 for ( ; iBeg < iEnd; ++iBeg )
5068 lEdges.push_back( data._edges[ iBeg ] );
5069 prepareEdgeToShrink( *data._edges[ iBeg ], F, helper, smDS );
5074 dumpFunction(SMESH_Comment("beforeShrinkFace")<<f2sd->first); // debug
5075 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
5076 while ( fIt->more() )
5077 if ( const SMDS_MeshElement* f = fIt->next() )
5078 dumpChangeNodes( f );
5080 // Replace source nodes by target nodes in mesh faces to shrink
5081 dumpFunction(SMESH_Comment("replNodesOnFace")<<f2sd->first); // debug
5082 const SMDS_MeshNode* nodes[20];
5083 for ( size_t i = 0; i < lEdges.size(); ++i )
5085 _LayerEdge& edge = *lEdges[i];
5086 const SMDS_MeshNode* srcNode = edge._nodes[0];
5087 const SMDS_MeshNode* tgtNode = edge._nodes.back();
5088 SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
5089 while ( fIt->more() )
5091 const SMDS_MeshElement* f = fIt->next();
5092 if ( !smDS->Contains( f ))
5094 SMDS_NodeIteratorPtr nIt = f->nodeIterator();
5095 for ( int iN = 0; nIt->more(); ++iN )
5097 const SMDS_MeshNode* n = nIt->next();
5098 nodes[iN] = ( n == srcNode ? tgtNode : n );
5100 helper.GetMeshDS()->ChangeElementNodes( f, nodes, f->NbNodes() );
5101 dumpChangeNodes( f );
5105 // find out if a FACE is concave
5106 const bool isConcaveFace = isConcave( F, helper );
5108 // Create _SmoothNode's on face F
5109 vector< _SmoothNode > nodesToSmooth( smoothNodes.size() );
5111 dumpFunction(SMESH_Comment("fixUVOnFace")<<f2sd->first); // debug
5112 const bool sortSimplices = isConcaveFace;
5113 for ( size_t i = 0; i < smoothNodes.size(); ++i )
5115 const SMDS_MeshNode* n = smoothNodes[i];
5116 nodesToSmooth[ i ]._node = n;
5117 // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices
5118 getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, sortSimplices );
5119 // fix up incorrect uv of nodes on the FACE
5120 helper.GetNodeUV( F, n, 0, &isOkUV);
5124 //if ( nodesToSmooth.empty() ) continue;
5126 // Find EDGE's to shrink and set simpices to LayerEdge's
5127 set< _Shrinker1D* > eShri1D;
5129 for ( size_t i = 0; i < lEdges.size(); ++i )
5131 _LayerEdge* edge = lEdges[i];
5132 if ( edge->_sWOL.ShapeType() == TopAbs_EDGE )
5134 TGeomID edgeIndex = getMeshDS()->ShapeToIndex( edge->_sWOL );
5135 _Shrinker1D& srinker = e2shrMap[ edgeIndex ];
5136 eShri1D.insert( & srinker );
5137 srinker.AddEdge( edge, helper );
5138 VISCOUS_3D::ToClearSubWithMain( _mesh->GetSubMesh( edge->_sWOL ), data._solid );
5139 // restore params of nodes on EGDE if the EDGE has been already
5140 // srinked while srinking another FACE
5141 srinker.RestoreParams();
5143 getSimplices( /*tgtNode=*/edge->_nodes.back(), edge->_simplices, ignoreShapes );
5147 bool toFixTria = false; // to improve quality of trias by diagonal swap
5148 if ( isConcaveFace )
5150 const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles();
5151 if ( hasTria != hasQuad ) {
5152 toFixTria = hasTria;
5155 set<int> nbNodesSet;
5156 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
5157 while ( fIt->more() && nbNodesSet.size() < 2 )
5158 nbNodesSet.insert( fIt->next()->NbCornerNodes() );
5159 toFixTria = ( *nbNodesSet.begin() == 3 );
5163 // ==================
5164 // Perform shrinking
5165 // ==================
5167 bool shrinked = true;
5168 int badNb, shriStep=0, smooStep=0;
5169 _SmoothNode::SmoothType smoothType
5170 = isConcaveFace ? _SmoothNode::ANGULAR : _SmoothNode::LAPLACIAN;
5174 // Move boundary nodes (actually just set new UV)
5175 // -----------------------------------------------
5176 dumpFunction(SMESH_Comment("moveBoundaryOnF")<<f2sd->first<<"_st"<<shriStep ); // debug
5178 for ( size_t i = 0; i < lEdges.size(); ++i )
5180 shrinked |= lEdges[i]->SetNewLength2d( surface,F,helper );
5184 // Move nodes on EDGE's
5185 // (XYZ is set as soon as a needed length reached in SetNewLength2d())
5186 set< _Shrinker1D* >::iterator shr = eShri1D.begin();
5187 for ( ; shr != eShri1D.end(); ++shr )
5188 (*shr)->Compute( /*set3D=*/false, helper );
5191 // -----------------
5192 int nbNoImpSteps = 0;
5195 while (( nbNoImpSteps < 5 && badNb > 0) && moved)
5197 dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
5199 int oldBadNb = badNb;
5202 for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5204 moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
5205 smoothType, /*set3D=*/isConcaveFace);
5207 if ( badNb < oldBadNb )
5215 return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first );
5216 if ( shriStep > 200 )
5217 return error(SMESH_Comment("Infinite loop at shrinking 2D mesh on face ") << f2sd->first );
5219 // Fix narrow triangles by swapping diagonals
5220 // ---------------------------------------
5223 set<const SMDS_MeshNode*> usedNodes;
5224 fixBadFaces( F, helper, /*is2D=*/true, shriStep, & usedNodes); // swap diagonals
5226 // update working data
5227 set<const SMDS_MeshNode*>::iterator n;
5228 for ( size_t i = 0; i < nodesToSmooth.size() && !usedNodes.empty(); ++i )
5230 n = usedNodes.find( nodesToSmooth[ i ]._node );
5231 if ( n != usedNodes.end())
5233 getSimplices( nodesToSmooth[ i ]._node,
5234 nodesToSmooth[ i ]._simplices,
5236 /*sortSimplices=*/ smoothType == _SmoothNode::ANGULAR );
5237 usedNodes.erase( n );
5240 for ( size_t i = 0; i < lEdges.size() && !usedNodes.empty(); ++i )
5242 n = usedNodes.find( /*tgtNode=*/ lEdges[i]->_nodes.back() );
5243 if ( n != usedNodes.end())
5245 getSimplices( lEdges[i]->_nodes.back(),
5246 lEdges[i]->_simplices,
5248 usedNodes.erase( n );
5252 // TODO: check effect of this additional smooth
5253 // additional laplacian smooth to increase allowed shrink step
5254 // for ( int st = 1; st; --st )
5256 // dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
5257 // for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5259 // nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
5260 // _SmoothNode::LAPLACIAN,/*set3D=*/false);
5263 } // while ( shrinked )
5265 // No wrongly shaped faces remain; final smooth. Set node XYZ.
5266 bool isStructuredFixed = false;
5267 if ( SMESH_2D_Algo* algo = dynamic_cast<SMESH_2D_Algo*>( sm->GetAlgo() ))
5268 isStructuredFixed = algo->FixInternalNodes( *data._proxyMesh, F );
5269 if ( !isStructuredFixed )
5271 if ( isConcaveFace ) // fix narrow faces by swapping diagonals
5272 fixBadFaces( F, helper, /*is2D=*/false, ++shriStep );
5274 for ( int st = 3; st; --st )
5277 case 1: smoothType = _SmoothNode::LAPLACIAN; break;
5278 case 2: smoothType = _SmoothNode::LAPLACIAN; break;
5279 case 3: smoothType = _SmoothNode::ANGULAR; break;
5281 dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
5282 for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5284 nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
5285 smoothType,/*set3D=*/st==1 );
5290 // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh
5291 VISCOUS_3D::ToClearSubWithMain( sm, data._solid );
5293 if ( !getMeshDS()->IsEmbeddedMode() )
5294 // Log node movement
5295 for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5297 SMESH_TNodeXYZ p ( nodesToSmooth[i]._node );
5298 getMeshDS()->MoveNode( nodesToSmooth[i]._node, p.X(), p.Y(), p.Z() );
5301 } // loop on FACES to srink mesh on
5304 // Replace source nodes by target nodes in shrinked mesh edges
5306 map< int, _Shrinker1D >::iterator e2shr = e2shrMap.begin();
5307 for ( ; e2shr != e2shrMap.end(); ++e2shr )
5308 e2shr->second.SwapSrcTgtNodes( getMeshDS() );
5313 //================================================================================
5315 * \brief Computes 2d shrink direction and finds nodes limiting shrinking
5317 //================================================================================
5319 bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge,
5320 const TopoDS_Face& F,
5321 SMESH_MesherHelper& helper,
5322 const SMESHDS_SubMesh* faceSubMesh)
5324 const SMDS_MeshNode* srcNode = edge._nodes[0];
5325 const SMDS_MeshNode* tgtNode = edge._nodes.back();
5327 if ( edge._sWOL.ShapeType() == TopAbs_FACE )
5329 gp_XY srcUV( edge._pos[0].X(), edge._pos[0].Y() );//helper.GetNodeUV( F, srcNode );
5330 gp_XY tgtUV = edge.LastUV( F ); //helper.GetNodeUV( F, tgtNode );
5331 gp_Vec2d uvDir( srcUV, tgtUV );
5332 double uvLen = uvDir.Magnitude();
5334 edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0 );
5337 edge._pos.resize(1);
5338 edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 );
5340 // set UV of source node to target node
5341 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
5342 pos->SetUParameter( srcUV.X() );
5343 pos->SetVParameter( srcUV.Y() );
5345 else // _sWOL is TopAbs_EDGE
5347 const TopoDS_Edge& E = TopoDS::Edge( edge._sWOL );
5348 SMESHDS_SubMesh* edgeSM = getMeshDS()->MeshElements( E );
5349 if ( !edgeSM || edgeSM->NbElements() == 0 )
5350 return error(SMESH_Comment("Not meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
5352 const SMDS_MeshNode* n2 = 0;
5353 SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
5354 while ( eIt->more() && !n2 )
5356 const SMDS_MeshElement* e = eIt->next();
5357 if ( !edgeSM->Contains(e)) continue;
5358 n2 = e->GetNode( 0 );
5359 if ( n2 == srcNode ) n2 = e->GetNode( 1 );
5362 return error(SMESH_Comment("Wrongly meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
5364 double uSrc = helper.GetNodeU( E, srcNode, n2 );
5365 double uTgt = helper.GetNodeU( E, tgtNode, srcNode );
5366 double u2 = helper.GetNodeU( E, n2, srcNode );
5370 if ( fabs( uSrc-uTgt ) < 0.99 * fabs( uSrc-u2 ))
5372 // tgtNode is located so that it does not make faces with wrong orientation
5375 edge._pos.resize(1);
5376 edge._pos[0].SetCoord( U_TGT, uTgt );
5377 edge._pos[0].SetCoord( U_SRC, uSrc );
5378 edge._pos[0].SetCoord( LEN_TGT, fabs( uSrc-uTgt ));
5380 edge._simplices.resize( 1 );
5381 edge._simplices[0]._nPrev = n2;
5383 // set U of source node to the target node
5384 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
5385 pos->SetUParameter( uSrc );
5390 //================================================================================
5392 * \brief Restore position of a sole node of a _LayerEdge based on _noShrinkShapes
5394 //================================================================================
5396 void _ViscousBuilder::restoreNoShrink( _LayerEdge& edge ) const
5398 if ( edge._nodes.size() == 1 )
5403 const SMDS_MeshNode* srcNode = edge._nodes[0];
5404 TopoDS_Shape S = SMESH_MesherHelper::GetSubShapeByNode( srcNode, getMeshDS() );
5405 if ( S.IsNull() ) return;
5409 switch ( S.ShapeType() )
5414 TopLoc_Location loc;
5415 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( S ), loc, f, l );
5416 if ( curve.IsNull() ) return;
5417 SMDS_EdgePosition* ePos = static_cast<SMDS_EdgePosition*>( srcNode->GetPosition() );
5418 p = curve->Value( ePos->GetUParameter() );
5423 p = BRep_Tool::Pnt( TopoDS::Vertex( S ));
5428 getMeshDS()->MoveNode( srcNode, p.X(), p.Y(), p.Z() );
5429 dumpMove( srcNode );
5433 //================================================================================
5435 * \brief Try to fix triangles with high aspect ratio by swaping diagonals
5437 //================================================================================
5439 void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F,
5440 SMESH_MesherHelper& helper,
5443 set<const SMDS_MeshNode*> * involvedNodes)
5445 SMESH::Controls::AspectRatio qualifier;
5446 SMESH::Controls::TSequenceOfXYZ points(3), points1(3), points2(3);
5447 const double maxAspectRatio = is2D ? 4. : 2;
5448 _NodeCoordHelper xyz( F, helper, is2D );
5450 // find bad triangles
5452 vector< const SMDS_MeshElement* > badTrias;
5453 vector< double > badAspects;
5454 SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( F );
5455 SMDS_ElemIteratorPtr fIt = sm->GetElements();
5456 while ( fIt->more() )
5458 const SMDS_MeshElement * f = fIt->next();
5459 if ( f->NbCornerNodes() != 3 ) continue;
5460 for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = xyz( f->GetNode(iP));
5461 double aspect = qualifier.GetValue( points );
5462 if ( aspect > maxAspectRatio )
5464 badTrias.push_back( f );
5465 badAspects.push_back( aspect );
5470 dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID());
5471 SMDS_ElemIteratorPtr fIt = sm->GetElements();
5472 while ( fIt->more() )
5474 const SMDS_MeshElement * f = fIt->next();
5475 if ( f->NbCornerNodes() == 3 )
5476 dumpChangeNodes( f );
5480 if ( badTrias.empty() )
5483 // find couples of faces to swap diagonal
5485 typedef pair < const SMDS_MeshElement* , const SMDS_MeshElement* > T2Trias;
5486 vector< T2Trias > triaCouples;
5488 TIDSortedElemSet involvedFaces, emptySet;
5489 for ( size_t iTia = 0; iTia < badTrias.size(); ++iTia )
5492 double aspRatio [3];
5495 if ( !involvedFaces.insert( badTrias[iTia] ).second )
5497 for ( int iP = 0; iP < 3; ++iP )
5498 points(iP+1) = xyz( badTrias[iTia]->GetNode(iP));
5500 // find triangles adjacent to badTrias[iTia] with better aspect ratio after diag-swaping
5501 int bestCouple = -1;
5502 for ( int iSide = 0; iSide < 3; ++iSide )
5504 const SMDS_MeshNode* n1 = badTrias[iTia]->GetNode( iSide );
5505 const SMDS_MeshNode* n2 = badTrias[iTia]->GetNode(( iSide+1 ) % 3 );
5506 trias [iSide].first = badTrias[iTia];
5507 trias [iSide].second = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, involvedFaces,
5509 if (( ! trias[iSide].second ) ||
5510 ( trias[iSide].second->NbCornerNodes() != 3 ) ||
5511 ( ! sm->Contains( trias[iSide].second )))
5514 // aspect ratio of an adjacent tria
5515 for ( int iP = 0; iP < 3; ++iP )
5516 points2(iP+1) = xyz( trias[iSide].second->GetNode(iP));
5517 double aspectInit = qualifier.GetValue( points2 );
5519 // arrange nodes as after diag-swaping
5520 if ( helper.WrapIndex( i1+1, 3 ) == i2 )
5521 i3 = helper.WrapIndex( i1-1, 3 );
5523 i3 = helper.WrapIndex( i1+1, 3 );
5525 points1( 1+ iSide ) = points2( 1+ i3 );
5526 points2( 1+ i2 ) = points1( 1+ ( iSide+2 ) % 3 );
5528 // aspect ratio after diag-swaping
5529 aspRatio[ iSide ] = qualifier.GetValue( points1 ) + qualifier.GetValue( points2 );
5530 if ( aspRatio[ iSide ] > aspectInit + badAspects[ iTia ] )
5533 // prevent inversion of a triangle
5534 gp_Vec norm1 = gp_Vec( points1(1), points1(3) ) ^ gp_Vec( points1(1), points1(2) );
5535 gp_Vec norm2 = gp_Vec( points2(1), points2(3) ) ^ gp_Vec( points2(1), points2(2) );
5536 if ( norm1 * norm2 < 0. && norm1.Angle( norm2 ) > 70./180.*M_PI )
5539 if ( bestCouple < 0 || aspRatio[ bestCouple ] > aspRatio[ iSide ] )
5543 if ( bestCouple >= 0 )
5545 triaCouples.push_back( trias[bestCouple] );
5546 involvedFaces.insert ( trias[bestCouple].second );
5550 involvedFaces.erase( badTrias[iTia] );
5553 if ( triaCouples.empty() )
5558 SMESH_MeshEditor editor( helper.GetMesh() );
5559 dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
5560 for ( size_t i = 0; i < triaCouples.size(); ++i )
5562 dumpChangeNodes( triaCouples[i].first );
5563 dumpChangeNodes( triaCouples[i].second );
5564 editor.InverseDiag( triaCouples[i].first, triaCouples[i].second );
5567 if ( involvedNodes )
5568 for ( size_t i = 0; i < triaCouples.size(); ++i )
5570 involvedNodes->insert( triaCouples[i].first->begin_nodes(),
5571 triaCouples[i].first->end_nodes() );
5572 involvedNodes->insert( triaCouples[i].second->begin_nodes(),
5573 triaCouples[i].second->end_nodes() );
5576 // just for debug dump resulting triangles
5577 dumpFunction(SMESH_Comment("swapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
5578 for ( size_t i = 0; i < triaCouples.size(); ++i )
5580 dumpChangeNodes( triaCouples[i].first );
5581 dumpChangeNodes( triaCouples[i].second );
5585 //================================================================================
5587 * \brief Move target node to it's final position on the FACE during shrinking
5589 //================================================================================
5591 bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
5592 const TopoDS_Face& F,
5593 SMESH_MesherHelper& helper )
5596 return false; // already at the target position
5598 SMDS_MeshNode* tgtNode = const_cast< SMDS_MeshNode*& >( _nodes.back() );
5600 if ( _sWOL.ShapeType() == TopAbs_FACE )
5602 gp_XY curUV = helper.GetNodeUV( F, tgtNode );
5603 gp_Pnt2d tgtUV( _pos[0].X(), _pos[0].Y() );
5604 gp_Vec2d uvDir( _normal.X(), _normal.Y() );
5605 const double uvLen = tgtUV.Distance( curUV );
5606 const double kSafe = Max( 0.5, 1. - 0.1 * _simplices.size() );
5608 // Select shrinking step such that not to make faces with wrong orientation.
5609 double stepSize = 1e100;
5610 for ( size_t i = 0; i < _simplices.size(); ++i )
5612 // find intersection of 2 lines: curUV-tgtUV and that connecting simplex nodes
5613 gp_XY uvN1 = helper.GetNodeUV( F, _simplices[i]._nPrev );
5614 gp_XY uvN2 = helper.GetNodeUV( F, _simplices[i]._nNext );
5615 gp_XY dirN = uvN2 - uvN1;
5616 double det = uvDir.Crossed( dirN );
5617 if ( Abs( det ) < std::numeric_limits<double>::min() ) continue;
5618 gp_XY dirN2Cur = curUV - uvN1;
5619 double step = dirN.Crossed( dirN2Cur ) / det;
5621 stepSize = Min( step, stepSize );
5624 if ( uvLen <= stepSize )
5629 else if ( stepSize > 0 )
5631 newUV = curUV + uvDir.XY() * stepSize * kSafe;
5637 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
5638 pos->SetUParameter( newUV.X() );
5639 pos->SetVParameter( newUV.Y() );
5642 gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
5643 tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
5644 dumpMove( tgtNode );
5647 else // _sWOL is TopAbs_EDGE
5649 const TopoDS_Edge& E = TopoDS::Edge( _sWOL );
5650 const SMDS_MeshNode* n2 = _simplices[0]._nPrev;
5651 SMDS_EdgePosition* tgtPos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
5653 const double u2 = helper.GetNodeU( E, n2, tgtNode );
5654 const double uSrc = _pos[0].Coord( U_SRC );
5655 const double lenTgt = _pos[0].Coord( LEN_TGT );
5657 double newU = _pos[0].Coord( U_TGT );
5658 if ( lenTgt < 0.99 * fabs( uSrc-u2 )) // n2 got out of src-tgt range
5664 newU = 0.1 * tgtPos->GetUParameter() + 0.9 * u2;
5666 tgtPos->SetUParameter( newU );
5668 gp_XY newUV = helper.GetNodeUV( F, tgtNode, _nodes[0]);
5669 gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
5670 tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
5671 dumpMove( tgtNode );
5677 //================================================================================
5679 * \brief Perform smooth on the FACE
5680 * \retval bool - true if the node has been moved
5682 //================================================================================
5684 bool _SmoothNode::Smooth(int& badNb,
5685 Handle(Geom_Surface)& surface,
5686 SMESH_MesherHelper& helper,
5687 const double refSign,
5691 const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() );
5693 // get uv of surrounding nodes
5694 vector<gp_XY> uv( _simplices.size() );
5695 for ( size_t i = 0; i < _simplices.size(); ++i )
5696 uv[i] = helper.GetNodeUV( face, _simplices[i]._nPrev, _node );
5698 // compute new UV for the node
5700 if ( how == TFI && _simplices.size() == 4 )
5703 for ( size_t i = 0; i < _simplices.size(); ++i )
5704 if ( _simplices[i]._nOpp )
5705 corners[i] = helper.GetNodeUV( face, _simplices[i]._nOpp, _node );
5707 throw SALOME_Exception(LOCALIZED("TFI smoothing: _Simplex::_nOpp not set!"));
5709 newPos = helper.calcTFI ( 0.5, 0.5,
5710 corners[0], corners[1], corners[2], corners[3],
5711 uv[1], uv[2], uv[3], uv[0] );
5713 else if ( how == ANGULAR )
5715 newPos = computeAngularPos( uv, helper.GetNodeUV( face, _node ), refSign );
5717 else if ( how == CENTROIDAL && _simplices.size() > 3 )
5719 // average centers of diagonals wieghted with their reciprocal lengths
5720 if ( _simplices.size() == 4 )
5722 double w1 = 1. / ( uv[2]-uv[0] ).SquareModulus();
5723 double w2 = 1. / ( uv[3]-uv[1] ).SquareModulus();
5724 newPos = ( w1 * ( uv[2]+uv[0] ) + w2 * ( uv[3]+uv[1] )) / ( w1+w2 ) / 2;
5728 double sumWeight = 0;
5729 int nb = _simplices.size() == 4 ? 2 : _simplices.size();
5730 for ( int i = 0; i < nb; ++i )
5733 int iTo = i + _simplices.size() - 1;
5734 for ( int j = iFrom; j < iTo; ++j )
5736 int i2 = SMESH_MesherHelper::WrapIndex( j, _simplices.size() );
5737 double w = 1. / ( uv[i]-uv[i2] ).SquareModulus();
5739 newPos += w * ( uv[i]+uv[i2] );
5742 newPos /= 2 * sumWeight; // 2 is to get a middle between uv's
5748 for ( size_t i = 0; i < _simplices.size(); ++i )
5750 newPos /= _simplices.size();
5753 // count quality metrics (orientation) of triangles around the node
5755 gp_XY tgtUV = helper.GetNodeUV( face, _node );
5756 for ( size_t i = 0; i < _simplices.size(); ++i )
5757 nbOkBefore += _simplices[i].IsForward( tgtUV, _node, face, helper, refSign );
5760 for ( size_t i = 0; i < _simplices.size(); ++i )
5761 nbOkAfter += _simplices[i].IsForward( newPos, _node, face, helper, refSign );
5763 if ( nbOkAfter < nbOkBefore )
5765 badNb += _simplices.size() - nbOkBefore;
5769 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( _node->GetPosition() );
5770 pos->SetUParameter( newPos.X() );
5771 pos->SetVParameter( newPos.Y() );
5778 gp_Pnt p = surface->Value( newPos.X(), newPos.Y() );
5779 const_cast< SMDS_MeshNode* >( _node )->setXYZ( p.X(), p.Y(), p.Z() );
5783 badNb += _simplices.size() - nbOkAfter;
5784 return ( (tgtUV-newPos).SquareModulus() > 1e-10 );
5787 //================================================================================
5789 * \brief Computes new UV using angle based smoothing technic
5791 //================================================================================
5793 gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
5794 const gp_XY& uvToFix,
5795 const double refSign)
5797 uv.push_back( uv.front() );
5799 vector< gp_XY > edgeDir ( uv.size() );
5800 vector< double > edgeSize( uv.size() );
5801 for ( size_t i = 1; i < edgeDir.size(); ++i )
5803 edgeDir [i-1] = uv[i] - uv[i-1];
5804 edgeSize[i-1] = edgeDir[i-1].Modulus();
5805 if ( edgeSize[i-1] < numeric_limits<double>::min() )
5806 edgeDir[i-1].SetX( 100 );
5808 edgeDir[i-1] /= edgeSize[i-1] * refSign;
5810 edgeDir.back() = edgeDir.front();
5811 edgeSize.back() = edgeSize.front();
5816 for ( size_t i = 1; i < edgeDir.size(); ++i )
5818 if ( edgeDir[i-1].X() > 1. ) continue;
5820 while ( edgeDir[i].X() > 1. && ++i < edgeDir.size() );
5821 if ( i == edgeDir.size() ) break;
5823 gp_XY norm1( -edgeDir[i1].Y(), edgeDir[i1].X() );
5824 gp_XY norm2( -edgeDir[i].Y(), edgeDir[i].X() );
5825 gp_XY bisec = norm1 + norm2;
5826 double bisecSize = bisec.Modulus();
5827 if ( bisecSize < numeric_limits<double>::min() )
5829 bisec = -edgeDir[i1] + edgeDir[i];
5830 bisecSize = bisec.Modulus();
5834 gp_XY dirToN = uvToFix - p;
5835 double distToN = dirToN.Modulus();
5836 if ( bisec * dirToN < 0 )
5839 newPos += ( p + bisec * distToN ) * ( edgeSize[i1] + edgeSize[i] );
5841 sumSize += edgeSize[i1] + edgeSize[i];
5843 newPos /= /*nbEdges * */sumSize;
5847 //================================================================================
5849 * \brief Delete _SolidData
5851 //================================================================================
5853 _SolidData::~_SolidData()
5855 for ( size_t i = 0; i < _edges.size(); ++i )
5857 if ( _edges[i] && _edges[i]->_2neibors )
5858 delete _edges[i]->_2neibors;
5863 //================================================================================
5865 * \brief Add a _LayerEdge inflated along the EDGE
5867 //================================================================================
5869 void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper )
5872 if ( _nodes.empty() )
5874 _edges[0] = _edges[1] = 0;
5878 if ( e == _edges[0] || e == _edges[1] )
5880 if ( e->_sWOL.IsNull() || e->_sWOL.ShapeType() != TopAbs_EDGE )
5881 throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
5882 if ( _edges[0] && _edges[0]->_sWOL != e->_sWOL )
5883 throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
5886 const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
5888 BRep_Tool::Range( E, f,l );
5889 double u = helper.GetNodeU( E, e->_nodes[0], e->_nodes.back());
5890 _edges[ u < 0.5*(f+l) ? 0 : 1 ] = e;
5894 const SMDS_MeshNode* tgtNode0 = _edges[0] ? _edges[0]->_nodes.back() : 0;
5895 const SMDS_MeshNode* tgtNode1 = _edges[1] ? _edges[1]->_nodes.back() : 0;
5897 if ( _nodes.empty() )
5899 SMESHDS_SubMesh * eSubMesh = helper.GetMeshDS()->MeshElements( E );
5900 if ( !eSubMesh || eSubMesh->NbNodes() < 1 )
5902 TopLoc_Location loc;
5903 Handle(Geom_Curve) C = BRep_Tool::Curve(E, loc, f,l);
5904 GeomAdaptor_Curve aCurve(C, f,l);
5905 const double totLen = GCPnts_AbscissaPoint::Length(aCurve, f, l);
5907 int nbExpectNodes = eSubMesh->NbNodes();
5908 _initU .reserve( nbExpectNodes );
5909 _normPar.reserve( nbExpectNodes );
5910 _nodes .reserve( nbExpectNodes );
5911 SMDS_NodeIteratorPtr nIt = eSubMesh->GetNodes();
5912 while ( nIt->more() )
5914 const SMDS_MeshNode* node = nIt->next();
5915 if ( node->NbInverseElements(SMDSAbs_Edge) == 0 ||
5916 node == tgtNode0 || node == tgtNode1 )
5917 continue; // refinement nodes
5918 _nodes.push_back( node );
5919 _initU.push_back( helper.GetNodeU( E, node ));
5920 double len = GCPnts_AbscissaPoint::Length(aCurve, f, _initU.back());
5921 _normPar.push_back( len / totLen );
5926 // remove target node of the _LayerEdge from _nodes
5928 for ( size_t i = 0; i < _nodes.size(); ++i )
5929 if ( !_nodes[i] || _nodes[i] == tgtNode0 || _nodes[i] == tgtNode1 )
5930 _nodes[i] = 0, nbFound++;
5931 if ( nbFound == _nodes.size() )
5936 //================================================================================
5938 * \brief Move nodes on EDGE from ends where _LayerEdge's are inflated
5940 //================================================================================
5942 void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
5944 if ( _done || _nodes.empty())
5946 const _LayerEdge* e = _edges[0];
5947 if ( !e ) e = _edges[1];
5950 _done = (( !_edges[0] || _edges[0]->_pos.empty() ) &&
5951 ( !_edges[1] || _edges[1]->_pos.empty() ));
5953 const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
5955 if ( set3D || _done )
5957 Handle(Geom_Curve) C = BRep_Tool::Curve(E, f,l);
5958 GeomAdaptor_Curve aCurve(C, f,l);
5961 f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
5963 l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
5964 double totLen = GCPnts_AbscissaPoint::Length( aCurve, f, l );
5966 for ( size_t i = 0; i < _nodes.size(); ++i )
5968 if ( !_nodes[i] ) continue;
5969 double len = totLen * _normPar[i];
5970 GCPnts_AbscissaPoint discret( aCurve, len, f );
5971 if ( !discret.IsDone() )
5972 return throw SALOME_Exception(LOCALIZED("GCPnts_AbscissaPoint failed"));
5973 double u = discret.Parameter();
5974 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
5975 pos->SetUParameter( u );
5976 gp_Pnt p = C->Value( u );
5977 const_cast< SMDS_MeshNode*>( _nodes[i] )->setXYZ( p.X(), p.Y(), p.Z() );
5982 BRep_Tool::Range( E, f,l );
5984 f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
5986 l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
5988 for ( size_t i = 0; i < _nodes.size(); ++i )
5990 if ( !_nodes[i] ) continue;
5991 double u = f * ( 1-_normPar[i] ) + l * _normPar[i];
5992 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
5993 pos->SetUParameter( u );
5998 //================================================================================
6000 * \brief Restore initial parameters of nodes on EDGE
6002 //================================================================================
6004 void _Shrinker1D::RestoreParams()
6007 for ( size_t i = 0; i < _nodes.size(); ++i )
6009 if ( !_nodes[i] ) continue;
6010 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
6011 pos->SetUParameter( _initU[i] );
6016 //================================================================================
6018 * \brief Replace source nodes by target nodes in shrinked mesh edges
6020 //================================================================================
6022 void _Shrinker1D::SwapSrcTgtNodes( SMESHDS_Mesh* mesh )
6024 const SMDS_MeshNode* nodes[3];
6025 for ( int i = 0; i < 2; ++i )
6027 if ( !_edges[i] ) continue;
6029 SMESHDS_SubMesh * eSubMesh = mesh->MeshElements( _edges[i]->_sWOL );
6030 if ( !eSubMesh ) return;
6031 const SMDS_MeshNode* srcNode = _edges[i]->_nodes[0];
6032 const SMDS_MeshNode* tgtNode = _edges[i]->_nodes.back();
6033 SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
6034 while ( eIt->more() )
6036 const SMDS_MeshElement* e = eIt->next();
6037 if ( !eSubMesh->Contains( e ))
6039 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
6040 for ( int iN = 0; iN < e->NbNodes(); ++iN )
6042 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
6043 nodes[iN] = ( n == srcNode ? tgtNode : n );
6045 mesh->ChangeElementNodes( e, nodes, e->NbNodes() );
6050 //================================================================================
6052 * \brief Creates 2D and 1D elements on boundaries of new prisms
6054 //================================================================================
6056 bool _ViscousBuilder::addBoundaryElements()
6058 SMESH_MesherHelper helper( *_mesh );
6060 for ( size_t i = 0; i < _sdVec.size(); ++i )
6062 _SolidData& data = _sdVec[i];
6063 TopTools_IndexedMapOfShape geomEdges;
6064 TopExp::MapShapes( data._solid, TopAbs_EDGE, geomEdges );
6065 for ( int iE = 1; iE <= geomEdges.Extent(); ++iE )
6067 const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE));
6068 if ( data._noShrinkShapes.count( getMeshDS()->ShapeToIndex( E )))
6071 // Get _LayerEdge's based on E
6073 map< double, const SMDS_MeshNode* > u2nodes;
6074 if ( !SMESH_Algo::GetSortedNodesOnEdge( getMeshDS(), E, /*ignoreMedium=*/false, u2nodes))
6077 vector< _LayerEdge* > ledges; ledges.reserve( u2nodes.size() );
6078 TNode2Edge & n2eMap = data._n2eMap;
6079 map< double, const SMDS_MeshNode* >::iterator u2n = u2nodes.begin();
6081 //check if 2D elements are needed on E
6082 TNode2Edge::iterator n2e = n2eMap.find( u2n->second );
6083 if ( n2e == n2eMap.end() ) continue; // no layers on vertex
6084 ledges.push_back( n2e->second );
6086 if (( n2e = n2eMap.find( u2n->second )) == n2eMap.end() )
6087 continue; // no layers on E
6088 ledges.push_back( n2eMap[ u2n->second ]);
6090 const SMDS_MeshNode* tgtN0 = ledges[0]->_nodes.back();
6091 const SMDS_MeshNode* tgtN1 = ledges[1]->_nodes.back();
6092 int nbSharedPyram = 0;
6093 SMDS_ElemIteratorPtr vIt = tgtN0->GetInverseElementIterator(SMDSAbs_Volume);
6094 while ( vIt->more() )
6096 const SMDS_MeshElement* v = vIt->next();
6097 nbSharedPyram += int( v->GetNodeIndex( tgtN1 ) >= 0 );
6099 if ( nbSharedPyram > 1 )
6100 continue; // not free border of the pyramid
6102 if ( getMeshDS()->FindFace( ledges[0]->_nodes[0], ledges[0]->_nodes[1],
6103 ledges[1]->_nodes[0], ledges[1]->_nodes[1]))
6104 continue; // faces already created
6106 for ( ++u2n; u2n != u2nodes.end(); ++u2n )
6107 ledges.push_back( n2eMap[ u2n->second ]);
6109 // Find out orientation and type of face to create
6111 bool reverse = false, isOnFace;
6113 map< TGeomID, TopoDS_Shape >::iterator e2f =
6114 data._shrinkShape2Shape.find( getMeshDS()->ShapeToIndex( E ));
6116 if (( isOnFace = ( e2f != data._shrinkShape2Shape.end() )))
6118 F = e2f->second.Oriented( TopAbs_FORWARD );
6119 reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED );
6120 if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED )
6121 reverse = !reverse, F.Reverse();
6122 if ( helper.IsReversedSubMesh( TopoDS::Face(F) ))
6127 // find FACE with layers sharing E
6128 PShapeIteratorPtr fIt = helper.GetAncestors( E, *_mesh, TopAbs_FACE );
6129 while ( fIt->more() && F.IsNull() )
6131 const TopoDS_Shape* pF = fIt->next();
6132 if ( helper.IsSubShape( *pF, data._solid) &&
6133 !data._ignoreFaceIds.count( e2f->first ))
6137 // Find the sub-mesh to add new faces
6138 SMESHDS_SubMesh* sm = 0;
6140 sm = getMeshDS()->MeshElements( F );
6142 sm = data._proxyMesh->getFaceSubM( TopoDS::Face(F), /*create=*/true );
6144 return error("error in addBoundaryElements()", data._index);
6147 const int dj1 = reverse ? 0 : 1;
6148 const int dj2 = reverse ? 1 : 0;
6149 for ( size_t j = 1; j < ledges.size(); ++j )
6151 vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes;
6152 vector< const SMDS_MeshNode*>& nn2 = ledges[j-dj2]->_nodes;
6153 if ( nn1.size() == nn2.size() )
6156 for ( size_t z = 1; z < nn1.size(); ++z )
6157 sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
6159 for ( size_t z = 1; z < nn1.size(); ++z )
6160 sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
6162 else if ( nn1.size() == 1 )
6165 for ( size_t z = 1; z < nn2.size(); ++z )
6166 sm->AddElement( getMeshDS()->AddFace( nn1[0], nn2[z-1], nn2[z] ));
6168 for ( size_t z = 1; z < nn2.size(); ++z )
6169 sm->AddElement( new SMDS_FaceOfNodes( nn1[0], nn2[z-1], nn2[z] ));
6174 for ( size_t z = 1; z < nn1.size(); ++z )
6175 sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[0], nn1[z] ));
6177 for ( size_t z = 1; z < nn1.size(); ++z )
6178 sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[0], nn2[z] ));
6183 for ( int isFirst = 0; isFirst < 2; ++isFirst )
6185 _LayerEdge* edge = isFirst ? ledges.front() : ledges.back();
6186 if ( !edge->_sWOL.IsNull() && edge->_sWOL.ShapeType() == TopAbs_EDGE )
6188 vector< const SMDS_MeshNode*>& nn = edge->_nodes;
6189 if ( nn.size() < 2 || nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() )
6191 helper.SetSubShape( edge->_sWOL );
6192 helper.SetElementsOnShape( true );
6193 for ( size_t z = 1; z < nn.size(); ++z )
6194 helper.AddEdge( nn[z-1], nn[z] );