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