1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 // SMESH SMESH : implementaion of SMESH idl descriptions
21 // File : StdMeshers_RadialQuadrangle_1D2D.cxx
24 #include "StdMeshers_RadialQuadrangle_1D2D.hxx"
26 #include "StdMeshers_NumberOfLayers.hxx"
27 #include "StdMeshers_LayerDistribution.hxx"
28 #include "StdMeshers_Regular_1D.hxx"
29 #include "StdMeshers_NumberOfSegments.hxx"
31 #include "SMDS_MeshNode.hxx"
32 #include "SMESHDS_Mesh.hxx"
33 #include "SMESHDS_SubMesh.hxx"
34 #include "SMESH_Gen.hxx"
35 #include "SMESH_HypoFilter.hxx"
36 #include "SMESH_Mesh.hxx"
37 #include "SMESH_MesherHelper.hxx"
38 #include "SMESH_subMesh.hxx"
39 #include "SMESH_subMeshEventListener.hxx"
41 #include "utilities.h"
43 #include <BRepAdaptor_Curve.hxx>
44 #include <BRepBuilderAPI_MakeEdge.hxx>
45 #include <BRep_Tool.hxx>
46 #include <GeomAPI_ProjectPointOnSurf.hxx>
47 #include <Geom_Circle.hxx>
48 #include <Geom_Line.hxx>
49 #include <Geom_TrimmedCurve.hxx>
50 #include <TColgp_SequenceOfPnt.hxx>
51 #include <TColgp_SequenceOfPnt2d.hxx>
53 #include <TopExp_Explorer.hxx>
54 #include <TopTools_ListIteratorOfListOfShape.hxx>
60 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
61 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
64 //=======================================================================
65 //function : StdMeshers_RadialQuadrangle_1D2D
67 //=======================================================================
69 StdMeshers_RadialQuadrangle_1D2D::StdMeshers_RadialQuadrangle_1D2D(int hypId,
72 :SMESH_2D_Algo(hypId, studyId, gen)
74 _name = "RadialQuadrangle_1D2D";
75 _shapeType = (1 << TopAbs_FACE); // 1 bit per shape type
77 _compatibleHypothesis.push_back("LayerDistribution2D");
78 _compatibleHypothesis.push_back("NumberOfLayers2D");
79 _requireDiscreteBoundary = false;
80 _supportSubmeshes = true;
81 _neededLowerHyps[ 1 ] = true; // suppress warning on hiding a global 1D algo
84 myDistributionHypo = 0;
88 //================================================================================
92 //================================================================================
94 StdMeshers_RadialQuadrangle_1D2D::~StdMeshers_RadialQuadrangle_1D2D()
98 //=======================================================================
99 //function : CheckHypothesis
101 //=======================================================================
103 bool StdMeshers_RadialQuadrangle_1D2D::CheckHypothesis
105 const TopoDS_Shape& aShape,
106 SMESH_Hypothesis::Hypothesis_Status& aStatus)
110 myDistributionHypo = 0;
112 list <const SMESHDS_Hypothesis * >::const_iterator itl;
114 const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
115 if ( hyps.size() == 0 ) {
116 aStatus = SMESH_Hypothesis::HYP_OK;
117 return true; // can work with no hypothesis
120 if ( hyps.size() > 1 ) {
121 aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
125 const SMESHDS_Hypothesis *theHyp = hyps.front();
127 string hypName = theHyp->GetName();
129 if (hypName == "NumberOfLayers2D") {
130 myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
131 aStatus = SMESH_Hypothesis::HYP_OK;
134 if (hypName == "LayerDistribution2D") {
135 myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
136 aStatus = SMESH_Hypothesis::HYP_OK;
139 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
145 // ------------------------------------------------------------------------------
147 * \brief Listener used to mark edges meshed by StdMeshers_RadialQuadrangle_1D2D
149 class TEdgeMarker : public SMESH_subMeshEventListener
151 TEdgeMarker(): SMESH_subMeshEventListener(/*isDeletable=*/false,
152 "StdMeshers_RadialQuadrangle_1D2D::TEdgeMarker") {}
154 //!< Return static listener
155 static SMESH_subMeshEventListener* getListener()
157 static TEdgeMarker theEdgeMarker;
158 return &theEdgeMarker;
160 //! Clear face sumbesh if something happens on edges
161 void ProcessEvent(const int event,
163 SMESH_subMesh* edgeSubMesh,
164 EventListenerData* data,
165 const SMESH_Hypothesis* /*hyp*/)
167 if ( data && !data->mySubMeshes.empty() && eventType == SMESH_subMesh::ALGO_EVENT)
169 ASSERT( data->mySubMeshes.front() != edgeSubMesh );
170 SMESH_subMesh* faceSubMesh = data->mySubMeshes.front();
171 faceSubMesh->ComputeStateEngine( SMESH_subMesh::CLEAN );
176 // ------------------------------------------------------------------------------
178 * \brief Mark an edge as computed by StdMeshers_RadialQuadrangle_1D2D
180 void markEdgeAsComputedByMe(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
182 if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
184 if ( !edgeSM->GetEventListenerData( TEdgeMarker::getListener() ))
185 faceSubMesh->SetEventListener( TEdgeMarker::getListener(),
186 SMESH_subMeshEventListenerData::MakeData(faceSubMesh),
190 // ------------------------------------------------------------------------------
192 * \brief Return true if a radial edge was meshed with StdMeshers_RadialQuadrangle_1D2D with
193 * the same radial distribution
195 // bool isEdgeCompatiballyMeshed(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
197 // if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
199 // if ( SMESH_subMeshEventListenerData* otherFaceData =
200 // edgeSM->GetEventListenerData( TEdgeMarker::getListener() ))
202 // // compare hypothesis aplied to two disk faces sharing radial edges
203 // SMESH_Mesh& mesh = *faceSubMesh->GetFather();
204 // SMESH_Algo* radialQuadAlgo = mesh.GetGen()->GetAlgo(mesh, faceSubMesh->GetSubShape() );
205 // SMESH_subMesh* otherFaceSubMesh = otherFaceData->mySubMeshes.front();
206 // list <const SMESHDS_Hypothesis *> hyps1 =
207 // radialQuadAlgo->GetUsedHypothesis( mesh, faceSubMesh->GetSubShape());
208 // list <const SMESHDS_Hypothesis *> hyps2 =
209 // radialQuadAlgo->GetUsedHypothesis( mesh, otherFaceSubMesh->GetSubShape());
210 // if( hyps1.empty() && hyps2.empty() )
211 // return true; // defaul hyps
212 // if ( hyps1.size() != hyps2.size() )
214 // return *hyps1.front() == *hyps2.front();
220 //================================================================================
222 * \brief Return base curve of the edge and extremum parameters
224 //================================================================================
226 Handle(Geom_Curve) getCurve(const TopoDS_Edge& edge, double* f=0, double* l=0)
228 Handle(Geom_Curve) C;
229 if ( !edge.IsNull() )
231 double first = 0., last = 0.;
232 C = BRep_Tool::Curve(edge, first, last);
235 Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C);
236 while( !tc.IsNull() ) {
237 C = tc->BasisCurve();
238 tc = Handle(Geom_TrimmedCurve)::DownCast(C);
247 //================================================================================
249 * \brief Return edges of the face
250 * \retval int - nb of edges
252 //================================================================================
254 int analyseFace(const TopoDS_Shape& face,
255 TopoDS_Edge& CircEdge,
256 TopoDS_Edge& LinEdge1,
257 TopoDS_Edge& LinEdge2)
259 CircEdge.Nullify(); LinEdge1.Nullify(); LinEdge2.Nullify();
262 for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next(), ++nbe )
264 const TopoDS_Edge& E = TopoDS::Edge( exp.Current() );
266 Handle(Geom_Curve) C = getCurve(E,&f,&l);
269 if ( C->IsKind( STANDARD_TYPE(Geom_Circle)))
271 if ( CircEdge.IsNull() )
276 else if ( LinEdge1.IsNull() )
284 //================================================================================
286 * \brief Checks if the common vertex between LinEdge's lies inside the circle
288 * \param [in] CircEdge -
289 * \param [in] LinEdge1 -
290 * \param [in] LinEdge2 -
291 * \return bool - false if there are 3 EDGEs and the corner is outside
293 //================================================================================
295 bool isCornerInsideCircle(const TopoDS_Edge& CircEdge,
296 const TopoDS_Edge& LinEdge1,
297 const TopoDS_Edge& LinEdge2)
299 if ( !CircEdge.IsNull() &&
300 !LinEdge1.IsNull() &&
303 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
304 TopoDS_Vertex aCommonV;
305 if ( !aCirc.IsNull() &&
306 TopExp::CommonVertex( LinEdge1, LinEdge2, aCommonV ))
308 gp_Pnt aCommonP = BRep_Tool::Pnt( aCommonV );
309 gp_Pnt aCenter = aCirc->Location();
310 double dist = aCenter.Distance( aCommonP );
311 return dist < 0.1 * aCirc->Radius();
317 //================================================================================
318 //================================================================================
320 * \brief Class computing layers distribution using data of
321 * StdMeshers_LayerDistribution hypothesis
323 //================================================================================
324 //================================================================================
326 class TNodeDistributor: public StdMeshers_Regular_1D
328 list <const SMESHDS_Hypothesis *> myUsedHyps;
330 // -----------------------------------------------------------------------------
331 static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
333 const int myID = -1001;
334 TNodeDistributor* myHyp = dynamic_cast<TNodeDistributor*>( aMesh.GetHypothesis( myID ));
336 myHyp = new TNodeDistributor( myID, 0, aMesh.GetGen() );
339 // -----------------------------------------------------------------------------
340 //! Computes distribution of nodes on a straight line ending at pIn and pOut
341 bool Compute( vector< double > & positions,
345 const SMESH_Hypothesis* hyp1d)
347 if ( !hyp1d ) return error( "Invalid LayerDistribution hypothesis");
349 double len = pIn.Distance( pOut );
350 if ( len <= DBL_MIN ) return error("Too close points of inner and outer shells");
353 myUsedHyps.push_back( hyp1d );
355 TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut );
356 SMESH_Hypothesis::Hypothesis_Status aStatus;
357 if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
358 return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
359 "with LayerDistribution hypothesis");
361 BRepAdaptor_Curve C3D(edge);
362 double f = C3D.FirstParameter(), l = C3D.LastParameter();
363 list< double > params;
364 if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false ))
365 return error("StdMeshers_Regular_1D failed to compute layers distribution");
368 positions.reserve( params.size() );
369 for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
370 positions.push_back( *itU / len );
373 // -----------------------------------------------------------------------------
374 //! Make mesh on an adge using assigned 1d hyp or defaut nb of segments
375 bool ComputeCircularEdge(SMESH_Mesh& aMesh,
376 const TopoDS_Edge& anEdge)
378 _gen->Compute( aMesh, anEdge);
379 SMESH_subMesh *sm = aMesh.GetSubMesh(anEdge);
380 if ( sm->GetComputeState() != SMESH_subMesh::COMPUTE_OK)
382 // find any 1d hyp assigned (there can be a hyp w/o algo)
383 myUsedHyps = SMESH_Algo::GetUsedHypothesis(aMesh, anEdge, /*ignoreAux=*/true);
384 Hypothesis_Status aStatus;
385 if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, anEdge, aStatus ))
387 // no valid 1d hyp assigned, use default nb of segments
388 _hypType = NB_SEGMENTS;
389 _ivalue[ DISTR_TYPE_IND ] = StdMeshers_NumberOfSegments::DT_Regular;
390 _ivalue[ NB_SEGMENTS_IND ] = _gen->GetDefaultNbSegments();
392 return StdMeshers_Regular_1D::Compute( aMesh, anEdge );
396 // -----------------------------------------------------------------------------
397 //! Make mesh on an adge using assigned 1d hyp or defaut nb of segments
398 bool EvaluateCircularEdge(SMESH_Mesh& aMesh,
399 const TopoDS_Edge& anEdge,
400 MapShapeNbElems& aResMap)
402 _gen->Evaluate( aMesh, anEdge, aResMap );
403 if ( aResMap.count( aMesh.GetSubMesh( anEdge )))
406 // find any 1d hyp assigned
407 myUsedHyps = SMESH_Algo::GetUsedHypothesis(aMesh, anEdge, /*ignoreAux=*/true);
408 Hypothesis_Status aStatus;
409 if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, anEdge, aStatus ))
411 // no valid 1d hyp assigned, use default nb of segments
412 _hypType = NB_SEGMENTS;
413 _ivalue[ DISTR_TYPE_IND ] = StdMeshers_NumberOfSegments::DT_Regular;
414 _ivalue[ NB_SEGMENTS_IND ] = _gen->GetDefaultNbSegments();
416 return StdMeshers_Regular_1D::Evaluate( aMesh, anEdge, aResMap );
419 // -----------------------------------------------------------------------------
420 TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen)
421 : StdMeshers_Regular_1D( hypId, studyId, gen)
424 // -----------------------------------------------------------------------------
425 virtual const list <const SMESHDS_Hypothesis *> &
426 GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
430 // -----------------------------------------------------------------------------
434 //=======================================================================
436 * \brief Allow algo to do something after persistent restoration
437 * \param subMesh - restored submesh
439 * call markEdgeAsComputedByMe()
441 //=======================================================================
443 void StdMeshers_RadialQuadrangle_1D2D::SubmeshRestored(SMESH_subMesh* faceSubMesh)
445 if ( !faceSubMesh->IsEmpty() )
447 TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
448 analyseFace( faceSubMesh->GetSubShape(), CircEdge, LinEdge1, LinEdge2 );
449 if ( !CircEdge.IsNull() ) markEdgeAsComputedByMe( CircEdge, faceSubMesh );
450 if ( !LinEdge1.IsNull() ) markEdgeAsComputedByMe( LinEdge1, faceSubMesh );
451 if ( !LinEdge2.IsNull() ) markEdgeAsComputedByMe( LinEdge2, faceSubMesh );
455 //=======================================================================
458 //=======================================================================
460 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
461 const TopoDS_Shape& aShape)
463 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
465 myHelper = new SMESH_MesherHelper( aMesh );
466 // to delete helper at exit from Compute()
467 SMESHUtils::Deleter<SMESH_MesherHelper> helperDeleter( myHelper );
469 TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
471 TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
472 int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 );
473 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
474 if( nbe > 3 || nbe < 1 || aCirc.IsNull() )
475 return error("The face must be a full circle or a part of circle (i.e. the number "
476 "of edges is less or equal to 3 and one of them is a circle curve)");
479 // points for rotation
480 TColgp_SequenceOfPnt Points;
481 // angles for rotation
482 TColStd_SequenceOfReal Angles;
483 // Nodes1 and Nodes2 - nodes along radiuses
484 // CNodes - nodes on circle edge
485 vector< const SMDS_MeshNode* > Nodes1, Nodes2, CNodes;
487 // parameters edge nodes on face
488 TColgp_SequenceOfPnt2d Pnts2d1;
491 int faceID = meshDS->ShapeToIndex(aShape);
492 TopoDS_Face F = TopoDS::Face(aShape);
493 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
498 if (!algo1d->ComputeCircularEdge( aMesh, CircEdge ))
499 return error( algo1d->GetComputeError() );
500 map< double, const SMDS_MeshNode* > theNodes;
501 if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes))
502 return error("Circular edge is incorrectly meshed");
504 myHelper->IsQuadraticSubMesh( aShape );
507 map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
508 const SMDS_MeshNode* NF = (*itn).second;
509 CNodes.push_back( (*itn).second );
510 double fang = (*itn).first;
511 if ( itn != theNodes.end() ) {
513 for(; itn != theNodes.end(); itn++ ) {
514 CNodes.push_back( (*itn).second );
515 double ang = (*itn).first - fang;
516 if( ang>M_PI ) ang = ang - 2.*M_PI;
517 if( ang<-M_PI ) ang = ang + 2.*M_PI;
518 Angles.Append( ang );
521 P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
522 P0 = aCirc->Location();
524 if ( !computeLayerPositions(P0,P1))
527 TopoDS_Vertex V1 = myHelper->IthVertex(0, CircEdge );
528 gp_Pnt2d p2dV = BRep_Tool::Parameters( V1, TopoDS::Face(aShape) );
530 NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
531 GeomAPI_ProjectPointOnSurf PPS(P0,S);
533 PPS.Parameters(1,U0,V0);
534 meshDS->SetNodeOnFace(NC, faceID, U0, V0);
535 PC = gp_Pnt2d(U0,V0);
538 gp_Vec2d aVec2d(PC,p2dV);
539 Nodes1.resize( myLayerPositions.size()+1 );
540 Nodes2.resize( myLayerPositions.size()+1 );
542 for ( ; i < myLayerPositions.size(); i++ ) {
543 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
544 P0.Y() + aVec.Y()*myLayerPositions[i],
545 P0.Z() + aVec.Z()*myLayerPositions[i] );
547 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
550 double U = PC.X() + aVec2d.X()*myLayerPositions[i];
551 double V = PC.Y() + aVec2d.Y()*myLayerPositions[i];
552 meshDS->SetNodeOnFace( node, faceID, U, V );
553 Pnts2d1.Append(gp_Pnt2d(U,V));
555 Nodes1[Nodes1.size()-1] = NF;
556 Nodes2[Nodes1.size()-1] = NF;
558 else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL )
560 // one curve must be a half of circle and other curve must be
563 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
564 if( fabs(fabs(lp-fp)-M_PI) > Precision::Confusion() ) {
565 // not half of circle
566 return error(COMPERR_BAD_SHAPE);
568 Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
569 if( aLine.IsNull() ) {
570 // other curve not line
571 return error(COMPERR_BAD_SHAPE);
574 if ( !algo1d->ComputeCircularEdge( aMesh, CircEdge ))
575 return error( algo1d->GetComputeError() );
576 map< double, const SMDS_MeshNode* > theNodes;
577 if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes) )
578 return error("Circular edge is incorrectly meshed");
580 myHelper->IsQuadraticSubMesh( aShape );
582 map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
584 CNodes.push_back( itn->second );
585 double fang = (*itn).first;
587 for(; itn != theNodes.end(); itn++ ) {
588 CNodes.push_back( (*itn).second );
589 double ang = (*itn).first - fang;
590 if( ang>M_PI ) ang = ang - 2.*M_PI;
591 if( ang<-M_PI ) ang = ang + 2.*M_PI;
592 Angles.Append( ang );
594 const SMDS_MeshNode* NF = theNodes.begin()->second;
595 const SMDS_MeshNode* NL = theNodes.rbegin()->second;
596 P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
597 gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
598 P0 = aCirc->Location();
600 bool linEdgeComputed;
601 if ( !computeLayerPositions(P0,P1,LinEdge1,&linEdgeComputed))
604 if ( linEdgeComputed )
606 if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
607 return error("Invalid mesh on a straight edge");
609 Nodes1.resize( myLayerPositions.size()+1 );
610 Nodes2.resize( myLayerPositions.size()+1 );
611 vector< const SMDS_MeshNode* > *pNodes1 = &Nodes1, *pNodes2 = &Nodes2;
612 bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
613 if ( !nodesFromP0ToP1 ) std::swap( pNodes1, pNodes2 );
615 map< double, const SMDS_MeshNode* >::reverse_iterator ritn = theNodes.rbegin();
616 itn = theNodes.begin();
617 for ( int i = Nodes1.size()-1; i > -1; ++itn, ++ritn, --i )
619 (*pNodes1)[i] = ritn->second;
620 (*pNodes2)[i] = itn->second;
621 Points.Prepend( gpXYZ( Nodes1[i]));
622 Pnts2d1.Prepend( myHelper->GetNodeUV( F, Nodes1[i]));
624 NC = const_cast<SMDS_MeshNode*>( itn->second );
625 Points.Remove( Nodes1.size() );
630 int edgeID = meshDS->ShapeToIndex(LinEdge1);
632 Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
636 if( P1.Distance(Ptmp) > Precision::Confusion() )
638 // get UV points for edge
640 BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
641 PC = gp_Pnt2d( (PF.X()+PL.X())/2, (PF.Y()+PL.Y())/2 );
643 if(ori) V2d = gp_Vec2d(PC,PF);
644 else V2d = gp_Vec2d(PC,PL);
646 double cp = (fp+lp)/2;
647 double dp2 = (lp-fp)/2;
648 NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
649 meshDS->SetNodeOnEdge(NC, edgeID, cp);
650 Nodes1.resize( myLayerPositions.size()+1 );
651 Nodes2.resize( myLayerPositions.size()+1 );
653 for(; i<myLayerPositions.size(); i++) {
654 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
655 P0.Y() + aVec.Y()*myLayerPositions[i],
656 P0.Z() + aVec.Z()*myLayerPositions[i] );
658 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
662 param = fp + dp2*(1-myLayerPositions[i]);
664 param = cp + dp2*myLayerPositions[i];
665 meshDS->SetNodeOnEdge(node, edgeID, param);
666 P = gp_Pnt( P0.X() - aVec.X()*myLayerPositions[i],
667 P0.Y() - aVec.Y()*myLayerPositions[i],
668 P0.Z() - aVec.Z()*myLayerPositions[i] );
669 node = meshDS->AddNode(P.X(), P.Y(), P.Z());
672 param = fp + dp2*(1-myLayerPositions[i]);
674 param = cp + dp2*myLayerPositions[i];
675 meshDS->SetNodeOnEdge(node, edgeID, param);
676 // parameters on face
677 gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
678 PC.Y() + V2d.Y()*myLayerPositions[i] );
681 Nodes1[ myLayerPositions.size() ] = NF;
682 Nodes2[ myLayerPositions.size() ] = NL;
683 // create 1D elements on edge
684 vector< const SMDS_MeshNode* > tmpNodes;
685 tmpNodes.resize(2*Nodes1.size()+1);
686 for(i=0; i<Nodes2.size(); i++)
687 tmpNodes[Nodes2.size()-i-1] = Nodes2[i];
688 tmpNodes[Nodes2.size()] = NC;
689 for(i=0; i<Nodes1.size(); i++)
690 tmpNodes[Nodes2.size()+1+i] = Nodes1[i];
691 for(i=1; i<tmpNodes.size(); i++) {
692 SMDS_MeshEdge* ME = myHelper->AddEdge( tmpNodes[i-1], tmpNodes[i] );
693 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
695 markEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
698 else // nbe==3 or ( nbe==2 && linEdge is INTERNAL )
700 if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
703 // one curve must be a part of circle and other curves must be
706 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
707 Handle(Geom_Line) aLine1 = Handle(Geom_Line )::DownCast( getCurve( LinEdge1 ));
708 Handle(Geom_Line) aLine2 = Handle(Geom_Line )::DownCast( getCurve( LinEdge2 ));
709 if ( aCirc.IsNull() || aLine1.IsNull() || aLine2.IsNull() )
710 return error(COMPERR_BAD_SHAPE);
711 if ( !isCornerInsideCircle( CircEdge, LinEdge1, LinEdge2 ))
712 return error(COMPERR_BAD_SHAPE);
714 if ( !algo1d->ComputeCircularEdge( aMesh, CircEdge ))
715 return error( algo1d->GetComputeError() );
716 map< double, const SMDS_MeshNode* > theNodes;
717 if ( !GetSortedNodesOnEdge( aMesh.GetMeshDS(), CircEdge, true, theNodes ))
718 return error("Circular edge is incorrectly meshed");
720 myHelper->IsQuadraticSubMesh( aShape );
722 const SMDS_MeshNode* NF = theNodes.begin()->second;
723 const SMDS_MeshNode* NL = theNodes.rbegin()->second;
725 CNodes.push_back( NF );
726 map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
727 double fang = (*itn).first;
729 for(; itn != theNodes.end(); itn++ ) {
730 CNodes.push_back( (*itn).second );
731 double ang = (*itn).first - fang;
732 if( ang>M_PI ) ang = ang - 2.*M_PI;
733 if( ang<-M_PI ) ang = ang + 2.*M_PI;
734 Angles.Append( ang );
736 P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
737 gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
738 P0 = aCirc->Location();
740 // make P1 belong to LinEdge1
741 TopoDS_Vertex V1 = myHelper->IthVertex( 0, LinEdge1 );
742 TopoDS_Vertex V2 = myHelper->IthVertex( 1, LinEdge1 );
743 gp_Pnt PE1 = BRep_Tool::Pnt(V1);
744 gp_Pnt PE2 = BRep_Tool::Pnt(V2);
745 if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
746 ( P1.Distance(PE2) > Precision::Confusion() ) )
747 std::swap( LinEdge1, LinEdge2 );
749 bool linEdge1Computed, linEdge2Computed;
750 if ( !computeLayerPositions(P0,P1,LinEdge1,&linEdge1Computed))
753 Nodes1.resize( myLayerPositions.size()+1 );
754 Nodes2.resize( myLayerPositions.size()+1 );
756 // check that both linear edges have same hypotheses
757 if ( !computeLayerPositions(P0,P2,LinEdge2, &linEdge2Computed))
759 if ( Nodes1.size() != myLayerPositions.size()+1 )
760 return error("Different hypotheses apply to radial edges");
762 // find the central vertex
763 TopoDS_Vertex VC = V2;
764 if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
765 ( P2.Distance(PE1) > Precision::Confusion() ) )
767 int vertID = meshDS->ShapeToIndex(VC);
770 if ( linEdge1Computed )
772 if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
773 return error("Invalid mesh on a straight edge");
775 bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
776 NC = const_cast<SMDS_MeshNode*>
777 ( nodesFromP0ToP1 ? theNodes.begin()->second : theNodes.rbegin()->second );
778 size_t i = 0, ir = Nodes1.size()-1;
779 size_t * pi = nodesFromP0ToP1 ? &i : &ir;
780 itn = theNodes.begin();
781 if ( nodesFromP0ToP1 ) ++itn;
782 for ( ; i < Nodes1.size(); ++i, --ir, ++itn )
784 Nodes1[*pi] = itn->second;
786 for ( i = 0; i < Nodes1.size()-1; ++i )
788 Points.Append( gpXYZ( Nodes1[i]));
789 Pnts2d1.Append( myHelper->GetNodeUV( F, Nodes1[i]));
794 int edgeID = meshDS->ShapeToIndex(LinEdge1);
797 Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
798 gp_Pnt Ptmp = Crv->Value(fp);
800 if( P1.Distance(Ptmp) > Precision::Confusion() )
802 // get UV points for edge
804 BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
807 V2d = gp_Vec2d(PF,PL);
811 V2d = gp_Vec2d(PL,PF);
814 NC = const_cast<SMDS_MeshNode*>( VertexNode( VC, meshDS ));
817 NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
818 meshDS->SetNodeOnVertex(NC, vertID);
822 for ( ; i < myLayerPositions.size(); i++ ) {
823 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
824 P0.Y() + aVec.Y()*myLayerPositions[i],
825 P0.Z() + aVec.Z()*myLayerPositions[i] );
827 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
831 param = fp + dp*(1-myLayerPositions[i]);
833 param = fp + dp*myLayerPositions[i];
834 meshDS->SetNodeOnEdge(node, edgeID, param);
835 // parameters on face
836 gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
837 PC.Y() + V2d.Y()*myLayerPositions[i] );
840 Nodes1[ myLayerPositions.size() ] = NF;
841 // create 1D elements on edge
842 SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes1[0] );
843 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
844 for ( i = 1; i < Nodes1.size(); i++ ) {
845 ME = myHelper->AddEdge( Nodes1[i-1], Nodes1[i] );
846 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
848 if ( nbe == 2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
851 markEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
854 if ( linEdge2Computed )
856 if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge2,true,theNodes))
857 return error("Invalid mesh on a straight edge");
859 bool nodesFromP0ToP2 = ( theNodes.rbegin()->second == NL );
860 size_t i = 0, ir = Nodes1.size()-1;
861 size_t * pi = nodesFromP0ToP2 ? &i : &ir;
862 itn = theNodes.begin();
863 if ( nodesFromP0ToP2 ) ++itn;
864 for ( ; i < Nodes2.size(); ++i, --ir, ++itn )
865 Nodes2[*pi] = itn->second;
869 int edgeID = meshDS->ShapeToIndex(LinEdge2);
870 gp_Vec aVec = gp_Vec(P0,P2);
872 Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge2,fp,lp);
873 gp_Pnt Ptmp = Crv->Value(fp);
875 if( P2.Distance(Ptmp) > Precision::Confusion() )
877 // get UV points for edge
879 BRep_Tool::UVPoints( LinEdge2, TopoDS::Face(aShape), PF, PL );
882 V2d = gp_Vec2d(PF,PL);
886 V2d = gp_Vec2d(PL,PF);
890 for ( size_t i = 0; i < myLayerPositions.size(); i++ ) {
891 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
892 P0.Y() + aVec.Y()*myLayerPositions[i],
893 P0.Z() + aVec.Z()*myLayerPositions[i] );
894 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
898 param = fp + dp*(1-myLayerPositions[i]);
900 param = fp + dp*myLayerPositions[i];
901 meshDS->SetNodeOnEdge(node, edgeID, param);
902 // parameters on face
903 gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
904 PC.Y() + V2d.Y()*myLayerPositions[i] );
906 Nodes2[ myLayerPositions.size() ] = NL;
907 // create 1D elements on edge
908 SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes2[0] );
909 if ( ME ) meshDS->SetMeshElementOnShape(ME, edgeID);
910 for ( size_t i = 1; i < Nodes2.size(); i++ ) {
911 ME = myHelper->AddEdge( Nodes2[i-1], Nodes2[i] );
912 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
915 markEdgeAsComputedByMe( LinEdge2, aMesh.GetSubMesh( F ));
917 markEdgeAsComputedByMe( CircEdge, aMesh.GetSubMesh( F ));
920 bool IsForward = ( CircEdge.Orientation()==TopAbs_FORWARD );
921 const double angleSign = ( F.Orientation() == TopAbs_REVERSED ? -1.0 : 1.0 );
923 // create nodes and mesh elements on face
924 // find axis of rotation
925 gp_Pnt P2 = gp_Pnt( CNodes[1]->X(), CNodes[1]->Y(), CNodes[1]->Z() );
928 gp_Vec Axis = Vec1.Crossed(Vec2);
931 //cout<<"Angles.Length() = "<<Angles.Length()<<" Points.Length() = "<<Points.Length()<<endl;
932 //cout<<"Nodes1.size() = "<<Nodes1.size()<<" Pnts2d1.Length() = "<<Pnts2d1.Length()<<endl;
933 for(; i<Angles.Length(); i++) {
934 vector< const SMDS_MeshNode* > tmpNodes;
936 gp_Ax1 theAxis(P0,gp_Dir(Axis));
937 aTrsf.SetRotation( theAxis, Angles.Value(i) );
939 aTrsf2d.SetRotation( PC, Angles.Value(i) * angleSign );
942 for(; j<=Points.Length(); j++) {
944 Points.Value(j).Coord( cx, cy, cz );
945 aTrsf.Transforms( cx, cy, cz );
946 SMDS_MeshNode* node = myHelper->AddNode( cx, cy, cz );
947 // find parameters on face
948 Pnts2d1.Value(j).Coord( cx, cy );
949 aTrsf2d.Transforms( cx, cy );
951 meshDS->SetNodeOnFace( node, faceID, cx, cy );
952 tmpNodes.push_back(node);
955 tmpNodes.push_back( CNodes[i] );
957 for ( j = 0; j < (int)Nodes1.size() - 1; j++ ) {
960 MF = myHelper->AddFace( tmpNodes[j], Nodes1[j],
961 Nodes1[j+1], tmpNodes[j+1] );
963 MF = myHelper->AddFace( tmpNodes[j], tmpNodes[j+1],
964 Nodes1[j+1], Nodes1[j] );
965 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
970 MF = myHelper->AddFace( NC, Nodes1[0], tmpNodes[0] );
972 MF = myHelper->AddFace( NC, tmpNodes[0], Nodes1[0] );
973 if ( MF ) meshDS->SetMeshElementOnShape(MF, faceID);
974 for ( j = 0; j < (int) Nodes1.size(); j++ ) {
975 Nodes1[j] = tmpNodes[j];
980 for ( i = 0; i < (int)Nodes1.size()-1; i++ ) {
983 MF = myHelper->AddFace( Nodes2[i], Nodes1[i],
984 Nodes1[i+1], Nodes2[i+1] );
986 MF = myHelper->AddFace( Nodes2[i], Nodes2[i+1],
987 Nodes1[i+1], Nodes1[i] );
988 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
993 MF = myHelper->AddFace( NC, Nodes1[0], Nodes2[0] );
995 MF = myHelper->AddFace( NC, Nodes2[0], Nodes1[0] );
996 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
1001 //================================================================================
1003 * \brief Compute positions of nodes on the radial edge
1004 * \retval bool - is a success
1006 //================================================================================
1008 bool StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(const gp_Pnt& p1,
1010 const TopoDS_Edge& linEdge,
1011 bool* linEdgeComputed)
1013 // First, try to compute positions of layers
1015 myLayerPositions.clear();
1017 SMESH_Mesh * mesh = myHelper->GetMesh();
1019 const SMESH_Hypothesis* hyp1D = myDistributionHypo ? myDistributionHypo->GetLayerDistribution() : 0;
1020 int nbLayers = myNbLayerHypo ? myNbLayerHypo->GetNumberOfLayers() : 0;
1022 if ( !hyp1D && !nbLayers )
1024 // No own algo hypotheses assigned, so first try to find any 1D hypothesis.
1025 // We need some edge
1026 TopoDS_Shape edge = linEdge;
1027 if ( edge.IsNull() && !myHelper->GetSubShape().IsNull())
1028 for ( TopExp_Explorer e(myHelper->GetSubShape(), TopAbs_EDGE); e.More(); e.Next())
1030 if ( !edge.IsNull() )
1032 // find a hyp usable by TNodeDistributor
1033 const SMESH_HypoFilter* hypKind =
1034 TNodeDistributor::GetDistributor(*mesh)->GetCompatibleHypoFilter(/*ignoreAux=*/true);
1035 hyp1D = mesh->GetHypothesis( edge, *hypKind, /*fromAncestors=*/true);
1038 if ( hyp1D ) // try to compute with hyp1D
1040 if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions,p1,p2,*mesh,hyp1D )) {
1041 if ( myDistributionHypo ) { // bad hyp assigned
1042 return error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
1045 // bad hyp found, its Ok, lets try with default nb of segnents
1050 if ( myLayerPositions.empty() ) // try to use nb of layers
1053 nbLayers = _gen->GetDefaultNbSegments();
1057 myLayerPositions.resize( nbLayers - 1 );
1058 for ( int z = 1; z < nbLayers; ++z )
1059 myLayerPositions[ z - 1 ] = double( z )/ double( nbLayers );
1063 // Second, check presence of a mesh built by other algo on linEdge
1064 // and mesh conformity to my hypothesis
1066 bool meshComputed = (!linEdge.IsNull() && !mesh->GetSubMesh(linEdge)->IsEmpty() );
1067 if ( linEdgeComputed ) *linEdgeComputed = meshComputed;
1071 vector< double > nodeParams;
1072 GetNodeParamOnEdge( mesh->GetMeshDS(), linEdge, nodeParams );
1074 // nb of present nodes must be different in cases of 1 and 2 straight edges
1076 TopoDS_Vertex VV[2];
1077 TopExp::Vertices( linEdge, VV[0], VV[1]);
1078 const gp_Pnt* points[] = { &p1, &p2 };
1079 gp_Pnt vPoints[] = { BRep_Tool::Pnt(VV[0]), BRep_Tool::Pnt(VV[1]) };
1080 const double tol[] = { BRep_Tool::Tolerance(VV[0]), BRep_Tool::Tolerance(VV[1]) };
1081 bool pointsAreOnVertices = true;
1082 for ( int iP = 0; iP < 2 && pointsAreOnVertices; ++iP )
1083 pointsAreOnVertices = ( points[iP]->Distance( vPoints[0] ) < tol[0] ||
1084 points[iP]->Distance( vPoints[1] ) < tol[1] );
1086 int nbNodes = nodeParams.size() - 2; // 2 straight edges
1087 if ( !pointsAreOnVertices )
1088 nbNodes = ( nodeParams.size() - 3 ) / 2; // 1 straight edge
1090 if ( myLayerPositions.empty() )
1092 myLayerPositions.resize( nbNodes );
1094 else if ( myDistributionHypo || myNbLayerHypo )
1096 // linEdge is computed by other algo. Check if there is a meshed face
1097 // using nodes on linEdge
1098 bool nodesAreUsed = false;
1099 TopTools_ListIteratorOfListOfShape ancestIt = mesh->GetAncestors( linEdge );
1100 for ( ; ancestIt.More() && !nodesAreUsed; ancestIt.Next() )
1101 if ( ancestIt.Value().ShapeType() == TopAbs_FACE )
1102 nodesAreUsed = (!mesh->GetSubMesh( ancestIt.Value() )->IsEmpty());
1103 if ( !nodesAreUsed ) {
1105 mesh->GetSubMesh( linEdge )->ComputeStateEngine( SMESH_subMesh::CLEAN );
1106 if ( linEdgeComputed ) *linEdgeComputed = false;
1110 if ((int) myLayerPositions.size() != nbNodes )
1111 return error("Radial edge is meshed by other algorithm");
1116 return !myLayerPositions.empty();
1120 //=======================================================================
1121 //function : Evaluate
1123 //=======================================================================
1125 bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
1126 const TopoDS_Shape& aShape,
1127 MapShapeNbElems& aResMap)
1129 if( aShape.ShapeType() != TopAbs_FACE ) {
1132 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1133 if( aResMap.count(sm) )
1136 vector<int>& aResVec =
1137 aResMap.insert( make_pair(sm, vector<int>(SMDSEntity_Last,0))).first->second;
1139 myHelper = new SMESH_MesherHelper( aMesh );
1140 myHelper->SetSubShape( aShape );
1141 auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
1143 TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
1145 TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
1146 int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 );
1147 if( nbe>3 || nbe < 1 || CircEdge.IsNull() )
1150 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
1151 if( aCirc.IsNull() )
1152 return error(COMPERR_BAD_SHAPE);
1154 gp_Pnt P0 = aCirc->Location();
1155 gp_Pnt P1 = aCirc->Value(0.);
1156 computeLayerPositions( P0, P1, LinEdge1 );
1158 int nb0d=0, nb2d_tria=0, nb2d_quad=0;
1159 bool isQuadratic = false, ok = true;
1162 // C1 must be a circle
1163 ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap );
1165 const vector<int>& aVec = aResMap[aMesh.GetSubMesh(CircEdge)];
1166 isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1169 nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1170 // radial medium nodes
1171 nb0d += (aVec[SMDSEntity_Node]+1) * (myLayerPositions.size()+1);
1172 // other medium nodes
1173 nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1176 nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1178 nb2d_tria = aVec[SMDSEntity_Node] + 1;
1182 else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL)
1184 // one curve must be a half of circle and other curve must be
1185 // a segment of line
1187 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
1188 if( fabs(fabs(lp-fp)-M_PI) > Precision::Confusion() ) {
1189 // not half of circle
1190 return error(COMPERR_BAD_SHAPE);
1192 Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
1193 if( aLine.IsNull() ) {
1194 // other curve not line
1195 return error(COMPERR_BAD_SHAPE);
1197 ok = !aResMap.count( aMesh.GetSubMesh(LinEdge1) );
1199 const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(LinEdge1) ];
1200 ok = ( aVec[SMDSEntity_Node] == (int) myLayerPositions.size() );
1203 ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap );
1206 const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(CircEdge) ];
1207 isQuadratic = aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge];
1210 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1211 // radial medium nodes
1212 nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1213 // other medium nodes
1214 nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1217 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1219 nb2d_tria = aVec[SMDSEntity_Node] + 1;
1220 nb2d_quad = nb2d_tria * myLayerPositions.size();
1221 // add evaluation for edges
1222 vector<int> aResVec(SMDSEntity_Last,0);
1224 aResVec[SMDSEntity_Node] = 4*myLayerPositions.size() + 3;
1225 aResVec[SMDSEntity_Quad_Edge] = 2*myLayerPositions.size() + 2;
1228 aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1229 aResVec[SMDSEntity_Edge] = 2*myLayerPositions.size() + 2;
1231 aResMap[ aMesh.GetSubMesh(LinEdge1) ] = aResVec;
1234 else // nbe==3 or ( nbe==2 && linEdge is INTERNAL )
1236 if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
1237 LinEdge2 = LinEdge1;
1239 // one curve must be a part of circle and other curves must be
1241 Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
1242 Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 ));
1243 if( aLine1.IsNull() || aLine2.IsNull() ) {
1244 // other curve not line
1245 return error(COMPERR_BAD_SHAPE);
1247 size_t nbLayers = myLayerPositions.size();
1248 computeLayerPositions( P0, P1, LinEdge2 );
1249 if ( nbLayers != myLayerPositions.size() )
1250 return error("Different hypotheses apply to radial edges");
1252 bool ok = !aResMap.count( aMesh.GetSubMesh(LinEdge1));
1254 if ( myDistributionHypo || myNbLayerHypo )
1255 ok = true; // override other 1d hyps
1257 const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(LinEdge1) ];
1258 ok = ( aVec[SMDSEntity_Node] == (int) myLayerPositions.size() );
1261 if( ok && aResMap.count( aMesh.GetSubMesh(LinEdge2) )) {
1262 if ( myDistributionHypo || myNbLayerHypo )
1263 ok = true; // override other 1d hyps
1265 const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(LinEdge2) ];
1266 ok = ( aVec[SMDSEntity_Node] == (int) myLayerPositions.size() );
1270 ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap );
1273 const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(CircEdge) ];
1274 isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1277 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1278 // radial medium nodes
1279 nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1280 // other medium nodes
1281 nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1284 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1286 nb2d_tria = aVec[SMDSEntity_Node] + 1;
1287 nb2d_quad = nb2d_tria * myLayerPositions.size();
1288 // add evaluation for edges
1289 vector<int> aResVec(SMDSEntity_Last, 0);
1291 aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1292 aResVec[SMDSEntity_Quad_Edge] = myLayerPositions.size() + 1;
1295 aResVec[SMDSEntity_Node] = myLayerPositions.size();
1296 aResVec[SMDSEntity_Edge] = myLayerPositions.size() + 1;
1298 sm = aMesh.GetSubMesh(LinEdge1);
1299 aResMap[sm] = aResVec;
1300 sm = aMesh.GetSubMesh(LinEdge2);
1301 aResMap[sm] = aResVec;
1308 aResVec[SMDSEntity_Quad_Triangle] = nb2d_tria;
1309 aResVec[SMDSEntity_Quad_Quadrangle] = nb2d_quad;
1312 aResVec[SMDSEntity_Triangle] = nb2d_tria;
1313 aResVec[SMDSEntity_Quadrangle] = nb2d_quad;
1319 sm = aMesh.GetSubMesh(aShape);
1320 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1321 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
1322 "Submesh can not be evaluated",this));
1327 //================================================================================
1329 * \brief Return true if the algorithm can compute mesh on this shape
1331 //================================================================================
1333 bool StdMeshers_RadialQuadrangle_1D2D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
1335 int nbFoundFaces = 0;
1336 for (TopExp_Explorer exp( aShape, TopAbs_FACE ); exp.More(); exp.Next(), ++nbFoundFaces )
1338 TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
1339 int nbe = analyseFace( exp.Current(), CircEdge, LinEdge1, LinEdge2 );
1340 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
1341 bool ok = ( nbe <= 3 && nbe >= 1 && !aCirc.IsNull() &&
1342 isCornerInsideCircle( CircEdge, LinEdge1, LinEdge2 ));
1343 if( toCheckAll && !ok ) return false;
1344 if( !toCheckAll && ok ) return true;
1346 if( toCheckAll && nbFoundFaces != 0 ) return true;