1 // Copyright (C) 2007-2019 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 // File : StdMeshers_RadialQuadrangle_1D2D.cxx
23 #include "StdMeshers_RadialQuadrangle_1D2D.hxx"
25 #include "StdMeshers_NumberOfLayers.hxx"
26 #include "StdMeshers_LayerDistribution.hxx"
27 #include "StdMeshers_Regular_1D.hxx"
28 #include "StdMeshers_NumberOfSegments.hxx"
29 #include "StdMeshers_FaceSide.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"
40 #include "SMESH_Block.hxx"
42 #include "utilities.h"
44 #include <BRepAdaptor_CompCurve.hxx>
45 #include <BRepAdaptor_Curve.hxx>
46 #include <BRepBuilderAPI_MakeEdge.hxx>
47 #include <BRep_Builder.hxx>
48 #include <BRep_Tool.hxx>
49 #include <GCPnts_AbscissaPoint.hxx>
50 #include <Geom2d_Line.hxx>
51 #include <Geom2d_TrimmedCurve.hxx>
52 #include <GeomAPI_ProjectPointOnSurf.hxx>
53 #include <Geom_Circle.hxx>
54 #include <Geom_Line.hxx>
55 #include <Geom_TrimmedCurve.hxx>
56 #include <ShapeFix_Edge.hxx>
57 #include <TColgp_SequenceOfPnt.hxx>
58 #include <TColgp_SequenceOfPnt2d.hxx>
60 #include <TopExp_Explorer.hxx>
61 #include <TopTools_ListIteratorOfListOfShape.hxx>
67 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
68 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
71 //=======================================================================
72 //function : StdMeshers_RadialQuadrangle_1D2D
74 //=======================================================================
76 StdMeshers_RadialQuadrangle_1D2D::StdMeshers_RadialQuadrangle_1D2D(int hypId,
78 :StdMeshers_Quadrangle_2D( hypId, gen )
80 _name = "RadialQuadrangle_1D2D";
81 _shapeType = (1 << TopAbs_FACE); // 1 bit per shape type
83 _compatibleHypothesis.push_back("LayerDistribution2D");
84 _compatibleHypothesis.push_back("NumberOfLayers2D");
85 _requireDiscreteBoundary = false;
86 _supportSubmeshes = true;
87 _neededLowerHyps[ 1 ] = true; // suppress warning on hiding a global 1D algo
90 myDistributionHypo = 0;
94 //================================================================================
98 //================================================================================
100 StdMeshers_RadialQuadrangle_1D2D::~StdMeshers_RadialQuadrangle_1D2D()
104 //=======================================================================
105 //function : CheckHypothesis
107 //=======================================================================
109 bool StdMeshers_RadialQuadrangle_1D2D::CheckHypothesis
111 const TopoDS_Shape& aShape,
112 SMESH_Hypothesis::Hypothesis_Status& aStatus)
116 myDistributionHypo = 0;
118 list <const SMESHDS_Hypothesis * >::const_iterator itl;
120 const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
121 if ( hyps.size() == 0 ) {
122 aStatus = SMESH_Hypothesis::HYP_OK;
123 return true; // can work with no hypothesis
126 if ( hyps.size() > 1 ) {
127 aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
131 const SMESHDS_Hypothesis *theHyp = hyps.front();
133 string hypName = theHyp->GetName();
135 if (hypName == "NumberOfLayers2D") {
136 myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
137 aStatus = SMESH_Hypothesis::HYP_OK;
140 if (hypName == "LayerDistribution2D") {
141 myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
142 aStatus = SMESH_Hypothesis::HYP_OK;
145 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
151 // ------------------------------------------------------------------------------
153 * \brief Listener used to mark edges meshed by StdMeshers_RadialQuadrangle_1D2D
155 class TEdgeMarker : public SMESH_subMeshEventListener
157 TEdgeMarker(): SMESH_subMeshEventListener(/*isDeletable=*/false,
158 "StdMeshers_RadialQuadrangle_1D2D::TEdgeMarker") {}
160 //!< Return static listener
161 static SMESH_subMeshEventListener* getListener()
163 static TEdgeMarker theEdgeMarker;
164 return &theEdgeMarker;
166 //! Clear edge sumbesh if something happens on face
167 void ProcessEvent(const int event,
169 SMESH_subMesh* faceSubMesh,
170 EventListenerData* edgesHolder,
171 const SMESH_Hypothesis* /*hyp*/)
173 if ( edgesHolder && eventType == SMESH_subMesh::ALGO_EVENT)
175 std::list<SMESH_subMesh*>::iterator smIt = edgesHolder->mySubMeshes.begin();
176 for ( ; smIt != edgesHolder->mySubMeshes.end(); ++smIt )
178 SMESH_subMesh* edgeSM = *smIt;
179 edgeSM->ComputeStateEngine( SMESH_subMesh::CLEAN );
183 //! Store edge SMESH_subMesh'es computed by the algo
184 static void markEdge( const TopoDS_Edge& edge, SMESH_subMesh* faceSM )
186 if ( SMESH_subMesh* edgeSM = faceSM->GetFather()->GetSubMeshContaining( edge ))
188 EventListenerData* edgesHolder = faceSM->GetEventListenerData( getListener() );
191 std::list<SMESH_subMesh*>::iterator smIt = std::find( edgesHolder->mySubMeshes.begin(),
192 edgesHolder->mySubMeshes.end(),
194 if ( smIt == edgesHolder->mySubMeshes.end() )
195 edgesHolder->mySubMeshes.push_back( edgeSM );
199 edgesHolder = SMESH_subMeshEventListenerData::MakeData( edgeSM );
200 faceSM->SetEventListener( TEdgeMarker::getListener(), edgesHolder, faceSM );
206 //================================================================================
208 * \brief Return sides of the face connected in the order: aCircEdge, aLinEdge1, aLinEdge2
209 * \retval int - nb of sides
211 //================================================================================
213 int analyseFace(const TopoDS_Shape& aShape,
215 StdMeshers_FaceSidePtr& aCircSide,
216 StdMeshers_FaceSidePtr& aLinSide1,
217 StdMeshers_FaceSidePtr& aLinSide2,
218 SMESH_MesherHelper* helper)
220 const TopoDS_Face& face = TopoDS::Face( aShape );
221 aCircSide.reset(); aLinSide1.reset(); aLinSide2.reset();
223 list< TopoDS_Edge > edges;
224 list< int > nbEdgesInWire;
225 int nbWire = SMESH_Block::GetOrderedEdges ( face, edges, nbEdgesInWire );
226 if ( nbWire > 2 || nbEdgesInWire.front() < 1 ) return 0;
228 // remove degenerated EDGEs
229 list<TopoDS_Edge>::iterator edge = edges.begin();
230 while ( edge != edges.end() )
231 if ( SMESH_Algo::isDegenerated( *edge ))
232 edge = edges.erase( edge );
235 int nbEdges = edges.size();
237 // find VERTEXes between continues EDGEs
238 TopTools_MapOfShape contVV;
241 TopoDS_Edge ePrev = edges.back();
242 for ( edge = edges.begin(); edge != edges.end(); ++edge )
244 if ( SMESH_Algo::IsContinuous( ePrev, *edge ))
245 contVV.Add( SMESH_MesherHelper::IthVertex( 0, *edge ));
249 // make edges start from a non-continues VERTEX
250 if ( 1 < contVV.Extent() && contVV.Extent() < nbEdges )
252 while ( contVV.Contains( SMESH_MesherHelper::IthVertex( 0, edges.front() )))
253 edges.splice( edges.end(), edges, edges.begin() );
258 while ( !edges.empty() )
260 list< TopoDS_Edge > sideEdges;
261 sideEdges.splice( sideEdges.end(), edges, edges.begin() );
262 while ( !edges.empty() &&
263 contVV.Contains( SMESH_MesherHelper::IthVertex( 0, edges.front() )))
264 sideEdges.splice( sideEdges.end(), edges, edges.begin() );
266 StdMeshers_FaceSidePtr side;
268 side = StdMeshers_FaceSide::New( face, sideEdges, aMesh,
269 /*isFwd=*/true, /*skipMedium=*/ true, helper );
270 sides.push_back( side );
273 if ( !aMesh ) // call from IsApplicable()
276 if ( sides.size() > 3 )
279 if ( nbWire == 2 && (( sides.size() != 2 ) ||
280 ( sides[0]->IsClosed() && sides[1]->IsClosed() ) ||
281 ( !sides[0]->IsClosed() && !sides[1]->IsClosed() )))
284 // detect an elliptic side
286 if ( sides.size() == 1 )
288 aCircSide = sides[0];
292 // sort sides by deviation from a straight line
293 multimap< double, int > deviation2sideInd;
294 const double nbSamples = 7;
295 for ( size_t iS = 0; iS < sides.size(); ++iS )
297 gp_Pnt pf = BRep_Tool::Pnt( sides[iS]->FirstVertex() );
298 gp_Pnt pl = BRep_Tool::Pnt( sides[iS]->LastVertex() );
300 double v1Len = v1.Magnitude();
301 if ( v1Len < std::numeric_limits< double >::min() )
303 deviation2sideInd.insert( make_pair( sides[iS]->Length(), iS )); // the side seems closed
307 for ( int i = 0; i < nbSamples; ++i )
309 gp_Pnt pi( sides[iS]->Value3d(( i + 1 ) / nbSamples ));
311 double h = 0.5 * v1.Crossed( vi ).Magnitude() / v1Len;
312 devia = Max( devia, h );
314 deviation2sideInd.insert( make_pair( devia, iS ));
316 double maxDevi = deviation2sideInd.rbegin()->first;
317 if ( maxDevi < 1e-7 && sides.size() == 3 )
319 // a triangle FACE; use a side with the most outstanding length as an elliptic one
320 deviation2sideInd.clear();
321 multimap< double, int > len2sideInd;
322 for ( size_t iS = 0; iS < sides.size(); ++iS )
323 len2sideInd.insert( make_pair( sides[iS]->Length(), iS ));
325 multimap< double, int >::iterator l2i = len2sideInd.begin();
326 double len0 = l2i->first;
327 double len1 = (++l2i)->first;
328 double len2 = (++l2i)->first;
329 if ( len1 - len0 > len2 - len1 )
330 deviation2sideInd.insert( make_pair( 0., len2sideInd.begin()->second ));
332 deviation2sideInd.insert( make_pair( 0., len2sideInd.rbegin()->second ));
335 int iCirc = deviation2sideInd.rbegin()->second;
336 aCircSide = sides[ iCirc ];
337 aLinSide1 = sides[( iCirc + 1 ) % sides.size() ];
338 if ( sides.size() > 2 )
340 aLinSide2 = sides[( iCirc + 2 ) % sides.size() ];
341 aLinSide2->Reverse(); // to be "parallel" to aLinSide1
344 if (( nbWire == 2 && aLinSide1 ) &&
345 ( aLinSide1->Edge(0).Orientation() == TopAbs_INTERNAL ) &&
346 ( aCircSide->IsClosed() ))
348 // assure that aCircSide starts at aLinSide1
349 TopoDS_Vertex v0 = aLinSide1->FirstVertex();
350 TopoDS_Vertex v1 = aLinSide1->LastVertex();
351 if ( ! aCircSide->FirstVertex().IsSame( v0 ) &&
352 ! aCircSide->FirstVertex().IsSame( v1 ))
355 for ( ; iE < aCircSide->NbEdges(); ++iE )
356 if ( aCircSide->FirstVertex(iE).IsSame( v0 ) ||
357 aCircSide->FirstVertex(iE).IsSame( v1 ))
359 if ( iE == aCircSide->NbEdges() )
363 for ( int i = 0; i < aCircSide->NbEdges(); ++i, ++iE )
364 edges.push_back( aCircSide->Edge( iE % aCircSide->NbEdges() ));
366 aCircSide = StdMeshers_FaceSide::New( face, edges, aMesh,
367 /*isFwd=*/true, /*skipMedium=*/ true, helper );
374 //================================================================================
376 * \brief Checks if the common vertex between LinSide's lies inside the circle
378 * \return bool - false if there are 3 EDGEs and the corner is outside
380 //================================================================================
382 bool isCornerInsideCircle(const StdMeshers_FaceSidePtr& CircSide,
383 const StdMeshers_FaceSidePtr& LinSide1,
384 const StdMeshers_FaceSidePtr& LinSide2)
386 // if ( CircSide && LinSide1 && LinSide2 )
388 // Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircSide ));
389 // TopoDS_Vertex aCommonV;
390 // if ( !aCirc.IsNull() &&
391 // TopExp::CommonVertex( LinSide1, LinSide2, aCommonV ))
393 // gp_Pnt aCommonP = BRep_Tool::Pnt( aCommonV );
394 // gp_Pnt aCenter = aCirc->Location();
395 // double dist = aCenter.Distance( aCommonP );
396 // return dist < 0.1 * aCirc->Radius();
402 //================================================================================
404 * \brief Create an EDGE connecting the ellipse center with the most distant point
406 * \param [in] circSide - the elliptic side
407 * \param [in] face - the FACE
408 * \param [out] circNode - a node on circSide most distant from the center
409 * \return TopoDS_Edge - the create EDGE
411 //================================================================================
413 TopoDS_Edge makeEdgeToCenter( StdMeshers_FaceSidePtr& circSide,
414 const TopoDS_Face& face,
415 const SMDS_MeshNode*& circNode)
417 // find the center and a point most distant from it
419 double maxDist = 0, normPar;
421 for ( int i = 0; i < 32; ++i )
423 double u = 0.5 * i / 32.;
424 gp_Pnt2d p1 = circSide->Value2d( u );
425 gp_Pnt2d p2 = circSide->Value2d( u + 0.5 );
426 double dist = p1.SquareDistance( p2 );
427 if ( dist > maxDist )
435 gp_XY center = 0.5 * ( uv1 + uv2 );
436 double len = 0.5 * Sqrt( maxDist );
437 bool isCirc = ( Abs( len - circSide->Value2d( 0 ).Distance( center )) < 1e-3 * len );
439 // find a node closest to the most distant point
442 const UVPtStructVec& circNodes = circSide->GetUVPtStruct();
445 double minDist = 1e100;
446 for ( size_t i = 0; i <= circNodes.size(); ++i )
448 double dist = Abs( circNodes[i].normParam - normPar );
449 if ( dist < minDist )
456 circNode = circNodes[iDist].node;
457 uv1 = circNodes[iDist].UV();
458 len = ( uv1 - center ).Modulus();
460 // make the EDGE between the most distant point and the center
462 Handle(Geom2d_Line) line = new Geom2d_Line( uv1, gp_Dir2d( center - uv1 ) );
463 Handle(Geom2d_Curve) pcu = new Geom2d_TrimmedCurve( line, 0, len );
465 Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
466 TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pcu, surface, 0, len );
468 BRep_Builder().UpdateEdge( edge, pcu, face, 1e-7 );
469 ShapeFix_Edge().FixAddCurve3d( edge );
472 // assure that circSide starts at circNode
473 if ( iDist != 0 && iDist != circNodes.size()-1 )
475 // create a new circSide
476 UVPtStructVec nodesNew;
477 nodesNew.reserve( circNodes.size() );
478 nodesNew.insert( nodesNew.end(), circNodes.begin() + iDist, circNodes.end() );
479 nodesNew.insert( nodesNew.end(), circNodes.begin() + 1, circNodes.begin() + iDist + 1 );
480 circSide = StdMeshers_FaceSide::New( nodesNew );
486 //================================================================================
488 * \brief Set nodes existing on a linSide to UVPtStructVec and create missing nodes
489 * corresponding to layerPositions
491 //================================================================================
493 void makeMissingMesh( StdMeshers_FaceSidePtr& linSide,
494 UVPtStructVec& nodes,
495 const vector< double >& layerPositions,
496 SMESH_MesherHelper* helper )
498 // tolerance to compare normParam
500 for ( size_t i = 1; i < layerPositions.size(); ++i )
501 tol = Min( tol, layerPositions[i] - layerPositions[i-1] );
504 // merge existing nodes with layerPositions into UVPtStructVec
505 // ------------------------------------------------------------
507 const UVPtStructVec& exiNodes = linSide->GetUVPtStruct();
509 nodes.reserve( layerPositions.size() + exiNodes.size() );
510 vector< double >::const_iterator pos = layerPositions.begin(), posEnd = layerPositions.end();
511 for ( size_t i = 0; i < exiNodes.size(); ++i )
513 switch ( exiNodes[i].node->GetPosition()->GetTypeOfPosition() )
515 case SMDS_TOP_VERTEX:
517 // allocate UVPtStruct's for non-existing nodes
518 while ( pos != posEnd && *pos < exiNodes[i].normParam - tol )
521 uvPS.normParam = *pos++;
522 nodes.push_back( uvPS );
524 // save existing node on a VERTEX
525 nodes.push_back( exiNodes[i] );
530 // save existing nodes on an EDGE
531 while ( i < exiNodes.size() && exiNodes[i].node->GetPosition()->GetDim() == 1 )
533 nodes.push_back( exiNodes[i++] );
535 // save existing node on a VERTEX
536 if ( i < exiNodes.size() && exiNodes[i].node->GetPosition()->GetDim() == 0 )
538 nodes.push_back( exiNodes[i] );
545 // skip layer positions covered by saved nodes
546 while ( pos != posEnd && *pos < nodes.back().normParam + tol )
551 // allocate UVPtStruct's for the rest non-existing nodes
552 while ( pos != posEnd )
555 uvPS.normParam = *pos++;
556 nodes.push_back( uvPS );
559 // create missing nodes
560 // ---------------------
562 SMESHDS_Mesh * meshDS = helper->GetMeshDS();
563 const TopoDS_Face& F = TopoDS::Face( helper->GetSubShape() );
564 Handle(ShapeAnalysis_Surface) surface = helper->GetSurface( F );
566 helper->SetElementsOnShape( false ); // we create nodes on EDGEs, not on the FACE
568 for ( size_t i = 0; i < nodes.size(); ++i )
570 if ( nodes[ i ].node ) continue;
572 gp_Pnt2d uv = linSide->Value2d( nodes[i].normParam );
573 gp_Pnt xyz = surface->Value( uv.X(), uv.Y() );
575 nodes[ i ].SetUV( uv.XY() );
576 nodes[ i ].node = helper->AddNode( xyz.X(), xyz.Y(), xyz.Z() );
579 // set nodes on VERTEXes
580 for ( int iE = 0; iE < linSide->NbEdges(); ++iE )
582 TopoDS_Vertex v = linSide->LastVertex( iE );
583 if ( SMESH_Algo::VertexNode( v, meshDS ))
586 double normPar = linSide->LastParameter( iE );
588 while ( nodes[ i ].normParam < normPar )
590 if (( nodes[ i ].normParam - normPar ) > ( normPar - nodes[ i-1 ].normParam ))
592 meshDS->SetNodeOnVertex( nodes[ i ].node, v );
595 // set nodes on EDGEs
597 for ( size_t i = 0; i < nodes.size(); ++i )
599 if ( nodes[ i ].node->getshapeId() > 0 ) continue;
601 double u = linSide->Parameter( nodes[i].normParam, edgeID );
602 meshDS->SetNodeOnEdge( nodes[ i ].node, edgeID, u );
606 for ( size_t i = 1; i < nodes.size(); ++i )
608 if ( meshDS->FindEdge( nodes[i].node, nodes[i-1].node )) continue;
610 const SMDS_MeshElement* seg = helper->AddEdge( nodes[i].node, nodes[i-1].node );
612 double normParam = 0.5 * ( nodes[i].normParam + nodes[i-1].normParam );
613 edgeID = linSide->EdgeID( linSide->EdgeIndex( normParam ));
614 meshDS->SetMeshElementOnShape( seg, edgeID );
617 helper->SetElementsOnShape( true );
619 } // end makeMissingMesh()
621 //================================================================================
622 //================================================================================
624 * \brief Class computing layers distribution using data of
625 * StdMeshers_LayerDistribution hypothesis
627 //================================================================================
628 //================================================================================
630 class TNodeDistributor: public StdMeshers_Regular_1D
632 list <const SMESHDS_Hypothesis *> myUsedHyps;
634 // -----------------------------------------------------------------------------
635 static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
637 const int myID = -1001;
638 TNodeDistributor* myHyp = dynamic_cast<TNodeDistributor*>( aMesh.GetHypothesis( myID ));
640 myHyp = new TNodeDistributor( myID, aMesh.GetGen() );
643 // -----------------------------------------------------------------------------
644 //! Computes distribution of nodes on a straight line ending at pIn and pOut
645 bool Compute( vector< double > & positions,
646 const TopoDS_Edge& edge,
647 Adaptor3d_Curve& curve,
651 const SMESH_Hypothesis* hyp1d)
653 if ( !hyp1d ) return error( "Invalid LayerDistribution hypothesis");
656 myUsedHyps.push_back( hyp1d );
658 SMESH_Hypothesis::Hypothesis_Status aStatus;
659 if ( !StdMeshers_Regular_1D::CheckHypothesis( mesh, edge, aStatus ))
660 return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
661 "with LayerDistribution hypothesis");
663 double len = GCPnts_AbscissaPoint::Length( curve, f, l );
665 list< double > params;
666 if ( !StdMeshers_Regular_1D::computeInternalParameters( mesh, curve, len, f, l, params, false ))
667 return error("StdMeshers_Regular_1D failed to compute layers distribution");
669 params.push_front( f );
670 params.push_back ( l );
672 positions.reserve( params.size() );
673 for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
674 positions.push_back(( *itU - f ) / ( l - f ));
677 // -----------------------------------------------------------------------------
678 //! Make mesh on an edge using assigned 1d hyp or default nb of segments
679 bool ComputeCircularEdge( SMESH_Mesh& aMesh,
680 const StdMeshers_FaceSidePtr& aSide )
683 for ( int i = 0; i < aSide->NbEdges(); ++i )
685 const TopoDS_Edge& edge = aSide->Edge( i );
686 _gen->Compute( aMesh, edge );
687 SMESH_subMesh *sm = aMesh.GetSubMesh( edge );
688 if ( sm->GetComputeState() != SMESH_subMesh::COMPUTE_OK)
690 // find any 1d hyp assigned (there can be a hyp w/o algo)
691 myUsedHyps = SMESH_Algo::GetUsedHypothesis( aMesh, edge, /*ignoreAux=*/true );
692 Hypothesis_Status aStatus;
693 if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
695 // no valid 1d hyp assigned, use default nb of segments
696 _hypType = NB_SEGMENTS;
697 _ivalue[ DISTR_TYPE_IND ] = StdMeshers_NumberOfSegments::DT_Regular;
698 _ivalue[ NB_SEGMENTS_IND ] = _gen->GetDefaultNbSegments();
700 ok &= StdMeshers_Regular_1D::Compute( aMesh, edge );
705 // -----------------------------------------------------------------------------
706 //! Make mesh on an edge using assigned 1d hyp or default nb of segments
707 bool EvaluateCircularEdge(SMESH_Mesh& aMesh,
708 const StdMeshers_FaceSidePtr aSide,
709 MapShapeNbElems& aResMap)
712 for ( int i = 0; i < aSide->NbEdges(); ++i )
714 const TopoDS_Edge& anEdge = aSide->Edge( i );
715 _gen->Evaluate( aMesh, anEdge, aResMap );
716 if ( aResMap.count( aMesh.GetSubMesh( anEdge )))
719 // find any 1d hyp assigned
720 myUsedHyps = SMESH_Algo::GetUsedHypothesis(aMesh, anEdge, /*ignoreAux=*/true);
721 Hypothesis_Status aStatus;
722 if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, anEdge, aStatus ))
724 // no valid 1d hyp assigned, use default nb of segments
725 _hypType = NB_SEGMENTS;
726 _ivalue[ DISTR_TYPE_IND ] = StdMeshers_NumberOfSegments::DT_Regular;
727 _ivalue[ NB_SEGMENTS_IND ] = _gen->GetDefaultNbSegments();
729 ok &= StdMeshers_Regular_1D::Evaluate( aMesh, anEdge, aResMap );
734 // -----------------------------------------------------------------------------
735 TNodeDistributor( int hypId, SMESH_Gen* gen)
736 : StdMeshers_Regular_1D( hypId, gen)
739 // -----------------------------------------------------------------------------
740 virtual const list <const SMESHDS_Hypothesis *> &
741 GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
745 // -----------------------------------------------------------------------------
749 //=======================================================================
751 * \brief Allow algo to do something after persistent restoration
752 * \param subMesh - restored submesh
754 * call TEdgeMarker::markEdge()
756 //=======================================================================
758 void StdMeshers_RadialQuadrangle_1D2D::SubmeshRestored(SMESH_subMesh* faceSubMesh)
760 if ( !faceSubMesh->IsEmpty() )
762 for ( TopExp_Explorer e( faceSubMesh->GetSubShape(), TopAbs_EDGE ); e.More(); e.Next() )
764 TEdgeMarker::markEdge( TopoDS::Edge( e.Current() ), faceSubMesh );
769 //=======================================================================
772 //=======================================================================
774 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
775 const TopoDS_Shape& aShape)
777 SMESH_MesherHelper helper( aMesh );
778 StdMeshers_Quadrangle_2D::myHelper = & helper;
779 StdMeshers_Quadrangle_2D::myNeedSmooth = false;
780 StdMeshers_Quadrangle_2D::myCheckOri = false;
781 StdMeshers_Quadrangle_2D::myQuadList.clear();
783 myHelper->SetSubShape( aShape );
784 myHelper->SetElementsOnShape( true );
786 StdMeshers_FaceSidePtr circSide, linSide1, linSide2;
787 int nbSides = analyseFace( aShape, &aMesh, circSide, linSide1, linSide2, myHelper );
788 if( nbSides > 3 || nbSides < 1 )
789 return error("The face must be a full ellipse or a part of ellipse (i.e. the number "
790 "of edges is less or equal to 3 and one of them is an ellipse curve)");
792 // get not yet computed EDGEs
793 list< TopoDS_Edge > emptyEdges;
794 for ( TopExp_Explorer e( aShape, TopAbs_EDGE ); e.More(); e.Next() )
796 if ( aMesh.GetSubMesh( e.Current() )->IsEmpty() )
797 emptyEdges.push_back( TopoDS::Edge( e.Current() ));
800 TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
802 if ( !algo1d->ComputeCircularEdge( aMesh, circSide ))
803 return error( algo1d->GetComputeError() );
806 TopoDS_Face F = TopoDS::Face(aShape);
807 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
809 myHelper->IsQuadraticSubMesh( aShape );
811 vector< double > layerPositions; // [0,1]
813 const SMDS_MeshNode* centerNode = 0;
814 gp_Pnt2d centerUV(0,0);
816 // ------------------------------------------------------------------------------------------
817 if ( nbSides == 1 ) // full ellipse
819 const SMDS_MeshNode* circNode;
820 TopoDS_Edge linEdge = makeEdgeToCenter( circSide, F, circNode );
822 StdMeshers_FaceSidePtr tmpSide =
823 StdMeshers_FaceSide::New( F, linEdge, &aMesh, /*isFrw=*/true, /*skipMedium=*/true, myHelper );
825 if ( !computeLayerPositions( tmpSide, layerPositions ))
828 UVPtStructVec nodes( layerPositions.size() );
829 nodes[0].node = circNode;
830 for ( size_t i = 0; i < layerPositions.size(); ++i )
832 gp_Pnt2d uv = tmpSide->Value2d( layerPositions[i] );
833 gp_Pnt xyz = S->Value( uv.X(), uv.Y() );
835 nodes[ i ].SetUV( uv.XY() );
836 nodes[ i ].normParam = layerPositions[i];
838 nodes[ i ].node = myHelper->AddNode( xyz.X(), xyz.Y(), xyz.Z(), 0, uv.X(), uv.Y() );
841 linSide1 = StdMeshers_FaceSide::New( nodes );
842 linSide2 = StdMeshers_FaceSide::New( nodes );
844 centerNode = nodes.back().node;
845 centerUV = nodes.back().UV();
847 // ------------------------------------------------------------------------------------------
848 else if ( nbSides == 2 && linSide1->Edge(0).Orientation() == TopAbs_INTERNAL )
850 // full ellipse with an internal radial side
852 // eliminate INTERNAL orientation
853 list< TopoDS_Edge > edges;
854 for ( int iE = 0; iE < linSide1->NbEdges(); ++iE )
856 edges.push_back( linSide1->Edge( iE ));
857 edges.back().Orientation( TopAbs_FORWARD );
860 // orient the internal side
861 bool isVIn0Shared = false;
862 TopoDS_Vertex vIn0 = myHelper->IthVertex( 0, edges.front() );
863 for ( int iE = 0; iE < circSide->NbEdges() && !isVIn0Shared; ++iE )
864 isVIn0Shared = vIn0.IsSame( circSide->FirstVertex( iE ));
866 linSide1 = StdMeshers_FaceSide::New( F, edges, &aMesh,
867 /*isFrw=*/isVIn0Shared, /*skipMedium=*/true, myHelper );
870 if ( !computeLayerPositions( linSide1, layerPositions, &nbMeshedEdges ))
873 // merge existing nodes with new nodes at layerPositions into a UVPtStructVec
875 if ( nbMeshedEdges != linSide1->NbEdges() )
876 makeMissingMesh( linSide1, nodes, layerPositions, myHelper );
878 nodes = linSide1->GetUVPtStruct();
880 linSide1 = StdMeshers_FaceSide::New( nodes );
881 linSide2 = StdMeshers_FaceSide::New( nodes );
883 centerNode = nodes.back().node;
884 centerUV = nodes.back().UV();
886 // ------------------------------------------------------------------------------------------
887 else if ( nbSides == 2 )
889 // find positions of layers for the first half of linSide1
891 if ( !computeLayerPositions( linSide1, layerPositions, &nbMeshedEdges, /*useHalf=*/true ))
894 // make positions for the whole linSide1
895 for ( size_t i = 0; i < layerPositions.size(); ++i )
897 layerPositions[i] *= 0.5;
899 layerPositions.reserve( layerPositions.size() * 2 );
900 for ( int nb = layerPositions.size()-1; nb > 0; --nb )
901 layerPositions.push_back( layerPositions.back() + layerPositions[nb] - layerPositions[nb-1] );
903 // merge existing nodes with new nodes at layerPositions into a UVPtStructVec
905 if ( nbMeshedEdges != linSide1->NbEdges() )
906 makeMissingMesh( linSide1, nodes, layerPositions, myHelper );
908 nodes = linSide1->GetUVPtStruct();
910 // find a central node
912 while ( nodes[ i ].normParam < 0.5 ) ++i;
913 if (( nodes[ i ].normParam - 0.5 ) > ( 0.5 - nodes[ i-1 ].normParam )) --i;
915 // distribute nodes between two linear sides
916 UVPtStructVec nodes2( nodes.rbegin(), nodes.rbegin() + nodes.size() - i );
917 nodes.resize( i + 1 );
919 linSide1 = StdMeshers_FaceSide::New( nodes );
920 linSide2 = StdMeshers_FaceSide::New( nodes2 );
922 centerNode = nodes.back().node;
923 centerUV = nodes.back().UV();
925 // ------------------------------------------------------------------------------------------
928 // one curve must be a part of ellipse and 2 other curves must be segments of line
930 int nbMeshedEdges1, nbMeshedEdges2;
931 vector< double > layerPositions2;
932 bool ok1 = computeLayerPositions( linSide1, layerPositions, &nbMeshedEdges1 );
933 bool ok2 = computeLayerPositions( linSide2, layerPositions2, &nbMeshedEdges2 );
937 bool linSide1Computed = ( nbMeshedEdges1 == linSide1->NbEdges() );
938 bool linSide2Computed = ( nbMeshedEdges2 == linSide2->NbEdges() );
942 if ( linSide1Computed && !linSide2Computed )
944 // use layer positions of linSide1 to mesh linSide2
945 makeMissingMesh( linSide2, nodes, layerPositions, myHelper );
946 linSide2 = StdMeshers_FaceSide::New( nodes );
948 else if ( linSide2Computed && !linSide1Computed )
950 // use layer positions of linSide2 to mesh linSide1
951 makeMissingMesh( linSide1, nodes, layerPositions2, myHelper );
952 linSide1 = StdMeshers_FaceSide::New( nodes );
954 else if ( !linSide2Computed && !linSide1Computed )
956 // use layer positions of a longer side to mesh the shorter side
957 vector< double >& lp =
958 ( linSide1->Length() > linSide2->Length() ) ? layerPositions : layerPositions2;
960 makeMissingMesh( linSide1, nodes, lp, myHelper );
961 linSide1 = StdMeshers_FaceSide::New( nodes );
962 makeMissingMesh( linSide2, nodes, lp, myHelper );
963 linSide2 = StdMeshers_FaceSide::New( nodes );
966 const UVPtStructVec& nodes2 = linSide2->GetUVPtStruct();
967 centerNode = nodes2.back().node;
968 centerUV = nodes2.back().UV();
971 list< TopoDS_Edge >::iterator ee = emptyEdges.begin();
972 for ( ; ee != emptyEdges.end(); ++ee )
973 TEdgeMarker::markEdge( *ee, aMesh.GetSubMesh( F ));
975 circSide->GetUVPtStruct(); // let sides take into account just computed nodes
976 linSide1->GetUVPtStruct();
977 linSide2->GetUVPtStruct();
979 FaceQuadStruct::Ptr quad( new FaceQuadStruct );
980 quad->side.resize( 4 );
981 quad->side[0] = circSide;
982 quad->side[1] = linSide1;
983 quad->side[2] = StdMeshers_FaceSide::New( circSide.get(), centerNode, ¢erUV );
984 quad->side[3] = linSide2;
986 myQuadList.push_back( quad );
988 // create quadrangles
990 if ( linSide1->NbPoints() == linSide2->NbPoints() )
991 ok = StdMeshers_Quadrangle_2D::computeQuadDominant( aMesh, F, quad );
993 ok = StdMeshers_Quadrangle_2D::computeTriangles( aMesh, F, quad );
995 StdMeshers_Quadrangle_2D::myHelper = 0;
1000 //================================================================================
1002 * \brief Compute nodes on the radial edge
1003 * \retval int - nb of segments on the linSide
1005 //================================================================================
1007 int StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(StdMeshers_FaceSidePtr linSide,
1008 vector< double > & positions,
1012 // First, try to compute positions of layers
1016 SMESH_Mesh * mesh = myHelper->GetMesh();
1018 const SMESH_Hypothesis* hyp1D = myDistributionHypo ? myDistributionHypo->GetLayerDistribution() : 0;
1019 int nbLayers = myNbLayerHypo ? myNbLayerHypo->GetNumberOfLayers() : 0;
1021 if ( !hyp1D && !nbLayers )
1023 // No own algo hypotheses assigned, so first try to find any 1D hypothesis.
1024 // find a hyp usable by TNodeDistributor
1025 TopoDS_Shape edge = linSide->Edge(0);
1026 const SMESH_HypoFilter* hypKind =
1027 TNodeDistributor::GetDistributor(*mesh)->GetCompatibleHypoFilter(/*ignoreAux=*/true);
1028 hyp1D = mesh->GetHypothesis( edge, *hypKind, /*fromAncestors=*/true);
1030 if ( hyp1D ) // try to compute with hyp1D
1032 BRepAdaptor_CompCurve* curve = linSide->GetCurve3d();
1033 SMESHUtils::Deleter< BRepAdaptor_CompCurve > delCurve( curve );
1034 double f = curve->FirstParameter();
1035 double l = curve->LastParameter();
1038 l = 0.5 * ( f + l ); // first part of linSide is used
1040 if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( positions, linSide->Edge(0),
1041 *curve, f, l, *mesh, hyp1D ))
1043 if ( myDistributionHypo ) { // bad hyp assigned
1044 return error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
1047 // bad hyp found, its Ok, lets try with default nb of segments
1052 if ( positions.empty() ) // try to use nb of layers
1055 nbLayers = _gen->GetDefaultNbSegments();
1059 positions.resize( nbLayers + 1 );
1060 for ( int z = 0; z < nbLayers; ++z )
1061 positions[ z ] = double( z )/ double( nbLayers );
1062 positions.back() = 1;
1066 // Second, check presence of a mesh built by other algo on linSide
1068 int nbEdgesComputed = 0;
1069 for ( int i = 0; i < linSide->NbEdges(); ++i )
1071 nbEdgesComputed += ( !mesh->GetSubMesh( linSide->Edge(i))->IsEmpty() );
1074 if ( nbEdgesComputed == linSide->NbEdges() )
1076 const UVPtStructVec& points = linSide->GetUVPtStruct();
1077 if ( points.size() >= 2 )
1079 positions.resize( points.size() );
1080 for ( size_t i = 0; i < points.size(); ++i )
1081 positions[ i ] = points[i].normParam;
1085 if ( nbMeshedEdges ) *nbMeshedEdges = nbEdgesComputed;
1087 return positions.size();
1091 //=======================================================================
1092 //function : Evaluate
1094 //=======================================================================
1096 bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
1097 const TopoDS_Shape& aShape,
1098 MapShapeNbElems& aResMap)
1100 if( aShape.ShapeType() != TopAbs_FACE ) {
1103 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1104 if( aResMap.count(sm) )
1107 vector<int>& aResVec =
1108 aResMap.insert( make_pair(sm, vector<int>(SMDSEntity_Last,0))).first->second;
1110 myHelper = new SMESH_MesherHelper( aMesh );
1111 myHelper->SetSubShape( aShape );
1112 SMESHUtils::Deleter<SMESH_MesherHelper> helperDeleter( myHelper );
1114 TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
1116 StdMeshers_FaceSidePtr circSide, linSide1, linSide2;
1117 int nbSides = analyseFace( aShape, &aMesh, circSide, linSide1, linSide2, myHelper );
1118 if( nbSides > 3 || nbSides < 1 )
1121 if ( algo1d->EvaluateCircularEdge( aMesh, circSide, aResMap ))
1124 vector< double > layerPositions; // [0,1]
1126 // ------------------------------------------------------------------------------------------
1129 const TopoDS_Face& F = TopoDS::Face( aShape );
1131 const SMDS_MeshNode* circNode;
1132 TopoDS_Edge linEdge = makeEdgeToCenter( circSide, F, circNode );
1134 StdMeshers_FaceSidePtr tmpSide =
1135 StdMeshers_FaceSide::New( F, linEdge, &aMesh, /*isFrw=*/true, /*skipMedium=*/true, myHelper );
1137 if ( !computeLayerPositions( tmpSide, layerPositions ))
1140 // ------------------------------------------------------------------------------------------
1141 else if ( nbSides == 2 && linSide1->Edge(0).Orientation() == TopAbs_INTERNAL )
1143 if ( !computeLayerPositions( linSide1, layerPositions ))
1146 // ------------------------------------------------------------------------------------------
1147 else if ( nbSides == 2 )
1149 // find positions of layers for the first half of linSide1
1150 if ( !computeLayerPositions( linSide1, layerPositions, 0, /*useHalf=*/true ))
1153 // ------------------------------------------------------------------------------------------
1154 else // nbSides == 3
1156 if ( !computeLayerPositions(( linSide1->Length() > linSide2->Length() ) ? linSide1 : linSide2,
1161 bool isQuadratic = false;
1162 for ( TopExp_Explorer edge( aShape, TopAbs_EDGE ); edge.More() && !isQuadratic ; edge.Next() )
1164 sm = aMesh.GetSubMesh( edge.Current() );
1165 vector<int>& nbElems = aResMap[ sm ];
1166 if ( SMDSEntity_Quad_Edge < (int) nbElems.size() )
1167 isQuadratic = nbElems[ SMDSEntity_Quad_Edge ];
1170 int nbCircSegments = 0;
1171 for ( int iE = 0; iE < circSide->NbEdges(); ++iE )
1173 sm = aMesh.GetSubMesh( circSide->Edge( iE ));
1174 vector<int>& nbElems = aResMap[ sm ];
1175 if ( SMDSEntity_Quad_Edge < (int) nbElems.size() )
1176 nbCircSegments += ( nbElems[ SMDSEntity_Edge ] + nbElems[ SMDSEntity_Quad_Edge ]);
1179 int nbQuads = nbCircSegments * ( layerPositions.size() - 1 );
1180 int nbTria = nbCircSegments;
1181 int nbNodes = ( nbCircSegments - 1 ) * ( layerPositions.size() - 2 );
1184 nbNodes += (( nbCircSegments - 1 ) * ( layerPositions.size() - 1 ) + // radial
1185 ( nbCircSegments ) * ( layerPositions.size() - 2 )); // circular
1186 aResVec[SMDSEntity_Quad_Triangle ] = nbTria;
1187 aResVec[SMDSEntity_Quad_Quadrangle] = nbQuads;
1191 aResVec[SMDSEntity_Triangle ] = nbTria;
1192 aResVec[SMDSEntity_Quadrangle] = nbQuads;
1194 aResVec[SMDSEntity_Node] = nbNodes;
1198 // evaluation for linSides
1199 vector<int> aResVec(SMDSEntity_Last, 0);
1200 if ( isQuadratic ) {
1201 aResVec[SMDSEntity_Node ] = 2 * ( layerPositions.size() - 1 ) + 1;
1202 aResVec[SMDSEntity_Quad_Edge] = layerPositions.size() - 1;
1205 aResVec[SMDSEntity_Node] = layerPositions.size() - 2;
1206 aResVec[SMDSEntity_Edge] = layerPositions.size() - 1;
1208 sm = aMesh.GetSubMesh( linSide1->Edge(0) );
1209 aResMap[sm] = aResVec;
1212 sm = aMesh.GetSubMesh( linSide2->Edge(0) );
1213 aResMap[sm] = aResVec;
1220 //================================================================================
1222 * \brief Return true if the algorithm can compute mesh on this shape
1224 //================================================================================
1226 bool StdMeshers_RadialQuadrangle_1D2D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
1228 int nbFoundFaces = 0;
1229 for (TopExp_Explorer exp( aShape, TopAbs_FACE ); exp.More(); exp.Next(), ++nbFoundFaces )
1231 StdMeshers_FaceSidePtr circSide, linSide1, linSide2;
1232 int nbSides = analyseFace( exp.Current(), NULL, circSide, linSide1, linSide2, NULL );
1233 bool ok = ( 0 < nbSides && nbSides <= 3 &&
1234 isCornerInsideCircle( circSide, linSide1, linSide2 ));
1235 if( toCheckAll && !ok ) return false;
1236 if( !toCheckAll && ok ) return true;
1238 if( toCheckAll && nbFoundFaces != 0 ) return true;