Salome HOME
Merge from V6_main 13/12/2012
[modules/smesh.git] / src / StdMeshers / StdMeshers_FaceSide.cxx
1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
5 //
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
10 //
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14 // Lesser General Public License for more details.
15 //
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
19 //
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 //
22
23 // File      : StdMeshers_FaceSide.hxx
24 // Created   : Wed Jan 31 18:41:25 2007
25 // Author    : Edward AGAPOV (eap)
26 // Module    : SMESH
27 //
28 #include "StdMeshers_FaceSide.hxx"
29
30 #include "SMDS_EdgePosition.hxx"
31 #include "SMDS_MeshNode.hxx"
32 #include "SMESHDS_Mesh.hxx"
33 #include "SMESHDS_SubMesh.hxx"
34 #include "SMESH_Algo.hxx"
35 #include "SMESH_Mesh.hxx"
36 #include "SMESH_MesherHelper.hxx"
37 #include "SMESH_ComputeError.hxx"
38 #include "SMESH_Block.hxx"
39
40 #include <Adaptor2d_Curve2d.hxx>
41 #include <BRepAdaptor_CompCurve.hxx>
42 #include <BRep_Builder.hxx>
43 #include <BRep_Tool.hxx>
44 #include <GCPnts_AbscissaPoint.hxx>
45 #include <Geom2dAdaptor_Curve.hxx>
46 #include <TopExp.hxx>
47 #include <TopExp_Explorer.hxx>
48 #include <TopoDS.hxx>
49 #include <TopoDS_Face.hxx>
50 #include <TopoDS_Vertex.hxx>
51 #include <TopoDS_Wire.hxx>
52
53 #include <map>
54
55 #include "utilities.h"
56
57 //================================================================================
58 /*!
59  * \brief Constructor of a side of one edge
60   * \param theFace - the face
61   * \param theEdge - the edge
62  */
63 //================================================================================
64
65 StdMeshers_FaceSide::StdMeshers_FaceSide(const TopoDS_Face&   theFace,
66                                          const TopoDS_Edge&   theEdge,
67                                          SMESH_Mesh*          theMesh,
68                                          const bool           theIsForward,
69                                          const bool           theIgnoreMediumNodes,
70                                          SMESH_ProxyMesh::Ptr theProxyMesh)
71 {
72   list<TopoDS_Edge> edges(1,theEdge);
73   *this = StdMeshers_FaceSide( theFace, edges, theMesh, theIsForward,
74                                theIgnoreMediumNodes, theProxyMesh );
75 }
76
77 //================================================================================
78 /*!
79  * \brief Constructor of a side of several edges
80   * \param theFace - the face
81   * \param theEdge - the edge
82  */
83 //================================================================================
84
85 StdMeshers_FaceSide::StdMeshers_FaceSide(const TopoDS_Face&   theFace,
86                                          list<TopoDS_Edge>&   theEdges,
87                                          SMESH_Mesh*          theMesh,
88                                          const bool           theIsForward,
89                                          const bool           theIgnoreMediumNodes,
90                                          SMESH_ProxyMesh::Ptr theProxyMesh)
91 {
92   int nbEdges = theEdges.size();
93   myEdge.resize      ( nbEdges );
94   myEdgeID.resize    ( nbEdges );
95   myC2d.resize       ( nbEdges );
96   myC3dAdaptor.resize( nbEdges );
97   myFirst.resize     ( nbEdges );
98   myLast.resize      ( nbEdges );
99   myNormPar.resize   ( nbEdges );
100   myEdgeLength.resize( nbEdges );
101   myIsUniform.resize ( nbEdges, true );
102   myLength             = 0;
103   myNbPonits           = myNbSegments = 0;
104   myProxyMesh          = theProxyMesh;
105   myMissingVertexNodes = false;
106   myIgnoreMediumNodes  = theIgnoreMediumNodes;
107   myDefaultPnt2d       = gp_Pnt2d( 1e+100, 1e+100 );
108   if ( !myProxyMesh ) myProxyMesh.reset( new SMESH_ProxyMesh( *theMesh ));
109   if ( nbEdges == 0 ) return;
110
111   SMESHDS_Mesh* meshDS = myProxyMesh->GetMeshDS();
112
113   int nbDegen = 0;
114   list<TopoDS_Edge>::iterator edge = theEdges.begin();
115   TopoDS_Iterator vExp;
116   for ( int index = 0; edge != theEdges.end(); ++index, ++edge )
117   {
118     int i = theIsForward ? index : nbEdges-index-1;
119     myEdgeLength[i] = SMESH_Algo::EdgeLength( *edge );
120     if ( myEdgeLength[i] < DBL_MIN ) nbDegen++;
121     myLength += myEdgeLength[i];
122     myEdge[i] = *edge;
123     myEdgeID[i] = meshDS->ShapeToIndex( *edge );
124     if ( !theIsForward ) myEdge[i].Reverse();
125
126     if ( theFace.IsNull() )
127       BRep_Tool::Range( *edge, myFirst[i], myLast[i] );
128     else
129       myC2d[i] = BRep_Tool::CurveOnSurface( *edge, theFace, myFirst[i], myLast[i] );
130     if ( myEdge[i].Orientation() == TopAbs_REVERSED )
131       std::swap( myFirst[i], myLast[i] );
132
133     if ( const SMESHDS_SubMesh* sm = myProxyMesh->GetSubMesh( *edge )) {
134       int nbN = sm->NbNodes();
135       if ( theIgnoreMediumNodes ) {
136         SMDS_ElemIteratorPtr elemIt = sm->GetElements();
137         if ( elemIt->more() && elemIt->next()->IsQuadratic() )
138           nbN -= sm->NbElements();
139       }
140       myNbPonits += nbN;
141       myNbSegments += sm->NbElements();
142     }
143
144     // TopExp::FirstVertex() and TopExp::LastVertex() return NULL from INTERNAL edge
145     vExp.Initialize( *edge );
146     if ( vExp.Value().Orientation() == TopAbs_REVERSED ) vExp.Next();
147     if ( SMESH_Algo::VertexNode( TopoDS::Vertex( vExp.Value()), meshDS ))
148       myNbPonits += 1; // for the first end
149     else
150       myMissingVertexNodes = true;
151
152     // check if the edge has a non-uniform parametrization (issue 0020705)
153     if ( !myC2d[i].IsNull() && myEdgeLength[i] > DBL_MIN)
154     {
155       Geom2dAdaptor_Curve A2dC( myC2d[i],
156                                 std::min( myFirst[i], myLast[i] ),
157                                 std::max( myFirst[i], myLast[i] ));
158       double p2 = myFirst[i]+(myLast[i]-myFirst[i])/2., p4 = myFirst[i]+(myLast[i]-myFirst[i])/4.;
159       double d2 = GCPnts_AbscissaPoint::Length( A2dC, myFirst[i], p2 );
160       double d4 = GCPnts_AbscissaPoint::Length( A2dC, myFirst[i], p4 );
161       //cout<<"len = "<<len<<"  d2 = "<<d2<<"  fabs(2*d2/len-1.0) = "<<fabs(2*d2/len-1.0)<<endl;
162       myIsUniform[i] = !( fabs(2*d2/myEdgeLength[i]-1.0) > 0.01 || fabs(2*d4/d2-1.0) > 0.01 );
163       if ( !myIsUniform[i] )
164       {
165         double fp,lp;
166         TopLoc_Location L;
167         Handle(Geom_Curve) C3d = BRep_Tool::Curve(myEdge[i],L,fp,lp);
168         myC3dAdaptor[i].Load( C3d, fp,lp );
169       }
170     }
171     // reverse a proxy submesh
172     if ( !theIsForward )
173       reverseProxySubmesh( myEdge[i] );
174
175   } // loop on edges
176
177   vExp.Initialize( theEdges.back() );
178   if ( vExp.Value().Orientation() != TopAbs_REVERSED ) vExp.Next();
179   if ( vExp.More() )
180   {
181     if ( SMESH_Algo::VertexNode( TopoDS::Vertex( vExp.Value()), meshDS ))
182       myNbPonits++; // for the last end
183     else
184       myMissingVertexNodes = true;
185   }
186   if ( nbEdges > 1 && myLength > DBL_MIN ) {
187     const double degenNormLen = 1.e-5;
188     double totLength = myLength;
189     if ( nbDegen )
190       totLength += myLength * degenNormLen * nbDegen;
191     double prevNormPar = 0;
192     for ( int i = 0; i < nbEdges; ++i ) {
193       if ( myEdgeLength[ i ] < DBL_MIN )
194         myEdgeLength[ i ] = myLength * degenNormLen;
195       myNormPar[ i ] = prevNormPar + myEdgeLength[i]/totLength;
196       prevNormPar = myNormPar[ i ];
197     }
198   }
199   myNormPar[nbEdges-1] = 1.;
200   //dump();
201 }
202
203 //================================================================================
204 /*!
205  * \brief Constructor of a side for vertex using data from other FaceSide
206   * \param theVertex - the vertex
207   * \param theSide - the side
208  */
209 //================================================================================
210
211 StdMeshers_FaceSide::StdMeshers_FaceSide(const SMDS_MeshNode* theNode,
212                                          const gp_Pnt2d thePnt2d,
213                                          const StdMeshers_FaceSide* theSide)
214 {
215   myC2d.resize(1);
216   myLength       = 0;
217   myProxyMesh    = theSide->myProxyMesh;
218   myDefaultPnt2d = thePnt2d;
219
220   myPoints = theSide->GetUVPtStruct();
221   myNbPonits = myNbSegments = myPoints.size();
222   std::vector<uvPtStruct>::iterator it = myPoints.begin();
223   for(; it!=myPoints.end(); it++) {
224     (*it).u = thePnt2d.X();
225     (*it).v = thePnt2d.Y();
226     (*it).y = 0.0;
227     (*it).node = theNode;
228   }
229 }
230
231 //================================================================================
232 /*!
233  * \brief Return info on nodes on the side
234   * \retval UVPtStruct* - array of data structures
235  */
236 //================================================================================
237
238 const vector<UVPtStruct>& StdMeshers_FaceSide::GetUVPtStruct(bool   isXConst,
239                                                              double constValue) const
240 {
241   if ( myPoints.empty() ) {
242
243     if ( NbEdges() == 0 ) return myPoints;
244
245     SMESHDS_Mesh* meshDS = myProxyMesh->GetMeshDS();
246     SMESH_MesherHelper helper(*myProxyMesh->GetMesh());
247     bool paramOK;
248     double eps = 1e-100;
249
250     // sort nodes of all edges putting them into a map
251
252     map< double, const SMDS_MeshNode*> u2node;
253     vector< const SMESH_ProxyMesh::SubMesh* > proxySubMesh( myEdge.size());
254     int nbProxyNodes = 0;
255     for ( size_t iE = 0; iE < myEdge.size(); ++iE )
256     {
257       proxySubMesh[iE] = myProxyMesh->GetProxySubMesh( myEdge[iE] );
258       if ( proxySubMesh[iE] )
259       {
260         if ( proxySubMesh[iE]->GetUVPtStructVec().empty() ) {
261           proxySubMesh[iE] = 0;
262         }
263         else {
264           nbProxyNodes += proxySubMesh[iE]->GetUVPtStructVec().size() - 1;
265           if ( iE+1 == myEdge.size() )
266             ++nbProxyNodes;
267           continue;
268         }
269       }
270       // Put 1st vertex node of a current edge
271       TopoDS_Vertex VV[2]; // TopExp::FirstVertex() returns NULL for INTERNAL edge
272       VV[0] = SMESH_MesherHelper::IthVertex( 0, myEdge[iE]);
273       VV[1] = SMESH_MesherHelper::IthVertex( 1, myEdge[iE]);
274       const SMDS_MeshNode* node = SMESH_Algo::VertexNode( VV[0], meshDS );
275       double prevNormPar = ( iE == 0 ? 0 : myNormPar[ iE-1 ]); // normalized param
276       if ( node ) { // nodes on internal vertices may be missing
277         u2node.insert( u2node.end(), make_pair( prevNormPar, node ));
278       }
279       else if ( iE == 0 ) {
280         MESSAGE(" NO NODE on VERTEX" );
281         return myPoints;
282       }
283
284       // Put internal nodes
285       if ( const SMESHDS_SubMesh* sm = myProxyMesh->GetSubMesh( myEdge[iE] ))
286       {
287         vector< pair< double, const SMDS_MeshNode*> > u2nodeVec;
288         u2nodeVec.reserve( sm->NbNodes() );
289         SMDS_NodeIteratorPtr nItr = sm->GetNodes();
290         double paramSize = myLast[iE] - myFirst[iE];
291         double r         = myNormPar[iE] - prevNormPar;
292         helper.SetSubShape( myEdge[iE] );
293         helper.ToFixNodeParameters( true );
294         if ( !myIsUniform[iE] )
295           while ( nItr->more() )
296           {
297             const SMDS_MeshNode* node = nItr->next();
298             if ( myIgnoreMediumNodes && SMESH_MeshEditor::IsMedium( node, SMDSAbs_Edge ))
299               continue;
300             double u = helper.GetNodeU( myEdge[iE], node, 0, &paramOK );
301             double aLenU = GCPnts_AbscissaPoint::Length
302               ( const_cast<GeomAdaptor_Curve&>( myC3dAdaptor[iE]), myFirst[iE], u );
303             if ( myEdgeLength[iE] < aLenU ) // nonregression test "3D_mesh_NETGEN/G6"
304             {
305               u2nodeVec.clear();
306               break;
307             }
308             double normPar = prevNormPar + r*aLenU/myEdgeLength[iE];
309             u2nodeVec.push_back( make_pair( normPar, node ));
310           }
311         nItr = sm->GetNodes();
312         if ( u2nodeVec.empty() )
313           while ( nItr->more() )
314           {
315             const SMDS_MeshNode* node = nItr->next();
316             if ( myIgnoreMediumNodes && SMESH_MeshEditor::IsMedium( node, SMDSAbs_Edge ))
317               continue;
318             double u = helper.GetNodeU( myEdge[iE], node, 0, &paramOK );
319
320             // paramSize is signed so orientation is taken into account
321             double normPar = prevNormPar + r * ( u - myFirst[iE] ) / paramSize;
322             u2nodeVec.push_back( make_pair( normPar, node ));
323           }
324         for ( size_t j = 0; j < u2nodeVec.size(); ++j )
325           u2node.insert( u2node.end(), u2nodeVec[j] );
326       }
327
328       // Put 2nd vertex node for a last edge
329       if ( iE+1 == myEdge.size() ) {
330         node = SMESH_Algo::VertexNode( VV[1], meshDS );
331         if ( !node ) {
332           MESSAGE(" NO NODE on VERTEX" );
333           return myPoints;
334         }
335         u2node.insert( u2node.end(), make_pair( 1., node ));
336       }
337     } // loop on myEdge's
338
339     if ( u2node.size() + nbProxyNodes != myNbPonits )
340     {
341       MESSAGE("Wrong node parameters on edges, u2node.size():"
342               <<u2node.size()<<" !=  myNbPonits:"<<myNbPonits);
343       return myPoints;
344     }
345
346     // fill array of UVPtStruct
347
348     UVPtStructVec& points = const_cast< UVPtStructVec& >( myPoints );
349     points.resize( myNbPonits );
350
351     int iPt = 0;
352     double prevNormPar = 0, paramSize = myNormPar[ 0 ];
353     map< double, const SMDS_MeshNode*>::iterator u_node = u2node.begin();
354     for ( size_t iE = 0; iE < myEdge.size(); ++iE )
355     {
356       if ( proxySubMesh[ iE ] ) // copy data from a proxy sub-mesh
357       {
358         const UVPtStructVec& edgeUVPtStruct = proxySubMesh[iE]->GetUVPtStructVec();
359         std::copy( edgeUVPtStruct.begin(), edgeUVPtStruct.end(), & points[iPt] );
360         // update normalized params
361         if ( myEdge.size() > 1 ) {
362           for ( size_t i = 0; i < edgeUVPtStruct.size(); ++i, ++iPt )
363           {
364             UVPtStruct & uvPt = points[iPt];
365             uvPt.normParam = prevNormPar + uvPt.normParam * paramSize;
366             uvPt.x = uvPt.y = uvPt.normParam;
367           }
368           --iPt; // to point to the 1st VERTEX of the next EDGE
369         }
370       }
371       else
372       {
373         for ( ; u_node != u2node.end(); ++u_node, ++iPt )
374         {
375           if ( myNormPar[ iE ]-eps < u_node->first )
376             break; // u_node is at VERTEX of the next EDGE 
377
378           UVPtStruct & uvPt = points[iPt];
379           uvPt.node       = u_node->second;
380           // -- normParam, x, y --------------------------------
381           uvPt.normParam  = u_node->first;
382           uvPt.x = uvPt.y = uvPt.normParam;
383           // -- U ----------------------------------------------
384           const SMDS_EdgePosition* epos =
385             dynamic_cast<const SMDS_EdgePosition*>(uvPt.node->GetPosition());
386           if ( epos ) {
387             uvPt.param = epos->GetUParameter();
388           }
389           else {
390             double r = ( uvPt.normParam - prevNormPar )/ paramSize;
391             uvPt.param = ( r > 0.5 ? myLast[iE] : myFirst[iE] );
392           }
393           // -- UV ---------------------------------------------
394           if ( !myC2d[ iE ].IsNull() ) {
395             gp_Pnt2d p = myC2d[ iE ]->Value( uvPt.param );
396             uvPt.u = p.X();
397             uvPt.v = p.Y();
398           }
399           else {
400             uvPt.u = uvPt.v = 1e+100;
401           }
402         }
403       }
404       // prepare for the next EDGE
405       if ( iE+1 < myEdge.size() )
406       {
407         prevNormPar = myNormPar[ iE ];
408         paramSize   = myNormPar[ iE+1 ] - prevNormPar;
409       }
410     } // loop on myEdge's
411
412     // set <constValue>
413     if ( isXConst )
414       for ( iPt = 0; iPt < points.size(); ++iPt ) points[ iPt ].x = constValue;
415     else
416       for ( iPt = 0; iPt < points.size(); ++iPt ) points[ iPt ].y = constValue;
417
418   } // if ( myPoints.empty())
419
420   return myPoints;
421 }
422
423 //================================================================================
424 /*!
425  * \brief Falsificate info on nodes
426  * \param nbSeg - nb of segments on the side
427  * \retval UVPtStruct* - array of data structures
428  */
429 //================================================================================
430
431 const vector<UVPtStruct>& StdMeshers_FaceSide::SimulateUVPtStruct(int    nbSeg,
432                                                                   bool   isXConst,
433                                                                   double constValue) const
434 {
435   if ( myFalsePoints.empty() ) {
436
437     if ( NbEdges() == 0 ) return myFalsePoints;
438
439     vector<uvPtStruct>* points = const_cast<vector<uvPtStruct>*>( &myFalsePoints );
440     points->resize( nbSeg+1 );
441
442     int EdgeIndex = 0;
443     double prevNormPar = 0, paramSize = myNormPar[ EdgeIndex ];
444     for (int i = 0 ; i < myFalsePoints.size(); ++i ) {
445       double normPar = double(i) / double(nbSeg);
446       UVPtStruct & uvPt = (*points)[i];
447       uvPt.node = 0;
448       uvPt.x = uvPt.y = uvPt.param = uvPt.normParam = normPar;
449       if ( isXConst ) uvPt.x = constValue;
450       else            uvPt.y = constValue;
451       if ( myNormPar[ EdgeIndex ] < normPar ) {
452         prevNormPar = myNormPar[ EdgeIndex ];
453         ++EdgeIndex;
454         paramSize = myNormPar[ EdgeIndex ] - prevNormPar;
455       }
456       double r = ( normPar - prevNormPar )/ paramSize;
457       uvPt.param = myFirst[EdgeIndex] * ( 1 - r ) + myLast[EdgeIndex] * r;
458       if ( !myC2d[ EdgeIndex ].IsNull() ) {
459         gp_Pnt2d p = myC2d[ EdgeIndex ]->Value( uvPt.param );
460         uvPt.u = p.X();
461         uvPt.v = p.Y();
462       }
463       else {
464         uvPt.u = uvPt.v = 1e+100;
465       }
466     }
467   }
468   return myFalsePoints;
469 }
470
471 //=======================================================================
472 //function : GetOrderedNodes
473 //purpose  : Return nodes in the order they encounter while walking along the side
474 //=======================================================================
475
476 std::vector<const SMDS_MeshNode*> StdMeshers_FaceSide::GetOrderedNodes() const
477 {
478   vector<const SMDS_MeshNode*> resultNodes;
479   if ( myPoints.empty() )
480   {
481     if ( NbEdges() == 0 ) return resultNodes;
482
483     SMESHDS_Mesh* meshDS = myProxyMesh->GetMeshDS();
484     SMESH_MesherHelper helper(*myProxyMesh->GetMesh());
485     bool paramOK;
486
487     // Sort nodes of all edges putting them into a map
488
489     map< double, const SMDS_MeshNode*> u2node;
490     for ( int i = 0; i < myEdge.size(); ++i )
491     {
492       // Put 1st vertex node of a current edge
493       TopoDS_Vertex VV[2]; // TopExp::FirstVertex() returns NULL for INTERNAL edge
494       VV[0] = SMESH_MesherHelper::IthVertex( 0, myEdge[i]);
495       VV[1] = SMESH_MesherHelper::IthVertex( 1, myEdge[i]);
496       const SMDS_MeshNode* node = SMESH_Algo::VertexNode( VV[0], meshDS );
497       double prevNormPar = ( i == 0 ? 0 : myNormPar[ i-1 ]); // normalized param
498       if ( node ) { // nodes on internal vertices may be missing
499         u2node.insert( make_pair( prevNormPar, node ));
500       }
501       else if ( i == 0 ) {
502         MESSAGE(" NO NODE on VERTEX" );
503         return resultNodes;
504       }
505
506       // Put internal nodes
507       if ( SMESHDS_SubMesh* sm = meshDS->MeshElements( myEdge[i] ))
508       {
509         SMDS_NodeIteratorPtr nItr = sm->GetNodes();
510         double paramSize = myLast[i] - myFirst[i];
511         double r         = myNormPar[i] - prevNormPar;
512         helper.SetSubShape( myEdge[i] );
513         helper.ToFixNodeParameters( true );
514         while ( nItr->more() )
515         {
516           const SMDS_MeshNode* node = nItr->next();
517           if ( myIgnoreMediumNodes && SMESH_MeshEditor::IsMedium( node, SMDSAbs_Edge ))
518             continue;
519           double u = helper.GetNodeU( myEdge[i], node, 0, &paramOK );
520
521           // paramSize is signed so orientation is taken into account
522           double normPar = prevNormPar + r * ( u - myFirst[i] ) / paramSize;
523           u2node.insert( u2node.end(), make_pair( normPar, node ));
524         }
525       }
526
527       // Put 2nd vertex node for a last edge
528       if ( i+1 == myEdge.size() ) {
529         node = SMESH_Algo::VertexNode( VV[1], meshDS );
530         if ( !node ) {
531           return resultNodes;
532         }
533         u2node.insert( u2node.end(), make_pair( 1., node ));
534       }
535     }
536
537     // Fill the result vector
538
539     if ( u2node.size() == myNbPonits )
540     {
541       resultNodes.reserve( u2node.size() );
542       map< double, const SMDS_MeshNode*>::iterator u2n = u2node.begin();
543       for ( ; u2n != u2node.end(); ++u2n )
544         resultNodes.push_back( u2n->second );
545     }
546   }
547   else
548   {
549     resultNodes.resize( myPoints.size() );
550     for ( size_t i = 0; i < myPoints.size(); ++i )
551       resultNodes[i] = myPoints[i].node;
552   }
553
554   return resultNodes;
555 }
556
557 //================================================================================
558 /*!
559  * \brief reverse order of vector elements
560   * \param vec - vector to reverse
561  */
562 //================================================================================
563
564 template <typename T > void reverse(vector<T> & vec)
565 {
566   for ( int f=0, r=vec.size()-1; f < r; ++f, --r )
567     std::swap( vec[f], vec[r] );
568 }
569
570 //================================================================================
571 /*!
572  * \brief Change orientation of side geometry
573  */
574 //================================================================================
575
576 void StdMeshers_FaceSide::Reverse()
577 {
578   int nbEdges = myEdge.size();
579   for ( int i = nbEdges-1; i >= 0; --i ) {
580     std::swap( myFirst[i], myLast[i] );
581     myEdge[i].Reverse();
582     if ( i > 0 ) // at the first loop 1. is overwritten
583       myNormPar[i] = 1 - myNormPar[i-1];
584   }
585   if ( nbEdges > 1 ) {
586     reverse( myEdge );
587     reverse( myEdgeID );
588     reverse( myC2d );
589     reverse( myC3dAdaptor );
590     reverse( myFirst );
591     reverse( myLast );
592     reverse( myNormPar );
593     reverse( myEdgeLength );
594     reverse( myIsUniform );
595   }
596   if ( nbEdges > 0 )
597   {
598     myNormPar[nbEdges-1]=1.;
599     myPoints.clear();
600     myFalsePoints.clear();
601     for ( size_t i = 0; i < myEdge.size(); ++i )
602       reverseProxySubmesh( myEdge[i] );
603   }
604 }
605
606 //================================================================================
607 /*!
608  * \brief Reverse UVPtStructVec if a proxy sub-mesh of E
609  */
610 //================================================================================
611
612 void StdMeshers_FaceSide::reverseProxySubmesh( const TopoDS_Edge& E )
613 {
614   if ( !myProxyMesh ) return;
615   if ( const SMESH_ProxyMesh::SubMesh* sm = myProxyMesh->GetProxySubMesh( E ))
616   {
617     UVPtStructVec& edgeUVPtStruct = (UVPtStructVec& ) sm->GetUVPtStructVec();
618     for ( size_t i = 0; i < edgeUVPtStruct.size(); ++i )
619     {
620       UVPtStruct & uvPt = edgeUVPtStruct[i];
621       uvPt.normParam = 1 - uvPt.normParam;
622       uvPt.x         = 1 - uvPt.x;
623       uvPt.y         = 1 - uvPt.y;
624     }
625     reverse( edgeUVPtStruct );
626   }
627 }
628
629 //================================================================================
630 /*!
631  * \brief Show side features
632  */
633 //================================================================================
634
635 void StdMeshers_FaceSide::dump(const char* msg) const
636 {
637 #ifdef _DEBUG_
638   if (msg) MESSAGE ( std::endl << msg );
639   MESSAGE_BEGIN ("NB EDGES: "<< myEdge.size() );
640   MESSAGE_ADD ( "nbPoints: "<<myNbPonits<<" vecSize: " << myPoints.size()<<" "<<myFalsePoints.size() );
641   for ( int i=0; i<myEdge.size(); ++i)
642   {
643     MESSAGE_ADD ( "\t"<<i+1 );
644     MESSAGE_ADD ( "\tEDGE: " );
645     if (myEdge[i].IsNull()) {
646       MESSAGE_ADD ( "NULL" );
647     }
648     else {
649       TopAbs::Print(myEdge[i].Orientation(),cout)<<" "<<myEdge[i].TShape().operator->()<<endl;
650       MESSAGE_ADD ( "\tV1: " << TopExp::FirstVertex( myEdge[i], 1).TShape().operator->()
651                  << "  V2: " << TopExp::LastVertex( myEdge[i], 1).TShape().operator->() );
652     }
653     MESSAGE_ADD ( "\tC2d: ");
654     
655     if (myC2d[i].IsNull()) {
656       MESSAGE_ADD ( "NULL" );
657     }
658     else {
659       MESSAGE_ADD ( myC2d[i].operator->() );
660     }
661       
662     MESSAGE_ADD ( "\tF: "<<myFirst[i]<< " L: "<< myLast[i] );
663     MESSAGE_END ( "\tnormPar: "<<myNormPar[i]<<endl );
664   }
665 #endif
666 }
667
668 //================================================================================
669 /*!
670  * \brief Creates a Adaptor2d_Curve2d to be used in SMESH_Block
671   * \retval Adaptor2d_Curve2d* - 
672  */
673 //================================================================================
674
675 struct Adaptor2dCurve2d : public Adaptor2d_Curve2d
676 {
677   const StdMeshers_FaceSide* mySide;
678   Adaptor2dCurve2d(const StdMeshers_FaceSide* faceSide):mySide(faceSide) {}
679   gp_Pnt2d Value(const Standard_Real U) const { return mySide->Value2d( U ); }
680   Standard_Real FirstParameter() const { return 0; }
681   Standard_Real LastParameter() const { return 1; }
682 };
683
684 Adaptor2d_Curve2d* StdMeshers_FaceSide::GetCurve2d() const
685 {
686   return new Adaptor2dCurve2d( this );
687 }
688
689 //================================================================================
690 /*!
691  * \brief Creates a fully functional Adaptor_Curve
692  */
693 //================================================================================
694
695 BRepAdaptor_CompCurve* StdMeshers_FaceSide::GetCurve3d() const
696 {
697   if ( myEdge.empty() )
698     return 0;
699
700   TopoDS_Wire aWire;
701   BRep_Builder aBuilder;
702   aBuilder.MakeWire(aWire);
703   for ( int i=0; i<myEdge.size(); ++i )
704     aBuilder.Add( aWire, myEdge[i] );
705
706   if ( myEdge.size() == 2 && FirstVertex().IsSame( LastVertex() ))
707     aWire.Closed(true); // issue 0021141
708
709   return new BRepAdaptor_CompCurve( aWire );
710 }
711
712
713 //================================================================================
714 /*!
715  * \brief Return 2D point by normalized parameter
716   * \param U - normalized parameter value
717   * \retval gp_Pnt2d - point
718  */
719 //================================================================================
720
721 gp_Pnt2d StdMeshers_FaceSide::Value2d(double U) const
722 {
723   if ( !myC2d[0].IsNull() ) {
724     int i = EdgeIndex( U );
725     double prevU = i ? myNormPar[ i-1 ] : 0;
726     double r = ( U - prevU )/ ( myNormPar[ i ] - prevU );
727
728     double par = myFirst[i] * ( 1 - r ) + myLast[i] * r;
729     
730     // check parametrization of curve
731     if( !myIsUniform[i] )
732     {
733       double aLen3dU = r * myEdgeLength[i] * ( myFirst[i]>myLast[i] ? -1. : 1.);
734       GCPnts_AbscissaPoint AbPnt
735         ( const_cast<GeomAdaptor_Curve&>( myC3dAdaptor[i]), aLen3dU, myFirst[i] );
736       if( AbPnt.IsDone() ) {
737         par = AbPnt.Parameter();
738       }
739     }
740     return myC2d[ i ]->Value(par);
741
742   }
743   return myDefaultPnt2d;
744 }
745
746 //================================================================================
747 /*!
748  * \brief Return wires of a face as StdMeshers_FaceSide's
749  */
750 //================================================================================
751
752 TSideVector StdMeshers_FaceSide::GetFaceWires(const TopoDS_Face&   theFace,
753                                               SMESH_Mesh &         theMesh,
754                                               const bool           theIgnoreMediumNodes,
755                                               TError &             theError,
756                                               SMESH_ProxyMesh::Ptr theProxyMesh)
757 {
758   TopoDS_Vertex V1;
759   list< TopoDS_Edge > edges, internalEdges;
760   list< int > nbEdgesInWires;
761   int nbWires = SMESH_Block::GetOrderedEdges (theFace, V1, edges, nbEdgesInWires);
762
763   // split list of all edges into separate wires
764   TSideVector wires( nbWires );
765   list< int >::iterator nbE = nbEdgesInWires.begin();
766   list< TopoDS_Edge >::iterator from = edges.begin(), to = from;
767   for ( int iW = 0; iW < nbWires; ++iW, ++nbE )
768   {
769     std::advance( to, *nbE );
770     if ( *nbE == 0 ) // Issue 0020676
771     {
772       --nbWires;
773       --iW;
774       wires.resize( nbWires );
775       continue;
776     }
777     list< TopoDS_Edge > wireEdges( from, to );
778     // assure that there is a node on the first vertex
779     // as StdMeshers_FaceSide::GetUVPtStruct() requires
780     if ( wireEdges.front().Orientation() != TopAbs_INTERNAL ) // Issue 0020676
781     {
782       while ( !SMESH_Algo::VertexNode( TopExp::FirstVertex( wireEdges.front(), true),
783                                        theMesh.GetMeshDS()))
784       {
785         wireEdges.splice(wireEdges.end(), wireEdges,
786                          wireEdges.begin(), ++wireEdges.begin());
787         if ( from->IsSame( wireEdges.front() )) {
788           theError = TError
789             ( new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,"No nodes on vertices"));
790           return TSideVector(0);
791         }
792       }
793     }
794     else if ( *nbE > 1 ) // Issue 0020676 (Face_pb_netgen.brep) - several internal edges in a wire
795     {
796       internalEdges.splice( internalEdges.end(), wireEdges, ++wireEdges.begin(), wireEdges.end());
797     }
798
799     StdMeshers_FaceSide* wire = new StdMeshers_FaceSide( theFace, wireEdges, &theMesh,
800                                                          /*isForward=*/true, theIgnoreMediumNodes,
801                                                          theProxyMesh );
802     wires[ iW ] = StdMeshers_FaceSidePtr( wire );
803     from = to;
804   }
805   while ( !internalEdges.empty() )
806   {
807     StdMeshers_FaceSide* wire = new StdMeshers_FaceSide( theFace, internalEdges.back(), &theMesh,
808                                                          /*isForward=*/true, theIgnoreMediumNodes,
809                                                          theProxyMesh );
810     wires.push_back( StdMeshers_FaceSidePtr( wire ));
811     internalEdges.pop_back();
812   }
813   return wires;
814 }
815
816 //================================================================================
817 /*!
818  * \brief Return 1st vertex of the i-the edge
819  */
820 //================================================================================
821
822 TopoDS_Vertex StdMeshers_FaceSide::FirstVertex(int i) const
823 {
824   TopoDS_Vertex v;
825   if ( i < NbEdges() )
826   {
827     v = myEdge[i].Orientation() <= TopAbs_REVERSED ? // FORWARD || REVERSED
828         TopExp::FirstVertex( myEdge[i], 1 )        :
829         TopoDS::Vertex( TopoDS_Iterator( myEdge[i] ).Value() );
830   }
831   return v;
832 }
833
834 //================================================================================
835 /*!
836  * \brief Return last vertex of the i-the edge
837  */
838 //================================================================================
839
840 TopoDS_Vertex StdMeshers_FaceSide::LastVertex(int i) const
841 {
842   TopoDS_Vertex v;
843   if ( i < NbEdges() )
844   {
845     const TopoDS_Edge& edge = i<0 ? myEdge[ NbEdges() + i ] : myEdge[i];
846     if ( edge.Orientation() <= TopAbs_REVERSED ) // FORWARD || REVERSED
847       v = TopExp::LastVertex( edge, 1 );
848     else
849       for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() )
850         v = TopoDS::Vertex( vIt.Value() );
851   }
852   return v;
853 }
854