]> SALOME platform Git repositories - modules/smesh.git/blob - src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx
Salome HOME
0021825: Error in the example of "Projection Algorithms" in the user's guide
[modules/smesh.git] / src / StdMeshers / StdMeshers_RadialQuadrangle_1D2D.cxx
1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
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.
7 //
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.
12 //
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
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19
20 //  SMESH SMESH : implementaion of SMESH idl descriptions
21 // File      : StdMeshers_RadialQuadrangle_1D2D.cxx
22 // Module    : SMESH
23
24 #include "StdMeshers_RadialQuadrangle_1D2D.hxx"
25
26 #include "StdMeshers_NumberOfLayers.hxx"
27 #include "StdMeshers_LayerDistribution.hxx"
28 #include "StdMeshers_Regular_1D.hxx"
29 #include "StdMeshers_NumberOfSegments.hxx"
30
31 #include "SMDS_MeshNode.hxx"
32 #include "SMESHDS_SubMesh.hxx"
33 #include "SMESH_Gen.hxx"
34 #include "SMESH_HypoFilter.hxx"
35 #include "SMESH_Mesh.hxx"
36 #include "SMESH_MesherHelper.hxx"
37 #include "SMESH_subMesh.hxx"
38 #include "SMESH_subMeshEventListener.hxx"
39
40 #include "utilities.h"
41
42 #include <BRepAdaptor_Curve.hxx>
43 #include <BRepBuilderAPI_MakeEdge.hxx>
44 #include <BRep_Tool.hxx>
45 #include <GeomAPI_ProjectPointOnSurf.hxx>
46 #include <Geom_Circle.hxx>
47 #include <Geom_Line.hxx>
48 #include <Geom_TrimmedCurve.hxx>
49 #include <TColgp_SequenceOfPnt.hxx>
50 #include <TColgp_SequenceOfPnt2d.hxx>
51 #include <TopExp.hxx>
52 #include <TopExp_Explorer.hxx>
53 #include <TopTools_ListIteratorOfListOfShape.hxx>
54 #include <TopoDS.hxx>
55
56
57 using namespace std;
58
59 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
60 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
61
62
63 //=======================================================================
64 //function : StdMeshers_RadialQuadrangle_1D2D
65 //purpose  : 
66 //=======================================================================
67
68 StdMeshers_RadialQuadrangle_1D2D::StdMeshers_RadialQuadrangle_1D2D(int hypId,
69                                                                    int studyId,
70                                                                    SMESH_Gen* gen)
71   :SMESH_2D_Algo(hypId, studyId, gen)
72 {
73   _name = "RadialQuadrangle_1D2D";
74   _shapeType = (1 << TopAbs_FACE);        // 1 bit per shape type
75
76   _compatibleHypothesis.push_back("LayerDistribution2D");
77   _compatibleHypothesis.push_back("NumberOfLayers2D");
78   myNbLayerHypo = 0;
79   myDistributionHypo = 0;
80   _requireDiscreteBoundary = false;
81   _supportSubmeshes = true;
82 }
83
84
85 //================================================================================
86 /*!
87  * \brief Destructor
88  */
89 //================================================================================
90
91 StdMeshers_RadialQuadrangle_1D2D::~StdMeshers_RadialQuadrangle_1D2D()
92 {}
93
94
95 //=======================================================================
96 //function : CheckHypothesis
97 //purpose  : 
98 //=======================================================================
99
100 bool StdMeshers_RadialQuadrangle_1D2D::CheckHypothesis
101                            (SMESH_Mesh&                          aMesh,
102                             const TopoDS_Shape&                  aShape,
103                             SMESH_Hypothesis::Hypothesis_Status& aStatus)
104 {
105   // check aShape 
106   myNbLayerHypo = 0;
107   myDistributionHypo = 0;
108
109   list <const SMESHDS_Hypothesis * >::const_iterator itl;
110
111   const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
112   if ( hyps.size() == 0 ) {
113     aStatus = SMESH_Hypothesis::HYP_OK;
114     return true;  // can work with no hypothesis
115   }
116
117   if ( hyps.size() > 1 ) {
118     aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
119     return false;
120   }
121
122   const SMESHDS_Hypothesis *theHyp = hyps.front();
123
124   string hypName = theHyp->GetName();
125
126   if (hypName == "NumberOfLayers2D") {
127     myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
128     aStatus = SMESH_Hypothesis::HYP_OK;
129     return true;
130   }
131   if (hypName == "LayerDistribution2D") {
132     myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
133     aStatus = SMESH_Hypothesis::HYP_OK;
134     return true;
135   }
136   aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
137   return true;
138 }
139
140 namespace
141 {
142   // ------------------------------------------------------------------------------
143   /*!
144    * \brief Listener used to mark edges meshed by StdMeshers_RadialQuadrangle_1D2D
145    */
146   class TEdgeMarker : public SMESH_subMeshEventListener
147   {
148     TEdgeMarker(): SMESH_subMeshEventListener(/*isDeletable=*/false,
149                                               "StdMeshers_RadialQuadrangle_1D2D::TEdgeMarker") {}
150   public:
151     //!<  Return static listener
152     static SMESH_subMeshEventListener* getListener()
153     {
154       static TEdgeMarker theEdgeMarker;
155       return &theEdgeMarker;
156     }
157     //! Clear face sumbesh if something happens on edges
158     void ProcessEvent(const int          event,
159                       const int          eventType,
160                       SMESH_subMesh*     edgeSubMesh,
161                       EventListenerData* data,
162                       const SMESH_Hypothesis*  /*hyp*/)
163     {
164       if ( data && !data->mySubMeshes.empty() && eventType == SMESH_subMesh::ALGO_EVENT)
165       {
166         ASSERT( data->mySubMeshes.front() != edgeSubMesh );
167         SMESH_subMesh* faceSubMesh = data->mySubMeshes.front();
168         faceSubMesh->ComputeStateEngine( SMESH_subMesh::CLEAN );
169       }
170     }
171   };
172
173   // ------------------------------------------------------------------------------
174   /*!
175    * \brief Mark an edge as computed by StdMeshers_RadialQuadrangle_1D2D
176    */
177   void markEdgeAsComputedByMe(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
178   {
179     if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
180     {
181       if ( !edgeSM->GetEventListenerData( TEdgeMarker::getListener() ))
182         faceSubMesh->SetEventListener( TEdgeMarker::getListener(),
183                                        SMESH_subMeshEventListenerData::MakeData(faceSubMesh),
184                                        edgeSM);
185     }
186   }
187   // ------------------------------------------------------------------------------
188   /*!
189    * \brief Return true if a radial edge was meshed with StdMeshers_RadialQuadrangle_1D2D with
190    * the same radial distribution
191    */
192 //   bool isEdgeCompatiballyMeshed(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
193 //   {
194 //     if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
195 //     {
196 //       if ( SMESH_subMeshEventListenerData* otherFaceData =
197 //            edgeSM->GetEventListenerData( TEdgeMarker::getListener() ))
198 //       {
199 //         // compare hypothesis aplied to two disk faces sharing radial edges
200 //         SMESH_Mesh& mesh = *faceSubMesh->GetFather();
201 //         SMESH_Algo* radialQuadAlgo = mesh.GetGen()->GetAlgo(mesh, faceSubMesh->GetSubShape() );
202 //         SMESH_subMesh* otherFaceSubMesh = otherFaceData->mySubMeshes.front();
203 //         list <const SMESHDS_Hypothesis *> hyps1 =
204 //           radialQuadAlgo->GetUsedHypothesis( mesh, faceSubMesh->GetSubShape());
205 //         list <const SMESHDS_Hypothesis *> hyps2 =
206 //           radialQuadAlgo->GetUsedHypothesis( mesh, otherFaceSubMesh->GetSubShape());
207 //         if( hyps1.empty() && hyps2.empty() )
208 //           return true; // defaul hyps
209 //         if ( hyps1.size() != hyps2.size() )
210 //           return false;
211 //         return *hyps1.front() == *hyps2.front();
212 //       }
213 //     }
214 //     return false;
215 //   }
216
217   //================================================================================
218   /*!
219    * \brief Return base curve of the edge and extremum parameters
220    */
221   //================================================================================
222
223   Handle(Geom_Curve) getCurve(const TopoDS_Edge& edge, double* f=0, double* l=0)
224   {
225     Handle(Geom_Curve) C;
226     if ( !edge.IsNull() )
227     {
228       double first = 0., last = 0.;
229       C = BRep_Tool::Curve(edge, first, last);
230       if ( !C.IsNull() )
231       {
232         Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C);
233         while( !tc.IsNull() ) {
234           C = tc->BasisCurve();
235           tc = Handle(Geom_TrimmedCurve)::DownCast(C);
236         }
237         if ( f ) *f = first;
238         if ( l ) *l = last;
239       }
240     }
241     return C;
242   }
243
244   //================================================================================
245   /*!
246    * \brief Return edges of the face
247    *  \retval int - nb of edges
248    */
249   //================================================================================
250
251   int analyseFace(const TopoDS_Shape& face,
252                   TopoDS_Edge&        CircEdge,
253                   TopoDS_Edge&        LinEdge1,
254                   TopoDS_Edge&        LinEdge2)
255   {
256     CircEdge.Nullify(); LinEdge1.Nullify(); LinEdge2.Nullify();
257     int nbe = 0;
258
259     for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next(), ++nbe )
260     {
261       const TopoDS_Edge& E = TopoDS::Edge( exp.Current() );
262       double f,l;
263       Handle(Geom_Curve) C = getCurve(E,&f,&l);
264       if ( !C.IsNull() )
265       {
266         if ( C->IsKind( STANDARD_TYPE(Geom_Circle)))
267         {
268           if ( CircEdge.IsNull() )
269             CircEdge = E;
270           else
271             return 0;
272         }
273         else if ( LinEdge1.IsNull() )
274           LinEdge1 = E;
275         else
276           LinEdge2 = E;
277       }
278     }
279     return nbe;
280   }
281
282 //================================================================================
283 //================================================================================
284 /*!
285  * \brief Class computing layers distribution using data of
286  *        StdMeshers_LayerDistribution hypothesis
287  */
288 //================================================================================
289 //================================================================================
290
291 class TNodeDistributor: public StdMeshers_Regular_1D
292 {
293   list <const SMESHDS_Hypothesis *> myUsedHyps;
294 public:
295   // -----------------------------------------------------------------------------
296   static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
297   {
298     const int myID = -1000;
299     map < int, SMESH_1D_Algo * > & algoMap = aMesh.GetGen()->_map1D_Algo;
300     map < int, SMESH_1D_Algo * >::iterator id_algo = algoMap.find( myID );
301     if ( id_algo == algoMap.end() )
302       return new TNodeDistributor( myID, 0, aMesh.GetGen() );
303     return static_cast< TNodeDistributor* >( id_algo->second );
304   }
305   // -----------------------------------------------------------------------------
306   //! Computes distribution of nodes on a straight line ending at pIn and pOut
307   bool Compute( vector< double > &      positions,
308                 gp_Pnt                  pIn,
309                 gp_Pnt                  pOut,
310                 SMESH_Mesh&             aMesh,
311                 const SMESH_Hypothesis* hyp1d)
312   {
313     if ( !hyp1d ) return error( "Invalid LayerDistribution hypothesis");
314
315     double len = pIn.Distance( pOut );
316     if ( len <= DBL_MIN ) return error("Too close points of inner and outer shells");
317
318     myUsedHyps.clear();
319     myUsedHyps.push_back( hyp1d );
320
321     TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut );
322     SMESH_Hypothesis::Hypothesis_Status aStatus;
323     if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
324       return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
325                     "with LayerDistribution hypothesis");
326
327     BRepAdaptor_Curve C3D(edge);
328     double f = C3D.FirstParameter(), l = C3D.LastParameter();
329     list< double > params;
330     if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false ))
331       return error("StdMeshers_Regular_1D failed to compute layers distribution");
332
333     positions.clear();
334     positions.reserve( params.size() );
335     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
336       positions.push_back( *itU / len );
337     return true;
338   }
339   // -----------------------------------------------------------------------------
340   //! Make mesh on an adge using assigned 1d hyp or defaut nb of segments
341   bool ComputeCircularEdge(SMESH_Mesh&         aMesh,
342                            const TopoDS_Edge& anEdge)
343   {
344     _gen->Compute( aMesh, anEdge);
345     SMESH_subMesh *sm = aMesh.GetSubMesh(anEdge);
346     if ( sm->GetComputeState() != SMESH_subMesh::COMPUTE_OK)
347     {
348       // find any 1d hyp assigned (there can be a hyp w/o algo)
349       myUsedHyps = SMESH_Algo::GetUsedHypothesis(aMesh, anEdge, /*ignoreAux=*/true);
350       Hypothesis_Status aStatus;
351       if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, anEdge, aStatus ))
352       {
353         // no valid 1d hyp assigned, use default nb of segments
354         _hypType                    = NB_SEGMENTS;
355         _ivalue[ DISTR_TYPE_IND ]   = StdMeshers_NumberOfSegments::DT_Regular;
356         _ivalue[ NB_SEGMENTS_IND  ] = _gen->GetDefaultNbSegments();
357       }
358       return StdMeshers_Regular_1D::Compute( aMesh, anEdge );
359     }
360     return true;
361   }
362   // -----------------------------------------------------------------------------
363   //! Make mesh on an adge using assigned 1d hyp or defaut nb of segments
364   bool EvaluateCircularEdge(SMESH_Mesh&        aMesh,
365                             const TopoDS_Edge& anEdge,
366                             MapShapeNbElems&   aResMap)
367   {
368     _gen->Evaluate( aMesh, anEdge, aResMap );
369     if ( aResMap.count( aMesh.GetSubMesh( anEdge )))
370       return true;
371
372     // find any 1d hyp assigned
373     myUsedHyps = SMESH_Algo::GetUsedHypothesis(aMesh, anEdge, /*ignoreAux=*/true);
374     Hypothesis_Status aStatus;
375     if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, anEdge, aStatus ))
376     {
377       // no valid 1d hyp assigned, use default nb of segments
378       _hypType                    = NB_SEGMENTS;
379       _ivalue[ DISTR_TYPE_IND ]   = StdMeshers_NumberOfSegments::DT_Regular;
380       _ivalue[ NB_SEGMENTS_IND  ] = _gen->GetDefaultNbSegments();
381     }
382     return StdMeshers_Regular_1D::Evaluate( aMesh, anEdge, aResMap );
383   }
384 protected:
385   // -----------------------------------------------------------------------------
386   TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen)
387     : StdMeshers_Regular_1D( hypId, studyId, gen)
388   {
389   }
390   // -----------------------------------------------------------------------------
391   virtual const list <const SMESHDS_Hypothesis *> &
392     GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
393   {
394     return myUsedHyps;
395   }
396   // -----------------------------------------------------------------------------
397 };
398 }
399
400 //=======================================================================
401 /*!
402  * \brief Allow algo to do something after persistent restoration
403  * \param subMesh - restored submesh
404  *
405  * call markEdgeAsComputedByMe()
406  */
407 //=======================================================================
408
409 void StdMeshers_RadialQuadrangle_1D2D::SubmeshRestored(SMESH_subMesh* faceSubMesh)
410 {
411   if ( !faceSubMesh->IsEmpty() )
412   {
413     TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
414     analyseFace( faceSubMesh->GetSubShape(), CircEdge, LinEdge1, LinEdge2 );
415     if ( !CircEdge.IsNull() ) markEdgeAsComputedByMe( CircEdge, faceSubMesh );
416     if ( !LinEdge1.IsNull() ) markEdgeAsComputedByMe( LinEdge1, faceSubMesh );
417     if ( !LinEdge2.IsNull() ) markEdgeAsComputedByMe( LinEdge2, faceSubMesh );
418   }
419 }
420
421 //=======================================================================
422 //function : Compute
423 //purpose  : 
424 //=======================================================================
425
426 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
427                                                const TopoDS_Shape& aShape)
428 {
429   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
430
431   myHelper = new SMESH_MesherHelper( aMesh );
432   myHelper->IsQuadraticSubMesh( aShape );
433   // to delete helper at exit from Compute()
434   auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
435
436   TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
437
438   TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
439   int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 );
440   Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
441   if( nbe>3 || nbe < 1 || aCirc.IsNull() )
442     return error("The face must be a full circle or a part of circle (i.e. the number of edges is less or equal to 3 and one of them is a circle curve)");
443   
444   gp_Pnt P0, P1;
445   // points for rotation
446   TColgp_SequenceOfPnt Points;
447   // angles for rotation
448   TColStd_SequenceOfReal Angles;
449   // Nodes1 and Nodes2 - nodes along radiuses
450   // CNodes - nodes on circle edge
451   vector< const SMDS_MeshNode* > Nodes1, Nodes2, CNodes;
452   SMDS_MeshNode * NC;
453   // parameters edge nodes on face
454   TColgp_SequenceOfPnt2d Pnts2d1;
455   gp_Pnt2d PC;
456
457   int faceID = meshDS->ShapeToIndex(aShape);
458   TopoDS_Face F = TopoDS::Face(aShape);
459   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
460
461
462   if(nbe==1)
463   {
464     if (!algo1d->ComputeCircularEdge( aMesh, CircEdge ))
465       return error( algo1d->GetComputeError() );
466     map< double, const SMDS_MeshNode* > theNodes;
467     if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes))
468       return error("Circular edge is incorrectly meshed");
469
470     CNodes.clear();
471     map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
472     const SMDS_MeshNode* NF = (*itn).second;
473     CNodes.push_back( (*itn).second );
474     double fang = (*itn).first;
475     if ( itn != theNodes.end() ) {
476       itn++;
477       for(; itn != theNodes.end(); itn++ ) {
478         CNodes.push_back( (*itn).second );
479         double ang = (*itn).first - fang;
480         if( ang>M_PI ) ang = ang - 2.*M_PI;
481         if( ang<-M_PI ) ang = ang + 2.*M_PI;
482         Angles.Append( ang ); 
483       }
484     }
485     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
486     P0 = aCirc->Location();
487
488     if ( !computeLayerPositions(P0,P1))
489       return false;
490
491     TopoDS_Vertex V1 = myHelper->IthVertex(0, CircEdge );
492     gp_Pnt2d p2dV = BRep_Tool::Parameters( V1, TopoDS::Face(aShape) );
493
494     NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
495     GeomAPI_ProjectPointOnSurf PPS(P0,S);
496     double U0,V0;
497     PPS.Parameters(1,U0,V0);
498     meshDS->SetNodeOnFace(NC, faceID, U0, V0);
499     PC = gp_Pnt2d(U0,V0);
500
501     gp_Vec aVec(P0,P1);
502     gp_Vec2d aVec2d(PC,p2dV);
503     Nodes1.resize( myLayerPositions.size()+1 );
504     Nodes2.resize( myLayerPositions.size()+1 );
505     int i = 0;
506     for(; i<myLayerPositions.size(); i++) {
507       gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
508                 P0.Y() + aVec.Y()*myLayerPositions[i],
509                 P0.Z() + aVec.Z()*myLayerPositions[i] );
510       Points.Append(P);
511       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
512       Nodes1[i] = node;
513       Nodes2[i] = node;
514       double U = PC.X() + aVec2d.X()*myLayerPositions[i];
515       double V = PC.Y() + aVec2d.Y()*myLayerPositions[i];
516       meshDS->SetNodeOnFace( node, faceID, U, V );
517       Pnts2d1.Append(gp_Pnt2d(U,V));
518     }
519     Nodes1[Nodes1.size()-1] = NF;
520     Nodes2[Nodes1.size()-1] = NF;
521   }
522   else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL )
523   {
524     // one curve must be a half of circle and other curve must be
525     // a segment of line
526     double fp, lp;
527     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
528     if( fabs(fabs(lp-fp)-M_PI) > Precision::Confusion() ) {
529       // not half of circle
530       return error(COMPERR_BAD_SHAPE);
531     }
532     Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
533     if( aLine.IsNull() ) {
534       // other curve not line
535       return error(COMPERR_BAD_SHAPE);
536     }
537
538     if ( !algo1d->ComputeCircularEdge( aMesh, CircEdge ))
539       return error( algo1d->GetComputeError() );
540     map< double, const SMDS_MeshNode* > theNodes;
541     if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes) )
542       return error("Circular edge is incorrectly meshed");
543
544     map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
545     CNodes.clear();
546     CNodes.push_back( itn->second );
547     double fang = (*itn).first;
548     itn++;
549     for(; itn != theNodes.end(); itn++ ) {
550       CNodes.push_back( (*itn).second );
551       double ang = (*itn).first - fang;
552       if( ang>M_PI ) ang = ang - 2.*M_PI;
553       if( ang<-M_PI ) ang = ang + 2.*M_PI;
554       Angles.Append( ang );
555     }
556     const SMDS_MeshNode* NF = theNodes.begin()->second;
557     const SMDS_MeshNode* NL = theNodes.rbegin()->second;
558     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
559     gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
560     P0 = aCirc->Location();
561
562     bool linEdgeComputed;
563     if ( !computeLayerPositions(P0,P1,LinEdge1,&linEdgeComputed))
564       return false;
565
566     if ( linEdgeComputed )
567     {
568       if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
569         return error("Invalid mesh on a straight edge");
570
571       Nodes1.resize( myLayerPositions.size()+1 );
572       Nodes2.resize( myLayerPositions.size()+1 );
573       vector< const SMDS_MeshNode* > *pNodes1 = &Nodes1, *pNodes2 = &Nodes2;
574       bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
575       if ( !nodesFromP0ToP1 ) std::swap( pNodes1, pNodes2 );
576
577       map< double, const SMDS_MeshNode* >::reverse_iterator ritn = theNodes.rbegin();
578       itn = theNodes.begin();
579       for ( int i = Nodes1.size()-1; i > -1; ++itn, ++ritn, --i )
580       {
581         (*pNodes1)[i] = ritn->second;
582         (*pNodes2)[i] =  itn->second;
583         Points.Prepend( gpXYZ( Nodes1[i]));
584         Pnts2d1.Prepend( myHelper->GetNodeUV( F, Nodes1[i]));
585       }
586       NC = const_cast<SMDS_MeshNode*>( itn->second );
587       Points.Remove( Nodes1.size() );
588     }
589     else
590     {
591       gp_Vec aVec(P0,P1);
592       int edgeID = meshDS->ShapeToIndex(LinEdge1);
593       // check orientation
594       Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
595       gp_Pnt Ptmp;
596       Crv->D0(fp,Ptmp);
597       bool ori = true;
598       if( P1.Distance(Ptmp) > Precision::Confusion() )
599         ori = false;
600       // get UV points for edge
601       gp_Pnt2d PF,PL;
602       BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
603       PC = gp_Pnt2d( (PF.X()+PL.X())/2, (PF.Y()+PL.Y())/2 );
604       gp_Vec2d V2d;
605       if(ori) V2d = gp_Vec2d(PC,PF);
606       else V2d = gp_Vec2d(PC,PL);
607       // add nodes on edge
608       double cp = (fp+lp)/2;
609       double dp2 = (lp-fp)/2;
610       NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
611       meshDS->SetNodeOnEdge(NC, edgeID, cp);
612       Nodes1.resize( myLayerPositions.size()+1 );
613       Nodes2.resize( myLayerPositions.size()+1 );
614       int i = 0;
615       for(; i<myLayerPositions.size(); i++) {
616         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
617                   P0.Y() + aVec.Y()*myLayerPositions[i],
618                   P0.Z() + aVec.Z()*myLayerPositions[i] );
619         Points.Append(P);
620         SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
621         Nodes1[i] = node;
622         double param;
623         if(ori)
624           param = fp + dp2*(1-myLayerPositions[i]);
625         else
626           param = cp + dp2*myLayerPositions[i];
627         meshDS->SetNodeOnEdge(node, edgeID, param);
628         P = gp_Pnt( P0.X() - aVec.X()*myLayerPositions[i],
629                     P0.Y() - aVec.Y()*myLayerPositions[i],
630                     P0.Z() - aVec.Z()*myLayerPositions[i] );
631         node = meshDS->AddNode(P.X(), P.Y(), P.Z());
632         Nodes2[i] = node;
633         if(!ori)
634           param = fp + dp2*(1-myLayerPositions[i]);
635         else
636           param = cp + dp2*myLayerPositions[i];
637         meshDS->SetNodeOnEdge(node, edgeID, param);
638         // parameters on face
639         gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
640                       PC.Y() + V2d.Y()*myLayerPositions[i] );
641         Pnts2d1.Append(P2d);
642       }
643       Nodes1[ myLayerPositions.size() ] = NF;
644       Nodes2[ myLayerPositions.size() ] = NL;
645       // create 1D elements on edge
646       vector< const SMDS_MeshNode* > tmpNodes;
647       tmpNodes.resize(2*Nodes1.size()+1);
648       for(i=0; i<Nodes2.size(); i++)
649         tmpNodes[Nodes2.size()-i-1] = Nodes2[i];
650       tmpNodes[Nodes2.size()] = NC;
651       for(i=0; i<Nodes1.size(); i++)
652         tmpNodes[Nodes2.size()+1+i] = Nodes1[i];
653       for(i=1; i<tmpNodes.size(); i++) {
654         SMDS_MeshEdge* ME = myHelper->AddEdge( tmpNodes[i-1], tmpNodes[i] );
655         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
656       }
657       markEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
658     }
659   }
660   else // nbe==3 or ( nbe==2 && linEdge is INTERNAL )
661   {
662     if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
663       LinEdge2 = LinEdge1;
664
665     // one curve must be a part of circle and other curves must be
666     // segments of line
667     double fp, lp;
668     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
669     Handle(Geom_Line)  aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
670     Handle(Geom_Line)  aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 ));
671     if( aCirc.IsNull() || aLine1.IsNull() || aLine2.IsNull() )
672       return error(COMPERR_BAD_SHAPE);
673
674     if ( !algo1d->ComputeCircularEdge( aMesh, CircEdge ))
675       return error( algo1d->GetComputeError() );
676     map< double, const SMDS_MeshNode* > theNodes;
677     if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes))
678       return error("Circular edge is incorrectly meshed");
679
680     const SMDS_MeshNode* NF = theNodes.begin()->second;
681     const SMDS_MeshNode* NL = theNodes.rbegin()->second;
682     CNodes.clear();
683     CNodes.push_back( NF );
684     map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
685     double fang = (*itn).first;
686     itn++;
687     for(; itn != theNodes.end(); itn++ ) {
688       CNodes.push_back( (*itn).second );
689       double ang = (*itn).first - fang;
690       if( ang>M_PI ) ang = ang - 2.*M_PI;
691       if( ang<-M_PI ) ang = ang + 2.*M_PI;
692       Angles.Append( ang );
693     }
694     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
695     gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
696     P0 = aCirc->Location();
697
698     // make P1 belong to LinEdge1
699     TopoDS_Vertex V1 = myHelper->IthVertex( 0, LinEdge1 );
700     TopoDS_Vertex V2 = myHelper->IthVertex( 1, LinEdge1 );
701     gp_Pnt PE1 = BRep_Tool::Pnt(V1);
702     gp_Pnt PE2 = BRep_Tool::Pnt(V2);
703     if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
704         ( P1.Distance(PE2) > Precision::Confusion() ) )
705       std::swap( LinEdge1, LinEdge2 );
706
707     bool linEdge1Computed, linEdge2Computed;
708     if ( !computeLayerPositions(P0,P1,LinEdge1,&linEdge1Computed))
709       return false;
710
711     Nodes1.resize( myLayerPositions.size()+1 );
712     Nodes2.resize( myLayerPositions.size()+1 );
713
714     // check that both linear edges have same hypotheses
715     if ( !computeLayerPositions(P0,P2,LinEdge2, &linEdge2Computed))
716          return false;
717     if ( Nodes1.size() != myLayerPositions.size()+1 )
718       return error("Different hypotheses apply to radial edges");
719
720     // find the central vertex
721     TopoDS_Vertex VC = V2;
722     if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
723         ( P2.Distance(PE1) > Precision::Confusion() ) )
724       VC = V1;
725     int vertID = meshDS->ShapeToIndex(VC);
726
727     // LinEdge1
728     if ( linEdge1Computed )
729     {
730       if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
731         return error("Invalid mesh on a straight edge");
732
733       bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
734       NC = const_cast<SMDS_MeshNode*>
735         ( nodesFromP0ToP1 ? theNodes.begin()->second : theNodes.rbegin()->second );
736       int i = 0, ir = Nodes1.size()-1;
737       int * pi = nodesFromP0ToP1 ? &i : &ir;
738       itn = theNodes.begin();
739       if ( nodesFromP0ToP1 ) ++itn;
740       for ( ; i < Nodes1.size(); ++i, --ir, ++itn )
741       {
742         Nodes1[*pi] = itn->second;
743       }
744       for ( i = 0; i < Nodes1.size()-1; ++i )
745       {
746         Points.Append( gpXYZ( Nodes1[i]));
747         Pnts2d1.Append( myHelper->GetNodeUV( F, Nodes1[i]));
748       }
749     }
750     else
751     {
752       int edgeID = meshDS->ShapeToIndex(LinEdge1);
753       gp_Vec aVec(P0,P1);
754       // check orientation
755       Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
756       gp_Pnt Ptmp = Crv->Value(fp);
757       bool ori = false;
758       if( P1.Distance(Ptmp) > Precision::Confusion() )
759         ori = true;
760       // get UV points for edge
761       gp_Pnt2d PF,PL;
762       BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
763       gp_Vec2d V2d;
764       if(ori) {
765         V2d = gp_Vec2d(PF,PL);
766         PC = PF;
767       }
768       else {
769         V2d = gp_Vec2d(PL,PF);
770         PC = PL;
771       }
772       NC = const_cast<SMDS_MeshNode*>( VertexNode( VC, meshDS ));
773       if ( !NC )
774       {
775         NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
776         meshDS->SetNodeOnVertex(NC, vertID);
777       }
778       double dp = lp-fp;
779       int i = 0;
780       for(; i<myLayerPositions.size(); i++) {
781         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
782                   P0.Y() + aVec.Y()*myLayerPositions[i],
783                   P0.Z() + aVec.Z()*myLayerPositions[i] );
784         Points.Append(P);
785         SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
786         Nodes1[i] = node;
787         double param;
788         if(!ori)
789           param = fp + dp*(1-myLayerPositions[i]);
790         else
791           param = fp + dp*myLayerPositions[i];
792         meshDS->SetNodeOnEdge(node, edgeID, param);
793         // parameters on face
794         gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
795                       PC.Y() + V2d.Y()*myLayerPositions[i] );
796         Pnts2d1.Append(P2d);
797       }
798       Nodes1[ myLayerPositions.size() ] = NF;
799       // create 1D elements on edge
800       SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes1[0] );
801       if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
802       for(i=1; i<Nodes1.size(); i++) {
803         ME = myHelper->AddEdge( Nodes1[i-1], Nodes1[i] );
804         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
805       }
806       if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
807         Nodes2 = Nodes1;
808     }
809     markEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
810
811     // LinEdge2
812     if ( linEdge2Computed )
813     {
814       if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge2,true,theNodes))
815         return error("Invalid mesh on a straight edge");
816
817       bool nodesFromP0ToP2 = ( theNodes.rbegin()->second == NL );
818       int i = 0, ir = Nodes1.size()-1;
819       int * pi = nodesFromP0ToP2 ? &i : &ir;
820       itn = theNodes.begin();
821       if ( nodesFromP0ToP2 ) ++itn;
822       for ( ; i < Nodes2.size(); ++i, --ir, ++itn )
823         Nodes2[*pi] = itn->second;
824     }
825     else
826     {
827       int edgeID = meshDS->ShapeToIndex(LinEdge2);
828       gp_Vec aVec = gp_Vec(P0,P2);
829       // check orientation
830       Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge2,fp,lp);
831       gp_Pnt Ptmp = Crv->Value(fp);
832       bool ori = false;
833       if( P2.Distance(Ptmp) > Precision::Confusion() )
834         ori = true;
835       // get UV points for edge
836       gp_Pnt2d PF,PL;
837       BRep_Tool::UVPoints( LinEdge2, TopoDS::Face(aShape), PF, PL );
838       gp_Vec2d V2d;
839       if(ori) {
840         V2d = gp_Vec2d(PF,PL);
841         PC = PF;
842       }
843       else {
844         V2d = gp_Vec2d(PL,PF);
845         PC = PL;
846       }
847       double dp = lp-fp;
848       for(int i=0; i<myLayerPositions.size(); i++) {
849         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
850                   P0.Y() + aVec.Y()*myLayerPositions[i],
851                   P0.Z() + aVec.Z()*myLayerPositions[i] );
852         SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
853         Nodes2[i] = node;
854         double param;
855         if(!ori)
856           param = fp + dp*(1-myLayerPositions[i]);
857         else
858           param = fp + dp*myLayerPositions[i];
859         meshDS->SetNodeOnEdge(node, edgeID, param);
860         // parameters on face
861         gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
862                       PC.Y() + V2d.Y()*myLayerPositions[i] );
863       }
864       Nodes2[ myLayerPositions.size() ] = NL;
865       // create 1D elements on edge
866       SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes2[0] );
867       if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
868       for(int i=1; i<Nodes2.size(); i++) {
869         ME = myHelper->AddEdge( Nodes2[i-1], Nodes2[i] );
870         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
871       }
872     }
873     markEdgeAsComputedByMe( LinEdge2, aMesh.GetSubMesh( F ));
874   }
875   markEdgeAsComputedByMe( CircEdge, aMesh.GetSubMesh( F ));
876
877   // orientation
878   bool IsForward = ( CircEdge.Orientation()==TopAbs_FORWARD );
879   const double angleSign = ( F.Orientation() == TopAbs_REVERSED ? -1.0 : 1.0 );
880
881   // create nodes and mesh elements on face
882   // find axis of rotation
883   gp_Pnt P2 = gp_Pnt( CNodes[1]->X(), CNodes[1]->Y(), CNodes[1]->Z() );
884   gp_Vec Vec1(P0,P1);
885   gp_Vec Vec2(P0,P2);
886   gp_Vec Axis = Vec1.Crossed(Vec2);
887   // create elements
888   int i = 1;
889   //cout<<"Angles.Length() = "<<Angles.Length()<<"   Points.Length() = "<<Points.Length()<<endl;
890   //cout<<"Nodes1.size() = "<<Nodes1.size()<<"   Pnts2d1.Length() = "<<Pnts2d1.Length()<<endl;
891   for(; i<Angles.Length(); i++) {
892     vector< const SMDS_MeshNode* > tmpNodes;
893     tmpNodes.reserve(Nodes1.size());
894     gp_Trsf aTrsf;
895     gp_Ax1 theAxis(P0,gp_Dir(Axis));
896     aTrsf.SetRotation( theAxis, Angles.Value(i) );
897     gp_Trsf2d aTrsf2d;
898     aTrsf2d.SetRotation( PC, Angles.Value(i) * angleSign );
899     // create nodes
900     int j = 1;
901     for(; j<=Points.Length(); j++) {
902       double cx,cy,cz;
903       Points.Value(j).Coord( cx, cy, cz );
904       aTrsf.Transforms( cx, cy, cz );
905       SMDS_MeshNode* node = myHelper->AddNode( cx, cy, cz );
906       // find parameters on face
907       Pnts2d1.Value(j).Coord( cx, cy );
908       aTrsf2d.Transforms( cx, cy );
909       // set node on face
910       meshDS->SetNodeOnFace( node, faceID, cx, cy );
911       tmpNodes[j-1] = node;
912     }
913     // create faces
914     tmpNodes[Points.Length()] = CNodes[i];
915     // quad
916     for(j=0; j<Nodes1.size()-1; j++) {
917       SMDS_MeshFace* MF;
918       if(IsForward)
919         MF = myHelper->AddFace( tmpNodes[j], Nodes1[j],
920                                 Nodes1[j+1], tmpNodes[j+1] );
921       else
922         MF = myHelper->AddFace( tmpNodes[j], tmpNodes[j+1],
923                                 Nodes1[j+1], Nodes1[j] );
924       if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
925     }
926     // tria
927     SMDS_MeshFace* MF;
928     if(IsForward)
929       MF = myHelper->AddFace( NC, Nodes1[0], tmpNodes[0] );
930     else
931       MF = myHelper->AddFace( NC, tmpNodes[0], Nodes1[0] );
932     if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
933     for(j=0; j<Nodes1.size(); j++) {
934       Nodes1[j] = tmpNodes[j];
935     }
936   }
937   // create last faces
938   // quad
939   for(i=0; i<Nodes1.size()-1; i++) {
940     SMDS_MeshFace* MF;
941     if(IsForward)
942       MF = myHelper->AddFace( Nodes2[i], Nodes1[i],
943                               Nodes1[i+1], Nodes2[i+1] );
944     else
945       MF = myHelper->AddFace( Nodes2[i],  Nodes2[i+1],
946                               Nodes1[i+1], Nodes1[i] );
947     if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
948   }
949   // tria
950   SMDS_MeshFace* MF;
951   if(IsForward)
952     MF = myHelper->AddFace( NC, Nodes1[0], Nodes2[0] );
953   else
954     MF = myHelper->AddFace( NC, Nodes2[0], Nodes1[0] );
955   if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
956
957   return true;
958 }
959
960 //================================================================================
961 /*!
962  * \brief Compute positions of nodes on the radial edge
963   * \retval bool - is a success
964  */
965 //================================================================================
966
967 bool StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(const gp_Pnt&      p1,
968                                                              const gp_Pnt&      p2,
969                                                              const TopoDS_Edge& linEdge,
970                                                              bool*              linEdgeComputed)
971 {
972   // First, try to compute positions of layers
973
974   myLayerPositions.clear();
975
976   SMESH_Mesh * mesh = myHelper->GetMesh();
977
978   const SMESH_Hypothesis* hyp1D = myDistributionHypo ? myDistributionHypo->GetLayerDistribution() : 0;
979   int                  nbLayers = myNbLayerHypo ? myNbLayerHypo->GetNumberOfLayers() : 0;
980
981   if ( !hyp1D && !nbLayers )
982   {
983     // No own algo hypotheses assigned, so first try to find any 1D hypothesis.
984     // We need some edge
985     TopoDS_Shape edge = linEdge;
986     if ( edge.IsNull() && !myHelper->GetSubShape().IsNull())
987       for ( TopExp_Explorer e(myHelper->GetSubShape(), TopAbs_EDGE); e.More(); e.Next())
988         edge = e.Current();
989     if ( !edge.IsNull() )
990     {
991       // find a hyp usable by TNodeDistributor
992       SMESH_HypoFilter hypKind;
993       TNodeDistributor::GetDistributor(*mesh)->InitCompatibleHypoFilter(hypKind,/*ignoreAux=*/1);
994       hyp1D = mesh->GetHypothesis( edge, hypKind, /*fromAncestors=*/true);
995     }
996   }
997   if ( hyp1D ) // try to compute with hyp1D
998   {
999     if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions,p1,p2,*mesh,hyp1D )) {
1000       if ( myDistributionHypo ) { // bad hyp assigned 
1001         return error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
1002       }
1003       else {
1004         // bad hyp found, its Ok, lets try with default nb of segnents
1005       }
1006     }
1007   }
1008   
1009   if ( myLayerPositions.empty() ) // try to use nb of layers
1010   {
1011     if ( !nbLayers )
1012       nbLayers = _gen->GetDefaultNbSegments();
1013
1014     if ( nbLayers )
1015     {
1016       myLayerPositions.resize( nbLayers - 1 );
1017       for ( int z = 1; z < nbLayers; ++z )
1018         myLayerPositions[ z - 1 ] = double( z )/ double( nbLayers );
1019     }
1020   }
1021
1022   // Second, check presence of a mesh built by other algo on linEdge
1023   // and mesh conformity to my hypothesis
1024
1025   bool meshComputed = (!linEdge.IsNull() && !mesh->GetSubMesh(linEdge)->IsEmpty() );
1026   if ( linEdgeComputed ) *linEdgeComputed = meshComputed;
1027
1028   if ( meshComputed )
1029   {
1030     vector< double > nodeParams;
1031     GetNodeParamOnEdge( mesh->GetMeshDS(), linEdge, nodeParams );
1032
1033     // nb of present nodes must be different in cases of 1 and 2 straight edges
1034
1035     TopoDS_Vertex VV[2];
1036     TopExp::Vertices( linEdge, VV[0], VV[1]);
1037     const gp_Pnt* points[] = { &p1, &p2 };
1038     gp_Pnt       vPoints[] = { BRep_Tool::Pnt(VV[0]), BRep_Tool::Pnt(VV[1]) };
1039     const double     tol[] = { BRep_Tool::Tolerance(VV[0]), BRep_Tool::Tolerance(VV[1]) };
1040     bool pointsAreOnVertices = true;
1041     for ( int iP = 0; iP < 2 && pointsAreOnVertices; ++iP )
1042       pointsAreOnVertices = ( points[iP]->Distance( vPoints[0] ) < tol[0] ||
1043                               points[iP]->Distance( vPoints[1] ) < tol[1] );
1044
1045     int nbNodes = nodeParams.size() - 2; // 2 straight edges
1046     if ( !pointsAreOnVertices )
1047       nbNodes = ( nodeParams.size() - 3 ) / 2; // 1 straight edge
1048
1049     if ( myLayerPositions.empty() )
1050     {
1051       myLayerPositions.resize( nbNodes );
1052     }
1053     else if ( myDistributionHypo || myNbLayerHypo )
1054     {
1055       // linEdge is computed by other algo. Check if there is a meshed face
1056       // using nodes on linEdge
1057       bool nodesAreUsed = false;
1058       TopTools_ListIteratorOfListOfShape ancestIt = mesh->GetAncestors( linEdge );
1059       for ( ; ancestIt.More() && !nodesAreUsed; ancestIt.Next() )
1060         if ( ancestIt.Value().ShapeType() == TopAbs_FACE )
1061           nodesAreUsed = (!mesh->GetSubMesh( ancestIt.Value() )->IsEmpty());
1062       if ( !nodesAreUsed ) {
1063         // rebuild them
1064         mesh->GetSubMesh( linEdge )->ComputeStateEngine( SMESH_subMesh::CLEAN );
1065         if ( linEdgeComputed ) *linEdgeComputed = false;
1066       }
1067       else {
1068         
1069         if ( myLayerPositions.size() != nbNodes )
1070           return error("Radial edge is meshed by other algorithm");
1071       }
1072     }
1073   }
1074
1075   return !myLayerPositions.empty();
1076 }
1077
1078
1079 //=======================================================================
1080 //function : Evaluate
1081 //purpose  : 
1082 //=======================================================================
1083
1084 bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
1085                                                 const TopoDS_Shape& aShape,
1086                                                 MapShapeNbElems& aResMap)
1087 {
1088   if( aShape.ShapeType() != TopAbs_FACE ) {
1089     return false;
1090   }
1091   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1092   if( aResMap.count(sm) )
1093     return false;
1094
1095   vector<int>& aResVec =
1096     aResMap.insert( make_pair(sm, vector<int>(SMDSEntity_Last,0))).first->second;
1097
1098   myHelper = new SMESH_MesherHelper( aMesh );
1099   myHelper->SetSubShape( aShape );
1100   auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
1101
1102   TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
1103
1104   TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
1105   int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 );
1106   if( nbe>3 || nbe < 1 || CircEdge.IsNull() )
1107     return false;
1108
1109   Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
1110   if( aCirc.IsNull() )
1111     return error(COMPERR_BAD_SHAPE);
1112
1113   gp_Pnt P0 = aCirc->Location();
1114   gp_Pnt P1 = aCirc->Value(0.);
1115   computeLayerPositions( P0, P1, LinEdge1 );
1116
1117   int nb0d=0, nb2d_tria=0, nb2d_quad=0;
1118   bool isQuadratic = false, ok = true;
1119   if(nbe==1)
1120   {
1121     // C1 must be a circle
1122     ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap );
1123     if(ok) {
1124       const vector<int>& aVec = aResMap[aMesh.GetSubMesh(CircEdge)];
1125       isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1126       if(isQuadratic) {
1127         // main nodes
1128         nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1129         // radial medium nodes
1130         nb0d += (aVec[SMDSEntity_Node]+1) * (myLayerPositions.size()+1);
1131         // other medium nodes
1132         nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1133       }
1134       else {
1135         nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1136       }
1137       nb2d_tria = aVec[SMDSEntity_Node] + 1;
1138       nb2d_quad = nb0d;
1139     }
1140   }
1141   else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL)
1142   {
1143     // one curve must be a half of circle and other curve must be
1144     // a segment of line
1145     double fp, lp;
1146     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
1147     if( fabs(fabs(lp-fp)-M_PI) > Precision::Confusion() ) {
1148       // not half of circle
1149       return error(COMPERR_BAD_SHAPE);
1150     }
1151     Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
1152     if( aLine.IsNull() ) {
1153       // other curve not line
1154       return error(COMPERR_BAD_SHAPE);
1155     }
1156     ok = !aResMap.count( aMesh.GetSubMesh(LinEdge1) );
1157     if ( !ok ) {
1158       const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(LinEdge1) ];
1159       ok = ( aVec[SMDSEntity_Node] == myLayerPositions.size() );
1160     }
1161     if(ok) {
1162       ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap );
1163     }
1164     if(ok) {
1165       const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(CircEdge) ];
1166       isQuadratic = aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge];
1167       if(isQuadratic) {
1168         // main nodes
1169         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1170         // radial medium nodes
1171         nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1172         // other medium nodes
1173         nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1174       }
1175       else {
1176         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1177       }
1178       nb2d_tria = aVec[SMDSEntity_Node] + 1;
1179       nb2d_quad = nb2d_tria * myLayerPositions.size();
1180       // add evaluation for edges
1181       vector<int> aResVec(SMDSEntity_Last,0);
1182       if(isQuadratic) {
1183         aResVec[SMDSEntity_Node] = 4*myLayerPositions.size() + 3;
1184         aResVec[SMDSEntity_Quad_Edge] = 2*myLayerPositions.size() + 2;
1185       }
1186       else {
1187         aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1188         aResVec[SMDSEntity_Edge] = 2*myLayerPositions.size() + 2;
1189       }
1190       aResMap[ aMesh.GetSubMesh(LinEdge1) ] = aResVec;
1191     }
1192   }
1193   else  // nbe==3 or ( nbe==2 && linEdge is INTERNAL )
1194   {
1195     if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
1196       LinEdge2 = LinEdge1;
1197
1198     // one curve must be a part of circle and other curves must be
1199     // segments of line
1200     Handle(Geom_Line)  aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
1201     Handle(Geom_Line)  aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 ));
1202     if( aLine1.IsNull() || aLine2.IsNull() ) {
1203       // other curve not line
1204       return error(COMPERR_BAD_SHAPE);
1205     }
1206     int nbLayers = myLayerPositions.size();
1207     computeLayerPositions( P0, P1, LinEdge2 );
1208     if ( nbLayers != myLayerPositions.size() )
1209       return error("Different hypotheses apply to radial edges");
1210       
1211     bool ok = !aResMap.count( aMesh.GetSubMesh(LinEdge1));
1212     if ( !ok ) {
1213       if ( myDistributionHypo || myNbLayerHypo )
1214         ok = true; // override other 1d hyps
1215       else {
1216         const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(LinEdge1) ];
1217         ok = ( aVec[SMDSEntity_Node] == myLayerPositions.size() );
1218       }
1219     }
1220     if( ok && aResMap.count( aMesh.GetSubMesh(LinEdge2) )) {
1221       if ( myDistributionHypo || myNbLayerHypo )
1222         ok = true; // override other 1d hyps
1223       else {
1224         const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(LinEdge2) ];
1225         ok = ( aVec[SMDSEntity_Node] == myLayerPositions.size() );
1226       }
1227     }
1228     if(ok) {
1229       ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap );
1230     }
1231     if(ok) {
1232       const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(CircEdge) ];
1233       isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1234       if(isQuadratic) {
1235         // main nodes
1236         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1237         // radial medium nodes
1238         nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1239         // other medium nodes
1240         nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1241       }
1242       else {
1243         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1244       }
1245       nb2d_tria = aVec[SMDSEntity_Node] + 1;
1246       nb2d_quad = nb2d_tria * myLayerPositions.size();
1247       // add evaluation for edges
1248       vector<int> aResVec(SMDSEntity_Last, 0);
1249       if(isQuadratic) {
1250         aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1251         aResVec[SMDSEntity_Quad_Edge] = myLayerPositions.size() + 1;
1252       }
1253       else {
1254         aResVec[SMDSEntity_Node] = myLayerPositions.size();
1255         aResVec[SMDSEntity_Edge] = myLayerPositions.size() + 1;
1256       }
1257       sm = aMesh.GetSubMesh(LinEdge1);
1258       aResMap[sm] = aResVec;
1259       sm = aMesh.GetSubMesh(LinEdge2);
1260       aResMap[sm] = aResVec;
1261     }
1262   }
1263
1264   if(nb0d>0) {
1265     aResVec[0] = nb0d;
1266     if(isQuadratic) {
1267       aResVec[SMDSEntity_Quad_Triangle] = nb2d_tria;
1268       aResVec[SMDSEntity_Quad_Quadrangle] = nb2d_quad;
1269     }
1270     else {
1271       aResVec[SMDSEntity_Triangle] = nb2d_tria;
1272       aResVec[SMDSEntity_Quadrangle] = nb2d_quad;
1273     }
1274     return true;
1275   }
1276
1277   // invalid case
1278   sm = aMesh.GetSubMesh(aShape);
1279   SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1280   smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
1281                                         "Submesh can not be evaluated",this));
1282   return false;
1283
1284 }