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;
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_MISSING;
117 return false; // can't 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 TLinEdgeMarker : public SMESH_subMeshEventListener
151 TLinEdgeMarker(): SMESH_subMeshEventListener(/*isDeletable=*/false) {}
153 static SMESH_subMeshEventListener* getListener()
155 static TLinEdgeMarker theEdgeMarker;
156 return &theEdgeMarker;
160 // ------------------------------------------------------------------------------
162 * \brief Mark an edge as computed by StdMeshers_RadialQuadrangle_1D2D
164 void markLinEdgeAsComputedByMe(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
166 if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
168 if ( !edgeSM->GetEventListenerData( TLinEdgeMarker::getListener() ))
169 faceSubMesh->SetEventListener( TLinEdgeMarker::getListener(),
170 SMESH_subMeshEventListenerData::MakeData(faceSubMesh),
174 // ------------------------------------------------------------------------------
176 * \brief Return true if a radial edge was meshed with StdMeshers_RadialQuadrangle_1D2D with
177 * the same radial distribution
179 bool isEdgeCompitaballyMeshed(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
181 if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
183 if ( SMESH_subMeshEventListenerData* otherFaceData =
184 edgeSM->GetEventListenerData( TLinEdgeMarker::getListener() ))
186 // compare hypothesis aplied to two disk faces sharing radial edges
187 SMESH_Mesh& mesh = *faceSubMesh->GetFather();
188 SMESH_Algo* radialQuadAlgo = mesh.GetGen()->GetAlgo(mesh, faceSubMesh->GetSubShape() );
189 SMESH_subMesh* otherFaceSubMesh = otherFaceData->mySubMeshes.front();
190 const list <const SMESHDS_Hypothesis *> & hyps1 =
191 radialQuadAlgo->GetUsedHypothesis( mesh, faceSubMesh->GetSubShape());
192 const list <const SMESHDS_Hypothesis *> & hyps2 =
193 radialQuadAlgo->GetUsedHypothesis( mesh, otherFaceSubMesh->GetSubShape());
194 if( hyps1.empty() && hyps2.empty() )
195 return true; // defaul hyps
196 if ( hyps1.size() != hyps2.size() ||
197 strcmp( hyps1.front()->GetName(), hyps2.front()->GetName() ))
199 ostringstream hypDump1, hypDump2;
200 list <const SMESHDS_Hypothesis*>::const_iterator hyp1 = hyps1.begin();
201 for ( ; hyp1 != hyps1.end(); ++hyp1 )
202 const_cast<SMESHDS_Hypothesis*>(*hyp1)->SaveTo( hypDump1 );
203 list <const SMESHDS_Hypothesis*>::const_iterator hyp2 = hyps2.begin();
204 for ( ; hyp2 != hyps2.end(); ++hyp2 )
205 const_cast<SMESHDS_Hypothesis*>(*hyp2)->SaveTo( hypDump2 );
206 return hypDump1.str() == hypDump2.str();
212 //================================================================================
214 * \brief Return base curve of the edge and extremum parameters
216 //================================================================================
218 Handle(Geom_Curve) getCurve(const TopoDS_Edge& edge, double* f=0, double* l=0)
220 Handle(Geom_Curve) C;
221 if ( !edge.IsNull() )
223 double first = 0., last = 0.;
224 C = BRep_Tool::Curve(edge, first, last);
227 Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C);
228 while( !tc.IsNull() ) {
229 C = tc->BasisCurve();
230 tc = Handle(Geom_TrimmedCurve)::DownCast(C);
239 //================================================================================
241 * \brief Return edges of the face
242 * \retval int - nb of edges
244 //================================================================================
246 int analyseFace(const TopoDS_Shape& face,
247 TopoDS_Edge& CircEdge,
248 TopoDS_Edge& LinEdge1,
249 TopoDS_Edge& LinEdge2)
251 CircEdge.Nullify(); LinEdge1.Nullify(); LinEdge2.Nullify();
254 for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next(), ++nbe )
256 const TopoDS_Edge& E = TopoDS::Edge( exp.Current() );
258 Handle(Geom_Curve) C = getCurve(E,&f,&l);
261 if ( C->IsKind( STANDARD_TYPE(Geom_Circle)))
263 if ( CircEdge.IsNull() )
268 else if ( LinEdge1.IsNull() )
278 //=======================================================================
280 * \brief Allow algo to do something after persistent restoration
281 * \param subMesh - restored submesh
283 * call markLinEdgeAsComputedByMe()
285 //=======================================================================
287 void StdMeshers_RadialQuadrangle_1D2D::SubmeshRestored(SMESH_subMesh* faceSubMesh)
289 if ( !faceSubMesh->IsEmpty() )
291 TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
292 analyseFace( faceSubMesh->GetSubShape(), CircEdge, LinEdge1, LinEdge2 );
293 if ( !LinEdge1.IsNull() ) markLinEdgeAsComputedByMe( LinEdge1, faceSubMesh );
294 if ( !LinEdge2.IsNull() ) markLinEdgeAsComputedByMe( LinEdge2, faceSubMesh );
298 //=======================================================================
301 //=======================================================================
303 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
304 const TopoDS_Shape& aShape)
307 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
309 myHelper = new SMESH_MesherHelper( aMesh );
310 myHelper->IsQuadraticSubMesh( aShape );
311 // to delete helper at exit from Compute()
312 auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
314 myLayerPositions.clear();
316 TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
317 int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 );
318 if( nbe>3 || nbe < 1 || CircEdge.IsNull() )
319 return error(COMPERR_BAD_SHAPE);
322 // points for rotation
323 TColgp_SequenceOfPnt Points;
324 // angles for rotation
325 TColStd_SequenceOfReal Angles;
326 // Nodes1 and Nodes2 - nodes along radiuses
327 // CNodes - nodes on circle edge
328 vector< const SMDS_MeshNode* > Nodes1, Nodes2, CNodes;
330 // parameters edge nodes on face
331 TColgp_SequenceOfPnt2d Pnts2d1;
334 int faceID = meshDS->ShapeToIndex(aShape);
335 TopoDS_Face F = TopoDS::Face(aShape);
336 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
340 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
342 bool ok = _gen->Compute( aMesh, CircEdge );
343 if( !ok ) return false;
344 map< double, const SMDS_MeshNode* > theNodes;
345 ok = GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
346 if( !ok ) return false;
349 map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
350 const SMDS_MeshNode* NF = (*itn).second;
351 CNodes.push_back( (*itn).second );
352 double fang = (*itn).first;
353 if ( itn != theNodes.end() ) {
355 for(; itn != theNodes.end(); itn++ ) {
356 CNodes.push_back( (*itn).second );
357 double ang = (*itn).first - fang;
358 if( ang>PI ) ang = ang - 2*PI;
359 if( ang<-PI ) ang = ang + 2*PI;
360 Angles.Append( ang );
363 P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
364 P0 = aCirc->Location();
366 myLayerPositions.clear();
367 computeLayerPositions(P0,P1);
369 exp.Init( CircEdge, TopAbs_VERTEX );
370 TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
371 gp_Pnt2d p2dV = BRep_Tool::Parameters( V1, TopoDS::Face(aShape) );
373 NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
374 GeomAPI_ProjectPointOnSurf PPS(P0,S);
376 PPS.Parameters(1,U0,V0);
377 meshDS->SetNodeOnFace(NC, faceID, U0, V0);
378 PC = gp_Pnt2d(U0,V0);
381 gp_Vec2d aVec2d(PC,p2dV);
382 Nodes1.resize( myLayerPositions.size()+1 );
383 Nodes2.resize( myLayerPositions.size()+1 );
385 for(; i<myLayerPositions.size(); i++) {
386 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
387 P0.Y() + aVec.Y()*myLayerPositions[i],
388 P0.Z() + aVec.Z()*myLayerPositions[i] );
390 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
393 double U = PC.X() + aVec2d.X()*myLayerPositions[i];
394 double V = PC.Y() + aVec2d.Y()*myLayerPositions[i];
395 meshDS->SetNodeOnFace( node, faceID, U, V );
396 Pnts2d1.Append(gp_Pnt2d(U,V));
398 Nodes1[Nodes1.size()-1] = NF;
399 Nodes2[Nodes1.size()-1] = NF;
401 else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL )
403 // one curve must be a half of circle and other curve must be
406 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
407 if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) {
408 // not half of circle
409 return error(COMPERR_BAD_SHAPE);
411 Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
412 if( aLine.IsNull() ) {
413 // other curve not line
414 return error(COMPERR_BAD_SHAPE);
416 bool linEdgeComputed = false;
417 if( SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1) ) {
418 if( !sm1->IsEmpty() )
419 if( isEdgeCompitaballyMeshed( LinEdge1, aMesh.GetSubMesh(F) ))
420 linEdgeComputed = true;
422 return error("Invalid set of hypotheses");
425 bool ok = _gen->Compute( aMesh, CircEdge );
426 if( !ok ) return false;
427 map< double, const SMDS_MeshNode* > theNodes;
428 GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
431 map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
432 double fang = (*itn).first;
434 for(; itn != theNodes.end(); itn++ ) {
435 CNodes.push_back( (*itn).second );
436 double ang = (*itn).first - fang;
437 if( ang>PI ) ang = ang - 2*PI;
438 if( ang<-PI ) ang = ang + 2*PI;
439 Angles.Append( ang );
441 const SMDS_MeshNode* NF = theNodes.begin()->second;
442 const SMDS_MeshNode* NL = theNodes.rbegin()->second;
443 CNodes.push_back( NF );
444 P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
445 gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
446 P0 = aCirc->Location();
448 myLayerPositions.clear();
449 computeLayerPositions(P0,P1);
451 if ( linEdgeComputed )
453 if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
454 return error("Invalid mesh on a straight edge");
456 vector< const SMDS_MeshNode* > *pNodes1 = &Nodes1, *pNodes2 = &Nodes2;
457 bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
458 if ( !nodesFromP0ToP1 ) std::swap( pNodes1, pNodes2 );
460 map< double, const SMDS_MeshNode* >::reverse_iterator ritn = theNodes.rbegin();
461 itn = theNodes.begin();
462 for ( int i = Nodes1.size()-1; i > -1; ++itn, ++ritn, --i )
464 (*pNodes1)[i] = ritn->second;
465 (*pNodes2)[i] = itn->second;
466 Points.Append( gpXYZ( Nodes1[i]));
467 Pnts2d1.Append( myHelper->GetNodeUV( F, Nodes1[i]));
469 NC = const_cast<SMDS_MeshNode*>( itn->second );
470 Points.Remove( Nodes1.size() );
475 int edgeID = meshDS->ShapeToIndex(LinEdge1);
477 Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
481 if( P1.Distance(Ptmp) > Precision::Confusion() )
483 // get UV points for edge
485 BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
486 PC = gp_Pnt2d( (PF.X()+PL.X())/2, (PF.Y()+PL.Y())/2 );
488 if(ori) V2d = gp_Vec2d(PC,PF);
489 else V2d = gp_Vec2d(PC,PL);
491 double cp = (fp+lp)/2;
492 double dp2 = (lp-fp)/2;
493 NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
494 meshDS->SetNodeOnEdge(NC, edgeID, cp);
495 Nodes1.resize( myLayerPositions.size()+1 );
496 Nodes2.resize( myLayerPositions.size()+1 );
498 for(; i<myLayerPositions.size(); i++) {
499 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
500 P0.Y() + aVec.Y()*myLayerPositions[i],
501 P0.Z() + aVec.Z()*myLayerPositions[i] );
503 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
507 param = fp + dp2*(1-myLayerPositions[i]);
509 param = cp + dp2*myLayerPositions[i];
510 meshDS->SetNodeOnEdge(node, edgeID, param);
511 P = gp_Pnt( P0.X() - aVec.X()*myLayerPositions[i],
512 P0.Y() - aVec.Y()*myLayerPositions[i],
513 P0.Z() - aVec.Z()*myLayerPositions[i] );
514 node = meshDS->AddNode(P.X(), P.Y(), P.Z());
517 param = fp + dp2*(1-myLayerPositions[i]);
519 param = cp + dp2*myLayerPositions[i];
520 meshDS->SetNodeOnEdge(node, edgeID, param);
521 // parameters on face
522 gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
523 PC.Y() + V2d.Y()*myLayerPositions[i] );
526 Nodes1[ myLayerPositions.size() ] = NF;
527 Nodes2[ myLayerPositions.size() ] = NL;
528 // create 1D elements on edge
529 vector< const SMDS_MeshNode* > tmpNodes;
530 tmpNodes.resize(2*Nodes1.size()+1);
531 for(i=0; i<Nodes2.size(); i++)
532 tmpNodes[Nodes2.size()-i-1] = Nodes2[i];
533 tmpNodes[Nodes2.size()] = NC;
534 for(i=0; i<Nodes1.size(); i++)
535 tmpNodes[Nodes2.size()+1+i] = Nodes1[i];
536 for(i=1; i<tmpNodes.size(); i++) {
537 SMDS_MeshEdge* ME = myHelper->AddEdge( tmpNodes[i-1], tmpNodes[i] );
538 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
540 markLinEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
543 else // nbe==3 or ( nbe==2 && linEdge is INTERNAL )
545 if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
548 // one curve must be a part of circle and other curves must be
551 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
552 Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
553 Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 ));
554 if( aLine1.IsNull() || aLine2.IsNull() ) {
555 // other curve not line
556 return error(COMPERR_BAD_SHAPE);
559 bool linEdge1Computed = false;
560 if ( SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1))
561 if( !sm1->IsEmpty() )
562 if( isEdgeCompitaballyMeshed( LinEdge1, aMesh.GetSubMesh(F) ))
563 linEdge1Computed = true;
565 return error("Invalid set of hypotheses");
567 bool linEdge2Computed = false;
568 if ( SMESH_subMesh* sm2 = aMesh.GetSubMesh(LinEdge2))
569 if( !sm2->IsEmpty() )
570 if( isEdgeCompitaballyMeshed( LinEdge2, aMesh.GetSubMesh(F) ))
571 linEdge2Computed = true;
573 return error("Invalid set of hypotheses");
575 bool ok = _gen->Compute( aMesh, CircEdge );
576 if( !ok ) return false;
577 map< double, const SMDS_MeshNode* > theNodes;
578 GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
580 const SMDS_MeshNode* NF = theNodes.begin()->second;
581 const SMDS_MeshNode* NL = theNodes.rbegin()->second;
583 CNodes.push_back( NF );
584 map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
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>PI ) ang = ang - 2*PI;
591 if( ang<-PI ) ang = ang + 2*PI;
592 Angles.Append( ang );
594 P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
595 gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
596 P0 = aCirc->Location();
598 myLayerPositions.clear();
599 computeLayerPositions(P0,P1);
601 Nodes1.resize( myLayerPositions.size()+1 );
602 Nodes2.resize( myLayerPositions.size()+1 );
604 exp.Init( LinEdge1, TopAbs_VERTEX );
605 TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
607 TopoDS_Vertex V2 = TopoDS::Vertex( exp.Current() );
608 gp_Pnt PE1 = BRep_Tool::Pnt(V1);
609 gp_Pnt PE2 = BRep_Tool::Pnt(V2);
610 if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
611 ( P1.Distance(PE2) > Precision::Confusion() ) )
613 std::swap( LinEdge1, LinEdge2 );
614 std::swap( linEdge1Computed, linEdge2Computed );
616 TopoDS_Vertex VC = V2;
617 if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
618 ( P2.Distance(PE1) > Precision::Confusion() ) )
620 int vertID = meshDS->ShapeToIndex(VC);
623 if ( linEdge1Computed )
625 if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
626 return error("Invalid mesh on a straight edge");
628 bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
629 NC = const_cast<SMDS_MeshNode*>
630 ( nodesFromP0ToP1 ? theNodes.begin()->second : theNodes.rbegin()->second );
631 int i = 0, ir = Nodes1.size()-1;
632 int * pi = nodesFromP0ToP1 ? &i : &ir;
633 itn = theNodes.begin();
634 if ( nodesFromP0ToP1 ) ++itn;
635 for ( ; i < Nodes1.size(); ++i, --ir, ++itn )
637 Nodes1[*pi] = itn->second;
639 for ( i = 0; i < Nodes1.size()-1; ++i )
641 Points.Append( gpXYZ( Nodes1[i]));
642 Pnts2d1.Append( myHelper->GetNodeUV( F, Nodes1[i]));
647 int edgeID = meshDS->ShapeToIndex(LinEdge1);
650 Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
651 gp_Pnt Ptmp = Crv->Value(fp);
653 if( P1.Distance(Ptmp) > Precision::Confusion() )
655 // get UV points for edge
657 BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
660 V2d = gp_Vec2d(PF,PL);
664 V2d = gp_Vec2d(PL,PF);
667 NC = const_cast<SMDS_MeshNode*>( VertexNode( VC, meshDS ));
670 NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
671 meshDS->SetNodeOnVertex(NC, vertID);
675 for(; i<myLayerPositions.size(); i++) {
676 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
677 P0.Y() + aVec.Y()*myLayerPositions[i],
678 P0.Z() + aVec.Z()*myLayerPositions[i] );
680 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
684 param = fp + dp*(1-myLayerPositions[i]);
686 param = fp + dp*myLayerPositions[i];
687 meshDS->SetNodeOnEdge(node, edgeID, param);
688 // parameters on face
689 gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
690 PC.Y() + V2d.Y()*myLayerPositions[i] );
693 Nodes1[ myLayerPositions.size() ] = NF;
694 // create 1D elements on edge
695 SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes1[0] );
696 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
697 for(i=1; i<Nodes1.size(); i++) {
698 ME = myHelper->AddEdge( Nodes1[i-1], Nodes1[i] );
699 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
701 if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
704 markLinEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
707 if ( linEdge2Computed )
709 if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge2,true,theNodes))
710 return error("Invalid mesh on a straight edge");
712 bool nodesFromP0ToP2 = ( theNodes.rbegin()->second == NL );
713 int i = 0, ir = Nodes1.size()-1;
714 int * pi = nodesFromP0ToP2 ? &i : &ir;
715 itn = theNodes.begin();
716 if ( nodesFromP0ToP2 ) ++itn;
717 for ( ; i < Nodes2.size(); ++i, --ir, ++itn )
718 Nodes2[*pi] = itn->second;
722 int edgeID = meshDS->ShapeToIndex(LinEdge2);
723 gp_Vec aVec = gp_Vec(P0,P2);
725 Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge2,fp,lp);
726 gp_Pnt Ptmp = Crv->Value(fp);
728 if( P2.Distance(Ptmp) > Precision::Confusion() )
730 // get UV points for edge
732 BRep_Tool::UVPoints( LinEdge2, TopoDS::Face(aShape), PF, PL );
735 V2d = gp_Vec2d(PF,PL);
739 V2d = gp_Vec2d(PL,PF);
743 for(int i=0; i<myLayerPositions.size(); i++) {
744 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
745 P0.Y() + aVec.Y()*myLayerPositions[i],
746 P0.Z() + aVec.Z()*myLayerPositions[i] );
747 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
751 param = fp + dp*(1-myLayerPositions[i]);
753 param = fp + dp*myLayerPositions[i];
754 meshDS->SetNodeOnEdge(node, edgeID, param);
755 // parameters on face
756 gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
757 PC.Y() + V2d.Y()*myLayerPositions[i] );
759 Nodes2[ myLayerPositions.size() ] = NL;
760 // create 1D elements on edge
761 SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes2[0] );
762 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
763 for(int i=1; i<Nodes2.size(); i++) {
764 ME = myHelper->AddEdge( Nodes2[i-1], Nodes2[i] );
765 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
768 markLinEdgeAsComputedByMe( LinEdge2, aMesh.GetSubMesh( F ));
772 bool IsForward = ( CircEdge.Orientation()==TopAbs_FORWARD );
774 // create nodes and mesh elements on face
775 // find axis of rotation
776 gp_Pnt P2 = gp_Pnt( CNodes[1]->X(), CNodes[1]->Y(), CNodes[1]->Z() );
779 gp_Vec Axis = Vec1.Crossed(Vec2);
782 //cout<<"Angles.Length() = "<<Angles.Length()<<" Points.Length() = "<<Points.Length()<<endl;
783 //cout<<"Nodes1.size() = "<<Nodes1.size()<<" Pnts2d1.Length() = "<<Pnts2d1.Length()<<endl;
784 for(; i<Angles.Length(); i++) {
785 vector< const SMDS_MeshNode* > tmpNodes;
786 tmpNodes.reserve(Nodes1.size());
788 gp_Ax1 theAxis(P0,gp_Dir(Axis));
789 aTrsf.SetRotation( theAxis, Angles.Value(i) );
791 aTrsf2d.SetRotation( PC, Angles.Value(i) );
794 for(; j<=Points.Length(); j++) {
796 Points.Value(j).Coord( cx, cy, cz );
797 aTrsf.Transforms( cx, cy, cz );
798 SMDS_MeshNode* node = myHelper->AddNode( cx, cy, cz );
799 // find parameters on face
800 Pnts2d1.Value(j).Coord( cx, cy );
801 aTrsf2d.Transforms( cx, cy );
803 meshDS->SetNodeOnFace( node, faceID, cx, cy );
804 tmpNodes[j-1] = node;
807 tmpNodes[Points.Length()] = CNodes[i];
809 for(j=0; j<Nodes1.size()-1; j++) {
812 MF = myHelper->AddFace( tmpNodes[j], Nodes1[j],
813 Nodes1[j+1], tmpNodes[j+1] );
815 MF = myHelper->AddFace( tmpNodes[j], tmpNodes[j+1],
816 Nodes1[j+1], Nodes1[j] );
817 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
822 MF = myHelper->AddFace( NC, Nodes1[0], tmpNodes[0] );
824 MF = myHelper->AddFace( NC, tmpNodes[0], Nodes1[0] );
825 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
826 for(j=0; j<Nodes1.size(); j++) {
827 Nodes1[j] = tmpNodes[j];
832 for(i=0; i<Nodes1.size()-1; i++) {
835 MF = myHelper->AddFace( Nodes2[i], Nodes1[i],
836 Nodes1[i+1], Nodes2[i+1] );
838 MF = myHelper->AddFace( Nodes2[i], Nodes2[i+1],
839 Nodes1[i+1], Nodes1[i] );
840 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
845 MF = myHelper->AddFace( NC, Nodes1[0], Nodes2[0] );
847 MF = myHelper->AddFace( NC, Nodes2[0], Nodes1[0] );
848 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
854 //================================================================================
855 //================================================================================
857 * \brief Class computing layers distribution using data of
858 * StdMeshers_LayerDistribution hypothesis
860 //================================================================================
861 //================================================================================
863 class TNodeDistributor: public StdMeshers_Regular_1D
865 list <const SMESHDS_Hypothesis *> myUsedHyps;
867 // -----------------------------------------------------------------------------
868 static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
870 const int myID = -1000;
871 map < int, SMESH_1D_Algo * > & algoMap = aMesh.GetGen()->_map1D_Algo;
872 map < int, SMESH_1D_Algo * >::iterator id_algo = algoMap.find( myID );
873 if ( id_algo == algoMap.end() )
874 return new TNodeDistributor( myID, 0, aMesh.GetGen() );
875 return static_cast< TNodeDistributor* >( id_algo->second );
877 // -----------------------------------------------------------------------------
878 bool Compute( vector< double > & positions,
882 const StdMeshers_LayerDistribution* hyp)
884 double len = pIn.Distance( pOut );
885 if ( len <= DBL_MIN ) return error("Too close points of inner and outer shells");
887 if ( !hyp || !hyp->GetLayerDistribution() )
888 return error( "Invalid LayerDistribution hypothesis");
890 myUsedHyps.push_back( hyp->GetLayerDistribution() );
892 TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut );
893 SMESH_Hypothesis::Hypothesis_Status aStatus;
894 if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
895 return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
896 "with LayerDistribution hypothesis");
898 BRepAdaptor_Curve C3D(edge);
899 double f = C3D.FirstParameter(), l = C3D.LastParameter();
900 list< double > params;
901 if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false ))
902 return error("StdMeshers_Regular_1D failed to compute layers distribution");
905 positions.reserve( params.size() );
906 for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
907 positions.push_back( *itU / len );
911 // -----------------------------------------------------------------------------
912 TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen)
913 : StdMeshers_Regular_1D( hypId, studyId, gen)
916 // -----------------------------------------------------------------------------
917 virtual const list <const SMESHDS_Hypothesis *> &
918 GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
922 // -----------------------------------------------------------------------------
925 //================================================================================
927 * \brief Compute positions of nodes between the internal and the external surfaces
928 * \retval bool - is a success
930 //================================================================================
932 bool StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(const gp_Pnt& pIn,
937 int nbSegments = myNbLayerHypo->GetNumberOfLayers();
938 myLayerPositions.resize( nbSegments - 1 );
939 for ( int z = 1; z < nbSegments; ++z )
940 myLayerPositions[ z - 1 ] = double( z )/ double( nbSegments );
943 if ( myDistributionHypo ) {
944 SMESH_Mesh * mesh = myHelper->GetMesh();
945 if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions, pIn, pOut,
946 *mesh, myDistributionHypo ))
948 error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
952 RETURN_BAD_RESULT("Bad hypothesis");
956 //=======================================================================
957 //function : Evaluate
959 //=======================================================================
961 bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
962 const TopoDS_Shape& aShape,
963 MapShapeNbElems& aResMap)
965 if( aShape.ShapeType() != TopAbs_FACE ) {
968 SMESH_subMesh * smf = aMesh.GetSubMesh(aShape);
969 MapShapeNbElemsItr anIt = aResMap.find(smf);
970 if( anIt != aResMap.end() ) {
974 myLayerPositions.clear();
977 computeLayerPositions(P0,P1);
979 TopoDS_Edge E1,E2,E3;
980 Handle(Geom_Curve) C1,C2,C3;
981 double f1,l1,f2,l2,f3,l3;
984 for ( exp.Init( aShape, TopAbs_EDGE ); exp.More(); exp.Next() ) {
986 TopoDS_Edge E = TopoDS::Edge( exp.Current() );
989 C1 = BRep_Tool::Curve(E,f1,l1);
993 C2 = BRep_Tool::Curve(E,f2,l2);
997 C3 = BRep_Tool::Curve(E,f3,l3);
1001 TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
1002 int nb0d=0, nb2d_tria=0, nb2d_quad=0;
1003 bool isQuadratic = false;
1005 // C1 must be a circle
1006 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
1007 if( !aCirc.IsNull() ) {
1008 bool ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
1010 SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
1011 MapShapeNbElemsItr anIt = aResMap.find(sm);
1012 vector<int> aVec = (*anIt).second;
1013 isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1016 nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1017 // radial medium nodes
1018 nb0d += (aVec[SMDSEntity_Node]+1) * (myLayerPositions.size()+1);
1019 // other medium nodes
1020 nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1023 nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1025 nb2d_tria = aVec[SMDSEntity_Node] + 1;
1031 // one curve must be a half of circle and other curve must be
1032 // a segment of line
1033 Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
1034 while( !tc.IsNull() ) {
1035 C1 = tc->BasisCurve();
1036 tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
1038 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1039 while( !tc.IsNull() ) {
1040 C2 = tc->BasisCurve();
1041 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1043 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
1044 Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast(C2);
1049 if( aCirc.IsNull() ) {
1050 aCirc = Handle(Geom_Circle)::DownCast(C2);
1055 aLine = Handle(Geom_Line)::DownCast(C3);
1057 bool ok = !aCirc.IsNull() && !aLine.IsNull();
1058 if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) {
1059 // not half of circle
1062 SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1);
1063 MapShapeNbElemsItr anIt = aResMap.find(sm1);
1064 if( anIt!=aResMap.end() ) {
1068 ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
1071 SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
1072 MapShapeNbElemsItr anIt = aResMap.find(sm);
1073 vector<int> aVec = (*anIt).second;
1074 isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1077 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1078 // radial medium nodes
1079 nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1080 // other medium nodes
1081 nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1084 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1086 nb2d_tria = aVec[SMDSEntity_Node] + 1;
1087 nb2d_quad = nb2d_tria * myLayerPositions.size();
1088 // add evaluation for edges
1089 vector<int> aResVec(SMDSEntity_Last);
1090 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1092 aResVec[SMDSEntity_Node] = 4*myLayerPositions.size() + 3;
1093 aResVec[SMDSEntity_Quad_Edge] = 2*myLayerPositions.size() + 2;
1096 aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1097 aResVec[SMDSEntity_Edge] = 2*myLayerPositions.size() + 2;
1099 sm = aMesh.GetSubMesh(LinEdge1);
1100 aResMap.insert(make_pair(sm,aResVec));
1104 // one curve must be a part of circle and other curves must be
1106 Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
1107 while( !tc.IsNull() ) {
1108 C1 = tc->BasisCurve();
1109 tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
1111 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1112 while( !tc.IsNull() ) {
1113 C2 = tc->BasisCurve();
1114 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1116 tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
1117 while( !tc.IsNull() ) {
1118 C3 = tc->BasisCurve();
1119 tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
1121 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
1122 Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast(C2);
1123 Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast(C3);
1129 if( aCirc.IsNull() ) {
1130 aCirc = Handle(Geom_Circle)::DownCast(C2);
1136 aLine1 = Handle(Geom_Line)::DownCast(C3);
1137 aLine2 = Handle(Geom_Line)::DownCast(C1);
1138 if( aCirc.IsNull() ) {
1139 aCirc = Handle(Geom_Circle)::DownCast(C3);
1145 aLine1 = Handle(Geom_Line)::DownCast(C1);
1146 aLine2 = Handle(Geom_Line)::DownCast(C2);
1149 bool ok = !aCirc.IsNull() && !aLine1.IsNull() && !aLine1.IsNull();
1150 SMESH_subMesh* sm = aMesh.GetSubMesh(LinEdge1);
1151 MapShapeNbElemsItr anIt = aResMap.find(sm);
1152 if( anIt!=aResMap.end() ) {
1155 sm = aMesh.GetSubMesh(LinEdge2);
1156 anIt = aResMap.find(sm);
1157 if( anIt!=aResMap.end() ) {
1161 ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
1164 SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
1165 MapShapeNbElemsItr anIt = aResMap.find(sm);
1166 vector<int> aVec = (*anIt).second;
1167 isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1170 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1171 // radial medium nodes
1172 nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1173 // other medium nodes
1174 nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1177 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1179 nb2d_tria = aVec[SMDSEntity_Node] + 1;
1180 nb2d_quad = nb2d_tria * myLayerPositions.size();
1181 // add evaluation for edges
1182 vector<int> aResVec(SMDSEntity_Last);
1183 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1185 aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1186 aResVec[SMDSEntity_Quad_Edge] = myLayerPositions.size() + 1;
1189 aResVec[SMDSEntity_Node] = myLayerPositions.size();
1190 aResVec[SMDSEntity_Edge] = myLayerPositions.size() + 1;
1192 sm = aMesh.GetSubMesh(LinEdge1);
1193 aResMap.insert(make_pair(sm,aResVec));
1194 sm = aMesh.GetSubMesh(LinEdge2);
1195 aResMap.insert(make_pair(sm,aResVec));
1199 vector<int> aResVec(SMDSEntity_Last);
1200 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1201 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1203 //cout<<"nb0d = "<<nb0d<<" nb2d_tria = "<<nb2d_tria<<" nb2d_quad = "<<nb2d_quad<<endl;
1207 aResVec[SMDSEntity_Quad_Triangle] = nb2d_tria;
1208 aResVec[SMDSEntity_Quad_Quadrangle] = nb2d_quad;
1211 aResVec[SMDSEntity_Triangle] = nb2d_tria;
1212 aResVec[SMDSEntity_Quadrangle] = nb2d_quad;
1214 aResMap.insert(make_pair(sm,aResVec));
1219 aResMap.insert(make_pair(sm,aResVec));
1220 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1221 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
1222 "Submesh can not be evaluated",this));