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_Vertex& fromV,
858 const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok,
860 //--------------------------------------------------------------------------------
861 gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE,
862 const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok)
865 Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l );
868 TopoDS_Vertex v = helper.IthVertex( 0, fromE );
869 return getFaceDir( F, v, node, helper, ok );
871 gp_XY uv = helper.GetNodeUV( F, node, 0, &ok );
872 Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
873 gp_Pnt p; gp_Vec du, dv, norm;
874 surface->D1( uv.X(),uv.Y(), p, du,dv );
877 double u = helper.GetNodeU( fromE, node, 0, &ok );
879 TopAbs_Orientation o = helper.GetSubShapeOri( F.Oriented(TopAbs_FORWARD), fromE);
880 if ( o == TopAbs_REVERSED )
883 gp_Vec dir = norm ^ du;
885 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX &&
886 helper.IsClosedEdge( fromE ))
888 if ( fabs(u-f) < fabs(u-l)) c->D1( l, p, dv );
889 else c->D1( f, p, dv );
890 if ( o == TopAbs_REVERSED )
892 gp_Vec dir2 = norm ^ dv;
893 dir = dir.Normalized() + dir2.Normalized();
897 //--------------------------------------------------------------------------------
898 gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
899 const SMDS_MeshNode* node, SMESH_MesherHelper& helper,
900 bool& ok, double* cosin)
902 TopoDS_Face faceFrw = F;
903 faceFrw.Orientation( TopAbs_FORWARD );
904 double f,l; TopLoc_Location loc;
905 TopoDS_Edge edges[2]; // sharing a vertex
909 TopExp_Explorer exp( faceFrw, TopAbs_EDGE );
910 for ( ; exp.More() && nbEdges < 2; exp.Next() )
912 const TopoDS_Edge& e = TopoDS::Edge( exp.Current() );
913 if ( SMESH_Algo::isDegenerated( e )) continue;
914 TopExp::Vertices( e, VV[0], VV[1], /*CumOri=*/true );
915 if ( VV[1].IsSame( fromV )) {
916 nbEdges += edges[ 0 ].IsNull();
919 else if ( VV[0].IsSame( fromV )) {
920 nbEdges += edges[ 1 ].IsNull();
925 gp_XYZ dir(0,0,0), edgeDir[2];
928 // get dirs of edges going fromV
930 for ( size_t i = 0; i < nbEdges && ok; ++i )
932 edgeDir[i] = getEdgeDir( edges[i], fromV );
933 double size2 = edgeDir[i].SquareModulus();
934 if (( ok = size2 > numeric_limits<double>::min() ))
935 edgeDir[i] /= sqrt( size2 );
937 if ( !ok ) return dir;
939 // get angle between the 2 edges
941 double angle = helper.GetAngle( edges[0], edges[1], faceFrw, fromV, &faceNormal );
942 if ( Abs( angle ) < 5 * M_PI/180 )
944 dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] );
948 dir = edgeDir[0] + edgeDir[1];
953 double angle = gp_Vec( edgeDir[0] ).Angle( dir );
954 *cosin = Cos( angle );
957 else if ( nbEdges == 1 )
959 dir = getFaceDir( faceFrw, edges[ edges[0].IsNull() ], node, helper, ok );
960 if ( cosin ) *cosin = 1.;
969 //================================================================================
971 * \brief Returns true if a FACE is bound by a concave EDGE
973 //================================================================================
975 bool isConcave( const TopoDS_Face& F, SMESH_MesherHelper& helper )
977 // if ( helper.Count( F, TopAbs_WIRE, /*useMap=*/false) > 1 )
981 TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE );
982 for ( ; eExp.More(); eExp.Next() )
984 const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
985 if ( SMESH_Algo::isDegenerated( E )) continue;
986 // check if 2D curve is concave
987 BRepAdaptor_Curve2d curve( E, F );
988 const int nbIntervals = curve.NbIntervals( GeomAbs_C2 );
989 TColStd_Array1OfReal intervals(1, nbIntervals + 1 );
990 curve.Intervals( intervals, GeomAbs_C2 );
991 bool isConvex = true;
992 for ( int i = 1; i <= nbIntervals && isConvex; ++i )
994 double u1 = intervals( i );
995 double u2 = intervals( i+1 );
996 curve.D2( 0.5*( u1+u2 ), p, drv1, drv2 );
997 double cross = drv2 ^ drv1;
998 if ( E.Orientation() == TopAbs_REVERSED )
1000 isConvex = ( cross > 0.1 ); //-1e-9 );
1004 //cout << "Concave FACE " << helper.GetMeshDS()->ShapeToIndex( F ) << endl;
1008 // check angles at VERTEXes
1010 TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *helper.GetMesh(), 0, error );
1011 for ( size_t iW = 0; iW < wires.size(); ++iW )
1013 const int nbEdges = wires[iW]->NbEdges();
1014 if ( nbEdges < 2 && SMESH_Algo::isDegenerated( wires[iW]->Edge(0)))
1016 for ( int iE1 = 0; iE1 < nbEdges; ++iE1 )
1018 if ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE1 ))) continue;
1019 int iE2 = ( iE1 + 1 ) % nbEdges;
1020 while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 )))
1021 iE2 = ( iE2 + 1 ) % nbEdges;
1022 double angle = helper.GetAngle( wires[iW]->Edge( iE1 ),
1023 wires[iW]->Edge( iE2 ), F,
1024 wires[iW]->FirstVertex( iE2 ));
1025 if ( angle < -5. * M_PI / 180. )
1031 //--------------------------------------------------------------------------------
1032 // DEBUG. Dump intermediate node positions into a python script
1033 // HOWTO use: run python commands written in a console to see
1034 // construction steps of viscous layers
1040 const char* fname = "/tmp/viscous.py";
1041 cout << "execfile('"<<fname<<"')"<<endl;
1042 py = new ofstream(fname);
1043 *py << "import SMESH" << endl
1044 << "from salome.smesh import smeshBuilder" << endl
1045 << "smesh = smeshBuilder.New(salome.myStudy)" << endl
1046 << "meshSO = smesh.GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
1047 << "mesh = smesh.Mesh( meshSO.GetObject() )"<<endl;
1052 *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Viscous Prisms',"
1053 "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA))"<<endl;
1054 *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Neg Volumes',"
1055 "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_Volume3D,'<',0))"<<endl;
1059 ~PyDump() { Finish(); cout << "NB FUNCTIONS: " << theNbFunc << endl; }
1061 #define dumpFunction(f) { _dumpFunction(f, __LINE__);}
1062 #define dumpMove(n) { _dumpMove(n, __LINE__);}
1063 #define dumpCmd(txt) { _dumpCmd(txt, __LINE__);}
1064 void _dumpFunction(const string& fun, int ln)
1065 { if (py) *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl; ++theNbFunc; }
1066 void _dumpMove(const SMDS_MeshNode* n, int ln)
1067 { if (py) *py<< " mesh.MoveNode( "<<n->GetID()<< ", "<< n->X()
1068 << ", "<<n->Y()<<", "<< n->Z()<< ")\t\t # "<< ln <<endl; }
1069 void _dumpCmd(const string& txt, int ln)
1070 { if (py) *py<< " "<<txt<<" # "<< ln <<endl; }
1071 void dumpFunctionEnd()
1072 { if (py) *py<< " return"<< endl; }
1073 void dumpChangeNodes( const SMDS_MeshElement* f )
1074 { if (py) { *py<< " mesh.ChangeElemNodes( " << f->GetID()<<", [";
1075 for ( int i=1; i < f->NbNodes(); ++i ) *py << f->GetNode(i-1)->GetID()<<", ";
1076 *py << f->GetNode( f->NbNodes()-1 )->GetID() << " ])"<< endl; }}
1077 #define debugMsg( txt ) { cout << txt << " (line: " << __LINE__ << ")" << endl; }
1079 struct PyDump { void Finish() {} };
1080 #define dumpFunction(f) f
1082 #define dumpCmd(txt)
1083 #define dumpFunctionEnd()
1084 #define dumpChangeNodes(f)
1085 #define debugMsg( txt ) {}
1089 using namespace VISCOUS_3D;
1091 //================================================================================
1093 * \brief Constructor of _ViscousBuilder
1095 //================================================================================
1097 _ViscousBuilder::_ViscousBuilder()
1099 _error = SMESH_ComputeError::New(COMPERR_OK);
1103 //================================================================================
1105 * \brief Stores error description and returns false
1107 //================================================================================
1109 bool _ViscousBuilder::error(const string& text, int solidId )
1111 const string prefix = string("Viscous layers builder: ");
1112 _error->myName = COMPERR_ALGO_FAILED;
1113 _error->myComment = prefix + text;
1116 SMESH_subMesh* sm = _mesh->GetSubMeshContaining( solidId );
1117 if ( !sm && !_sdVec.empty() )
1118 sm = _mesh->GetSubMeshContaining( solidId = _sdVec[0]._index );
1119 if ( sm && sm->GetSubShape().ShapeType() == TopAbs_SOLID )
1121 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1122 if ( smError && smError->myAlgo )
1123 _error->myAlgo = smError->myAlgo;
1125 sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1127 // set KO to all solids
1128 for ( size_t i = 0; i < _sdVec.size(); ++i )
1130 if ( _sdVec[i]._index == solidId )
1132 sm = _mesh->GetSubMesh( _sdVec[i]._solid );
1133 if ( !sm->IsEmpty() )
1135 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1136 if ( !smError || smError->IsOK() )
1138 smError = SMESH_ComputeError::New( COMPERR_ALGO_FAILED, prefix + "failed");
1139 sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1143 makeGroupOfLE(); // debug
1148 //================================================================================
1150 * \brief At study restoration, restore event listeners used to clear an inferior
1151 * dim sub-mesh modified by viscous layers
1153 //================================================================================
1155 void _ViscousBuilder::RestoreListeners()
1160 //================================================================================
1162 * \brief computes SMESH_ProxyMesh::SubMesh::_n2n
1164 //================================================================================
1166 bool _ViscousBuilder::MakeN2NMap( _MeshOfSolid* pm )
1168 SMESH_subMesh* solidSM = pm->mySubMeshes.front();
1169 TopExp_Explorer fExp( solidSM->GetSubShape(), TopAbs_FACE );
1170 for ( ; fExp.More(); fExp.Next() )
1172 SMESHDS_SubMesh* srcSmDS = pm->GetMeshDS()->MeshElements( fExp.Current() );
1173 const SMESH_ProxyMesh::SubMesh* prxSmDS = pm->GetProxySubMesh( fExp.Current() );
1175 if ( !srcSmDS || !prxSmDS || !srcSmDS->NbElements() || !prxSmDS->NbElements() )
1177 if ( srcSmDS->GetElements()->next() == prxSmDS->GetElements()->next())
1180 if ( srcSmDS->NbElements() != prxSmDS->NbElements() )
1181 return error( "Different nb elements in a source and a proxy sub-mesh", solidSM->GetId());
1183 SMDS_ElemIteratorPtr srcIt = srcSmDS->GetElements();
1184 SMDS_ElemIteratorPtr prxIt = prxSmDS->GetElements();
1185 while( prxIt->more() )
1187 const SMDS_MeshElement* fSrc = srcIt->next();
1188 const SMDS_MeshElement* fPrx = prxIt->next();
1189 if ( fSrc->NbNodes() != fPrx->NbNodes())
1190 return error( "Different elements in a source and a proxy sub-mesh", solidSM->GetId());
1191 for ( int i = 0 ; i < fPrx->NbNodes(); ++i )
1192 pm->setNode2Node( fSrc->GetNode(i), fPrx->GetNode(i), prxSmDS );
1195 pm->_n2nMapComputed = true;
1199 //================================================================================
1201 * \brief Does its job
1203 //================================================================================
1205 SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh,
1206 const TopoDS_Shape& theShape)
1208 // TODO: set priority of solids during Gen::Compute()
1212 // check if proxy mesh already computed
1213 TopExp_Explorer exp( theShape, TopAbs_SOLID );
1215 return error("No SOLID's in theShape"), _error;
1217 if ( _ViscousListener::GetSolidMesh( _mesh, exp.Current(), /*toCreate=*/false))
1218 return SMESH_ComputeErrorPtr(); // everything already computed
1222 // TODO: ignore already computed SOLIDs
1223 if ( !findSolidsWithLayers())
1226 if ( !findFacesWithLayers() )
1229 for ( size_t i = 0; i < _sdVec.size(); ++i )
1231 if ( ! makeLayer(_sdVec[i]) )
1234 if ( _sdVec[i]._edges.size() == 0 )
1237 if ( ! inflate(_sdVec[i]) )
1240 if ( ! refine(_sdVec[i]) )
1246 addBoundaryElements();
1248 makeGroupOfLE(); // debug
1254 //================================================================================
1256 * \brief Finds SOLIDs to compute using viscous layers. Fills _sdVec
1258 //================================================================================
1260 bool _ViscousBuilder::findSolidsWithLayers()
1263 TopTools_IndexedMapOfShape allSolids;
1264 TopExp::MapShapes( _mesh->GetShapeToMesh(), TopAbs_SOLID, allSolids );
1265 _sdVec.reserve( allSolids.Extent());
1267 SMESH_Gen* gen = _mesh->GetGen();
1268 SMESH_HypoFilter filter;
1269 for ( int i = 1; i <= allSolids.Extent(); ++i )
1271 // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
1272 SMESH_Algo* algo = gen->GetAlgo( *_mesh, allSolids(i) );
1273 if ( !algo ) continue;
1274 // TODO: check if algo is hidden
1275 const list <const SMESHDS_Hypothesis *> & allHyps =
1276 algo->GetUsedHypothesis(*_mesh, allSolids(i), /*ignoreAuxiliary=*/false);
1277 list< const SMESHDS_Hypothesis *>::const_iterator hyp = allHyps.begin();
1278 const StdMeshers_ViscousLayers* viscHyp = 0;
1279 for ( ; hyp != allHyps.end() && !viscHyp; ++hyp )
1280 viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp );
1283 TopoDS_Shape hypShape;
1284 filter.Init( filter.Is( viscHyp ));
1285 _mesh->GetHypothesis( allSolids(i), filter, true, &hypShape );
1287 _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
1290 _sdVec.push_back( _SolidData( allSolids(i), viscHyp, hypShape, proxyMesh ));
1291 _sdVec.back()._index = getMeshDS()->ShapeToIndex( allSolids(i));
1294 if ( _sdVec.empty() )
1296 ( SMESH_Comment(StdMeshers_ViscousLayers::GetHypType()) << " hypothesis not found",0);
1301 //================================================================================
1305 //================================================================================
1307 bool _ViscousBuilder::findFacesWithLayers()
1309 SMESH_MesherHelper helper( *_mesh );
1310 TopExp_Explorer exp;
1311 TopTools_IndexedMapOfShape solids;
1313 // collect all faces to ignore defined by hyp
1314 for ( size_t i = 0; i < _sdVec.size(); ++i )
1316 solids.Add( _sdVec[i]._solid );
1318 vector<TGeomID> ids = _sdVec[i]._hyp->GetBndShapes();
1319 if ( _sdVec[i]._hyp->IsToIgnoreShapes() ) // FACEs to ignore are given
1321 for ( size_t ii = 0; ii < ids.size(); ++ii )
1323 const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[ii] );
1324 if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
1325 _sdVec[i]._ignoreFaceIds.insert( ids[ii] );
1328 else // FACEs with layers are given
1330 exp.Init( _sdVec[i]._solid, TopAbs_FACE );
1331 for ( ; exp.More(); exp.Next() )
1333 TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
1334 if ( find( ids.begin(), ids.end(), faceInd ) == ids.end() )
1335 _sdVec[i]._ignoreFaceIds.insert( faceInd );
1339 // ignore internal FACEs if inlets and outlets are specified
1341 TopTools_IndexedDataMapOfShapeListOfShape solidsOfFace;
1342 if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
1343 TopExp::MapShapesAndAncestors( _sdVec[i]._hypShape,
1344 TopAbs_FACE, TopAbs_SOLID, solidsOfFace);
1346 exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
1347 for ( ; exp.More(); exp.Next() )
1349 const TopoDS_Face& face = TopoDS::Face( exp.Current() );
1350 if ( helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) < 2 )
1353 const TGeomID faceInd = getMeshDS()->ShapeToIndex( face );
1354 if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
1356 int nbSolids = solidsOfFace.FindFromKey( face ).Extent();
1358 _sdVec[i]._ignoreFaceIds.insert( faceInd );
1361 if ( helper.IsReversedSubMesh( face ))
1363 _sdVec[i]._reversedFaceIds.insert( faceInd );
1369 // Find faces to shrink mesh on (solution 2 in issue 0020832);
1370 TopTools_IndexedMapOfShape shapes;
1371 for ( size_t i = 0; i < _sdVec.size(); ++i )
1374 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_EDGE, shapes);
1375 for ( int iE = 1; iE <= shapes.Extent(); ++iE )
1377 const TopoDS_Shape& edge = shapes(iE);
1378 // find 2 faces sharing an edge
1380 PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE);
1381 while ( fIt->more())
1383 const TopoDS_Shape* f = fIt->next();
1384 if ( helper.IsSubShape( *f, _sdVec[i]._solid))
1385 FF[ int( !FF[0].IsNull()) ] = *f;
1387 if( FF[1].IsNull() ) continue; // seam edge can be shared by 1 FACE only
1388 // check presence of layers on them
1390 for ( int j = 0; j < 2; ++j )
1391 ignore[j] = _sdVec[i]._ignoreFaceIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
1392 if ( ignore[0] == ignore[1] )
1393 continue; // nothing interesting
1394 TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ];
1395 // check presence of layers on fWOL within an adjacent SOLID
1396 bool collision = false;
1397 PShapeIteratorPtr sIt = helper.GetAncestors( fWOL, *_mesh, TopAbs_SOLID );
1398 while ( const TopoDS_Shape* solid = sIt->next() )
1399 if ( !solid->IsSame( _sdVec[i]._solid ))
1401 int iSolid = solids.FindIndex( *solid );
1402 int iFace = getMeshDS()->ShapeToIndex( fWOL );
1403 if ( iSolid > 0 && !_sdVec[ iSolid-1 ]._ignoreFaceIds.count( iFace ))
1405 //_sdVec[i]._noShrinkShapes.insert( iFace );
1411 if ( !fWOL.IsNull())
1413 TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
1414 _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
1417 // _shrinkShape2Shape will be used to temporary inflate _LayerEdge's based
1418 // on the edge but shrink won't be performed
1419 _sdVec[i]._noShrinkShapes.insert( edgeInd );
1424 // Exclude from _shrinkShape2Shape FACE's that can't be shrinked since
1425 // the algo of the SOLID sharing the FACE does not support it
1426 set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D");
1427 for ( size_t i = 0; i < _sdVec.size(); ++i )
1429 map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
1430 for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
1432 const TopoDS_Shape& fWOL = e2f->second;
1433 const TGeomID edgeID = e2f->first;
1434 bool notShrinkFace = false;
1435 PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
1436 while ( soIt->more() )
1438 const TopoDS_Shape* solid = soIt->next();
1439 if ( _sdVec[i]._solid.IsSame( *solid )) continue;
1440 SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid );
1441 if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue;
1442 notShrinkFace = true;
1444 for ( ; iSolid < _sdVec.size(); ++iSolid )
1446 if ( _sdVec[iSolid]._solid.IsSame( *solid ) ) {
1447 if ( _sdVec[iSolid]._shrinkShape2Shape.count( edgeID ))
1448 notShrinkFace = false;
1452 if ( notShrinkFace )
1454 _sdVec[i]._noShrinkShapes.insert( edgeID );
1456 // add VERTEXes of the edge in _noShrinkShapes
1457 TopoDS_Shape edge = getMeshDS()->IndexToShape( edgeID );
1458 for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() )
1459 _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( vIt.Value() ));
1461 // check if there is a collision with to-shrink-from EDGEs in iSolid
1462 if ( iSolid == _sdVec.size() )
1463 continue; // no VL in the solid
1465 TopExp::MapShapes( fWOL, TopAbs_EDGE, shapes);
1466 for ( int iE = 1; iE <= shapes.Extent(); ++iE )
1468 const TopoDS_Edge& E = TopoDS::Edge( shapes( iE ));
1469 const TGeomID eID = getMeshDS()->ShapeToIndex( E );
1470 if ( eID == edgeID ||
1471 !_sdVec[iSolid]._shrinkShape2Shape.count( eID ) ||
1472 _sdVec[i]._noShrinkShapes.count( eID ))
1474 for ( int is1st = 0; is1st < 2; ++is1st )
1476 TopoDS_Vertex V = helper.IthVertex( is1st, E );
1477 if ( _sdVec[i]._noShrinkShapes.count( getMeshDS()->ShapeToIndex( V ) ))
1479 // _sdVec[i]._noShrinkShapes.insert( eID );
1480 // V = helper.IthVertex( !is1st, E );
1481 // _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( V ));
1482 //iE = 0; // re-start the loop on EDGEs of fWOL
1483 return error("No way to make a conformal mesh with "
1484 "the given set of faces with layers", _sdVec[i]._index);
1490 } // while ( soIt->more() )
1491 } // loop on _sdVec[i]._shrinkShape2Shape
1492 } // loop on _sdVec to fill in _SolidData::_noShrinkShapes
1494 // Find the SHAPE along which to inflate _LayerEdge based on VERTEX
1496 for ( size_t i = 0; i < _sdVec.size(); ++i )
1499 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_VERTEX, shapes);
1500 for ( int iV = 1; iV <= shapes.Extent(); ++iV )
1502 const TopoDS_Shape& vertex = shapes(iV);
1503 // find faces WOL sharing the vertex
1504 vector< TopoDS_Shape > facesWOL;
1505 int totalNbFaces = 0;
1506 PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE);
1507 while ( fIt->more())
1509 const TopoDS_Shape* f = fIt->next();
1510 if ( helper.IsSubShape( *f, _sdVec[i]._solid ) )
1513 const int fID = getMeshDS()->ShapeToIndex( *f );
1514 if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&&
1515 !_sdVec[i]._noShrinkShapes.count( fID )*/)
1516 facesWOL.push_back( *f );
1519 if ( facesWOL.size() == totalNbFaces || facesWOL.empty() )
1520 continue; // no layers at this vertex or no WOL
1521 TGeomID vInd = getMeshDS()->ShapeToIndex( vertex );
1522 switch ( facesWOL.size() )
1526 helper.SetSubShape( facesWOL[0] );
1527 if ( helper.IsRealSeam( vInd )) // inflate along a seam edge?
1529 TopoDS_Shape seamEdge;
1530 PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1531 while ( eIt->more() && seamEdge.IsNull() )
1533 const TopoDS_Shape* e = eIt->next();
1534 if ( helper.IsRealSeam( *e ) )
1537 if ( !seamEdge.IsNull() )
1539 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge ));
1543 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] ));
1548 // find an edge shared by 2 faces
1549 PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1550 while ( eIt->more())
1552 const TopoDS_Shape* e = eIt->next();
1553 if ( helper.IsSubShape( *e, facesWOL[0]) &&
1554 helper.IsSubShape( *e, facesWOL[1]))
1556 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break;
1562 return error("Not yet supported case", _sdVec[i]._index);
1567 // add FACEs of other SOLIDs to _ignoreFaceIds
1568 for ( size_t i = 0; i < _sdVec.size(); ++i )
1571 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_FACE, shapes);
1573 for ( exp.Init( _mesh->GetShapeToMesh(), TopAbs_FACE ); exp.More(); exp.Next() )
1575 if ( !shapes.Contains( exp.Current() ))
1576 _sdVec[i]._ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( exp.Current() ));
1583 //================================================================================
1585 * \brief Create the inner surface of the viscous layer and prepare data for infation
1587 //================================================================================
1589 bool _ViscousBuilder::makeLayer(_SolidData& data)
1591 // get all sub-shapes to make layers on
1592 set<TGeomID> subIds, faceIds;
1593 subIds = data._noShrinkShapes;
1594 TopExp_Explorer exp( data._solid, TopAbs_FACE );
1595 for ( ; exp.More(); exp.Next() )
1597 SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
1598 if ( ! data._ignoreFaceIds.count( fSubM->GetId() ))
1599 faceIds.insert( fSubM->GetId() );
1600 SMESH_subMeshIteratorPtr subIt = fSubM->getDependsOnIterator(/*includeSelf=*/true);
1601 while ( subIt->more() )
1602 subIds.insert( subIt->next()->GetId() );
1605 // make a map to find new nodes on sub-shapes shared with other SOLID
1606 map< TGeomID, TNode2Edge* >::iterator s2ne;
1607 map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
1608 for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
1610 TGeomID shapeInd = s2s->first;
1611 for ( size_t i = 0; i < _sdVec.size(); ++i )
1613 if ( _sdVec[i]._index == data._index ) continue;
1614 map< TGeomID, TopoDS_Shape >::iterator s2s2 = _sdVec[i]._shrinkShape2Shape.find( shapeInd );
1615 if ( s2s2 != _sdVec[i]._shrinkShape2Shape.end() &&
1616 *s2s == *s2s2 && !_sdVec[i]._n2eMap.empty() )
1618 data._s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap ));
1624 // Create temporary faces and _LayerEdge's
1626 dumpFunction(SMESH_Comment("makeLayers_")<<data._index);
1628 data._stepSize = Precision::Infinite();
1629 data._stepSizeNodes[0] = 0;
1631 SMESH_MesherHelper helper( *_mesh );
1632 helper.SetSubShape( data._solid );
1633 helper.SetElementsOnShape( true );
1635 vector< const SMDS_MeshNode*> newNodes; // of a mesh face
1636 TNode2Edge::iterator n2e2;
1638 // collect _LayerEdge's of shapes they are based on
1639 const int nbShapes = getMeshDS()->MaxShapeIndex();
1640 vector< vector<_LayerEdge*> > edgesByGeom( nbShapes+1 );
1642 for ( set<TGeomID>::iterator id = faceIds.begin(); id != faceIds.end(); ++id )
1644 SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( *id );
1645 if ( !smDS ) return error(SMESH_Comment("Not meshed face ") << *id, data._index );
1647 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( *id ));
1648 SMESH_ProxyMesh::SubMesh* proxySub =
1649 data._proxyMesh->getFaceSubM( F, /*create=*/true);
1651 SMDS_ElemIteratorPtr eIt = smDS->GetElements();
1652 while ( eIt->more() )
1654 const SMDS_MeshElement* face = eIt->next();
1655 newNodes.resize( face->NbCornerNodes() );
1656 double faceMaxCosin = -1;
1657 _LayerEdge* maxCosinEdge = 0;
1658 for ( int i = 0 ; i < face->NbCornerNodes(); ++i )
1660 const SMDS_MeshNode* n = face->GetNode(i);
1661 TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first;
1662 if ( !(*n2e).second )
1665 _LayerEdge* edge = new _LayerEdge();
1667 edge->_nodes.push_back( n );
1668 const int shapeID = n->getshapeId();
1669 const bool noShrink = data._noShrinkShapes.count( shapeID );
1670 edgesByGeom[ shapeID ].push_back( edge );
1672 SMESH_TNodeXYZ xyz( n );
1674 // set edge data or find already refined _LayerEdge and get data from it
1675 if (( !noShrink ) &&
1676 ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE ) &&
1677 (( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() ) &&
1678 (( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end() ))
1680 _LayerEdge* foundEdge = (*n2e2).second;
1681 gp_XYZ lastPos = edge->Copy( *foundEdge, helper );
1682 foundEdge->_pos.push_back( lastPos );
1683 // location of the last node is modified and we restore it by foundEdge->_pos.back()
1684 const_cast< SMDS_MeshNode* >
1685 ( edge->_nodes.back() )->setXYZ( xyz.X(), xyz.Y(), xyz.Z() );
1691 edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ));
1693 if ( !setEdgeData( *edge, subIds, helper, data ))
1696 dumpMove(edge->_nodes.back());
1698 if ( edge->_cosin > faceMaxCosin )
1700 faceMaxCosin = edge->_cosin;
1701 maxCosinEdge = edge;
1704 newNodes[ i ] = n2e->second->_nodes.back();
1706 // create a temporary face
1707 const SMDS_MeshElement* newFace =
1708 new _TmpMeshFace( newNodes, --_tmpFaceID, face->getshapeId() );
1709 proxySub->AddElement( newFace );
1711 // compute inflation step size by min size of element on a convex surface
1712 if ( faceMaxCosin > theMinSmoothCosin )
1713 limitStepSize( data, face, maxCosinEdge );
1714 } // loop on 2D elements on a FACE
1715 } // loop on FACEs of a SOLID
1717 data._epsilon = 1e-7;
1718 if ( data._stepSize < 1. )
1719 data._epsilon *= data._stepSize;
1721 // Put _LayerEdge's into the vector data._edges
1722 if ( !sortEdges( data, edgesByGeom ))
1725 // limit data._stepSize depending on surface curvature and fill data._convexFaces
1726 limitStepSizeByCurvature( data ); // !!! it must be before node substitution in _Simplex
1728 // Set target nodes into _Simplex and _LayerEdge's to _2NearEdges
1729 TNode2Edge::iterator n2e;
1730 const SMDS_MeshNode* nn[2];
1731 for ( size_t i = 0; i < data._edges.size(); ++i )
1733 if ( data._edges[i]->IsOnEdge() )
1735 // get neighbor nodes
1736 bool hasData = ( data._edges[i]->_2neibors->_edges[0] );
1737 if ( hasData ) // _LayerEdge is a copy of another one
1739 nn[0] = data._edges[i]->_2neibors->srcNode(0);
1740 nn[1] = data._edges[i]->_2neibors->srcNode(1);
1742 else if ( !findNeiborsOnEdge( data._edges[i], nn[0],nn[1], data ))
1746 // set neighbor _LayerEdge's
1747 for ( int j = 0; j < 2; ++j )
1749 if (( n2e = data._n2eMap.find( nn[j] )) == data._n2eMap.end() )
1750 return error("_LayerEdge not found by src node", data._index);
1751 data._edges[i]->_2neibors->_edges[j] = n2e->second;
1754 data._edges[i]->SetDataByNeighbors( nn[0], nn[1], helper);
1757 for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j )
1759 _Simplex& s = data._edges[i]->_simplices[j];
1760 s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
1761 s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
1769 //================================================================================
1771 * \brief Compute inflation step size by min size of element on a convex surface
1773 //================================================================================
1775 void _ViscousBuilder::limitStepSize( _SolidData& data,
1776 const SMDS_MeshElement* face,
1777 const _LayerEdge* maxCosinEdge )
1780 double minSize = 10 * data._stepSize;
1781 const int nbNodes = face->NbCornerNodes();
1782 for ( int i = 0; i < nbNodes; ++i )
1784 const SMDS_MeshNode* nextN = face->GetNode( SMESH_MesherHelper::WrapIndex( i+1, nbNodes ));
1785 const SMDS_MeshNode* curN = face->GetNode( i );
1786 if ( nextN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ||
1787 curN-> GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
1789 double dist = SMESH_TNodeXYZ( curN ).Distance( nextN );
1790 if ( dist < minSize )
1791 minSize = dist, iN = i;
1794 double newStep = 0.8 * minSize / maxCosinEdge->_lenFactor;
1795 if ( newStep < data._stepSize )
1797 data._stepSize = newStep;
1798 data._stepSizeCoeff = 0.8 / maxCosinEdge->_lenFactor;
1799 data._stepSizeNodes[0] = face->GetNode( iN );
1800 data._stepSizeNodes[1] = face->GetNode( SMESH_MesherHelper::WrapIndex( iN+1, nbNodes ));
1804 //================================================================================
1806 * \brief Compute inflation step size by min size of element on a convex surface
1808 //================================================================================
1810 void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize )
1812 if ( minSize < data._stepSize )
1814 data._stepSize = minSize;
1815 if ( data._stepSizeNodes[0] )
1818 SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
1819 data._stepSizeCoeff = data._stepSize / dist;
1824 //================================================================================
1826 * \brief Limit data._stepSize by evaluating curvature of shapes and fill data._convexFaces
1828 //================================================================================
1830 void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
1832 const int nbTestPnt = 5; // on a FACE sub-shape
1833 const double minCurvature = 0.9 / data._hyp->GetTotalThickness();
1835 BRepLProp_SLProps surfProp( 2, 1e-6 );
1836 SMESH_MesherHelper helper( *_mesh );
1838 data._convexFaces.clear();
1840 TopExp_Explorer face( data._solid, TopAbs_FACE );
1841 for ( ; face.More(); face.Next() )
1843 const TopoDS_Face& F = TopoDS::Face( face.Current() );
1844 SMESH_subMesh * sm = _mesh->GetSubMesh( F );
1845 const TGeomID faceID = sm->GetId();
1846 if ( data._ignoreFaceIds.count( faceID )) continue;
1848 BRepAdaptor_Surface surface( F, false );
1849 surfProp.SetSurface( surface );
1851 bool isTooCurved = false;
1854 _ConvexFace cnvFace;
1855 const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. );
1856 SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
1857 while ( smIt->more() )
1860 const TGeomID subID = sm->GetId();
1861 // find _LayerEdge's of a sub-shape
1863 if ( data.GetShapeEdges( subID, edgesEnd, &iBeg, &iEnd ))
1864 cnvFace._subIdToEdgeEnd.insert( make_pair( subID, edgesEnd ));
1867 // check concavity and curvature and limit data._stepSize
1868 int nbLEdges = iEnd - iBeg;
1869 int iStep = Max( 1, nbLEdges / nbTestPnt );
1870 for ( ; iBeg < iEnd; iBeg += iStep )
1872 gp_XY uv = helper.GetNodeUV( F, data._edges[ iBeg ]->_nodes[0] );
1873 surfProp.SetParameters( uv.X(), uv.Y() );
1874 if ( !surfProp.IsCurvatureDefined() )
1876 if ( surfProp.MaxCurvature() * oriFactor > minCurvature )
1878 limitStepSize( data, 0.9 / surfProp.MaxCurvature() * oriFactor );
1881 if ( surfProp.MinCurvature() * oriFactor > minCurvature )
1883 limitStepSize( data, 0.9 / surfProp.MinCurvature() * oriFactor );
1887 } // loop on sub-shapes of the FACE
1889 if ( !isTooCurved ) continue;
1891 _ConvexFace & convFace =
1892 data._convexFaces.insert( make_pair( faceID, cnvFace )).first->second;
1895 convFace._normalsFixed = false;
1897 // Fill _ConvexFace::_simplexTestEdges. These _LayerEdge's are used to detect
1898 // prism distortion.
1899 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.find( faceID );
1900 if ( id2end != convFace._subIdToEdgeEnd.end() )
1902 // there are _LayerEdge's on the FACE it-self;
1903 // select _LayerEdge's near EDGEs
1904 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
1905 for ( ; iBeg < iEnd; ++iBeg )
1907 _LayerEdge* ledge = data._edges[ iBeg ];
1908 for ( size_t j = 0; j < ledge->_simplices.size(); ++j )
1909 if ( ledge->_simplices[j]._nNext->GetPosition()->GetDim() < 2 )
1911 convFace._simplexTestEdges.push_back( ledge );
1918 // where there are no _LayerEdge's on a _ConvexFace,
1919 // as e.g. on a fillet surface with no internal nodes - issue 22580,
1920 // so that collision of viscous internal faces is not detected by check of
1921 // intersection of _LayerEdge's with the viscous internal faces.
1923 set< const SMDS_MeshNode* > usedNodes;
1925 // look for _LayerEdge's with null _sWOL
1926 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.begin();
1927 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
1929 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
1930 if ( iBeg >= iEnd || !data._edges[ iBeg ]->_sWOL.IsNull() )
1932 for ( ; iBeg < iEnd; ++iBeg )
1934 _LayerEdge* ledge = data._edges[ iBeg ];
1935 const SMDS_MeshNode* srcNode = ledge->_nodes[0];
1936 if ( !usedNodes.insert( srcNode ).second ) continue;
1938 getSimplices( srcNode, ledge->_simplices, data._ignoreFaceIds, &data );
1939 for ( size_t i = 0; i < ledge->_simplices.size(); ++i )
1941 usedNodes.insert( ledge->_simplices[i]._nPrev );
1942 usedNodes.insert( ledge->_simplices[i]._nNext );
1944 convFace._simplexTestEdges.push_back( ledge );
1948 } // loop on FACEs of data._solid
1951 //================================================================================
1953 * \brief Separate shapes (and _LayerEdge's on them) to smooth from the rest ones
1955 //================================================================================
1957 bool _ViscousBuilder::sortEdges( _SolidData& data,
1958 vector< vector<_LayerEdge*> >& edgesByGeom)
1960 // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's
1961 // boundry inclined at a sharp angle to the shape
1963 list< TGeomID > shapesToSmooth;
1965 SMESH_MesherHelper helper( *_mesh );
1968 for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
1970 vector<_LayerEdge*>& eS = edgesByGeom[iS];
1971 if ( eS.empty() ) continue;
1972 const TopoDS_Shape& S = getMeshDS()->IndexToShape( iS );
1973 bool needSmooth = false;
1974 switch ( S.ShapeType() )
1978 bool isShrinkEdge = !eS[0]->_sWOL.IsNull();
1979 for ( TopoDS_Iterator vIt( S ); vIt.More() && !needSmooth; vIt.Next() )
1981 TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
1982 vector<_LayerEdge*>& eV = edgesByGeom[ iV ];
1983 if ( eV.empty() ) continue;
1984 // double cosin = eV[0]->_cosin;
1986 // ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge));
1989 // gp_Vec dir1, dir2;
1990 // if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE )
1991 // dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ));
1993 // dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ),
1994 // eV[0]->_nodes[0], helper, ok);
1995 // dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
1996 // double angle = dir1.Angle( dir2 );
1997 // cosin = cos( angle );
1999 gp_Vec eDir = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
2000 double angle = eDir.Angle( eV[0]->_normal );
2001 double cosin = Cos( angle );
2002 needSmooth = ( cosin > theMinSmoothCosin );
2008 for ( TopExp_Explorer eExp( S, TopAbs_EDGE ); eExp.More() && !needSmooth; eExp.Next() )
2010 TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
2011 vector<_LayerEdge*>& eE = edgesByGeom[ iE ];
2012 if ( eE.empty() ) continue;
2013 if ( eE[0]->_sWOL.IsNull() )
2015 for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
2016 needSmooth = ( eE[i]->_cosin > theMinSmoothCosin );
2020 const TopoDS_Face& F1 = TopoDS::Face( S );
2021 const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
2022 const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
2023 for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
2025 gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
2026 gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
2027 double angle = dir1.Angle( dir2 );
2028 double cosin = cos( angle );
2029 needSmooth = ( cosin > theMinSmoothCosin );
2041 if ( S.ShapeType() == TopAbs_EDGE ) shapesToSmooth.push_front( iS );
2042 else shapesToSmooth.push_back ( iS );
2045 } // loop on edgesByGeom
2047 data._edges.reserve( data._n2eMap.size() );
2048 data._endEdgeOnShape.clear();
2050 // first we put _LayerEdge's on shapes to smooth
2051 data._nbShapesToSmooth = 0;
2052 list< TGeomID >::iterator gIt = shapesToSmooth.begin();
2053 for ( ; gIt != shapesToSmooth.end(); ++gIt )
2055 vector<_LayerEdge*>& eVec = edgesByGeom[ *gIt ];
2056 if ( eVec.empty() ) continue;
2057 data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
2058 data._endEdgeOnShape.push_back( data._edges.size() );
2059 data._nbShapesToSmooth++;
2063 // then the rest _LayerEdge's
2064 for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
2066 vector<_LayerEdge*>& eVec = edgesByGeom[iS];
2067 if ( eVec.empty() ) continue;
2068 data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
2069 data._endEdgeOnShape.push_back( data._edges.size() );
2076 //================================================================================
2078 * \brief Set data of _LayerEdge needed for smoothing
2079 * \param subIds - ids of sub-shapes of a SOLID to take into account faces from
2081 //================================================================================
2083 bool _ViscousBuilder::setEdgeData(_LayerEdge& edge,
2084 const set<TGeomID>& subIds,
2085 SMESH_MesherHelper& helper,
2088 SMESH_MeshEditor editor(_mesh);
2090 const SMDS_MeshNode* node = edge._nodes[0]; // source node
2091 SMDS_TypeOfPosition posType = node->GetPosition()->GetTypeOfPosition();
2095 edge._curvature = 0;
2097 // --------------------------
2098 // Compute _normal and _cosin
2099 // --------------------------
2102 edge._normal.SetCoord(0,0,0);
2104 int totalNbFaces = 0;
2108 const TGeomID shapeInd = node->getshapeId();
2109 map< TGeomID, TopoDS_Shape >::const_iterator s2s = data._shrinkShape2Shape.find( shapeInd );
2110 const bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
2112 if ( onShrinkShape ) // one of faces the node is on has no layers
2114 TopoDS_Shape vertEdge = getMeshDS()->IndexToShape( s2s->first ); // vertex or edge
2115 if ( s2s->second.ShapeType() == TopAbs_EDGE )
2117 // inflate from VERTEX along EDGE
2118 edge._normal = getEdgeDir( TopoDS::Edge( s2s->second ), TopoDS::Vertex( vertEdge ));
2120 else if ( vertEdge.ShapeType() == TopAbs_VERTEX )
2122 // inflate from VERTEX along FACE
2123 edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Vertex( vertEdge ),
2124 node, helper, normOK, &edge._cosin);
2128 // inflate from EDGE along FACE
2129 edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Edge( vertEdge ),
2130 node, helper, normOK);
2133 else // layers are on all faces of SOLID the node is on
2135 // find indices of geom faces the node lies on
2136 set<TGeomID> faceIds;
2137 if ( posType == SMDS_TOP_FACE )
2139 faceIds.insert( node->getshapeId() );
2143 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2144 while ( fIt->more() )
2145 faceIds.insert( editor.FindShape(fIt->next()));
2148 set<TGeomID>::iterator id = faceIds.begin();
2150 std::pair< TGeomID, gp_XYZ > id2Norm[20];
2151 for ( ; id != faceIds.end(); ++id )
2153 const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id );
2154 if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
2156 F = TopoDS::Face( s );
2157 geomNorm = getFaceNormal( node, F, helper, normOK );
2158 if ( !normOK ) continue;
2160 if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
2162 id2Norm[ totalNbFaces ].first = *id;
2163 id2Norm[ totalNbFaces ].second = geomNorm.XYZ();
2165 edge._normal += geomNorm.XYZ();
2167 if ( totalNbFaces == 0 )
2168 return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index);
2170 if ( normOK && edge._normal.Modulus() < 1e-3 && totalNbFaces > 1 )
2172 // opposite normals, re-get normals at shifted positions (IPAL 52426)
2173 edge._normal.SetCoord( 0,0,0 );
2174 for ( int i = 0; i < totalNbFaces; ++i )
2176 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( id2Norm[i].first ));
2177 geomNorm = getFaceNormal( node, F, helper, normOK, /*shiftInside=*/true );
2178 if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
2181 id2Norm[ i ].second = geomNorm.XYZ();
2182 edge._normal += id2Norm[ i ].second;
2186 if ( totalNbFaces < 3 )
2188 //edge._normal /= totalNbFaces;
2192 edge._normal = getWeigthedNormal( node, id2Norm, totalNbFaces );
2198 edge._cosin = 0; break;
2200 case SMDS_TOP_EDGE: {
2201 TopoDS_Edge E = TopoDS::Edge( helper.GetSubShapeByNode( node, getMeshDS()));
2202 gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK );
2203 double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
2204 edge._cosin = cos( angle );
2205 //cout << "Cosin on EDGE " << edge._cosin << " node " << node->GetID() << endl;
2208 case SMDS_TOP_VERTEX: {
2209 TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS()));
2210 gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK );
2211 double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
2212 edge._cosin = cos( angle );
2213 //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl;
2217 return error(SMESH_Comment("Invalid shape position of node ")<<node, data._index);
2219 } // case _sWOL.IsNull()
2221 double normSize = edge._normal.SquareModulus();
2222 if ( normSize < numeric_limits<double>::min() )
2223 return error(SMESH_Comment("Bad normal at node ")<< node->GetID(), data._index );
2225 edge._normal /= sqrt( normSize );
2227 // TODO: if ( !normOK ) then get normal by mesh faces
2229 // Set the rest data
2230 // --------------------
2231 if ( onShrinkShape )
2233 edge._sWOL = (*s2s).second;
2235 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( edge._nodes.back() );
2236 if ( SMESHDS_SubMesh* sm = getMeshDS()->MeshElements( data._solid ))
2237 sm->RemoveNode( tgtNode , /*isNodeDeleted=*/false );
2239 // set initial position which is parameters on _sWOL in this case
2240 if ( edge._sWOL.ShapeType() == TopAbs_EDGE )
2242 double u = helper.GetNodeU( TopoDS::Edge( edge._sWOL ), node, 0, &normOK );
2243 edge._pos.push_back( gp_XYZ( u, 0, 0 ));
2244 if ( edge._nodes.size() > 1 )
2245 getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( edge._sWOL ), u );
2249 gp_XY uv = helper.GetNodeUV( TopoDS::Face( edge._sWOL ), node, 0, &normOK );
2250 edge._pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
2251 if ( edge._nodes.size() > 1 )
2252 getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( edge._sWOL ), uv.X(), uv.Y() );
2257 edge._pos.push_back( SMESH_TNodeXYZ( node ));
2259 if ( posType == SMDS_TOP_FACE )
2261 getSimplices( node, edge._simplices, data._ignoreFaceIds, &data );
2262 double avgNormProj = 0, avgLen = 0;
2263 for ( size_t i = 0; i < edge._simplices.size(); ++i )
2265 gp_XYZ vec = edge._pos.back() - SMESH_TNodeXYZ( edge._simplices[i]._nPrev );
2266 avgNormProj += edge._normal * vec;
2267 avgLen += vec.Modulus();
2269 avgNormProj /= edge._simplices.size();
2270 avgLen /= edge._simplices.size();
2271 edge._curvature = _Curvature::New( avgNormProj, avgLen );
2275 // Set neighbour nodes for a _LayerEdge based on EDGE
2277 if ( posType == SMDS_TOP_EDGE /*||
2278 ( onShrinkShape && posType == SMDS_TOP_VERTEX && fabs( edge._cosin ) < 1e-10 )*/)
2280 edge._2neibors = new _2NearEdges;
2281 // target node instead of source ones will be set later
2282 // if ( ! findNeiborsOnEdge( &edge,
2283 // edge._2neibors->_nodes[0],
2284 // edge._2neibors->_nodes[1],
2287 // edge.SetDataByNeighbors( edge._2neibors->_nodes[0],
2288 // edge._2neibors->_nodes[1],
2292 edge.SetCosin( edge._cosin ); // to update edge._lenFactor
2297 //================================================================================
2299 * \brief Return normal to a FACE at a node
2300 * \param [in] n - node
2301 * \param [in] face - FACE
2302 * \param [in] helper - helper
2303 * \param [out] isOK - true or false
2304 * \param [in] shiftInside - to find normal at a position shifted inside the face
2305 * \return gp_XYZ - normal
2307 //================================================================================
2309 gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
2310 const TopoDS_Face& face,
2311 SMESH_MesherHelper& helper,
2318 // get a shifted position
2319 gp_Pnt p = SMESH_TNodeXYZ( node );
2320 gp_XYZ shift( 0,0,0 );
2321 TopoDS_Shape S = helper.GetSubShapeByNode( node, helper.GetMeshDS() );
2322 switch ( S.ShapeType() ) {
2325 shift = getFaceDir( face, TopoDS::Vertex( S ), node, helper, isOK );
2330 shift = getFaceDir( face, TopoDS::Edge( S ), node, helper, isOK );
2338 p.Translate( shift * 1e-5 );
2340 TopLoc_Location loc;
2341 GeomAPI_ProjectPointOnSurf& projector = helper.GetProjector( face, loc, 1e-7 );
2343 if ( !loc.IsIdentity() ) p.Transform( loc.Transformation().Inverted() );
2345 projector.Perform( p );
2346 if ( !projector.IsDone() || projector.NbPoints() < 1 )
2351 Quantity_Parameter U,V;
2352 projector.LowerDistanceParameters(U,V);
2357 uv = helper.GetNodeUV( face, node, 0, &isOK );
2363 Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
2364 int pointKind = GeomLib::NormEstim( surface, uv, 1e-5, normal );
2365 enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE };
2367 if ( pointKind == IMPOSSIBLE &&
2368 node->GetPosition()->GetDim() == 2 ) // node inside the FACE
2370 // probably NormEstim() failed due to a too high tolerance
2371 pointKind = GeomLib::NormEstim( surface, uv, 1e-20, normal );
2372 isOK = ( pointKind < IMPOSSIBLE );
2374 if ( pointKind < IMPOSSIBLE )
2376 if ( pointKind != REGULAR &&
2378 node->GetPosition()->GetDim() < 2 ) // FACE boundary
2380 gp_XYZ normShift = getFaceNormal( node, face, helper, isOK, /*shiftInside=*/true );
2381 if ( normShift * normal.XYZ() < 0. )
2387 if ( !isOK ) // hard singularity, to call with shiftInside=true ?
2389 const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face );
2391 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2392 while ( fIt->more() )
2394 const SMDS_MeshElement* f = fIt->next();
2395 if ( f->getshapeId() == faceID )
2397 isOK = SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) normal.XYZ(), /*normalized=*/true );
2400 TopoDS_Face ff = face;
2401 ff.Orientation( TopAbs_FORWARD );
2402 if ( helper.IsReversedSubMesh( ff ))
2409 return normal.XYZ();
2412 //================================================================================
2414 * \brief Return a normal at a node weighted with angles taken by FACEs
2415 * \param [in] n - the node
2416 * \param [in] fId2Normal - FACE ids and normals
2417 * \param [in] nbFaces - nb of FACEs meeting at the node
2418 * \return gp_XYZ - computed normal
2420 //================================================================================
2422 gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode* n,
2423 std::pair< TGeomID, gp_XYZ > fId2Normal[],
2426 gp_XYZ resNorm(0,0,0);
2427 TopoDS_Shape V = SMESH_MesherHelper::GetSubShapeByNode( n, getMeshDS() );
2428 if ( V.ShapeType() != TopAbs_VERTEX )
2430 for ( int i = 0; i < nbFaces; ++i )
2431 resNorm += fId2Normal[i].second / nbFaces ;
2436 for ( int i = 0; i < nbFaces; ++i )
2438 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( fId2Normal[i].first ));
2440 // look for two EDGEs shared by F and other FACEs within fId2Normal
2443 PShapeIteratorPtr eIt = SMESH_MesherHelper::GetAncestors( V, *_mesh, TopAbs_EDGE );
2444 while ( const TopoDS_Shape* E = eIt->next() )
2446 if ( !SMESH_MesherHelper::IsSubShape( *E, F ))
2448 bool isSharedEdge = false;
2449 for ( int j = 0; j < nbFaces && !isSharedEdge; ++j )
2451 if ( i == j ) continue;
2452 const TopoDS_Shape& otherF = getMeshDS()->IndexToShape( fId2Normal[j].first );
2453 isSharedEdge = SMESH_MesherHelper::IsSubShape( *E, otherF );
2455 if ( !isSharedEdge )
2457 ee[ nbE ] = TopoDS::Edge( *E );
2458 ee[ nbE ].Orientation( SMESH_MesherHelper::GetSubShapeOri( F, *E ));
2463 // get an angle between the two EDGEs
2465 if ( nbE < 1 ) continue;
2472 if ( !V.IsSame( SMESH_MesherHelper::IthVertex( 0, ee[ 1 ] )))
2473 std::swap( ee[0], ee[1] );
2475 angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F, TopoDS::Vertex( V ));
2478 // compute a weighted normal
2479 double sumAngle = 0;
2480 for ( int i = 0; i < nbFaces; ++i )
2482 angles[i] = ( angles[i] > 2*M_PI ) ? 0 : M_PI - angles[i];
2483 sumAngle += angles[i];
2485 for ( int i = 0; i < nbFaces; ++i )
2486 resNorm += angles[i] / sumAngle * fId2Normal[i].second;
2491 //================================================================================
2493 * \brief Find 2 neigbor nodes of a node on EDGE
2495 //================================================================================
2497 bool _ViscousBuilder::findNeiborsOnEdge(const _LayerEdge* edge,
2498 const SMDS_MeshNode*& n1,
2499 const SMDS_MeshNode*& n2,
2502 const SMDS_MeshNode* node = edge->_nodes[0];
2503 const int shapeInd = node->getshapeId();
2504 SMESHDS_SubMesh* edgeSM = 0;
2505 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE )
2507 edgeSM = getMeshDS()->MeshElements( shapeInd );
2508 if ( !edgeSM || edgeSM->NbElements() == 0 )
2509 return error(SMESH_Comment("Not meshed EDGE ") << shapeInd, data._index);
2513 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Edge);
2514 while ( eIt->more() && !n2 )
2516 const SMDS_MeshElement* e = eIt->next();
2517 const SMDS_MeshNode* nNeibor = e->GetNode( 0 );
2518 if ( nNeibor == node ) nNeibor = e->GetNode( 1 );
2521 if (!edgeSM->Contains(e)) continue;
2525 TopoDS_Shape s = SMESH_MesherHelper::GetSubShapeByNode(nNeibor, getMeshDS() );
2526 if ( !SMESH_MesherHelper::IsSubShape( s, edge->_sWOL )) continue;
2528 ( iN++ ? n2 : n1 ) = nNeibor;
2531 return error(SMESH_Comment("Wrongly meshed EDGE ") << shapeInd, data._index);
2535 //================================================================================
2537 * \brief Set _curvature and _2neibors->_plnNorm by 2 neigbor nodes residing the same EDGE
2539 //================================================================================
2541 void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
2542 const SMDS_MeshNode* n2,
2543 SMESH_MesherHelper& helper)
2545 if ( _nodes[0]->GetPosition()->GetTypeOfPosition() != SMDS_TOP_EDGE )
2548 gp_XYZ pos = SMESH_TNodeXYZ( _nodes[0] );
2549 gp_XYZ vec1 = pos - SMESH_TNodeXYZ( n1 );
2550 gp_XYZ vec2 = pos - SMESH_TNodeXYZ( n2 );
2554 double sumLen = vec1.Modulus() + vec2.Modulus();
2555 _2neibors->_wgt[0] = 1 - vec1.Modulus() / sumLen;
2556 _2neibors->_wgt[1] = 1 - vec2.Modulus() / sumLen;
2557 double avgNormProj = 0.5 * ( _normal * vec1 + _normal * vec2 );
2558 double avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() );
2559 if ( _curvature ) delete _curvature;
2560 _curvature = _Curvature::New( avgNormProj, avgLen );
2561 // if ( _curvature )
2562 // debugMsg( _nodes[0]->GetID()
2563 // << " CURV r,k: " << _curvature->_r<<","<<_curvature->_k
2564 // << " proj = "<<avgNormProj<< " len = " << avgLen << "| lenDelta(0) = "
2565 // << _curvature->lenDelta(0) );
2569 if ( _sWOL.IsNull() )
2571 TopoDS_Shape S = helper.GetSubShapeByNode( _nodes[0], helper.GetMeshDS() );
2572 gp_XYZ dirE = getEdgeDir( TopoDS::Edge( S ), _nodes[0], helper );
2573 gp_XYZ plnNorm = dirE ^ _normal;
2574 double proj0 = plnNorm * vec1;
2575 double proj1 = plnNorm * vec2;
2576 if ( fabs( proj0 ) > 1e-10 || fabs( proj1 ) > 1e-10 )
2578 if ( _2neibors->_plnNorm ) delete _2neibors->_plnNorm;
2579 _2neibors->_plnNorm = new gp_XYZ( plnNorm.Normalized() );
2584 //================================================================================
2586 * \brief Copy data from a _LayerEdge of other SOLID and based on the same node;
2587 * this and other _LayerEdge's are inflated along a FACE or an EDGE
2589 //================================================================================
2591 gp_XYZ _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
2593 _nodes = other._nodes;
2594 _normal = other._normal;
2596 _lenFactor = other._lenFactor;
2597 _cosin = other._cosin;
2598 _sWOL = other._sWOL;
2599 _2neibors = other._2neibors;
2600 _curvature = 0; std::swap( _curvature, other._curvature );
2601 _2neibors = 0; std::swap( _2neibors, other._2neibors );
2603 gp_XYZ lastPos( 0,0,0 );
2604 if ( _sWOL.ShapeType() == TopAbs_EDGE )
2606 double u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes[0] );
2607 _pos.push_back( gp_XYZ( u, 0, 0));
2609 u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes.back() );
2614 gp_XY uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes[0]);
2615 _pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
2617 uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes.back() );
2618 lastPos.SetX( uv.X() );
2619 lastPos.SetY( uv.Y() );
2624 //================================================================================
2626 * \brief Set _cosin and _lenFactor
2628 //================================================================================
2630 void _LayerEdge::SetCosin( double cosin )
2633 cosin = Abs( _cosin );
2634 _lenFactor = ( /*0.1 < cosin &&*/ cosin < 1-1e-12 ) ? 1./sqrt(1-cosin*cosin) : 1.0;
2637 //================================================================================
2639 * \brief Fills a vector<_Simplex >
2641 //================================================================================
2643 void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node,
2644 vector<_Simplex>& simplices,
2645 const set<TGeomID>& ingnoreShapes,
2646 const _SolidData* dataToCheckOri,
2650 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2651 while ( fIt->more() )
2653 const SMDS_MeshElement* f = fIt->next();
2654 const TGeomID shapeInd = f->getshapeId();
2655 if ( ingnoreShapes.count( shapeInd )) continue;
2656 const int nbNodes = f->NbCornerNodes();
2657 const int srcInd = f->GetNodeIndex( node );
2658 const SMDS_MeshNode* nPrev = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd-1, nbNodes ));
2659 const SMDS_MeshNode* nNext = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+1, nbNodes ));
2660 const SMDS_MeshNode* nOpp = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+2, nbNodes ));
2661 if ( dataToCheckOri && dataToCheckOri->_reversedFaceIds.count( shapeInd ))
2662 std::swap( nPrev, nNext );
2663 simplices.push_back( _Simplex( nPrev, nNext, nOpp ));
2668 vector<_Simplex> sortedSimplices( simplices.size() );
2669 sortedSimplices[0] = simplices[0];
2671 for ( size_t i = 1; i < simplices.size(); ++i )
2673 for ( size_t j = 1; j < simplices.size(); ++j )
2674 if ( sortedSimplices[i-1]._nNext == simplices[j]._nPrev )
2676 sortedSimplices[i] = simplices[j];
2681 if ( nbFound == simplices.size() - 1 )
2682 simplices.swap( sortedSimplices );
2686 //================================================================================
2688 * \brief DEBUG. Create groups contating temorary data of _LayerEdge's
2690 //================================================================================
2692 void _ViscousBuilder::makeGroupOfLE()
2695 for ( size_t i = 0 ; i < _sdVec.size(); ++i )
2697 if ( _sdVec[i]._edges.empty() ) continue;
2699 dumpFunction( SMESH_Comment("make_LayerEdge_") << i );
2700 for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
2702 _LayerEdge* le = _sdVec[i]._edges[j];
2703 for ( size_t iN = 1; iN < le->_nodes.size(); ++iN )
2704 dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<le->_nodes[iN-1]->GetID()
2705 << ", " << le->_nodes[iN]->GetID() <<"])");
2709 dumpFunction( SMESH_Comment("makeNormals") << i );
2710 for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
2712 _LayerEdge& edge = *_sdVec[i]._edges[j];
2713 SMESH_TNodeXYZ nXYZ( edge._nodes[0] );
2714 nXYZ += edge._normal * _sdVec[i]._stepSize;
2715 dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<edge._nodes[0]->GetID()
2716 << ", mesh.AddNode( " << nXYZ.X()<<","<< nXYZ.Y()<<","<< nXYZ.Z()<<")])");
2720 dumpFunction( SMESH_Comment("makeTmpFaces_") << i );
2721 TopExp_Explorer fExp( _sdVec[i]._solid, TopAbs_FACE );
2722 for ( ; fExp.More(); fExp.Next() )
2724 if (const SMESHDS_SubMesh* sm = _sdVec[i]._proxyMesh->GetProxySubMesh( fExp.Current()))
2726 SMDS_ElemIteratorPtr fIt = sm->GetElements();
2727 while ( fIt->more())
2729 const SMDS_MeshElement* e = fIt->next();
2730 SMESH_Comment cmd("mesh.AddFace([");
2731 for ( int j=0; j < e->NbCornerNodes(); ++j )
2732 cmd << e->GetNode(j)->GetID() << (j+1<e->NbCornerNodes() ? ",": "])");
2742 //================================================================================
2744 * \brief Increase length of _LayerEdge's to reach the required thickness of layers
2746 //================================================================================
2748 bool _ViscousBuilder::inflate(_SolidData& data)
2750 SMESH_MesherHelper helper( *_mesh );
2752 // Limit inflation step size by geometry size found by itersecting
2753 // normals of _LayerEdge's with mesh faces
2754 double geomSize = Precision::Infinite(), intersecDist;
2755 auto_ptr<SMESH_ElementSearcher> searcher
2756 ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
2757 data._proxyMesh->GetFaces( data._solid )) );
2758 for ( size_t i = 0; i < data._edges.size(); ++i )
2760 if ( data._edges[i]->IsOnEdge() ) continue;
2761 data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon );
2762 if ( geomSize > intersecDist && intersecDist > 0 )
2763 geomSize = intersecDist;
2765 if ( data._stepSize > 0.3 * geomSize )
2766 limitStepSize( data, 0.3 * geomSize );
2768 const double tgtThick = data._hyp->GetTotalThickness();
2769 if ( data._stepSize > tgtThick )
2770 limitStepSize( data, tgtThick );
2772 if ( data._stepSize < 1. )
2773 data._epsilon = data._stepSize * 1e-7;
2775 debugMsg( "-- geomSize = " << geomSize << ", stepSize = " << data._stepSize );
2777 double avgThick = 0, curThick = 0, distToIntersection = Precision::Infinite();
2778 int nbSteps = 0, nbRepeats = 0;
2779 while ( 1.01 * avgThick < tgtThick )
2781 // new target length
2782 curThick += data._stepSize;
2783 if ( curThick > tgtThick )
2785 curThick = tgtThick + ( tgtThick-avgThick ) * nbRepeats;
2789 // Elongate _LayerEdge's
2790 dumpFunction(SMESH_Comment("inflate")<<data._index<<"_step"<<nbSteps); // debug
2791 for ( size_t i = 0; i < data._edges.size(); ++i )
2793 data._edges[i]->SetNewLength( curThick, helper );
2797 if ( !updateNormals( data, helper, nbSteps ))
2800 // Improve and check quality
2801 if ( !smoothAndCheck( data, nbSteps, distToIntersection ))
2805 dumpFunction(SMESH_Comment("invalidate")<<data._index<<"_step"<<nbSteps); // debug
2806 for ( size_t i = 0; i < data._edges.size(); ++i )
2808 data._edges[i]->InvalidateStep( nbSteps+1 );
2812 break; // no more inflating possible
2816 // Evaluate achieved thickness
2818 for ( size_t i = 0; i < data._edges.size(); ++i )
2819 avgThick += data._edges[i]->_len;
2820 avgThick /= data._edges.size();
2821 debugMsg( "-- Thickness " << avgThick << " reached" );
2823 if ( distToIntersection < avgThick*1.5 )
2825 debugMsg( "-- Stop inflation since "
2826 << " distToIntersection( "<<distToIntersection<<" ) < avgThick( "
2827 << avgThick << " ) * 1.5" );
2831 limitStepSize( data, 0.25 * distToIntersection );
2832 if ( data._stepSizeNodes[0] )
2833 data._stepSize = data._stepSizeCoeff *
2834 SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
2836 } // while ( 1.01 * avgThick < tgtThick )
2839 return error("failed at the very first inflation step", data._index);
2841 if ( 1.01 * avgThick < tgtThick )
2842 if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( data._index ))
2844 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
2845 if ( !smError || smError->IsOK() )
2847 ( new SMESH_ComputeError (COMPERR_WARNING,
2848 SMESH_Comment("Thickness ") << tgtThick <<
2849 " of viscous layers not reached,"
2850 " average reached thickness is " << avgThick ));
2854 // Restore position of src nodes moved by infaltion on _noShrinkShapes
2855 dumpFunction(SMESH_Comment("restoNoShrink_So")<<data._index); // debug
2857 for ( int iS = 0; iS < data._endEdgeOnShape.size(); ++iS )
2860 iEnd = data._endEdgeOnShape[ iS ];
2861 if ( data._edges[ iBeg ]->_nodes.size() == 1 )
2862 for ( ; iBeg < iEnd; ++iBeg )
2864 restoreNoShrink( *data._edges[ iBeg ] );
2872 //================================================================================
2874 * \brief Improve quality of layer inner surface and check intersection
2876 //================================================================================
2878 bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
2880 double & distToIntersection)
2882 if ( data._nbShapesToSmooth == 0 )
2883 return true; // no shapes needing smoothing
2885 bool moved, improved;
2887 SMESH_MesherHelper helper(*_mesh);
2888 Handle(Geom_Surface) surface;
2892 for ( int iS = 0; iS < data._nbShapesToSmooth; ++iS )
2895 iEnd = data._endEdgeOnShape[ iS ];
2897 if ( !data._edges[ iBeg ]->_sWOL.IsNull() &&
2898 data._edges[ iBeg ]->_sWOL.ShapeType() == TopAbs_FACE )
2900 if ( !F.IsSame( data._edges[ iBeg ]->_sWOL )) {
2901 F = TopoDS::Face( data._edges[ iBeg ]->_sWOL );
2902 helper.SetSubShape( F );
2903 surface = BRep_Tool::Surface( F );
2908 F.Nullify(); surface.Nullify();
2910 TGeomID sInd = data._edges[ iBeg ]->_nodes[0]->getshapeId();
2912 if ( data._edges[ iBeg ]->IsOnEdge() )
2914 dumpFunction(SMESH_Comment("smooth")<<data._index << "_Ed"<<sInd <<"_InfStep"<<nbSteps);
2916 // try a simple solution on an analytic EDGE
2917 if ( !smoothAnalyticEdge( data, iBeg, iEnd, surface, F, helper ))
2923 for ( int i = iBeg; i < iEnd; ++i )
2925 moved |= data._edges[i]->SmoothOnEdge(surface, F, helper);
2927 dumpCmd( SMESH_Comment("# end step ")<<step);
2929 while ( moved && step++ < 5 );
2936 int step = 0, stepLimit = 5, badNb = 0; moved = true;
2937 while (( ++step <= stepLimit && moved ) || improved )
2939 dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
2940 <<"_InfStep"<<nbSteps<<"_"<<step); // debug
2941 int oldBadNb = badNb;
2945 for ( int i = iBeg; i < iEnd; ++i )
2946 moved |= data._edges[i]->Smooth(badNb);
2948 for ( int i = iEnd-1; i >= iBeg; --i )
2949 moved |= data._edges[i]->Smooth(badNb);
2950 improved = ( badNb < oldBadNb );
2952 // issue 22576 -- no bad faces but still there are intersections to fix
2953 if ( improved && badNb == 0 )
2954 stepLimit = step + 3;
2961 for ( int i = iBeg; i < iEnd; ++i )
2963 _LayerEdge* edge = data._edges[i];
2964 SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
2965 for ( size_t j = 0; j < edge->_simplices.size(); ++j )
2966 if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
2968 cout << "Bad simplex ( " << edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
2969 << " "<< edge->_simplices[j]._nPrev->GetID()
2970 << " "<< edge->_simplices[j]._nNext->GetID() << " )" << endl;
2978 } // loop on shapes to smooth
2980 // Check orientation of simplices of _ConvexFace::_simplexTestEdges
2981 map< TGeomID, _ConvexFace >::iterator id2face = data._convexFaces.begin();
2982 for ( ; id2face != data._convexFaces.end(); ++id2face )
2984 _ConvexFace & convFace = (*id2face).second;
2985 if ( !convFace._simplexTestEdges.empty() &&
2986 convFace._simplexTestEdges[0]->_nodes[0]->GetPosition()->GetDim() == 2 )
2987 continue; // _simplexTestEdges are based on FACE -- already checked while smoothing
2989 if ( !convFace.CheckPrisms() )
2993 // Check if the last segments of _LayerEdge intersects 2D elements;
2994 // checked elements are either temporary faces or faces on surfaces w/o the layers
2996 auto_ptr<SMESH_ElementSearcher> searcher
2997 ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
2998 data._proxyMesh->GetFaces( data._solid )) );
3000 distToIntersection = Precision::Infinite();
3002 const SMDS_MeshElement* intFace = 0;
3003 const SMDS_MeshElement* closestFace = 0;
3005 for ( size_t i = 0; i < data._edges.size(); ++i )
3007 if ( !data._edges[i]->_sWOL.IsNull() )
3009 if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace ))
3011 if ( distToIntersection > dist )
3013 // ignore intersection of a _LayerEdge based on a _ConvexFace with a face
3014 // lying on this _ConvexFace
3015 if ( _ConvexFace* convFace = data.GetConvexFace( intFace->getshapeId() ))
3016 if ( convFace->_subIdToEdgeEnd.count ( data._edges[i]->_nodes[0]->getshapeId() ))
3019 distToIntersection = dist;
3021 closestFace = intFace;
3027 SMDS_MeshElement::iterator nIt = closestFace->begin_nodes();
3028 cout << "Shortest distance: _LayerEdge nodes: tgt " << data._edges[iLE]->_nodes.back()->GetID()
3029 << " src " << data._edges[iLE]->_nodes[0]->GetID()<< ", intersection with face ("
3030 << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
3031 << ") distance = " << distToIntersection<< endl;
3038 //================================================================================
3040 * \brief Return a curve of the EDGE to be used for smoothing and arrange
3041 * _LayerEdge's to be in a consequent order
3043 //================================================================================
3045 Handle(Geom_Curve) _SolidData::CurveForSmooth( const TopoDS_Edge& E,
3048 Handle(Geom_Surface)& surface,
3049 const TopoDS_Face& F,
3050 SMESH_MesherHelper& helper)
3052 TGeomID eIndex = helper.GetMeshDS()->ShapeToIndex( E );
3054 map< TGeomID, Handle(Geom_Curve)>::iterator i2curve = _edge2curve.find( eIndex );
3056 if ( i2curve == _edge2curve.end() )
3058 // sort _LayerEdge's by position on the EDGE
3059 SortOnEdge( E, iFrom, iTo, helper );
3061 SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( eIndex );
3063 TopLoc_Location loc; double f,l;
3065 Handle(Geom_Line) line;
3066 Handle(Geom_Circle) circle;
3067 bool isLine, isCirc;
3068 if ( F.IsNull() ) // 3D case
3070 // check if the EDGE is a line
3071 Handle(Geom_Curve) curve = BRep_Tool::Curve( E, loc, f, l);
3072 if ( curve->IsKind( STANDARD_TYPE( Geom_TrimmedCurve )))
3073 curve = Handle(Geom_TrimmedCurve)::DownCast( curve )->BasisCurve();
3075 line = Handle(Geom_Line)::DownCast( curve );
3076 circle = Handle(Geom_Circle)::DownCast( curve );
3077 isLine = (!line.IsNull());
3078 isCirc = (!circle.IsNull());
3080 if ( !isLine && !isCirc ) // Check if the EDGE is close to a line
3083 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
3084 while ( nIt->more() )
3085 bndBox.Add( SMESH_TNodeXYZ( nIt->next() ));
3086 gp_XYZ size = bndBox.CornerMax() - bndBox.CornerMin();
3088 SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->tgtNode(0) );
3089 SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->tgtNode(1) );
3090 const double lineTol = 1e-2 * ( p0 - p1 ).Modulus();
3091 for ( int i = 0; i < 3 && !isLine; ++i )
3092 isLine = ( size.Coord( i+1 ) <= lineTol );
3094 if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle
3101 // check if the EDGE is a line
3102 Handle(Geom2d_Curve) curve = BRep_Tool::CurveOnSurface( E, F, f, l);
3103 if ( curve->IsKind( STANDARD_TYPE( Geom2d_TrimmedCurve )))
3104 curve = Handle(Geom2d_TrimmedCurve)::DownCast( curve )->BasisCurve();
3106 Handle(Geom2d_Line) line2d = Handle(Geom2d_Line)::DownCast( curve );
3107 Handle(Geom2d_Circle) circle2d = Handle(Geom2d_Circle)::DownCast( curve );
3108 isLine = (!line2d.IsNull());
3109 isCirc = (!circle2d.IsNull());
3111 if ( !isLine && !isCirc) // Check if the EDGE is close to a line
3114 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
3115 while ( nIt->more() )
3116 bndBox.Add( helper.GetNodeUV( F, nIt->next() ));
3117 gp_XY size = bndBox.CornerMax() - bndBox.CornerMin();
3119 const double lineTol = 1e-2 * sqrt( bndBox.SquareExtent() );
3120 for ( int i = 0; i < 2 && !isLine; ++i )
3121 isLine = ( size.Coord( i+1 ) <= lineTol );
3123 if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle
3129 line = new Geom_Line( gp::OX() ); // only type does matter
3133 gp_Pnt2d p = circle2d->Location();
3134 gp_Ax2 ax( gp_Pnt( p.X(), p.Y(), 0), gp::DX());
3135 circle = new Geom_Circle( ax, 1.); // only center position does matter
3139 Handle(Geom_Curve)& res = _edge2curve[ eIndex ];
3147 return i2curve->second;
3150 //================================================================================
3152 * \brief Sort _LayerEdge's by a parameter on a given EDGE
3154 //================================================================================
3156 void _SolidData::SortOnEdge( const TopoDS_Edge& E,
3159 SMESH_MesherHelper& helper)
3161 map< double, _LayerEdge* > u2edge;
3162 for ( int i = iFrom; i < iTo; ++i )
3163 u2edge.insert( make_pair( helper.GetNodeU( E, _edges[i]->_nodes[0] ), _edges[i] ));
3165 ASSERT( u2edge.size() == iTo - iFrom );
3166 map< double, _LayerEdge* >::iterator u2e = u2edge.begin();
3167 for ( int i = iFrom; i < iTo; ++i, ++u2e )
3168 _edges[i] = u2e->second;
3170 // set _2neibors according to the new order
3171 for ( int i = iFrom; i < iTo-1; ++i )
3172 if ( _edges[i]->_2neibors->tgtNode(1) != _edges[i+1]->_nodes.back() )
3173 _edges[i]->_2neibors->reverse();
3174 if ( u2edge.size() > 1 &&
3175 _edges[iTo-1]->_2neibors->tgtNode(0) != _edges[iTo-2]->_nodes.back() )
3176 _edges[iTo-1]->_2neibors->reverse();
3179 //================================================================================
3181 * \brief Return index corresponding to the shape in _endEdgeOnShape
3183 //================================================================================
3185 bool _SolidData::GetShapeEdges(const TGeomID shapeID,
3190 int beg = 0, end = 0;
3191 for ( edgesEnd = 0; edgesEnd < _endEdgeOnShape.size(); ++edgesEnd )
3193 end = _endEdgeOnShape[ edgesEnd ];
3194 TGeomID sID = _edges[ beg ]->_nodes[0]->getshapeId();
3195 if ( sID == shapeID )
3197 if ( iBeg ) *iBeg = beg;
3198 if ( iEnd ) *iEnd = end;
3206 //================================================================================
3208 * \brief Add faces for smoothing
3210 //================================================================================
3212 void _SolidData::AddShapesToSmooth( const set< TGeomID >& faceIDs )
3214 // convert faceIDs to indices in _endEdgeOnShape
3215 set< size_t > iEnds;
3217 set< TGeomID >::const_iterator fId = faceIDs.begin();
3218 for ( ; fId != faceIDs.end(); ++fId )
3219 if ( GetShapeEdges( *fId, end ) && end >= _nbShapesToSmooth )
3220 iEnds.insert( end );
3222 set< size_t >::iterator endsIt = iEnds.begin();
3224 // "add" by move of _nbShapesToSmooth
3225 int nbFacesToAdd = iEnds.size();
3226 while ( endsIt != iEnds.end() && *endsIt == _nbShapesToSmooth )
3229 ++_nbShapesToSmooth;
3232 if ( endsIt == iEnds.end() )
3235 // Move _LayerEdge's on FACEs just after _nbShapesToSmooth
3237 vector< _LayerEdge* > nonSmoothLE, smoothLE;
3238 size_t lastSmooth = *iEnds.rbegin();
3240 for ( size_t i = _nbShapesToSmooth; i <= lastSmooth; ++i )
3242 vector< _LayerEdge* > & edgesVec = iEnds.count(i) ? smoothLE : nonSmoothLE;
3243 iBeg = i ? _endEdgeOnShape[ i-1 ] : 0;
3244 iEnd = _endEdgeOnShape[ i ];
3245 edgesVec.insert( edgesVec.end(), _edges.begin() + iBeg, _edges.begin() + iEnd );
3248 iBeg = _nbShapesToSmooth ? _endEdgeOnShape[ _nbShapesToSmooth-1 ] : 0;
3249 std::copy( smoothLE.begin(), smoothLE.end(), &_edges[ iBeg ] );
3250 std::copy( nonSmoothLE.begin(), nonSmoothLE.end(), &_edges[ iBeg + smoothLE.size()]);
3252 // update _endEdgeOnShape
3253 for ( size_t i = _nbShapesToSmooth; i < _endEdgeOnShape.size(); ++i )
3255 TGeomID curShape = _edges[ iBeg ]->_nodes[0]->getshapeId();
3256 while ( ++iBeg < _edges.size() &&
3257 curShape == _edges[ iBeg ]->_nodes[0]->getshapeId() );
3259 _endEdgeOnShape[ i ] = iBeg;
3262 _nbShapesToSmooth += nbFacesToAdd;
3265 //================================================================================
3267 * \brief smooth _LayerEdge's on a staight EDGE or circular EDGE
3269 //================================================================================
3271 bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data,
3274 Handle(Geom_Surface)& surface,
3275 const TopoDS_Face& F,
3276 SMESH_MesherHelper& helper)
3278 TopoDS_Shape S = helper.GetSubShapeByNode( data._edges[ iFrom ]->_nodes[0],
3279 helper.GetMeshDS());
3280 TopoDS_Edge E = TopoDS::Edge( S );
3282 Handle(Geom_Curve) curve = data.CurveForSmooth( E, iFrom, iTo, surface, F, helper );
3283 if ( curve.IsNull() ) return false;
3285 // compute a relative length of segments
3286 vector< double > len( iTo-iFrom+1 );
3288 double curLen, prevLen = len[0] = 1.0;
3289 for ( int i = iFrom; i < iTo; ++i )
3291 curLen = prevLen * data._edges[i]->_2neibors->_wgt[0] / data._edges[i]->_2neibors->_wgt[1];
3292 len[i-iFrom+1] = len[i-iFrom] + curLen;
3297 if ( curve->IsKind( STANDARD_TYPE( Geom_Line )))
3299 if ( F.IsNull() ) // 3D
3301 SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->tgtNode(0));
3302 SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->tgtNode(1));
3303 for ( int i = iFrom; i < iTo; ++i )
3305 double r = len[i-iFrom] / len.back();
3306 gp_XYZ newPos = p0 * ( 1. - r ) + p1 * r;
3307 data._edges[i]->_pos.back() = newPos;
3308 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3309 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3310 dumpMove( tgtNode );
3315 // gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->tgtNode(0));
3316 // gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->tgtNode(1));
3317 gp_XY uv0 = data._edges[iFrom]->_2neibors->_edges[0]->LastUV( F );
3318 gp_XY uv1 = data._edges[iTo-1]->_2neibors->_edges[1]->LastUV( F );
3319 if ( data._edges[iFrom]->_2neibors->tgtNode(0) ==
3320 data._edges[iTo-1]->_2neibors->tgtNode(1) ) // closed edge
3322 int iPeriodic = helper.GetPeriodicIndex();
3323 if ( iPeriodic == 1 || iPeriodic == 2 )
3325 uv1.SetCoord( iPeriodic, helper.GetOtherParam( uv1.Coord( iPeriodic )));
3326 if ( uv0.Coord( iPeriodic ) > uv1.Coord( iPeriodic ))
3327 std::swap( uv0, uv1 );
3330 const gp_XY rangeUV = uv1 - uv0;
3331 for ( int i = iFrom; i < iTo; ++i )
3333 double r = len[i-iFrom] / len.back();
3334 gp_XY newUV = uv0 + r * rangeUV;
3335 data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
3337 gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
3338 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3339 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3340 dumpMove( tgtNode );
3342 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
3343 pos->SetUParameter( newUV.X() );
3344 pos->SetVParameter( newUV.Y() );
3350 if ( curve->IsKind( STANDARD_TYPE( Geom_Circle )))
3352 Handle(Geom_Circle) circle = Handle(Geom_Circle)::DownCast( curve );
3353 gp_Pnt center3D = circle->Location();
3355 if ( F.IsNull() ) // 3D
3357 if ( data._edges[iFrom]->_2neibors->tgtNode(0) ==
3358 data._edges[iTo-1]->_2neibors->tgtNode(1) )
3359 return true; // closed EDGE - nothing to do
3361 return false; // TODO ???
3365 const gp_XY center( center3D.X(), center3D.Y() );
3367 gp_XY uv0 = data._edges[iFrom]->_2neibors->_edges[0]->LastUV( F );
3368 gp_XY uvM = data._edges[iFrom]->LastUV( F );
3369 gp_XY uv1 = data._edges[iTo-1]->_2neibors->_edges[1]->LastUV( F );
3370 // gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->tgtNode(0));
3371 // gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back());
3372 // gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->tgtNode(1));
3373 gp_Vec2d vec0( center, uv0 );
3374 gp_Vec2d vecM( center, uvM );
3375 gp_Vec2d vec1( center, uv1 );
3376 double uLast = vec0.Angle( vec1 ); // -PI - +PI
3377 double uMidl = vec0.Angle( vecM );
3378 if ( uLast * uMidl <= 0. )
3379 uLast += ( uMidl > 0 ? +2. : -2. ) * M_PI;
3380 const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() );
3382 gp_Ax2d axis( center, vec0 );
3383 gp_Circ2d circ( axis, radius );
3384 for ( int i = iFrom; i < iTo; ++i )
3386 double newU = uLast * len[i-iFrom] / len.back();
3387 gp_Pnt2d newUV = ElCLib::Value( newU, circ );
3388 data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
3390 gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
3391 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3392 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3393 dumpMove( tgtNode );
3395 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
3396 pos->SetUParameter( newUV.X() );
3397 pos->SetVParameter( newUV.Y() );
3406 //================================================================================
3408 * \brief Modify normals of _LayerEdge's on EDGE's to avoid intersection with
3409 * _LayerEdge's on neighbor EDGE's
3411 //================================================================================
3413 bool _ViscousBuilder::updateNormals( _SolidData& data,
3414 SMESH_MesherHelper& helper,
3418 return updateNormalsOfConvexFaces( data, helper, stepNb );
3420 // make temporary quadrangles got by extrusion of
3421 // mesh edges along _LayerEdge._normal's
3423 vector< const SMDS_MeshElement* > tmpFaces;
3425 set< SMESH_TLink > extrudedLinks; // contains target nodes
3426 vector< const SMDS_MeshNode*> nodes(4); // of a tmp mesh face
3428 dumpFunction(SMESH_Comment("makeTmpFacesOnEdges")<<data._index);
3429 for ( size_t i = 0; i < data._edges.size(); ++i )
3431 _LayerEdge* edge = data._edges[i];
3432 if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
3433 const SMDS_MeshNode* tgt1 = edge->_nodes.back();
3434 for ( int j = 0; j < 2; ++j ) // loop on _2NearEdges
3436 const SMDS_MeshNode* tgt2 = edge->_2neibors->tgtNode(j);
3437 pair< set< SMESH_TLink >::iterator, bool > link_isnew =
3438 extrudedLinks.insert( SMESH_TLink( tgt1, tgt2 ));
3439 if ( !link_isnew.second )
3441 extrudedLinks.erase( link_isnew.first );
3442 continue; // already extruded and will no more encounter
3444 // a _LayerEdge containg tgt2
3445 _LayerEdge* neiborEdge = edge->_2neibors->_edges[j];
3447 _TmpMeshFaceOnEdge* f = new _TmpMeshFaceOnEdge( edge, neiborEdge, --_tmpFaceID );
3448 tmpFaces.push_back( f );
3450 dumpCmd(SMESH_Comment("mesh.AddFace([ ")
3451 <<f->_nn[0]->GetID()<<", "<<f->_nn[1]->GetID()<<", "
3452 <<f->_nn[2]->GetID()<<", "<<f->_nn[3]->GetID()<<" ])");
3457 // Check if _LayerEdge's based on EDGE's intersects tmpFaces.
3458 // Perform two loops on _LayerEdge on EDGE's:
3459 // 1) to find and fix intersection
3460 // 2) to check that no new intersection appears as result of 1)
3462 SMDS_ElemIteratorPtr fIt( new SMDS_ElementVectorIterator( tmpFaces.begin(),
3464 auto_ptr<SMESH_ElementSearcher> searcher
3465 ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), fIt ));
3467 // 1) Find intersections
3469 const SMDS_MeshElement* face;
3470 typedef map< _LayerEdge*, set< _LayerEdge*, _LayerEdgeCmp >, _LayerEdgeCmp > TLEdge2LEdgeSet;
3471 TLEdge2LEdgeSet edge2CloseEdge;
3473 const double eps = data._epsilon * data._epsilon;
3474 for ( size_t i = 0; i < data._edges.size(); ++i )
3476 _LayerEdge* edge = data._edges[i];
3477 if (( !edge->IsOnEdge() ) &&
3478 ( edge->_sWOL.IsNull() || edge->_sWOL.ShapeType() != TopAbs_FACE ))
3480 if ( edge->FindIntersection( *searcher, dist, eps, &face ))
3482 const _TmpMeshFaceOnEdge* f = (const _TmpMeshFaceOnEdge*) face;
3483 set< _LayerEdge*, _LayerEdgeCmp > & ee = edge2CloseEdge[ edge ];
3484 ee.insert( f->_le1 );
3485 ee.insert( f->_le2 );
3486 if ( f->_le1->IsOnEdge() && f->_le1->_sWOL.IsNull() )
3487 edge2CloseEdge[ f->_le1 ].insert( edge );
3488 if ( f->_le2->IsOnEdge() && f->_le2->_sWOL.IsNull() )
3489 edge2CloseEdge[ f->_le2 ].insert( edge );
3493 // Set _LayerEdge._normal
3495 if ( !edge2CloseEdge.empty() )
3497 dumpFunction(SMESH_Comment("updateNormals")<<data._index);
3499 set< TGeomID > shapesToSmooth;
3501 // vector to store new _normal and _cosin for each edge in edge2CloseEdge
3502 vector< pair< _LayerEdge*, _LayerEdge > > edge2newEdge( edge2CloseEdge.size() );
3504 TLEdge2LEdgeSet::iterator e2ee = edge2CloseEdge.begin();
3505 for ( size_t iE = 0; e2ee != edge2CloseEdge.end(); ++e2ee, ++iE )
3507 _LayerEdge* edge1 = e2ee->first;
3508 _LayerEdge* edge2 = 0;
3509 set< _LayerEdge*, _LayerEdgeCmp >& ee = e2ee->second;
3511 edge2newEdge[ iE ].first = NULL;
3513 // find EDGEs the edges reside
3514 // TopoDS_Edge E1, E2;
3515 // TopoDS_Shape S = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
3516 // if ( S.ShapeType() != TopAbs_EDGE )
3517 // continue; // TODO: find EDGE by VERTEX
3518 // E1 = TopoDS::Edge( S );
3519 set< _LayerEdge*, _LayerEdgeCmp >::iterator eIt = ee.begin();
3520 for ( ; !edge2 && eIt != ee.end(); ++eIt )
3522 if ( edge1->_sWOL == (*eIt)->_sWOL )
3525 if ( !edge2 ) continue;
3527 edge2newEdge[ iE ].first = edge1;
3528 _LayerEdge& newEdge = edge2newEdge[ iE ].second;
3529 // while ( E2.IsNull() && eIt != ee.end())
3531 // _LayerEdge* e2 = *eIt++;
3532 // TopoDS_Shape S = helper.GetSubShapeByNode( e2->_nodes[0], getMeshDS() );
3533 // if ( S.ShapeType() == TopAbs_EDGE )
3534 // E2 = TopoDS::Edge( S ), edge2 = e2;
3536 // if ( E2.IsNull() ) continue; // TODO: find EDGE by VERTEX
3538 // find 3 FACEs sharing 2 EDGEs
3540 // TopoDS_Face FF1[2], FF2[2];
3541 // PShapeIteratorPtr fIt = helper.GetAncestors(E1, *_mesh, TopAbs_FACE);
3542 // while ( fIt->more() && FF1[1].IsNull() )
3544 // const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
3545 // if ( helper.IsSubShape( *F, data._solid))
3546 // FF1[ FF1[0].IsNull() ? 0 : 1 ] = *F;
3548 // fIt = helper.GetAncestors(E2, *_mesh, TopAbs_FACE);
3549 // while ( fIt->more() && FF2[1].IsNull())
3551 // const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
3552 // if ( helper.IsSubShape( *F, data._solid))
3553 // FF2[ FF2[0].IsNull() ? 0 : 1 ] = *F;
3555 // // exclude a FACE common to E1 and E2 (put it to FFn[1] )
3556 // if ( FF1[0].IsSame( FF2[0]) || FF1[0].IsSame( FF2[1]))
3557 // std::swap( FF1[0], FF1[1] );
3558 // if ( FF2[0].IsSame( FF1[0]) )
3559 // std::swap( FF2[0], FF2[1] );
3560 // if ( FF1[0].IsNull() || FF2[0].IsNull() )
3563 // get a new normal for edge1
3565 gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal;
3566 // if ( edge1->_cosin < 0 )
3567 // dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized();
3568 // if ( edge2->_cosin < 0 )
3569 // dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized();
3571 double cos1 = Abs( edge1->_cosin ), cos2 = Abs( edge2->_cosin );
3572 double wgt1 = ( cos1 + 0.001 ) / ( cos1 + cos2 + 0.002 );
3573 double wgt2 = ( cos2 + 0.001 ) / ( cos1 + cos2 + 0.002 );
3574 newEdge._normal = ( wgt1 * dir1 + wgt2 * dir2 ).XYZ();
3575 newEdge._normal.Normalize();
3577 // cout << edge1->_nodes[0]->GetID() << " "
3578 // << edge2->_nodes[0]->GetID() << " NORM: "
3579 // << newEdge._normal.X() << ", " << newEdge._normal.Y() << ", " << newEdge._normal.Z() << endl;
3582 if ( cos1 < theMinSmoothCosin )
3584 newEdge._cosin = edge2->_cosin;
3586 else if ( cos2 > theMinSmoothCosin ) // both cos1 and cos2 > theMinSmoothCosin
3588 // gp_Vec dirInFace;
3589 // if ( edge1->_cosin < 0 )
3590 // dirInFace = dir1;
3592 // dirInFace = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
3593 // double angle = dirInFace.Angle( edge1->_normal ); // [0,PI]
3594 // edge1->SetCosin( Cos( angle ));
3595 //newEdge._cosin = 0; // ???????????
3596 newEdge._cosin = ( wgt1 * cos1 + wgt2 * cos2 ) * edge1->_cosin / cos1;
3600 newEdge._cosin = edge1->_cosin;
3603 // find shapes that need smoothing due to change of _normal
3604 if ( edge1->_cosin < theMinSmoothCosin &&
3605 newEdge._cosin > theMinSmoothCosin )
3607 if ( edge1->_sWOL.IsNull() )
3609 SMDS_ElemIteratorPtr fIt = edge1->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
3610 while ( fIt->more() )
3611 shapesToSmooth.insert( fIt->next()->getshapeId() );
3612 //limitStepSize( data, fIt->next(), edge1->_cosin ); // too late
3614 else // edge1 inflates along a FACE
3616 TopoDS_Shape V = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
3617 PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE );
3618 while ( const TopoDS_Shape* E = eIt->next() )
3620 if ( !helper.IsSubShape( *E, /*FACE=*/edge1->_sWOL ))
3622 gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V ));
3623 double angle = edgeDir.Angle( newEdge._normal ); // [0,PI]
3624 if ( angle < M_PI / 2 )
3625 shapesToSmooth.insert( getMeshDS()->ShapeToIndex( *E ));
3631 data.AddShapesToSmooth( shapesToSmooth );
3633 // Update data of edges depending on a new _normal
3635 for ( size_t iE = 0; iE < edge2newEdge.size(); ++iE )
3637 _LayerEdge* edge1 = edge2newEdge[ iE ].first;
3638 _LayerEdge& newEdge = edge2newEdge[ iE ].second;
3639 if ( !edge1 ) continue;
3641 edge1->_normal = newEdge._normal;
3642 edge1->SetCosin( newEdge._cosin );
3643 edge1->InvalidateStep( 1 );
3645 edge1->SetNewLength( data._stepSize, helper );
3646 if ( edge1->IsOnEdge() )
3648 const SMDS_MeshNode * n1 = edge1->_2neibors->srcNode(0);
3649 const SMDS_MeshNode * n2 = edge1->_2neibors->srcNode(1);
3650 edge1->SetDataByNeighbors( n1, n2, helper );
3653 // Update normals and other dependent data of not intersecting _LayerEdge's
3654 // neighboring the intersecting ones
3656 if ( !edge1->_2neibors )
3658 for ( int j = 0; j < 2; ++j ) // loop on 2 neighbors
3660 _LayerEdge* neighbor = edge1->_2neibors->_edges[j];
3661 if ( edge2CloseEdge.count ( neighbor ))
3662 continue; // j-th neighbor is also intersected
3663 _LayerEdge* prevEdge = edge1;
3664 const int nbSteps = 10;
3665 for ( int step = nbSteps; step; --step ) // step from edge1 in j-th direction
3667 if ( !neighbor->_2neibors )
3668 break; // neighbor is on VERTEX
3670 _LayerEdge* nextEdge = neighbor->_2neibors->_edges[iNext];
3671 if ( nextEdge == prevEdge )
3672 nextEdge = neighbor->_2neibors->_edges[ ++iNext ];
3673 double r = double(step-1)/nbSteps;
3674 if ( !nextEdge->_2neibors )
3677 gp_XYZ newNorm = prevEdge->_normal * r + nextEdge->_normal * (1-r);
3678 newNorm.Normalize();
3680 neighbor->_normal = newNorm;
3681 neighbor->SetCosin( prevEdge->_cosin * r + nextEdge->_cosin * (1-r) );
3682 neighbor->SetDataByNeighbors( prevEdge->_nodes[0], nextEdge->_nodes[0], helper );
3684 neighbor->InvalidateStep( 1 );
3686 neighbor->SetNewLength( data._stepSize, helper );
3688 // goto the next neighbor
3689 prevEdge = neighbor;
3690 neighbor = nextEdge;
3696 // 2) Check absence of intersections
3699 for ( size_t i = 0 ; i < tmpFaces.size(); ++i )
3705 //================================================================================
3707 * \brief Modify normals of _LayerEdge's on _ConvexFace's
3709 //================================================================================
3711 bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData& data,
3712 SMESH_MesherHelper& helper,
3715 SMESHDS_Mesh* meshDS = helper.GetMeshDS();
3718 map< TGeomID, _ConvexFace >::iterator id2face = data._convexFaces.begin();
3719 for ( ; id2face != data._convexFaces.end(); ++id2face )
3721 _ConvexFace & convFace = (*id2face).second;
3722 if ( convFace._normalsFixed )
3723 continue; // already fixed
3724 if ( convFace.CheckPrisms() )
3725 continue; // nothing to fix
3727 convFace._normalsFixed = true;
3729 BRepAdaptor_Surface surface ( convFace._face, false );
3730 BRepLProp_SLProps surfProp( surface, 2, 1e-6 );
3732 // check if the convex FACE is of spherical shape
3734 Bnd_B3d centersBox; // bbox of centers of curvature of _LayerEdge's on VERTEXes
3739 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.begin();
3740 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3742 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3744 if ( meshDS->IndexToShape( id2end->first ).ShapeType() == TopAbs_VERTEX )
3746 _LayerEdge* ledge = data._edges[ iBeg ];
3747 if ( convFace.GetCenterOfCurvature( ledge, surfProp, helper, center ))
3748 centersBox.Add( center );
3750 for ( ; iBeg < iEnd; ++iBeg )
3751 nodesBox.Add( SMESH_TNodeXYZ( data._edges[ iBeg ]->_nodes[0] ));
3753 if ( centersBox.IsVoid() )
3755 debugMsg( "Error: centersBox.IsVoid()" );
3758 const bool isSpherical =
3759 ( centersBox.SquareExtent() < 1e-6 * nodesBox.SquareExtent() );
3761 int nbEdges = helper.Count( convFace._face, TopAbs_EDGE, /*ignoreSame=*/false );
3762 vector < _CentralCurveOnEdge > centerCurves( nbEdges );
3766 // set _LayerEdge::_normal as average of all normals
3768 // WARNING: different density of nodes on EDGEs is not taken into account that
3769 // can lead to an improper new normal
3771 gp_XYZ avgNormal( 0,0,0 );
3773 id2end = convFace._subIdToEdgeEnd.begin();
3774 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3776 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3777 // set data of _CentralCurveOnEdge
3778 const TopoDS_Shape& S = meshDS->IndexToShape( id2end->first );
3779 if ( S.ShapeType() == TopAbs_EDGE )
3781 _CentralCurveOnEdge& ceCurve = centerCurves[ nbEdges++ ];
3782 ceCurve.SetShapes( TopoDS::Edge(S), convFace, data, helper );
3783 if ( !data._edges[ iBeg ]->_sWOL.IsNull() )
3784 ceCurve._adjFace.Nullify();
3786 ceCurve._ledges.insert( ceCurve._ledges.end(),
3787 &data._edges[ iBeg ], &data._edges[ iEnd ]);
3789 // summarize normals
3790 for ( ; iBeg < iEnd; ++iBeg )
3791 avgNormal += data._edges[ iBeg ]->_normal;
3793 double normSize = avgNormal.SquareModulus();
3794 if ( normSize < 1e-200 )
3796 debugMsg( "updateNormalsOfConvexFaces(): zero avgNormal" );
3799 avgNormal /= Sqrt( normSize );
3801 // compute new _LayerEdge::_cosin on EDGEs
3802 double avgCosin = 0;
3805 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3807 _CentralCurveOnEdge& ceCurve = centerCurves[ iE ];
3808 if ( ceCurve._adjFace.IsNull() )
3810 for ( size_t iLE = 0; iLE < ceCurve._ledges.size(); ++iLE )
3812 const SMDS_MeshNode* node = ceCurve._ledges[ iLE ]->_nodes[0];
3813 inFaceDir = getFaceDir( ceCurve._adjFace, ceCurve._edge, node, helper, isOK );
3816 double angle = inFaceDir.Angle( avgNormal ); // [0,PI]
3817 ceCurve._ledges[ iLE ]->_cosin = Cos( angle );
3818 avgCosin += ceCurve._ledges[ iLE ]->_cosin;
3824 avgCosin /= nbCosin;
3826 // set _LayerEdge::_normal = avgNormal
3827 id2end = convFace._subIdToEdgeEnd.begin();
3828 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3830 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3831 const TopoDS_Shape& S = meshDS->IndexToShape( id2end->first );
3832 if ( S.ShapeType() != TopAbs_EDGE )
3833 for ( int i = iBeg; i < iEnd; ++i )
3834 data._edges[ i ]->_cosin = avgCosin;
3836 for ( ; iBeg < iEnd; ++iBeg )
3837 data._edges[ iBeg ]->_normal = avgNormal;
3840 else // if ( isSpherical )
3842 // We suppose that centers of curvature at all points of the FACE
3843 // lie on some curve, let's call it "central curve". For all _LayerEdge's
3844 // having a common center of curvature we define the same new normal
3845 // as a sum of normals of _LayerEdge's on EDGEs among them.
3847 // get all centers of curvature for each EDGE
3849 helper.SetSubShape( convFace._face );
3850 _LayerEdge* vertexLEdges[2], **edgeLEdge, **edgeLEdgeEnd;
3852 TopExp_Explorer edgeExp( convFace._face, TopAbs_EDGE );
3853 for ( int iE = 0; edgeExp.More(); edgeExp.Next(), ++iE )
3855 const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
3857 // set adjacent FACE
3858 centerCurves[ iE ].SetShapes( edge, convFace, data, helper );
3860 // get _LayerEdge's of the EDGE
3861 TGeomID edgeID = meshDS->ShapeToIndex( edge );
3862 id2end = convFace._subIdToEdgeEnd.find( edgeID );
3863 if ( id2end == convFace._subIdToEdgeEnd.end() )
3865 // no _LayerEdge's on EDGE, use _LayerEdge's on VERTEXes
3866 for ( int iV = 0; iV < 2; ++iV )
3868 TopoDS_Vertex v = helper.IthVertex( iV, edge );
3869 TGeomID vID = meshDS->ShapeToIndex( v );
3870 int end = convFace._subIdToEdgeEnd[ vID ];
3871 int iBeg = end > 0 ? data._endEdgeOnShape[ end-1 ] : 0;
3872 vertexLEdges[ iV ] = data._edges[ iBeg ];
3874 edgeLEdge = &vertexLEdges[0];
3875 edgeLEdgeEnd = edgeLEdge + 2;
3877 centerCurves[ iE ]._adjFace.Nullify();
3881 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3882 if ( id2end->second >= data._nbShapesToSmooth )
3883 data.SortOnEdge( edge, iBeg, iEnd, helper );
3884 edgeLEdge = &data._edges[ iBeg ];
3885 edgeLEdgeEnd = edgeLEdge + iEnd - iBeg;
3886 vertexLEdges[0] = data._edges[ iBeg ]->_2neibors->_edges[0];
3887 vertexLEdges[1] = data._edges[ iEnd-1 ]->_2neibors->_edges[1];
3889 if ( ! data._edges[ iBeg ]->_sWOL.IsNull() )
3890 centerCurves[ iE ]._adjFace.Nullify();
3893 // Get curvature centers
3897 if ( edgeLEdge[0]->IsOnEdge() &&
3898 convFace.GetCenterOfCurvature( vertexLEdges[0], surfProp, helper, center ))
3900 centerCurves[ iE ].Append( center, vertexLEdges[0] );
3901 centersBox.Add( center );
3903 for ( ; edgeLEdge < edgeLEdgeEnd; ++edgeLEdge )
3904 if ( convFace.GetCenterOfCurvature( *edgeLEdge, surfProp, helper, center ))
3905 { // EDGE or VERTEXes
3906 centerCurves[ iE ].Append( center, *edgeLEdge );
3907 centersBox.Add( center );
3909 if ( edgeLEdge[-1]->IsOnEdge() &&
3910 convFace.GetCenterOfCurvature( vertexLEdges[1], surfProp, helper, center ))
3912 centerCurves[ iE ].Append( center, vertexLEdges[1] );
3913 centersBox.Add( center );
3915 centerCurves[ iE ]._isDegenerated =
3916 ( centersBox.IsVoid() || centersBox.SquareExtent() < 1e-6 * nodesBox.SquareExtent() );
3918 } // loop on EDGES of convFace._face to set up data of centerCurves
3920 // Compute new normals for _LayerEdge's on EDGEs
3922 double avgCosin = 0;
3925 for ( size_t iE1 = 0; iE1 < centerCurves.size(); ++iE1 )
3927 _CentralCurveOnEdge& ceCurve = centerCurves[ iE1 ];
3928 if ( ceCurve._isDegenerated )
3930 const vector< gp_Pnt >& centers = ceCurve._curvaCenters;
3931 vector< gp_XYZ > & newNormals = ceCurve._normals;
3932 for ( size_t iC1 = 0; iC1 < centers.size(); ++iC1 )
3935 for ( size_t iE2 = 0; iE2 < centerCurves.size() && !isOK; ++iE2 )
3938 isOK = centerCurves[ iE2 ].FindNewNormal( centers[ iC1 ], newNormals[ iC1 ]);
3940 if ( isOK && !ceCurve._adjFace.IsNull() )
3942 // compute new _LayerEdge::_cosin
3943 const SMDS_MeshNode* node = ceCurve._ledges[ iC1 ]->_nodes[0];
3944 inFaceDir = getFaceDir( ceCurve._adjFace, ceCurve._edge, node, helper, isOK );
3947 double angle = inFaceDir.Angle( newNormals[ iC1 ] ); // [0,PI]
3948 ceCurve._ledges[ iC1 ]->_cosin = Cos( angle );
3949 avgCosin += ceCurve._ledges[ iC1 ]->_cosin;
3955 // set new normals to _LayerEdge's of NOT degenerated central curves
3956 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3958 if ( centerCurves[ iE ]._isDegenerated )
3960 for ( size_t iLE = 0; iLE < centerCurves[ iE ]._ledges.size(); ++iLE )
3961 centerCurves[ iE ]._ledges[ iLE ]->_normal = centerCurves[ iE ]._normals[ iLE ];
3963 // set new normals to _LayerEdge's of degenerated central curves
3964 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3966 if ( !centerCurves[ iE ]._isDegenerated ||
3967 centerCurves[ iE ]._ledges.size() < 3 )
3969 // new normal is an average of new normals at VERTEXes that
3970 // was computed on non-degenerated _CentralCurveOnEdge's
3971 gp_XYZ newNorm = ( centerCurves[ iE ]._ledges.front()->_normal +
3972 centerCurves[ iE ]._ledges.back ()->_normal );
3973 double sz = newNorm.Modulus();
3977 double newCosin = ( 0.5 * centerCurves[ iE ]._ledges.front()->_cosin +
3978 0.5 * centerCurves[ iE ]._ledges.back ()->_cosin );
3979 for ( size_t iLE = 1, nb = centerCurves[ iE ]._ledges.size() - 1; iLE < nb; ++iLE )
3981 centerCurves[ iE ]._ledges[ iLE ]->_normal = newNorm;
3982 centerCurves[ iE ]._ledges[ iLE ]->_cosin = newCosin;
3986 // Find new normals for _LayerEdge's based on FACE
3989 avgCosin /= nbCosin;
3990 const TGeomID faceID = meshDS->ShapeToIndex( convFace._face );
3991 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.find( faceID );
3992 if ( id2end != convFace._subIdToEdgeEnd.end() )
3996 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3997 for ( ; iBeg < iEnd; ++iBeg )
3999 _LayerEdge* ledge = data._edges[ iBeg ];
4000 if ( !convFace.GetCenterOfCurvature( ledge, surfProp, helper, center ))
4002 for ( size_t i = 0; i < centerCurves.size(); ++i, ++iE )
4004 iE = iE % centerCurves.size();
4005 if ( centerCurves[ iE ]._isDegenerated )
4007 newNorm.SetCoord( 0,0,0 );
4008 if ( centerCurves[ iE ].FindNewNormal( center, newNorm ))
4010 ledge->_normal = newNorm;
4011 ledge->_cosin = avgCosin;
4018 } // not a quasi-spherical FACE
4020 // Update _LayerEdge's data according to a new normal
4022 dumpFunction(SMESH_Comment("updateNormalsOfConvexFaces")<<data._index
4023 <<"_F"<<meshDS->ShapeToIndex( convFace._face ));
4025 id2end = convFace._subIdToEdgeEnd.begin();
4026 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
4028 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
4029 for ( ; iBeg < iEnd; ++iBeg )
4031 _LayerEdge* & ledge = data._edges[ iBeg ];
4032 double len = ledge->_len;
4033 ledge->InvalidateStep( stepNb + 1, /*restoreLength=*/true );
4034 ledge->SetCosin( ledge->_cosin );
4035 ledge->SetNewLength( len, helper );
4038 } // loop on sub-shapes of convFace._face
4040 // Find FACEs adjacent to convFace._face that got necessity to smooth
4041 // as a result of normals modification
4043 set< TGeomID > adjFacesToSmooth;
4044 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
4046 if ( centerCurves[ iE ]._adjFace.IsNull() ||
4047 centerCurves[ iE ]._adjFaceToSmooth )
4049 for ( size_t iLE = 0; iLE < centerCurves[ iE ]._ledges.size(); ++iLE )
4051 if ( centerCurves[ iE ]._ledges[ iLE ]->_cosin > theMinSmoothCosin )
4053 adjFacesToSmooth.insert( meshDS->ShapeToIndex( centerCurves[ iE ]._adjFace ));
4058 data.AddShapesToSmooth( adjFacesToSmooth );
4063 } // loop on data._convexFaces
4068 //================================================================================
4070 * \brief Finds a center of curvature of a surface at a _LayerEdge
4072 //================================================================================
4074 bool _ConvexFace::GetCenterOfCurvature( _LayerEdge* ledge,
4075 BRepLProp_SLProps& surfProp,
4076 SMESH_MesherHelper& helper,
4077 gp_Pnt & center ) const
4079 gp_XY uv = helper.GetNodeUV( _face, ledge->_nodes[0] );
4080 surfProp.SetParameters( uv.X(), uv.Y() );
4081 if ( !surfProp.IsCurvatureDefined() )
4084 const double oriFactor = ( _face.Orientation() == TopAbs_REVERSED ? +1. : -1. );
4085 double surfCurvatureMax = surfProp.MaxCurvature() * oriFactor;
4086 double surfCurvatureMin = surfProp.MinCurvature() * oriFactor;
4087 if ( surfCurvatureMin > surfCurvatureMax )
4088 center = surfProp.Value().Translated( surfProp.Normal().XYZ() / surfCurvatureMin * oriFactor );
4090 center = surfProp.Value().Translated( surfProp.Normal().XYZ() / surfCurvatureMax * oriFactor );
4095 //================================================================================
4097 * \brief Check that prisms are not distorted
4099 //================================================================================
4101 bool _ConvexFace::CheckPrisms() const
4103 for ( size_t i = 0; i < _simplexTestEdges.size(); ++i )
4105 const _LayerEdge* edge = _simplexTestEdges[i];
4106 SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
4107 for ( size_t j = 0; j < edge->_simplices.size(); ++j )
4108 if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
4110 debugMsg( "Bad simplex of _simplexTestEdges ("
4111 << " "<< edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
4112 << " "<< edge->_simplices[j]._nPrev->GetID()
4113 << " "<< edge->_simplices[j]._nNext->GetID() << " )" );
4120 //================================================================================
4122 * \brief Try to compute a new normal by interpolating normals of _LayerEdge's
4123 * stored in this _CentralCurveOnEdge.
4124 * \param [in] center - curvature center of a point of another _CentralCurveOnEdge.
4125 * \param [in,out] newNormal - current normal at this point, to be redefined
4126 * \return bool - true if succeeded.
4128 //================================================================================
4130 bool _CentralCurveOnEdge::FindNewNormal( const gp_Pnt& center, gp_XYZ& newNormal )
4132 if ( this->_isDegenerated )
4135 // find two centers the given one lies between
4137 for ( size_t i = 0, nb = _curvaCenters.size()-1; i < nb; ++i )
4139 double sl2 = 1.001 * _segLength2[ i ];
4141 double d1 = center.SquareDistance( _curvaCenters[ i ]);
4145 double d2 = center.SquareDistance( _curvaCenters[ i+1 ]);
4146 if ( d2 > sl2 || d2 + d1 < 1e-100 )
4151 double r = d1 / ( d1 + d2 );
4152 gp_XYZ norm = (( 1. - r ) * _ledges[ i ]->_normal +
4153 ( r ) * _ledges[ i+1 ]->_normal );
4157 double sz = newNormal.Modulus();
4166 //================================================================================
4168 * \brief Set shape members
4170 //================================================================================
4172 void _CentralCurveOnEdge::SetShapes( const TopoDS_Edge& edge,
4173 const _ConvexFace& convFace,
4174 const _SolidData& data,
4175 SMESH_MesherHelper& helper)
4179 PShapeIteratorPtr fIt = helper.GetAncestors( edge, *helper.GetMesh(), TopAbs_FACE );
4180 while ( const TopoDS_Shape* F = fIt->next())
4181 if ( !convFace._face.IsSame( *F ))
4183 _adjFace = TopoDS::Face( *F );
4184 _adjFaceToSmooth = false;
4185 // _adjFace already in a smoothing queue ?
4187 TGeomID adjFaceID = helper.GetMeshDS()->ShapeToIndex( *F );
4188 if ( data.GetShapeEdges( adjFaceID, end ))
4189 _adjFaceToSmooth = ( end < data._nbShapesToSmooth );
4194 //================================================================================
4196 * \brief Looks for intersection of it's last segment with faces
4197 * \param distance - returns shortest distance from the last node to intersection
4199 //================================================================================
4201 bool _LayerEdge::FindIntersection( SMESH_ElementSearcher& searcher,
4203 const double& epsilon,
4204 const SMDS_MeshElement** face)
4206 vector< const SMDS_MeshElement* > suspectFaces;
4208 gp_Ax1 lastSegment = LastSegment(segLen);
4209 searcher.GetElementsNearLine( lastSegment, SMDSAbs_Face, suspectFaces );
4211 bool segmentIntersected = false;
4212 distance = Precision::Infinite();
4213 int iFace = -1; // intersected face
4214 for ( size_t j = 0 ; j < suspectFaces.size() /*&& !segmentIntersected*/; ++j )
4216 const SMDS_MeshElement* face = suspectFaces[j];
4217 if ( face->GetNodeIndex( _nodes.back() ) >= 0 ||
4218 face->GetNodeIndex( _nodes[0] ) >= 0 )
4219 continue; // face sharing _LayerEdge node
4220 const int nbNodes = face->NbCornerNodes();
4221 bool intFound = false;
4223 SMDS_MeshElement::iterator nIt = face->begin_nodes();
4226 intFound = SegTriaInter( lastSegment, *nIt++, *nIt++, *nIt++, dist, epsilon );
4230 const SMDS_MeshNode* tria[3];
4233 for ( int n2 = 2; n2 < nbNodes && !intFound; ++n2 )
4236 intFound = SegTriaInter(lastSegment, tria[0], tria[1], tria[2], dist, epsilon );
4242 if ( dist < segLen*(1.01) && dist > -(_len*_lenFactor-segLen) )
4243 segmentIntersected = true;
4244 if ( distance > dist )
4245 distance = dist, iFace = j;
4248 if ( iFace != -1 && face ) *face = suspectFaces[iFace];
4250 if ( segmentIntersected )
4253 SMDS_MeshElement::iterator nIt = suspectFaces[iFace]->begin_nodes();
4254 gp_XYZ intP( lastSegment.Location().XYZ() + lastSegment.Direction().XYZ() * distance );
4255 cout << "nodes: tgt " << _nodes.back()->GetID() << " src " << _nodes[0]->GetID()
4256 << ", intersection with face ("
4257 << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
4258 << ") at point (" << intP.X() << ", " << intP.Y() << ", " << intP.Z()
4259 << ") distance = " << distance - segLen<< endl;
4265 return segmentIntersected;
4268 //================================================================================
4270 * \brief Returns size and direction of the last segment
4272 //================================================================================
4274 gp_Ax1 _LayerEdge::LastSegment(double& segLen) const
4276 // find two non-coincident positions
4277 gp_XYZ orig = _pos.back();
4279 int iPrev = _pos.size() - 2;
4280 while ( iPrev >= 0 )
4282 dir = orig - _pos[iPrev];
4283 if ( dir.SquareModulus() > 1e-100 )
4293 segDir.SetLocation( SMESH_TNodeXYZ( _nodes[0] ));
4294 segDir.SetDirection( _normal );
4299 gp_Pnt pPrev = _pos[ iPrev ];
4300 if ( !_sWOL.IsNull() )
4302 TopLoc_Location loc;
4303 if ( _sWOL.ShapeType() == TopAbs_EDGE )
4306 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
4307 pPrev = curve->Value( pPrev.X() ).Transformed( loc );
4311 Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
4312 pPrev = surface->Value( pPrev.X(), pPrev.Y() ).Transformed( loc );
4314 dir = SMESH_TNodeXYZ( _nodes.back() ) - pPrev.XYZ();
4316 segDir.SetLocation( pPrev );
4317 segDir.SetDirection( dir );
4318 segLen = dir.Modulus();
4324 //================================================================================
4326 * \brief Return the last position of the target node on a FACE.
4327 * \param [in] F - the FACE this _LayerEdge is inflated along
4328 * \return gp_XY - result UV
4330 //================================================================================
4332 gp_XY _LayerEdge::LastUV( const TopoDS_Face& F ) const
4334 if ( F.IsSame( _sWOL )) // F is my FACE
4335 return gp_XY( _pos.back().X(), _pos.back().Y() );
4337 if ( _sWOL.IsNull() || _sWOL.ShapeType() != TopAbs_EDGE ) // wrong call
4338 return gp_XY( 1e100, 1e100 );
4340 // _sWOL is EDGE of F; _pos.back().X() is the last U on the EDGE
4341 double f, l, u = _pos.back().X();
4342 Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( TopoDS::Edge(_sWOL), F, f,l);
4343 if ( !C2d.IsNull() && f <= u && u <= l )
4344 return C2d->Value( u ).XY();
4346 return gp_XY( 1e100, 1e100 );
4349 //================================================================================
4351 * \brief Test intersection of the last segment with a given triangle
4352 * using Moller-Trumbore algorithm
4353 * Intersection is detected if distance to intersection is less than _LayerEdge._len
4355 //================================================================================
4357 bool _LayerEdge::SegTriaInter( const gp_Ax1& lastSegment,
4358 const SMDS_MeshNode* n0,
4359 const SMDS_MeshNode* n1,
4360 const SMDS_MeshNode* n2,
4362 const double& EPSILON) const
4364 //const double EPSILON = 1e-6;
4366 gp_XYZ orig = lastSegment.Location().XYZ();
4367 gp_XYZ dir = lastSegment.Direction().XYZ();
4369 SMESH_TNodeXYZ vert0( n0 );
4370 SMESH_TNodeXYZ vert1( n1 );
4371 SMESH_TNodeXYZ vert2( n2 );
4373 /* calculate distance from vert0 to ray origin */
4374 gp_XYZ tvec = orig - vert0;
4376 //if ( tvec * dir > EPSILON )
4377 // intersected face is at back side of the temporary face this _LayerEdge belongs to
4380 gp_XYZ edge1 = vert1 - vert0;
4381 gp_XYZ edge2 = vert2 - vert0;
4383 /* begin calculating determinant - also used to calculate U parameter */
4384 gp_XYZ pvec = dir ^ edge2;
4386 /* if determinant is near zero, ray lies in plane of triangle */
4387 double det = edge1 * pvec;
4389 if (det > -EPSILON && det < EPSILON)
4391 double inv_det = 1.0 / det;
4393 /* calculate U parameter and test bounds */
4394 double u = ( tvec * pvec ) * inv_det;
4395 //if (u < 0.0 || u > 1.0)
4396 if (u < -EPSILON || u > 1.0 + EPSILON)
4399 /* prepare to test V parameter */
4400 gp_XYZ qvec = tvec ^ edge1;
4402 /* calculate V parameter and test bounds */
4403 double v = (dir * qvec) * inv_det;
4404 //if ( v < 0.0 || u + v > 1.0 )
4405 if ( v < -EPSILON || u + v > 1.0 + EPSILON)
4408 /* calculate t, ray intersects triangle */
4409 t = (edge2 * qvec) * inv_det;
4415 //================================================================================
4417 * \brief Perform smooth of _LayerEdge's based on EDGE's
4418 * \retval bool - true if node has been moved
4420 //================================================================================
4422 bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface,
4423 const TopoDS_Face& F,
4424 SMESH_MesherHelper& helper)
4426 ASSERT( IsOnEdge() );
4428 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( _nodes.back() );
4429 SMESH_TNodeXYZ oldPos( tgtNode );
4430 double dist01, distNewOld;
4432 SMESH_TNodeXYZ p0( _2neibors->tgtNode(0));
4433 SMESH_TNodeXYZ p1( _2neibors->tgtNode(1));
4434 dist01 = p0.Distance( _2neibors->tgtNode(1) );
4436 gp_Pnt newPos = p0 * _2neibors->_wgt[0] + p1 * _2neibors->_wgt[1];
4437 double lenDelta = 0;
4440 //lenDelta = _curvature->lenDelta( _len );
4441 lenDelta = _curvature->lenDeltaByDist( dist01 );
4442 newPos.ChangeCoord() += _normal * lenDelta;
4445 distNewOld = newPos.Distance( oldPos );
4449 if ( _2neibors->_plnNorm )
4451 // put newPos on the plane defined by source node and _plnNorm
4452 gp_XYZ new2src = SMESH_TNodeXYZ( _nodes[0] ) - newPos.XYZ();
4453 double new2srcProj = (*_2neibors->_plnNorm) * new2src;
4454 newPos.ChangeCoord() += (*_2neibors->_plnNorm) * new2srcProj;
4456 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4457 _pos.back() = newPos.XYZ();
4461 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4462 gp_XY uv( Precision::Infinite(), 0 );
4463 helper.CheckNodeUV( F, tgtNode, uv, 1e-10, /*force=*/true );
4464 _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
4466 newPos = surface->Value( uv.X(), uv.Y() );
4467 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4470 if ( _curvature && lenDelta < 0 )
4472 gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
4473 _len -= prevPos.Distance( oldPos );
4474 _len += prevPos.Distance( newPos );
4476 bool moved = distNewOld > dist01/50;
4478 dumpMove( tgtNode ); // debug
4483 //================================================================================
4485 * \brief Perform laplacian smooth in 3D of nodes inflated from FACE
4486 * \retval bool - true if _tgtNode has been moved
4488 //================================================================================
4490 bool _LayerEdge::Smooth(int& badNb)
4492 if ( _simplices.size() < 2 )
4493 return false; // _LayerEdge inflated along EDGE or FACE
4495 // compute new position for the last _pos
4496 gp_XYZ newPos (0,0,0);
4497 for ( size_t i = 0; i < _simplices.size(); ++i )
4498 newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev );
4499 newPos /= _simplices.size();
4501 const gp_XYZ& curPos ( _pos.back() );
4502 const gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
4505 double delta = _curvature->lenDelta( _len );
4507 newPos += _normal * delta;
4510 double segLen = _normal * ( newPos - prevPos.XYZ() );
4511 if ( segLen + delta > 0 )
4512 newPos += _normal * delta;
4514 // double segLenChange = _normal * ( curPos - newPos );
4515 // newPos += 0.5 * _normal * segLenChange;
4518 // count quality metrics (orientation) of tetras around _tgtNode
4520 for ( size_t i = 0; i < _simplices.size(); ++i )
4521 nbOkBefore += _simplices[i].IsForward( _nodes[0], &curPos );
4524 for ( size_t i = 0; i < _simplices.size(); ++i )
4525 nbOkAfter += _simplices[i].IsForward( _nodes[0], &newPos );
4527 if ( nbOkAfter < nbOkBefore )
4530 SMDS_MeshNode* n = const_cast< SMDS_MeshNode* >( _nodes.back() );
4532 _len -= prevPos.Distance(SMESH_TNodeXYZ( n ));
4533 _len += prevPos.Distance(newPos);
4535 n->setXYZ( newPos.X(), newPos.Y(), newPos.Z());
4536 _pos.back() = newPos;
4538 badNb += _simplices.size() - nbOkAfter;
4545 //================================================================================
4547 * \brief Add a new segment to _LayerEdge during inflation
4549 //================================================================================
4551 void _LayerEdge::SetNewLength( double len, SMESH_MesherHelper& helper )
4553 if ( _len - len > -1e-6 )
4555 _pos.push_back( _pos.back() );
4559 SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
4560 SMESH_TNodeXYZ oldXYZ( n );
4561 gp_XYZ nXYZ = oldXYZ + _normal * ( len - _len ) * _lenFactor;
4562 n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
4564 _pos.push_back( nXYZ );
4566 if ( !_sWOL.IsNull() )
4569 if ( _sWOL.ShapeType() == TopAbs_EDGE )
4571 double u = Precision::Infinite(); // to force projection w/o distance check
4572 helper.CheckNodeU( TopoDS::Edge( _sWOL ), n, u, 1e-10, /*force=*/true, distXYZ );
4573 _pos.back().SetCoord( u, 0, 0 );
4574 if ( _nodes.size() > 1 )
4576 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
4577 pos->SetUParameter( u );
4582 gp_XY uv( Precision::Infinite(), 0 );
4583 helper.CheckNodeUV( TopoDS::Face( _sWOL ), n, uv, 1e-10, /*force=*/true, distXYZ );
4584 _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
4585 if ( _nodes.size() > 1 )
4587 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
4588 pos->SetUParameter( uv.X() );
4589 pos->SetVParameter( uv.Y() );
4592 n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]);
4594 dumpMove( n ); //debug
4597 //================================================================================
4599 * \brief Remove last inflation step
4601 //================================================================================
4603 void _LayerEdge::InvalidateStep( int curStep, bool restoreLength )
4605 if ( _pos.size() > curStep )
4607 if ( restoreLength )
4608 _len -= ( _pos[ curStep-1 ] - _pos.back() ).Modulus();
4610 _pos.resize( curStep );
4611 gp_Pnt nXYZ = _pos.back();
4612 SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
4613 if ( !_sWOL.IsNull() )
4615 TopLoc_Location loc;
4616 if ( _sWOL.ShapeType() == TopAbs_EDGE )
4618 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
4619 pos->SetUParameter( nXYZ.X() );
4621 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
4622 nXYZ = curve->Value( nXYZ.X() ).Transformed( loc );
4626 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
4627 pos->SetUParameter( nXYZ.X() );
4628 pos->SetVParameter( nXYZ.Y() );
4629 Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
4630 nXYZ = surface->Value( nXYZ.X(), nXYZ.Y() ).Transformed( loc );
4633 n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
4638 //================================================================================
4640 * \brief Create layers of prisms
4642 //================================================================================
4644 bool _ViscousBuilder::refine(_SolidData& data)
4646 SMESH_MesherHelper helper( *_mesh );
4647 helper.SetSubShape( data._solid );
4648 helper.SetElementsOnShape(false);
4650 Handle(Geom_Curve) curve;
4651 Handle(Geom_Surface) surface;
4652 TopoDS_Edge geomEdge;
4653 TopoDS_Face geomFace;
4654 TopoDS_Shape prevSWOL;
4655 TopLoc_Location loc;
4659 TGeomID prevBaseId = -1;
4660 TNode2Edge* n2eMap = 0;
4661 TNode2Edge::iterator n2e;
4663 // Create intermediate nodes on each _LayerEdge
4665 for ( size_t i = 0; i < data._edges.size(); ++i )
4667 _LayerEdge& edge = *data._edges[i];
4669 if ( edge._nodes.size() < 2 )
4670 continue; // on _noShrinkShapes
4672 // get accumulated length of segments
4673 vector< double > segLen( edge._pos.size() );
4675 for ( size_t j = 1; j < edge._pos.size(); ++j )
4676 segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus();
4678 // allocate memory for new nodes if it is not yet refined
4679 const SMDS_MeshNode* tgtNode = edge._nodes.back();
4680 if ( edge._nodes.size() == 2 )
4682 edge._nodes.resize( data._hyp->GetNumberLayers() + 1, 0 );
4684 edge._nodes.back() = tgtNode;
4686 // get data of a shrink shape
4687 if ( !edge._sWOL.IsNull() && edge._sWOL != prevSWOL )
4689 isOnEdge = ( edge._sWOL.ShapeType() == TopAbs_EDGE );
4692 geomEdge = TopoDS::Edge( edge._sWOL );
4693 curve = BRep_Tool::Curve( geomEdge, loc, f,l);
4697 geomFace = TopoDS::Face( edge._sWOL );
4698 surface = BRep_Tool::Surface( geomFace, loc );
4700 prevSWOL = edge._sWOL;
4702 // restore shapePos of the last node by already treated _LayerEdge of another _SolidData
4703 const TGeomID baseShapeId = edge._nodes[0]->getshapeId();
4704 if ( baseShapeId != prevBaseId )
4706 map< TGeomID, TNode2Edge* >::iterator s2ne = data._s2neMap.find( baseShapeId );
4707 n2eMap = ( s2ne == data._s2neMap.end() ) ? 0 : n2eMap = s2ne->second;
4708 prevBaseId = baseShapeId;
4710 _LayerEdge* edgeOnSameNode = 0;
4711 if ( n2eMap && (( n2e = n2eMap->find( edge._nodes[0] )) != n2eMap->end() ))
4713 edgeOnSameNode = n2e->second;
4714 const gp_XYZ& otherTgtPos = edgeOnSameNode->_pos.back();
4715 SMDS_PositionPtr lastPos = tgtNode->GetPosition();
4718 SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( lastPos );
4719 epos->SetUParameter( otherTgtPos.X() );
4723 SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( lastPos );
4724 fpos->SetUParameter( otherTgtPos.X() );
4725 fpos->SetVParameter( otherTgtPos.Y() );
4728 // calculate height of the first layer
4730 const double T = segLen.back(); //data._hyp.GetTotalThickness();
4731 const double f = data._hyp->GetStretchFactor();
4732 const int N = data._hyp->GetNumberLayers();
4733 const double fPowN = pow( f, N );
4734 if ( fPowN - 1 <= numeric_limits<double>::min() )
4737 h0 = T * ( f - 1 )/( fPowN - 1 );
4739 const double zeroLen = std::numeric_limits<double>::min();
4741 // create intermediate nodes
4742 double hSum = 0, hi = h0/f;
4744 for ( size_t iStep = 1; iStep < edge._nodes.size(); ++iStep )
4746 // compute an intermediate position
4749 while ( hSum > segLen[iSeg] && iSeg < segLen.size()-1)
4751 int iPrevSeg = iSeg-1;
4752 while ( fabs( segLen[iPrevSeg] - segLen[iSeg]) <= zeroLen && iPrevSeg > 0 )
4754 double r = ( segLen[iSeg] - hSum ) / ( segLen[iSeg] - segLen[iPrevSeg] );
4755 gp_Pnt pos = r * edge._pos[iPrevSeg] + (1-r) * edge._pos[iSeg];
4757 SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >( edge._nodes[ iStep ]);
4758 if ( !edge._sWOL.IsNull() )
4760 // compute XYZ by parameters <pos>
4765 pos = curve->Value( u ).Transformed(loc);
4769 uv.SetCoord( pos.X(), pos.Y() );
4771 pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc);
4774 // create or update the node
4777 node = helper.AddNode( pos.X(), pos.Y(), pos.Z());
4778 if ( !edge._sWOL.IsNull() )
4781 getMeshDS()->SetNodeOnEdge( node, geomEdge, u );
4783 getMeshDS()->SetNodeOnFace( node, geomFace, uv.X(), uv.Y() );
4787 getMeshDS()->SetNodeInVolume( node, helper.GetSubShapeID() );
4792 if ( !edge._sWOL.IsNull() )
4794 // make average pos from new and current parameters
4797 u = 0.5 * ( u + helper.GetNodeU( geomEdge, node ));
4798 pos = curve->Value( u ).Transformed(loc);
4800 SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( node->GetPosition() );
4801 epos->SetUParameter( u );
4805 uv = 0.5 * ( uv + helper.GetNodeUV( geomFace, node ));
4806 pos = surface->Value( uv.X(), uv.Y()).Transformed(loc);
4808 SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( node->GetPosition() );
4809 fpos->SetUParameter( uv.X() );
4810 fpos->SetVParameter( uv.Y() );
4813 node->setXYZ( pos.X(), pos.Y(), pos.Z() );
4815 } // loop on edge._nodes
4817 if ( !edge._sWOL.IsNull() ) // prepare for shrink()
4820 edge._pos.back().SetCoord( u, 0,0);
4822 edge._pos.back().SetCoord( uv.X(), uv.Y() ,0);
4824 if ( edgeOnSameNode )
4825 edgeOnSameNode->_pos.back() = edge._pos.back();
4828 } // loop on data._edges to create nodes
4830 if ( !getMeshDS()->IsEmbeddedMode() )
4831 // Log node movement
4832 for ( size_t i = 0; i < data._edges.size(); ++i )
4834 _LayerEdge& edge = *data._edges[i];
4835 SMESH_TNodeXYZ p ( edge._nodes.back() );
4836 getMeshDS()->MoveNode( p._node, p.X(), p.Y(), p.Z() );
4841 helper.SetElementsOnShape(true);
4843 vector< vector<const SMDS_MeshNode*>* > nnVec;
4844 set< int > degenEdgeInd;
4845 vector<const SMDS_MeshElement*> degenVols;
4847 TopExp_Explorer exp( data._solid, TopAbs_FACE );
4848 for ( ; exp.More(); exp.Next() )
4850 if ( data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
4852 SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
4853 SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
4854 while ( fIt->more() )
4856 const SMDS_MeshElement* face = fIt->next();
4857 const int nbNodes = face->NbCornerNodes();
4858 nnVec.resize( nbNodes );
4859 degenEdgeInd.clear();
4861 SMDS_ElemIteratorPtr nIt = face->nodesIterator();
4862 for ( int iN = 0; iN < nbNodes; ++iN )
4864 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
4865 nnVec[ iN ] = & data._n2eMap[ n ]->_nodes;
4866 if ( nnVec[ iN ]->size() < 2 )
4867 degenEdgeInd.insert( iN );
4869 nbZ = nnVec[ iN ]->size();
4877 switch ( degenEdgeInd.size() )
4881 for ( int iZ = 1; iZ < nbZ; ++iZ )
4882 helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
4883 (*nnVec[0])[iZ], (*nnVec[1])[iZ], (*nnVec[2])[iZ]);
4888 int i2 = *degenEdgeInd.begin();
4889 int i0 = helper.WrapIndex( i2 - 1, nbNodes );
4890 int i1 = helper.WrapIndex( i2 + 1, nbNodes );
4891 for ( int iZ = 1; iZ < nbZ; ++iZ )
4892 helper.AddVolume( (*nnVec[i0])[iZ-1], (*nnVec[i1])[iZ-1],
4893 (*nnVec[i1])[iZ], (*nnVec[i0])[iZ], (*nnVec[i2])[0]);
4898 int i3 = !degenEdgeInd.count(0) ? 0 : !degenEdgeInd.count(1) ? 1 : 2;
4899 for ( int iZ = 1; iZ < nbZ; ++iZ )
4900 helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
4908 switch ( degenEdgeInd.size() )
4912 for ( int iZ = 1; iZ < nbZ; ++iZ )
4913 helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1],
4914 (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1],
4915 (*nnVec[0])[iZ], (*nnVec[1])[iZ],
4916 (*nnVec[2])[iZ], (*nnVec[3])[iZ]);
4921 int i2 = *degenEdgeInd.begin();
4922 int i3 = *degenEdgeInd.rbegin();
4923 bool ok = ( i3 - i2 == 1 );
4924 if ( i2 == 0 && i3 == 3 ) { i2 = 3; i3 = 0; ok = true; }
4925 int i0 = helper.WrapIndex( i3 + 1, nbNodes );
4926 int i1 = helper.WrapIndex( i0 + 1, nbNodes );
4927 for ( int iZ = 1; iZ < nbZ; ++iZ )
4929 const SMDS_MeshElement* vol =
4930 helper.AddVolume( (*nnVec[i3])[0], (*nnVec[i0])[iZ], (*nnVec[i0])[iZ-1],
4931 (*nnVec[i2])[0], (*nnVec[i1])[iZ], (*nnVec[i1])[iZ-1]);
4933 degenVols.push_back( vol );
4937 case 3: // degen HEX
4939 const SMDS_MeshNode* nn[8];
4940 for ( int iZ = 1; iZ < nbZ; ++iZ )
4942 const SMDS_MeshElement* vol =
4943 helper.AddVolume( nnVec[0]->size() > 1 ? (*nnVec[0])[iZ-1] : (*nnVec[0])[0],
4944 nnVec[1]->size() > 1 ? (*nnVec[1])[iZ-1] : (*nnVec[1])[0],
4945 nnVec[2]->size() > 1 ? (*nnVec[2])[iZ-1] : (*nnVec[2])[0],
4946 nnVec[3]->size() > 1 ? (*nnVec[3])[iZ-1] : (*nnVec[3])[0],
4947 nnVec[0]->size() > 1 ? (*nnVec[0])[iZ] : (*nnVec[0])[0],
4948 nnVec[1]->size() > 1 ? (*nnVec[1])[iZ] : (*nnVec[1])[0],
4949 nnVec[2]->size() > 1 ? (*nnVec[2])[iZ] : (*nnVec[2])[0],
4950 nnVec[3]->size() > 1 ? (*nnVec[3])[iZ] : (*nnVec[3])[0]);
4951 degenVols.push_back( vol );
4959 return error("Not supported type of element", data._index);
4961 } // switch ( nbNodes )
4962 } // while ( fIt->more() )
4965 if ( !degenVols.empty() )
4967 SMESH_ComputeErrorPtr& err = _mesh->GetSubMesh( data._solid )->GetComputeError();
4968 if ( !err || err->IsOK() )
4970 err.reset( new SMESH_ComputeError( COMPERR_WARNING,
4971 "Degenerated volumes created" ));
4972 err->myBadElements.insert( err->myBadElements.end(),
4973 degenVols.begin(),degenVols.end() );
4980 //================================================================================
4982 * \brief Shrink 2D mesh on faces to let space for inflated layers
4984 //================================================================================
4986 bool _ViscousBuilder::shrink()
4988 // make map of (ids of FACEs to shrink mesh on) to (_SolidData containing _LayerEdge's
4989 // inflated along FACE or EDGE)
4990 map< TGeomID, _SolidData* > f2sdMap;
4991 for ( size_t i = 0 ; i < _sdVec.size(); ++i )
4993 _SolidData& data = _sdVec[i];
4994 TopTools_MapOfShape FFMap;
4995 map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
4996 for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
4997 if ( s2s->second.ShapeType() == TopAbs_FACE )
4999 f2sdMap.insert( make_pair( getMeshDS()->ShapeToIndex( s2s->second ), &data ));
5001 if ( FFMap.Add( (*s2s).second ))
5002 // Put mesh faces on the shrinked FACE to the proxy sub-mesh to avoid
5003 // usage of mesh faces made in addBoundaryElements() by the 3D algo or
5004 // by StdMeshers_QuadToTriaAdaptor
5005 if ( SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( s2s->second ))
5007 SMESH_ProxyMesh::SubMesh* proxySub =
5008 data._proxyMesh->getFaceSubM( TopoDS::Face( s2s->second ), /*create=*/true);
5009 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
5010 while ( fIt->more() )
5011 proxySub->AddElement( fIt->next() );
5012 // as a result 3D algo will use elements from proxySub and not from smDS
5017 SMESH_MesherHelper helper( *_mesh );
5018 helper.ToFixNodeParameters( true );
5021 map< TGeomID, _Shrinker1D > e2shrMap;
5022 vector< _LayerEdge* > lEdges;
5024 // loop on FACES to srink mesh on
5025 map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin();
5026 for ( ; f2sd != f2sdMap.end(); ++f2sd )
5028 _SolidData& data = *f2sd->second;
5029 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first ));
5030 SMESH_subMesh* sm = _mesh->GetSubMesh( F );
5031 SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
5033 Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
5035 helper.SetSubShape(F);
5037 // ===========================
5038 // Prepare data for shrinking
5039 // ===========================
5041 // Collect nodes to smooth, as src nodes are not yet replaced by tgt ones
5042 // and thus all nodes on a FACE connected to 2d elements are to be smoothed
5043 vector < const SMDS_MeshNode* > smoothNodes;
5045 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
5046 while ( nIt->more() )
5048 const SMDS_MeshNode* n = nIt->next();
5049 if ( n->NbInverseElements( SMDSAbs_Face ) > 0 )
5050 smoothNodes.push_back( n );
5053 // Find out face orientation
5055 const set<TGeomID> ignoreShapes;
5057 if ( !smoothNodes.empty() )
5059 vector<_Simplex> simplices;
5060 getSimplices( smoothNodes[0], simplices, ignoreShapes );
5061 helper.GetNodeUV( F, simplices[0]._nPrev, 0, &isOkUV ); // fix UV of silpmex nodes
5062 helper.GetNodeUV( F, simplices[0]._nNext, 0, &isOkUV );
5063 gp_XY uv = helper.GetNodeUV( F, smoothNodes[0], 0, &isOkUV );
5064 if ( !simplices[0].IsForward(uv, smoothNodes[0], F, helper,refSign) )
5068 // Find _LayerEdge's inflated along F
5071 set< TGeomID > subIDs;
5072 SMESH_subMeshIteratorPtr subIt = sm->getDependsOnIterator(/*includeSelf=*/false);
5073 while ( subIt->more() )
5074 subIDs.insert( subIt->next()->GetId() );
5077 for ( int iS = 0; iS < data._endEdgeOnShape.size() && !subIDs.empty(); ++iS )
5080 iEnd = data._endEdgeOnShape[ iS ];
5081 TGeomID shapeID = data._edges[ iBeg ]->_nodes[0]->getshapeId();
5082 set< TGeomID >::iterator idIt = subIDs.find( shapeID );
5083 if ( idIt == subIDs.end() ||
5084 data._edges[ iBeg ]->_sWOL.IsNull() ) continue;
5085 subIDs.erase( idIt );
5087 if ( !data._noShrinkShapes.count( shapeID ))
5088 for ( ; iBeg < iEnd; ++iBeg )
5090 lEdges.push_back( data._edges[ iBeg ] );
5091 prepareEdgeToShrink( *data._edges[ iBeg ], F, helper, smDS );
5096 dumpFunction(SMESH_Comment("beforeShrinkFace")<<f2sd->first); // debug
5097 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
5098 while ( fIt->more() )
5099 if ( const SMDS_MeshElement* f = fIt->next() )
5100 dumpChangeNodes( f );
5102 // Replace source nodes by target nodes in mesh faces to shrink
5103 dumpFunction(SMESH_Comment("replNodesOnFace")<<f2sd->first); // debug
5104 const SMDS_MeshNode* nodes[20];
5105 for ( size_t i = 0; i < lEdges.size(); ++i )
5107 _LayerEdge& edge = *lEdges[i];
5108 const SMDS_MeshNode* srcNode = edge._nodes[0];
5109 const SMDS_MeshNode* tgtNode = edge._nodes.back();
5110 SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
5111 while ( fIt->more() )
5113 const SMDS_MeshElement* f = fIt->next();
5114 if ( !smDS->Contains( f ))
5116 SMDS_NodeIteratorPtr nIt = f->nodeIterator();
5117 for ( int iN = 0; nIt->more(); ++iN )
5119 const SMDS_MeshNode* n = nIt->next();
5120 nodes[iN] = ( n == srcNode ? tgtNode : n );
5122 helper.GetMeshDS()->ChangeElementNodes( f, nodes, f->NbNodes() );
5123 dumpChangeNodes( f );
5127 // find out if a FACE is concave
5128 const bool isConcaveFace = isConcave( F, helper );
5130 // Create _SmoothNode's on face F
5131 vector< _SmoothNode > nodesToSmooth( smoothNodes.size() );
5133 dumpFunction(SMESH_Comment("fixUVOnFace")<<f2sd->first); // debug
5134 const bool sortSimplices = isConcaveFace;
5135 for ( size_t i = 0; i < smoothNodes.size(); ++i )
5137 const SMDS_MeshNode* n = smoothNodes[i];
5138 nodesToSmooth[ i ]._node = n;
5139 // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices
5140 getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, sortSimplices );
5141 // fix up incorrect uv of nodes on the FACE
5142 helper.GetNodeUV( F, n, 0, &isOkUV);
5146 //if ( nodesToSmooth.empty() ) continue;
5148 // Find EDGE's to shrink and set simpices to LayerEdge's
5149 set< _Shrinker1D* > eShri1D;
5151 for ( size_t i = 0; i < lEdges.size(); ++i )
5153 _LayerEdge* edge = lEdges[i];
5154 if ( edge->_sWOL.ShapeType() == TopAbs_EDGE )
5156 TGeomID edgeIndex = getMeshDS()->ShapeToIndex( edge->_sWOL );
5157 _Shrinker1D& srinker = e2shrMap[ edgeIndex ];
5158 eShri1D.insert( & srinker );
5159 srinker.AddEdge( edge, helper );
5160 VISCOUS_3D::ToClearSubWithMain( _mesh->GetSubMesh( edge->_sWOL ), data._solid );
5161 // restore params of nodes on EGDE if the EDGE has been already
5162 // srinked while srinking another FACE
5163 srinker.RestoreParams();
5165 getSimplices( /*tgtNode=*/edge->_nodes.back(), edge->_simplices, ignoreShapes );
5169 bool toFixTria = false; // to improve quality of trias by diagonal swap
5170 if ( isConcaveFace )
5172 const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles();
5173 if ( hasTria != hasQuad ) {
5174 toFixTria = hasTria;
5177 set<int> nbNodesSet;
5178 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
5179 while ( fIt->more() && nbNodesSet.size() < 2 )
5180 nbNodesSet.insert( fIt->next()->NbCornerNodes() );
5181 toFixTria = ( *nbNodesSet.begin() == 3 );
5185 // ==================
5186 // Perform shrinking
5187 // ==================
5189 bool shrinked = true;
5190 int badNb, shriStep=0, smooStep=0;
5191 _SmoothNode::SmoothType smoothType
5192 = isConcaveFace ? _SmoothNode::ANGULAR : _SmoothNode::LAPLACIAN;
5196 // Move boundary nodes (actually just set new UV)
5197 // -----------------------------------------------
5198 dumpFunction(SMESH_Comment("moveBoundaryOnF")<<f2sd->first<<"_st"<<shriStep ); // debug
5200 for ( size_t i = 0; i < lEdges.size(); ++i )
5202 shrinked |= lEdges[i]->SetNewLength2d( surface,F,helper );
5206 // Move nodes on EDGE's
5207 // (XYZ is set as soon as a needed length reached in SetNewLength2d())
5208 set< _Shrinker1D* >::iterator shr = eShri1D.begin();
5209 for ( ; shr != eShri1D.end(); ++shr )
5210 (*shr)->Compute( /*set3D=*/false, helper );
5213 // -----------------
5214 int nbNoImpSteps = 0;
5217 while (( nbNoImpSteps < 5 && badNb > 0) && moved)
5219 dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
5221 int oldBadNb = badNb;
5224 for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5226 moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
5227 smoothType, /*set3D=*/isConcaveFace);
5229 if ( badNb < oldBadNb )
5237 return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first );
5238 if ( shriStep > 200 )
5239 return error(SMESH_Comment("Infinite loop at shrinking 2D mesh on face ") << f2sd->first );
5241 // Fix narrow triangles by swapping diagonals
5242 // ---------------------------------------
5245 set<const SMDS_MeshNode*> usedNodes;
5246 fixBadFaces( F, helper, /*is2D=*/true, shriStep, & usedNodes); // swap diagonals
5248 // update working data
5249 set<const SMDS_MeshNode*>::iterator n;
5250 for ( size_t i = 0; i < nodesToSmooth.size() && !usedNodes.empty(); ++i )
5252 n = usedNodes.find( nodesToSmooth[ i ]._node );
5253 if ( n != usedNodes.end())
5255 getSimplices( nodesToSmooth[ i ]._node,
5256 nodesToSmooth[ i ]._simplices,
5258 /*sortSimplices=*/ smoothType == _SmoothNode::ANGULAR );
5259 usedNodes.erase( n );
5262 for ( size_t i = 0; i < lEdges.size() && !usedNodes.empty(); ++i )
5264 n = usedNodes.find( /*tgtNode=*/ lEdges[i]->_nodes.back() );
5265 if ( n != usedNodes.end())
5267 getSimplices( lEdges[i]->_nodes.back(),
5268 lEdges[i]->_simplices,
5270 usedNodes.erase( n );
5274 // TODO: check effect of this additional smooth
5275 // additional laplacian smooth to increase allowed shrink step
5276 // for ( int st = 1; st; --st )
5278 // dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
5279 // for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5281 // nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
5282 // _SmoothNode::LAPLACIAN,/*set3D=*/false);
5285 } // while ( shrinked )
5287 // No wrongly shaped faces remain; final smooth. Set node XYZ.
5288 bool isStructuredFixed = false;
5289 if ( SMESH_2D_Algo* algo = dynamic_cast<SMESH_2D_Algo*>( sm->GetAlgo() ))
5290 isStructuredFixed = algo->FixInternalNodes( *data._proxyMesh, F );
5291 if ( !isStructuredFixed )
5293 if ( isConcaveFace ) // fix narrow faces by swapping diagonals
5294 fixBadFaces( F, helper, /*is2D=*/false, ++shriStep );
5296 for ( int st = 3; st; --st )
5299 case 1: smoothType = _SmoothNode::LAPLACIAN; break;
5300 case 2: smoothType = _SmoothNode::LAPLACIAN; break;
5301 case 3: smoothType = _SmoothNode::ANGULAR; break;
5303 dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
5304 for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5306 nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
5307 smoothType,/*set3D=*/st==1 );
5312 // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh
5313 VISCOUS_3D::ToClearSubWithMain( sm, data._solid );
5315 if ( !getMeshDS()->IsEmbeddedMode() )
5316 // Log node movement
5317 for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5319 SMESH_TNodeXYZ p ( nodesToSmooth[i]._node );
5320 getMeshDS()->MoveNode( nodesToSmooth[i]._node, p.X(), p.Y(), p.Z() );
5323 } // loop on FACES to srink mesh on
5326 // Replace source nodes by target nodes in shrinked mesh edges
5328 map< int, _Shrinker1D >::iterator e2shr = e2shrMap.begin();
5329 for ( ; e2shr != e2shrMap.end(); ++e2shr )
5330 e2shr->second.SwapSrcTgtNodes( getMeshDS() );
5335 //================================================================================
5337 * \brief Computes 2d shrink direction and finds nodes limiting shrinking
5339 //================================================================================
5341 bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge,
5342 const TopoDS_Face& F,
5343 SMESH_MesherHelper& helper,
5344 const SMESHDS_SubMesh* faceSubMesh)
5346 const SMDS_MeshNode* srcNode = edge._nodes[0];
5347 const SMDS_MeshNode* tgtNode = edge._nodes.back();
5349 if ( edge._sWOL.ShapeType() == TopAbs_FACE )
5351 gp_XY srcUV( edge._pos[0].X(), edge._pos[0].Y() );//helper.GetNodeUV( F, srcNode );
5352 gp_XY tgtUV = edge.LastUV( F ); //helper.GetNodeUV( F, tgtNode );
5353 gp_Vec2d uvDir( srcUV, tgtUV );
5354 double uvLen = uvDir.Magnitude();
5356 edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0 );
5359 edge._pos.resize(1);
5360 edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 );
5362 // set UV of source node to target node
5363 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
5364 pos->SetUParameter( srcUV.X() );
5365 pos->SetVParameter( srcUV.Y() );
5367 else // _sWOL is TopAbs_EDGE
5369 const TopoDS_Edge& E = TopoDS::Edge( edge._sWOL );
5370 SMESHDS_SubMesh* edgeSM = getMeshDS()->MeshElements( E );
5371 if ( !edgeSM || edgeSM->NbElements() == 0 )
5372 return error(SMESH_Comment("Not meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
5374 const SMDS_MeshNode* n2 = 0;
5375 SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
5376 while ( eIt->more() && !n2 )
5378 const SMDS_MeshElement* e = eIt->next();
5379 if ( !edgeSM->Contains(e)) continue;
5380 n2 = e->GetNode( 0 );
5381 if ( n2 == srcNode ) n2 = e->GetNode( 1 );
5384 return error(SMESH_Comment("Wrongly meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
5386 double uSrc = helper.GetNodeU( E, srcNode, n2 );
5387 double uTgt = helper.GetNodeU( E, tgtNode, srcNode );
5388 double u2 = helper.GetNodeU( E, n2, srcNode );
5392 if ( fabs( uSrc-uTgt ) < 0.99 * fabs( uSrc-u2 ))
5394 // tgtNode is located so that it does not make faces with wrong orientation
5397 edge._pos.resize(1);
5398 edge._pos[0].SetCoord( U_TGT, uTgt );
5399 edge._pos[0].SetCoord( U_SRC, uSrc );
5400 edge._pos[0].SetCoord( LEN_TGT, fabs( uSrc-uTgt ));
5402 edge._simplices.resize( 1 );
5403 edge._simplices[0]._nPrev = n2;
5405 // set U of source node to the target node
5406 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
5407 pos->SetUParameter( uSrc );
5412 //================================================================================
5414 * \brief Restore position of a sole node of a _LayerEdge based on _noShrinkShapes
5416 //================================================================================
5418 void _ViscousBuilder::restoreNoShrink( _LayerEdge& edge ) const
5420 if ( edge._nodes.size() == 1 )
5425 const SMDS_MeshNode* srcNode = edge._nodes[0];
5426 TopoDS_Shape S = SMESH_MesherHelper::GetSubShapeByNode( srcNode, getMeshDS() );
5427 if ( S.IsNull() ) return;
5431 switch ( S.ShapeType() )
5436 TopLoc_Location loc;
5437 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( S ), loc, f, l );
5438 if ( curve.IsNull() ) return;
5439 SMDS_EdgePosition* ePos = static_cast<SMDS_EdgePosition*>( srcNode->GetPosition() );
5440 p = curve->Value( ePos->GetUParameter() );
5445 p = BRep_Tool::Pnt( TopoDS::Vertex( S ));
5450 getMeshDS()->MoveNode( srcNode, p.X(), p.Y(), p.Z() );
5451 dumpMove( srcNode );
5455 //================================================================================
5457 * \brief Try to fix triangles with high aspect ratio by swaping diagonals
5459 //================================================================================
5461 void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F,
5462 SMESH_MesherHelper& helper,
5465 set<const SMDS_MeshNode*> * involvedNodes)
5467 SMESH::Controls::AspectRatio qualifier;
5468 SMESH::Controls::TSequenceOfXYZ points(3), points1(3), points2(3);
5469 const double maxAspectRatio = is2D ? 4. : 2;
5470 _NodeCoordHelper xyz( F, helper, is2D );
5472 // find bad triangles
5474 vector< const SMDS_MeshElement* > badTrias;
5475 vector< double > badAspects;
5476 SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( F );
5477 SMDS_ElemIteratorPtr fIt = sm->GetElements();
5478 while ( fIt->more() )
5480 const SMDS_MeshElement * f = fIt->next();
5481 if ( f->NbCornerNodes() != 3 ) continue;
5482 for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = xyz( f->GetNode(iP));
5483 double aspect = qualifier.GetValue( points );
5484 if ( aspect > maxAspectRatio )
5486 badTrias.push_back( f );
5487 badAspects.push_back( aspect );
5492 dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID());
5493 SMDS_ElemIteratorPtr fIt = sm->GetElements();
5494 while ( fIt->more() )
5496 const SMDS_MeshElement * f = fIt->next();
5497 if ( f->NbCornerNodes() == 3 )
5498 dumpChangeNodes( f );
5502 if ( badTrias.empty() )
5505 // find couples of faces to swap diagonal
5507 typedef pair < const SMDS_MeshElement* , const SMDS_MeshElement* > T2Trias;
5508 vector< T2Trias > triaCouples;
5510 TIDSortedElemSet involvedFaces, emptySet;
5511 for ( size_t iTia = 0; iTia < badTrias.size(); ++iTia )
5514 double aspRatio [3];
5517 if ( !involvedFaces.insert( badTrias[iTia] ).second )
5519 for ( int iP = 0; iP < 3; ++iP )
5520 points(iP+1) = xyz( badTrias[iTia]->GetNode(iP));
5522 // find triangles adjacent to badTrias[iTia] with better aspect ratio after diag-swaping
5523 int bestCouple = -1;
5524 for ( int iSide = 0; iSide < 3; ++iSide )
5526 const SMDS_MeshNode* n1 = badTrias[iTia]->GetNode( iSide );
5527 const SMDS_MeshNode* n2 = badTrias[iTia]->GetNode(( iSide+1 ) % 3 );
5528 trias [iSide].first = badTrias[iTia];
5529 trias [iSide].second = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, involvedFaces,
5531 if (( ! trias[iSide].second ) ||
5532 ( trias[iSide].second->NbCornerNodes() != 3 ) ||
5533 ( ! sm->Contains( trias[iSide].second )))
5536 // aspect ratio of an adjacent tria
5537 for ( int iP = 0; iP < 3; ++iP )
5538 points2(iP+1) = xyz( trias[iSide].second->GetNode(iP));
5539 double aspectInit = qualifier.GetValue( points2 );
5541 // arrange nodes as after diag-swaping
5542 if ( helper.WrapIndex( i1+1, 3 ) == i2 )
5543 i3 = helper.WrapIndex( i1-1, 3 );
5545 i3 = helper.WrapIndex( i1+1, 3 );
5547 points1( 1+ iSide ) = points2( 1+ i3 );
5548 points2( 1+ i2 ) = points1( 1+ ( iSide+2 ) % 3 );
5550 // aspect ratio after diag-swaping
5551 aspRatio[ iSide ] = qualifier.GetValue( points1 ) + qualifier.GetValue( points2 );
5552 if ( aspRatio[ iSide ] > aspectInit + badAspects[ iTia ] )
5555 // prevent inversion of a triangle
5556 gp_Vec norm1 = gp_Vec( points1(1), points1(3) ) ^ gp_Vec( points1(1), points1(2) );
5557 gp_Vec norm2 = gp_Vec( points2(1), points2(3) ) ^ gp_Vec( points2(1), points2(2) );
5558 if ( norm1 * norm2 < 0. && norm1.Angle( norm2 ) > 70./180.*M_PI )
5561 if ( bestCouple < 0 || aspRatio[ bestCouple ] > aspRatio[ iSide ] )
5565 if ( bestCouple >= 0 )
5567 triaCouples.push_back( trias[bestCouple] );
5568 involvedFaces.insert ( trias[bestCouple].second );
5572 involvedFaces.erase( badTrias[iTia] );
5575 if ( triaCouples.empty() )
5580 SMESH_MeshEditor editor( helper.GetMesh() );
5581 dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
5582 for ( size_t i = 0; i < triaCouples.size(); ++i )
5584 dumpChangeNodes( triaCouples[i].first );
5585 dumpChangeNodes( triaCouples[i].second );
5586 editor.InverseDiag( triaCouples[i].first, triaCouples[i].second );
5589 if ( involvedNodes )
5590 for ( size_t i = 0; i < triaCouples.size(); ++i )
5592 involvedNodes->insert( triaCouples[i].first->begin_nodes(),
5593 triaCouples[i].first->end_nodes() );
5594 involvedNodes->insert( triaCouples[i].second->begin_nodes(),
5595 triaCouples[i].second->end_nodes() );
5598 // just for debug dump resulting triangles
5599 dumpFunction(SMESH_Comment("swapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
5600 for ( size_t i = 0; i < triaCouples.size(); ++i )
5602 dumpChangeNodes( triaCouples[i].first );
5603 dumpChangeNodes( triaCouples[i].second );
5607 //================================================================================
5609 * \brief Move target node to it's final position on the FACE during shrinking
5611 //================================================================================
5613 bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
5614 const TopoDS_Face& F,
5615 SMESH_MesherHelper& helper )
5618 return false; // already at the target position
5620 SMDS_MeshNode* tgtNode = const_cast< SMDS_MeshNode*& >( _nodes.back() );
5622 if ( _sWOL.ShapeType() == TopAbs_FACE )
5624 gp_XY curUV = helper.GetNodeUV( F, tgtNode );
5625 gp_Pnt2d tgtUV( _pos[0].X(), _pos[0].Y() );
5626 gp_Vec2d uvDir( _normal.X(), _normal.Y() );
5627 const double uvLen = tgtUV.Distance( curUV );
5628 const double kSafe = Max( 0.5, 1. - 0.1 * _simplices.size() );
5630 // Select shrinking step such that not to make faces with wrong orientation.
5631 double stepSize = 1e100;
5632 for ( size_t i = 0; i < _simplices.size(); ++i )
5634 // find intersection of 2 lines: curUV-tgtUV and that connecting simplex nodes
5635 gp_XY uvN1 = helper.GetNodeUV( F, _simplices[i]._nPrev );
5636 gp_XY uvN2 = helper.GetNodeUV( F, _simplices[i]._nNext );
5637 gp_XY dirN = uvN2 - uvN1;
5638 double det = uvDir.Crossed( dirN );
5639 if ( Abs( det ) < std::numeric_limits<double>::min() ) continue;
5640 gp_XY dirN2Cur = curUV - uvN1;
5641 double step = dirN.Crossed( dirN2Cur ) / det;
5643 stepSize = Min( step, stepSize );
5646 if ( uvLen <= stepSize )
5651 else if ( stepSize > 0 )
5653 newUV = curUV + uvDir.XY() * stepSize * kSafe;
5659 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
5660 pos->SetUParameter( newUV.X() );
5661 pos->SetVParameter( newUV.Y() );
5664 gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
5665 tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
5666 dumpMove( tgtNode );
5669 else // _sWOL is TopAbs_EDGE
5671 const TopoDS_Edge& E = TopoDS::Edge( _sWOL );
5672 const SMDS_MeshNode* n2 = _simplices[0]._nPrev;
5673 SMDS_EdgePosition* tgtPos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
5675 const double u2 = helper.GetNodeU( E, n2, tgtNode );
5676 const double uSrc = _pos[0].Coord( U_SRC );
5677 const double lenTgt = _pos[0].Coord( LEN_TGT );
5679 double newU = _pos[0].Coord( U_TGT );
5680 if ( lenTgt < 0.99 * fabs( uSrc-u2 )) // n2 got out of src-tgt range
5686 newU = 0.1 * tgtPos->GetUParameter() + 0.9 * u2;
5688 tgtPos->SetUParameter( newU );
5690 gp_XY newUV = helper.GetNodeUV( F, tgtNode, _nodes[0]);
5691 gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
5692 tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
5693 dumpMove( tgtNode );
5699 //================================================================================
5701 * \brief Perform smooth on the FACE
5702 * \retval bool - true if the node has been moved
5704 //================================================================================
5706 bool _SmoothNode::Smooth(int& badNb,
5707 Handle(Geom_Surface)& surface,
5708 SMESH_MesherHelper& helper,
5709 const double refSign,
5713 const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() );
5715 // get uv of surrounding nodes
5716 vector<gp_XY> uv( _simplices.size() );
5717 for ( size_t i = 0; i < _simplices.size(); ++i )
5718 uv[i] = helper.GetNodeUV( face, _simplices[i]._nPrev, _node );
5720 // compute new UV for the node
5722 if ( how == TFI && _simplices.size() == 4 )
5725 for ( size_t i = 0; i < _simplices.size(); ++i )
5726 if ( _simplices[i]._nOpp )
5727 corners[i] = helper.GetNodeUV( face, _simplices[i]._nOpp, _node );
5729 throw SALOME_Exception(LOCALIZED("TFI smoothing: _Simplex::_nOpp not set!"));
5731 newPos = helper.calcTFI ( 0.5, 0.5,
5732 corners[0], corners[1], corners[2], corners[3],
5733 uv[1], uv[2], uv[3], uv[0] );
5735 else if ( how == ANGULAR )
5737 newPos = computeAngularPos( uv, helper.GetNodeUV( face, _node ), refSign );
5739 else if ( how == CENTROIDAL && _simplices.size() > 3 )
5741 // average centers of diagonals wieghted with their reciprocal lengths
5742 if ( _simplices.size() == 4 )
5744 double w1 = 1. / ( uv[2]-uv[0] ).SquareModulus();
5745 double w2 = 1. / ( uv[3]-uv[1] ).SquareModulus();
5746 newPos = ( w1 * ( uv[2]+uv[0] ) + w2 * ( uv[3]+uv[1] )) / ( w1+w2 ) / 2;
5750 double sumWeight = 0;
5751 int nb = _simplices.size() == 4 ? 2 : _simplices.size();
5752 for ( int i = 0; i < nb; ++i )
5755 int iTo = i + _simplices.size() - 1;
5756 for ( int j = iFrom; j < iTo; ++j )
5758 int i2 = SMESH_MesherHelper::WrapIndex( j, _simplices.size() );
5759 double w = 1. / ( uv[i]-uv[i2] ).SquareModulus();
5761 newPos += w * ( uv[i]+uv[i2] );
5764 newPos /= 2 * sumWeight; // 2 is to get a middle between uv's
5770 for ( size_t i = 0; i < _simplices.size(); ++i )
5772 newPos /= _simplices.size();
5775 // count quality metrics (orientation) of triangles around the node
5777 gp_XY tgtUV = helper.GetNodeUV( face, _node );
5778 for ( size_t i = 0; i < _simplices.size(); ++i )
5779 nbOkBefore += _simplices[i].IsForward( tgtUV, _node, face, helper, refSign );
5782 for ( size_t i = 0; i < _simplices.size(); ++i )
5783 nbOkAfter += _simplices[i].IsForward( newPos, _node, face, helper, refSign );
5785 if ( nbOkAfter < nbOkBefore )
5787 badNb += _simplices.size() - nbOkBefore;
5791 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( _node->GetPosition() );
5792 pos->SetUParameter( newPos.X() );
5793 pos->SetVParameter( newPos.Y() );
5800 gp_Pnt p = surface->Value( newPos.X(), newPos.Y() );
5801 const_cast< SMDS_MeshNode* >( _node )->setXYZ( p.X(), p.Y(), p.Z() );
5805 badNb += _simplices.size() - nbOkAfter;
5806 return ( (tgtUV-newPos).SquareModulus() > 1e-10 );
5809 //================================================================================
5811 * \brief Computes new UV using angle based smoothing technic
5813 //================================================================================
5815 gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
5816 const gp_XY& uvToFix,
5817 const double refSign)
5819 uv.push_back( uv.front() );
5821 vector< gp_XY > edgeDir ( uv.size() );
5822 vector< double > edgeSize( uv.size() );
5823 for ( size_t i = 1; i < edgeDir.size(); ++i )
5825 edgeDir [i-1] = uv[i] - uv[i-1];
5826 edgeSize[i-1] = edgeDir[i-1].Modulus();
5827 if ( edgeSize[i-1] < numeric_limits<double>::min() )
5828 edgeDir[i-1].SetX( 100 );
5830 edgeDir[i-1] /= edgeSize[i-1] * refSign;
5832 edgeDir.back() = edgeDir.front();
5833 edgeSize.back() = edgeSize.front();
5838 for ( size_t i = 1; i < edgeDir.size(); ++i )
5840 if ( edgeDir[i-1].X() > 1. ) continue;
5842 while ( edgeDir[i].X() > 1. && ++i < edgeDir.size() );
5843 if ( i == edgeDir.size() ) break;
5845 gp_XY norm1( -edgeDir[i1].Y(), edgeDir[i1].X() );
5846 gp_XY norm2( -edgeDir[i].Y(), edgeDir[i].X() );
5847 gp_XY bisec = norm1 + norm2;
5848 double bisecSize = bisec.Modulus();
5849 if ( bisecSize < numeric_limits<double>::min() )
5851 bisec = -edgeDir[i1] + edgeDir[i];
5852 bisecSize = bisec.Modulus();
5856 gp_XY dirToN = uvToFix - p;
5857 double distToN = dirToN.Modulus();
5858 if ( bisec * dirToN < 0 )
5861 newPos += ( p + bisec * distToN ) * ( edgeSize[i1] + edgeSize[i] );
5863 sumSize += edgeSize[i1] + edgeSize[i];
5865 newPos /= /*nbEdges * */sumSize;
5869 //================================================================================
5871 * \brief Delete _SolidData
5873 //================================================================================
5875 _SolidData::~_SolidData()
5877 for ( size_t i = 0; i < _edges.size(); ++i )
5879 if ( _edges[i] && _edges[i]->_2neibors )
5880 delete _edges[i]->_2neibors;
5885 //================================================================================
5887 * \brief Add a _LayerEdge inflated along the EDGE
5889 //================================================================================
5891 void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper )
5894 if ( _nodes.empty() )
5896 _edges[0] = _edges[1] = 0;
5900 if ( e == _edges[0] || e == _edges[1] )
5902 if ( e->_sWOL.IsNull() || e->_sWOL.ShapeType() != TopAbs_EDGE )
5903 throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
5904 if ( _edges[0] && _edges[0]->_sWOL != e->_sWOL )
5905 throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
5908 const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
5910 BRep_Tool::Range( E, f,l );
5911 double u = helper.GetNodeU( E, e->_nodes[0], e->_nodes.back());
5912 _edges[ u < 0.5*(f+l) ? 0 : 1 ] = e;
5916 const SMDS_MeshNode* tgtNode0 = _edges[0] ? _edges[0]->_nodes.back() : 0;
5917 const SMDS_MeshNode* tgtNode1 = _edges[1] ? _edges[1]->_nodes.back() : 0;
5919 if ( _nodes.empty() )
5921 SMESHDS_SubMesh * eSubMesh = helper.GetMeshDS()->MeshElements( E );
5922 if ( !eSubMesh || eSubMesh->NbNodes() < 1 )
5924 TopLoc_Location loc;
5925 Handle(Geom_Curve) C = BRep_Tool::Curve(E, loc, f,l);
5926 GeomAdaptor_Curve aCurve(C, f,l);
5927 const double totLen = GCPnts_AbscissaPoint::Length(aCurve, f, l);
5929 int nbExpectNodes = eSubMesh->NbNodes();
5930 _initU .reserve( nbExpectNodes );
5931 _normPar.reserve( nbExpectNodes );
5932 _nodes .reserve( nbExpectNodes );
5933 SMDS_NodeIteratorPtr nIt = eSubMesh->GetNodes();
5934 while ( nIt->more() )
5936 const SMDS_MeshNode* node = nIt->next();
5937 if ( node->NbInverseElements(SMDSAbs_Edge) == 0 ||
5938 node == tgtNode0 || node == tgtNode1 )
5939 continue; // refinement nodes
5940 _nodes.push_back( node );
5941 _initU.push_back( helper.GetNodeU( E, node ));
5942 double len = GCPnts_AbscissaPoint::Length(aCurve, f, _initU.back());
5943 _normPar.push_back( len / totLen );
5948 // remove target node of the _LayerEdge from _nodes
5950 for ( size_t i = 0; i < _nodes.size(); ++i )
5951 if ( !_nodes[i] || _nodes[i] == tgtNode0 || _nodes[i] == tgtNode1 )
5952 _nodes[i] = 0, nbFound++;
5953 if ( nbFound == _nodes.size() )
5958 //================================================================================
5960 * \brief Move nodes on EDGE from ends where _LayerEdge's are inflated
5962 //================================================================================
5964 void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
5966 if ( _done || _nodes.empty())
5968 const _LayerEdge* e = _edges[0];
5969 if ( !e ) e = _edges[1];
5972 _done = (( !_edges[0] || _edges[0]->_pos.empty() ) &&
5973 ( !_edges[1] || _edges[1]->_pos.empty() ));
5975 const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
5977 if ( set3D || _done )
5979 Handle(Geom_Curve) C = BRep_Tool::Curve(E, f,l);
5980 GeomAdaptor_Curve aCurve(C, f,l);
5983 f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
5985 l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
5986 double totLen = GCPnts_AbscissaPoint::Length( aCurve, f, l );
5988 for ( size_t i = 0; i < _nodes.size(); ++i )
5990 if ( !_nodes[i] ) continue;
5991 double len = totLen * _normPar[i];
5992 GCPnts_AbscissaPoint discret( aCurve, len, f );
5993 if ( !discret.IsDone() )
5994 return throw SALOME_Exception(LOCALIZED("GCPnts_AbscissaPoint failed"));
5995 double u = discret.Parameter();
5996 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
5997 pos->SetUParameter( u );
5998 gp_Pnt p = C->Value( u );
5999 const_cast< SMDS_MeshNode*>( _nodes[i] )->setXYZ( p.X(), p.Y(), p.Z() );
6004 BRep_Tool::Range( E, f,l );
6006 f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
6008 l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
6010 for ( size_t i = 0; i < _nodes.size(); ++i )
6012 if ( !_nodes[i] ) continue;
6013 double u = f * ( 1-_normPar[i] ) + l * _normPar[i];
6014 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
6015 pos->SetUParameter( u );
6020 //================================================================================
6022 * \brief Restore initial parameters of nodes on EDGE
6024 //================================================================================
6026 void _Shrinker1D::RestoreParams()
6029 for ( size_t i = 0; i < _nodes.size(); ++i )
6031 if ( !_nodes[i] ) continue;
6032 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
6033 pos->SetUParameter( _initU[i] );
6038 //================================================================================
6040 * \brief Replace source nodes by target nodes in shrinked mesh edges
6042 //================================================================================
6044 void _Shrinker1D::SwapSrcTgtNodes( SMESHDS_Mesh* mesh )
6046 const SMDS_MeshNode* nodes[3];
6047 for ( int i = 0; i < 2; ++i )
6049 if ( !_edges[i] ) continue;
6051 SMESHDS_SubMesh * eSubMesh = mesh->MeshElements( _edges[i]->_sWOL );
6052 if ( !eSubMesh ) return;
6053 const SMDS_MeshNode* srcNode = _edges[i]->_nodes[0];
6054 const SMDS_MeshNode* tgtNode = _edges[i]->_nodes.back();
6055 SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
6056 while ( eIt->more() )
6058 const SMDS_MeshElement* e = eIt->next();
6059 if ( !eSubMesh->Contains( e ))
6061 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
6062 for ( int iN = 0; iN < e->NbNodes(); ++iN )
6064 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
6065 nodes[iN] = ( n == srcNode ? tgtNode : n );
6067 mesh->ChangeElementNodes( e, nodes, e->NbNodes() );
6072 //================================================================================
6074 * \brief Creates 2D and 1D elements on boundaries of new prisms
6076 //================================================================================
6078 bool _ViscousBuilder::addBoundaryElements()
6080 SMESH_MesherHelper helper( *_mesh );
6082 for ( size_t i = 0; i < _sdVec.size(); ++i )
6084 _SolidData& data = _sdVec[i];
6085 TopTools_IndexedMapOfShape geomEdges;
6086 TopExp::MapShapes( data._solid, TopAbs_EDGE, geomEdges );
6087 for ( int iE = 1; iE <= geomEdges.Extent(); ++iE )
6089 const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE));
6090 if ( data._noShrinkShapes.count( getMeshDS()->ShapeToIndex( E )))
6093 // Get _LayerEdge's based on E
6095 map< double, const SMDS_MeshNode* > u2nodes;
6096 if ( !SMESH_Algo::GetSortedNodesOnEdge( getMeshDS(), E, /*ignoreMedium=*/false, u2nodes))
6099 vector< _LayerEdge* > ledges; ledges.reserve( u2nodes.size() );
6100 TNode2Edge & n2eMap = data._n2eMap;
6101 map< double, const SMDS_MeshNode* >::iterator u2n = u2nodes.begin();
6103 //check if 2D elements are needed on E
6104 TNode2Edge::iterator n2e = n2eMap.find( u2n->second );
6105 if ( n2e == n2eMap.end() ) continue; // no layers on vertex
6106 ledges.push_back( n2e->second );
6108 if (( n2e = n2eMap.find( u2n->second )) == n2eMap.end() )
6109 continue; // no layers on E
6110 ledges.push_back( n2eMap[ u2n->second ]);
6112 const SMDS_MeshNode* tgtN0 = ledges[0]->_nodes.back();
6113 const SMDS_MeshNode* tgtN1 = ledges[1]->_nodes.back();
6114 int nbSharedPyram = 0;
6115 SMDS_ElemIteratorPtr vIt = tgtN0->GetInverseElementIterator(SMDSAbs_Volume);
6116 while ( vIt->more() )
6118 const SMDS_MeshElement* v = vIt->next();
6119 nbSharedPyram += int( v->GetNodeIndex( tgtN1 ) >= 0 );
6121 if ( nbSharedPyram > 1 )
6122 continue; // not free border of the pyramid
6124 if ( getMeshDS()->FindFace( ledges[0]->_nodes[0], ledges[0]->_nodes[1],
6125 ledges[1]->_nodes[0], ledges[1]->_nodes[1]))
6126 continue; // faces already created
6128 for ( ++u2n; u2n != u2nodes.end(); ++u2n )
6129 ledges.push_back( n2eMap[ u2n->second ]);
6131 // Find out orientation and type of face to create
6133 bool reverse = false, isOnFace;
6135 map< TGeomID, TopoDS_Shape >::iterator e2f =
6136 data._shrinkShape2Shape.find( getMeshDS()->ShapeToIndex( E ));
6138 if (( isOnFace = ( e2f != data._shrinkShape2Shape.end() )))
6140 F = e2f->second.Oriented( TopAbs_FORWARD );
6141 reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED );
6142 if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED )
6143 reverse = !reverse, F.Reverse();
6144 if ( helper.IsReversedSubMesh( TopoDS::Face(F) ))
6149 // find FACE with layers sharing E
6150 PShapeIteratorPtr fIt = helper.GetAncestors( E, *_mesh, TopAbs_FACE );
6151 while ( fIt->more() && F.IsNull() )
6153 const TopoDS_Shape* pF = fIt->next();
6154 if ( helper.IsSubShape( *pF, data._solid) &&
6155 !data._ignoreFaceIds.count( e2f->first ))
6159 // Find the sub-mesh to add new faces
6160 SMESHDS_SubMesh* sm = 0;
6162 sm = getMeshDS()->MeshElements( F );
6164 sm = data._proxyMesh->getFaceSubM( TopoDS::Face(F), /*create=*/true );
6166 return error("error in addBoundaryElements()", data._index);
6169 const int dj1 = reverse ? 0 : 1;
6170 const int dj2 = reverse ? 1 : 0;
6171 for ( size_t j = 1; j < ledges.size(); ++j )
6173 vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes;
6174 vector< const SMDS_MeshNode*>& nn2 = ledges[j-dj2]->_nodes;
6175 if ( nn1.size() == nn2.size() )
6178 for ( size_t z = 1; z < nn1.size(); ++z )
6179 sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
6181 for ( size_t z = 1; z < nn1.size(); ++z )
6182 sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
6184 else if ( nn1.size() == 1 )
6187 for ( size_t z = 1; z < nn2.size(); ++z )
6188 sm->AddElement( getMeshDS()->AddFace( nn1[0], nn2[z-1], nn2[z] ));
6190 for ( size_t z = 1; z < nn2.size(); ++z )
6191 sm->AddElement( new SMDS_FaceOfNodes( nn1[0], nn2[z-1], nn2[z] ));
6196 for ( size_t z = 1; z < nn1.size(); ++z )
6197 sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[0], nn1[z] ));
6199 for ( size_t z = 1; z < nn1.size(); ++z )
6200 sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[0], nn2[z] ));
6205 for ( int isFirst = 0; isFirst < 2; ++isFirst )
6207 _LayerEdge* edge = isFirst ? ledges.front() : ledges.back();
6208 if ( !edge->_sWOL.IsNull() && edge->_sWOL.ShapeType() == TopAbs_EDGE )
6210 vector< const SMDS_MeshNode*>& nn = edge->_nodes;
6211 if ( nn.size() < 2 || nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() )
6213 helper.SetSubShape( edge->_sWOL );
6214 helper.SetElementsOnShape( true );
6215 for ( size_t z = 1; z < nn.size(); ++z )
6216 helper.AddEdge( nn[z-1], nn[z] );