Back to index

salome-smesh  6.5.0
StdMeshers_RadialQuadrangle_1D2D.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 //  SMESH SMESH : implementaion of SMESH idl descriptions
00021 // File      : StdMeshers_RadialQuadrangle_1D2D.cxx
00022 // Module    : SMESH
00023 // Created   : Fri Oct 20 11:37:07 2006
00024 // Author    : Edward AGAPOV (eap)
00025 
00026 #include "StdMeshers_RadialQuadrangle_1D2D.hxx"
00027 
00028 #include "StdMeshers_NumberOfLayers.hxx"
00029 #include "StdMeshers_LayerDistribution.hxx"
00030 #include "StdMeshers_Regular_1D.hxx"
00031 #include "StdMeshers_NumberOfSegments.hxx"
00032 
00033 #include "SMDS_MeshNode.hxx"
00034 #include "SMESHDS_SubMesh.hxx"
00035 #include "SMESH_Gen.hxx"
00036 #include "SMESH_HypoFilter.hxx"
00037 #include "SMESH_Mesh.hxx"
00038 #include "SMESH_MesherHelper.hxx"
00039 #include "SMESH_subMesh.hxx"
00040 #include "SMESH_subMeshEventListener.hxx"
00041 
00042 #include "utilities.h"
00043 
00044 #include <BRepAdaptor_Curve.hxx>
00045 #include <BRepBuilderAPI_MakeEdge.hxx>
00046 #include <BRep_Tool.hxx>
00047 #include <GeomAPI_ProjectPointOnSurf.hxx>
00048 #include <Geom_Circle.hxx>
00049 #include <Geom_Line.hxx>
00050 #include <Geom_TrimmedCurve.hxx>
00051 #include <TColgp_SequenceOfPnt.hxx>
00052 #include <TColgp_SequenceOfPnt2d.hxx>
00053 #include <TopExp.hxx>
00054 #include <TopExp_Explorer.hxx>
00055 #include <TopTools_ListIteratorOfListOfShape.hxx>
00056 #include <TopoDS.hxx>
00057 
00058 
00059 using namespace std;
00060 
00061 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
00062 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
00063 
00064 
00065 //=======================================================================
00066 //function : StdMeshers_RadialQuadrangle_1D2D
00067 //purpose  : 
00068 //=======================================================================
00069 
00070 StdMeshers_RadialQuadrangle_1D2D::StdMeshers_RadialQuadrangle_1D2D(int hypId,
00071                                                                    int studyId,
00072                                                                    SMESH_Gen* gen)
00073   :SMESH_2D_Algo(hypId, studyId, gen)
00074 {
00075   _name = "RadialQuadrangle_1D2D";
00076   _shapeType = (1 << TopAbs_FACE);        // 1 bit per shape type
00077 
00078   _compatibleHypothesis.push_back("LayerDistribution2D");
00079   _compatibleHypothesis.push_back("NumberOfLayers2D");
00080   myNbLayerHypo = 0;
00081   myDistributionHypo = 0;
00082   _requireDiscreteBoundary = false;
00083   _supportSubmeshes = true;
00084 }
00085 
00086 
00087 //================================================================================
00091 //================================================================================
00092 
00093 StdMeshers_RadialQuadrangle_1D2D::~StdMeshers_RadialQuadrangle_1D2D()
00094 {}
00095 
00096 
00097 //=======================================================================
00098 //function : CheckHypothesis
00099 //purpose  : 
00100 //=======================================================================
00101 
00102 bool StdMeshers_RadialQuadrangle_1D2D::CheckHypothesis
00103                            (SMESH_Mesh&                          aMesh,
00104                             const TopoDS_Shape&                  aShape,
00105                             SMESH_Hypothesis::Hypothesis_Status& aStatus)
00106 {
00107   // check aShape 
00108   myNbLayerHypo = 0;
00109   myDistributionHypo = 0;
00110 
00111   list <const SMESHDS_Hypothesis * >::const_iterator itl;
00112 
00113   const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
00114   if ( hyps.size() == 0 ) {
00115     aStatus = SMESH_Hypothesis::HYP_OK;
00116     return true;  // can work with no hypothesis
00117   }
00118 
00119   if ( hyps.size() > 1 ) {
00120     aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
00121     return false;
00122   }
00123 
00124   const SMESHDS_Hypothesis *theHyp = hyps.front();
00125 
00126   string hypName = theHyp->GetName();
00127 
00128   if (hypName == "NumberOfLayers2D") {
00129     myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
00130     aStatus = SMESH_Hypothesis::HYP_OK;
00131     return true;
00132   }
00133   if (hypName == "LayerDistribution2D") {
00134     myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
00135     aStatus = SMESH_Hypothesis::HYP_OK;
00136     return true;
00137   }
00138   aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
00139   return true;
00140 }
00141 
00142 namespace
00143 {
00144   // ------------------------------------------------------------------------------
00148   class TEdgeMarker : public SMESH_subMeshEventListener
00149   {
00150     TEdgeMarker(): SMESH_subMeshEventListener(/*isDeletable=*/false,
00151                                               "StdMeshers_RadialQuadrangle_1D2D::TEdgeMarker") {}
00152   public:
00154     static SMESH_subMeshEventListener* getListener()
00155     {
00156       static TEdgeMarker theEdgeMarker;
00157       return &theEdgeMarker;
00158     }
00160     void ProcessEvent(const int          event,
00161                       const int          eventType,
00162                       SMESH_subMesh*     edgeSubMesh,
00163                       EventListenerData* data,
00164                       const SMESH_Hypothesis*  /*hyp*/)
00165     {
00166       if ( data && !data->mySubMeshes.empty() && eventType == SMESH_subMesh::ALGO_EVENT)
00167       {
00168         ASSERT( data->mySubMeshes.front() != edgeSubMesh );
00169         SMESH_subMesh* faceSubMesh = data->mySubMeshes.front();
00170         faceSubMesh->ComputeStateEngine( SMESH_subMesh::CLEAN );
00171       }
00172     }
00173   };
00174 
00175   // ------------------------------------------------------------------------------
00179   void markEdgeAsComputedByMe(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
00180   {
00181     if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
00182     {
00183       if ( !edgeSM->GetEventListenerData( TEdgeMarker::getListener() ))
00184         faceSubMesh->SetEventListener( TEdgeMarker::getListener(),
00185                                        SMESH_subMeshEventListenerData::MakeData(faceSubMesh),
00186                                        edgeSM);
00187     }
00188   }
00189   // ------------------------------------------------------------------------------
00194 //   bool isEdgeCompatiballyMeshed(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
00195 //   {
00196 //     if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
00197 //     {
00198 //       if ( SMESH_subMeshEventListenerData* otherFaceData =
00199 //            edgeSM->GetEventListenerData( TEdgeMarker::getListener() ))
00200 //       {
00201 //         // compare hypothesis aplied to two disk faces sharing radial edges
00202 //         SMESH_Mesh& mesh = *faceSubMesh->GetFather();
00203 //         SMESH_Algo* radialQuadAlgo = mesh.GetGen()->GetAlgo(mesh, faceSubMesh->GetSubShape() );
00204 //         SMESH_subMesh* otherFaceSubMesh = otherFaceData->mySubMeshes.front();
00205 //         list <const SMESHDS_Hypothesis *> hyps1 =
00206 //           radialQuadAlgo->GetUsedHypothesis( mesh, faceSubMesh->GetSubShape());
00207 //         list <const SMESHDS_Hypothesis *> hyps2 =
00208 //           radialQuadAlgo->GetUsedHypothesis( mesh, otherFaceSubMesh->GetSubShape());
00209 //         if( hyps1.empty() && hyps2.empty() )
00210 //           return true; // defaul hyps
00211 //         if ( hyps1.size() != hyps2.size() )
00212 //           return false;
00213 //         return *hyps1.front() == *hyps2.front();
00214 //       }
00215 //     }
00216 //     return false;
00217 //   }
00218 
00219   //================================================================================
00223   //================================================================================
00224 
00225   Handle(Geom_Curve) getCurve(const TopoDS_Edge& edge, double* f=0, double* l=0)
00226   {
00227     Handle(Geom_Curve) C;
00228     if ( !edge.IsNull() )
00229     {
00230       double first = 0., last = 0.;
00231       C = BRep_Tool::Curve(edge, first, last);
00232       if ( !C.IsNull() )
00233       {
00234         Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C);
00235         while( !tc.IsNull() ) {
00236           C = tc->BasisCurve();
00237           tc = Handle(Geom_TrimmedCurve)::DownCast(C);
00238         }
00239         if ( f ) *f = first;
00240         if ( l ) *l = last;
00241       }
00242     }
00243     return C;
00244   }
00245 
00246   //================================================================================
00251   //================================================================================
00252 
00253   int analyseFace(const TopoDS_Shape& face,
00254                   TopoDS_Edge&        CircEdge,
00255                   TopoDS_Edge&        LinEdge1,
00256                   TopoDS_Edge&        LinEdge2)
00257   {
00258     CircEdge.Nullify(); LinEdge1.Nullify(); LinEdge2.Nullify();
00259     int nbe = 0;
00260 
00261     for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next(), ++nbe )
00262     {
00263       const TopoDS_Edge& E = TopoDS::Edge( exp.Current() );
00264       double f,l;
00265       Handle(Geom_Curve) C = getCurve(E,&f,&l);
00266       if ( !C.IsNull() )
00267       {
00268         if ( C->IsKind( STANDARD_TYPE(Geom_Circle)))
00269         {
00270           if ( CircEdge.IsNull() )
00271             CircEdge = E;
00272           else
00273             return 0;
00274         }
00275         else if ( LinEdge1.IsNull() )
00276           LinEdge1 = E;
00277         else
00278           LinEdge2 = E;
00279       }
00280     }
00281     return nbe;
00282   }
00283 
00284 //================================================================================
00285 //================================================================================
00290 //================================================================================
00291 //================================================================================
00292 
00293 class TNodeDistributor: public StdMeshers_Regular_1D
00294 {
00295   list <const SMESHDS_Hypothesis *> myUsedHyps;
00296 public:
00297   // -----------------------------------------------------------------------------
00298   static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
00299   {
00300     const int myID = -1000;
00301     map < int, SMESH_1D_Algo * > & algoMap = aMesh.GetGen()->_map1D_Algo;
00302     map < int, SMESH_1D_Algo * >::iterator id_algo = algoMap.find( myID );
00303     if ( id_algo == algoMap.end() )
00304       return new TNodeDistributor( myID, 0, aMesh.GetGen() );
00305     return static_cast< TNodeDistributor* >( id_algo->second );
00306   }
00307   // -----------------------------------------------------------------------------
00309   bool Compute( vector< double > &      positions,
00310                 gp_Pnt                  pIn,
00311                 gp_Pnt                  pOut,
00312                 SMESH_Mesh&             aMesh,
00313                 const SMESH_Hypothesis* hyp1d)
00314   {
00315     if ( !hyp1d ) return error( "Invalid LayerDistribution hypothesis");
00316 
00317     double len = pIn.Distance( pOut );
00318     if ( len <= DBL_MIN ) return error("Too close points of inner and outer shells");
00319 
00320     myUsedHyps.clear();
00321     myUsedHyps.push_back( hyp1d );
00322 
00323     TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut );
00324     SMESH_Hypothesis::Hypothesis_Status aStatus;
00325     if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
00326       return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
00327                     "with LayerDistribution hypothesis");
00328 
00329     BRepAdaptor_Curve C3D(edge);
00330     double f = C3D.FirstParameter(), l = C3D.LastParameter();
00331     list< double > params;
00332     if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false ))
00333       return error("StdMeshers_Regular_1D failed to compute layers distribution");
00334 
00335     positions.clear();
00336     positions.reserve( params.size() );
00337     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
00338       positions.push_back( *itU / len );
00339     return true;
00340   }
00341   // -----------------------------------------------------------------------------
00343   bool ComputeCircularEdge(SMESH_Mesh&         aMesh,
00344                            const TopoDS_Edge& anEdge)
00345   {
00346     _gen->Compute( aMesh, anEdge);
00347     SMESH_subMesh *sm = aMesh.GetSubMesh(anEdge);
00348     if ( sm->GetComputeState() != SMESH_subMesh::COMPUTE_OK)
00349     {
00350       // find any 1d hyp assigned (there can be a hyp w/o algo)
00351       myUsedHyps = SMESH_Algo::GetUsedHypothesis(aMesh, anEdge, /*ignoreAux=*/true);
00352       Hypothesis_Status aStatus;
00353       if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, anEdge, aStatus ))
00354       {
00355         // no valid 1d hyp assigned, use default nb of segments
00356         _hypType                    = NB_SEGMENTS;
00357         _ivalue[ DISTR_TYPE_IND ]   = StdMeshers_NumberOfSegments::DT_Regular;
00358         _ivalue[ NB_SEGMENTS_IND  ] = _gen->GetDefaultNbSegments();
00359       }
00360       return StdMeshers_Regular_1D::Compute( aMesh, anEdge );
00361     }
00362     return true;
00363   }
00364   // -----------------------------------------------------------------------------
00366   bool EvaluateCircularEdge(SMESH_Mesh&        aMesh,
00367                             const TopoDS_Edge& anEdge,
00368                             MapShapeNbElems&   aResMap)
00369   {
00370     _gen->Evaluate( aMesh, anEdge, aResMap );
00371     if ( aResMap.count( aMesh.GetSubMesh( anEdge )))
00372       return true;
00373 
00374     // find any 1d hyp assigned
00375     myUsedHyps = SMESH_Algo::GetUsedHypothesis(aMesh, anEdge, /*ignoreAux=*/true);
00376     Hypothesis_Status aStatus;
00377     if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, anEdge, aStatus ))
00378     {
00379       // no valid 1d hyp assigned, use default nb of segments
00380       _hypType                    = NB_SEGMENTS;
00381       _ivalue[ DISTR_TYPE_IND ]   = StdMeshers_NumberOfSegments::DT_Regular;
00382       _ivalue[ NB_SEGMENTS_IND  ] = _gen->GetDefaultNbSegments();
00383     }
00384     return StdMeshers_Regular_1D::Evaluate( aMesh, anEdge, aResMap );
00385   }
00386 protected:
00387   // -----------------------------------------------------------------------------
00388   TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen)
00389     : StdMeshers_Regular_1D( hypId, studyId, gen)
00390   {
00391   }
00392   // -----------------------------------------------------------------------------
00393   virtual const list <const SMESHDS_Hypothesis *> &
00394     GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
00395   {
00396     return myUsedHyps;
00397   }
00398   // -----------------------------------------------------------------------------
00399 };
00400 }
00401 
00402 //=======================================================================
00409 //=======================================================================
00410 
00411 void StdMeshers_RadialQuadrangle_1D2D::SubmeshRestored(SMESH_subMesh* faceSubMesh)
00412 {
00413   if ( !faceSubMesh->IsEmpty() )
00414   {
00415     TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
00416     analyseFace( faceSubMesh->GetSubShape(), CircEdge, LinEdge1, LinEdge2 );
00417     if ( !CircEdge.IsNull() ) markEdgeAsComputedByMe( CircEdge, faceSubMesh );
00418     if ( !LinEdge1.IsNull() ) markEdgeAsComputedByMe( LinEdge1, faceSubMesh );
00419     if ( !LinEdge2.IsNull() ) markEdgeAsComputedByMe( LinEdge2, faceSubMesh );
00420   }
00421 }
00422 
00423 //=======================================================================
00424 //function : Compute
00425 //purpose  : 
00426 //=======================================================================
00427 
00428 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
00429                                                const TopoDS_Shape& aShape)
00430 {
00431   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
00432 
00433   myHelper = new SMESH_MesherHelper( aMesh );
00434   myHelper->IsQuadraticSubMesh( aShape );
00435   // to delete helper at exit from Compute()
00436   auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
00437 
00438   TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
00439 
00440   TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
00441   int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 );
00442   Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
00443   if( nbe>3 || nbe < 1 || aCirc.IsNull() )
00444     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)");
00445   
00446   gp_Pnt P0, P1;
00447   // points for rotation
00448   TColgp_SequenceOfPnt Points;
00449   // angles for rotation
00450   TColStd_SequenceOfReal Angles;
00451   // Nodes1 and Nodes2 - nodes along radiuses
00452   // CNodes - nodes on circle edge
00453   vector< const SMDS_MeshNode* > Nodes1, Nodes2, CNodes;
00454   SMDS_MeshNode * NC;
00455   // parameters edge nodes on face
00456   TColgp_SequenceOfPnt2d Pnts2d1;
00457   gp_Pnt2d PC;
00458 
00459   int faceID = meshDS->ShapeToIndex(aShape);
00460   TopoDS_Face F = TopoDS::Face(aShape);
00461   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
00462 
00463 
00464   if(nbe==1)
00465   {
00466     if (!algo1d->ComputeCircularEdge( aMesh, CircEdge ))
00467       return error( algo1d->GetComputeError() );
00468     map< double, const SMDS_MeshNode* > theNodes;
00469     if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes))
00470       return error("Circular edge is incorrectly meshed");
00471 
00472     CNodes.clear();
00473     map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
00474     const SMDS_MeshNode* NF = (*itn).second;
00475     CNodes.push_back( (*itn).second );
00476     double fang = (*itn).first;
00477     if ( itn != theNodes.end() ) {
00478       itn++;
00479       for(; itn != theNodes.end(); itn++ ) {
00480         CNodes.push_back( (*itn).second );
00481         double ang = (*itn).first - fang;
00482         if( ang>M_PI ) ang = ang - 2.*M_PI;
00483         if( ang<-M_PI ) ang = ang + 2.*M_PI;
00484         Angles.Append( ang ); 
00485       }
00486     }
00487     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
00488     P0 = aCirc->Location();
00489 
00490     if ( !computeLayerPositions(P0,P1))
00491       return false;
00492 
00493     TopoDS_Vertex V1 = myHelper->IthVertex(0, CircEdge );
00494     gp_Pnt2d p2dV = BRep_Tool::Parameters( V1, TopoDS::Face(aShape) );
00495 
00496     NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
00497     GeomAPI_ProjectPointOnSurf PPS(P0,S);
00498     double U0,V0;
00499     PPS.Parameters(1,U0,V0);
00500     meshDS->SetNodeOnFace(NC, faceID, U0, V0);
00501     PC = gp_Pnt2d(U0,V0);
00502 
00503     gp_Vec aVec(P0,P1);
00504     gp_Vec2d aVec2d(PC,p2dV);
00505     Nodes1.resize( myLayerPositions.size()+1 );
00506     Nodes2.resize( myLayerPositions.size()+1 );
00507     int i = 0;
00508     for(; i<myLayerPositions.size(); i++) {
00509       gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
00510                 P0.Y() + aVec.Y()*myLayerPositions[i],
00511                 P0.Z() + aVec.Z()*myLayerPositions[i] );
00512       Points.Append(P);
00513       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
00514       Nodes1[i] = node;
00515       Nodes2[i] = node;
00516       double U = PC.X() + aVec2d.X()*myLayerPositions[i];
00517       double V = PC.Y() + aVec2d.Y()*myLayerPositions[i];
00518       meshDS->SetNodeOnFace( node, faceID, U, V );
00519       Pnts2d1.Append(gp_Pnt2d(U,V));
00520     }
00521     Nodes1[Nodes1.size()-1] = NF;
00522     Nodes2[Nodes1.size()-1] = NF;
00523   }
00524   else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL )
00525   {
00526     // one curve must be a half of circle and other curve must be
00527     // a segment of line
00528     double fp, lp;
00529     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
00530     if( fabs(fabs(lp-fp)-M_PI) > Precision::Confusion() ) {
00531       // not half of circle
00532       return error(COMPERR_BAD_SHAPE);
00533     }
00534     Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
00535     if( aLine.IsNull() ) {
00536       // other curve not line
00537       return error(COMPERR_BAD_SHAPE);
00538     }
00539 
00540     if ( !algo1d->ComputeCircularEdge( aMesh, CircEdge ))
00541       return error( algo1d->GetComputeError() );
00542     map< double, const SMDS_MeshNode* > theNodes;
00543     if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes) )
00544       return error("Circular edge is incorrectly meshed");
00545 
00546     map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
00547     CNodes.clear();
00548     CNodes.push_back( itn->second );
00549     double fang = (*itn).first;
00550     itn++;
00551     for(; itn != theNodes.end(); itn++ ) {
00552       CNodes.push_back( (*itn).second );
00553       double ang = (*itn).first - fang;
00554       if( ang>M_PI ) ang = ang - 2.*M_PI;
00555       if( ang<-M_PI ) ang = ang + 2.*M_PI;
00556       Angles.Append( ang );
00557     }
00558     const SMDS_MeshNode* NF = theNodes.begin()->second;
00559     const SMDS_MeshNode* NL = theNodes.rbegin()->second;
00560     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
00561     gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
00562     P0 = aCirc->Location();
00563 
00564     bool linEdgeComputed;
00565     if ( !computeLayerPositions(P0,P1,LinEdge1,&linEdgeComputed))
00566       return false;
00567 
00568     if ( linEdgeComputed )
00569     {
00570       if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
00571         return error("Invalid mesh on a straight edge");
00572 
00573       Nodes1.resize( myLayerPositions.size()+1 );
00574       Nodes2.resize( myLayerPositions.size()+1 );
00575       vector< const SMDS_MeshNode* > *pNodes1 = &Nodes1, *pNodes2 = &Nodes2;
00576       bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
00577       if ( !nodesFromP0ToP1 ) std::swap( pNodes1, pNodes2 );
00578 
00579       map< double, const SMDS_MeshNode* >::reverse_iterator ritn = theNodes.rbegin();
00580       itn = theNodes.begin();
00581       for ( int i = Nodes1.size()-1; i > -1; ++itn, ++ritn, --i )
00582       {
00583         (*pNodes1)[i] = ritn->second;
00584         (*pNodes2)[i] =  itn->second;
00585         Points.Prepend( gpXYZ( Nodes1[i]));
00586         Pnts2d1.Prepend( myHelper->GetNodeUV( F, Nodes1[i]));
00587       }
00588       NC = const_cast<SMDS_MeshNode*>( itn->second );
00589       Points.Remove( Nodes1.size() );
00590     }
00591     else
00592     {
00593       gp_Vec aVec(P0,P1);
00594       int edgeID = meshDS->ShapeToIndex(LinEdge1);
00595       // check orientation
00596       Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
00597       gp_Pnt Ptmp;
00598       Crv->D0(fp,Ptmp);
00599       bool ori = true;
00600       if( P1.Distance(Ptmp) > Precision::Confusion() )
00601         ori = false;
00602       // get UV points for edge
00603       gp_Pnt2d PF,PL;
00604       BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
00605       PC = gp_Pnt2d( (PF.X()+PL.X())/2, (PF.Y()+PL.Y())/2 );
00606       gp_Vec2d V2d;
00607       if(ori) V2d = gp_Vec2d(PC,PF);
00608       else V2d = gp_Vec2d(PC,PL);
00609       // add nodes on edge
00610       double cp = (fp+lp)/2;
00611       double dp2 = (lp-fp)/2;
00612       NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
00613       meshDS->SetNodeOnEdge(NC, edgeID, cp);
00614       Nodes1.resize( myLayerPositions.size()+1 );
00615       Nodes2.resize( myLayerPositions.size()+1 );
00616       int i = 0;
00617       for(; i<myLayerPositions.size(); i++) {
00618         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
00619                   P0.Y() + aVec.Y()*myLayerPositions[i],
00620                   P0.Z() + aVec.Z()*myLayerPositions[i] );
00621         Points.Append(P);
00622         SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
00623         Nodes1[i] = node;
00624         double param;
00625         if(ori)
00626           param = fp + dp2*(1-myLayerPositions[i]);
00627         else
00628           param = cp + dp2*myLayerPositions[i];
00629         meshDS->SetNodeOnEdge(node, edgeID, param);
00630         P = gp_Pnt( P0.X() - aVec.X()*myLayerPositions[i],
00631                     P0.Y() - aVec.Y()*myLayerPositions[i],
00632                     P0.Z() - aVec.Z()*myLayerPositions[i] );
00633         node = meshDS->AddNode(P.X(), P.Y(), P.Z());
00634         Nodes2[i] = node;
00635         if(!ori)
00636           param = fp + dp2*(1-myLayerPositions[i]);
00637         else
00638           param = cp + dp2*myLayerPositions[i];
00639         meshDS->SetNodeOnEdge(node, edgeID, param);
00640         // parameters on face
00641         gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
00642                       PC.Y() + V2d.Y()*myLayerPositions[i] );
00643         Pnts2d1.Append(P2d);
00644       }
00645       Nodes1[ myLayerPositions.size() ] = NF;
00646       Nodes2[ myLayerPositions.size() ] = NL;
00647       // create 1D elements on edge
00648       vector< const SMDS_MeshNode* > tmpNodes;
00649       tmpNodes.resize(2*Nodes1.size()+1);
00650       for(i=0; i<Nodes2.size(); i++)
00651         tmpNodes[Nodes2.size()-i-1] = Nodes2[i];
00652       tmpNodes[Nodes2.size()] = NC;
00653       for(i=0; i<Nodes1.size(); i++)
00654         tmpNodes[Nodes2.size()+1+i] = Nodes1[i];
00655       for(i=1; i<tmpNodes.size(); i++) {
00656         SMDS_MeshEdge* ME = myHelper->AddEdge( tmpNodes[i-1], tmpNodes[i] );
00657         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
00658       }
00659       markEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
00660     }
00661   }
00662   else // nbe==3 or ( nbe==2 && linEdge is INTERNAL )
00663   {
00664     if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
00665       LinEdge2 = LinEdge1;
00666 
00667     // one curve must be a part of circle and other curves must be
00668     // segments of line
00669     double fp, lp;
00670     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
00671     Handle(Geom_Line)  aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
00672     Handle(Geom_Line)  aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 ));
00673     if( aCirc.IsNull() || aLine1.IsNull() || aLine2.IsNull() )
00674       return error(COMPERR_BAD_SHAPE);
00675 
00676     if ( !algo1d->ComputeCircularEdge( aMesh, CircEdge ))
00677       return error( algo1d->GetComputeError() );
00678     map< double, const SMDS_MeshNode* > theNodes;
00679     if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes))
00680       return error("Circular edge is incorrectly meshed");
00681 
00682     const SMDS_MeshNode* NF = theNodes.begin()->second;
00683     const SMDS_MeshNode* NL = theNodes.rbegin()->second;
00684     CNodes.clear();
00685     CNodes.push_back( NF );
00686     map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
00687     double fang = (*itn).first;
00688     itn++;
00689     for(; itn != theNodes.end(); itn++ ) {
00690       CNodes.push_back( (*itn).second );
00691       double ang = (*itn).first - fang;
00692       if( ang>M_PI ) ang = ang - 2.*M_PI;
00693       if( ang<-M_PI ) ang = ang + 2.*M_PI;
00694       Angles.Append( ang );
00695     }
00696     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
00697     gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
00698     P0 = aCirc->Location();
00699 
00700     // make P1 belong to LinEdge1
00701     TopoDS_Vertex V1 = myHelper->IthVertex( 0, LinEdge1 );
00702     TopoDS_Vertex V2 = myHelper->IthVertex( 1, LinEdge1 );
00703     gp_Pnt PE1 = BRep_Tool::Pnt(V1);
00704     gp_Pnt PE2 = BRep_Tool::Pnt(V2);
00705     if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
00706         ( P1.Distance(PE2) > Precision::Confusion() ) )
00707       std::swap( LinEdge1, LinEdge2 );
00708 
00709     bool linEdge1Computed, linEdge2Computed;
00710     if ( !computeLayerPositions(P0,P1,LinEdge1,&linEdge1Computed))
00711       return false;
00712 
00713     Nodes1.resize( myLayerPositions.size()+1 );
00714     Nodes2.resize( myLayerPositions.size()+1 );
00715 
00716     // check that both linear edges have same hypotheses
00717     if ( !computeLayerPositions(P0,P2,LinEdge2, &linEdge2Computed))
00718          return false;
00719     if ( Nodes1.size() != myLayerPositions.size()+1 )
00720       return error("Different hypotheses apply to radial edges");
00721 
00722     // find the central vertex
00723     TopoDS_Vertex VC = V2;
00724     if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
00725         ( P2.Distance(PE1) > Precision::Confusion() ) )
00726       VC = V1;
00727     int vertID = meshDS->ShapeToIndex(VC);
00728 
00729     // LinEdge1
00730     if ( linEdge1Computed )
00731     {
00732       if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
00733         return error("Invalid mesh on a straight edge");
00734 
00735       bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
00736       NC = const_cast<SMDS_MeshNode*>
00737         ( nodesFromP0ToP1 ? theNodes.begin()->second : theNodes.rbegin()->second );
00738       int i = 0, ir = Nodes1.size()-1;
00739       int * pi = nodesFromP0ToP1 ? &i : &ir;
00740       itn = theNodes.begin();
00741       if ( nodesFromP0ToP1 ) ++itn;
00742       for ( ; i < Nodes1.size(); ++i, --ir, ++itn )
00743       {
00744         Nodes1[*pi] = itn->second;
00745       }
00746       for ( i = 0; i < Nodes1.size()-1; ++i )
00747       {
00748         Points.Append( gpXYZ( Nodes1[i]));
00749         Pnts2d1.Append( myHelper->GetNodeUV( F, Nodes1[i]));
00750       }
00751     }
00752     else
00753     {
00754       int edgeID = meshDS->ShapeToIndex(LinEdge1);
00755       gp_Vec aVec(P0,P1);
00756       // check orientation
00757       Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
00758       gp_Pnt Ptmp = Crv->Value(fp);
00759       bool ori = false;
00760       if( P1.Distance(Ptmp) > Precision::Confusion() )
00761         ori = true;
00762       // get UV points for edge
00763       gp_Pnt2d PF,PL;
00764       BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
00765       gp_Vec2d V2d;
00766       if(ori) {
00767         V2d = gp_Vec2d(PF,PL);
00768         PC = PF;
00769       }
00770       else {
00771         V2d = gp_Vec2d(PL,PF);
00772         PC = PL;
00773       }
00774       NC = const_cast<SMDS_MeshNode*>( VertexNode( VC, meshDS ));
00775       if ( !NC )
00776       {
00777         NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
00778         meshDS->SetNodeOnVertex(NC, vertID);
00779       }
00780       double dp = lp-fp;
00781       int i = 0;
00782       for(; i<myLayerPositions.size(); i++) {
00783         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
00784                   P0.Y() + aVec.Y()*myLayerPositions[i],
00785                   P0.Z() + aVec.Z()*myLayerPositions[i] );
00786         Points.Append(P);
00787         SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
00788         Nodes1[i] = node;
00789         double param;
00790         if(!ori)
00791           param = fp + dp*(1-myLayerPositions[i]);
00792         else
00793           param = fp + dp*myLayerPositions[i];
00794         meshDS->SetNodeOnEdge(node, edgeID, param);
00795         // parameters on face
00796         gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
00797                       PC.Y() + V2d.Y()*myLayerPositions[i] );
00798         Pnts2d1.Append(P2d);
00799       }
00800       Nodes1[ myLayerPositions.size() ] = NF;
00801       // create 1D elements on edge
00802       SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes1[0] );
00803       if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
00804       for(i=1; i<Nodes1.size(); i++) {
00805         ME = myHelper->AddEdge( Nodes1[i-1], Nodes1[i] );
00806         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
00807       }
00808       if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
00809         Nodes2 = Nodes1;
00810     }
00811     markEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
00812 
00813     // LinEdge2
00814     if ( linEdge2Computed )
00815     {
00816       if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge2,true,theNodes))
00817         return error("Invalid mesh on a straight edge");
00818 
00819       bool nodesFromP0ToP2 = ( theNodes.rbegin()->second == NL );
00820       int i = 0, ir = Nodes1.size()-1;
00821       int * pi = nodesFromP0ToP2 ? &i : &ir;
00822       itn = theNodes.begin();
00823       if ( nodesFromP0ToP2 ) ++itn;
00824       for ( ; i < Nodes2.size(); ++i, --ir, ++itn )
00825         Nodes2[*pi] = itn->second;
00826     }
00827     else
00828     {
00829       int edgeID = meshDS->ShapeToIndex(LinEdge2);
00830       gp_Vec aVec = gp_Vec(P0,P2);
00831       // check orientation
00832       Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge2,fp,lp);
00833       gp_Pnt Ptmp = Crv->Value(fp);
00834       bool ori = false;
00835       if( P2.Distance(Ptmp) > Precision::Confusion() )
00836         ori = true;
00837       // get UV points for edge
00838       gp_Pnt2d PF,PL;
00839       BRep_Tool::UVPoints( LinEdge2, TopoDS::Face(aShape), PF, PL );
00840       gp_Vec2d V2d;
00841       if(ori) {
00842         V2d = gp_Vec2d(PF,PL);
00843         PC = PF;
00844       }
00845       else {
00846         V2d = gp_Vec2d(PL,PF);
00847         PC = PL;
00848       }
00849       double dp = lp-fp;
00850       for(int i=0; i<myLayerPositions.size(); i++) {
00851         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
00852                   P0.Y() + aVec.Y()*myLayerPositions[i],
00853                   P0.Z() + aVec.Z()*myLayerPositions[i] );
00854         SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
00855         Nodes2[i] = node;
00856         double param;
00857         if(!ori)
00858           param = fp + dp*(1-myLayerPositions[i]);
00859         else
00860           param = fp + dp*myLayerPositions[i];
00861         meshDS->SetNodeOnEdge(node, edgeID, param);
00862         // parameters on face
00863         gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
00864                       PC.Y() + V2d.Y()*myLayerPositions[i] );
00865       }
00866       Nodes2[ myLayerPositions.size() ] = NL;
00867       // create 1D elements on edge
00868       SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes2[0] );
00869       if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
00870       for(int i=1; i<Nodes2.size(); i++) {
00871         ME = myHelper->AddEdge( Nodes2[i-1], Nodes2[i] );
00872         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
00873       }
00874     }
00875     markEdgeAsComputedByMe( LinEdge2, aMesh.GetSubMesh( F ));
00876   }
00877   markEdgeAsComputedByMe( CircEdge, aMesh.GetSubMesh( F ));
00878 
00879   // orientation
00880   bool IsForward = ( CircEdge.Orientation()==TopAbs_FORWARD );
00881 
00882   // create nodes and mesh elements on face
00883   // find axis of rotation
00884   gp_Pnt P2 = gp_Pnt( CNodes[1]->X(), CNodes[1]->Y(), CNodes[1]->Z() );
00885   gp_Vec Vec1(P0,P1);
00886   gp_Vec Vec2(P0,P2);
00887   gp_Vec Axis = Vec1.Crossed(Vec2);
00888   // create elements
00889   int i = 1;
00890   //cout<<"Angles.Length() = "<<Angles.Length()<<"   Points.Length() = "<<Points.Length()<<endl;
00891   //cout<<"Nodes1.size() = "<<Nodes1.size()<<"   Pnts2d1.Length() = "<<Pnts2d1.Length()<<endl;
00892   for(; i<Angles.Length(); i++) {
00893     vector< const SMDS_MeshNode* > tmpNodes;
00894     tmpNodes.reserve(Nodes1.size());
00895     gp_Trsf aTrsf;
00896     gp_Ax1 theAxis(P0,gp_Dir(Axis));
00897     aTrsf.SetRotation( theAxis, Angles.Value(i) );
00898     gp_Trsf2d aTrsf2d;
00899     aTrsf2d.SetRotation( PC, Angles.Value(i) );
00900     // create nodes
00901     int j = 1;
00902     for(; j<=Points.Length(); j++) {
00903       double cx,cy,cz;
00904       Points.Value(j).Coord( cx, cy, cz );
00905       aTrsf.Transforms( cx, cy, cz );
00906       SMDS_MeshNode* node = myHelper->AddNode( cx, cy, cz );
00907       // find parameters on face
00908       Pnts2d1.Value(j).Coord( cx, cy );
00909       aTrsf2d.Transforms( cx, cy );
00910       // set node on face
00911       meshDS->SetNodeOnFace( node, faceID, cx, cy );
00912       tmpNodes[j-1] = node;
00913     }
00914     // create faces
00915     tmpNodes[Points.Length()] = CNodes[i];
00916     // quad
00917     for(j=0; j<Nodes1.size()-1; j++) {
00918       SMDS_MeshFace* MF;
00919       if(IsForward)
00920         MF = myHelper->AddFace( tmpNodes[j], Nodes1[j],
00921                                 Nodes1[j+1], tmpNodes[j+1] );
00922       else
00923         MF = myHelper->AddFace( tmpNodes[j], tmpNodes[j+1],
00924                                 Nodes1[j+1], Nodes1[j] );
00925       if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
00926     }
00927     // tria
00928     SMDS_MeshFace* MF;
00929     if(IsForward)
00930       MF = myHelper->AddFace( NC, Nodes1[0], tmpNodes[0] );
00931     else
00932       MF = myHelper->AddFace( NC, tmpNodes[0], Nodes1[0] );
00933     if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
00934     for(j=0; j<Nodes1.size(); j++) {
00935       Nodes1[j] = tmpNodes[j];
00936     }
00937   }
00938   // create last faces
00939   // quad
00940   for(i=0; i<Nodes1.size()-1; i++) {
00941     SMDS_MeshFace* MF;
00942     if(IsForward)
00943       MF = myHelper->AddFace( Nodes2[i], Nodes1[i],
00944                               Nodes1[i+1], Nodes2[i+1] );
00945     else
00946       MF = myHelper->AddFace( Nodes2[i],  Nodes2[i+1],
00947                               Nodes1[i+1], Nodes1[i] );
00948     if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
00949   }
00950   // tria
00951   SMDS_MeshFace* MF;
00952   if(IsForward)
00953     MF = myHelper->AddFace( NC, Nodes1[0], Nodes2[0] );
00954   else
00955     MF = myHelper->AddFace( NC, Nodes2[0], Nodes1[0] );
00956   if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
00957 
00958   return true;
00959 }
00960 
00961 //================================================================================
00966 //================================================================================
00967 
00968 bool StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(const gp_Pnt&      p1,
00969                                                              const gp_Pnt&      p2,
00970                                                              const TopoDS_Edge& linEdge,
00971                                                              bool*              linEdgeComputed)
00972 {
00973   // First, try to compute positions of layers
00974 
00975   myLayerPositions.clear();
00976 
00977   SMESH_Mesh * mesh = myHelper->GetMesh();
00978 
00979   const SMESH_Hypothesis* hyp1D = myDistributionHypo ? myDistributionHypo->GetLayerDistribution() : 0;
00980   int                  nbLayers = myNbLayerHypo ? myNbLayerHypo->GetNumberOfLayers() : 0;
00981 
00982   if ( !hyp1D && !nbLayers )
00983   {
00984     // No own algo hypotheses assigned, so first try to find any 1D hypothesis.
00985     // We need some edge
00986     TopoDS_Shape edge = linEdge;
00987     if ( edge.IsNull() && !myHelper->GetSubShape().IsNull())
00988       for ( TopExp_Explorer e(myHelper->GetSubShape(), TopAbs_EDGE); e.More(); e.Next())
00989         edge = e.Current();
00990     if ( !edge.IsNull() )
00991     {
00992       // find a hyp usable by TNodeDistributor
00993       SMESH_HypoFilter hypKind;
00994       TNodeDistributor::GetDistributor(*mesh)->InitCompatibleHypoFilter(hypKind,/*ignoreAux=*/1);
00995       hyp1D = mesh->GetHypothesis( edge, hypKind, /*fromAncestors=*/true);
00996     }
00997   }
00998   if ( hyp1D ) // try to compute with hyp1D
00999   {
01000     if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions,p1,p2,*mesh,hyp1D )) {
01001       if ( myDistributionHypo ) { // bad hyp assigned 
01002         return error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
01003       }
01004       else {
01005         // bad hyp found, its Ok, lets try with default nb of segnents
01006       }
01007     }
01008   }
01009   
01010   if ( myLayerPositions.empty() ) // try to use nb of layers
01011   {
01012     if ( !nbLayers )
01013       nbLayers = _gen->GetDefaultNbSegments();
01014 
01015     if ( nbLayers )
01016     {
01017       myLayerPositions.resize( nbLayers - 1 );
01018       for ( int z = 1; z < nbLayers; ++z )
01019         myLayerPositions[ z - 1 ] = double( z )/ double( nbLayers );
01020     }
01021   }
01022 
01023   // Second, check presence of a mesh built by other algo on linEdge
01024   // and mesh conformity to my hypothesis
01025 
01026   bool meshComputed = (!linEdge.IsNull() && !mesh->GetSubMesh(linEdge)->IsEmpty() );
01027   if ( linEdgeComputed ) *linEdgeComputed = meshComputed;
01028 
01029   if ( meshComputed )
01030   {
01031     vector< double > nodeParams;
01032     GetNodeParamOnEdge( mesh->GetMeshDS(), linEdge, nodeParams );
01033 
01034     // nb of present nodes must be different in cases of 1 and 2 straight edges
01035 
01036     TopoDS_Vertex VV[2];
01037     TopExp::Vertices( linEdge, VV[0], VV[1]);
01038     const gp_Pnt* points[] = { &p1, &p2 };
01039     gp_Pnt       vPoints[] = { BRep_Tool::Pnt(VV[0]), BRep_Tool::Pnt(VV[1]) };
01040     const double     tol[] = { BRep_Tool::Tolerance(VV[0]), BRep_Tool::Tolerance(VV[1]) };
01041     bool pointsAreOnVertices = true;
01042     for ( int iP = 0; iP < 2 && pointsAreOnVertices; ++iP )
01043       pointsAreOnVertices = ( points[iP]->Distance( vPoints[0] ) < tol[0] ||
01044                               points[iP]->Distance( vPoints[1] ) < tol[1] );
01045 
01046     int nbNodes = nodeParams.size() - 2; // 2 straight edges
01047     if ( !pointsAreOnVertices )
01048       nbNodes = ( nodeParams.size() - 3 ) / 2; // 1 straight edge
01049 
01050     if ( myLayerPositions.empty() )
01051     {
01052       myLayerPositions.resize( nbNodes );
01053     }
01054     else if ( myDistributionHypo || myNbLayerHypo )
01055     {
01056       // linEdge is computed by other algo. Check if there is a meshed face
01057       // using nodes on linEdge
01058       bool nodesAreUsed = false;
01059       TopTools_ListIteratorOfListOfShape ancestIt = mesh->GetAncestors( linEdge );
01060       for ( ; ancestIt.More() && !nodesAreUsed; ancestIt.Next() )
01061         if ( ancestIt.Value().ShapeType() == TopAbs_FACE )
01062           nodesAreUsed = (!mesh->GetSubMesh( ancestIt.Value() )->IsEmpty());
01063       if ( !nodesAreUsed ) {
01064         // rebuild them
01065         mesh->GetSubMesh( linEdge )->ComputeStateEngine( SMESH_subMesh::CLEAN );
01066         if ( linEdgeComputed ) *linEdgeComputed = false;
01067       }
01068       else {
01069         
01070         if ( myLayerPositions.size() != nbNodes )
01071           return error("Radial edge is meshed by other algorithm");
01072       }
01073     }
01074   }
01075 
01076   return !myLayerPositions.empty();
01077 }
01078 
01079 
01080 //=======================================================================
01081 //function : Evaluate
01082 //purpose  : 
01083 //=======================================================================
01084 
01085 bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
01086                                                 const TopoDS_Shape& aShape,
01087                                                 MapShapeNbElems& aResMap)
01088 {
01089   if( aShape.ShapeType() != TopAbs_FACE ) {
01090     return false;
01091   }
01092   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
01093   if( aResMap.count(sm) )
01094     return false;
01095 
01096   vector<int>& aResVec =
01097     aResMap.insert( make_pair(sm, vector<int>(SMDSEntity_Last,0))).first->second;
01098 
01099   myHelper = new SMESH_MesherHelper( aMesh );
01100   myHelper->SetSubShape( aShape );
01101   auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
01102 
01103   TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
01104 
01105   TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
01106   int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 );
01107   if( nbe>3 || nbe < 1 || CircEdge.IsNull() )
01108     return false;
01109 
01110   Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
01111   if( aCirc.IsNull() )
01112     return error(COMPERR_BAD_SHAPE);
01113 
01114   gp_Pnt P0 = aCirc->Location();
01115   gp_Pnt P1 = aCirc->Value(0.);
01116   computeLayerPositions( P0, P1, LinEdge1 );
01117 
01118   int nb0d=0, nb2d_tria=0, nb2d_quad=0;
01119   bool isQuadratic = false, ok = true;
01120   if(nbe==1)
01121   {
01122     // C1 must be a circle
01123     ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap );
01124     if(ok) {
01125       const vector<int>& aVec = aResMap[aMesh.GetSubMesh(CircEdge)];
01126       isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
01127       if(isQuadratic) {
01128         // main nodes
01129         nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
01130         // radial medium nodes
01131         nb0d += (aVec[SMDSEntity_Node]+1) * (myLayerPositions.size()+1);
01132         // other medium nodes
01133         nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
01134       }
01135       else {
01136         nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
01137       }
01138       nb2d_tria = aVec[SMDSEntity_Node] + 1;
01139       nb2d_quad = nb0d;
01140     }
01141   }
01142   else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL)
01143   {
01144     // one curve must be a half of circle and other curve must be
01145     // a segment of line
01146     double fp, lp;
01147     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
01148     if( fabs(fabs(lp-fp)-M_PI) > Precision::Confusion() ) {
01149       // not half of circle
01150       return error(COMPERR_BAD_SHAPE);
01151     }
01152     Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
01153     if( aLine.IsNull() ) {
01154       // other curve not line
01155       return error(COMPERR_BAD_SHAPE);
01156     }
01157     ok = !aResMap.count( aMesh.GetSubMesh(LinEdge1) );
01158     if ( !ok ) {
01159       const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(LinEdge1) ];
01160       ok = ( aVec[SMDSEntity_Node] == myLayerPositions.size() );
01161     }
01162     if(ok) {
01163       ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap );
01164     }
01165     if(ok) {
01166       const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(CircEdge) ];
01167       isQuadratic = aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge];
01168       if(isQuadratic) {
01169         // main nodes
01170         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
01171         // radial medium nodes
01172         nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
01173         // other medium nodes
01174         nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
01175       }
01176       else {
01177         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
01178       }
01179       nb2d_tria = aVec[SMDSEntity_Node] + 1;
01180       nb2d_quad = nb2d_tria * myLayerPositions.size();
01181       // add evaluation for edges
01182       vector<int> aResVec(SMDSEntity_Last,0);
01183       if(isQuadratic) {
01184         aResVec[SMDSEntity_Node] = 4*myLayerPositions.size() + 3;
01185         aResVec[SMDSEntity_Quad_Edge] = 2*myLayerPositions.size() + 2;
01186       }
01187       else {
01188         aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
01189         aResVec[SMDSEntity_Edge] = 2*myLayerPositions.size() + 2;
01190       }
01191       aResMap[ aMesh.GetSubMesh(LinEdge1) ] = aResVec;
01192     }
01193   }
01194   else  // nbe==3 or ( nbe==2 && linEdge is INTERNAL )
01195   {
01196     if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
01197       LinEdge2 = LinEdge1;
01198 
01199     // one curve must be a part of circle and other curves must be
01200     // segments of line
01201     Handle(Geom_Line)  aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
01202     Handle(Geom_Line)  aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 ));
01203     if( aLine1.IsNull() || aLine2.IsNull() ) {
01204       // other curve not line
01205       return error(COMPERR_BAD_SHAPE);
01206     }
01207     int nbLayers = myLayerPositions.size();
01208     computeLayerPositions( P0, P1, LinEdge2 );
01209     if ( nbLayers != myLayerPositions.size() )
01210       return error("Different hypotheses apply to radial edges");
01211       
01212     bool ok = !aResMap.count( aMesh.GetSubMesh(LinEdge1));
01213     if ( !ok ) {
01214       if ( myDistributionHypo || myNbLayerHypo )
01215         ok = true; // override other 1d hyps
01216       else {
01217         const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(LinEdge1) ];
01218         ok = ( aVec[SMDSEntity_Node] == myLayerPositions.size() );
01219       }
01220     }
01221     if( ok && aResMap.count( aMesh.GetSubMesh(LinEdge2) )) {
01222       if ( myDistributionHypo || myNbLayerHypo )
01223         ok = true; // override other 1d hyps
01224       else {
01225         const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(LinEdge2) ];
01226         ok = ( aVec[SMDSEntity_Node] == myLayerPositions.size() );
01227       }
01228     }
01229     if(ok) {
01230       ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap );
01231     }
01232     if(ok) {
01233       const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(CircEdge) ];
01234       isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
01235       if(isQuadratic) {
01236         // main nodes
01237         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
01238         // radial medium nodes
01239         nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
01240         // other medium nodes
01241         nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
01242       }
01243       else {
01244         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
01245       }
01246       nb2d_tria = aVec[SMDSEntity_Node] + 1;
01247       nb2d_quad = nb2d_tria * myLayerPositions.size();
01248       // add evaluation for edges
01249       vector<int> aResVec(SMDSEntity_Last, 0);
01250       if(isQuadratic) {
01251         aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
01252         aResVec[SMDSEntity_Quad_Edge] = myLayerPositions.size() + 1;
01253       }
01254       else {
01255         aResVec[SMDSEntity_Node] = myLayerPositions.size();
01256         aResVec[SMDSEntity_Edge] = myLayerPositions.size() + 1;
01257       }
01258       sm = aMesh.GetSubMesh(LinEdge1);
01259       aResMap[sm] = aResVec;
01260       sm = aMesh.GetSubMesh(LinEdge2);
01261       aResMap[sm] = aResVec;
01262     }
01263   }
01264 
01265   if(nb0d>0) {
01266     aResVec[0] = nb0d;
01267     if(isQuadratic) {
01268       aResVec[SMDSEntity_Quad_Triangle] = nb2d_tria;
01269       aResVec[SMDSEntity_Quad_Quadrangle] = nb2d_quad;
01270     }
01271     else {
01272       aResVec[SMDSEntity_Triangle] = nb2d_tria;
01273       aResVec[SMDSEntity_Quadrangle] = nb2d_quad;
01274     }
01275     return true;
01276   }
01277 
01278   // invalid case
01279   sm = aMesh.GetSubMesh(aShape);
01280   SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
01281   smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
01282                                         "Submesh can not be evaluated",this));
01283   return false;
01284 
01285 }