1 // Copyright (C) 2007-2008 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 // SMESH SMESH : implementaion of SMESH idl descriptions
23 // File : StdMeshers_RadialQuadrangle_1D2D.cxx
25 // Created : Fri Oct 20 11:37:07 2006
26 // Author : Edward AGAPOV (eap)
28 #include "StdMeshers_RadialQuadrangle_1D2D.hxx"
30 //#include "StdMeshers_ProjectionUtils.hxx"
31 #include "StdMeshers_NumberOfLayers.hxx"
32 #include "StdMeshers_LayerDistribution.hxx"
33 //#include "StdMeshers_Prism_3D.hxx"
34 #include "StdMeshers_Regular_1D.hxx"
36 #include "SMDS_MeshNode.hxx"
37 #include "SMESHDS_SubMesh.hxx"
38 #include "SMESH_Gen.hxx"
39 #include "SMESH_Mesh.hxx"
40 #include "SMESH_MesherHelper.hxx"
41 #include "SMESH_subMesh.hxx"
42 #include "SMESH_subMeshEventListener.hxx"
44 #include "utilities.h"
46 #include <BRepAdaptor_Curve.hxx>
47 #include <BRepBuilderAPI_MakeEdge.hxx>
48 #include <BRep_Tool.hxx>
49 #include <GeomAPI_ProjectPointOnSurf.hxx>
50 #include <Geom_Circle.hxx>
51 #include <Geom_Line.hxx>
52 #include <Geom_TrimmedCurve.hxx>
53 #include <TColgp_SequenceOfPnt.hxx>
54 #include <TColgp_SequenceOfPnt2d.hxx>
55 #include <TopExp_Explorer.hxx>
61 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
62 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
64 //typedef StdMeshers_ProjectionUtils TAssocTool;
67 //=======================================================================
68 //function : StdMeshers_RadialQuadrangle_1D2D
70 //=======================================================================
72 StdMeshers_RadialQuadrangle_1D2D::StdMeshers_RadialQuadrangle_1D2D(int hypId,
75 :SMESH_2D_Algo(hypId, studyId, gen)
77 _name = "RadialQuadrangle_1D2D";
78 _shapeType = (1 << TopAbs_FACE); // 1 bit per shape type
80 _compatibleHypothesis.push_back("LayerDistribution2D");
81 _compatibleHypothesis.push_back("NumberOfLayers2D");
83 myDistributionHypo = 0;
84 _requireDescretBoundary = false;
85 _supportSubmeshes = true;
89 //================================================================================
93 //================================================================================
95 StdMeshers_RadialQuadrangle_1D2D::~StdMeshers_RadialQuadrangle_1D2D()
99 //=======================================================================
100 //function : CheckHypothesis
102 //=======================================================================
104 bool StdMeshers_RadialQuadrangle_1D2D::CheckHypothesis
106 const TopoDS_Shape& aShape,
107 SMESH_Hypothesis::Hypothesis_Status& aStatus)
111 myDistributionHypo = 0;
113 list <const SMESHDS_Hypothesis * >::const_iterator itl;
115 const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
116 if ( hyps.size() == 0 ) {
117 aStatus = SMESH_Hypothesis::HYP_MISSING;
118 return false; // can't work with no hypothesis
121 if ( hyps.size() > 1 ) {
122 aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
126 const SMESHDS_Hypothesis *theHyp = hyps.front();
128 string hypName = theHyp->GetName();
130 if (hypName == "NumberOfLayers2D") {
131 myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
132 aStatus = SMESH_Hypothesis::HYP_OK;
135 if (hypName == "LayerDistribution2D") {
136 myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
137 aStatus = SMESH_Hypothesis::HYP_OK;
140 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
146 // ------------------------------------------------------------------------------
148 * \brief Listener used to mark edges meshed by StdMeshers_RadialQuadrangle_1D2D
150 class TLinEdgeMarker : public SMESH_subMeshEventListener
152 TLinEdgeMarker(): SMESH_subMeshEventListener(/*isDeletable=*/false) {}
154 static SMESH_subMeshEventListener* getListener()
156 static TLinEdgeMarker theEdgeMarker;
157 return &theEdgeMarker;
161 // ------------------------------------------------------------------------------
163 * \brief Mark an edge as computed by StdMeshers_RadialQuadrangle_1D2D
165 void markLinEdgeAsComputedByMe(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
167 if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
169 if ( !edgeSM->GetEventListenerData( TLinEdgeMarker::getListener() ))
170 faceSubMesh->SetEventListener( TLinEdgeMarker::getListener(),
171 SMESH_subMeshEventListenerData::MakeData(faceSubMesh),
175 // ------------------------------------------------------------------------------
177 * \brief Return true if a radial edge was meshed with StdMeshers_RadialQuadrangle_1D2D with
178 * the same radial distribution
180 bool isEdgeCompitaballyMeshed(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
182 if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
184 if ( SMESH_subMeshEventListenerData* otherFaceData =
185 edgeSM->GetEventListenerData( TLinEdgeMarker::getListener() ))
187 // compare hypothesis aplied to two disk faces sharing radial edges
188 SMESH_Mesh& mesh = *faceSubMesh->GetFather();
189 SMESH_Algo* radialQuadAlgo = mesh.GetGen()->GetAlgo(mesh, faceSubMesh->GetSubShape() );
190 SMESH_subMesh* otherFaceSubMesh = otherFaceData->mySubMeshes.front();
191 const list <const SMESHDS_Hypothesis *> & hyps1 =
192 radialQuadAlgo->GetUsedHypothesis( mesh, faceSubMesh->GetSubShape());
193 const list <const SMESHDS_Hypothesis *> & hyps2 =
194 radialQuadAlgo->GetUsedHypothesis( mesh, otherFaceSubMesh->GetSubShape());
195 if( hyps1.empty() && hyps2.empty() )
196 return true; // defaul hyps
197 if ( hyps1.size() != hyps2.size() ||
198 strcmp( hyps1.front()->GetName(), hyps2.front()->GetName() ))
200 ostringstream hypDump1, hypDump2;
201 list <const SMESHDS_Hypothesis*>::const_iterator hyp1 = hyps1.begin();
202 for ( ; hyp1 != hyps1.end(); ++hyp1 )
203 const_cast<SMESHDS_Hypothesis*>(*hyp1)->SaveTo( hypDump1 );
204 list <const SMESHDS_Hypothesis*>::const_iterator hyp2 = hyps2.begin();
205 for ( ; hyp2 != hyps2.end(); ++hyp2 )
206 const_cast<SMESHDS_Hypothesis*>(*hyp2)->SaveTo( hypDump2 );
207 return hypDump1.str() == hypDump2.str();
213 //================================================================================
215 * \brief Return base curve of the edge and extremum parameters
217 //================================================================================
219 Handle(Geom_Curve) getCurve(const TopoDS_Edge& edge, double* f=0, double* l=0)
221 Handle(Geom_Curve) C;
222 if ( !edge.IsNull() )
224 double first = 0., last = 0.;
225 C = BRep_Tool::Curve(edge, first, last);
228 Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C);
229 while( !tc.IsNull() ) {
230 C = tc->BasisCurve();
231 tc = Handle(Geom_TrimmedCurve)::DownCast(C);
240 //================================================================================
242 * \brief Return edges of the face
243 * \retval int - nb of edges
245 //================================================================================
247 int analyseFace(const TopoDS_Shape& face,
248 TopoDS_Edge& CircEdge,
249 TopoDS_Edge& LinEdge1,
250 TopoDS_Edge& LinEdge2)
252 CircEdge.Nullify(); LinEdge1.Nullify(); LinEdge2.Nullify();
255 for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next(), ++nbe )
257 const TopoDS_Edge& E = TopoDS::Edge( exp.Current() );
259 Handle(Geom_Curve) C = getCurve(E,&f,&l);
262 if ( C->IsKind( STANDARD_TYPE(Geom_Circle)))
264 if ( CircEdge.IsNull() )
269 else if ( LinEdge1.IsNull() )
279 //=======================================================================
281 * \brief Allow algo to do something after persistent restoration
282 * \param subMesh - restored submesh
284 * call markLinEdgeAsComputedByMe()
286 //=======================================================================
288 void StdMeshers_RadialQuadrangle_1D2D::SubmeshRestored(SMESH_subMesh* faceSubMesh)
290 if ( !faceSubMesh->IsEmpty() )
292 TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
293 analyseFace( faceSubMesh->GetSubShape(), CircEdge, LinEdge1, LinEdge2 );
294 if ( !LinEdge1.IsNull() ) markLinEdgeAsComputedByMe( LinEdge1, faceSubMesh );
295 if ( !LinEdge2.IsNull() ) markLinEdgeAsComputedByMe( LinEdge2, faceSubMesh );
299 //=======================================================================
302 //=======================================================================
304 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
305 const TopoDS_Shape& aShape)
308 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
310 myHelper = new SMESH_MesherHelper( aMesh );
311 myHelper->IsQuadraticSubMesh( aShape );
312 // to delete helper at exit from Compute()
313 auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
315 myLayerPositions.clear();
317 TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
318 int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 );
319 if( nbe>3 || nbe < 1 || CircEdge.IsNull() )
320 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)");
323 // points for rotation
324 TColgp_SequenceOfPnt Points;
325 // angles for rotation
326 TColStd_SequenceOfReal Angles;
327 // Nodes1 and Nodes2 - nodes along radiuses
328 // CNodes - nodes on circle edge
329 vector< const SMDS_MeshNode* > Nodes1, Nodes2, CNodes;
331 // parameters edge nodes on face
332 TColgp_SequenceOfPnt2d Pnts2d1;
335 int faceID = meshDS->ShapeToIndex(aShape);
336 TopoDS_Face F = TopoDS::Face(aShape);
337 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
341 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
343 bool ok = _gen->Compute( aMesh, CircEdge );
344 if( !ok ) return false;
345 map< double, const SMDS_MeshNode* > theNodes;
346 ok = GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
347 if( !ok ) return false;
350 map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
351 const SMDS_MeshNode* NF = (*itn).second;
352 CNodes.push_back( (*itn).second );
353 double fang = (*itn).first;
354 if ( itn != theNodes.end() ) {
356 for(; itn != theNodes.end(); itn++ ) {
357 CNodes.push_back( (*itn).second );
358 double ang = (*itn).first - fang;
359 if( ang>PI ) ang = ang - 2*PI;
360 if( ang<-PI ) ang = ang + 2*PI;
361 Angles.Append( ang );
364 P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
365 P0 = aCirc->Location();
367 myLayerPositions.clear();
368 computeLayerPositions(P0,P1);
370 exp.Init( CircEdge, TopAbs_VERTEX );
371 TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
372 gp_Pnt2d p2dV = BRep_Tool::Parameters( V1, TopoDS::Face(aShape) );
374 NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
375 GeomAPI_ProjectPointOnSurf PPS(P0,S);
377 PPS.Parameters(1,U0,V0);
378 meshDS->SetNodeOnFace(NC, faceID, U0, V0);
379 PC = gp_Pnt2d(U0,V0);
382 gp_Vec2d aVec2d(PC,p2dV);
383 Nodes1.resize( myLayerPositions.size()+1 );
384 Nodes2.resize( myLayerPositions.size()+1 );
386 for(; i<myLayerPositions.size(); i++) {
387 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
388 P0.Y() + aVec.Y()*myLayerPositions[i],
389 P0.Z() + aVec.Z()*myLayerPositions[i] );
391 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
394 double U = PC.X() + aVec2d.X()*myLayerPositions[i];
395 double V = PC.Y() + aVec2d.Y()*myLayerPositions[i];
396 meshDS->SetNodeOnFace( node, faceID, U, V );
397 Pnts2d1.Append(gp_Pnt2d(U,V));
399 Nodes1[Nodes1.size()-1] = NF;
400 Nodes2[Nodes1.size()-1] = NF;
402 else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL )
404 // one curve must be a half of circle and other curve must be
407 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
408 if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) {
409 // not half of circle
410 return error(COMPERR_BAD_SHAPE);
412 Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
413 if( aLine.IsNull() ) {
414 // other curve not line
415 return error(COMPERR_BAD_SHAPE);
417 bool linEdgeComputed = false;
418 if( SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1) ) {
419 if( !sm1->IsEmpty() )
420 if( isEdgeCompitaballyMeshed( LinEdge1, aMesh.GetSubMesh(F) ))
421 linEdgeComputed = true;
423 return error("Invalid set of hypotheses");
426 bool ok = _gen->Compute( aMesh, CircEdge );
427 if( !ok ) return false;
428 map< double, const SMDS_MeshNode* > theNodes;
429 GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
432 map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
433 double fang = (*itn).first;
435 for(; itn != theNodes.end(); itn++ ) {
436 CNodes.push_back( (*itn).second );
437 double ang = (*itn).first - fang;
438 if( ang>PI ) ang = ang - 2*PI;
439 if( ang<-PI ) ang = ang + 2*PI;
440 Angles.Append( ang );
442 const SMDS_MeshNode* NF = theNodes.begin()->second;
443 const SMDS_MeshNode* NL = theNodes.rbegin()->second;
444 CNodes.push_back( NF );
445 P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
446 gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
447 P0 = aCirc->Location();
449 myLayerPositions.clear();
450 computeLayerPositions(P0,P1);
452 if ( linEdgeComputed )
454 if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
455 return error("Invalid mesh on a straight edge");
457 vector< const SMDS_MeshNode* > *pNodes1 = &Nodes1, *pNodes2 = &Nodes2;
458 bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
459 if ( !nodesFromP0ToP1 ) std::swap( pNodes1, pNodes2 );
461 map< double, const SMDS_MeshNode* >::reverse_iterator ritn = theNodes.rbegin();
462 itn = theNodes.begin();
463 for ( int i = Nodes1.size()-1; i > -1; ++itn, ++ritn, --i )
465 (*pNodes1)[i] = ritn->second;
466 (*pNodes2)[i] = itn->second;
467 Points.Append( gpXYZ( Nodes1[i]));
468 Pnts2d1.Append( myHelper->GetNodeUV( F, Nodes1[i]));
470 NC = const_cast<SMDS_MeshNode*>( itn->second );
471 Points.Remove( Nodes1.size() );
476 int edgeID = meshDS->ShapeToIndex(LinEdge1);
478 Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
482 if( P1.Distance(Ptmp) > Precision::Confusion() )
484 // get UV points for edge
486 BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
487 PC = gp_Pnt2d( (PF.X()+PL.X())/2, (PF.Y()+PL.Y())/2 );
489 if(ori) V2d = gp_Vec2d(PC,PF);
490 else V2d = gp_Vec2d(PC,PL);
492 double cp = (fp+lp)/2;
493 double dp2 = (lp-fp)/2;
494 NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
495 meshDS->SetNodeOnEdge(NC, edgeID, cp);
496 Nodes1.resize( myLayerPositions.size()+1 );
497 Nodes2.resize( myLayerPositions.size()+1 );
499 for(; i<myLayerPositions.size(); i++) {
500 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
501 P0.Y() + aVec.Y()*myLayerPositions[i],
502 P0.Z() + aVec.Z()*myLayerPositions[i] );
504 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
508 param = fp + dp2*(1-myLayerPositions[i]);
510 param = cp + dp2*myLayerPositions[i];
511 meshDS->SetNodeOnEdge(node, edgeID, param);
512 P = gp_Pnt( P0.X() - aVec.X()*myLayerPositions[i],
513 P0.Y() - aVec.Y()*myLayerPositions[i],
514 P0.Z() - aVec.Z()*myLayerPositions[i] );
515 node = meshDS->AddNode(P.X(), P.Y(), P.Z());
518 param = fp + dp2*(1-myLayerPositions[i]);
520 param = cp + dp2*myLayerPositions[i];
521 meshDS->SetNodeOnEdge(node, edgeID, param);
522 // parameters on face
523 gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
524 PC.Y() + V2d.Y()*myLayerPositions[i] );
527 Nodes1[ myLayerPositions.size() ] = NF;
528 Nodes2[ myLayerPositions.size() ] = NL;
529 // create 1D elements on edge
530 vector< const SMDS_MeshNode* > tmpNodes;
531 tmpNodes.resize(2*Nodes1.size()+1);
532 for(i=0; i<Nodes2.size(); i++)
533 tmpNodes[Nodes2.size()-i-1] = Nodes2[i];
534 tmpNodes[Nodes2.size()] = NC;
535 for(i=0; i<Nodes1.size(); i++)
536 tmpNodes[Nodes2.size()+1+i] = Nodes1[i];
537 for(i=1; i<tmpNodes.size(); i++) {
538 SMDS_MeshEdge* ME = myHelper->AddEdge( tmpNodes[i-1], tmpNodes[i] );
539 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
541 markLinEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
544 else // nbe==3 or ( nbe==2 && linEdge is INTERNAL )
546 if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
549 // one curve must be a part of circle and other curves must be
552 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
553 Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
554 Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 ));
555 if( aLine1.IsNull() || aLine2.IsNull() ) {
556 // other curve not line
557 return error(COMPERR_BAD_SHAPE);
560 bool linEdge1Computed = false;
561 if ( SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1))
562 if( !sm1->IsEmpty() )
563 if( isEdgeCompitaballyMeshed( LinEdge1, aMesh.GetSubMesh(F) ))
564 linEdge1Computed = true;
566 return error("Invalid set of hypotheses");
568 bool linEdge2Computed = false;
569 if ( SMESH_subMesh* sm2 = aMesh.GetSubMesh(LinEdge2))
570 if( !sm2->IsEmpty() )
571 if( isEdgeCompitaballyMeshed( LinEdge2, aMesh.GetSubMesh(F) ))
572 linEdge2Computed = true;
574 return error("Invalid set of hypotheses");
576 bool ok = _gen->Compute( aMesh, CircEdge );
577 if( !ok ) return false;
578 map< double, const SMDS_MeshNode* > theNodes;
579 GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
581 const SMDS_MeshNode* NF = theNodes.begin()->second;
582 const SMDS_MeshNode* NL = theNodes.rbegin()->second;
584 CNodes.push_back( NF );
585 map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
586 double fang = (*itn).first;
588 for(; itn != theNodes.end(); itn++ ) {
589 CNodes.push_back( (*itn).second );
590 double ang = (*itn).first - fang;
591 if( ang>PI ) ang = ang - 2*PI;
592 if( ang<-PI ) ang = ang + 2*PI;
593 Angles.Append( ang );
595 P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
596 gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
597 P0 = aCirc->Location();
599 myLayerPositions.clear();
600 computeLayerPositions(P0,P1);
602 Nodes1.resize( myLayerPositions.size()+1 );
603 Nodes2.resize( myLayerPositions.size()+1 );
605 exp.Init( LinEdge1, TopAbs_VERTEX );
606 TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
608 TopoDS_Vertex V2 = TopoDS::Vertex( exp.Current() );
609 gp_Pnt PE1 = BRep_Tool::Pnt(V1);
610 gp_Pnt PE2 = BRep_Tool::Pnt(V2);
611 if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
612 ( P1.Distance(PE2) > Precision::Confusion() ) )
614 std::swap( LinEdge1, LinEdge2 );
615 std::swap( linEdge1Computed, linEdge2Computed );
617 TopoDS_Vertex VC = V2;
618 if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
619 ( P2.Distance(PE1) > Precision::Confusion() ) )
621 int vertID = meshDS->ShapeToIndex(VC);
624 if ( linEdge1Computed )
626 if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
627 return error("Invalid mesh on a straight edge");
629 bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
630 NC = const_cast<SMDS_MeshNode*>
631 ( nodesFromP0ToP1 ? theNodes.begin()->second : theNodes.rbegin()->second );
632 int i = 0, ir = Nodes1.size()-1;
633 int * pi = nodesFromP0ToP1 ? &i : &ir;
634 itn = theNodes.begin();
635 if ( nodesFromP0ToP1 ) ++itn;
636 for ( ; i < Nodes1.size(); ++i, --ir, ++itn )
638 Nodes1[*pi] = itn->second;
640 for ( i = 0; i < Nodes1.size()-1; ++i )
642 Points.Append( gpXYZ( Nodes1[i]));
643 Pnts2d1.Append( myHelper->GetNodeUV( F, Nodes1[i]));
648 int edgeID = meshDS->ShapeToIndex(LinEdge1);
651 Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
652 gp_Pnt Ptmp = Crv->Value(fp);
654 if( P1.Distance(Ptmp) > Precision::Confusion() )
656 // get UV points for edge
658 BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
661 V2d = gp_Vec2d(PF,PL);
665 V2d = gp_Vec2d(PL,PF);
668 NC = const_cast<SMDS_MeshNode*>( VertexNode( VC, meshDS ));
671 NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
672 meshDS->SetNodeOnVertex(NC, vertID);
676 for(; i<myLayerPositions.size(); i++) {
677 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
678 P0.Y() + aVec.Y()*myLayerPositions[i],
679 P0.Z() + aVec.Z()*myLayerPositions[i] );
681 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
685 param = fp + dp*(1-myLayerPositions[i]);
687 param = fp + dp*myLayerPositions[i];
688 meshDS->SetNodeOnEdge(node, edgeID, param);
689 // parameters on face
690 gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
691 PC.Y() + V2d.Y()*myLayerPositions[i] );
694 Nodes1[ myLayerPositions.size() ] = NF;
695 // create 1D elements on edge
696 SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes1[0] );
697 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
698 for(i=1; i<Nodes1.size(); i++) {
699 ME = myHelper->AddEdge( Nodes1[i-1], Nodes1[i] );
700 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
702 if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
705 markLinEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
708 if ( linEdge2Computed )
710 if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge2,true,theNodes))
711 return error("Invalid mesh on a straight edge");
713 bool nodesFromP0ToP2 = ( theNodes.rbegin()->second == NL );
714 int i = 0, ir = Nodes1.size()-1;
715 int * pi = nodesFromP0ToP2 ? &i : &ir;
716 itn = theNodes.begin();
717 if ( nodesFromP0ToP2 ) ++itn;
718 for ( ; i < Nodes2.size(); ++i, --ir, ++itn )
719 Nodes2[*pi] = itn->second;
723 int edgeID = meshDS->ShapeToIndex(LinEdge2);
724 gp_Vec aVec = gp_Vec(P0,P2);
726 Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge2,fp,lp);
727 gp_Pnt Ptmp = Crv->Value(fp);
729 if( P2.Distance(Ptmp) > Precision::Confusion() )
731 // get UV points for edge
733 BRep_Tool::UVPoints( LinEdge2, TopoDS::Face(aShape), PF, PL );
736 V2d = gp_Vec2d(PF,PL);
740 V2d = gp_Vec2d(PL,PF);
744 for(int i=0; i<myLayerPositions.size(); i++) {
745 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
746 P0.Y() + aVec.Y()*myLayerPositions[i],
747 P0.Z() + aVec.Z()*myLayerPositions[i] );
748 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
752 param = fp + dp*(1-myLayerPositions[i]);
754 param = fp + dp*myLayerPositions[i];
755 meshDS->SetNodeOnEdge(node, edgeID, param);
756 // parameters on face
757 gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
758 PC.Y() + V2d.Y()*myLayerPositions[i] );
760 Nodes2[ myLayerPositions.size() ] = NL;
761 // create 1D elements on edge
762 SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes2[0] );
763 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
764 for(int i=1; i<Nodes2.size(); i++) {
765 ME = myHelper->AddEdge( Nodes2[i-1], Nodes2[i] );
766 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
769 markLinEdgeAsComputedByMe( LinEdge2, aMesh.GetSubMesh( F ));
773 bool IsForward = ( CircEdge.Orientation()==TopAbs_FORWARD );
775 // create nodes and mesh elements on face
776 // find axis of rotation
777 gp_Pnt P2 = gp_Pnt( CNodes[1]->X(), CNodes[1]->Y(), CNodes[1]->Z() );
780 gp_Vec Axis = Vec1.Crossed(Vec2);
783 //cout<<"Angles.Length() = "<<Angles.Length()<<" Points.Length() = "<<Points.Length()<<endl;
784 //cout<<"Nodes1.size() = "<<Nodes1.size()<<" Pnts2d1.Length() = "<<Pnts2d1.Length()<<endl;
785 for(; i<Angles.Length(); i++) {
786 vector< const SMDS_MeshNode* > tmpNodes;
787 tmpNodes.reserve(Nodes1.size());
789 gp_Ax1 theAxis(P0,gp_Dir(Axis));
790 aTrsf.SetRotation( theAxis, Angles.Value(i) );
792 aTrsf2d.SetRotation( PC, Angles.Value(i) );
795 for(; j<=Points.Length(); j++) {
797 Points.Value(j).Coord( cx, cy, cz );
798 aTrsf.Transforms( cx, cy, cz );
799 SMDS_MeshNode* node = myHelper->AddNode( cx, cy, cz );
800 // find parameters on face
801 Pnts2d1.Value(j).Coord( cx, cy );
802 aTrsf2d.Transforms( cx, cy );
804 meshDS->SetNodeOnFace( node, faceID, cx, cy );
805 tmpNodes[j-1] = node;
808 tmpNodes[Points.Length()] = CNodes[i];
810 for(j=0; j<Nodes1.size()-1; j++) {
813 MF = myHelper->AddFace( tmpNodes[j], Nodes1[j],
814 Nodes1[j+1], tmpNodes[j+1] );
816 MF = myHelper->AddFace( tmpNodes[j], tmpNodes[j+1],
817 Nodes1[j+1], Nodes1[j] );
818 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
823 MF = myHelper->AddFace( NC, Nodes1[0], tmpNodes[0] );
825 MF = myHelper->AddFace( NC, tmpNodes[0], Nodes1[0] );
826 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
827 for(j=0; j<Nodes1.size(); j++) {
828 Nodes1[j] = tmpNodes[j];
833 for(i=0; i<Nodes1.size()-1; i++) {
836 MF = myHelper->AddFace( Nodes2[i], Nodes1[i],
837 Nodes1[i+1], Nodes2[i+1] );
839 MF = myHelper->AddFace( Nodes2[i], Nodes2[i+1],
840 Nodes1[i+1], Nodes1[i] );
841 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
846 MF = myHelper->AddFace( NC, Nodes1[0], Nodes2[0] );
848 MF = myHelper->AddFace( NC, Nodes2[0], Nodes1[0] );
849 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
855 //================================================================================
856 //================================================================================
858 * \brief Class computing layers distribution using data of
859 * StdMeshers_LayerDistribution hypothesis
861 //================================================================================
862 //================================================================================
864 class TNodeDistributor: public StdMeshers_Regular_1D
866 list <const SMESHDS_Hypothesis *> myUsedHyps;
868 // -----------------------------------------------------------------------------
869 static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
871 const int myID = -1000;
872 map < int, SMESH_1D_Algo * > & algoMap = aMesh.GetGen()->_map1D_Algo;
873 map < int, SMESH_1D_Algo * >::iterator id_algo = algoMap.find( myID );
874 if ( id_algo == algoMap.end() )
875 return new TNodeDistributor( myID, 0, aMesh.GetGen() );
876 return static_cast< TNodeDistributor* >( id_algo->second );
878 // -----------------------------------------------------------------------------
879 bool Compute( vector< double > & positions,
883 const StdMeshers_LayerDistribution* hyp)
885 double len = pIn.Distance( pOut );
886 if ( len <= DBL_MIN ) return error("Too close points of inner and outer shells");
888 if ( !hyp || !hyp->GetLayerDistribution() )
889 return error( "Invalid LayerDistribution hypothesis");
891 myUsedHyps.push_back( hyp->GetLayerDistribution() );
893 TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut );
894 SMESH_Hypothesis::Hypothesis_Status aStatus;
895 if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
896 return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
897 "with LayerDistribution hypothesis");
899 BRepAdaptor_Curve C3D(edge);
900 double f = C3D.FirstParameter(), l = C3D.LastParameter();
901 list< double > params;
902 if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false ))
903 return error("StdMeshers_Regular_1D failed to compute layers distribution");
906 positions.reserve( params.size() );
907 for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
908 positions.push_back( *itU / len );
912 // -----------------------------------------------------------------------------
913 TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen)
914 : StdMeshers_Regular_1D( hypId, studyId, gen)
917 // -----------------------------------------------------------------------------
918 virtual const list <const SMESHDS_Hypothesis *> &
919 GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
923 // -----------------------------------------------------------------------------
926 //================================================================================
928 * \brief Compute positions of nodes between the internal and the external surfaces
929 * \retval bool - is a success
931 //================================================================================
933 bool StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(const gp_Pnt& pIn,
938 int nbSegments = myNbLayerHypo->GetNumberOfLayers();
939 myLayerPositions.resize( nbSegments - 1 );
940 for ( int z = 1; z < nbSegments; ++z )
941 myLayerPositions[ z - 1 ] = double( z )/ double( nbSegments );
944 if ( myDistributionHypo ) {
945 SMESH_Mesh * mesh = myHelper->GetMesh();
946 if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions, pIn, pOut,
947 *mesh, myDistributionHypo ))
949 error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
953 RETURN_BAD_RESULT("Bad hypothesis");
957 //=======================================================================
958 //function : Evaluate
960 //=======================================================================
962 bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
963 const TopoDS_Shape& aShape,
964 MapShapeNbElems& aResMap)
966 if( aShape.ShapeType() != TopAbs_FACE ) {
969 SMESH_subMesh * smf = aMesh.GetSubMesh(aShape);
970 MapShapeNbElemsItr anIt = aResMap.find(smf);
971 if( anIt != aResMap.end() ) {
975 myLayerPositions.clear();
978 computeLayerPositions(P0,P1);
980 TopoDS_Edge E1,E2,E3;
981 Handle(Geom_Curve) C1,C2,C3;
982 double f1,l1,f2,l2,f3,l3;
985 for ( exp.Init( aShape, TopAbs_EDGE ); exp.More(); exp.Next() ) {
987 TopoDS_Edge E = TopoDS::Edge( exp.Current() );
990 C1 = BRep_Tool::Curve(E,f1,l1);
994 C2 = BRep_Tool::Curve(E,f2,l2);
998 C3 = BRep_Tool::Curve(E,f3,l3);
1002 TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
1003 int nb0d=0, nb2d_tria=0, nb2d_quad=0;
1004 bool isQuadratic = false;
1006 // C1 must be a circle
1007 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
1008 if( !aCirc.IsNull() ) {
1009 bool ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
1011 SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
1012 MapShapeNbElemsItr anIt = aResMap.find(sm);
1013 vector<int> aVec = (*anIt).second;
1014 isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1017 nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1018 // radial medium nodes
1019 nb0d += (aVec[SMDSEntity_Node]+1) * (myLayerPositions.size()+1);
1020 // other medium nodes
1021 nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1024 nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1026 nb2d_tria = aVec[SMDSEntity_Node] + 1;
1032 // one curve must be a half of circle and other curve must be
1033 // a segment of line
1034 Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
1035 while( !tc.IsNull() ) {
1036 C1 = tc->BasisCurve();
1037 tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
1039 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1040 while( !tc.IsNull() ) {
1041 C2 = tc->BasisCurve();
1042 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1044 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
1045 Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast(C2);
1050 if( aCirc.IsNull() ) {
1051 aCirc = Handle(Geom_Circle)::DownCast(C2);
1056 aLine = Handle(Geom_Line)::DownCast(C3);
1058 bool ok = !aCirc.IsNull() && !aLine.IsNull();
1059 if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) {
1060 // not half of circle
1063 SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1);
1064 MapShapeNbElemsItr anIt = aResMap.find(sm1);
1065 if( anIt!=aResMap.end() ) {
1069 ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
1072 SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
1073 MapShapeNbElemsItr anIt = aResMap.find(sm);
1074 vector<int> aVec = (*anIt).second;
1075 isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1078 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1079 // radial medium nodes
1080 nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1081 // other medium nodes
1082 nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1085 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1087 nb2d_tria = aVec[SMDSEntity_Node] + 1;
1088 nb2d_quad = nb2d_tria * myLayerPositions.size();
1089 // add evaluation for edges
1090 vector<int> aResVec(SMDSEntity_Last);
1091 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1093 aResVec[SMDSEntity_Node] = 4*myLayerPositions.size() + 3;
1094 aResVec[SMDSEntity_Quad_Edge] = 2*myLayerPositions.size() + 2;
1097 aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1098 aResVec[SMDSEntity_Edge] = 2*myLayerPositions.size() + 2;
1100 sm = aMesh.GetSubMesh(LinEdge1);
1101 aResMap.insert(make_pair(sm,aResVec));
1105 // one curve must be a part of circle and other curves must be
1107 Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
1108 while( !tc.IsNull() ) {
1109 C1 = tc->BasisCurve();
1110 tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
1112 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1113 while( !tc.IsNull() ) {
1114 C2 = tc->BasisCurve();
1115 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1117 tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
1118 while( !tc.IsNull() ) {
1119 C3 = tc->BasisCurve();
1120 tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
1122 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
1123 Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast(C2);
1124 Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast(C3);
1130 if( aCirc.IsNull() ) {
1131 aCirc = Handle(Geom_Circle)::DownCast(C2);
1137 aLine1 = Handle(Geom_Line)::DownCast(C3);
1138 aLine2 = Handle(Geom_Line)::DownCast(C1);
1139 if( aCirc.IsNull() ) {
1140 aCirc = Handle(Geom_Circle)::DownCast(C3);
1146 aLine1 = Handle(Geom_Line)::DownCast(C1);
1147 aLine2 = Handle(Geom_Line)::DownCast(C2);
1150 bool ok = !aCirc.IsNull() && !aLine1.IsNull() && !aLine1.IsNull();
1151 SMESH_subMesh* sm = aMesh.GetSubMesh(LinEdge1);
1152 MapShapeNbElemsItr anIt = aResMap.find(sm);
1153 if( anIt!=aResMap.end() ) {
1156 sm = aMesh.GetSubMesh(LinEdge2);
1157 anIt = aResMap.find(sm);
1158 if( anIt!=aResMap.end() ) {
1162 ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
1165 SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
1166 MapShapeNbElemsItr anIt = aResMap.find(sm);
1167 vector<int> aVec = (*anIt).second;
1168 isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1171 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1172 // radial medium nodes
1173 nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1174 // other medium nodes
1175 nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1178 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1180 nb2d_tria = aVec[SMDSEntity_Node] + 1;
1181 nb2d_quad = nb2d_tria * myLayerPositions.size();
1182 // add evaluation for edges
1183 vector<int> aResVec(SMDSEntity_Last);
1184 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1186 aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1187 aResVec[SMDSEntity_Quad_Edge] = myLayerPositions.size() + 1;
1190 aResVec[SMDSEntity_Node] = myLayerPositions.size();
1191 aResVec[SMDSEntity_Edge] = myLayerPositions.size() + 1;
1193 sm = aMesh.GetSubMesh(LinEdge1);
1194 aResMap.insert(make_pair(sm,aResVec));
1195 sm = aMesh.GetSubMesh(LinEdge2);
1196 aResMap.insert(make_pair(sm,aResVec));
1200 vector<int> aResVec(SMDSEntity_Last);
1201 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1202 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1204 //cout<<"nb0d = "<<nb0d<<" nb2d_tria = "<<nb2d_tria<<" nb2d_quad = "<<nb2d_quad<<endl;
1208 aResVec[SMDSEntity_Quad_Triangle] = nb2d_tria;
1209 aResVec[SMDSEntity_Quad_Quadrangle] = nb2d_quad;
1212 aResVec[SMDSEntity_Triangle] = nb2d_tria;
1213 aResVec[SMDSEntity_Quadrangle] = nb2d_quad;
1215 aResMap.insert(make_pair(sm,aResVec));
1220 aResMap.insert(make_pair(sm,aResVec));
1221 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1222 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
1223 "Submesh can not be evaluated",this));