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 <BRepTools.hxx>
49 #include <BRep_Tool.hxx>
50 #include <TopExp_Explorer.hxx>
52 //#include <TopoDS_Shell.hxx>
53 //#include <TopoDS_Solid.hxx>
54 //#include <TopTools_MapOfShape.hxx>
56 //#include <gp_Pnt.hxx>
59 #include <Geom_TrimmedCurve.hxx>
60 #include <Geom_Circle.hxx>
61 #include <Geom_Line.hxx>
62 #include <TColgp_SequenceOfPnt.hxx>
63 #include <TColgp_SequenceOfPnt2d.hxx>
64 #include <GeomAPI_ProjectPointOnSurf.hxx>
69 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
70 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
72 //typedef StdMeshers_ProjectionUtils TAssocTool;
75 //=======================================================================
76 //function : StdMeshers_RadialQuadrangle_1D2D
78 //=======================================================================
80 StdMeshers_RadialQuadrangle_1D2D::StdMeshers_RadialQuadrangle_1D2D(int hypId,
83 :SMESH_2D_Algo(hypId, studyId, gen)
85 _name = "RadialQuadrangle_1D2D";
86 _shapeType = (1 << TopAbs_FACE); // 1 bit per shape type
88 _compatibleHypothesis.push_back("LayerDistribution2D");
89 _compatibleHypothesis.push_back("NumberOfLayers2D");
91 myDistributionHypo = 0;
92 _requireDescretBoundary = false;
96 //================================================================================
100 //================================================================================
102 StdMeshers_RadialQuadrangle_1D2D::~StdMeshers_RadialQuadrangle_1D2D()
106 //=======================================================================
107 //function : CheckHypothesis
109 //=======================================================================
111 bool StdMeshers_RadialQuadrangle_1D2D::CheckHypothesis
113 const TopoDS_Shape& aShape,
114 SMESH_Hypothesis::Hypothesis_Status& aStatus)
118 myDistributionHypo = 0;
120 list <const SMESHDS_Hypothesis * >::const_iterator itl;
122 const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
123 if ( hyps.size() == 0 ) {
124 aStatus = SMESH_Hypothesis::HYP_MISSING;
125 return false; // can't work with no hypothesis
128 if ( hyps.size() > 1 ) {
129 aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
133 const SMESHDS_Hypothesis *theHyp = hyps.front();
135 string hypName = theHyp->GetName();
137 if (hypName == "NumberOfLayers2D") {
138 myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
139 aStatus = SMESH_Hypothesis::HYP_OK;
142 if (hypName == "LayerDistribution2D") {
143 myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
144 aStatus = SMESH_Hypothesis::HYP_OK;
147 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
152 //=======================================================================
155 //=======================================================================
157 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
158 const TopoDS_Shape& aShape)
161 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
163 myHelper = new SMESH_MesherHelper( aMesh );
164 myHelper->IsQuadraticSubMesh( aShape );
166 myLayerPositions.clear();
168 TopoDS_Edge E1,E2,E3;
169 Handle(Geom_Curve) C1,C2,C3;
170 double f1,l1,f2,l2,f3,l3;
172 for ( exp.Init( aShape, TopAbs_EDGE ); exp.More(); exp.Next() ) {
174 TopoDS_Edge E = TopoDS::Edge( exp.Current() );
177 C1 = BRep_Tool::Curve(E,f1,l1);
181 C2 = BRep_Tool::Curve(E,f2,l2);
185 C3 = BRep_Tool::Curve(E,f3,l3);
190 return error(COMPERR_BAD_SHAPE);
193 // points for rotation
194 TColgp_SequenceOfPnt Points;
195 // angles for rotation
196 TColStd_SequenceOfReal Angles;
197 // Nodes1 and Nodes2 - nodes along radiuses
198 // CNodes - nodes on circle edge
199 std::vector< const SMDS_MeshNode* > Nodes1, Nodes2, CNodes;
201 // parameters edge nodes on face
202 TColgp_SequenceOfPnt2d Pnts2d1, Pnts2d2;
205 int faceID = meshDS->ShapeToIndex(aShape);
206 TopoDS_Face F = TopoDS::Face(aShape);
207 Handle(Geom_Surface) S = BRep_Tool::Surface(F);
210 bool IsForward = F.Orientation()==TopAbs_FORWARD;
212 //cout<<"RadialQuadrangle_1D2D::Compute nbe = "<<nbe<<endl;
213 TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
215 // C1 must be a circle
216 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
218 return error(COMPERR_BAD_SHAPE);
221 bool ok = _gen->Compute( aMesh, CircEdge, false, MeshDim_1D );
222 if( !ok ) return false;
223 std::map< double, const SMDS_MeshNode* > theNodes;
224 ok = GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
225 if( !ok ) return false;
228 std::map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
229 const SMDS_MeshNode* NF = (*itn).second;
230 CNodes.push_back( (*itn).second );
231 double fang = (*itn).first;
232 if ( itn != theNodes.end() ) {
234 for(; itn != theNodes.end(); itn++ ) {
235 CNodes.push_back( (*itn).second );
236 double ang = (*itn).first - fang;
237 if( ang>PI ) ang = ang - 2*PI;
238 if( ang<-PI ) ang = ang + 2*PI;
239 Angles.Append( ang );
242 P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
243 P0 = aCirc->Location();
245 myLayerPositions.clear();
246 computeLayerPositions(P0,P1);
248 exp.Init( CircEdge, TopAbs_VERTEX );
249 TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
250 gp_Pnt2d p2dV = BRep_Tool::Parameters( V1, TopoDS::Face(aShape) );
252 NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
253 GeomAPI_ProjectPointOnSurf PPS(P0,S);
255 PPS.Parameters(1,U0,V0);
256 meshDS->SetNodeOnFace(NC, faceID, U0, V0);
257 PC = gp_Pnt2d(U0,V0);
260 gp_Vec2d aVec2d(PC,p2dV);
261 Nodes1.resize( myLayerPositions.size()+1 );
262 Nodes2.resize( myLayerPositions.size()+1 );
264 for(; i<myLayerPositions.size(); i++) {
265 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
266 P0.Y() + aVec.Y()*myLayerPositions[i],
267 P0.Z() + aVec.Z()*myLayerPositions[i] );
269 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
272 double U = PC.X() + aVec2d.X()*myLayerPositions[i];
273 double V = PC.Y() + aVec2d.Y()*myLayerPositions[i];
274 meshDS->SetNodeOnFace( node, faceID, U, V );
275 Pnts2d1.Append(gp_Pnt2d(U,V));
276 Pnts2d2.Append(gp_Pnt2d(U,V));
278 Nodes1[Nodes1.size()-1] = NF;
279 Nodes2[Nodes1.size()-1] = NF;
282 // one curve must be a half of circle and other curve must be
284 Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
285 while( !tc.IsNull() ) {
286 C1 = tc->BasisCurve();
287 tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
289 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
290 while( !tc.IsNull() ) {
291 C2 = tc->BasisCurve();
292 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
295 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
296 Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast(C2);
301 if( aCirc.IsNull() ) {
302 aCirc = Handle(Geom_Circle)::DownCast(C2);
307 aLine = Handle(Geom_Line)::DownCast(C3);
309 if( aCirc.IsNull() ) {
311 return error(COMPERR_BAD_SHAPE);
313 if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) {
314 // not half of circle
315 return error(COMPERR_BAD_SHAPE);
317 if( aLine.IsNull() ) {
318 // other curve not line
319 return error(COMPERR_BAD_SHAPE);
321 SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1);
323 SMESHDS_SubMesh* sdssm1 = sm1->GetSubMeshDS();
325 if( sm1->GetSubMeshDS()->NbNodes()>0 ) {
326 SMESH_subMesh* sm = aMesh.GetSubMesh(F);
327 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
328 smError.reset(new SMESH_ComputeError(COMPERR_ALGO_FAILED,
329 "Invalid set of hypothesises",this));
335 bool ok = _gen->Compute( aMesh, CircEdge, false, MeshDim_1D );
336 if( !ok ) return false;
337 std::map< double, const SMDS_MeshNode* > theNodes;
338 GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
341 std::map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
342 const SMDS_MeshNode* NF = (*itn).second;
343 CNodes.push_back( (*itn).second );
344 double fang = (*itn).first;
346 const SMDS_MeshNode* NL;
348 for(; itn != theNodes.end(); itn++ ) {
350 if( nbn == theNodes.size() )
352 CNodes.push_back( (*itn).second );
353 double ang = (*itn).first - fang;
354 if( ang>PI ) ang = ang - 2*PI;
355 if( ang<-PI ) ang = ang + 2*PI;
356 Angles.Append( ang );
358 P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
359 gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
360 P0 = aCirc->Location();
362 myLayerPositions.clear();
363 computeLayerPositions(P0,P1);
366 int edgeID = meshDS->ShapeToIndex(LinEdge1);
368 Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
372 if( P1.Distance(Ptmp) > Precision::Confusion() )
374 // get UV points for edge
376 BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
377 PC = gp_Pnt2d( (PF.X()+PL.X())/2, (PF.Y()+PL.Y())/2 );
379 if(ori) V2d = gp_Vec2d(PC,PF);
380 else V2d = gp_Vec2d(PC,PL);
382 double cp = (fp+lp)/2;
383 double dp2 = (lp-fp)/2;
384 NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
385 meshDS->SetNodeOnEdge(NC, edgeID, cp);
386 Nodes1.resize( myLayerPositions.size()+1 );
387 Nodes2.resize( myLayerPositions.size()+1 );
389 for(; i<myLayerPositions.size(); i++) {
390 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
391 P0.Y() + aVec.Y()*myLayerPositions[i],
392 P0.Z() + aVec.Z()*myLayerPositions[i] );
394 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
398 param = fp + dp2*(1-myLayerPositions[i]);
400 param = cp + dp2*myLayerPositions[i];
401 meshDS->SetNodeOnEdge(node, edgeID, param);
402 P = gp_Pnt( P0.X() - aVec.X()*myLayerPositions[i],
403 P0.Y() - aVec.Y()*myLayerPositions[i],
404 P0.Z() - aVec.Z()*myLayerPositions[i] );
405 node = meshDS->AddNode(P.X(), P.Y(), P.Z());
408 param = fp + dp2*(1-myLayerPositions[i]);
410 param = cp + dp2*myLayerPositions[i];
411 meshDS->SetNodeOnEdge(node, edgeID, param);
412 // parameters on face
413 gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
414 PC.Y() + V2d.Y()*myLayerPositions[i] );
416 P2d = gp_Pnt2d( PC.X() - V2d.X()*myLayerPositions[i],
417 PC.Y() - V2d.Y()*myLayerPositions[i] );
420 Nodes1[ myLayerPositions.size() ] = NF;
421 Nodes2[ myLayerPositions.size() ] = NL;
422 // create 1D elements on edge
423 std::vector< const SMDS_MeshNode* > tmpNodes;
424 tmpNodes.resize(2*Nodes1.size()+1);
425 for(i=0; i<Nodes2.size(); i++)
426 tmpNodes[Nodes2.size()-i-1] = Nodes2[i];
427 tmpNodes[Nodes2.size()] = NC;
428 for(i=0; i<Nodes1.size(); i++)
429 tmpNodes[Nodes2.size()+1+i] = Nodes1[i];
430 for(i=1; i<tmpNodes.size(); i++) {
431 SMDS_MeshEdge* ME = myHelper->AddEdge( tmpNodes[i-1], tmpNodes[i] );
432 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
436 // one curve must be a part of circle and other curves must be
438 Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
439 while( !tc.IsNull() ) {
440 C1 = tc->BasisCurve();
441 tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
443 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
444 while( !tc.IsNull() ) {
445 C2 = tc->BasisCurve();
446 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
448 tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
449 while( !tc.IsNull() ) {
450 C3 = tc->BasisCurve();
451 tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
454 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
455 Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast(C2);
456 Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast(C3);
462 if( aCirc.IsNull() ) {
463 aCirc = Handle(Geom_Circle)::DownCast(C2);
469 aLine1 = Handle(Geom_Line)::DownCast(C3);
470 aLine2 = Handle(Geom_Line)::DownCast(C1);
471 if( aCirc.IsNull() ) {
472 aCirc = Handle(Geom_Circle)::DownCast(C3);
478 aLine1 = Handle(Geom_Line)::DownCast(C1);
479 aLine2 = Handle(Geom_Line)::DownCast(C2);
482 if( aCirc.IsNull() ) {
484 return error(COMPERR_BAD_SHAPE);
486 if( aLine1.IsNull() || aLine2.IsNull() ) {
487 // other curve not line
488 return error(COMPERR_BAD_SHAPE);
490 SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1);
491 SMESH_subMesh* sm2 = aMesh.GetSubMesh(LinEdge2);
493 SMESHDS_SubMesh* sdssm1 = sm1->GetSubMeshDS();
494 SMESHDS_SubMesh* sdssm2 = sm2->GetSubMeshDS();
495 if( sdssm1 && sdssm2 ) {
496 if( sm1->GetSubMeshDS()->NbNodes()>0 || sm2->GetSubMeshDS()->NbNodes()>0 ) {
497 SMESH_subMesh* sm = aMesh.GetSubMesh(F);
498 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
499 smError.reset(new SMESH_ComputeError(COMPERR_ALGO_FAILED,
500 "Invalid set of hypothesises",this));
506 bool ok = _gen->Compute( aMesh, CircEdge, false, MeshDim_1D );
507 if( !ok ) return false;
508 std::map< double, const SMDS_MeshNode* > theNodes;
509 GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
512 std::map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
513 const SMDS_MeshNode* NF = (*itn).second;
514 CNodes.push_back( (*itn).second );
515 double fang = (*itn).first;
517 const SMDS_MeshNode* NL;
519 for(; itn != theNodes.end(); itn++ ) {
521 if( nbn == theNodes.size() )
523 CNodes.push_back( (*itn).second );
524 double ang = (*itn).first - fang;
525 if( ang>PI ) ang = ang - 2*PI;
526 if( ang<-PI ) ang = ang + 2*PI;
527 Angles.Append( ang );
529 P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
530 gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
531 P0 = aCirc->Location();
533 myLayerPositions.clear();
534 computeLayerPositions(P0,P1);
536 exp.Init( LinEdge1, TopAbs_VERTEX );
537 TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
539 TopoDS_Vertex V2 = TopoDS::Vertex( exp.Current() );
540 gp_Pnt PE1 = BRep_Tool::Pnt(V1);
541 gp_Pnt PE2 = BRep_Tool::Pnt(V2);
542 if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
543 ( P1.Distance(PE2) > Precision::Confusion() ) ) {
544 TopoDS_Edge E = LinEdge1;
549 if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
550 ( P2.Distance(PE1) > Precision::Confusion() ) ) {
554 int vertID = meshDS->ShapeToIndex(VC);
556 int edgeID = meshDS->ShapeToIndex(LinEdge1);
559 Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
563 if( P1.Distance(Ptmp) > Precision::Confusion() )
565 // get UV points for edge
567 BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
570 V2d = gp_Vec2d(PF,PL);
574 V2d = gp_Vec2d(PL,PF);
577 NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
578 meshDS->SetNodeOnVertex(NC, vertID);
580 Nodes1.resize( myLayerPositions.size()+1 );
582 for(; i<myLayerPositions.size(); i++) {
583 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
584 P0.Y() + aVec.Y()*myLayerPositions[i],
585 P0.Z() + aVec.Z()*myLayerPositions[i] );
587 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
591 param = fp + dp*(1-myLayerPositions[i]);
593 param = fp + dp*myLayerPositions[i];
594 meshDS->SetNodeOnEdge(node, edgeID, param);
595 // parameters on face
596 gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
597 PC.Y() + V2d.Y()*myLayerPositions[i] );
600 Nodes1[ myLayerPositions.size() ] = NF;
601 // create 1D elements on edge
602 SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes1[0] );
603 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
604 for(i=1; i<Nodes1.size(); i++) {
605 SMDS_MeshEdge* ME = myHelper->AddEdge( Nodes1[i-1], Nodes1[i] );
606 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
609 edgeID = meshDS->ShapeToIndex(LinEdge2);
610 aVec = gp_Vec(P0,P2);
612 Crv = BRep_Tool::Curve(LinEdge2,fp,lp);
615 if( P2.Distance(Ptmp) > Precision::Confusion() )
617 // get UV points for edge
618 BRep_Tool::UVPoints( LinEdge2, TopoDS::Face(aShape), PF, PL );
620 V2d = gp_Vec2d(PF,PL);
624 V2d = gp_Vec2d(PL,PF);
628 Nodes2.resize( myLayerPositions.size()+1 );
629 for(i=0; i<myLayerPositions.size(); i++) {
630 gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
631 P0.Y() + aVec.Y()*myLayerPositions[i],
632 P0.Z() + aVec.Z()*myLayerPositions[i] );
633 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
637 param = fp + dp*(1-myLayerPositions[i]);
639 param = fp + dp*myLayerPositions[i];
640 meshDS->SetNodeOnEdge(node, edgeID, param);
641 // parameters on face
642 gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
643 PC.Y() + V2d.Y()*myLayerPositions[i] );
646 Nodes2[ myLayerPositions.size() ] = NL;
647 // create 1D elements on edge
648 ME = myHelper->AddEdge( NC, Nodes2[0] );
649 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
650 for(i=1; i<Nodes2.size(); i++) {
651 SMDS_MeshEdge* ME = myHelper->AddEdge( Nodes2[i-1], Nodes2[i] );
652 if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
656 // create nodes and mesh elements on face
657 // find axis of rotation
658 gp_Pnt P2 = gp_Pnt( CNodes[1]->X(), CNodes[1]->Y(), CNodes[1]->Z() );
661 gp_Vec Axis = Vec1.Crossed(Vec2);
664 //cout<<"Angles.Length() = "<<Angles.Length()<<" Points.Length() = "<<Points.Length()<<endl;
665 //cout<<"Nodes1.size() = "<<Nodes1.size()<<" Pnts2d1.Length() = "<<Pnts2d1.Length()<<endl;
666 for(; i<Angles.Length(); i++) {
667 std::vector< const SMDS_MeshNode* > tmpNodes;
668 tmpNodes.reserve(Nodes1.size());
670 gp_Ax1 theAxis(P0,gp_Dir(Axis));
671 aTrsf.SetRotation( theAxis, Angles.Value(i) );
673 aTrsf2d.SetRotation( PC, Angles.Value(i) );
676 for(; j<=Points.Length(); j++) {
678 Points.Value(j).Coord( cx, cy, cz );
679 aTrsf.Transforms( cx, cy, cz );
680 SMDS_MeshNode* node = myHelper->AddNode( cx, cy, cz );
681 // find parameters on face
682 Pnts2d1.Value(j).Coord( cx, cy );
683 aTrsf2d.Transforms( cx, cy );
685 meshDS->SetNodeOnFace( node, faceID, cx, cy );
686 tmpNodes[j-1] = node;
689 tmpNodes[Points.Length()] = CNodes[i];
691 for(j=0; j<Nodes1.size()-1; j++) {
694 MF = myHelper->AddFace( tmpNodes[j], Nodes1[j],
695 Nodes1[j+1], tmpNodes[j+1] );
697 MF = myHelper->AddFace( tmpNodes[j], tmpNodes[j+1],
698 Nodes1[j+1], Nodes1[j] );
699 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
704 MF = myHelper->AddFace( NC, Nodes1[0], tmpNodes[0] );
706 MF = myHelper->AddFace( NC, tmpNodes[0], Nodes1[0] );
707 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
708 for(j=0; j<Nodes1.size(); j++) {
709 Nodes1[j] = tmpNodes[j];
714 for(i=0; i<Nodes1.size()-1; i++) {
717 MF = myHelper->AddFace( Nodes2[i], Nodes1[i],
718 Nodes1[i+1], Nodes2[i+1] );
720 MF = myHelper->AddFace( Nodes2[i], Nodes2[i+1],
721 Nodes1[i+1], Nodes1[i] );
722 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
727 MF = myHelper->AddFace( NC, Nodes1[0], Nodes2[0] );
729 MF = myHelper->AddFace( NC, Nodes2[0], Nodes1[0] );
730 if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
733 // to delete helper at exit from Compute()
734 std::auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
740 //================================================================================
741 //================================================================================
743 * \brief Class computing layers distribution using data of
744 * StdMeshers_LayerDistribution hypothesis
746 //================================================================================
747 //================================================================================
749 class TNodeDistributor: public StdMeshers_Regular_1D
751 list <const SMESHDS_Hypothesis *> myUsedHyps;
753 // -----------------------------------------------------------------------------
754 static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
756 const int myID = -1000;
757 map < int, SMESH_1D_Algo * > & algoMap = aMesh.GetGen()->_map1D_Algo;
758 map < int, SMESH_1D_Algo * >::iterator id_algo = algoMap.find( myID );
759 if ( id_algo == algoMap.end() )
760 return new TNodeDistributor( myID, 0, aMesh.GetGen() );
761 return static_cast< TNodeDistributor* >( id_algo->second );
763 // -----------------------------------------------------------------------------
764 bool Compute( vector< double > & positions,
768 const StdMeshers_LayerDistribution* hyp)
770 double len = pIn.Distance( pOut );
771 if ( len <= DBL_MIN ) return error("Too close points of inner and outer shells");
773 if ( !hyp || !hyp->GetLayerDistribution() )
774 return error( "Invalid LayerDistribution hypothesis");
776 myUsedHyps.push_back( hyp->GetLayerDistribution() );
778 TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut );
779 SMESH_Hypothesis::Hypothesis_Status aStatus;
780 if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
781 return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
782 "with LayerDistribution hypothesis");
784 BRepAdaptor_Curve C3D(edge);
785 double f = C3D.FirstParameter(), l = C3D.LastParameter();
786 list< double > params;
787 if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false ))
788 return error("StdMeshers_Regular_1D failed to compute layers distribution");
791 positions.reserve( params.size() );
792 for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
793 positions.push_back( *itU / len );
797 // -----------------------------------------------------------------------------
798 TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen)
799 : StdMeshers_Regular_1D( hypId, studyId, gen)
802 // -----------------------------------------------------------------------------
803 virtual const list <const SMESHDS_Hypothesis *> &
804 GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
808 // -----------------------------------------------------------------------------
811 //================================================================================
813 * \brief Compute positions of nodes between the internal and the external surfaces
814 * \retval bool - is a success
816 //================================================================================
818 bool StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(const gp_Pnt& pIn,
823 int nbSegments = myNbLayerHypo->GetNumberOfLayers();
824 myLayerPositions.resize( nbSegments - 1 );
825 for ( int z = 1; z < nbSegments; ++z )
826 myLayerPositions[ z - 1 ] = double( z )/ double( nbSegments );
829 if ( myDistributionHypo ) {
830 SMESH_Mesh * mesh = myHelper->GetMesh();
831 if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions, pIn, pOut,
832 *mesh, myDistributionHypo ))
834 error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
838 RETURN_BAD_RESULT("Bad hypothesis");
842 //=======================================================================
843 //function : Evaluate
845 //=======================================================================
847 bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
848 const TopoDS_Shape& aShape,
849 MapShapeNbElems& aResMap)
851 if( aShape.ShapeType() != TopAbs_FACE ) {
854 SMESH_subMesh * smf = aMesh.GetSubMesh(aShape);
855 MapShapeNbElemsItr anIt = aResMap.find(smf);
856 if( anIt != aResMap.end() ) {
860 myLayerPositions.clear();
863 computeLayerPositions(P0,P1);
865 TopoDS_Edge E1,E2,E3;
866 Handle(Geom_Curve) C1,C2,C3;
867 double f1,l1,f2,l2,f3,l3;
870 for ( exp.Init( aShape, TopAbs_EDGE ); exp.More(); exp.Next() ) {
872 TopoDS_Edge E = TopoDS::Edge( exp.Current() );
875 C1 = BRep_Tool::Curve(E,f1,l1);
879 C2 = BRep_Tool::Curve(E,f2,l2);
883 C3 = BRep_Tool::Curve(E,f3,l3);
887 TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
888 int nb0d=0, nb2d_tria=0, nb2d_quad=0;
889 bool isQuadratic = false;
891 // C1 must be a circle
892 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
893 if( !aCirc.IsNull() ) {
894 bool ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
896 SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
897 MapShapeNbElemsItr anIt = aResMap.find(sm);
898 std::vector<int> aVec = (*anIt).second;
899 isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
902 nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
903 // radial medium nodes
904 nb0d += (aVec[SMDSEntity_Node]+1) * (myLayerPositions.size()+1);
905 // other medium nodes
906 nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
909 nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
911 nb2d_tria = aVec[SMDSEntity_Node] + 1;
917 // one curve must be a half of circle and other curve must be
919 Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
920 while( !tc.IsNull() ) {
921 C1 = tc->BasisCurve();
922 tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
924 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
925 while( !tc.IsNull() ) {
926 C2 = tc->BasisCurve();
927 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
929 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
930 Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast(C2);
935 if( aCirc.IsNull() ) {
936 aCirc = Handle(Geom_Circle)::DownCast(C2);
941 aLine = Handle(Geom_Line)::DownCast(C3);
943 bool ok = !aCirc.IsNull() && !aLine.IsNull();
944 if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) {
945 // not half of circle
948 SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1);
949 MapShapeNbElemsItr anIt = aResMap.find(sm1);
950 if( anIt!=aResMap.end() ) {
954 ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
957 SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
958 MapShapeNbElemsItr anIt = aResMap.find(sm);
959 std::vector<int> aVec = (*anIt).second;
960 isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
963 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
964 // radial medium nodes
965 nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
966 // other medium nodes
967 nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
970 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
972 nb2d_tria = aVec[SMDSEntity_Node] + 1;
973 nb2d_quad = nb2d_tria * myLayerPositions.size();
974 // add evaluation for edges
975 std::vector<int> aResVec(SMDSEntity_Last);
976 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
978 aResVec[SMDSEntity_Node] = 4*myLayerPositions.size() + 3;
979 aResVec[SMDSEntity_Quad_Edge] = 2*myLayerPositions.size() + 2;
982 aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
983 aResVec[SMDSEntity_Edge] = 2*myLayerPositions.size() + 2;
985 sm = aMesh.GetSubMesh(LinEdge1);
986 aResMap.insert(std::make_pair(sm,aResVec));
990 // one curve must be a part of circle and other curves must be
992 Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
993 while( !tc.IsNull() ) {
994 C1 = tc->BasisCurve();
995 tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
997 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
998 while( !tc.IsNull() ) {
999 C2 = tc->BasisCurve();
1000 tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1002 tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
1003 while( !tc.IsNull() ) {
1004 C3 = tc->BasisCurve();
1005 tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
1007 Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
1008 Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast(C2);
1009 Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast(C3);
1015 if( aCirc.IsNull() ) {
1016 aCirc = Handle(Geom_Circle)::DownCast(C2);
1022 aLine1 = Handle(Geom_Line)::DownCast(C3);
1023 aLine2 = Handle(Geom_Line)::DownCast(C1);
1024 if( aCirc.IsNull() ) {
1025 aCirc = Handle(Geom_Circle)::DownCast(C3);
1031 aLine1 = Handle(Geom_Line)::DownCast(C1);
1032 aLine2 = Handle(Geom_Line)::DownCast(C2);
1035 bool ok = !aCirc.IsNull() && !aLine1.IsNull() && !aLine1.IsNull();
1036 SMESH_subMesh* sm = aMesh.GetSubMesh(LinEdge1);
1037 MapShapeNbElemsItr anIt = aResMap.find(sm);
1038 if( anIt!=aResMap.end() ) {
1041 sm = aMesh.GetSubMesh(LinEdge2);
1042 anIt = aResMap.find(sm);
1043 if( anIt!=aResMap.end() ) {
1047 ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
1050 SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
1051 MapShapeNbElemsItr anIt = aResMap.find(sm);
1052 std::vector<int> aVec = (*anIt).second;
1053 isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1056 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1057 // radial medium nodes
1058 nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1059 // other medium nodes
1060 nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1063 nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1065 nb2d_tria = aVec[SMDSEntity_Node] + 1;
1066 nb2d_quad = nb2d_tria * myLayerPositions.size();
1067 // add evaluation for edges
1068 std::vector<int> aResVec(SMDSEntity_Last);
1069 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1071 aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1072 aResVec[SMDSEntity_Quad_Edge] = myLayerPositions.size() + 1;
1075 aResVec[SMDSEntity_Node] = myLayerPositions.size();
1076 aResVec[SMDSEntity_Edge] = myLayerPositions.size() + 1;
1078 sm = aMesh.GetSubMesh(LinEdge1);
1079 aResMap.insert(std::make_pair(sm,aResVec));
1080 sm = aMesh.GetSubMesh(LinEdge2);
1081 aResMap.insert(std::make_pair(sm,aResVec));
1085 std::vector<int> aResVec(SMDSEntity_Last);
1086 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1087 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1089 //cout<<"nb0d = "<<nb0d<<" nb2d_tria = "<<nb2d_tria<<" nb2d_quad = "<<nb2d_quad<<endl;
1093 aResVec[SMDSEntity_Quad_Triangle] = nb2d_tria;
1094 aResVec[SMDSEntity_Quad_Quadrangle] = nb2d_quad;
1097 aResVec[SMDSEntity_Triangle] = nb2d_tria;
1098 aResVec[SMDSEntity_Quadrangle] = nb2d_quad;
1100 aResMap.insert(std::make_pair(sm,aResVec));
1105 aResMap.insert(std::make_pair(sm,aResVec));
1106 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1107 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
1108 "Submesh can not be evaluated",this));