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 * Structure used to smooth a _LayerEdge (master) based on an EDGE.
305 // target nodes of 2 neighbour _LayerEdge's based on the same EDGE
306 const SMDS_MeshNode* _nodes[2];
307 // vectors from source nodes of 2 _LayerEdge's to the source node of master _LayerEdge
309 double _wgt[2]; // weights of _nodes
310 _LayerEdge* _edges[2];
312 // normal to plane passing through _LayerEdge._normal and tangent of EDGE
315 _2NearEdges() { _nodes[0]=_nodes[1]=0; _plnNorm = 0; }
317 std::swap( _nodes[0], _nodes[1] );
318 std::swap( _wgt [0], _wgt [1] );
319 std::swap( _edges[0], _edges[1] );
322 //--------------------------------------------------------------------------------
324 * \brief Edge normal to surface, connecting a node on solid surface (_nodes[0])
325 * and a node of the most internal layer (_nodes.back())
329 vector< const SMDS_MeshNode*> _nodes;
331 gp_XYZ _normal; // to solid surface
332 vector<gp_XYZ> _pos; // points computed during inflation
333 double _len; // length achived with the last inflation step
334 double _cosin; // of angle (_normal ^ surface)
335 double _lenFactor; // to compute _len taking _cosin into account
337 // face or edge w/o layer along or near which _LayerEdge is inflated
339 // simplices connected to the source node (_nodes[0]);
340 // used for smoothing and quality check of _LayerEdge's based on the FACE
341 vector<_Simplex> _simplices;
342 // data for smoothing of _LayerEdge's based on the EDGE
343 _2NearEdges* _2neibors;
345 _Curvature* _curvature;
346 // TODO:: detele _Curvature, _plnNorm
348 void SetNewLength( double len, SMESH_MesherHelper& helper );
349 bool SetNewLength2d( Handle(Geom_Surface)& surface,
350 const TopoDS_Face& F,
351 SMESH_MesherHelper& helper );
352 void SetDataByNeighbors( const SMDS_MeshNode* n1,
353 const SMDS_MeshNode* n2,
354 SMESH_MesherHelper& helper);
355 void InvalidateStep( int curStep, bool restoreLength=false );
356 bool Smooth(int& badNb);
357 bool SmoothOnEdge(Handle(Geom_Surface)& surface,
358 const TopoDS_Face& F,
359 SMESH_MesherHelper& helper);
360 bool FindIntersection( SMESH_ElementSearcher& searcher,
362 const double& epsilon,
363 const SMDS_MeshElement** face = 0);
364 bool SegTriaInter( const gp_Ax1& lastSegment,
365 const SMDS_MeshNode* n0,
366 const SMDS_MeshNode* n1,
367 const SMDS_MeshNode* n2,
369 const double& epsilon) const;
370 gp_Ax1 LastSegment(double& segLen) const;
371 bool IsOnEdge() const { return _2neibors; }
372 gp_XYZ Copy( _LayerEdge& other, SMESH_MesherHelper& helper );
373 void SetCosin( double cosin );
377 bool operator () (const _LayerEdge* e1, const _LayerEdge* e2) const
379 const bool cmpNodes = ( e1 && e2 && e1->_nodes.size() && e2->_nodes.size() );
380 return cmpNodes ? ( e1->_nodes[0]->GetID() < e2->_nodes[0]->GetID()) : ( e1 < e2 );
383 //--------------------------------------------------------------------------------
385 * \brief Convex FACE whose radius of curvature is less than the thickness of
386 * layers. It is used to detect distortion of prisms based on a convex
387 * FACE and to update normals to enable further increasing the thickness
393 // edges whose _simplices are used to detect prism destorsion
394 vector< _LayerEdge* > _simplexTestEdges;
396 // map a sub-shape to it's index in _SolidData::_endEdgeOnShape vector
397 map< TGeomID, int > _subIdToEdgeEnd;
401 bool GetCenterOfCurvature( _LayerEdge* ledge,
402 BRepLProp_SLProps& surfProp,
403 SMESH_MesherHelper& helper,
404 gp_Pnt & center ) const;
405 bool CheckPrisms() const;
408 //--------------------------------------------------------------------------------
410 typedef map< const SMDS_MeshNode*, _LayerEdge*, TIDCompare > TNode2Edge;
412 //--------------------------------------------------------------------------------
414 * \brief Data of a SOLID
419 const StdMeshers_ViscousLayers* _hyp;
420 TopoDS_Shape _hypShape;
421 _MeshOfSolid* _proxyMesh;
422 set<TGeomID> _reversedFaceIds;
423 set<TGeomID> _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDS
425 double _stepSize, _stepSizeCoeff;
426 const SMDS_MeshNode* _stepSizeNodes[2];
429 // map to find _n2eMap of another _SolidData by a shrink shape shared by two _SolidData's
430 map< TGeomID, TNode2Edge* > _s2neMap;
431 // edges of _n2eMap. We keep same data in two containers because
432 // iteration over the map is 5 time longer than over the vector
433 vector< _LayerEdge* > _edges;
435 // key: an id of shape (EDGE or VERTEX) shared by a FACE with
436 // layers and a FACE w/o layers
437 // value: the shape (FACE or EDGE) to shrink mesh on.
438 // _LayerEdge's basing on nodes on key shape are inflated along the value shape
439 map< TGeomID, TopoDS_Shape > _shrinkShape2Shape;
441 // Convex FACEs whose radius of curvature is less than the thickness of layers
442 map< TGeomID, _ConvexFace > _convexFaces;
444 // FACE's WOL, srink on which is forbiden due to algo on the adjacent SOLID
445 set< TGeomID > _noShrinkFaces;
447 // <EDGE to smooth on> to <it's curve> -- for analytic smooth
448 map< TGeomID,Handle(Geom_Curve)> _edge2curve;
450 // end indices in _edges of _LayerEdge on each shape, first go shapes to smooth
451 vector< int > _endEdgeOnShape;
452 int _nbShapesToSmooth;
454 double _epsilon; // precision for SegTriaInter()
456 int _index; // for debug
458 _SolidData(const TopoDS_Shape& s=TopoDS_Shape(),
459 const StdMeshers_ViscousLayers* h=0,
460 const TopoDS_Shape& hs=TopoDS_Shape(),
462 :_solid(s), _hyp(h), _hypShape(hs), _proxyMesh(m) {}
465 Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge& E,
468 Handle(Geom_Surface)& surface,
469 const TopoDS_Face& F,
470 SMESH_MesherHelper& helper);
472 void SortOnEdge( const TopoDS_Edge& E,
475 SMESH_MesherHelper& helper);
477 _ConvexFace* GetConvexFace( const TGeomID faceID )
479 map< TGeomID, _ConvexFace >::iterator id2face = _convexFaces.find( faceID );
480 return id2face == _convexFaces.end() ? 0 : & id2face->second;
482 void GetEdgesOnShape( size_t end, int & iBeg, int & iEnd )
484 iBeg = end > 0 ? _endEdgeOnShape[ end-1 ] : 0;
485 iEnd = _endEdgeOnShape[ end ];
488 bool GetShapeEdges(const TGeomID shapeID, size_t& edgeEnd, int* iBeg=0, int* iEnd=0 ) const;
490 void AddShapesToSmooth( const set< TGeomID >& shapeIDs );
492 //--------------------------------------------------------------------------------
494 * \brief Container of centers of curvature at nodes on an EDGE bounding _ConvexFace
496 struct _CentralCurveOnEdge
499 vector< gp_Pnt > _curvaCenters;
500 vector< _LayerEdge* > _ledges;
501 vector< gp_XYZ > _normals; // new normal for each of _ledges
502 vector< double > _segLength2;
505 TopoDS_Face _adjFace;
506 bool _adjFaceToSmooth;
508 void Append( const gp_Pnt& center, _LayerEdge* ledge )
510 if ( _curvaCenters.size() > 0 )
511 _segLength2.push_back( center.SquareDistance( _curvaCenters.back() ));
512 _curvaCenters.push_back( center );
513 _ledges.push_back( ledge );
514 _normals.push_back( ledge->_normal );
516 bool FindNewNormal( const gp_Pnt& center, gp_XYZ& newNormal );
517 void SetShapes( const TopoDS_Edge& edge,
518 const _ConvexFace& convFace,
519 const _SolidData& data,
520 SMESH_MesherHelper& helper);
522 //--------------------------------------------------------------------------------
524 * \brief Data of node on a shrinked FACE
528 const SMDS_MeshNode* _node;
529 vector<_Simplex> _simplices; // for quality check
531 enum SmoothType { LAPLACIAN, CENTROIDAL, ANGULAR, TFI };
533 bool Smooth(int& badNb,
534 Handle(Geom_Surface)& surface,
535 SMESH_MesherHelper& helper,
536 const double refSign,
540 gp_XY computeAngularPos(vector<gp_XY>& uv,
541 const gp_XY& uvToFix,
542 const double refSign );
544 //--------------------------------------------------------------------------------
546 * \brief Builder of viscous layers
548 class _ViscousBuilder
553 SMESH_ComputeErrorPtr Compute(SMESH_Mesh& mesh,
554 const TopoDS_Shape& shape);
556 // restore event listeners used to clear an inferior dim sub-mesh modified by viscous layers
557 void RestoreListeners();
559 // computes SMESH_ProxyMesh::SubMesh::_n2n;
560 bool MakeN2NMap( _MeshOfSolid* pm );
564 bool findSolidsWithLayers();
565 bool findFacesWithLayers();
566 bool makeLayer(_SolidData& data);
567 bool setEdgeData(_LayerEdge& edge, const set<TGeomID>& subIds,
568 SMESH_MesherHelper& helper, _SolidData& data);
569 gp_XYZ getFaceNormal(const SMDS_MeshNode* n,
570 const TopoDS_Face& face,
571 SMESH_MesherHelper& helper,
573 bool shiftInside=false);
574 gp_XYZ getWeigthedNormal( const SMDS_MeshNode* n,
575 std::pair< TGeomID, gp_XYZ > fId2Normal[],
577 bool findNeiborsOnEdge(const _LayerEdge* edge,
578 const SMDS_MeshNode*& n1,
579 const SMDS_MeshNode*& n2,
581 void getSimplices( const SMDS_MeshNode* node, vector<_Simplex>& simplices,
582 const set<TGeomID>& ingnoreShapes,
583 const _SolidData* dataToCheckOri = 0,
584 const bool toSort = false);
585 void findSimplexTestEdges( _SolidData& data,
586 vector< vector<_LayerEdge*> >& edgesByGeom);
587 bool sortEdges( _SolidData& data,
588 vector< vector<_LayerEdge*> >& edgesByGeom);
589 void limitStepSizeByCurvature( _SolidData& data );
590 void limitStepSize( _SolidData& data,
591 const SMDS_MeshElement* face,
592 const _LayerEdge* maxCosinEdge );
593 void limitStepSize( _SolidData& data, const double minSize);
594 bool inflate(_SolidData& data);
595 bool smoothAndCheck(_SolidData& data, const int nbSteps, double & distToIntersection);
596 bool smoothAnalyticEdge( _SolidData& data,
599 Handle(Geom_Surface)& surface,
600 const TopoDS_Face& F,
601 SMESH_MesherHelper& helper);
602 bool updateNormals( _SolidData& data, SMESH_MesherHelper& helper, int stepNb );
603 bool updateNormalsOfConvexFaces( _SolidData& data,
604 SMESH_MesherHelper& helper,
606 bool refine(_SolidData& data);
608 bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F,
609 SMESH_MesherHelper& helper,
610 const SMESHDS_SubMesh* faceSubMesh );
611 void fixBadFaces(const TopoDS_Face& F,
612 SMESH_MesherHelper& helper,
615 set<const SMDS_MeshNode*> * involvedNodes=NULL);
616 bool addBoundaryElements();
618 bool error( const string& text, int solidID=-1 );
619 SMESHDS_Mesh* getMeshDS() { return _mesh->GetMeshDS(); }
622 void makeGroupOfLE();
625 SMESH_ComputeErrorPtr _error;
627 vector< _SolidData > _sdVec;
630 //--------------------------------------------------------------------------------
632 * \brief Shrinker of nodes on the EDGE
636 vector<double> _initU;
637 vector<double> _normPar;
638 vector<const SMDS_MeshNode*> _nodes;
639 const _LayerEdge* _edges[2];
642 void AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper );
643 void Compute(bool set3D, SMESH_MesherHelper& helper);
644 void RestoreParams();
645 void SwapSrcTgtNodes(SMESHDS_Mesh* mesh);
647 //--------------------------------------------------------------------------------
649 * \brief Class of temporary mesh face.
650 * We can't use SMDS_FaceOfNodes since it's impossible to set it's ID which is
651 * needed because SMESH_ElementSearcher internaly uses set of elements sorted by ID
653 struct _TmpMeshFace : public SMDS_MeshElement
655 vector<const SMDS_MeshNode* > _nn;
656 _TmpMeshFace( const vector<const SMDS_MeshNode*>& nodes, int id, int faceID=-1):
657 SMDS_MeshElement(id), _nn(nodes) { setShapeId(faceID); }
658 virtual const SMDS_MeshNode* GetNode(const int ind) const { return _nn[ind]; }
659 virtual SMDSAbs_ElementType GetType() const { return SMDSAbs_Face; }
660 virtual vtkIdType GetVtkType() const { return -1; }
661 virtual SMDSAbs_EntityType GetEntityType() const { return SMDSEntity_Last; }
662 virtual SMDSAbs_GeometryType GetGeomType() const { return SMDSGeom_TRIANGLE; }
663 virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType) const
664 { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));}
666 //--------------------------------------------------------------------------------
668 * \brief Class of temporary mesh face storing _LayerEdge it's based on
670 struct _TmpMeshFaceOnEdge : public _TmpMeshFace
672 _LayerEdge *_le1, *_le2;
673 _TmpMeshFaceOnEdge( _LayerEdge* le1, _LayerEdge* le2, int ID ):
674 _TmpMeshFace( vector<const SMDS_MeshNode*>(4), ID ), _le1(le1), _le2(le2)
676 _nn[0]=_le1->_nodes[0];
677 _nn[1]=_le1->_nodes.back();
678 _nn[2]=_le2->_nodes.back();
679 _nn[3]=_le2->_nodes[0];
682 //--------------------------------------------------------------------------------
684 * \brief Retriever of node coordinates either directly of from a surface by node UV.
685 * \warning Location of a surface is ignored
687 struct _NodeCoordHelper
689 SMESH_MesherHelper& _helper;
690 const TopoDS_Face& _face;
691 Handle(Geom_Surface) _surface;
692 gp_XYZ (_NodeCoordHelper::* _fun)(const SMDS_MeshNode* n) const;
694 _NodeCoordHelper(const TopoDS_Face& F, SMESH_MesherHelper& helper, bool is2D)
695 : _helper( helper ), _face( F )
700 _surface = BRep_Tool::Surface( _face, loc );
702 if ( _surface.IsNull() )
703 _fun = & _NodeCoordHelper::direct;
705 _fun = & _NodeCoordHelper::byUV;
707 gp_XYZ operator()(const SMDS_MeshNode* n) const { return (this->*_fun)( n ); }
710 gp_XYZ direct(const SMDS_MeshNode* n) const
712 return SMESH_TNodeXYZ( n );
714 gp_XYZ byUV (const SMDS_MeshNode* n) const
716 gp_XY uv = _helper.GetNodeUV( _face, n );
717 return _surface->Value( uv.X(), uv.Y() ).XYZ();
720 } // namespace VISCOUS_3D
724 //================================================================================
725 // StdMeshers_ViscousLayers hypothesis
727 StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen)
728 :SMESH_Hypothesis(hypId, studyId, gen),
729 _isToIgnoreShapes(1), _nbLayers(1), _thickness(1), _stretchFactor(1)
731 _name = StdMeshers_ViscousLayers::GetHypType();
732 _param_algo_dim = -3; // auxiliary hyp used by 3D algos
733 } // --------------------------------------------------------------------------------
734 void StdMeshers_ViscousLayers::SetBndShapes(const std::vector<int>& faceIds, bool toIgnore)
736 if ( faceIds != _shapeIds )
737 _shapeIds = faceIds, NotifySubMeshesHypothesisModification();
738 if ( _isToIgnoreShapes != toIgnore )
739 _isToIgnoreShapes = toIgnore, NotifySubMeshesHypothesisModification();
740 } // --------------------------------------------------------------------------------
741 void StdMeshers_ViscousLayers::SetTotalThickness(double thickness)
743 if ( thickness != _thickness )
744 _thickness = thickness, NotifySubMeshesHypothesisModification();
745 } // --------------------------------------------------------------------------------
746 void StdMeshers_ViscousLayers::SetNumberLayers(int nb)
748 if ( _nbLayers != nb )
749 _nbLayers = nb, NotifySubMeshesHypothesisModification();
750 } // --------------------------------------------------------------------------------
751 void StdMeshers_ViscousLayers::SetStretchFactor(double factor)
753 if ( _stretchFactor != factor )
754 _stretchFactor = factor, NotifySubMeshesHypothesisModification();
755 } // --------------------------------------------------------------------------------
757 StdMeshers_ViscousLayers::Compute(SMESH_Mesh& theMesh,
758 const TopoDS_Shape& theShape,
759 const bool toMakeN2NMap) const
761 using namespace VISCOUS_3D;
762 _ViscousBuilder bulder;
763 SMESH_ComputeErrorPtr err = bulder.Compute( theMesh, theShape );
764 if ( err && !err->IsOK() )
765 return SMESH_ProxyMesh::Ptr();
767 vector<SMESH_ProxyMesh::Ptr> components;
768 TopExp_Explorer exp( theShape, TopAbs_SOLID );
769 for ( ; exp.More(); exp.Next() )
771 if ( _MeshOfSolid* pm =
772 _ViscousListener::GetSolidMesh( &theMesh, exp.Current(), /*toCreate=*/false))
774 if ( toMakeN2NMap && !pm->_n2nMapComputed )
775 if ( !bulder.MakeN2NMap( pm ))
776 return SMESH_ProxyMesh::Ptr();
777 components.push_back( SMESH_ProxyMesh::Ptr( pm ));
778 pm->myIsDeletable = false; // it will de deleted by boost::shared_ptr
780 _ViscousListener::RemoveSolidMesh ( &theMesh, exp.Current() );
782 switch ( components.size() )
786 case 1: return components[0];
788 default: return SMESH_ProxyMesh::Ptr( new SMESH_ProxyMesh( components ));
790 return SMESH_ProxyMesh::Ptr();
791 } // --------------------------------------------------------------------------------
792 std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save)
794 save << " " << _nbLayers
796 << " " << _stretchFactor
797 << " " << _shapeIds.size();
798 for ( size_t i = 0; i < _shapeIds.size(); ++i )
799 save << " " << _shapeIds[i];
800 save << " " << !_isToIgnoreShapes; // negate to keep the behavior in old studies.
802 } // --------------------------------------------------------------------------------
803 std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load)
805 int nbFaces, faceID, shapeToTreat;
806 load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces;
807 while ( _shapeIds.size() < nbFaces && load >> faceID )
808 _shapeIds.push_back( faceID );
809 if ( load >> shapeToTreat )
810 _isToIgnoreShapes = !shapeToTreat;
812 _isToIgnoreShapes = true; // old behavior
814 } // --------------------------------------------------------------------------------
815 bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh* theMesh,
816 const TopoDS_Shape& theShape)
821 // END StdMeshers_ViscousLayers hypothesis
822 //================================================================================
826 gp_XYZ getEdgeDir( const TopoDS_Edge& E, const TopoDS_Vertex& fromV )
830 Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
831 gp_Pnt p = BRep_Tool::Pnt( fromV );
832 double distF = p.SquareDistance( c->Value( f ));
833 double distL = p.SquareDistance( c->Value( l ));
834 c->D1(( distF < distL ? f : l), p, dir );
835 if ( distL < distF ) dir.Reverse();
838 //--------------------------------------------------------------------------------
839 gp_XYZ getEdgeDir( const TopoDS_Edge& E, const SMDS_MeshNode* atNode,
840 SMESH_MesherHelper& helper)
843 double f,l; gp_Pnt p;
844 Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l );
845 double u = helper.GetNodeU( E, atNode );
849 //--------------------------------------------------------------------------------
850 gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE,
851 const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok)
853 gp_XY uv = helper.GetNodeUV( F, node, 0, &ok );
854 Handle(Geom_Surface) surface = BRep_Tool::Surface( F );
855 gp_Pnt p; gp_Vec du, dv, norm;
856 surface->D1( uv.X(),uv.Y(), p, du,dv );
860 Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l );
861 double u = helper.GetNodeU( fromE, node, 0, &ok );
863 TopAbs_Orientation o = helper.GetSubShapeOri( F.Oriented(TopAbs_FORWARD), fromE);
864 if ( o == TopAbs_REVERSED )
867 gp_Vec dir = norm ^ du;
869 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX &&
870 helper.IsClosedEdge( fromE ))
872 if ( fabs(u-f) < fabs(u-l)) c->D1( l, p, dv );
873 else c->D1( f, p, dv );
874 if ( o == TopAbs_REVERSED )
876 gp_Vec dir2 = norm ^ dv;
877 dir = dir.Normalized() + dir2.Normalized();
881 //--------------------------------------------------------------------------------
882 gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV,
883 const SMDS_MeshNode* node, SMESH_MesherHelper& helper,
884 bool& ok, double* cosin=0)
886 TopoDS_Face faceFrw = F;
887 faceFrw.Orientation( TopAbs_FORWARD );
888 double f,l; TopLoc_Location loc;
889 TopoDS_Edge edges[2]; // sharing a vertex
893 TopExp_Explorer exp( faceFrw, TopAbs_EDGE );
894 for ( ; exp.More() && nbEdges < 2; exp.Next() )
896 const TopoDS_Edge& e = TopoDS::Edge( exp.Current() );
897 if ( SMESH_Algo::isDegenerated( e )) continue;
898 TopExp::Vertices( e, VV[0], VV[1], /*CumOri=*/true );
899 if ( VV[1].IsSame( fromV )) {
900 nbEdges += edges[ 0 ].IsNull();
903 else if ( VV[0].IsSame( fromV )) {
904 nbEdges += edges[ 1 ].IsNull();
909 gp_XYZ dir(0,0,0), edgeDir[2];
912 // get dirs of edges going fromV
914 for ( size_t i = 0; i < nbEdges && ok; ++i )
916 edgeDir[i] = getEdgeDir( edges[i], fromV );
917 double size2 = edgeDir[i].SquareModulus();
918 if (( ok = size2 > numeric_limits<double>::min() ))
919 edgeDir[i] /= sqrt( size2 );
921 if ( !ok ) return dir;
923 // get angle between the 2 edges
925 double angle = helper.GetAngle( edges[0], edges[1], faceFrw, &faceNormal );
926 if ( Abs( angle ) < 5 * M_PI/180 )
928 dir = ( faceNormal.XYZ() ^ edgeDir[0].Reversed()) + ( faceNormal.XYZ() ^ edgeDir[1] );
932 dir = edgeDir[0] + edgeDir[1];
937 double angle = gp_Vec( edgeDir[0] ).Angle( dir );
938 *cosin = Cos( angle );
941 else if ( nbEdges == 1 )
943 dir = getFaceDir( faceFrw, edges[ edges[0].IsNull() ], node, helper, ok );
944 if ( cosin ) *cosin = 1.;
953 //================================================================================
955 * \brief Returns true if a FACE is bound by a concave EDGE
957 //================================================================================
959 bool isConcave( const TopoDS_Face& F, SMESH_MesherHelper& helper )
961 // if ( helper.Count( F, TopAbs_WIRE, /*useMap=*/false) > 1 )
965 TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE );
966 for ( ; eExp.More(); eExp.Next() )
968 const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
969 if ( SMESH_Algo::isDegenerated( E )) continue;
970 // check if 2D curve is concave
971 BRepAdaptor_Curve2d curve( E, F );
972 const int nbIntervals = curve.NbIntervals( GeomAbs_C2 );
973 TColStd_Array1OfReal intervals(1, nbIntervals + 1 );
974 curve.Intervals( intervals, GeomAbs_C2 );
975 bool isConvex = true;
976 for ( int i = 1; i <= nbIntervals && isConvex; ++i )
978 double u1 = intervals( i );
979 double u2 = intervals( i+1 );
980 curve.D2( 0.5*( u1+u2 ), p, drv1, drv2 );
981 double cross = drv2 ^ drv1;
982 if ( E.Orientation() == TopAbs_REVERSED )
984 isConvex = ( cross > 0.1 ); //-1e-9 );
988 //cout << "Concave FACE " << helper.GetMeshDS()->ShapeToIndex( F ) << endl;
992 // check angles at VERTEXes
994 TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *helper.GetMesh(), 0, error );
995 for ( size_t iW = 0; iW < wires.size(); ++iW )
997 const int nbEdges = wires[iW]->NbEdges();
998 if ( nbEdges < 2 && SMESH_Algo::isDegenerated( wires[iW]->Edge(0)))
1000 for ( int iE1 = 0; iE1 < nbEdges; ++iE1 )
1002 if ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE1 ))) continue;
1003 int iE2 = ( iE1 + 1 ) % nbEdges;
1004 while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 )))
1005 iE2 = ( iE2 + 1 ) % nbEdges;
1006 double angle = helper.GetAngle( wires[iW]->Edge( iE1 ),
1007 wires[iW]->Edge( iE2 ), F );
1008 if ( angle < -5. * M_PI / 180. )
1014 //--------------------------------------------------------------------------------
1015 // DEBUG. Dump intermediate node positions into a python script
1016 // HOWTO use: run python commands written in a console to see
1017 // construction steps of viscous layers
1023 const char* fname = "/tmp/viscous.py";
1024 cout << "execfile('"<<fname<<"')"<<endl;
1025 py = new ofstream(fname);
1026 *py << "import SMESH" << endl
1027 << "from salome.smesh import smeshBuilder" << endl
1028 << "smesh = smeshBuilder.New(salome.myStudy)" << endl
1029 << "meshSO = smesh.GetCurrentStudy().FindObjectID('0:1:2:3')" << endl
1030 << "mesh = smesh.Mesh( meshSO.GetObject() )"<<endl;
1035 *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Viscous Prisms',"
1036 "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_ElemGeomType,'=',SMESH.Geom_PENTA))"<<endl;
1037 *py << "mesh.GroupOnFilter(SMESH.VOLUME,'Neg Volumes',"
1038 "smesh.GetFilter(SMESH.VOLUME,SMESH.FT_Volume3D,'<',0))"<<endl;
1042 ~PyDump() { Finish(); cout << "NB FUNCTIONS: " << theNbFunc << endl; }
1044 #define dumpFunction(f) { _dumpFunction(f, __LINE__);}
1045 #define dumpMove(n) { _dumpMove(n, __LINE__);}
1046 #define dumpCmd(txt) { _dumpCmd(txt, __LINE__);}
1047 void _dumpFunction(const string& fun, int ln)
1048 { if (py) *py<< "def "<<fun<<"(): # "<< ln <<endl; cout<<fun<<"()"<<endl; ++theNbFunc; }
1049 void _dumpMove(const SMDS_MeshNode* n, int ln)
1050 { if (py) *py<< " mesh.MoveNode( "<<n->GetID()<< ", "<< n->X()
1051 << ", "<<n->Y()<<", "<< n->Z()<< ")\t\t # "<< ln <<endl; }
1052 void _dumpCmd(const string& txt, int ln)
1053 { if (py) *py<< " "<<txt<<" # "<< ln <<endl; }
1054 void dumpFunctionEnd()
1055 { if (py) *py<< " return"<< endl; }
1056 void dumpChangeNodes( const SMDS_MeshElement* f )
1057 { if (py) { *py<< " mesh.ChangeElemNodes( " << f->GetID()<<", [";
1058 for ( int i=1; i < f->NbNodes(); ++i ) *py << f->GetNode(i-1)->GetID()<<", ";
1059 *py << f->GetNode( f->NbNodes()-1 )->GetID() << " ])"<< endl; }}
1060 #define debugMsg( txt ) { cout << txt << " (line: " << __LINE__ << ")" << endl; }
1062 struct PyDump { void Finish() {} };
1063 #define dumpFunction(f) f
1065 #define dumpCmd(txt)
1066 #define dumpFunctionEnd()
1067 #define dumpChangeNodes(f)
1068 #define debugMsg( txt ) {}
1072 using namespace VISCOUS_3D;
1074 //================================================================================
1076 * \brief Constructor of _ViscousBuilder
1078 //================================================================================
1080 _ViscousBuilder::_ViscousBuilder()
1082 _error = SMESH_ComputeError::New(COMPERR_OK);
1086 //================================================================================
1088 * \brief Stores error description and returns false
1090 //================================================================================
1092 bool _ViscousBuilder::error(const string& text, int solidId )
1094 _error->myName = COMPERR_ALGO_FAILED;
1095 _error->myComment = string("Viscous layers builder: ") + text;
1098 SMESH_subMesh* sm = _mesh->GetSubMeshContaining( solidId );
1099 if ( !sm && !_sdVec.empty() )
1100 sm = _mesh->GetSubMeshContaining( _sdVec[0]._index );
1101 if ( sm && sm->GetSubShape().ShapeType() == TopAbs_SOLID )
1103 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1104 if ( smError && smError->myAlgo )
1105 _error->myAlgo = smError->myAlgo;
1109 makeGroupOfLE(); // debug
1114 //================================================================================
1116 * \brief At study restoration, restore event listeners used to clear an inferior
1117 * dim sub-mesh modified by viscous layers
1119 //================================================================================
1121 void _ViscousBuilder::RestoreListeners()
1126 //================================================================================
1128 * \brief computes SMESH_ProxyMesh::SubMesh::_n2n
1130 //================================================================================
1132 bool _ViscousBuilder::MakeN2NMap( _MeshOfSolid* pm )
1134 SMESH_subMesh* solidSM = pm->mySubMeshes.front();
1135 TopExp_Explorer fExp( solidSM->GetSubShape(), TopAbs_FACE );
1136 for ( ; fExp.More(); fExp.Next() )
1138 SMESHDS_SubMesh* srcSmDS = pm->GetMeshDS()->MeshElements( fExp.Current() );
1139 const SMESH_ProxyMesh::SubMesh* prxSmDS = pm->GetProxySubMesh( fExp.Current() );
1141 if ( !srcSmDS || !prxSmDS || !srcSmDS->NbElements() || !prxSmDS->NbElements() )
1143 if ( srcSmDS->GetElements()->next() == prxSmDS->GetElements()->next())
1146 if ( srcSmDS->NbElements() != prxSmDS->NbElements() )
1147 return error( "Different nb elements in a source and a proxy sub-mesh", solidSM->GetId());
1149 SMDS_ElemIteratorPtr srcIt = srcSmDS->GetElements();
1150 SMDS_ElemIteratorPtr prxIt = prxSmDS->GetElements();
1151 while( prxIt->more() )
1153 const SMDS_MeshElement* fSrc = srcIt->next();
1154 const SMDS_MeshElement* fPrx = prxIt->next();
1155 if ( fSrc->NbNodes() != fPrx->NbNodes())
1156 return error( "Different elements in a source and a proxy sub-mesh", solidSM->GetId());
1157 for ( int i = 0 ; i < fPrx->NbNodes(); ++i )
1158 pm->setNode2Node( fSrc->GetNode(i), fPrx->GetNode(i), prxSmDS );
1161 pm->_n2nMapComputed = true;
1165 //================================================================================
1167 * \brief Does its job
1169 //================================================================================
1171 SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh,
1172 const TopoDS_Shape& theShape)
1174 // TODO: set priority of solids during Gen::Compute()
1178 // check if proxy mesh already computed
1179 TopExp_Explorer exp( theShape, TopAbs_SOLID );
1181 return error("No SOLID's in theShape"), _error;
1183 if ( _ViscousListener::GetSolidMesh( _mesh, exp.Current(), /*toCreate=*/false))
1184 return SMESH_ComputeErrorPtr(); // everything already computed
1188 // TODO: ignore already computed SOLIDs
1189 if ( !findSolidsWithLayers())
1192 if ( !findFacesWithLayers() )
1195 for ( size_t i = 0; i < _sdVec.size(); ++i )
1197 if ( ! makeLayer(_sdVec[i]) )
1200 if ( _sdVec[i]._edges.size() == 0 )
1203 if ( ! inflate(_sdVec[i]) )
1206 if ( ! refine(_sdVec[i]) )
1212 addBoundaryElements();
1214 makeGroupOfLE(); // debug
1220 //================================================================================
1222 * \brief Finds SOLIDs to compute using viscous layers. Fills _sdVec
1224 //================================================================================
1226 bool _ViscousBuilder::findSolidsWithLayers()
1229 TopTools_IndexedMapOfShape allSolids;
1230 TopExp::MapShapes( _mesh->GetShapeToMesh(), TopAbs_SOLID, allSolids );
1231 _sdVec.reserve( allSolids.Extent());
1233 SMESH_Gen* gen = _mesh->GetGen();
1234 SMESH_HypoFilter filter;
1235 for ( int i = 1; i <= allSolids.Extent(); ++i )
1237 // find StdMeshers_ViscousLayers hyp assigned to the i-th solid
1238 SMESH_Algo* algo = gen->GetAlgo( *_mesh, allSolids(i) );
1239 if ( !algo ) continue;
1240 // TODO: check if algo is hidden
1241 const list <const SMESHDS_Hypothesis *> & allHyps =
1242 algo->GetUsedHypothesis(*_mesh, allSolids(i), /*ignoreAuxiliary=*/false);
1243 list< const SMESHDS_Hypothesis *>::const_iterator hyp = allHyps.begin();
1244 const StdMeshers_ViscousLayers* viscHyp = 0;
1245 for ( ; hyp != allHyps.end() && !viscHyp; ++hyp )
1246 viscHyp = dynamic_cast<const StdMeshers_ViscousLayers*>( *hyp );
1249 TopoDS_Shape hypShape;
1250 filter.Init( filter.Is( viscHyp ));
1251 _mesh->GetHypothesis( allSolids(i), filter, true, &hypShape );
1253 _MeshOfSolid* proxyMesh = _ViscousListener::GetSolidMesh( _mesh,
1256 _sdVec.push_back( _SolidData( allSolids(i), viscHyp, hypShape, proxyMesh ));
1257 _sdVec.back()._index = getMeshDS()->ShapeToIndex( allSolids(i));
1260 if ( _sdVec.empty() )
1262 ( SMESH_Comment(StdMeshers_ViscousLayers::GetHypType()) << " hypothesis not found",0);
1267 //================================================================================
1271 //================================================================================
1273 bool _ViscousBuilder::findFacesWithLayers()
1275 SMESH_MesherHelper helper( *_mesh );
1276 TopExp_Explorer exp;
1277 TopTools_IndexedMapOfShape solids;
1279 // collect all faces to ignore defined by hyp
1280 for ( size_t i = 0; i < _sdVec.size(); ++i )
1282 solids.Add( _sdVec[i]._solid );
1284 vector<TGeomID> ids = _sdVec[i]._hyp->GetBndShapes();
1285 if ( _sdVec[i]._hyp->IsToIgnoreShapes() ) // FACEs to ignore are given
1287 for ( size_t ii = 0; ii < ids.size(); ++ii )
1289 const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[ii] );
1290 if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE )
1291 _sdVec[i]._ignoreFaceIds.insert( ids[ii] );
1294 else // FACEs with layers are given
1296 exp.Init( _sdVec[i]._solid, TopAbs_FACE );
1297 for ( ; exp.More(); exp.Next() )
1299 TGeomID faceInd = getMeshDS()->ShapeToIndex( exp.Current() );
1300 if ( find( ids.begin(), ids.end(), faceInd ) == ids.end() )
1301 _sdVec[i]._ignoreFaceIds.insert( faceInd );
1305 // ignore internal FACEs if inlets and outlets are specified
1307 TopTools_IndexedDataMapOfShapeListOfShape solidsOfFace;
1308 if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
1309 TopExp::MapShapesAndAncestors( _sdVec[i]._hypShape,
1310 TopAbs_FACE, TopAbs_SOLID, solidsOfFace);
1312 exp.Init( _sdVec[i]._solid.Oriented( TopAbs_FORWARD ), TopAbs_FACE );
1313 for ( ; exp.More(); exp.Next() )
1315 const TopoDS_Face& face = TopoDS::Face( exp.Current() );
1316 if ( helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) < 2 )
1319 const TGeomID faceInd = getMeshDS()->ShapeToIndex( face );
1320 if ( _sdVec[i]._hyp->IsToIgnoreShapes() )
1322 int nbSolids = solidsOfFace.FindFromKey( face ).Extent();
1324 _sdVec[i]._ignoreFaceIds.insert( faceInd );
1327 if ( helper.IsReversedSubMesh( face ))
1329 _sdVec[i]._reversedFaceIds.insert( faceInd );
1335 // Find faces to shrink mesh on (solution 2 in issue 0020832);
1336 TopTools_IndexedMapOfShape shapes;
1337 for ( size_t i = 0; i < _sdVec.size(); ++i )
1340 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_EDGE, shapes);
1341 for ( int iE = 1; iE <= shapes.Extent(); ++iE )
1343 const TopoDS_Shape& edge = shapes(iE);
1344 // find 2 faces sharing an edge
1346 PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE);
1347 while ( fIt->more())
1349 const TopoDS_Shape* f = fIt->next();
1350 if ( helper.IsSubShape( *f, _sdVec[i]._solid))
1351 FF[ int( !FF[0].IsNull()) ] = *f;
1353 if( FF[1].IsNull() ) continue; // seam edge can be shared by 1 FACE only
1354 // check presence of layers on them
1356 for ( int j = 0; j < 2; ++j )
1357 ignore[j] = _sdVec[i]._ignoreFaceIds.count ( getMeshDS()->ShapeToIndex( FF[j] ));
1358 if ( ignore[0] == ignore[1] )
1359 continue; // nothing interesting
1360 TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ];
1361 // check presence of layers on fWOL within an adjacent SOLID
1362 PShapeIteratorPtr sIt = helper.GetAncestors( fWOL, *_mesh, TopAbs_SOLID );
1363 while ( const TopoDS_Shape* solid = sIt->next() )
1364 if ( !solid->IsSame( _sdVec[i]._solid ))
1366 int iSolid = solids.FindIndex( *solid );
1367 int iFace = getMeshDS()->ShapeToIndex( fWOL );
1368 if ( iSolid > 0 && !_sdVec[ iSolid-1 ]._ignoreFaceIds.count( iFace ))
1370 _sdVec[i]._noShrinkFaces.insert( iFace );
1375 if ( !fWOL.IsNull())
1377 TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge );
1378 _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL ));
1382 // Exclude from _shrinkShape2Shape FACE's that can't be shrinked since
1383 // the algo of the SOLID sharing the FACE does not support it
1384 set< string > notSupportAlgos; notSupportAlgos.insert("Hexa_3D");
1385 for ( size_t i = 0; i < _sdVec.size(); ++i )
1387 TopTools_MapOfShape noShrinkVertices;
1388 map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin();
1389 for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f )
1391 const TopoDS_Shape& fWOL = e2f->second;
1392 TGeomID edgeID = e2f->first;
1393 bool notShrinkFace = false;
1394 PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID);
1395 while ( soIt->more())
1397 const TopoDS_Shape* solid = soIt->next();
1398 if ( _sdVec[i]._solid.IsSame( *solid )) continue;
1399 SMESH_Algo* algo = _mesh->GetGen()->GetAlgo( *_mesh, *solid );
1400 if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue;
1401 notShrinkFace = true;
1402 for ( size_t j = 0; j < _sdVec.size(); ++j )
1404 if ( _sdVec[j]._solid.IsSame( *solid ) )
1405 if ( _sdVec[j]._shrinkShape2Shape.count( edgeID ))
1406 notShrinkFace = false;
1409 if ( notShrinkFace )
1411 _sdVec[i]._noShrinkFaces.insert( getMeshDS()->ShapeToIndex( fWOL ));
1412 for ( TopExp_Explorer vExp( fWOL, TopAbs_VERTEX ); vExp.More(); vExp.Next() )
1413 noShrinkVertices.Add( vExp.Current() );
1416 // erase from _shrinkShape2Shape all srink EDGE's of a SOLID connected
1417 // to the found not shrinked fWOL's
1418 e2f = _sdVec[i]._shrinkShape2Shape.begin();
1419 for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); )
1421 TGeomID edgeID = e2f->first;
1422 TopoDS_Vertex VV[2];
1423 TopExp::Vertices( TopoDS::Edge( getMeshDS()->IndexToShape( edgeID )),VV[0],VV[1]);
1424 if ( noShrinkVertices.Contains( VV[0] ) || noShrinkVertices.Contains( VV[1] ))
1426 _sdVec[i]._noShrinkFaces.insert( getMeshDS()->ShapeToIndex( e2f->second ));
1427 _sdVec[i]._shrinkShape2Shape.erase( e2f++ );
1436 // Find the SHAPE along which to inflate _LayerEdge based on VERTEX
1438 for ( size_t i = 0; i < _sdVec.size(); ++i )
1441 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_VERTEX, shapes);
1442 for ( int iV = 1; iV <= shapes.Extent(); ++iV )
1444 const TopoDS_Shape& vertex = shapes(iV);
1445 // find faces WOL sharing the vertex
1446 vector< TopoDS_Shape > facesWOL;
1447 int totalNbFaces = 0;
1448 PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE);
1449 while ( fIt->more())
1451 const TopoDS_Shape* f = fIt->next();
1452 if ( helper.IsSubShape( *f, _sdVec[i]._solid ) )
1455 const int fID = getMeshDS()->ShapeToIndex( *f );
1456 if ( _sdVec[i]._ignoreFaceIds.count ( fID ) &&
1457 !_sdVec[i]._noShrinkFaces.count( fID ))
1458 facesWOL.push_back( *f );
1461 if ( facesWOL.size() == totalNbFaces || facesWOL.empty() )
1462 continue; // no layers at this vertex or no WOL
1463 TGeomID vInd = getMeshDS()->ShapeToIndex( vertex );
1464 switch ( facesWOL.size() )
1468 helper.SetSubShape( facesWOL[0] );
1469 if ( helper.IsRealSeam( vInd )) // inflate along a seam edge?
1471 TopoDS_Shape seamEdge;
1472 PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1473 while ( eIt->more() && seamEdge.IsNull() )
1475 const TopoDS_Shape* e = eIt->next();
1476 if ( helper.IsRealSeam( *e ) )
1479 if ( !seamEdge.IsNull() )
1481 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, seamEdge ));
1485 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, facesWOL[0] ));
1490 // find an edge shared by 2 faces
1491 PShapeIteratorPtr eIt = helper.GetAncestors(vertex, *_mesh, TopAbs_EDGE);
1492 while ( eIt->more())
1494 const TopoDS_Shape* e = eIt->next();
1495 if ( helper.IsSubShape( *e, facesWOL[0]) &&
1496 helper.IsSubShape( *e, facesWOL[1]))
1498 _sdVec[i]._shrinkShape2Shape.insert( make_pair( vInd, *e )); break;
1504 return error("Not yet supported case", _sdVec[i]._index);
1509 // add FACEs of other SOLIDs to _ignoreFaceIds
1510 for ( size_t i = 0; i < _sdVec.size(); ++i )
1513 TopExp::MapShapes(_sdVec[i]._solid, TopAbs_FACE, shapes);
1515 for ( exp.Init( _mesh->GetShapeToMesh(), TopAbs_FACE ); exp.More(); exp.Next() )
1517 if ( !shapes.Contains( exp.Current() ))
1518 _sdVec[i]._ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( exp.Current() ));
1525 //================================================================================
1527 * \brief Create the inner surface of the viscous layer and prepare data for infation
1529 //================================================================================
1531 bool _ViscousBuilder::makeLayer(_SolidData& data)
1533 // get all sub-shapes to make layers on
1534 set<TGeomID> subIds, faceIds;
1535 subIds = data._noShrinkFaces;
1536 TopExp_Explorer exp( data._solid, TopAbs_FACE );
1537 for ( ; exp.More(); exp.Next() )
1539 SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() );
1540 if ( ! data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
1541 faceIds.insert( fSubM->GetId() );
1542 SMESH_subMeshIteratorPtr subIt =
1543 fSubM->getDependsOnIterator(/*includeSelf=*/true, /*complexShapeFirst=*/false);
1544 while ( subIt->more() )
1545 subIds.insert( subIt->next()->GetId() );
1548 // make a map to find new nodes on sub-shapes shared with other SOLID
1549 map< TGeomID, TNode2Edge* >::iterator s2ne;
1550 map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
1551 for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
1553 TGeomID shapeInd = s2s->first;
1554 for ( size_t i = 0; i < _sdVec.size(); ++i )
1556 if ( _sdVec[i]._index == data._index ) continue;
1557 map< TGeomID, TopoDS_Shape >::iterator s2s2 = _sdVec[i]._shrinkShape2Shape.find( shapeInd );
1558 if ( s2s2 != _sdVec[i]._shrinkShape2Shape.end() &&
1559 *s2s == *s2s2 && !_sdVec[i]._n2eMap.empty() )
1561 data._s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap ));
1567 // Create temporary faces and _LayerEdge's
1569 dumpFunction(SMESH_Comment("makeLayers_")<<data._index);
1571 data._stepSize = Precision::Infinite();
1572 data._stepSizeNodes[0] = 0;
1574 SMESH_MesherHelper helper( *_mesh );
1575 helper.SetSubShape( data._solid );
1576 helper.SetElementsOnShape(true);
1578 vector< const SMDS_MeshNode*> newNodes; // of a mesh face
1579 TNode2Edge::iterator n2e2;
1581 // collect _LayerEdge's of shapes they are based on
1582 const int nbShapes = getMeshDS()->MaxShapeIndex();
1583 vector< vector<_LayerEdge*> > edgesByGeom( nbShapes+1 );
1585 for ( set<TGeomID>::iterator id = faceIds.begin(); id != faceIds.end(); ++id )
1587 SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( *id );
1588 if ( !smDS ) return error(SMESH_Comment("Not meshed face ") << *id, data._index );
1590 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( *id ));
1591 SMESH_ProxyMesh::SubMesh* proxySub =
1592 data._proxyMesh->getFaceSubM( F, /*create=*/true);
1594 SMDS_ElemIteratorPtr eIt = smDS->GetElements();
1595 while ( eIt->more() )
1597 const SMDS_MeshElement* face = eIt->next();
1598 newNodes.resize( face->NbCornerNodes() );
1599 double faceMaxCosin = -1;
1600 _LayerEdge* maxCosinEdge = 0;
1601 for ( int i = 0 ; i < face->NbCornerNodes(); ++i )
1603 const SMDS_MeshNode* n = face->GetNode(i);
1604 TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first;
1605 if ( !(*n2e).second )
1608 _LayerEdge* edge = new _LayerEdge();
1610 edge->_nodes.push_back( n );
1611 const int shapeID = n->getshapeId();
1612 edgesByGeom[ shapeID ].push_back( edge );
1614 SMESH_TNodeXYZ xyz( n );
1616 // set edge data or find already refined _LayerEdge and get data from it
1617 if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
1618 ( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() &&
1619 ( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end())
1621 _LayerEdge* foundEdge = (*n2e2).second;
1622 gp_XYZ lastPos = edge->Copy( *foundEdge, helper );
1623 foundEdge->_pos.push_back( lastPos );
1624 // location of the last node is modified and we restore it by foundEdge->_pos.back()
1625 const_cast< SMDS_MeshNode* >
1626 ( edge->_nodes.back() )->setXYZ( xyz.X(), xyz.Y(), xyz.Z() );
1630 edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ));
1631 if ( !setEdgeData( *edge, subIds, helper, data ))
1634 dumpMove(edge->_nodes.back());
1636 if ( edge->_cosin > faceMaxCosin )
1638 faceMaxCosin = edge->_cosin;
1639 maxCosinEdge = edge;
1642 newNodes[ i ] = n2e->second->_nodes.back();
1644 // create a temporary face
1645 const SMDS_MeshElement* newFace =
1646 new _TmpMeshFace( newNodes, --_tmpFaceID, face->getshapeId() );
1647 proxySub->AddElement( newFace );
1649 // compute inflation step size by min size of element on a convex surface
1650 if ( faceMaxCosin > theMinSmoothCosin )
1651 limitStepSize( data, face, maxCosinEdge );
1652 } // loop on 2D elements on a FACE
1653 } // loop on FACEs of a SOLID
1655 data._epsilon = 1e-7;
1656 if ( data._stepSize < 1. )
1657 data._epsilon *= data._stepSize;
1659 // Put _LayerEdge's into the vector data._edges
1660 if ( !sortEdges( data, edgesByGeom ))
1663 // limit data._stepSize depending on surface curvature and fill data._convexFaces
1664 limitStepSizeByCurvature( data ); // !!! it must be before node substitution in _Simplex
1666 // Set target nodes into _Simplex and _2NearEdges of _LayerEdge's
1667 TNode2Edge::iterator n2e;
1668 for ( size_t i = 0; i < data._edges.size(); ++i )
1670 if ( data._edges[i]->IsOnEdge())
1671 for ( int j = 0; j < 2; ++j )
1673 if ( data._edges[i]->_nodes.back()->NbInverseElements(SMDSAbs_Volume) > 0 )
1674 break; // _LayerEdge is shared by two _SolidData's
1675 const SMDS_MeshNode* & n = data._edges[i]->_2neibors->_nodes[j];
1676 if (( n2e = data._n2eMap.find( n )) == data._n2eMap.end() )
1677 return error("_LayerEdge not found by src node", data._index);
1678 n = (*n2e).second->_nodes.back();
1679 data._edges[i]->_2neibors->_edges[j] = n2e->second;
1682 for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j )
1684 _Simplex& s = data._edges[i]->_simplices[j];
1685 s._nNext = data._n2eMap[ s._nNext ]->_nodes.back();
1686 s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back();
1694 //================================================================================
1696 * \brief Compute inflation step size by min size of element on a convex surface
1698 //================================================================================
1700 void _ViscousBuilder::limitStepSize( _SolidData& data,
1701 const SMDS_MeshElement* face,
1702 const _LayerEdge* maxCosinEdge )
1705 double minSize = 10 * data._stepSize;
1706 const int nbNodes = face->NbCornerNodes();
1707 for ( int i = 0; i < nbNodes; ++i )
1709 const SMDS_MeshNode* nextN = face->GetNode( SMESH_MesherHelper::WrapIndex( i+1, nbNodes ));
1710 const SMDS_MeshNode* curN = face->GetNode( i );
1711 if ( nextN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ||
1712 curN-> GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE )
1714 double dist = SMESH_TNodeXYZ( curN ).Distance( nextN );
1715 if ( dist < minSize )
1716 minSize = dist, iN = i;
1719 double newStep = 0.8 * minSize / maxCosinEdge->_lenFactor;
1720 if ( newStep < data._stepSize )
1722 data._stepSize = newStep;
1723 data._stepSizeCoeff = 0.8 / maxCosinEdge->_lenFactor;
1724 data._stepSizeNodes[0] = face->GetNode( iN );
1725 data._stepSizeNodes[1] = face->GetNode( SMESH_MesherHelper::WrapIndex( iN+1, nbNodes ));
1729 //================================================================================
1731 * \brief Compute inflation step size by min size of element on a convex surface
1733 //================================================================================
1735 void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize )
1737 if ( minSize < data._stepSize )
1739 data._stepSize = minSize;
1740 if ( data._stepSizeNodes[0] )
1743 SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
1744 data._stepSizeCoeff = data._stepSize / dist;
1749 //================================================================================
1751 * \brief Limit data._stepSize by evaluating curvature of shapes and fill data._convexFaces
1753 //================================================================================
1755 void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data )
1757 const int nbTestPnt = 5; // on a FACE sub-shape
1758 const double minCurvature = 0.9 / data._hyp->GetTotalThickness();
1760 BRepLProp_SLProps surfProp( 2, 1e-6 );
1761 SMESH_MesherHelper helper( *_mesh );
1763 data._convexFaces.clear();
1765 TopExp_Explorer face( data._solid, TopAbs_FACE );
1766 for ( ; face.More(); face.Next() )
1768 const TopoDS_Face& F = TopoDS::Face( face.Current() );
1769 BRepAdaptor_Surface surface( F, false );
1770 surfProp.SetSurface( surface );
1772 bool isTooCurved = false;
1775 _ConvexFace cnvFace;
1776 SMESH_subMesh * sm = _mesh->GetSubMesh( F );
1777 const TGeomID faceID = sm->GetId();
1778 if ( data._ignoreFaceIds.count( faceID )) continue;
1779 const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. );
1780 SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
1781 while ( smIt->more() )
1784 const TGeomID subID = sm->GetId();
1785 // find _LayerEdge's of a sub-shape
1787 if ( data.GetShapeEdges( subID, edgesEnd, &iBeg, &iEnd ))
1788 cnvFace._subIdToEdgeEnd.insert( make_pair( subID, edgesEnd ));
1791 // check concavity and curvature and limit data._stepSize
1792 int nbLEdges = iEnd - iBeg;
1793 int step = Max( 1, nbLEdges / nbTestPnt );
1794 for ( ; iBeg < iEnd; iBeg += step )
1796 gp_XY uv = helper.GetNodeUV( F, data._edges[ iBeg ]->_nodes[0] );
1797 surfProp.SetParameters( uv.X(), uv.Y() );
1798 if ( !surfProp.IsCurvatureDefined() )
1800 if ( surfProp.MaxCurvature() * oriFactor > minCurvature )
1802 limitStepSize( data, 0.9 / surfProp.MaxCurvature() * oriFactor );
1805 if ( surfProp.MinCurvature() * oriFactor > minCurvature )
1807 limitStepSize( data, 0.9 / surfProp.MinCurvature() * oriFactor );
1811 } // loop on sub-shapes of the FACE
1813 if ( !isTooCurved ) continue;
1815 _ConvexFace & convFace =
1816 data._convexFaces.insert( make_pair( faceID, cnvFace )).first->second;
1819 convFace._normalsFixed = false;
1821 // Fill _ConvexFace::_simplexTestEdges. These _LayerEdge's are used to detect
1822 // prism distortion.
1823 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.find( faceID );
1824 if ( id2end != convFace._subIdToEdgeEnd.end() )
1826 // there are _LayerEdge's on the FACE it-self;
1827 // select _LayerEdge's near EDGEs
1828 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
1829 for ( ; iBeg < iEnd; ++iBeg )
1831 _LayerEdge* ledge = data._edges[ iBeg ];
1832 for ( size_t j = 0; j < ledge->_simplices.size(); ++j )
1833 if ( ledge->_simplices[j]._nNext->GetPosition()->GetDim() < 2 )
1835 convFace._simplexTestEdges.push_back( ledge );
1842 // where there are no _LayerEdge's on a _ConvexFace,
1843 // as e.g. on a fillet surface with no internal nodes - issue 22580,
1844 // so that collision of viscous internal faces is not detected by check of
1845 // intersection of _LayerEdge's with the viscous internal faces.
1847 set< const SMDS_MeshNode* > usedNodes;
1849 // look for _LayerEdge's with null _sWOL
1850 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.begin();
1851 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
1853 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
1854 if ( iBeg >= iEnd || !data._edges[ iBeg ]->_sWOL.IsNull() )
1856 for ( ; iBeg < iEnd; ++iBeg )
1858 _LayerEdge* ledge = data._edges[ iBeg ];
1859 const SMDS_MeshNode* srcNode = ledge->_nodes[0];
1860 if ( !usedNodes.insert( srcNode ).second ) continue;
1862 getSimplices( srcNode, ledge->_simplices, data._ignoreFaceIds, &data );
1863 for ( size_t i = 0; i < ledge->_simplices.size(); ++i )
1865 usedNodes.insert( ledge->_simplices[i]._nPrev );
1866 usedNodes.insert( ledge->_simplices[i]._nNext );
1868 convFace._simplexTestEdges.push_back( ledge );
1872 } // loop on FACEs of data._solid
1875 //================================================================================
1877 * \brief Separate shapes (and _LayerEdge's on them) to smooth from the rest ones
1879 //================================================================================
1881 bool _ViscousBuilder::sortEdges( _SolidData& data,
1882 vector< vector<_LayerEdge*> >& edgesByGeom)
1884 // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's
1885 // boundry inclined at a sharp angle to the shape
1887 list< TGeomID > shapesToSmooth;
1889 SMESH_MesherHelper helper( *_mesh );
1892 for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
1894 vector<_LayerEdge*>& eS = edgesByGeom[iS];
1895 if ( eS.empty() ) continue;
1896 const TopoDS_Shape& S = getMeshDS()->IndexToShape( iS );
1897 bool needSmooth = false;
1898 switch ( S.ShapeType() )
1902 bool isShrinkEdge = !eS[0]->_sWOL.IsNull();
1903 for ( TopoDS_Iterator vIt( S ); vIt.More() && !needSmooth; vIt.Next() )
1905 TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() );
1906 vector<_LayerEdge*>& eV = edgesByGeom[ iV ];
1907 if ( eV.empty() ) continue;
1908 // double cosin = eV[0]->_cosin;
1910 // ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge));
1913 // gp_Vec dir1, dir2;
1914 // if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE )
1915 // dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ));
1917 // dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ),
1918 // eV[0]->_nodes[0], helper, ok);
1919 // dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
1920 // double angle = dir1.Angle( dir2 );
1921 // cosin = cos( angle );
1923 gp_Vec eDir = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() ));
1924 double angle = eDir.Angle( eV[0]->_normal );
1925 double cosin = Cos( angle );
1926 needSmooth = ( cosin > theMinSmoothCosin );
1932 for ( TopExp_Explorer eExp( S, TopAbs_EDGE ); eExp.More() && !needSmooth; eExp.Next() )
1934 TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() );
1935 vector<_LayerEdge*>& eE = edgesByGeom[ iE ];
1936 if ( eE.empty() ) continue;
1937 if ( eE[0]->_sWOL.IsNull() )
1939 for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
1940 needSmooth = ( eE[i]->_cosin > theMinSmoothCosin );
1944 const TopoDS_Face& F1 = TopoDS::Face( S );
1945 const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL );
1946 const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() );
1947 for ( size_t i = 0; i < eE.size() && !needSmooth; ++i )
1949 gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok );
1950 gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok );
1951 double angle = dir1.Angle( dir2 );
1952 double cosin = cos( angle );
1953 needSmooth = ( cosin > theMinSmoothCosin );
1965 if ( S.ShapeType() == TopAbs_EDGE ) shapesToSmooth.push_front( iS );
1966 else shapesToSmooth.push_back ( iS );
1969 } // loop on edgesByGeom
1971 data._edges.reserve( data._n2eMap.size() );
1972 data._endEdgeOnShape.clear();
1974 // first we put _LayerEdge's on shapes to smooth
1975 data._nbShapesToSmooth = 0;
1976 list< TGeomID >::iterator gIt = shapesToSmooth.begin();
1977 for ( ; gIt != shapesToSmooth.end(); ++gIt )
1979 vector<_LayerEdge*>& eVec = edgesByGeom[ *gIt ];
1980 if ( eVec.empty() ) continue;
1981 data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
1982 data._endEdgeOnShape.push_back( data._edges.size() );
1983 data._nbShapesToSmooth++;
1987 // then the rest _LayerEdge's
1988 for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS )
1990 vector<_LayerEdge*>& eVec = edgesByGeom[iS];
1991 if ( eVec.empty() ) continue;
1992 data._edges.insert( data._edges.end(), eVec.begin(), eVec.end() );
1993 data._endEdgeOnShape.push_back( data._edges.size() );
2000 //================================================================================
2002 * \brief Set data of _LayerEdge needed for smoothing
2003 * \param subIds - ids of sub-shapes of a SOLID to take into account faces from
2005 //================================================================================
2007 bool _ViscousBuilder::setEdgeData(_LayerEdge& edge,
2008 const set<TGeomID>& subIds,
2009 SMESH_MesherHelper& helper,
2012 SMESH_MeshEditor editor(_mesh);
2014 const SMDS_MeshNode* node = edge._nodes[0]; // source node
2015 SMDS_TypeOfPosition posType = node->GetPosition()->GetTypeOfPosition();
2019 edge._curvature = 0;
2021 // --------------------------
2022 // Compute _normal and _cosin
2023 // --------------------------
2026 edge._normal.SetCoord(0,0,0);
2028 int totalNbFaces = 0;
2032 const TGeomID shapeInd = node->getshapeId();
2033 map< TGeomID, TopoDS_Shape >::const_iterator s2s = data._shrinkShape2Shape.find( shapeInd );
2034 const bool onShrinkShape ( s2s != data._shrinkShape2Shape.end() );
2036 if ( onShrinkShape ) // one of faces the node is on has no layers
2038 TopoDS_Shape vertEdge = getMeshDS()->IndexToShape( s2s->first ); // vertex or edge
2039 if ( s2s->second.ShapeType() == TopAbs_EDGE )
2041 // inflate from VERTEX along EDGE
2042 edge._normal = getEdgeDir( TopoDS::Edge( s2s->second ), TopoDS::Vertex( vertEdge ));
2044 else if ( vertEdge.ShapeType() == TopAbs_VERTEX )
2046 // inflate from VERTEX along FACE
2047 edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Vertex( vertEdge ),
2048 node, helper, normOK, &edge._cosin);
2052 // inflate from EDGE along FACE
2053 edge._normal = getFaceDir( TopoDS::Face( s2s->second ), TopoDS::Edge( vertEdge ),
2054 node, helper, normOK);
2057 else // layers are on all faces of SOLID the node is on
2059 // find indices of geom faces the node lies on
2060 set<TGeomID> faceIds;
2061 if ( posType == SMDS_TOP_FACE )
2063 faceIds.insert( node->getshapeId() );
2067 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2068 while ( fIt->more() )
2069 faceIds.insert( editor.FindShape(fIt->next()));
2072 set<TGeomID>::iterator id = faceIds.begin();
2074 std::pair< TGeomID, gp_XYZ > id2Norm[20];
2075 for ( ; id != faceIds.end(); ++id )
2077 const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id );
2078 if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id ))
2080 F = TopoDS::Face( s );
2081 geomNorm = getFaceNormal( node, F, helper, normOK );
2082 if ( !normOK ) continue;
2084 if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
2086 id2Norm[ totalNbFaces ].first = *id;
2087 id2Norm[ totalNbFaces ].second = geomNorm.XYZ();
2089 edge._normal += geomNorm.XYZ();
2091 if ( totalNbFaces == 0 )
2092 return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index);
2094 if ( normOK && edge._normal.Modulus() < 1e-3 && totalNbFaces > 1 )
2096 // opposite normals, re-get normals at shifted positions (IPAL 52426)
2097 edge._normal.SetCoord( 0,0,0 );
2098 for ( int i = 0; i < totalNbFaces; ++i )
2100 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( id2Norm[i].first ));
2101 geomNorm = getFaceNormal( node, F, helper, normOK, /*shiftInside=*/true );
2102 if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED )
2105 id2Norm[ i ].second = geomNorm.XYZ();
2106 edge._normal += id2Norm[ i ].second;
2110 if ( totalNbFaces < 3 )
2112 //edge._normal /= totalNbFaces;
2116 edge._normal = getWeigthedNormal( node, id2Norm, totalNbFaces );
2122 edge._cosin = 0; break;
2124 case SMDS_TOP_EDGE: {
2125 TopoDS_Edge E = TopoDS::Edge( helper.GetSubShapeByNode( node, getMeshDS()));
2126 gp_Vec inFaceDir = getFaceDir( F, E, node, helper, normOK );
2127 double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
2128 edge._cosin = cos( angle );
2129 //cout << "Cosin on EDGE " << edge._cosin << " node " << node->GetID() << endl;
2132 case SMDS_TOP_VERTEX: {
2133 TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS()));
2134 gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK );
2135 double angle = inFaceDir.Angle( edge._normal ); // [0,PI]
2136 edge._cosin = cos( angle );
2137 //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl;
2141 return error(SMESH_Comment("Invalid shape position of node ")<<node, data._index);
2143 } // case _sWOL.IsNull()
2145 double normSize = edge._normal.SquareModulus();
2146 if ( normSize < numeric_limits<double>::min() )
2147 return error(SMESH_Comment("Bad normal at node ")<< node->GetID(), data._index );
2149 edge._normal /= sqrt( normSize );
2151 // TODO: if ( !normOK ) then get normal by mesh faces
2153 // Set the rest data
2154 // --------------------
2155 if ( onShrinkShape )
2157 edge._sWOL = (*s2s).second;
2159 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( edge._nodes.back() );
2160 if ( SMESHDS_SubMesh* sm = getMeshDS()->MeshElements( data._solid ))
2161 sm->RemoveNode( tgtNode , /*isNodeDeleted=*/false );
2163 // set initial position which is parameters on _sWOL in this case
2164 if ( edge._sWOL.ShapeType() == TopAbs_EDGE )
2166 double u = helper.GetNodeU( TopoDS::Edge( edge._sWOL ), node, 0, &normOK );
2167 edge._pos.push_back( gp_XYZ( u, 0, 0));
2168 getMeshDS()->SetNodeOnEdge( tgtNode, TopoDS::Edge( edge._sWOL ), u );
2172 gp_XY uv = helper.GetNodeUV( TopoDS::Face( edge._sWOL ), node, 0, &normOK );
2173 edge._pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
2174 getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( edge._sWOL ), uv.X(), uv.Y() );
2179 edge._pos.push_back( SMESH_TNodeXYZ( node ));
2181 if ( posType == SMDS_TOP_FACE )
2183 getSimplices( node, edge._simplices, data._ignoreFaceIds, &data );
2184 double avgNormProj = 0, avgLen = 0;
2185 for ( size_t i = 0; i < edge._simplices.size(); ++i )
2187 gp_XYZ vec = edge._pos.back() - SMESH_TNodeXYZ( edge._simplices[i]._nPrev );
2188 avgNormProj += edge._normal * vec;
2189 avgLen += vec.Modulus();
2191 avgNormProj /= edge._simplices.size();
2192 avgLen /= edge._simplices.size();
2193 edge._curvature = _Curvature::New( avgNormProj, avgLen );
2197 // Set neighbour nodes for a _LayerEdge based on EDGE
2199 if ( posType == SMDS_TOP_EDGE /*||
2200 ( onShrinkShape && posType == SMDS_TOP_VERTEX && fabs( edge._cosin ) < 1e-10 )*/)
2202 edge._2neibors = new _2NearEdges;
2203 // target node instead of source ones will be set later
2204 if ( ! findNeiborsOnEdge( &edge,
2205 edge._2neibors->_nodes[0],
2206 edge._2neibors->_nodes[1],
2209 edge.SetDataByNeighbors( edge._2neibors->_nodes[0],
2210 edge._2neibors->_nodes[1],
2214 edge.SetCosin( edge._cosin ); // to update edge._lenFactor
2219 //================================================================================
2221 * \brief Return normal to a FACE at a node
2222 * \param [in] n - node
2223 * \param [in] face - FACE
2224 * \param [in] helper - helper
2225 * \param [out] isOK - true or false
2226 * \param [in] shiftInside - to find normal at a position shifted inside the face
2227 * \return gp_XYZ - normal
2229 //================================================================================
2231 gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node,
2232 const TopoDS_Face& face,
2233 SMESH_MesherHelper& helper,
2240 // get a shifted position
2241 gp_Pnt p = SMESH_TNodeXYZ( node );
2242 gp_XYZ shift( 0,0,0 );
2243 TopoDS_Shape S = helper.GetSubShapeByNode( node, helper.GetMeshDS() );
2244 switch ( S.ShapeType() ) {
2247 shift = getFaceDir( face, TopoDS::Vertex( S ), node, helper, isOK );
2252 shift = getFaceDir( face, TopoDS::Edge( S ), node, helper, isOK );
2260 p.Translate( shift * 1e-5 );
2262 TopLoc_Location loc;
2263 GeomAPI_ProjectPointOnSurf& projector = helper.GetProjector( face, loc, 1e-7 );
2265 if ( !loc.IsIdentity() ) p.Transform( loc.Transformation().Inverted() );
2267 projector.Perform( p );
2268 if ( !projector.IsDone() || projector.NbPoints() < 1 )
2273 Quantity_Parameter U,V;
2274 projector.LowerDistanceParameters(U,V);
2279 uv = helper.GetNodeUV( face, node, 0, &isOK );
2285 Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
2286 if ( GeomLib::NormEstim( surface, uv, 1e-10, normal ) < 3 )
2291 else // hard singularity
2293 const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face );
2295 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2296 while ( fIt->more() )
2298 const SMDS_MeshElement* f = fIt->next();
2299 if ( f->getshapeId() == faceID )
2301 isOK = SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) normal.XYZ(), /*normalized=*/true );
2304 if ( helper.IsReversedSubMesh( face ))
2311 return normal.XYZ();
2314 //================================================================================
2316 * \brief Return a normal at a node weighted with angles taken by FACEs
2317 * \param [in] n - the node
2318 * \param [in] fId2Normal - FACE ids and normals
2319 * \param [in] nbFaces - nb of FACEs meeting at the node
2320 * \return gp_XYZ - computed normal
2322 //================================================================================
2324 gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode* n,
2325 std::pair< TGeomID, gp_XYZ > fId2Normal[],
2328 gp_XYZ resNorm(0,0,0);
2329 TopoDS_Shape V = SMESH_MesherHelper::GetSubShapeByNode( n, getMeshDS() );
2330 if ( V.ShapeType() != TopAbs_VERTEX )
2332 for ( int i = 0; i < nbFaces; ++i )
2333 resNorm += fId2Normal[i].second / nbFaces ;
2338 for ( int i = 0; i < nbFaces; ++i )
2340 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( fId2Normal[i].first ));
2342 // look for two EDGEs shared by F and other FACEs within fId2Normal
2345 PShapeIteratorPtr eIt = SMESH_MesherHelper::GetAncestors( V, *_mesh, TopAbs_EDGE );
2346 while ( const TopoDS_Shape* E = eIt->next() )
2348 if ( !SMESH_MesherHelper::IsSubShape( *E, F ))
2350 bool isSharedEdge = false;
2351 for ( int j = 0; j < nbFaces && !isSharedEdge; ++j )
2353 if ( i == j ) continue;
2354 const TopoDS_Shape& otherF = getMeshDS()->IndexToShape( fId2Normal[j].first );
2355 isSharedEdge = SMESH_MesherHelper::IsSubShape( *E, otherF );
2357 if ( !isSharedEdge )
2359 ee[ nbE ] = TopoDS::Edge( *E );
2360 ee[ nbE ].Orientation( SMESH_MesherHelper::GetSubShapeOri( F, *E ));
2365 // get an angle between the two EDGEs
2367 if ( nbE < 1 ) continue;
2374 TopoDS_Vertex v10 = SMESH_MesherHelper::IthVertex( 1, ee[ 0 ]);
2375 TopoDS_Vertex v01 = SMESH_MesherHelper::IthVertex( 0, ee[ 1 ]);
2376 if ( !v10.IsSame( v01 ))
2377 std::swap( ee[0], ee[1] );
2379 angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F );
2382 // compute a weighted normal
2383 double sumAngle = 0;
2384 for ( int i = 0; i < nbFaces; ++i )
2386 angles[i] = ( angles[i] > 2*M_PI ) ? 0 : M_PI - angles[i];
2387 sumAngle += angles[i];
2389 for ( int i = 0; i < nbFaces; ++i )
2390 resNorm += angles[i] / sumAngle * fId2Normal[i].second;
2395 //================================================================================
2397 * \brief Find 2 neigbor nodes of a node on EDGE
2399 //================================================================================
2401 bool _ViscousBuilder::findNeiborsOnEdge(const _LayerEdge* edge,
2402 const SMDS_MeshNode*& n1,
2403 const SMDS_MeshNode*& n2,
2406 const SMDS_MeshNode* node = edge->_nodes[0];
2407 const int shapeInd = node->getshapeId();
2408 SMESHDS_SubMesh* edgeSM = 0;
2409 if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE )
2412 edgeSM = getMeshDS()->MeshElements( shapeInd );
2413 if ( !edgeSM || edgeSM->NbElements() == 0 )
2414 return error(SMESH_Comment("Not meshed EDGE ") << shapeInd, data._index);
2418 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(SMDSAbs_Edge);
2419 while ( eIt->more() && !n2 )
2421 const SMDS_MeshElement* e = eIt->next();
2422 const SMDS_MeshNode* nNeibor = e->GetNode( 0 );
2423 if ( nNeibor == node ) nNeibor = e->GetNode( 1 );
2426 if (!edgeSM->Contains(e)) continue;
2430 TopoDS_Shape s = SMESH_MesherHelper::GetSubShapeByNode(nNeibor, getMeshDS() );
2431 if ( !SMESH_MesherHelper::IsSubShape( s, edge->_sWOL )) continue;
2433 ( iN++ ? n2 : n1 ) = nNeibor;
2436 return error(SMESH_Comment("Wrongly meshed EDGE ") << shapeInd, data._index);
2440 //================================================================================
2442 * \brief Set _curvature and _2neibors->_plnNorm by 2 neigbor nodes residing the same EDGE
2444 //================================================================================
2446 void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1,
2447 const SMDS_MeshNode* n2,
2448 SMESH_MesherHelper& helper)
2450 if ( _nodes[0]->GetPosition()->GetTypeOfPosition() != SMDS_TOP_EDGE )
2453 gp_XYZ pos = SMESH_TNodeXYZ( _nodes[0] );
2454 gp_XYZ vec1 = pos - SMESH_TNodeXYZ( n1 );
2455 gp_XYZ vec2 = pos - SMESH_TNodeXYZ( n2 );
2459 double sumLen = vec1.Modulus() + vec2.Modulus();
2460 _2neibors->_wgt[0] = 1 - vec1.Modulus() / sumLen;
2461 _2neibors->_wgt[1] = 1 - vec2.Modulus() / sumLen;
2462 double avgNormProj = 0.5 * ( _normal * vec1 + _normal * vec2 );
2463 double avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() );
2464 if ( _curvature ) delete _curvature;
2465 _curvature = _Curvature::New( avgNormProj, avgLen );
2466 // if ( _curvature )
2467 // debugMsg( _nodes[0]->GetID()
2468 // << " CURV r,k: " << _curvature->_r<<","<<_curvature->_k
2469 // << " proj = "<<avgNormProj<< " len = " << avgLen << "| lenDelta(0) = "
2470 // << _curvature->lenDelta(0) );
2474 if ( _sWOL.IsNull() )
2476 TopoDS_Shape S = helper.GetSubShapeByNode( _nodes[0], helper.GetMeshDS() );
2477 gp_XYZ dirE = getEdgeDir( TopoDS::Edge( S ), _nodes[0], helper );
2478 gp_XYZ plnNorm = dirE ^ _normal;
2479 double proj0 = plnNorm * vec1;
2480 double proj1 = plnNorm * vec2;
2481 if ( fabs( proj0 ) > 1e-10 || fabs( proj1 ) > 1e-10 )
2483 if ( _2neibors->_plnNorm ) delete _2neibors->_plnNorm;
2484 _2neibors->_plnNorm = new gp_XYZ( plnNorm.Normalized() );
2489 //================================================================================
2491 * \brief Copy data from a _LayerEdge of other SOLID and based on the same node;
2492 * this and other _LayerEdge's are inflated along a FACE or an EDGE
2494 //================================================================================
2496 gp_XYZ _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper )
2498 _nodes = other._nodes;
2499 _normal = other._normal;
2501 _lenFactor = other._lenFactor;
2502 _cosin = other._cosin;
2503 _sWOL = other._sWOL;
2504 _2neibors = other._2neibors;
2505 _curvature = 0; std::swap( _curvature, other._curvature );
2506 _2neibors = 0; std::swap( _2neibors, other._2neibors );
2508 gp_XYZ lastPos( 0,0,0 );
2509 if ( _sWOL.ShapeType() == TopAbs_EDGE )
2511 double u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes[0] );
2512 _pos.push_back( gp_XYZ( u, 0, 0));
2514 u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes.back() );
2519 gp_XY uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes[0]);
2520 _pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0));
2522 uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes.back() );
2523 lastPos.SetX( uv.X() );
2524 lastPos.SetY( uv.Y() );
2529 //================================================================================
2531 * \brief Set _cosin and _lenFactor
2533 //================================================================================
2535 void _LayerEdge::SetCosin( double cosin )
2538 cosin = Abs( _cosin );
2539 _lenFactor = ( /*0.1 < cosin &&*/ cosin < 1-1e-12 ) ? 1./sqrt(1-cosin*cosin) : 1.0;
2542 //================================================================================
2544 * \brief Fills a vector<_Simplex >
2546 //================================================================================
2548 void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node,
2549 vector<_Simplex>& simplices,
2550 const set<TGeomID>& ingnoreShapes,
2551 const _SolidData* dataToCheckOri,
2555 SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face);
2556 while ( fIt->more() )
2558 const SMDS_MeshElement* f = fIt->next();
2559 const TGeomID shapeInd = f->getshapeId();
2560 if ( ingnoreShapes.count( shapeInd )) continue;
2561 const int nbNodes = f->NbCornerNodes();
2562 const int srcInd = f->GetNodeIndex( node );
2563 const SMDS_MeshNode* nPrev = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd-1, nbNodes ));
2564 const SMDS_MeshNode* nNext = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+1, nbNodes ));
2565 const SMDS_MeshNode* nOpp = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+2, nbNodes ));
2566 if ( dataToCheckOri && dataToCheckOri->_reversedFaceIds.count( shapeInd ))
2567 std::swap( nPrev, nNext );
2568 simplices.push_back( _Simplex( nPrev, nNext, nOpp ));
2573 vector<_Simplex> sortedSimplices( simplices.size() );
2574 sortedSimplices[0] = simplices[0];
2576 for ( size_t i = 1; i < simplices.size(); ++i )
2578 for ( size_t j = 1; j < simplices.size(); ++j )
2579 if ( sortedSimplices[i-1]._nNext == simplices[j]._nPrev )
2581 sortedSimplices[i] = simplices[j];
2586 if ( nbFound == simplices.size() - 1 )
2587 simplices.swap( sortedSimplices );
2591 //================================================================================
2593 * \brief DEBUG. Create groups contating temorary data of _LayerEdge's
2595 //================================================================================
2597 void _ViscousBuilder::makeGroupOfLE()
2600 for ( size_t i = 0 ; i < _sdVec.size(); ++i )
2602 if ( _sdVec[i]._edges.empty() ) continue;
2604 dumpFunction( SMESH_Comment("make_LayerEdge_") << i );
2605 for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
2607 _LayerEdge* le = _sdVec[i]._edges[j];
2608 for ( size_t iN = 1; iN < le->_nodes.size(); ++iN )
2609 dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<le->_nodes[iN-1]->GetID()
2610 << ", " << le->_nodes[iN]->GetID() <<"])");
2614 dumpFunction( SMESH_Comment("makeNormals") << i );
2615 for ( size_t j = 0 ; j < _sdVec[i]._edges.size(); ++j )
2617 _LayerEdge& edge = *_sdVec[i]._edges[j];
2618 SMESH_TNodeXYZ nXYZ( edge._nodes[0] );
2619 nXYZ += edge._normal * _sdVec[i]._stepSize;
2620 dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <<edge._nodes[0]->GetID()
2621 << ", mesh.AddNode( " << nXYZ.X()<<","<< nXYZ.Y()<<","<< nXYZ.Z()<<")])");
2625 dumpFunction( SMESH_Comment("makeTmpFaces_") << i );
2626 TopExp_Explorer fExp( _sdVec[i]._solid, TopAbs_FACE );
2627 for ( ; fExp.More(); fExp.Next() )
2629 if (const SMESHDS_SubMesh* sm = _sdVec[i]._proxyMesh->GetProxySubMesh( fExp.Current()))
2631 SMDS_ElemIteratorPtr fIt = sm->GetElements();
2632 while ( fIt->more())
2634 const SMDS_MeshElement* e = fIt->next();
2635 SMESH_Comment cmd("mesh.AddFace([");
2636 for ( int j=0; j < e->NbCornerNodes(); ++j )
2637 cmd << e->GetNode(j)->GetID() << (j+1<e->NbCornerNodes() ? ",": "])");
2647 //================================================================================
2649 * \brief Increase length of _LayerEdge's to reach the required thickness of layers
2651 //================================================================================
2653 bool _ViscousBuilder::inflate(_SolidData& data)
2655 SMESH_MesherHelper helper( *_mesh );
2657 // Limit inflation step size by geometry size found by itersecting
2658 // normals of _LayerEdge's with mesh faces
2659 double geomSize = Precision::Infinite(), intersecDist;
2660 auto_ptr<SMESH_ElementSearcher> searcher
2661 ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
2662 data._proxyMesh->GetFaces( data._solid )) );
2663 for ( size_t i = 0; i < data._edges.size(); ++i )
2665 if ( data._edges[i]->IsOnEdge() ) continue;
2666 data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon );
2667 if ( geomSize > intersecDist && intersecDist > 0 )
2668 geomSize = intersecDist;
2670 if ( data._stepSize > 0.3 * geomSize )
2671 limitStepSize( data, 0.3 * geomSize );
2673 const double tgtThick = data._hyp->GetTotalThickness();
2674 if ( data._stepSize > tgtThick )
2675 limitStepSize( data, tgtThick );
2677 if ( data._stepSize < 1. )
2678 data._epsilon = data._stepSize * 1e-7;
2680 debugMsg( "-- geomSize = " << geomSize << ", stepSize = " << data._stepSize );
2682 double avgThick = 0, curThick = 0, distToIntersection = Precision::Infinite();
2683 int nbSteps = 0, nbRepeats = 0;
2684 while ( 1.01 * avgThick < tgtThick )
2686 // new target length
2687 curThick += data._stepSize;
2688 if ( curThick > tgtThick )
2690 curThick = tgtThick + ( tgtThick-avgThick ) * nbRepeats;
2694 // Elongate _LayerEdge's
2695 dumpFunction(SMESH_Comment("inflate")<<data._index<<"_step"<<nbSteps); // debug
2696 for ( size_t i = 0; i < data._edges.size(); ++i )
2698 data._edges[i]->SetNewLength( curThick, helper );
2702 if ( !updateNormals( data, helper, nbSteps ))
2705 // Improve and check quality
2706 if ( !smoothAndCheck( data, nbSteps, distToIntersection ))
2710 dumpFunction(SMESH_Comment("invalidate")<<data._index<<"_step"<<nbSteps); // debug
2711 for ( size_t i = 0; i < data._edges.size(); ++i )
2713 data._edges[i]->InvalidateStep( nbSteps+1 );
2717 break; // no more inflating possible
2721 // Evaluate achieved thickness
2723 for ( size_t i = 0; i < data._edges.size(); ++i )
2724 avgThick += data._edges[i]->_len;
2725 avgThick /= data._edges.size();
2726 debugMsg( "-- Thickness " << avgThick << " reached" );
2728 if ( distToIntersection < avgThick*1.5 )
2730 debugMsg( "-- Stop inflation since "
2731 << " distToIntersection( "<<distToIntersection<<" ) < avgThick( "
2732 << avgThick << " ) * 1.5" );
2736 limitStepSize( data, 0.25 * distToIntersection );
2737 if ( data._stepSizeNodes[0] )
2738 data._stepSize = data._stepSizeCoeff *
2739 SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]);
2741 } // while ( 1.01 * avgThick < tgtThick )
2744 return error("failed at the very first inflation step", data._index);
2746 if ( 1.01 * avgThick < tgtThick )
2747 if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( data._index ))
2749 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
2750 if ( !smError || smError->IsOK() )
2752 ( new SMESH_ComputeError (COMPERR_WARNING,
2753 SMESH_Comment("Thickness ") << tgtThick <<
2754 " of viscous layers not reached,"
2755 " average reached thickness is " << avgThick ));
2761 //================================================================================
2763 * \brief Improve quality of layer inner surface and check intersection
2765 //================================================================================
2767 bool _ViscousBuilder::smoothAndCheck(_SolidData& data,
2769 double & distToIntersection)
2771 if ( data._nbShapesToSmooth == 0 )
2772 return true; // no shapes needing smoothing
2774 bool moved, improved;
2776 SMESH_MesherHelper helper(*_mesh);
2777 Handle(Geom_Surface) surface;
2781 for ( int iS = 0; iS < data._nbShapesToSmooth; ++iS )
2784 iEnd = data._endEdgeOnShape[ iS ];
2786 if ( !data._edges[ iBeg ]->_sWOL.IsNull() &&
2787 data._edges[ iBeg ]->_sWOL.ShapeType() == TopAbs_FACE )
2789 if ( !F.IsSame( data._edges[ iBeg ]->_sWOL )) {
2790 F = TopoDS::Face( data._edges[ iBeg ]->_sWOL );
2791 helper.SetSubShape( F );
2792 surface = BRep_Tool::Surface( F );
2797 F.Nullify(); surface.Nullify();
2799 TGeomID sInd = data._edges[ iBeg ]->_nodes[0]->getshapeId();
2801 if ( data._edges[ iBeg ]->IsOnEdge() )
2803 dumpFunction(SMESH_Comment("smooth")<<data._index << "_Ed"<<sInd <<"_InfStep"<<nbSteps);
2805 // try a simple solution on an analytic EDGE
2806 if ( !smoothAnalyticEdge( data, iBeg, iEnd, surface, F, helper ))
2812 for ( int i = iBeg; i < iEnd; ++i )
2814 moved |= data._edges[i]->SmoothOnEdge(surface, F, helper);
2816 dumpCmd( SMESH_Comment("# end step ")<<step);
2818 while ( moved && step++ < 5 );
2825 int step = 0, stepLimit = 5, badNb = 0; moved = true;
2826 while (( ++step <= stepLimit && moved ) || improved )
2828 dumpFunction(SMESH_Comment("smooth")<<data._index<<"_Fa"<<sInd
2829 <<"_InfStep"<<nbSteps<<"_"<<step); // debug
2830 int oldBadNb = badNb;
2834 for ( int i = iBeg; i < iEnd; ++i )
2835 moved |= data._edges[i]->Smooth(badNb);
2837 for ( int i = iEnd-1; i >= iBeg; --i )
2838 moved |= data._edges[i]->Smooth(badNb);
2839 improved = ( badNb < oldBadNb );
2841 // issue 22576 -- no bad faces but still there are intersections to fix
2842 if ( improved && badNb == 0 )
2843 stepLimit = step + 3;
2850 for ( int i = iBeg; i < iEnd; ++i )
2852 _LayerEdge* edge = data._edges[i];
2853 SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
2854 for ( size_t j = 0; j < edge->_simplices.size(); ++j )
2855 if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
2857 cout << "Bad simplex ( " << edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
2858 << " "<< edge->_simplices[j]._nPrev->GetID()
2859 << " "<< edge->_simplices[j]._nNext->GetID() << " )" << endl;
2867 } // loop on shapes to smooth
2869 // Check orientation of simplices of _ConvexFace::_simplexTestEdges
2870 map< TGeomID, _ConvexFace >::iterator id2face = data._convexFaces.begin();
2871 for ( ; id2face != data._convexFaces.end(); ++id2face )
2873 _ConvexFace & convFace = (*id2face).second;
2874 if ( !convFace._simplexTestEdges.empty() &&
2875 convFace._simplexTestEdges[0]->_nodes[0]->GetPosition()->GetDim() == 2 )
2876 continue; // _simplexTestEdges are based on FACE -- already checked while smoothing
2878 if ( !convFace.CheckPrisms() )
2882 // Check if the last segments of _LayerEdge intersects 2D elements;
2883 // checked elements are either temporary faces or faces on surfaces w/o the layers
2885 auto_ptr<SMESH_ElementSearcher> searcher
2886 ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(),
2887 data._proxyMesh->GetFaces( data._solid )) );
2889 distToIntersection = Precision::Infinite();
2891 const SMDS_MeshElement* intFace = 0;
2892 const SMDS_MeshElement* closestFace = 0;
2894 for ( size_t i = 0; i < data._edges.size(); ++i )
2896 if ( !data._edges[i]->_sWOL.IsNull() )
2898 if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace ))
2900 if ( distToIntersection > dist )
2902 // ignore intersection of a _LayerEdge based on a _ConvexFace with a face
2903 // lying on this _ConvexFace
2904 if ( _ConvexFace* convFace = data.GetConvexFace( intFace->getshapeId() ))
2905 if ( convFace->_subIdToEdgeEnd.count ( data._edges[i]->_nodes[0]->getshapeId() ))
2908 distToIntersection = dist;
2910 closestFace = intFace;
2916 SMDS_MeshElement::iterator nIt = closestFace->begin_nodes();
2917 cout << "Shortest distance: _LayerEdge nodes: tgt " << data._edges[iLE]->_nodes.back()->GetID()
2918 << " src " << data._edges[iLE]->_nodes[0]->GetID()<< ", intersection with face ("
2919 << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
2920 << ") distance = " << distToIntersection<< endl;
2927 //================================================================================
2929 * \brief Return a curve of the EDGE to be used for smoothing and arrange
2930 * _LayerEdge's to be in a consequent order
2932 //================================================================================
2934 Handle(Geom_Curve) _SolidData::CurveForSmooth( const TopoDS_Edge& E,
2937 Handle(Geom_Surface)& surface,
2938 const TopoDS_Face& F,
2939 SMESH_MesherHelper& helper)
2941 TGeomID eIndex = helper.GetMeshDS()->ShapeToIndex( E );
2943 map< TGeomID, Handle(Geom_Curve)>::iterator i2curve = _edge2curve.find( eIndex );
2945 if ( i2curve == _edge2curve.end() )
2947 // sort _LayerEdge's by position on the EDGE
2948 SortOnEdge( E, iFrom, iTo, helper );
2950 SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( eIndex );
2952 TopLoc_Location loc; double f,l;
2954 Handle(Geom_Line) line;
2955 Handle(Geom_Circle) circle;
2956 bool isLine, isCirc;
2957 if ( F.IsNull() ) // 3D case
2959 // check if the EDGE is a line
2960 Handle(Geom_Curve) curve = BRep_Tool::Curve( E, loc, f, l);
2961 if ( curve->IsKind( STANDARD_TYPE( Geom_TrimmedCurve )))
2962 curve = Handle(Geom_TrimmedCurve)::DownCast( curve )->BasisCurve();
2964 line = Handle(Geom_Line)::DownCast( curve );
2965 circle = Handle(Geom_Circle)::DownCast( curve );
2966 isLine = (!line.IsNull());
2967 isCirc = (!circle.IsNull());
2969 if ( !isLine && !isCirc ) // Check if the EDGE is close to a line
2972 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
2973 while ( nIt->more() )
2974 bndBox.Add( SMESH_TNodeXYZ( nIt->next() ));
2975 gp_XYZ size = bndBox.CornerMax() - bndBox.CornerMin();
2977 SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->_nodes[0] );
2978 SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->_nodes[1] );
2979 const double lineTol = 1e-2 * ( p0 - p1 ).Modulus();
2980 for ( int i = 0; i < 3 && !isLine; ++i )
2981 isLine = ( size.Coord( i+1 ) <= lineTol );
2983 if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle
2990 // check if the EDGE is a line
2991 Handle(Geom2d_Curve) curve = BRep_Tool::CurveOnSurface( E, F, f, l);
2992 if ( curve->IsKind( STANDARD_TYPE( Geom2d_TrimmedCurve )))
2993 curve = Handle(Geom2d_TrimmedCurve)::DownCast( curve )->BasisCurve();
2995 Handle(Geom2d_Line) line2d = Handle(Geom2d_Line)::DownCast( curve );
2996 Handle(Geom2d_Circle) circle2d = Handle(Geom2d_Circle)::DownCast( curve );
2997 isLine = (!line2d.IsNull());
2998 isCirc = (!circle2d.IsNull());
3000 if ( !isLine && !isCirc) // Check if the EDGE is close to a line
3003 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
3004 while ( nIt->more() )
3005 bndBox.Add( helper.GetNodeUV( F, nIt->next() ));
3006 gp_XY size = bndBox.CornerMax() - bndBox.CornerMin();
3008 const double lineTol = 1e-2 * sqrt( bndBox.SquareExtent() );
3009 for ( int i = 0; i < 2 && !isLine; ++i )
3010 isLine = ( size.Coord( i+1 ) <= lineTol );
3012 if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle
3018 line = new Geom_Line( gp::OX() ); // only type does matter
3022 gp_Pnt2d p = circle2d->Location();
3023 gp_Ax2 ax( gp_Pnt( p.X(), p.Y(), 0), gp::DX());
3024 circle = new Geom_Circle( ax, 1.); // only center position does matter
3028 Handle(Geom_Curve)& res = _edge2curve[ eIndex ];
3036 return i2curve->second;
3039 //================================================================================
3041 * \brief Sort _LayerEdge's by a parameter on a given EDGE
3043 //================================================================================
3045 void _SolidData::SortOnEdge( const TopoDS_Edge& E,
3048 SMESH_MesherHelper& helper)
3050 map< double, _LayerEdge* > u2edge;
3051 for ( int i = iFrom; i < iTo; ++i )
3052 u2edge.insert( make_pair( helper.GetNodeU( E, _edges[i]->_nodes[0] ), _edges[i] ));
3054 ASSERT( u2edge.size() == iTo - iFrom );
3055 map< double, _LayerEdge* >::iterator u2e = u2edge.begin();
3056 for ( int i = iFrom; i < iTo; ++i, ++u2e )
3057 _edges[i] = u2e->second;
3059 // set _2neibors according to the new order
3060 for ( int i = iFrom; i < iTo-1; ++i )
3061 if ( _edges[i]->_2neibors->_nodes[1] != _edges[i+1]->_nodes.back() )
3062 _edges[i]->_2neibors->reverse();
3063 if ( u2edge.size() > 1 &&
3064 _edges[iTo-1]->_2neibors->_nodes[0] != _edges[iTo-2]->_nodes.back() )
3065 _edges[iTo-1]->_2neibors->reverse();
3068 //================================================================================
3070 * \brief Return index corresponding to the shape in _endEdgeOnShape
3072 //================================================================================
3074 bool _SolidData::GetShapeEdges(const TGeomID shapeID,
3079 int beg = 0, end = 0;
3080 for ( edgesEnd = 0; edgesEnd < _endEdgeOnShape.size(); ++edgesEnd )
3082 end = _endEdgeOnShape[ edgesEnd ];
3083 TGeomID sID = _edges[ beg ]->_nodes[0]->getshapeId();
3084 if ( sID == shapeID )
3086 if ( iBeg ) *iBeg = beg;
3087 if ( iEnd ) *iEnd = end;
3095 //================================================================================
3097 * \brief Add faces for smoothing
3099 //================================================================================
3101 void _SolidData::AddShapesToSmooth( const set< TGeomID >& faceIDs )
3103 // convert faceIDs to indices in _endEdgeOnShape
3104 set< size_t > iEnds;
3106 set< TGeomID >::const_iterator fId = faceIDs.begin();
3107 for ( ; fId != faceIDs.end(); ++fId )
3108 if ( GetShapeEdges( *fId, end ) && end >= _nbShapesToSmooth )
3109 iEnds.insert( end );
3111 set< size_t >::iterator endsIt = iEnds.begin();
3113 // "add" by move of _nbShapesToSmooth
3114 int nbFacesToAdd = iEnds.size();
3115 while ( endsIt != iEnds.end() && *endsIt == _nbShapesToSmooth )
3118 ++_nbShapesToSmooth;
3121 if ( endsIt == iEnds.end() )
3124 // Move _LayerEdge's on FACEs just after _nbShapesToSmooth
3126 vector< _LayerEdge* > nonSmoothLE, smoothLE;
3127 size_t lastSmooth = *iEnds.rbegin();
3129 for ( size_t i = _nbShapesToSmooth; i <= lastSmooth; ++i )
3131 vector< _LayerEdge* > & edgesVec = iEnds.count(i) ? smoothLE : nonSmoothLE;
3132 iBeg = i ? _endEdgeOnShape[ i-1 ] : 0;
3133 iEnd = _endEdgeOnShape[ i ];
3134 edgesVec.insert( edgesVec.end(), _edges.begin() + iBeg, _edges.begin() + iEnd );
3137 iBeg = _nbShapesToSmooth ? _endEdgeOnShape[ _nbShapesToSmooth-1 ] : 0;
3138 std::copy( smoothLE.begin(), smoothLE.end(), &_edges[ iBeg ] );
3139 std::copy( nonSmoothLE.begin(), nonSmoothLE.end(), &_edges[ iBeg + smoothLE.size()]);
3141 // update _endEdgeOnShape
3142 for ( size_t i = _nbShapesToSmooth; i < _endEdgeOnShape.size(); ++i )
3144 TGeomID curShape = _edges[ iBeg ]->_nodes[0]->getshapeId();
3145 while ( ++iBeg < _edges.size() &&
3146 curShape == _edges[ iBeg ]->_nodes[0]->getshapeId() );
3148 _endEdgeOnShape[ i ] = iBeg;
3151 _nbShapesToSmooth += nbFacesToAdd;
3154 //================================================================================
3156 * \brief smooth _LayerEdge's on a staight EDGE or circular EDGE
3158 //================================================================================
3160 bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data,
3163 Handle(Geom_Surface)& surface,
3164 const TopoDS_Face& F,
3165 SMESH_MesherHelper& helper)
3167 TopoDS_Shape S = helper.GetSubShapeByNode( data._edges[ iFrom ]->_nodes[0],
3168 helper.GetMeshDS());
3169 TopoDS_Edge E = TopoDS::Edge( S );
3171 Handle(Geom_Curve) curve = data.CurveForSmooth( E, iFrom, iTo, surface, F, helper );
3172 if ( curve.IsNull() ) return false;
3174 // compute a relative length of segments
3175 vector< double > len( iTo-iFrom+1 );
3177 double curLen, prevLen = len[0] = 1.0;
3178 for ( int i = iFrom; i < iTo; ++i )
3180 curLen = prevLen * data._edges[i]->_2neibors->_wgt[0] / data._edges[i]->_2neibors->_wgt[1];
3181 len[i-iFrom+1] = len[i-iFrom] + curLen;
3186 if ( curve->IsKind( STANDARD_TYPE( Geom_Line )))
3188 if ( F.IsNull() ) // 3D
3190 SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->_nodes[0]);
3191 SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->_nodes[1]);
3192 for ( int i = iFrom; i < iTo; ++i )
3194 double r = len[i-iFrom] / len.back();
3195 gp_XYZ newPos = p0 * ( 1. - r ) + p1 * r;
3196 data._edges[i]->_pos.back() = newPos;
3197 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3198 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3199 dumpMove( tgtNode );
3204 gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]);
3205 gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]);
3206 if ( data._edges[iFrom]->_2neibors->_nodes[0] ==
3207 data._edges[iTo-1]->_2neibors->_nodes[1] ) // closed edge
3209 int iPeriodic = helper.GetPeriodicIndex();
3210 if ( iPeriodic == 1 || iPeriodic == 2 )
3212 uv1.SetCoord( iPeriodic, helper.GetOtherParam( uv1.Coord( iPeriodic )));
3213 if ( uv0.Coord( iPeriodic ) > uv1.Coord( iPeriodic ))
3214 std::swap( uv0, uv1 );
3217 const gp_XY rangeUV = uv1 - uv0;
3218 for ( int i = iFrom; i < iTo; ++i )
3220 double r = len[i-iFrom] / len.back();
3221 gp_XY newUV = uv0 + r * rangeUV;
3222 data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
3224 gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
3225 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3226 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3227 dumpMove( tgtNode );
3229 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
3230 pos->SetUParameter( newUV.X() );
3231 pos->SetVParameter( newUV.Y() );
3237 if ( curve->IsKind( STANDARD_TYPE( Geom_Circle )))
3239 Handle(Geom_Circle) circle = Handle(Geom_Circle)::DownCast( curve );
3240 gp_Pnt center3D = circle->Location();
3242 if ( F.IsNull() ) // 3D
3244 if ( data._edges[iFrom]->_2neibors->_nodes[0] ==
3245 data._edges[iTo-1]->_2neibors->_nodes[1] )
3246 return true; // closed EDGE - nothing to do
3248 return false; // TODO ???
3252 const gp_XY center( center3D.X(), center3D.Y() );
3254 gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]);
3255 gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back());
3256 gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]);
3257 gp_Vec2d vec0( center, uv0 );
3258 gp_Vec2d vecM( center, uvM );
3259 gp_Vec2d vec1( center, uv1 );
3260 double uLast = vec0.Angle( vec1 ); // -PI - +PI
3261 double uMidl = vec0.Angle( vecM );
3262 if ( uLast * uMidl <= 0. )
3263 uLast += ( uMidl > 0 ? +2. : -2. ) * M_PI;
3264 const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() );
3266 gp_Ax2d axis( center, vec0 );
3267 gp_Circ2d circ( axis, radius );
3268 for ( int i = iFrom; i < iTo; ++i )
3270 double newU = uLast * len[i-iFrom] / len.back();
3271 gp_Pnt2d newUV = ElCLib::Value( newU, circ );
3272 data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 );
3274 gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() );
3275 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( data._edges[i]->_nodes.back() );
3276 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
3277 dumpMove( tgtNode );
3279 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
3280 pos->SetUParameter( newUV.X() );
3281 pos->SetVParameter( newUV.Y() );
3290 //================================================================================
3292 * \brief Modify normals of _LayerEdge's on EDGE's to avoid intersection with
3293 * _LayerEdge's on neighbor EDGE's
3295 //================================================================================
3297 bool _ViscousBuilder::updateNormals( _SolidData& data,
3298 SMESH_MesherHelper& helper,
3302 return updateNormalsOfConvexFaces( data, helper, stepNb );
3304 // make temporary quadrangles got by extrusion of
3305 // mesh edges along _LayerEdge._normal's
3307 vector< const SMDS_MeshElement* > tmpFaces;
3309 set< SMESH_TLink > extrudedLinks; // contains target nodes
3310 vector< const SMDS_MeshNode*> nodes(4); // of a tmp mesh face
3312 dumpFunction(SMESH_Comment("makeTmpFacesOnEdges")<<data._index);
3313 for ( size_t i = 0; i < data._edges.size(); ++i )
3315 _LayerEdge* edge = data._edges[i];
3316 if ( !edge->IsOnEdge() || !edge->_sWOL.IsNull() ) continue;
3317 const SMDS_MeshNode* tgt1 = edge->_nodes.back();
3318 for ( int j = 0; j < 2; ++j ) // loop on _2NearEdges
3320 const SMDS_MeshNode* tgt2 = edge->_2neibors->_nodes[j];
3321 pair< set< SMESH_TLink >::iterator, bool > link_isnew =
3322 extrudedLinks.insert( SMESH_TLink( tgt1, tgt2 ));
3323 if ( !link_isnew.second )
3325 extrudedLinks.erase( link_isnew.first );
3326 continue; // already extruded and will no more encounter
3328 // a _LayerEdge containg tgt2
3329 _LayerEdge* neiborEdge = edge->_2neibors->_edges[j];
3331 _TmpMeshFaceOnEdge* f = new _TmpMeshFaceOnEdge( edge, neiborEdge, --_tmpFaceID );
3332 tmpFaces.push_back( f );
3334 dumpCmd(SMESH_Comment("mesh.AddFace([ ")
3335 <<f->_nn[0]->GetID()<<", "<<f->_nn[1]->GetID()<<", "
3336 <<f->_nn[2]->GetID()<<", "<<f->_nn[3]->GetID()<<" ])");
3341 // Check if _LayerEdge's based on EDGE's intersects tmpFaces.
3342 // Perform two loops on _LayerEdge on EDGE's:
3343 // 1) to find and fix intersection
3344 // 2) to check that no new intersection appears as result of 1)
3346 SMDS_ElemIteratorPtr fIt( new SMDS_ElementVectorIterator( tmpFaces.begin(),
3348 auto_ptr<SMESH_ElementSearcher> searcher
3349 ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), fIt ));
3351 // 1) Find intersections
3353 const SMDS_MeshElement* face;
3354 typedef map< _LayerEdge*, set< _LayerEdge*, _LayerEdgeCmp >, _LayerEdgeCmp > TLEdge2LEdgeSet;
3355 TLEdge2LEdgeSet edge2CloseEdge;
3357 const double eps = data._epsilon * data._epsilon;
3358 for ( size_t i = 0; i < data._edges.size(); ++i )
3360 _LayerEdge* edge = data._edges[i];
3361 if (( !edge->IsOnEdge() ) &&
3362 ( edge->_sWOL.IsNull() || edge->_sWOL.ShapeType() != TopAbs_FACE ))
3364 if ( edge->FindIntersection( *searcher, dist, eps, &face ))
3366 const _TmpMeshFaceOnEdge* f = (const _TmpMeshFaceOnEdge*) face;
3367 set< _LayerEdge*, _LayerEdgeCmp > & ee = edge2CloseEdge[ edge ];
3368 ee.insert( f->_le1 );
3369 ee.insert( f->_le2 );
3370 if ( f->_le1->IsOnEdge() && f->_le1->_sWOL.IsNull() )
3371 edge2CloseEdge[ f->_le1 ].insert( edge );
3372 if ( f->_le2->IsOnEdge() && f->_le2->_sWOL.IsNull() )
3373 edge2CloseEdge[ f->_le2 ].insert( edge );
3377 // Set _LayerEdge._normal
3379 if ( !edge2CloseEdge.empty() )
3381 dumpFunction(SMESH_Comment("updateNormals")<<data._index);
3383 set< TGeomID > shapesToSmooth;
3385 // vector to store new _normal and _cosin for each edge in edge2CloseEdge
3386 vector< pair< _LayerEdge*, _LayerEdge > > edge2newEdge( edge2CloseEdge.size() );
3388 TLEdge2LEdgeSet::iterator e2ee = edge2CloseEdge.begin();
3389 for ( size_t iE = 0; e2ee != edge2CloseEdge.end(); ++e2ee, ++iE )
3391 _LayerEdge* edge1 = e2ee->first;
3392 _LayerEdge* edge2 = 0;
3393 set< _LayerEdge*, _LayerEdgeCmp >& ee = e2ee->second;
3395 edge2newEdge[ iE ].first = NULL;
3397 // find EDGEs the edges reside
3398 // TopoDS_Edge E1, E2;
3399 // TopoDS_Shape S = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
3400 // if ( S.ShapeType() != TopAbs_EDGE )
3401 // continue; // TODO: find EDGE by VERTEX
3402 // E1 = TopoDS::Edge( S );
3403 set< _LayerEdge*, _LayerEdgeCmp >::iterator eIt = ee.begin();
3404 for ( ; !edge2 && eIt != ee.end(); ++eIt )
3406 if ( edge1->_sWOL == (*eIt)->_sWOL )
3409 if ( !edge2 ) continue;
3411 edge2newEdge[ iE ].first = edge1;
3412 _LayerEdge& newEdge = edge2newEdge[ iE ].second;
3413 // while ( E2.IsNull() && eIt != ee.end())
3415 // _LayerEdge* e2 = *eIt++;
3416 // TopoDS_Shape S = helper.GetSubShapeByNode( e2->_nodes[0], getMeshDS() );
3417 // if ( S.ShapeType() == TopAbs_EDGE )
3418 // E2 = TopoDS::Edge( S ), edge2 = e2;
3420 // if ( E2.IsNull() ) continue; // TODO: find EDGE by VERTEX
3422 // find 3 FACEs sharing 2 EDGEs
3424 // TopoDS_Face FF1[2], FF2[2];
3425 // PShapeIteratorPtr fIt = helper.GetAncestors(E1, *_mesh, TopAbs_FACE);
3426 // while ( fIt->more() && FF1[1].IsNull() )
3428 // const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
3429 // if ( helper.IsSubShape( *F, data._solid))
3430 // FF1[ FF1[0].IsNull() ? 0 : 1 ] = *F;
3432 // fIt = helper.GetAncestors(E2, *_mesh, TopAbs_FACE);
3433 // while ( fIt->more() && FF2[1].IsNull())
3435 // const TopoDS_Face *F = (const TopoDS_Face*) fIt->next();
3436 // if ( helper.IsSubShape( *F, data._solid))
3437 // FF2[ FF2[0].IsNull() ? 0 : 1 ] = *F;
3439 // // exclude a FACE common to E1 and E2 (put it to FFn[1] )
3440 // if ( FF1[0].IsSame( FF2[0]) || FF1[0].IsSame( FF2[1]))
3441 // std::swap( FF1[0], FF1[1] );
3442 // if ( FF2[0].IsSame( FF1[0]) )
3443 // std::swap( FF2[0], FF2[1] );
3444 // if ( FF1[0].IsNull() || FF2[0].IsNull() )
3447 // get a new normal for edge1
3449 gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal;
3450 // if ( edge1->_cosin < 0 )
3451 // dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized();
3452 // if ( edge2->_cosin < 0 )
3453 // dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized();
3455 double cos1 = Abs( edge1->_cosin ), cos2 = Abs( edge2->_cosin );
3456 double wgt1 = ( cos1 + 0.001 ) / ( cos1 + cos2 + 0.002 );
3457 double wgt2 = ( cos2 + 0.001 ) / ( cos1 + cos2 + 0.002 );
3458 newEdge._normal = ( wgt1 * dir1 + wgt2 * dir2 ).XYZ();
3459 newEdge._normal.Normalize();
3461 // cout << edge1->_nodes[0]->GetID() << " "
3462 // << edge2->_nodes[0]->GetID() << " NORM: "
3463 // << newEdge._normal.X() << ", " << newEdge._normal.Y() << ", " << newEdge._normal.Z() << endl;
3466 if ( cos1 < theMinSmoothCosin )
3468 newEdge._cosin = edge2->_cosin;
3470 else if ( cos2 > theMinSmoothCosin ) // both cos1 and cos2 > theMinSmoothCosin
3472 // gp_Vec dirInFace;
3473 // if ( edge1->_cosin < 0 )
3474 // dirInFace = dir1;
3476 // dirInFace = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok );
3477 // double angle = dirInFace.Angle( edge1->_normal ); // [0,PI]
3478 // edge1->SetCosin( Cos( angle ));
3479 //newEdge._cosin = 0; // ???????????
3480 newEdge._cosin = ( wgt1 * cos1 + wgt2 * cos2 ) * edge1->_cosin / cos1;
3484 newEdge._cosin = edge1->_cosin;
3487 // find shapes that need smoothing due to change of _normal
3488 if ( edge1->_cosin < theMinSmoothCosin &&
3489 newEdge._cosin > theMinSmoothCosin )
3491 if ( edge1->_sWOL.IsNull() )
3493 SMDS_ElemIteratorPtr fIt = edge1->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
3494 while ( fIt->more() )
3495 shapesToSmooth.insert( fIt->next()->getshapeId() );
3496 //limitStepSize( data, fIt->next(), edge1->_cosin ); // too late
3498 else // edge1 inflates along a FACE
3500 TopoDS_Shape V = helper.GetSubShapeByNode( edge1->_nodes[0], getMeshDS() );
3501 PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE );
3502 while ( const TopoDS_Shape* E = eIt->next() )
3504 if ( !helper.IsSubShape( *E, /*FACE=*/edge1->_sWOL ))
3506 gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V ));
3507 double angle = edgeDir.Angle( newEdge._normal ); // [0,PI]
3508 if ( angle < M_PI / 2 )
3509 shapesToSmooth.insert( getMeshDS()->ShapeToIndex( *E ));
3515 data.AddShapesToSmooth( shapesToSmooth );
3517 // Update data of edges depending on a new _normal
3519 for ( size_t iE = 0; iE < edge2newEdge.size(); ++iE )
3521 _LayerEdge* edge1 = edge2newEdge[ iE ].first;
3522 _LayerEdge& newEdge = edge2newEdge[ iE ].second;
3523 if ( !edge1 ) continue;
3525 edge1->_normal = newEdge._normal;
3526 edge1->SetCosin( newEdge._cosin );
3527 edge1->InvalidateStep( 1 );
3529 edge1->SetNewLength( data._stepSize, helper );
3530 if ( edge1->IsOnEdge() )
3532 const SMDS_MeshNode * n1 = edge1->_2neibors->_edges[0]->_nodes[0];
3533 const SMDS_MeshNode * n2 = edge1->_2neibors->_edges[1]->_nodes[0];
3534 edge1->SetDataByNeighbors( n1, n2, helper );
3537 // Update normals and other dependent data of not intersecting _LayerEdge's
3538 // neighboring the intersecting ones
3540 if ( !edge1->_2neibors )
3542 for ( int j = 0; j < 2; ++j ) // loop on 2 neighbors
3544 _LayerEdge* neighbor = edge1->_2neibors->_edges[j];
3545 if ( edge2CloseEdge.count ( neighbor ))
3546 continue; // j-th neighbor is also intersected
3547 _LayerEdge* prevEdge = edge1;
3548 const int nbSteps = 10;
3549 for ( int step = nbSteps; step; --step ) // step from edge1 in j-th direction
3551 if ( !neighbor->_2neibors )
3552 break; // neighbor is on VERTEX
3554 _LayerEdge* nextEdge = neighbor->_2neibors->_edges[iNext];
3555 if ( nextEdge == prevEdge )
3556 nextEdge = neighbor->_2neibors->_edges[ ++iNext ];
3557 double r = double(step-1)/nbSteps;
3558 if ( !nextEdge->_2neibors )
3561 gp_XYZ newNorm = prevEdge->_normal * r + nextEdge->_normal * (1-r);
3562 newNorm.Normalize();
3564 neighbor->_normal = newNorm;
3565 neighbor->SetCosin( prevEdge->_cosin * r + nextEdge->_cosin * (1-r) );
3566 neighbor->SetDataByNeighbors( prevEdge->_nodes[0], nextEdge->_nodes[0], helper );
3568 neighbor->InvalidateStep( 1 );
3570 neighbor->SetNewLength( data._stepSize, helper );
3572 // goto the next neighbor
3573 prevEdge = neighbor;
3574 neighbor = nextEdge;
3580 // 2) Check absence of intersections
3583 for ( size_t i = 0 ; i < tmpFaces.size(); ++i )
3589 //================================================================================
3591 * \brief Modify normals of _LayerEdge's on _ConvexFace's
3593 //================================================================================
3595 bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData& data,
3596 SMESH_MesherHelper& helper,
3599 SMESHDS_Mesh* meshDS = helper.GetMeshDS();
3602 map< TGeomID, _ConvexFace >::iterator id2face = data._convexFaces.begin();
3603 for ( ; id2face != data._convexFaces.end(); ++id2face )
3605 _ConvexFace & convFace = (*id2face).second;
3606 if ( convFace._normalsFixed )
3607 continue; // already fixed
3608 if ( convFace.CheckPrisms() )
3609 continue; // nothing to fix
3611 convFace._normalsFixed = true;
3613 BRepAdaptor_Surface surface ( convFace._face, false );
3614 BRepLProp_SLProps surfProp( surface, 2, 1e-6 );
3616 // check if the convex FACE is of spherical shape
3618 Bnd_B3d centersBox; // bbox of centers of curvature of _LayerEdge's on VERTEXes
3623 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.begin();
3624 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3626 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3628 if ( meshDS->IndexToShape( id2end->first ).ShapeType() == TopAbs_VERTEX )
3630 _LayerEdge* ledge = data._edges[ iBeg ];
3631 if ( convFace.GetCenterOfCurvature( ledge, surfProp, helper, center ))
3632 centersBox.Add( center );
3634 for ( ; iBeg < iEnd; ++iBeg )
3635 nodesBox.Add( SMESH_TNodeXYZ( data._edges[ iBeg ]->_nodes[0] ));
3637 if ( centersBox.IsVoid() )
3639 debugMsg( "Error: centersBox.IsVoid()" );
3642 const bool isSpherical =
3643 ( centersBox.SquareExtent() < 1e-6 * nodesBox.SquareExtent() );
3645 int nbEdges = helper.Count( convFace._face, TopAbs_EDGE, /*ignoreSame=*/false );
3646 vector < _CentralCurveOnEdge > centerCurves( nbEdges );
3650 // set _LayerEdge::_normal as average of all normals
3652 // WARNING: different density of nodes on EDGEs is not taken into account that
3653 // can lead to an improper new normal
3655 gp_XYZ avgNormal( 0,0,0 );
3657 id2end = convFace._subIdToEdgeEnd.begin();
3658 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3660 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3661 // set data of _CentralCurveOnEdge
3662 const TopoDS_Shape& S = meshDS->IndexToShape( id2end->first );
3663 if ( S.ShapeType() == TopAbs_EDGE )
3665 _CentralCurveOnEdge& ceCurve = centerCurves[ nbEdges++ ];
3666 ceCurve.SetShapes( TopoDS::Edge(S), convFace, data, helper );
3667 if ( !data._edges[ iBeg ]->_sWOL.IsNull() )
3668 ceCurve._adjFace.Nullify();
3670 ceCurve._ledges.insert( ceCurve._ledges.end(),
3671 &data._edges[ iBeg ], &data._edges[ iEnd ]);
3673 // summarize normals
3674 for ( ; iBeg < iEnd; ++iBeg )
3675 avgNormal += data._edges[ iBeg ]->_normal;
3677 double normSize = avgNormal.SquareModulus();
3678 if ( normSize < 1e-200 )
3680 debugMsg( "updateNormalsOfConvexFaces(): zero avgNormal" );
3683 avgNormal /= Sqrt( normSize );
3685 // compute new _LayerEdge::_cosin on EDGEs
3686 double avgCosin = 0;
3689 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3691 _CentralCurveOnEdge& ceCurve = centerCurves[ iE ];
3692 if ( ceCurve._adjFace.IsNull() )
3694 for ( size_t iLE = 0; iLE < ceCurve._ledges.size(); ++iLE )
3696 const SMDS_MeshNode* node = ceCurve._ledges[ iLE ]->_nodes[0];
3697 inFaceDir = getFaceDir( ceCurve._adjFace, ceCurve._edge, node, helper, isOK );
3700 double angle = inFaceDir.Angle( avgNormal ); // [0,PI]
3701 ceCurve._ledges[ iLE ]->_cosin = Cos( angle );
3702 avgCosin += ceCurve._ledges[ iLE ]->_cosin;
3708 avgCosin /= nbCosin;
3710 // set _LayerEdge::_normal = avgNormal
3711 id2end = convFace._subIdToEdgeEnd.begin();
3712 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3714 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3715 const TopoDS_Shape& S = meshDS->IndexToShape( id2end->first );
3716 if ( S.ShapeType() != TopAbs_EDGE )
3717 for ( int i = iBeg; i < iEnd; ++i )
3718 data._edges[ i ]->_cosin = avgCosin;
3720 for ( ; iBeg < iEnd; ++iBeg )
3721 data._edges[ iBeg ]->_normal = avgNormal;
3724 else // if ( isSpherical )
3726 // We suppose that centers of curvature at all points of the FACE
3727 // lie on some curve, let's call it "central curve". For all _LayerEdge's
3728 // having a common center of curvature we define the same new normal
3729 // as a sum of normals of _LayerEdge's on EDGEs among them.
3731 // get all centers of curvature for each EDGE
3733 helper.SetSubShape( convFace._face );
3734 _LayerEdge* vertexLEdges[2], **edgeLEdge, **edgeLEdgeEnd;
3736 TopExp_Explorer edgeExp( convFace._face, TopAbs_EDGE );
3737 for ( int iE = 0; edgeExp.More(); edgeExp.Next(), ++iE )
3739 const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
3741 // set adjacent FACE
3742 centerCurves[ iE ].SetShapes( edge, convFace, data, helper );
3744 // get _LayerEdge's of the EDGE
3745 TGeomID edgeID = meshDS->ShapeToIndex( edge );
3746 id2end = convFace._subIdToEdgeEnd.find( edgeID );
3747 if ( id2end == convFace._subIdToEdgeEnd.end() )
3749 // no _LayerEdge's on EDGE, use _LayerEdge's on VERTEXes
3750 for ( int iV = 0; iV < 2; ++iV )
3752 TopoDS_Vertex v = helper.IthVertex( iV, edge );
3753 TGeomID vID = meshDS->ShapeToIndex( v );
3754 int end = convFace._subIdToEdgeEnd[ vID ];
3755 int iBeg = end > 0 ? data._endEdgeOnShape[ end-1 ] : 0;
3756 vertexLEdges[ iV ] = data._edges[ iBeg ];
3758 edgeLEdge = &vertexLEdges[0];
3759 edgeLEdgeEnd = edgeLEdge + 2;
3761 centerCurves[ iE ]._adjFace.Nullify();
3765 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3766 if ( id2end->second >= data._nbShapesToSmooth )
3767 data.SortOnEdge( edge, iBeg, iEnd, helper );
3768 edgeLEdge = &data._edges[ iBeg ];
3769 edgeLEdgeEnd = edgeLEdge + iEnd - iBeg;
3770 vertexLEdges[0] = data._edges[ iBeg ]->_2neibors->_edges[0];
3771 vertexLEdges[1] = data._edges[ iEnd-1 ]->_2neibors->_edges[1];
3773 if ( ! data._edges[ iBeg ]->_sWOL.IsNull() )
3774 centerCurves[ iE ]._adjFace.Nullify();
3777 // Get curvature centers
3781 if ( edgeLEdge[0]->IsOnEdge() &&
3782 convFace.GetCenterOfCurvature( vertexLEdges[0], surfProp, helper, center ))
3784 centerCurves[ iE ].Append( center, vertexLEdges[0] );
3785 centersBox.Add( center );
3787 for ( ; edgeLEdge < edgeLEdgeEnd; ++edgeLEdge )
3788 if ( convFace.GetCenterOfCurvature( *edgeLEdge, surfProp, helper, center ))
3789 { // EDGE or VERTEXes
3790 centerCurves[ iE ].Append( center, *edgeLEdge );
3791 centersBox.Add( center );
3793 if ( edgeLEdge[-1]->IsOnEdge() &&
3794 convFace.GetCenterOfCurvature( vertexLEdges[1], surfProp, helper, center ))
3796 centerCurves[ iE ].Append( center, vertexLEdges[1] );
3797 centersBox.Add( center );
3799 centerCurves[ iE ]._isDegenerated =
3800 ( centersBox.IsVoid() || centersBox.SquareExtent() < 1e-6 * nodesBox.SquareExtent() );
3802 } // loop on EDGES of convFace._face to set up data of centerCurves
3804 // Compute new normals for _LayerEdge's on EDGEs
3806 double avgCosin = 0;
3809 for ( size_t iE1 = 0; iE1 < centerCurves.size(); ++iE1 )
3811 _CentralCurveOnEdge& ceCurve = centerCurves[ iE1 ];
3812 if ( ceCurve._isDegenerated )
3814 const vector< gp_Pnt >& centers = ceCurve._curvaCenters;
3815 vector< gp_XYZ > & newNormals = ceCurve._normals;
3816 for ( size_t iC1 = 0; iC1 < centers.size(); ++iC1 )
3819 for ( size_t iE2 = 0; iE2 < centerCurves.size() && !isOK; ++iE2 )
3822 isOK = centerCurves[ iE2 ].FindNewNormal( centers[ iC1 ], newNormals[ iC1 ]);
3824 if ( isOK && !ceCurve._adjFace.IsNull() )
3826 // compute new _LayerEdge::_cosin
3827 const SMDS_MeshNode* node = ceCurve._ledges[ iC1 ]->_nodes[0];
3828 inFaceDir = getFaceDir( ceCurve._adjFace, ceCurve._edge, node, helper, isOK );
3831 double angle = inFaceDir.Angle( newNormals[ iC1 ] ); // [0,PI]
3832 ceCurve._ledges[ iC1 ]->_cosin = Cos( angle );
3833 avgCosin += ceCurve._ledges[ iC1 ]->_cosin;
3839 // set new normals to _LayerEdge's of NOT degenerated central curves
3840 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3842 if ( centerCurves[ iE ]._isDegenerated )
3844 for ( size_t iLE = 0; iLE < centerCurves[ iE ]._ledges.size(); ++iLE )
3845 centerCurves[ iE ]._ledges[ iLE ]->_normal = centerCurves[ iE ]._normals[ iLE ];
3847 // set new normals to _LayerEdge's of degenerated central curves
3848 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3850 if ( !centerCurves[ iE ]._isDegenerated ||
3851 centerCurves[ iE ]._ledges.size() < 3 )
3853 // new normal is an average of new normals at VERTEXes that
3854 // was computed on non-degenerated _CentralCurveOnEdge's
3855 gp_XYZ newNorm = ( centerCurves[ iE ]._ledges.front()->_normal +
3856 centerCurves[ iE ]._ledges.back ()->_normal );
3857 double sz = newNorm.Modulus();
3861 double newCosin = ( 0.5 * centerCurves[ iE ]._ledges.front()->_cosin +
3862 0.5 * centerCurves[ iE ]._ledges.back ()->_cosin );
3863 for ( size_t iLE = 1, nb = centerCurves[ iE ]._ledges.size() - 1; iLE < nb; ++iLE )
3865 centerCurves[ iE ]._ledges[ iLE ]->_normal = newNorm;
3866 centerCurves[ iE ]._ledges[ iLE ]->_cosin = newCosin;
3870 // Find new normals for _LayerEdge's based on FACE
3873 avgCosin /= nbCosin;
3874 const TGeomID faceID = meshDS->ShapeToIndex( convFace._face );
3875 map< TGeomID, int >::iterator id2end = convFace._subIdToEdgeEnd.find( faceID );
3876 if ( id2end != convFace._subIdToEdgeEnd.end() )
3880 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3881 for ( ; iBeg < iEnd; ++iBeg )
3883 _LayerEdge* ledge = data._edges[ iBeg ];
3884 if ( !convFace.GetCenterOfCurvature( ledge, surfProp, helper, center ))
3886 for ( size_t i = 0; i < centerCurves.size(); ++i, ++iE )
3888 iE = iE % centerCurves.size();
3889 if ( centerCurves[ iE ]._isDegenerated )
3891 newNorm.SetCoord( 0,0,0 );
3892 if ( centerCurves[ iE ].FindNewNormal( center, newNorm ))
3894 ledge->_normal = newNorm;
3895 ledge->_cosin = avgCosin;
3902 } // not a quasi-spherical FACE
3904 // Update _LayerEdge's data according to a new normal
3906 dumpFunction(SMESH_Comment("updateNormalsOfConvexFaces")<<data._index
3907 <<"_F"<<meshDS->ShapeToIndex( convFace._face ));
3909 id2end = convFace._subIdToEdgeEnd.begin();
3910 for ( ; id2end != convFace._subIdToEdgeEnd.end(); ++id2end )
3912 data.GetEdgesOnShape( id2end->second, iBeg, iEnd );
3913 for ( ; iBeg < iEnd; ++iBeg )
3915 _LayerEdge* & ledge = data._edges[ iBeg ];
3916 double len = ledge->_len;
3917 ledge->InvalidateStep( stepNb + 1, /*restoreLength=*/true );
3918 ledge->SetCosin( ledge->_cosin );
3919 ledge->SetNewLength( len, helper );
3922 } // loop on sub-shapes of convFace._face
3924 // Find FACEs adjacent to convFace._face that got necessity to smooth
3925 // as a result of normals modification
3927 set< TGeomID > adjFacesToSmooth;
3928 for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
3930 if ( centerCurves[ iE ]._adjFace.IsNull() ||
3931 centerCurves[ iE ]._adjFaceToSmooth )
3933 for ( size_t iLE = 0; iLE < centerCurves[ iE ]._ledges.size(); ++iLE )
3935 if ( centerCurves[ iE ]._ledges[ iLE ]->_cosin > theMinSmoothCosin )
3937 adjFacesToSmooth.insert( meshDS->ShapeToIndex( centerCurves[ iE ]._adjFace ));
3942 data.AddShapesToSmooth( adjFacesToSmooth );
3947 } // loop on data._convexFaces
3952 //================================================================================
3954 * \brief Finds a center of curvature of a surface at a _LayerEdge
3956 //================================================================================
3958 bool _ConvexFace::GetCenterOfCurvature( _LayerEdge* ledge,
3959 BRepLProp_SLProps& surfProp,
3960 SMESH_MesherHelper& helper,
3961 gp_Pnt & center ) const
3963 gp_XY uv = helper.GetNodeUV( _face, ledge->_nodes[0] );
3964 surfProp.SetParameters( uv.X(), uv.Y() );
3965 if ( !surfProp.IsCurvatureDefined() )
3968 const double oriFactor = ( _face.Orientation() == TopAbs_REVERSED ? +1. : -1. );
3969 double surfCurvatureMax = surfProp.MaxCurvature() * oriFactor;
3970 double surfCurvatureMin = surfProp.MinCurvature() * oriFactor;
3971 if ( surfCurvatureMin > surfCurvatureMax )
3972 center = surfProp.Value().Translated( surfProp.Normal().XYZ() / surfCurvatureMin * oriFactor );
3974 center = surfProp.Value().Translated( surfProp.Normal().XYZ() / surfCurvatureMax * oriFactor );
3979 //================================================================================
3981 * \brief Check that prisms are not distorted
3983 //================================================================================
3985 bool _ConvexFace::CheckPrisms() const
3987 for ( size_t i = 0; i < _simplexTestEdges.size(); ++i )
3989 const _LayerEdge* edge = _simplexTestEdges[i];
3990 SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() );
3991 for ( size_t j = 0; j < edge->_simplices.size(); ++j )
3992 if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ ))
3994 debugMsg( "Bad simplex of _simplexTestEdges ("
3995 << " "<< edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
3996 << " "<< edge->_simplices[j]._nPrev->GetID()
3997 << " "<< edge->_simplices[j]._nNext->GetID() << " )" );
4004 //================================================================================
4006 * \brief Try to compute a new normal by interpolating normals of _LayerEdge's
4007 * stored in this _CentralCurveOnEdge.
4008 * \param [in] center - curvature center of a point of another _CentralCurveOnEdge.
4009 * \param [in,out] newNormal - current normal at this point, to be redefined
4010 * \return bool - true if succeeded.
4012 //================================================================================
4014 bool _CentralCurveOnEdge::FindNewNormal( const gp_Pnt& center, gp_XYZ& newNormal )
4016 if ( this->_isDegenerated )
4019 // find two centers the given one lies between
4021 for ( size_t i = 0, nb = _curvaCenters.size()-1; i < nb; ++i )
4023 double sl2 = 1.001 * _segLength2[ i ];
4025 double d1 = center.SquareDistance( _curvaCenters[ i ]);
4029 double d2 = center.SquareDistance( _curvaCenters[ i+1 ]);
4030 if ( d2 > sl2 || d2 + d1 < 1e-100 )
4035 double r = d1 / ( d1 + d2 );
4036 gp_XYZ norm = (( 1. - r ) * _ledges[ i ]->_normal +
4037 ( r ) * _ledges[ i+1 ]->_normal );
4041 double sz = newNormal.Modulus();
4050 //================================================================================
4052 * \brief Set shape members
4054 //================================================================================
4056 void _CentralCurveOnEdge::SetShapes( const TopoDS_Edge& edge,
4057 const _ConvexFace& convFace,
4058 const _SolidData& data,
4059 SMESH_MesherHelper& helper)
4063 PShapeIteratorPtr fIt = helper.GetAncestors( edge, *helper.GetMesh(), TopAbs_FACE );
4064 while ( const TopoDS_Shape* F = fIt->next())
4065 if ( !convFace._face.IsSame( *F ))
4067 _adjFace = TopoDS::Face( *F );
4068 _adjFaceToSmooth = false;
4069 // _adjFace already in a smoothing queue ?
4071 TGeomID adjFaceID = helper.GetMeshDS()->ShapeToIndex( *F );
4072 if ( data.GetShapeEdges( adjFaceID, end ))
4073 _adjFaceToSmooth = ( end < data._nbShapesToSmooth );
4078 //================================================================================
4080 * \brief Looks for intersection of it's last segment with faces
4081 * \param distance - returns shortest distance from the last node to intersection
4083 //================================================================================
4085 bool _LayerEdge::FindIntersection( SMESH_ElementSearcher& searcher,
4087 const double& epsilon,
4088 const SMDS_MeshElement** face)
4090 vector< const SMDS_MeshElement* > suspectFaces;
4092 gp_Ax1 lastSegment = LastSegment(segLen);
4093 searcher.GetElementsNearLine( lastSegment, SMDSAbs_Face, suspectFaces );
4095 bool segmentIntersected = false;
4096 distance = Precision::Infinite();
4097 int iFace = -1; // intersected face
4098 for ( size_t j = 0 ; j < suspectFaces.size() /*&& !segmentIntersected*/; ++j )
4100 const SMDS_MeshElement* face = suspectFaces[j];
4101 if ( face->GetNodeIndex( _nodes.back() ) >= 0 ||
4102 face->GetNodeIndex( _nodes[0] ) >= 0 )
4103 continue; // face sharing _LayerEdge node
4104 const int nbNodes = face->NbCornerNodes();
4105 bool intFound = false;
4107 SMDS_MeshElement::iterator nIt = face->begin_nodes();
4110 intFound = SegTriaInter( lastSegment, *nIt++, *nIt++, *nIt++, dist, epsilon );
4114 const SMDS_MeshNode* tria[3];
4117 for ( int n2 = 2; n2 < nbNodes && !intFound; ++n2 )
4120 intFound = SegTriaInter(lastSegment, tria[0], tria[1], tria[2], dist, epsilon );
4126 if ( dist < segLen*(1.01) && dist > -(_len*_lenFactor-segLen) )
4127 segmentIntersected = true;
4128 if ( distance > dist )
4129 distance = dist, iFace = j;
4132 if ( iFace != -1 && face ) *face = suspectFaces[iFace];
4134 if ( segmentIntersected )
4137 SMDS_MeshElement::iterator nIt = suspectFaces[iFace]->begin_nodes();
4138 gp_XYZ intP( lastSegment.Location().XYZ() + lastSegment.Direction().XYZ() * distance );
4139 cout << "nodes: tgt " << _nodes.back()->GetID() << " src " << _nodes[0]->GetID()
4140 << ", intersection with face ("
4141 << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()
4142 << ") at point (" << intP.X() << ", " << intP.Y() << ", " << intP.Z()
4143 << ") distance = " << distance - segLen<< endl;
4149 return segmentIntersected;
4152 //================================================================================
4154 * \brief Returns size and direction of the last segment
4156 //================================================================================
4158 gp_Ax1 _LayerEdge::LastSegment(double& segLen) const
4160 // find two non-coincident positions
4161 gp_XYZ orig = _pos.back();
4163 int iPrev = _pos.size() - 2;
4164 while ( iPrev >= 0 )
4166 dir = orig - _pos[iPrev];
4167 if ( dir.SquareModulus() > 1e-100 )
4177 segDir.SetLocation( SMESH_TNodeXYZ( _nodes[0] ));
4178 segDir.SetDirection( _normal );
4183 gp_Pnt pPrev = _pos[ iPrev ];
4184 if ( !_sWOL.IsNull() )
4186 TopLoc_Location loc;
4187 if ( _sWOL.ShapeType() == TopAbs_EDGE )
4190 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
4191 pPrev = curve->Value( pPrev.X() ).Transformed( loc );
4195 Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
4196 pPrev = surface->Value( pPrev.X(), pPrev.Y() ).Transformed( loc );
4198 dir = SMESH_TNodeXYZ( _nodes.back() ) - pPrev.XYZ();
4200 segDir.SetLocation( pPrev );
4201 segDir.SetDirection( dir );
4202 segLen = dir.Modulus();
4208 //================================================================================
4210 * \brief Test intersection of the last segment with a given triangle
4211 * using Moller-Trumbore algorithm
4212 * Intersection is detected if distance to intersection is less than _LayerEdge._len
4214 //================================================================================
4216 bool _LayerEdge::SegTriaInter( const gp_Ax1& lastSegment,
4217 const SMDS_MeshNode* n0,
4218 const SMDS_MeshNode* n1,
4219 const SMDS_MeshNode* n2,
4221 const double& EPSILON) const
4223 //const double EPSILON = 1e-6;
4225 gp_XYZ orig = lastSegment.Location().XYZ();
4226 gp_XYZ dir = lastSegment.Direction().XYZ();
4228 SMESH_TNodeXYZ vert0( n0 );
4229 SMESH_TNodeXYZ vert1( n1 );
4230 SMESH_TNodeXYZ vert2( n2 );
4232 /* calculate distance from vert0 to ray origin */
4233 gp_XYZ tvec = orig - vert0;
4235 //if ( tvec * dir > EPSILON )
4236 // intersected face is at back side of the temporary face this _LayerEdge belongs to
4239 gp_XYZ edge1 = vert1 - vert0;
4240 gp_XYZ edge2 = vert2 - vert0;
4242 /* begin calculating determinant - also used to calculate U parameter */
4243 gp_XYZ pvec = dir ^ edge2;
4245 /* if determinant is near zero, ray lies in plane of triangle */
4246 double det = edge1 * pvec;
4248 if (det > -EPSILON && det < EPSILON)
4250 double inv_det = 1.0 / det;
4252 /* calculate U parameter and test bounds */
4253 double u = ( tvec * pvec ) * inv_det;
4254 //if (u < 0.0 || u > 1.0)
4255 if (u < -EPSILON || u > 1.0 + EPSILON)
4258 /* prepare to test V parameter */
4259 gp_XYZ qvec = tvec ^ edge1;
4261 /* calculate V parameter and test bounds */
4262 double v = (dir * qvec) * inv_det;
4263 //if ( v < 0.0 || u + v > 1.0 )
4264 if ( v < -EPSILON || u + v > 1.0 + EPSILON)
4267 /* calculate t, ray intersects triangle */
4268 t = (edge2 * qvec) * inv_det;
4274 //================================================================================
4276 * \brief Perform smooth of _LayerEdge's based on EDGE's
4277 * \retval bool - true if node has been moved
4279 //================================================================================
4281 bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface,
4282 const TopoDS_Face& F,
4283 SMESH_MesherHelper& helper)
4285 ASSERT( IsOnEdge() );
4287 SMDS_MeshNode* tgtNode = const_cast<SMDS_MeshNode*>( _nodes.back() );
4288 SMESH_TNodeXYZ oldPos( tgtNode );
4289 double dist01, distNewOld;
4291 SMESH_TNodeXYZ p0( _2neibors->_nodes[0]);
4292 SMESH_TNodeXYZ p1( _2neibors->_nodes[1]);
4293 dist01 = p0.Distance( _2neibors->_nodes[1] );
4295 gp_Pnt newPos = p0 * _2neibors->_wgt[0] + p1 * _2neibors->_wgt[1];
4296 double lenDelta = 0;
4299 //lenDelta = _curvature->lenDelta( _len );
4300 lenDelta = _curvature->lenDeltaByDist( dist01 );
4301 newPos.ChangeCoord() += _normal * lenDelta;
4304 distNewOld = newPos.Distance( oldPos );
4308 if ( _2neibors->_plnNorm )
4310 // put newPos on the plane defined by source node and _plnNorm
4311 gp_XYZ new2src = SMESH_TNodeXYZ( _nodes[0] ) - newPos.XYZ();
4312 double new2srcProj = (*_2neibors->_plnNorm) * new2src;
4313 newPos.ChangeCoord() += (*_2neibors->_plnNorm) * new2srcProj;
4315 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4316 _pos.back() = newPos.XYZ();
4320 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4321 gp_XY uv( Precision::Infinite(), 0 );
4322 helper.CheckNodeUV( F, tgtNode, uv, 1e-10, /*force=*/true );
4323 _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
4325 newPos = surface->Value( uv.X(), uv.Y() );
4326 tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() );
4329 if ( _curvature && lenDelta < 0 )
4331 gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
4332 _len -= prevPos.Distance( oldPos );
4333 _len += prevPos.Distance( newPos );
4335 bool moved = distNewOld > dist01/50;
4337 dumpMove( tgtNode ); // debug
4342 //================================================================================
4344 * \brief Perform laplacian smooth in 3D of nodes inflated from FACE
4345 * \retval bool - true if _tgtNode has been moved
4347 //================================================================================
4349 bool _LayerEdge::Smooth(int& badNb)
4351 if ( _simplices.size() < 2 )
4352 return false; // _LayerEdge inflated along EDGE or FACE
4354 // compute new position for the last _pos
4355 gp_XYZ newPos (0,0,0);
4356 for ( size_t i = 0; i < _simplices.size(); ++i )
4357 newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev );
4358 newPos /= _simplices.size();
4360 const gp_XYZ& curPos ( _pos.back() );
4361 const gp_Pnt prevPos( _pos[ _pos.size()-2 ]);
4364 double delta = _curvature->lenDelta( _len );
4366 newPos += _normal * delta;
4369 double segLen = _normal * ( newPos - prevPos.XYZ() );
4370 if ( segLen + delta > 0 )
4371 newPos += _normal * delta;
4373 // double segLenChange = _normal * ( curPos - newPos );
4374 // newPos += 0.5 * _normal * segLenChange;
4377 // count quality metrics (orientation) of tetras around _tgtNode
4379 for ( size_t i = 0; i < _simplices.size(); ++i )
4380 nbOkBefore += _simplices[i].IsForward( _nodes[0], &curPos );
4383 for ( size_t i = 0; i < _simplices.size(); ++i )
4384 nbOkAfter += _simplices[i].IsForward( _nodes[0], &newPos );
4386 if ( nbOkAfter < nbOkBefore )
4389 SMDS_MeshNode* n = const_cast< SMDS_MeshNode* >( _nodes.back() );
4391 _len -= prevPos.Distance(SMESH_TNodeXYZ( n ));
4392 _len += prevPos.Distance(newPos);
4394 n->setXYZ( newPos.X(), newPos.Y(), newPos.Z());
4395 _pos.back() = newPos;
4397 badNb += _simplices.size() - nbOkAfter;
4404 //================================================================================
4406 * \brief Add a new segment to _LayerEdge during inflation
4408 //================================================================================
4410 void _LayerEdge::SetNewLength( double len, SMESH_MesherHelper& helper )
4412 if ( _len - len > -1e-6 )
4414 _pos.push_back( _pos.back() );
4418 SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
4419 SMESH_TNodeXYZ oldXYZ( n );
4420 gp_XYZ nXYZ = oldXYZ + _normal * ( len - _len ) * _lenFactor;
4421 n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
4423 _pos.push_back( nXYZ );
4425 if ( !_sWOL.IsNull() )
4428 if ( _sWOL.ShapeType() == TopAbs_EDGE )
4430 double u = Precision::Infinite(); // to force projection w/o distance check
4431 helper.CheckNodeU( TopoDS::Edge( _sWOL ), n, u, 1e-10, /*force=*/true, distXYZ );
4432 _pos.back().SetCoord( u, 0, 0 );
4433 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
4434 pos->SetUParameter( u );
4438 gp_XY uv( Precision::Infinite(), 0 );
4439 helper.CheckNodeUV( TopoDS::Face( _sWOL ), n, uv, 1e-10, /*force=*/true, distXYZ );
4440 _pos.back().SetCoord( uv.X(), uv.Y(), 0 );
4441 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
4442 pos->SetUParameter( uv.X() );
4443 pos->SetVParameter( uv.Y() );
4445 n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]);
4447 dumpMove( n ); //debug
4450 //================================================================================
4452 * \brief Remove last inflation step
4454 //================================================================================
4456 void _LayerEdge::InvalidateStep( int curStep, bool restoreLength )
4458 if ( _pos.size() > curStep )
4460 if ( restoreLength )
4461 _len -= ( _pos[ curStep-1 ] - _pos.back() ).Modulus();
4463 _pos.resize( curStep );
4464 gp_Pnt nXYZ = _pos.back();
4465 SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() );
4466 if ( !_sWOL.IsNull() )
4468 TopLoc_Location loc;
4469 if ( _sWOL.ShapeType() == TopAbs_EDGE )
4471 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( n->GetPosition() );
4472 pos->SetUParameter( nXYZ.X() );
4474 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( _sWOL ), loc, f,l);
4475 nXYZ = curve->Value( nXYZ.X() ).Transformed( loc );
4479 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( n->GetPosition() );
4480 pos->SetUParameter( nXYZ.X() );
4481 pos->SetVParameter( nXYZ.Y() );
4482 Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc );
4483 nXYZ = surface->Value( nXYZ.X(), nXYZ.Y() ).Transformed( loc );
4486 n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() );
4491 //================================================================================
4493 * \brief Create layers of prisms
4495 //================================================================================
4497 bool _ViscousBuilder::refine(_SolidData& data)
4499 SMESH_MesherHelper helper( *_mesh );
4500 helper.SetSubShape( data._solid );
4501 helper.SetElementsOnShape(false);
4503 Handle(Geom_Curve) curve;
4504 Handle(Geom_Surface) surface;
4505 TopoDS_Edge geomEdge;
4506 TopoDS_Face geomFace;
4507 TopoDS_Shape prevSWOL;
4508 TopLoc_Location loc;
4512 TGeomID prevBaseId = -1;
4513 TNode2Edge* n2eMap = 0;
4514 TNode2Edge::iterator n2e;
4516 for ( size_t i = 0; i < data._edges.size(); ++i )
4518 _LayerEdge& edge = *data._edges[i];
4520 // get accumulated length of segments
4521 vector< double > segLen( edge._pos.size() );
4523 for ( size_t j = 1; j < edge._pos.size(); ++j )
4524 segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus();
4526 // allocate memory for new nodes if it is not yet refined
4527 const SMDS_MeshNode* tgtNode = edge._nodes.back();
4528 if ( edge._nodes.size() == 2 )
4530 edge._nodes.resize( data._hyp->GetNumberLayers() + 1, 0 );
4532 edge._nodes.back() = tgtNode;
4534 // get data of a shrink shape
4535 if ( !edge._sWOL.IsNull() && edge._sWOL != prevSWOL )
4537 isOnEdge = ( edge._sWOL.ShapeType() == TopAbs_EDGE );
4540 geomEdge = TopoDS::Edge( edge._sWOL );
4541 curve = BRep_Tool::Curve( geomEdge, loc, f,l);
4545 geomFace = TopoDS::Face( edge._sWOL );
4546 surface = BRep_Tool::Surface( geomFace, loc );
4548 prevSWOL = edge._sWOL;
4550 // restore shapePos of the last node by already treated _LayerEdge of another _SolidData
4551 const TGeomID baseShapeId = edge._nodes[0]->getshapeId();
4552 if ( baseShapeId != prevBaseId )
4554 map< TGeomID, TNode2Edge* >::iterator s2ne = data._s2neMap.find( baseShapeId );
4555 n2eMap = ( s2ne == data._s2neMap.end() ) ? 0 : n2eMap = s2ne->second;
4556 prevBaseId = baseShapeId;
4558 if ( n2eMap && (( n2e = n2eMap->find( edge._nodes[0] )) != n2eMap->end() ))
4560 _LayerEdge* foundEdge = n2e->second;
4561 const gp_XYZ& foundPos = foundEdge->_pos.back();
4562 SMDS_PositionPtr lastPos = tgtNode->GetPosition();
4565 SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( lastPos );
4566 epos->SetUParameter( foundPos.X() );
4570 SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( lastPos );
4571 fpos->SetUParameter( foundPos.X() );
4572 fpos->SetVParameter( foundPos.Y() );
4575 // calculate height of the first layer
4577 const double T = segLen.back(); //data._hyp.GetTotalThickness();
4578 const double f = data._hyp->GetStretchFactor();
4579 const int N = data._hyp->GetNumberLayers();
4580 const double fPowN = pow( f, N );
4581 if ( fPowN - 1 <= numeric_limits<double>::min() )
4584 h0 = T * ( f - 1 )/( fPowN - 1 );
4586 const double zeroLen = std::numeric_limits<double>::min();
4588 // create intermediate nodes
4589 double hSum = 0, hi = h0/f;
4591 for ( size_t iStep = 1; iStep < edge._nodes.size(); ++iStep )
4593 // compute an intermediate position
4596 while ( hSum > segLen[iSeg] && iSeg < segLen.size()-1)
4598 int iPrevSeg = iSeg-1;
4599 while ( fabs( segLen[iPrevSeg] - segLen[iSeg]) <= zeroLen && iPrevSeg > 0 )
4601 double r = ( segLen[iSeg] - hSum ) / ( segLen[iSeg] - segLen[iPrevSeg] );
4602 gp_Pnt pos = r * edge._pos[iPrevSeg] + (1-r) * edge._pos[iSeg];
4604 SMDS_MeshNode*& node = const_cast< SMDS_MeshNode*& >(edge._nodes[ iStep ]);
4605 if ( !edge._sWOL.IsNull() )
4607 // compute XYZ by parameters <pos>
4612 pos = curve->Value( u ).Transformed(loc);
4616 uv.SetCoord( pos.X(), pos.Y() );
4618 pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc);
4621 // create or update the node
4624 node = helper.AddNode( pos.X(), pos.Y(), pos.Z());
4625 if ( !edge._sWOL.IsNull() )
4628 getMeshDS()->SetNodeOnEdge( node, geomEdge, u );
4630 getMeshDS()->SetNodeOnFace( node, geomFace, uv.X(), uv.Y() );
4634 getMeshDS()->SetNodeInVolume( node, helper.GetSubShapeID() );
4639 if ( !edge._sWOL.IsNull() )
4641 // make average pos from new and current parameters
4644 u = 0.5 * ( u + helper.GetNodeU( geomEdge, node ));
4645 pos = curve->Value( u ).Transformed(loc);
4647 SMDS_EdgePosition* epos = static_cast<SMDS_EdgePosition*>( node->GetPosition() );
4648 epos->SetUParameter( u );
4652 uv = 0.5 * ( uv + helper.GetNodeUV( geomFace, node ));
4653 pos = surface->Value( uv.X(), uv.Y()).Transformed(loc);
4655 SMDS_FacePosition* fpos = static_cast<SMDS_FacePosition*>( node->GetPosition() );
4656 fpos->SetUParameter( uv.X() );
4657 fpos->SetVParameter( uv.Y() );
4660 node->setXYZ( pos.X(), pos.Y(), pos.Z() );
4665 if ( !getMeshDS()->IsEmbeddedMode() )
4666 // Log node movement
4667 for ( size_t i = 0; i < data._edges.size(); ++i )
4669 _LayerEdge& edge = *data._edges[i];
4670 SMESH_TNodeXYZ p ( edge._nodes.back() );
4671 getMeshDS()->MoveNode( p._node, p.X(), p.Y(), p.Z() );
4674 // TODO: make quadratic prisms and polyhedrons(?)
4676 helper.SetElementsOnShape(true);
4678 TopExp_Explorer exp( data._solid, TopAbs_FACE );
4679 for ( ; exp.More(); exp.Next() )
4681 if ( data._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( exp.Current() )))
4683 SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() );
4684 SMDS_ElemIteratorPtr fIt = fSubM->GetElements();
4685 vector< vector<const SMDS_MeshNode*>* > nnVec;
4686 while ( fIt->more() )
4688 const SMDS_MeshElement* face = fIt->next();
4689 int nbNodes = face->NbCornerNodes();
4690 nnVec.resize( nbNodes );
4691 SMDS_ElemIteratorPtr nIt = face->nodesIterator();
4692 for ( int iN = 0; iN < nbNodes; ++iN )
4694 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
4695 nnVec[ iN ] = & data._n2eMap[ n ]->_nodes;
4698 int nbZ = nnVec[0]->size();
4702 for ( int iZ = 1; iZ < nbZ; ++iZ )
4703 helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1],
4704 (*nnVec[0])[iZ], (*nnVec[1])[iZ], (*nnVec[2])[iZ]);
4707 for ( int iZ = 1; iZ < nbZ; ++iZ )
4708 helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1],
4709 (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1],
4710 (*nnVec[0])[iZ], (*nnVec[1])[iZ],
4711 (*nnVec[2])[iZ], (*nnVec[3])[iZ]);
4714 return error("Not supported type of element", data._index);
4721 //================================================================================
4723 * \brief Shrink 2D mesh on faces to let space for inflated layers
4725 //================================================================================
4727 bool _ViscousBuilder::shrink()
4729 // make map of (ids of FACEs to shrink mesh on) to (_SolidData containing _LayerEdge's
4730 // inflated along FACE or EDGE)
4731 map< TGeomID, _SolidData* > f2sdMap;
4732 for ( size_t i = 0 ; i < _sdVec.size(); ++i )
4734 _SolidData& data = _sdVec[i];
4735 TopTools_MapOfShape FFMap;
4736 map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin();
4737 for (; s2s != data._shrinkShape2Shape.end(); ++s2s )
4738 if ( s2s->second.ShapeType() == TopAbs_FACE )
4740 f2sdMap.insert( make_pair( getMeshDS()->ShapeToIndex( s2s->second ), &data ));
4742 if ( FFMap.Add( (*s2s).second ))
4743 // Put mesh faces on the shrinked FACE to the proxy sub-mesh to avoid
4744 // usage of mesh faces made in addBoundaryElements() by the 3D algo or
4745 // by StdMeshers_QuadToTriaAdaptor
4746 if ( SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( s2s->second ))
4748 SMESH_ProxyMesh::SubMesh* proxySub =
4749 data._proxyMesh->getFaceSubM( TopoDS::Face( s2s->second ), /*create=*/true);
4750 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
4751 while ( fIt->more() )
4752 proxySub->AddElement( fIt->next() );
4753 // as a result 3D algo will use elements from proxySub and not from smDS
4758 SMESH_MesherHelper helper( *_mesh );
4759 helper.ToFixNodeParameters( true );
4762 map< TGeomID, _Shrinker1D > e2shrMap;
4764 // loop on FACES to srink mesh on
4765 map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin();
4766 for ( ; f2sd != f2sdMap.end(); ++f2sd )
4768 _SolidData& data = *f2sd->second;
4769 TNode2Edge& n2eMap = data._n2eMap;
4770 const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first ));
4772 Handle(Geom_Surface) surface = BRep_Tool::Surface(F);
4774 SMESH_subMesh* sm = _mesh->GetSubMesh( F );
4775 SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
4777 helper.SetSubShape(F);
4779 // ===========================
4780 // Prepare data for shrinking
4781 // ===========================
4783 // Collect nodes to smooth, as src nodes are not yet replaced by tgt ones
4784 // and thus all nodes on a FACE connected to 2d elements are to be smoothed
4785 vector < const SMDS_MeshNode* > smoothNodes;
4787 SMDS_NodeIteratorPtr nIt = smDS->GetNodes();
4788 while ( nIt->more() )
4790 const SMDS_MeshNode* n = nIt->next();
4791 if ( n->NbInverseElements( SMDSAbs_Face ) > 0 )
4792 smoothNodes.push_back( n );
4795 // Find out face orientation
4797 const set<TGeomID> ignoreShapes;
4799 if ( !smoothNodes.empty() )
4801 vector<_Simplex> simplices;
4802 getSimplices( smoothNodes[0], simplices, ignoreShapes );
4803 helper.GetNodeUV( F, simplices[0]._nPrev, 0, &isOkUV ); // fix UV of silpmex nodes
4804 helper.GetNodeUV( F, simplices[0]._nNext, 0, &isOkUV );
4805 gp_XY uv = helper.GetNodeUV( F, smoothNodes[0], 0, &isOkUV );
4806 if ( !simplices[0].IsForward(uv, smoothNodes[0], F, helper,refSign) )
4810 // Find _LayerEdge's inflated along F
4811 vector< _LayerEdge* > lEdges;
4813 SMESH_subMeshIteratorPtr subIt =
4814 sm->getDependsOnIterator(/*includeSelf=*/false, /*complexShapeFirst=*/false);
4815 while ( subIt->more() )
4817 SMESH_subMesh* sub = subIt->next();
4818 SMESHDS_SubMesh* subDS = sub->GetSubMeshDS();
4819 if ( subDS->NbNodes() == 0 || !n2eMap.count( subDS->GetNodes()->next() ))
4821 SMDS_NodeIteratorPtr nIt = subDS->GetNodes();
4822 while ( nIt->more() )
4824 _LayerEdge* edge = n2eMap[ nIt->next() ];
4825 lEdges.push_back( edge );
4826 prepareEdgeToShrink( *edge, F, helper, smDS );
4831 dumpFunction(SMESH_Comment("beforeShrinkFace")<<f2sd->first); // debug
4832 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
4833 while ( fIt->more() )
4834 if ( const SMDS_MeshElement* f = fIt->next() )
4835 dumpChangeNodes( f );
4837 // Replace source nodes by target nodes in mesh faces to shrink
4838 dumpFunction(SMESH_Comment("replNodesOnFace")<<f2sd->first); // debug
4839 const SMDS_MeshNode* nodes[20];
4840 for ( size_t i = 0; i < lEdges.size(); ++i )
4842 _LayerEdge& edge = *lEdges[i];
4843 const SMDS_MeshNode* srcNode = edge._nodes[0];
4844 const SMDS_MeshNode* tgtNode = edge._nodes.back();
4845 SMDS_ElemIteratorPtr fIt = srcNode->GetInverseElementIterator(SMDSAbs_Face);
4846 while ( fIt->more() )
4848 const SMDS_MeshElement* f = fIt->next();
4849 if ( !smDS->Contains( f ))
4851 SMDS_NodeIteratorPtr nIt = f->nodeIterator();
4852 for ( int iN = 0; nIt->more(); ++iN )
4854 const SMDS_MeshNode* n = nIt->next();
4855 nodes[iN] = ( n == srcNode ? tgtNode : n );
4857 helper.GetMeshDS()->ChangeElementNodes( f, nodes, f->NbNodes() );
4858 dumpChangeNodes( f );
4862 // find out if a FACE is concave
4863 const bool isConcaveFace = isConcave( F, helper );
4865 // Create _SmoothNode's on face F
4866 vector< _SmoothNode > nodesToSmooth( smoothNodes.size() );
4868 dumpFunction(SMESH_Comment("fixUVOnFace")<<f2sd->first); // debug
4869 const bool sortSimplices = isConcaveFace;
4870 for ( size_t i = 0; i < smoothNodes.size(); ++i )
4872 const SMDS_MeshNode* n = smoothNodes[i];
4873 nodesToSmooth[ i ]._node = n;
4874 // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices
4875 getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, sortSimplices );
4876 // fix up incorrect uv of nodes on the FACE
4877 helper.GetNodeUV( F, n, 0, &isOkUV);
4881 //if ( nodesToSmooth.empty() ) continue;
4883 // Find EDGE's to shrink and set simpices to LayerEdge's
4884 set< _Shrinker1D* > eShri1D;
4886 for ( size_t i = 0; i < lEdges.size(); ++i )
4888 _LayerEdge* edge = lEdges[i];
4889 if ( edge->_sWOL.ShapeType() == TopAbs_EDGE )
4891 TGeomID edgeIndex = getMeshDS()->ShapeToIndex( edge->_sWOL );
4892 _Shrinker1D& srinker = e2shrMap[ edgeIndex ];
4893 eShri1D.insert( & srinker );
4894 srinker.AddEdge( edge, helper );
4895 VISCOUS_3D::ToClearSubWithMain( _mesh->GetSubMesh( edge->_sWOL ), data._solid );
4896 // restore params of nodes on EGDE if the EDGE has been already
4897 // srinked while srinking another FACE
4898 srinker.RestoreParams();
4900 getSimplices( /*tgtNode=*/edge->_nodes.back(), edge->_simplices, ignoreShapes );
4904 bool toFixTria = false; // to improve quality of trias by diagonal swap
4905 if ( isConcaveFace )
4907 const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles();
4908 if ( hasTria != hasQuad ) {
4909 toFixTria = hasTria;
4912 set<int> nbNodesSet;
4913 SMDS_ElemIteratorPtr fIt = smDS->GetElements();
4914 while ( fIt->more() && nbNodesSet.size() < 2 )
4915 nbNodesSet.insert( fIt->next()->NbCornerNodes() );
4916 toFixTria = ( *nbNodesSet.begin() == 3 );
4920 // ==================
4921 // Perform shrinking
4922 // ==================
4924 bool shrinked = true;
4925 int badNb, shriStep=0, smooStep=0;
4926 _SmoothNode::SmoothType smoothType
4927 = isConcaveFace ? _SmoothNode::ANGULAR : _SmoothNode::LAPLACIAN;
4931 // Move boundary nodes (actually just set new UV)
4932 // -----------------------------------------------
4933 dumpFunction(SMESH_Comment("moveBoundaryOnF")<<f2sd->first<<"_st"<<shriStep ); // debug
4935 for ( size_t i = 0; i < lEdges.size(); ++i )
4937 shrinked |= lEdges[i]->SetNewLength2d( surface,F,helper );
4941 // Move nodes on EDGE's
4942 // (XYZ is set as soon as a needed length reached in SetNewLength2d())
4943 set< _Shrinker1D* >::iterator shr = eShri1D.begin();
4944 for ( ; shr != eShri1D.end(); ++shr )
4945 (*shr)->Compute( /*set3D=*/false, helper );
4948 // -----------------
4949 int nbNoImpSteps = 0;
4952 while (( nbNoImpSteps < 5 && badNb > 0) && moved)
4954 dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
4956 int oldBadNb = badNb;
4959 for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
4961 moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
4962 smoothType, /*set3D=*/isConcaveFace);
4964 if ( badNb < oldBadNb )
4972 return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first );
4973 if ( shriStep > 200 )
4974 return error(SMESH_Comment("Infinite loop at shrinking 2D mesh on face ") << f2sd->first );
4976 // Fix narrow triangles by swapping diagonals
4977 // ---------------------------------------
4980 set<const SMDS_MeshNode*> usedNodes;
4981 fixBadFaces( F, helper, /*is2D=*/true, shriStep, & usedNodes); // swap diagonals
4983 // update working data
4984 set<const SMDS_MeshNode*>::iterator n;
4985 for ( size_t i = 0; i < nodesToSmooth.size() && !usedNodes.empty(); ++i )
4987 n = usedNodes.find( nodesToSmooth[ i ]._node );
4988 if ( n != usedNodes.end())
4990 getSimplices( nodesToSmooth[ i ]._node,
4991 nodesToSmooth[ i ]._simplices,
4993 /*sortSimplices=*/ smoothType == _SmoothNode::ANGULAR );
4994 usedNodes.erase( n );
4997 for ( size_t i = 0; i < lEdges.size() && !usedNodes.empty(); ++i )
4999 n = usedNodes.find( /*tgtNode=*/ lEdges[i]->_nodes.back() );
5000 if ( n != usedNodes.end())
5002 getSimplices( lEdges[i]->_nodes.back(),
5003 lEdges[i]->_simplices,
5005 usedNodes.erase( n );
5009 // TODO: check effect of this additional smooth
5010 // additional laplacian smooth to increase allowed shrink step
5011 // for ( int st = 1; st; --st )
5013 // dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
5014 // for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5016 // nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
5017 // _SmoothNode::LAPLACIAN,/*set3D=*/false);
5020 } // while ( shrinked )
5022 // No wrongly shaped faces remain; final smooth. Set node XYZ.
5023 bool isStructuredFixed = false;
5024 if ( SMESH_2D_Algo* algo = dynamic_cast<SMESH_2D_Algo*>( sm->GetAlgo() ))
5025 isStructuredFixed = algo->FixInternalNodes( *data._proxyMesh, F );
5026 if ( !isStructuredFixed )
5028 if ( isConcaveFace ) // fix narrow faces by swapping diagonals
5029 fixBadFaces( F, helper, /*is2D=*/false, ++shriStep );
5031 for ( int st = 3; st; --st )
5034 case 1: smoothType = _SmoothNode::LAPLACIAN; break;
5035 case 2: smoothType = _SmoothNode::LAPLACIAN; break;
5036 case 3: smoothType = _SmoothNode::ANGULAR; break;
5038 dumpFunction(SMESH_Comment("shrinkFace")<<f2sd->first<<"_st"<<++smooStep); // debug
5039 for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5041 nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,
5042 smoothType,/*set3D=*/st==1 );
5047 // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh
5048 VISCOUS_3D::ToClearSubWithMain( sm, data._solid );
5050 if ( !getMeshDS()->IsEmbeddedMode() )
5051 // Log node movement
5052 for ( size_t i = 0; i < nodesToSmooth.size(); ++i )
5054 SMESH_TNodeXYZ p ( nodesToSmooth[i]._node );
5055 getMeshDS()->MoveNode( nodesToSmooth[i]._node, p.X(), p.Y(), p.Z() );
5058 } // loop on FACES to srink mesh on
5061 // Replace source nodes by target nodes in shrinked mesh edges
5063 map< int, _Shrinker1D >::iterator e2shr = e2shrMap.begin();
5064 for ( ; e2shr != e2shrMap.end(); ++e2shr )
5065 e2shr->second.SwapSrcTgtNodes( getMeshDS() );
5070 //================================================================================
5072 * \brief Computes 2d shrink direction and finds nodes limiting shrinking
5074 //================================================================================
5076 bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge,
5077 const TopoDS_Face& F,
5078 SMESH_MesherHelper& helper,
5079 const SMESHDS_SubMesh* faceSubMesh)
5081 const SMDS_MeshNode* srcNode = edge._nodes[0];
5082 const SMDS_MeshNode* tgtNode = edge._nodes.back();
5086 if ( edge._sWOL.ShapeType() == TopAbs_FACE )
5088 gp_XY srcUV = helper.GetNodeUV( F, srcNode );
5089 gp_XY tgtUV = helper.GetNodeUV( F, tgtNode );
5090 gp_Vec2d uvDir( srcUV, tgtUV );
5091 double uvLen = uvDir.Magnitude();
5093 edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0 );
5096 edge._pos.resize(1);
5097 edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 );
5099 // set UV of source node to target node
5100 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
5101 pos->SetUParameter( srcUV.X() );
5102 pos->SetVParameter( srcUV.Y() );
5104 else // _sWOL is TopAbs_EDGE
5106 TopoDS_Edge E = TopoDS::Edge( edge._sWOL);
5107 SMESHDS_SubMesh* edgeSM = getMeshDS()->MeshElements( E );
5108 if ( !edgeSM || edgeSM->NbElements() == 0 )
5109 return error(SMESH_Comment("Not meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
5111 const SMDS_MeshNode* n2 = 0;
5112 SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
5113 while ( eIt->more() && !n2 )
5115 const SMDS_MeshElement* e = eIt->next();
5116 if ( !edgeSM->Contains(e)) continue;
5117 n2 = e->GetNode( 0 );
5118 if ( n2 == srcNode ) n2 = e->GetNode( 1 );
5121 return error(SMESH_Comment("Wrongly meshed EDGE ") << getMeshDS()->ShapeToIndex( E ));
5123 double uSrc = helper.GetNodeU( E, srcNode, n2 );
5124 double uTgt = helper.GetNodeU( E, tgtNode, srcNode );
5125 double u2 = helper.GetNodeU( E, n2, srcNode );
5127 if ( fabs( uSrc-uTgt ) < 0.99 * fabs( uSrc-u2 ))
5129 // tgtNode is located so that it does not make faces with wrong orientation
5132 edge._pos.resize(1);
5133 edge._pos[0].SetCoord( U_TGT, uTgt );
5134 edge._pos[0].SetCoord( U_SRC, uSrc );
5135 edge._pos[0].SetCoord( LEN_TGT, fabs( uSrc-uTgt ));
5137 edge._simplices.resize( 1 );
5138 edge._simplices[0]._nPrev = n2;
5140 // set UV of source node to target node
5141 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
5142 pos->SetUParameter( uSrc );
5147 //================================================================================
5149 * \brief Try to fix triangles with high aspect ratio by swaping diagonals
5151 //================================================================================
5153 void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F,
5154 SMESH_MesherHelper& helper,
5157 set<const SMDS_MeshNode*> * involvedNodes)
5159 SMESH::Controls::AspectRatio qualifier;
5160 SMESH::Controls::TSequenceOfXYZ points(3), points1(3), points2(3);
5161 const double maxAspectRatio = is2D ? 4. : 2;
5162 _NodeCoordHelper xyz( F, helper, is2D );
5164 // find bad triangles
5166 vector< const SMDS_MeshElement* > badTrias;
5167 vector< double > badAspects;
5168 SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( F );
5169 SMDS_ElemIteratorPtr fIt = sm->GetElements();
5170 while ( fIt->more() )
5172 const SMDS_MeshElement * f = fIt->next();
5173 if ( f->NbCornerNodes() != 3 ) continue;
5174 for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = xyz( f->GetNode(iP));
5175 double aspect = qualifier.GetValue( points );
5176 if ( aspect > maxAspectRatio )
5178 badTrias.push_back( f );
5179 badAspects.push_back( aspect );
5184 dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID());
5185 SMDS_ElemIteratorPtr fIt = sm->GetElements();
5186 while ( fIt->more() )
5188 const SMDS_MeshElement * f = fIt->next();
5189 if ( f->NbCornerNodes() == 3 )
5190 dumpChangeNodes( f );
5194 if ( badTrias.empty() )
5197 // find couples of faces to swap diagonal
5199 typedef pair < const SMDS_MeshElement* , const SMDS_MeshElement* > T2Trias;
5200 vector< T2Trias > triaCouples;
5202 TIDSortedElemSet involvedFaces, emptySet;
5203 for ( size_t iTia = 0; iTia < badTrias.size(); ++iTia )
5206 double aspRatio [3];
5209 if ( !involvedFaces.insert( badTrias[iTia] ).second )
5211 for ( int iP = 0; iP < 3; ++iP )
5212 points(iP+1) = xyz( badTrias[iTia]->GetNode(iP));
5214 // find triangles adjacent to badTrias[iTia] with better aspect ratio after diag-swaping
5215 int bestCouple = -1;
5216 for ( int iSide = 0; iSide < 3; ++iSide )
5218 const SMDS_MeshNode* n1 = badTrias[iTia]->GetNode( iSide );
5219 const SMDS_MeshNode* n2 = badTrias[iTia]->GetNode(( iSide+1 ) % 3 );
5220 trias [iSide].first = badTrias[iTia];
5221 trias [iSide].second = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, involvedFaces,
5223 if (( ! trias[iSide].second ) ||
5224 ( trias[iSide].second->NbCornerNodes() != 3 ) ||
5225 ( ! sm->Contains( trias[iSide].second )))
5228 // aspect ratio of an adjacent tria
5229 for ( int iP = 0; iP < 3; ++iP )
5230 points2(iP+1) = xyz( trias[iSide].second->GetNode(iP));
5231 double aspectInit = qualifier.GetValue( points2 );
5233 // arrange nodes as after diag-swaping
5234 if ( helper.WrapIndex( i1+1, 3 ) == i2 )
5235 i3 = helper.WrapIndex( i1-1, 3 );
5237 i3 = helper.WrapIndex( i1+1, 3 );
5239 points1( 1+ iSide ) = points2( 1+ i3 );
5240 points2( 1+ i2 ) = points1( 1+ ( iSide+2 ) % 3 );
5242 // aspect ratio after diag-swaping
5243 aspRatio[ iSide ] = qualifier.GetValue( points1 ) + qualifier.GetValue( points2 );
5244 if ( aspRatio[ iSide ] > aspectInit + badAspects[ iTia ] )
5247 // prevent inversion of a triangle
5248 gp_Vec norm1 = gp_Vec( points1(1), points1(3) ) ^ gp_Vec( points1(1), points1(2) );
5249 gp_Vec norm2 = gp_Vec( points2(1), points2(3) ) ^ gp_Vec( points2(1), points2(2) );
5250 if ( norm1 * norm2 < 0. && norm1.Angle( norm2 ) > 70./180.*M_PI )
5253 if ( bestCouple < 0 || aspRatio[ bestCouple ] > aspRatio[ iSide ] )
5257 if ( bestCouple >= 0 )
5259 triaCouples.push_back( trias[bestCouple] );
5260 involvedFaces.insert ( trias[bestCouple].second );
5264 involvedFaces.erase( badTrias[iTia] );
5267 if ( triaCouples.empty() )
5272 SMESH_MeshEditor editor( helper.GetMesh() );
5273 dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
5274 for ( size_t i = 0; i < triaCouples.size(); ++i )
5276 dumpChangeNodes( triaCouples[i].first );
5277 dumpChangeNodes( triaCouples[i].second );
5278 editor.InverseDiag( triaCouples[i].first, triaCouples[i].second );
5281 if ( involvedNodes )
5282 for ( size_t i = 0; i < triaCouples.size(); ++i )
5284 involvedNodes->insert( triaCouples[i].first->begin_nodes(),
5285 triaCouples[i].first->end_nodes() );
5286 involvedNodes->insert( triaCouples[i].second->begin_nodes(),
5287 triaCouples[i].second->end_nodes() );
5290 // just for debug dump resulting triangles
5291 dumpFunction(SMESH_Comment("swapDiagonals_F")<<helper.GetSubShapeID()<<"_"<<step);
5292 for ( size_t i = 0; i < triaCouples.size(); ++i )
5294 dumpChangeNodes( triaCouples[i].first );
5295 dumpChangeNodes( triaCouples[i].second );
5299 //================================================================================
5301 * \brief Move target node to it's final position on the FACE during shrinking
5303 //================================================================================
5305 bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface,
5306 const TopoDS_Face& F,
5307 SMESH_MesherHelper& helper )
5310 return false; // already at the target position
5312 SMDS_MeshNode* tgtNode = const_cast< SMDS_MeshNode*& >( _nodes.back() );
5314 if ( _sWOL.ShapeType() == TopAbs_FACE )
5316 gp_XY curUV = helper.GetNodeUV( F, tgtNode );
5317 gp_Pnt2d tgtUV( _pos[0].X(), _pos[0].Y() );
5318 gp_Vec2d uvDir( _normal.X(), _normal.Y() );
5319 const double uvLen = tgtUV.Distance( curUV );
5320 const double kSafe = Max( 0.5, 1. - 0.1 * _simplices.size() );
5322 // Select shrinking step such that not to make faces with wrong orientation.
5323 double stepSize = 1e100;
5324 for ( size_t i = 0; i < _simplices.size(); ++i )
5326 // find intersection of 2 lines: curUV-tgtUV and that connecting simplex nodes
5327 gp_XY uvN1 = helper.GetNodeUV( F, _simplices[i]._nPrev );
5328 gp_XY uvN2 = helper.GetNodeUV( F, _simplices[i]._nNext );
5329 gp_XY dirN = uvN2 - uvN1;
5330 double det = uvDir.Crossed( dirN );
5331 if ( Abs( det ) < std::numeric_limits<double>::min() ) continue;
5332 gp_XY dirN2Cur = curUV - uvN1;
5333 double step = dirN.Crossed( dirN2Cur ) / det;
5335 stepSize = Min( step, stepSize );
5338 if ( uvLen <= stepSize )
5343 else if ( stepSize > 0 )
5345 newUV = curUV + uvDir.XY() * stepSize * kSafe;
5351 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( tgtNode->GetPosition() );
5352 pos->SetUParameter( newUV.X() );
5353 pos->SetVParameter( newUV.Y() );
5356 gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
5357 tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
5358 dumpMove( tgtNode );
5361 else // _sWOL is TopAbs_EDGE
5363 TopoDS_Edge E = TopoDS::Edge( _sWOL );
5364 const SMDS_MeshNode* n2 = _simplices[0]._nPrev;
5365 SMDS_EdgePosition* tgtPos = static_cast<SMDS_EdgePosition*>( tgtNode->GetPosition() );
5367 const double u2 = helper.GetNodeU( E, n2, tgtNode );
5368 const double uSrc = _pos[0].Coord( U_SRC );
5369 const double lenTgt = _pos[0].Coord( LEN_TGT );
5371 double newU = _pos[0].Coord( U_TGT );
5372 if ( lenTgt < 0.99 * fabs( uSrc-u2 )) // n2 got out of src-tgt range
5378 newU = 0.1 * tgtPos->GetUParameter() + 0.9 * u2;
5380 tgtPos->SetUParameter( newU );
5382 gp_XY newUV = helper.GetNodeUV( F, tgtNode, _nodes[0]);
5383 gp_Pnt p = surface->Value( newUV.X(), newUV.Y() );
5384 tgtNode->setXYZ( p.X(), p.Y(), p.Z() );
5385 dumpMove( tgtNode );
5391 //================================================================================
5393 * \brief Perform smooth on the FACE
5394 * \retval bool - true if the node has been moved
5396 //================================================================================
5398 bool _SmoothNode::Smooth(int& badNb,
5399 Handle(Geom_Surface)& surface,
5400 SMESH_MesherHelper& helper,
5401 const double refSign,
5405 const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() );
5407 // get uv of surrounding nodes
5408 vector<gp_XY> uv( _simplices.size() );
5409 for ( size_t i = 0; i < _simplices.size(); ++i )
5410 uv[i] = helper.GetNodeUV( face, _simplices[i]._nPrev, _node );
5412 // compute new UV for the node
5414 if ( how == TFI && _simplices.size() == 4 )
5417 for ( size_t i = 0; i < _simplices.size(); ++i )
5418 if ( _simplices[i]._nOpp )
5419 corners[i] = helper.GetNodeUV( face, _simplices[i]._nOpp, _node );
5421 throw SALOME_Exception(LOCALIZED("TFI smoothing: _Simplex::_nOpp not set!"));
5423 newPos = helper.calcTFI ( 0.5, 0.5,
5424 corners[0], corners[1], corners[2], corners[3],
5425 uv[1], uv[2], uv[3], uv[0] );
5427 else if ( how == ANGULAR )
5429 newPos = computeAngularPos( uv, helper.GetNodeUV( face, _node ), refSign );
5431 else if ( how == CENTROIDAL && _simplices.size() > 3 )
5433 // average centers of diagonals wieghted with their reciprocal lengths
5434 if ( _simplices.size() == 4 )
5436 double w1 = 1. / ( uv[2]-uv[0] ).SquareModulus();
5437 double w2 = 1. / ( uv[3]-uv[1] ).SquareModulus();
5438 newPos = ( w1 * ( uv[2]+uv[0] ) + w2 * ( uv[3]+uv[1] )) / ( w1+w2 ) / 2;
5442 double sumWeight = 0;
5443 int nb = _simplices.size() == 4 ? 2 : _simplices.size();
5444 for ( int i = 0; i < nb; ++i )
5447 int iTo = i + _simplices.size() - 1;
5448 for ( int j = iFrom; j < iTo; ++j )
5450 int i2 = SMESH_MesherHelper::WrapIndex( j, _simplices.size() );
5451 double w = 1. / ( uv[i]-uv[i2] ).SquareModulus();
5453 newPos += w * ( uv[i]+uv[i2] );
5456 newPos /= 2 * sumWeight; // 2 is to get a middle between uv's
5462 for ( size_t i = 0; i < _simplices.size(); ++i )
5464 newPos /= _simplices.size();
5467 // count quality metrics (orientation) of triangles around the node
5469 gp_XY tgtUV = helper.GetNodeUV( face, _node );
5470 for ( size_t i = 0; i < _simplices.size(); ++i )
5471 nbOkBefore += _simplices[i].IsForward( tgtUV, _node, face, helper, refSign );
5474 for ( size_t i = 0; i < _simplices.size(); ++i )
5475 nbOkAfter += _simplices[i].IsForward( newPos, _node, face, helper, refSign );
5477 if ( nbOkAfter < nbOkBefore )
5479 badNb += _simplices.size() - nbOkBefore;
5483 SMDS_FacePosition* pos = static_cast<SMDS_FacePosition*>( _node->GetPosition() );
5484 pos->SetUParameter( newPos.X() );
5485 pos->SetVParameter( newPos.Y() );
5492 gp_Pnt p = surface->Value( newPos.X(), newPos.Y() );
5493 const_cast< SMDS_MeshNode* >( _node )->setXYZ( p.X(), p.Y(), p.Z() );
5497 badNb += _simplices.size() - nbOkAfter;
5498 return ( (tgtUV-newPos).SquareModulus() > 1e-10 );
5501 //================================================================================
5503 * \brief Computes new UV using angle based smoothing technic
5505 //================================================================================
5507 gp_XY _SmoothNode::computeAngularPos(vector<gp_XY>& uv,
5508 const gp_XY& uvToFix,
5509 const double refSign)
5511 uv.push_back( uv.front() );
5513 vector< gp_XY > edgeDir ( uv.size() );
5514 vector< double > edgeSize( uv.size() );
5515 for ( size_t i = 1; i < edgeDir.size(); ++i )
5517 edgeDir [i-1] = uv[i] - uv[i-1];
5518 edgeSize[i-1] = edgeDir[i-1].Modulus();
5519 if ( edgeSize[i-1] < numeric_limits<double>::min() )
5520 edgeDir[i-1].SetX( 100 );
5522 edgeDir[i-1] /= edgeSize[i-1] * refSign;
5524 edgeDir.back() = edgeDir.front();
5525 edgeSize.back() = edgeSize.front();
5530 for ( size_t i = 1; i < edgeDir.size(); ++i )
5532 if ( edgeDir[i-1].X() > 1. ) continue;
5534 while ( edgeDir[i].X() > 1. && ++i < edgeDir.size() );
5535 if ( i == edgeDir.size() ) break;
5537 gp_XY norm1( -edgeDir[i1].Y(), edgeDir[i1].X() );
5538 gp_XY norm2( -edgeDir[i].Y(), edgeDir[i].X() );
5539 gp_XY bisec = norm1 + norm2;
5540 double bisecSize = bisec.Modulus();
5541 if ( bisecSize < numeric_limits<double>::min() )
5543 bisec = -edgeDir[i1] + edgeDir[i];
5544 bisecSize = bisec.Modulus();
5548 gp_XY dirToN = uvToFix - p;
5549 double distToN = dirToN.Modulus();
5550 if ( bisec * dirToN < 0 )
5553 newPos += ( p + bisec * distToN ) * ( edgeSize[i1] + edgeSize[i] );
5555 sumSize += edgeSize[i1] + edgeSize[i];
5557 newPos /= /*nbEdges * */sumSize;
5561 //================================================================================
5563 * \brief Delete _SolidData
5565 //================================================================================
5567 _SolidData::~_SolidData()
5569 for ( size_t i = 0; i < _edges.size(); ++i )
5571 if ( _edges[i] && _edges[i]->_2neibors )
5572 delete _edges[i]->_2neibors;
5577 //================================================================================
5579 * \brief Add a _LayerEdge inflated along the EDGE
5581 //================================================================================
5583 void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper )
5586 if ( _nodes.empty() )
5588 _edges[0] = _edges[1] = 0;
5592 if ( e == _edges[0] || e == _edges[1] )
5594 if ( e->_sWOL.IsNull() || e->_sWOL.ShapeType() != TopAbs_EDGE )
5595 throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
5596 if ( _edges[0] && _edges[0]->_sWOL != e->_sWOL )
5597 throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added"));
5600 const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
5602 BRep_Tool::Range( E, f,l );
5603 double u = helper.GetNodeU( E, e->_nodes[0], e->_nodes.back());
5604 _edges[ u < 0.5*(f+l) ? 0 : 1 ] = e;
5608 const SMDS_MeshNode* tgtNode0 = _edges[0] ? _edges[0]->_nodes.back() : 0;
5609 const SMDS_MeshNode* tgtNode1 = _edges[1] ? _edges[1]->_nodes.back() : 0;
5611 if ( _nodes.empty() )
5613 SMESHDS_SubMesh * eSubMesh = helper.GetMeshDS()->MeshElements( E );
5614 if ( !eSubMesh || eSubMesh->NbNodes() < 1 )
5616 TopLoc_Location loc;
5617 Handle(Geom_Curve) C = BRep_Tool::Curve(E, loc, f,l);
5618 GeomAdaptor_Curve aCurve(C, f,l);
5619 const double totLen = GCPnts_AbscissaPoint::Length(aCurve, f, l);
5621 int nbExpectNodes = eSubMesh->NbNodes();
5622 _initU .reserve( nbExpectNodes );
5623 _normPar.reserve( nbExpectNodes );
5624 _nodes .reserve( nbExpectNodes );
5625 SMDS_NodeIteratorPtr nIt = eSubMesh->GetNodes();
5626 while ( nIt->more() )
5628 const SMDS_MeshNode* node = nIt->next();
5629 if ( node->NbInverseElements(SMDSAbs_Edge) == 0 ||
5630 node == tgtNode0 || node == tgtNode1 )
5631 continue; // refinement nodes
5632 _nodes.push_back( node );
5633 _initU.push_back( helper.GetNodeU( E, node ));
5634 double len = GCPnts_AbscissaPoint::Length(aCurve, f, _initU.back());
5635 _normPar.push_back( len / totLen );
5640 // remove target node of the _LayerEdge from _nodes
5642 for ( size_t i = 0; i < _nodes.size(); ++i )
5643 if ( !_nodes[i] || _nodes[i] == tgtNode0 || _nodes[i] == tgtNode1 )
5644 _nodes[i] = 0, nbFound++;
5645 if ( nbFound == _nodes.size() )
5650 //================================================================================
5652 * \brief Move nodes on EDGE from ends where _LayerEdge's are inflated
5654 //================================================================================
5656 void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper)
5658 if ( _done || _nodes.empty())
5660 const _LayerEdge* e = _edges[0];
5661 if ( !e ) e = _edges[1];
5664 _done = (( !_edges[0] || _edges[0]->_pos.empty() ) &&
5665 ( !_edges[1] || _edges[1]->_pos.empty() ));
5667 const TopoDS_Edge& E = TopoDS::Edge( e->_sWOL );
5669 if ( set3D || _done )
5671 Handle(Geom_Curve) C = BRep_Tool::Curve(E, f,l);
5672 GeomAdaptor_Curve aCurve(C, f,l);
5675 f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
5677 l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
5678 double totLen = GCPnts_AbscissaPoint::Length( aCurve, f, l );
5680 for ( size_t i = 0; i < _nodes.size(); ++i )
5682 if ( !_nodes[i] ) continue;
5683 double len = totLen * _normPar[i];
5684 GCPnts_AbscissaPoint discret( aCurve, len, f );
5685 if ( !discret.IsDone() )
5686 return throw SALOME_Exception(LOCALIZED("GCPnts_AbscissaPoint failed"));
5687 double u = discret.Parameter();
5688 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
5689 pos->SetUParameter( u );
5690 gp_Pnt p = C->Value( u );
5691 const_cast< SMDS_MeshNode*>( _nodes[i] )->setXYZ( p.X(), p.Y(), p.Z() );
5696 BRep_Tool::Range( E, f,l );
5698 f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] );
5700 l = helper.GetNodeU( E, _edges[1]->_nodes.back(), _nodes.back() );
5702 for ( size_t i = 0; i < _nodes.size(); ++i )
5704 if ( !_nodes[i] ) continue;
5705 double u = f * ( 1-_normPar[i] ) + l * _normPar[i];
5706 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
5707 pos->SetUParameter( u );
5712 //================================================================================
5714 * \brief Restore initial parameters of nodes on EDGE
5716 //================================================================================
5718 void _Shrinker1D::RestoreParams()
5721 for ( size_t i = 0; i < _nodes.size(); ++i )
5723 if ( !_nodes[i] ) continue;
5724 SMDS_EdgePosition* pos = static_cast<SMDS_EdgePosition*>( _nodes[i]->GetPosition() );
5725 pos->SetUParameter( _initU[i] );
5730 //================================================================================
5732 * \brief Replace source nodes by target nodes in shrinked mesh edges
5734 //================================================================================
5736 void _Shrinker1D::SwapSrcTgtNodes( SMESHDS_Mesh* mesh )
5738 const SMDS_MeshNode* nodes[3];
5739 for ( int i = 0; i < 2; ++i )
5741 if ( !_edges[i] ) continue;
5743 SMESHDS_SubMesh * eSubMesh = mesh->MeshElements( _edges[i]->_sWOL );
5744 if ( !eSubMesh ) return;
5745 const SMDS_MeshNode* srcNode = _edges[i]->_nodes[0];
5746 const SMDS_MeshNode* tgtNode = _edges[i]->_nodes.back();
5747 SMDS_ElemIteratorPtr eIt = srcNode->GetInverseElementIterator(SMDSAbs_Edge);
5748 while ( eIt->more() )
5750 const SMDS_MeshElement* e = eIt->next();
5751 if ( !eSubMesh->Contains( e ))
5753 SMDS_ElemIteratorPtr nIt = e->nodesIterator();
5754 for ( int iN = 0; iN < e->NbNodes(); ++iN )
5756 const SMDS_MeshNode* n = static_cast<const SMDS_MeshNode*>( nIt->next() );
5757 nodes[iN] = ( n == srcNode ? tgtNode : n );
5759 mesh->ChangeElementNodes( e, nodes, e->NbNodes() );
5764 //================================================================================
5766 * \brief Creates 2D and 1D elements on boundaries of new prisms
5768 //================================================================================
5770 bool _ViscousBuilder::addBoundaryElements()
5772 SMESH_MesherHelper helper( *_mesh );
5774 for ( size_t i = 0; i < _sdVec.size(); ++i )
5776 _SolidData& data = _sdVec[i];
5777 TopTools_IndexedMapOfShape geomEdges;
5778 TopExp::MapShapes( data._solid, TopAbs_EDGE, geomEdges );
5779 for ( int iE = 1; iE <= geomEdges.Extent(); ++iE )
5781 const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE));
5783 // Get _LayerEdge's based on E
5785 map< double, const SMDS_MeshNode* > u2nodes;
5786 if ( !SMESH_Algo::GetSortedNodesOnEdge( getMeshDS(), E, /*ignoreMedium=*/false, u2nodes))
5789 vector< _LayerEdge* > ledges; ledges.reserve( u2nodes.size() );
5790 TNode2Edge & n2eMap = data._n2eMap;
5791 map< double, const SMDS_MeshNode* >::iterator u2n = u2nodes.begin();
5793 //check if 2D elements are needed on E
5794 TNode2Edge::iterator n2e = n2eMap.find( u2n->second );
5795 if ( n2e == n2eMap.end() ) continue; // no layers on vertex
5796 ledges.push_back( n2e->second );
5798 if (( n2e = n2eMap.find( u2n->second )) == n2eMap.end() )
5799 continue; // no layers on E
5800 ledges.push_back( n2eMap[ u2n->second ]);
5802 const SMDS_MeshNode* tgtN0 = ledges[0]->_nodes.back();
5803 const SMDS_MeshNode* tgtN1 = ledges[1]->_nodes.back();
5804 int nbSharedPyram = 0;
5805 SMDS_ElemIteratorPtr vIt = tgtN0->GetInverseElementIterator(SMDSAbs_Volume);
5806 while ( vIt->more() )
5808 const SMDS_MeshElement* v = vIt->next();
5809 nbSharedPyram += int( v->GetNodeIndex( tgtN1 ) >= 0 );
5811 if ( nbSharedPyram > 1 )
5812 continue; // not free border of the pyramid
5814 if ( getMeshDS()->FindFace( ledges[0]->_nodes[0], ledges[0]->_nodes[1],
5815 ledges[1]->_nodes[0], ledges[1]->_nodes[1]))
5816 continue; // faces already created
5818 for ( ++u2n; u2n != u2nodes.end(); ++u2n )
5819 ledges.push_back( n2eMap[ u2n->second ]);
5821 // Find out orientation and type of face to create
5823 bool reverse = false, isOnFace;
5825 map< TGeomID, TopoDS_Shape >::iterator e2f =
5826 data._shrinkShape2Shape.find( getMeshDS()->ShapeToIndex( E ));
5828 if (( isOnFace = ( e2f != data._shrinkShape2Shape.end() )))
5830 F = e2f->second.Oriented( TopAbs_FORWARD );
5831 reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED );
5832 if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED )
5833 reverse = !reverse, F.Reverse();
5834 if ( helper.IsReversedSubMesh( TopoDS::Face(F) ))
5839 // find FACE with layers sharing E
5840 PShapeIteratorPtr fIt = helper.GetAncestors( E, *_mesh, TopAbs_FACE );
5841 while ( fIt->more() && F.IsNull() )
5843 const TopoDS_Shape* pF = fIt->next();
5844 if ( helper.IsSubShape( *pF, data._solid) &&
5845 !data._ignoreFaceIds.count( e2f->first ))
5849 // Find the sub-mesh to add new faces
5850 SMESHDS_SubMesh* sm = 0;
5852 sm = getMeshDS()->MeshElements( F );
5854 sm = data._proxyMesh->getFaceSubM( TopoDS::Face(F), /*create=*/true );
5856 return error("error in addBoundaryElements()", data._index);
5859 const int dj1 = reverse ? 0 : 1;
5860 const int dj2 = reverse ? 1 : 0;
5861 for ( size_t j = 1; j < ledges.size(); ++j )
5863 vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes;
5864 vector< const SMDS_MeshNode*>& nn2 = ledges[j-dj2]->_nodes;
5866 for ( size_t z = 1; z < nn1.size(); ++z )
5867 sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] ));
5869 for ( size_t z = 1; z < nn1.size(); ++z )
5870 sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z]));
5874 for ( int isFirst = 0; isFirst < 2; ++isFirst )
5876 _LayerEdge* edge = isFirst ? ledges.front() : ledges.back();
5877 if ( !edge->_sWOL.IsNull() && edge->_sWOL.ShapeType() == TopAbs_EDGE )
5879 vector< const SMDS_MeshNode*>& nn = edge->_nodes;
5880 if ( nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() )
5882 helper.SetSubShape( edge->_sWOL );
5883 helper.SetElementsOnShape( true );
5884 for ( size_t z = 1; z < nn.size(); ++z )
5885 helper.AddEdge( nn[z-1], nn[z] );