// File: GeomPlate_BuildPlateSurface.cxx // Created: Mon May 5 17:07:03 1997 // Author: Jerome LEMONIER // // Modification de l'interface // Amelioration de l'aglo de remplissage // 29-09-97 ; jct; correction du calcul des points doublons (tol 2d et non 3d) // 16-10-97 ; jct; on remet un Raise au lieu du cout sous debug (sinon plantage dans k1fab) // 03-11-97 ; jct; plus de reference a BRepAdaptor et BRepTools #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include // pour la verif G2 #include // pour les intersection entre les courbes #include #include #include #include // projection #include #include #include #include #include #include #include #include #include #include // pour mes tests #ifdef DEB #include #endif #ifdef DRAW static Standard_Integer Affich=0; // 0 : Pas de display // 1 : Display des Geometries et controle intermediaire // 2 : Display du nombre de contrainte par courbe + Intersection // 3 : Dump des contraintes dans Plate static Standard_Integer NbPlan = 0; //static Standard_Integer NbCurv2d = 0; static Standard_Integer NbMark = 0; static Standard_Integer NbProj = 0; #endif #ifdef DRAW #include #include #include #include #endif #include #include #include #include #include #include //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// // ========================================================= // C O N S T R U C T E U R S // ========================================================= //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// //--------------------------------------------------------- // Constructeur compatible avec l'ancienne version //--------------------------------------------------------- GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface ( const Handle(TColStd_HArray1OfInteger)& NPoints, const Handle(GeomPlate_HArray1OfHCurveOnSurface)& TabCurve, const Handle(TColStd_HArray1OfInteger)& Tang, const Standard_Integer Degree, const Standard_Integer NbIter, const Standard_Real Tol2d, const Standard_Real Tol3d, const Standard_Real TolAng, const Standard_Real TolCurv, const Standard_Boolean Anisotropie ) : myAnisotropie(Anisotropie), myDegree(Degree), myNbIter(NbIter), myProj(), myTol2d(Tol2d), myTol3d(Tol3d), myTolAng(TolAng), myTolCurv(TolCurv), myNbBounds(0) { Standard_Integer NTCurve=TabCurve->Length();// Nombre de contraintes lineaires myNbPtsOnCur = 0; // Debrayage du calcul du nombre de points // en fonction de la longueur myLinCont = new GeomPlate_HSequenceOfCurveConstraint; myPntCont = new GeomPlate_HSequenceOfPointConstraint; if (myNbIter<1) Standard_ConstructionError::Raise("GeomPlate : Number of iteration must be >= 1"); if (NTCurve==0) Standard_ConstructionError::Raise("GeomPlate : the bounds Array is null"); if (Tang->Length()==0) Standard_ConstructionError::Raise("GeomPlate : the constraints Array is null"); Standard_Integer nbp = 0; Standard_Integer i ; for ( i=1;i<=NTCurve;i++) { nbp+=NPoints->Value(i); } if (nbp==0) Standard_ConstructionError::Raise("GeomPlate : the resolution is impossible if the number of constraints points is 0"); if (myDegree<2) Standard_ConstructionError::Raise("GeomPlate ; the degree resolution must be upper of 2"); // Remplissage des champs passage de l'ancien constructeur au nouveau for(i=1;i<=NTCurve;i++) { Handle(GeomPlate_CurveConstraint) Cont = new GeomPlate_CurveConstraint ( TabCurve->Value(i), Tang->Value(i), NPoints->Value(i)); myLinCont->Append(Cont); } mySurfInitIsGive=Standard_False; myIsLinear = Standard_True; myFree = Standard_False; } //------------------------------------------------------------------ // Constructeur avec surface d'init et degre de resolution de plate //------------------------------------------------------------------ GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface ( const Handle(Geom_Surface)& Surf, const Standard_Integer Degree, const Standard_Integer NbPtsOnCur, const Standard_Integer NbIter, const Standard_Real Tol2d, const Standard_Real Tol3d, const Standard_Real TolAng, const Standard_Real TolCurv, const Standard_Boolean Anisotropie ) : mySurfInit(Surf), myAnisotropie(Anisotropie), myDegree(Degree), myNbPtsOnCur(NbPtsOnCur), myNbIter(NbIter), myProj(), myTol2d(Tol2d), myTol3d(Tol3d), myTolAng(TolAng), myTolCurv(TolCurv), myNbBounds(0) { if (myNbIter<1) Standard_ConstructionError::Raise("GeomPlate : Number of iteration must be >= 1"); if (myDegree<2) Standard_ConstructionError::Raise("GeomPlate : the degree resolution must be upper of 2"); myLinCont = new GeomPlate_HSequenceOfCurveConstraint; myPntCont = new GeomPlate_HSequenceOfPointConstraint; mySurfInitIsGive=Standard_True; myIsLinear = Standard_True; myFree = Standard_False; } //--------------------------------------------------------- // Constructeur avec degre de resolution de plate //--------------------------------------------------------- GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface ( const Standard_Integer Degree, const Standard_Integer NbPtsOnCur, const Standard_Integer NbIter, const Standard_Real Tol2d, const Standard_Real Tol3d, const Standard_Real TolAng, const Standard_Real TolCurv, const Standard_Boolean Anisotropie ) : myAnisotropie(Anisotropie), myDegree(Degree), myNbPtsOnCur(NbPtsOnCur), myNbIter(NbIter), myProj(), myTol2d(Tol2d), myTol3d(Tol3d), myTolAng(TolAng), myTolCurv(TolCurv), myNbBounds(0) { if (myNbIter<1) Standard_ConstructionError::Raise("GeomPlate : Number of iteration must be >= 1"); if (myDegree<2) Standard_ConstructionError::Raise("GeomPlate : the degree resolution must be upper of 2"); myLinCont = new GeomPlate_HSequenceOfCurveConstraint; myPntCont = new GeomPlate_HSequenceOfPointConstraint; mySurfInitIsGive=Standard_False; myIsLinear = Standard_True; myFree = Standard_False; } //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// // ========================================================= // M E T H O D E S P U B L I Q U E S // ========================================================= //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// //------------------------------------------------------------------------- // Fonction : TrierTab, Fonction interne, ne fait partie de la classe //------------------------------------------------------------------------- // Reordonne le tableau des permutations // Apres l'appel a CourbeJointive l'ordre des courbes est modifie // Ex : ordre init des courbe ==> A B C D E F // Dans TabInit on note ==> 1 2 3 4 5 6 // apres CourbeJointive ==> A E C B D F // TabInit ==> 1 5 3 2 4 6 // apres TrierTab le Tableau contient ==> 1 4 3 5 2 6 // On peut ainsi acceder a la 2eme courbe en prenant TabInit[2] // c'est a dire la 4eme du tableau des courbes classees //------------------------------------------------------------------------- void TrierTab(Handle(TColStd_HArray1OfInteger)& Tab) { // Trie le tableau des permutations pour retrouver l'ordre initial Standard_Integer Nb=Tab->Length(); TColStd_Array1OfInteger TabTri(1,Nb); for (Standard_Integer i=1;i<=Nb;i++) TabTri.SetValue(Tab->Value(i),i); Tab->ChangeArray1()=TabTri; } //--------------------------------------------------------- // Fonction : ProjectCurve //--------------------------------------------------------- Handle(Geom2d_Curve) GeomPlate_BuildPlateSurface::ProjectCurve(const Handle(Adaptor3d_HCurve)& Curv) { // Projection d'une courbe sur un plan Handle(Geom2d_Curve) Curve2d ; Handle(GeomAdaptor_HSurface) hsur = new GeomAdaptor_HSurface(mySurfInit); gp_Pnt2d P2d; ProjLib_CompProjectedCurve Projector(hsur, Curv, myTol3d/10, myTol3d/10); Standard_Real Udeb, Ufin, ProjUdeb, ProjUfin; Udeb = Curv->FirstParameter(); Ufin = Curv->LastParameter(); Projector.Bounds( 1, ProjUdeb, ProjUfin ); if (Projector.NbCurves() != 1 || Abs( Udeb-ProjUdeb ) > Precision::PConfusion() || Abs( Ufin-ProjUfin ) > Precision::PConfusion()) { if (Projector.IsSinglePnt(1, P2d)) { //solution ponctuelles TColgp_Array1OfPnt2d poles(1, 2); poles.Init(P2d); Curve2d = new (Geom2d_BezierCurve) (poles); } else { Curve2d.Nullify(); // Pas de solution continue #if PLATE_DEB cout << "BuildPlateSurace :: Pas de projection continue" << endl; #endif } } else { GeomAbs_Shape Continuity = GeomAbs_C1; Standard_Integer MaxDegree = 10, MaxSeg; Standard_Real Udeb, Ufin; Handle(ProjLib_HCompProjectedCurve) HProjector = new ProjLib_HCompProjectedCurve(); HProjector->Set(Projector); Projector.Bounds(1, Udeb, Ufin); MaxSeg = 20 + HProjector->NbIntervals(GeomAbs_C3); Approx_CurveOnSurface appr(HProjector, hsur, Udeb, Ufin, myTol3d, Continuity, MaxDegree, MaxSeg, Standard_False, Standard_True); Curve2d = appr.Curve2d(); } #if DRAW if (Affich) { char* name = new char[100]; sprintf(name,"proj_%d",++NbProj); DrawTrSurf::Set(name, Curve2d); } #endif return Curve2d; } //--------------------------------------------------------- // Fonction : ProjectedCurve //--------------------------------------------------------- Handle(Adaptor2d_HCurve2d) GeomPlate_BuildPlateSurface::ProjectedCurve( Handle(Adaptor3d_HCurve)& Curv) { // Projection d'une courbe sur la surface d'init Handle(GeomAdaptor_HSurface) hsur = new GeomAdaptor_HSurface(mySurfInit); ProjLib_CompProjectedCurve Projector(hsur, Curv, myTolU/10, myTolV/10); Handle(ProjLib_HCompProjectedCurve) HProjector = new ProjLib_HCompProjectedCurve(); if (Projector.NbCurves() != 1) { HProjector.Nullify(); // Pas de solution continue #if PLATE_DEB cout << "BuildPlateSurace :: Pas de projection continue" << endl; #endif } else { Standard_Real First1,Last1,First2,Last2; First1=Curv->FirstParameter(); Last1 =Curv->LastParameter(); Projector.Bounds(1,First2,Last2); if (Abs(First1-First2) <= Max(myTolU,myTolV) && Abs(Last1-Last2) <= Max(myTolU,myTolV)) { HProjector->Set(Projector); HProjector = Handle(ProjLib_HCompProjectedCurve):: DownCast(HProjector->Trim(First2,Last2,Precision::PConfusion())); } else { HProjector.Nullify(); // Pas de solution continue #if PLATE_DEB cout << "BuildPlateSurace :: Pas de projection complete" << endl; #endif } } return HProjector; } //--------------------------------------------------------- // Fonction : ProjectPoint //--------------------------------------------------------- // Projete une point sur la surface d'init //--------------------------------------------------------- gp_Pnt2d GeomPlate_BuildPlateSurface::ProjectPoint(const gp_Pnt &p3d) { Extrema_POnSurf P; myProj.Perform(p3d); Standard_Integer nearest = 1; if( myProj.NbExt() > 1) { Standard_Real dist2mini = myProj.SquareDistance(1); for (Standard_Integer i=2; i<=myProj.NbExt();i++) { if (myProj.SquareDistance(i) < dist2mini) { dist2mini = myProj.SquareDistance(i); nearest = i; } } } P=myProj.Point(nearest); Standard_Real u,v; P.Parameter(u,v); gp_Pnt2d p2d; p2d.SetCoord(u,v); return p2d; } //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// // ========================================================= // M E T H O D E S P U B L I Q U E S // ========================================================= //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// //--------------------------------------------------------- // Fonction : Init //--------------------------------------------------------- // Initialise les contraintes lineaires et ponctuelles //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::Init() { myLinCont->Clear(); myPntCont->Clear(); myPntCont = new GeomPlate_HSequenceOfPointConstraint; myLinCont = new GeomPlate_HSequenceOfCurveConstraint; } //--------------------------------------------------------- // Fonction : LoadInitSurface //--------------------------------------------------------- // Charge la surface initiale //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::LoadInitSurface(const Handle(Geom_Surface)& Surf) { mySurfInit = Surf; mySurfInitIsGive=Standard_True; } //--------------------------------------------------------- //fonction : Add //--------------------------------------------------------- //--------------------------------------------------------- // Ajout d'une contrainte lineaire //--------------------------------------------------------- void GeomPlate_BuildPlateSurface:: Add(const Handle(GeomPlate_CurveConstraint)& Cont) { myLinCont->Append(Cont); } void GeomPlate_BuildPlateSurface::SetNbBounds( const Standard_Integer NbBounds ) { myNbBounds = NbBounds; } //--------------------------------------------------------- //fonction : Add //--------------------------------------------------------- //--------------------------------------------------------- // Ajout d'une contrainte ponctuelle //--------------------------------------------------------- void GeomPlate_BuildPlateSurface:: Add(const Handle(GeomPlate_PointConstraint)& Cont) { myPntCont->Append(Cont); } //--------------------------------------------------------- //fonction : Perform // Calcul la surface de remplissage avec les contraintes chargees //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::Perform() { #ifdef DEB // Chronmetrage OSD_Chronometer Chrono; Chrono.Reset(); Chrono.Start(); #endif if (myNbBounds == 0) myNbBounds = myLinCont->Length(); myPlate.Init(); //===================================================================== //Declaration des variables. //===================================================================== Standard_Integer NTLinCont = myLinCont->Length(), NTPntCont = myPntCont->Length(), NbBoucle=0; // La variable NTPoint peut etre enlevee Standard_Boolean Fini=Standard_True; if ((NTLinCont+NTPntCont)==0) Standard_RangeError::Raise("GeomPlate : The number of constraints is null."); //====================================================================== // Surface Initiale //====================================================================== if (!mySurfInitIsGive) ComputeSurfInit(); else { if (NTLinCont>=2) { // Tableau des permutations pour conserver l'ordre initial voir TrierTab myInitOrder = new TColStd_HArray1OfInteger(1,NTLinCont); for (Standard_Integer l=1;l<=NTLinCont;l++) myInitOrder->SetValue(l,l); if (!CourbeJointive(myTol3d)) {// Standard_Failure::Raise("Curves are not joined"); #ifdef PLATE_DEB cout<<"WARNING : Courbes non jointives a "< 0)//Patch { mySense = new TColStd_HArray1OfInteger( 1, NTLinCont, 0 ); myInitOrder = new TColStd_HArray1OfInteger( 1, NTLinCont, 1 ); } } Standard_Real u1,v1,u2,v2; mySurfInit->Bounds(u1,v1,u2,v2); GeomAdaptor_Surface Surf(mySurfInit); myTolU = Surf.UResolution(myTol3d); myTolV = Surf.VResolution(myTol3d); myProj.Initialize(Surf,u1,v1,u2,v2, myTolU,myTolV); //====================================================================== // Projection des courbes //====================================================================== Standard_Integer i; Standard_Boolean Ok = Standard_True; for (i = 1; i <= NTLinCont; i++) if(myLinCont->Value(i)->Curve2dOnSurf().IsNull()) { Handle( Geom2d_Curve ) Curve2d = ProjectCurve( myLinCont->Value(i)->Curve3d() ); if (Curve2d.IsNull()) { Ok = Standard_False; break; } myLinCont->ChangeValue(i)->SetCurve2dOnSurf( Curve2d ); } if (!Ok) { GeomPlate_MakeApprox App(myGeomPlateSurface, myTol3d, 1, 3, 15*myTol3d, -1, GeomAbs_C0, 1.3); mySurfInit = App.Surface(); mySurfInit->Bounds(u1,v1,u2,v2); GeomAdaptor_Surface Surf(mySurfInit); myTolU = Surf.UResolution(myTol3d); myTolV = Surf.VResolution(myTol3d); myProj.Initialize(Surf,u1,v1,u2,v2, myTolU,myTolV); Ok = Standard_True; for (i = 1; i <= NTLinCont; i++) { Handle( Geom2d_Curve ) Curve2d = ProjectCurve( myLinCont->Value(i)->Curve3d() ); if (Curve2d.IsNull()) { Ok = Standard_False; break; } myLinCont->ChangeValue(i)->SetCurve2dOnSurf( Curve2d ); } if (!Ok) { mySurfInit = myPlanarSurfInit; mySurfInit->Bounds(u1,v1,u2,v2); GeomAdaptor_Surface Surf(mySurfInit); myTolU = Surf.UResolution(myTol3d); myTolV = Surf.VResolution(myTol3d); myProj.Initialize(Surf,u1,v1,u2,v2, myTolU,myTolV); for (i = 1; i <= NTLinCont; i++) myLinCont->ChangeValue(i)-> SetCurve2dOnSurf(ProjectCurve( myLinCont->Value(i)->Curve3d() ) ); } else { // Project the points for ( i=1;i<=NTPntCont;i++) { gp_Pnt P; myPntCont->Value(i)->D0(P); myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P)); } } } //====================================================================== // Projection des points //====================================================================== for ( i=1;i<=NTPntCont;i++) { if (! myPntCont->Value(i)->HasPnt2dOnSurf()) { gp_Pnt P; myPntCont->Value(i)->D0(P); myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P)); } } //====================================================================== // Nbre de point par courbe //====================================================================== if ((NTLinCont !=0)&&(myNbPtsOnCur!=0)) CalculNbPtsInit(); //====================================================================== // Gestion des incompatibilites entre courbes //====================================================================== Handle(GeomPlate_HArray1OfSequenceOfReal) PntInter; Handle(GeomPlate_HArray1OfSequenceOfReal) PntG1G1; if (NTLinCont != 0) { PntInter = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont); PntG1G1 = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont); Intersect(PntInter, PntG1G1); } //====================================================================== // Boucle pour obtenir une meilleur surface //====================================================================== myFree = !myIsLinear; do { #if PLATE_DEB if (Affich && NbBoucle) { cout<<"Resultats boucle"<< NbBoucle << endl; cout<<"DistMax="<SetBounds(Umin,Umax,Vmin,Vmax); Fini = VerifSurface(NbBoucle); if ((NbBoucle >= myNbIter)&&(!Fini)) { #ifdef PLATE_DEB cout<<"Warning objectif non atteint"<SetBounds(Umin,Umax,Vmin,Vmax); Fini = Standard_True; Standard_Real di,an,cu; VerifPoints(di,an,cu); } } while (!Fini); // Fin boucle pour meilleur surface #ifdef PLATE_DEB if (NTLinCont != 0) { cout<<"======== Resultats globaux ==========="<Value(c).Length(); Standard_Real U; if (NbPt<3) NbPt=4; else NbPt=myParCont->Value(c).Length(); gp_Vec v1i, v1f, v2i, v2f, v3i, v3f; gp_Pnt Pi, Pf; gp_Pnt2d P2d;Standard_Integer i; Handle(GeomPlate_CurveConstraint) LinCont = myLinCont->Value(c); switch (LinCont->Order()) { case 0 : for (i=1; iValue(c).Value(i) + myParCont->Value(c).Value(i+1) )/2; LinCont->D0(U, Pi); if (!LinCont->ProjectedCurve().IsNull()) P2d = LinCont->ProjectedCurve()->Value(U); else { if (!LinCont->Curve2dOnSurf().IsNull()) P2d = LinCont->Curve2dOnSurf()->Value(U); else P2d = ProjectPoint(Pi); } myGeomPlateSurface->D0( P2d.Coord(1), P2d.Coord(2), Pf); an->Init(0); courb->Init(0); d->ChangeValue(i) = Pf.Distance(Pi); } break; case 1 : for (i=1; iValue(c).Value(i) + myParCont->Value(c).Value(i+1) )/2; LinCont->D1(U, Pi, v1i, v2i); if (!LinCont->ProjectedCurve().IsNull()) P2d = LinCont->ProjectedCurve()->Value(U); else { if (!LinCont->Curve2dOnSurf().IsNull()) P2d = LinCont->Curve2dOnSurf()->Value(U); else P2d = ProjectPoint(Pi); } myGeomPlateSurface->D1( P2d.Coord(1), P2d.Coord(2), Pf, v1f, v2f); d->ChangeValue(i) = Pf.Distance(Pi); v3i = v1i^v2i; v3f=v1f^v2f; Standard_Real angle=v3f.Angle(v3i); if (angle>(PI/2)) an->ChangeValue(i) = PI -angle; else an->ChangeValue(i) = angle; courb->Init(0); } break; case 2 : Handle(Geom_Surface) Splate; Splate = Handle(Geom_Surface)::DownCast(myGeomPlateSurface); LocalAnalysis_SurfaceContinuity CG2; for (i=1; iValue(c).Value(i) + myParCont->Value(c).Value(i+1) )/2; LinCont->D0(U, Pi); if (!LinCont->ProjectedCurve().IsNull()) P2d = LinCont->ProjectedCurve()->Value(U); else { if (!LinCont->Curve2dOnSurf().IsNull()) P2d = LinCont->Curve2dOnSurf()->Value(U); else P2d = ProjectPoint(Pi); } GeomLProp_SLProps Prop (Splate, P2d.Coord(1), P2d.Coord(2), 2, 0.001); CG2.ComputeAnalysis(Prop, myLinCont->Value(c)->LPropSurf(U), GeomAbs_G2); d->ChangeValue(i)=CG2.C0Value(); an->ChangeValue(i)=CG2.G1Angle(); courb->ChangeValue(i)=CG2.G2CurvatureGap(); } break; } } //--------------------------------------------------------- // fonction : Disc2dContour //--------------------------------------------------------- void GeomPlate_BuildPlateSurface:: Disc2dContour ( const Standard_Integer nbp, TColgp_SequenceOfXY& Seq2d) { #ifdef PLATE_DEB if (nbp!=4) cout<<"nbp doit etre egal a 4 pour Disc2dContour"<Length(); Standard_Integer NTPntCont = myPntCont->Length(); gp_Pnt2d P2d; gp_XY UV; gp_Pnt PP; Standard_Real u1,v1,u2,v2; Standard_Integer i ; mySurfInit->Bounds(u1,v1,u2,v2); GeomAdaptor_Surface Surf(mySurfInit); myProj.Initialize(Surf,u1,v1,u2,v2, myTolU,myTolV); for( i=1; i<=NTPntCont; i++) if (myPntCont->Value(i)->Order()!=-1) { P2d = myPntCont->Value(i)->Pnt2dOnSurf(); UV.SetX(P2d.Coord(1)); UV.SetY(P2d.Coord(2)); Seq2d.Append(UV); } for(i=1; i<=NTCurve; i++) { Handle(GeomPlate_CurveConstraint) LinCont = myLinCont->Value(i); if (LinCont->Order()!=-1) { Standard_Integer NbPt=myParCont->Value(i).Length(); // premier point de contrainte (j=0) if (!LinCont->ProjectedCurve().IsNull()) P2d = LinCont->ProjectedCurve()->Value(myParCont->Value(i).Value(1)); else { if (!LinCont->Curve2dOnSurf().IsNull()) P2d = LinCont->Curve2dOnSurf()->Value(myParCont->Value(i).Value(1)); else { LinCont->D0(myParCont->Value(i).Value(1),PP); P2d = ProjectPoint(PP); } } UV.SetX(P2d.Coord(1)); UV.SetY(P2d.Coord(2)); Seq2d.Append(UV); for (Standard_Integer j=2; jValue(i).Value(j), Ujp1=myParCont->Value(i).Value(j+1); if (!LinCont->ProjectedCurve().IsNull()) P2d = LinCont->ProjectedCurve()->Value((Ujp1+3*Uj)/4); else { if (!LinCont->Curve2dOnSurf().IsNull()) P2d = LinCont->Curve2dOnSurf()->Value((Ujp1+3*Uj)/4); else { LinCont->D0((Ujp1+3*Uj)/4,PP); P2d = ProjectPoint(PP); } } UV.SetX(P2d.Coord(1)); UV.SetY(P2d.Coord(2)); Seq2d.Append(UV); // point milieu precedent if (!LinCont->ProjectedCurve().IsNull()) P2d = LinCont->ProjectedCurve()->Value((Ujp1+Uj)/2); else { if (!LinCont->Curve2dOnSurf().IsNull()) P2d = LinCont->Curve2dOnSurf()->Value((Ujp1+Uj)/2); else { LinCont->D0((Ujp1+Uj)/2,PP); P2d = ProjectPoint(PP); } } UV.SetX(P2d.Coord(1)); UV.SetY(P2d.Coord(2)); Seq2d.Append(UV); // point 3/4 precedent if (!LinCont->ProjectedCurve().IsNull()) P2d = LinCont->ProjectedCurve()->Value((3*Ujp1+Uj)/4); else { if (!LinCont->Curve2dOnSurf().IsNull()) P2d = LinCont->Curve2dOnSurf()->Value((3*Ujp1+Uj)/4); else { LinCont->D0((3*Ujp1+Uj)/4,PP); P2d = ProjectPoint(PP); } } UV.SetX(P2d.Coord(1)); UV.SetY(P2d.Coord(2)); Seq2d.Append(UV); // point de contrainte courant if (!LinCont->ProjectedCurve().IsNull()) P2d = LinCont->ProjectedCurve()->Value(Ujp1); else { if (!LinCont->Curve2dOnSurf().IsNull()) P2d = LinCont->Curve2dOnSurf()->Value(Ujp1); else { LinCont->D0(Ujp1,PP); P2d = ProjectPoint(PP); } } UV.SetX(P2d.Coord(1)); UV.SetY(P2d.Coord(2)); Seq2d.Append(UV); } } } } //--------------------------------------------------------- // fonction : Disc3dContour //--------------------------------------------------------- void GeomPlate_BuildPlateSurface:: Disc3dContour ( const Standard_Integer nbp, const Standard_Integer iordre, TColgp_SequenceOfXYZ& Seq3d) { #ifdef PLATE_DEB if (nbp!=4) cout<<"nbp doit etre egal a 4 pour Disc3dContour"<Bounds(u1,v1,u2,v2); GeomAdaptor_Surface Surf(mySurfInit); myProj.Initialize(Surf,u1,v1,u2,v2, Surf.UResolution(myTol3d), Surf.VResolution(myTol3d)); Standard_Integer NTCurve = myLinCont->Length(); Standard_Integer NTPntCont = myPntCont->Length(); // gp_Pnt2d P2d; gp_Pnt P3d; gp_Vec v1h,v2h,v3h; gp_XYZ Pos; Standard_Integer i ; for( i=1; i<=NTPntCont; i++) if (myPntCont->Value(i)->Order()!=-1) { if (iordre==0) { myPntCont->Value(i)->D0(P3d); Pos.SetX(P3d.X()); Pos.SetY(P3d.Y()); Pos.SetZ(P3d.Z()); Seq3d.Append(Pos); } else { myPntCont->Value(i)->D1(P3d,v1h,v2h); v3h=v1h^v2h; Pos.SetX(v3h.X()); Pos.SetY(v3h.Y()); Pos.SetZ(v3h.Z()); Seq3d.Append(Pos); } } for(i=1; i<=NTCurve; i++) if (myLinCont->Value(i)->Order()!=-1) { Standard_Integer NbPt=myParCont->Value(i).Length();; // premier point de contrainte (j=0) // Standard_Integer NbPt=myParCont->Length(); if (iordre==0) { myLinCont->Value(i)->D0(myParCont->Value(i).Value(1),P3d); Pos.SetX(P3d.X()); Pos.SetY(P3d.Y()); Pos.SetZ(P3d.Z()); Seq3d.Append(Pos); } else { myLinCont->Value(i)->D1(myParCont->Value(i).Value(1),P3d,v1h,v2h); v3h=v1h^v2h; Pos.SetX(v3h.X()); Pos.SetY(v3h.Y()); Pos.SetZ(v3h.Z()); Seq3d.Append(Pos); } for (Standard_Integer j=2; jValue(i).Value(j), Ujp1=myParCont->Value(i).Value(j+1); if (iordre==0) { // point 1/4 precedent myLinCont->Value(i)->D0((Ujp1+3*Uj)/4,P3d); Pos.SetX(P3d.X()); Pos.SetY(P3d.Y()); Pos.SetZ(P3d.Z()); Seq3d.Append(Pos); // point milieu precedent myLinCont->Value(i)->D0((Ujp1+Uj)/2,P3d); Pos.SetX(P3d.X()); Pos.SetY(P3d.Y()); Pos.SetZ(P3d.Z()); Seq3d.Append(Pos); // point 3/4 precedent myLinCont->Value(i)->D0((3*Ujp1+Uj)/4,P3d); Pos.SetX(P3d.X()); Pos.SetY(P3d.Y()); Pos.SetZ(P3d.Z()); Seq3d.Append(Pos); // point de contrainte courant myLinCont->Value(i)->D0(Ujp1,P3d); Pos.SetX(P3d.X()); Pos.SetY(P3d.Y()); Pos.SetZ(P3d.Z()); Seq3d.Append(Pos); } else { // point 1/4 precedent myLinCont->Value(i)->D1((Ujp1+3*Uj)/4,P3d,v1h,v2h); v3h=v1h^v2h; Pos.SetX(v3h.X()); Pos.SetY(v3h.Y()); Pos.SetZ(v3h.Z()); Seq3d.Append(Pos); // point milieu precedent myLinCont->Value(i)->D1((Ujp1+Uj)/2,P3d,v1h,v2h); v3h=v1h^v2h; Pos.SetX(v3h.X()); Pos.SetY(v3h.Y()); Pos.SetZ(v3h.Z()); Seq3d.Append(Pos); // point 3/4 precedent myLinCont->Value(i)->D1((3*Ujp1+Uj)/4,P3d,v1h,v2h); v3h=v1h^v2h; Pos.SetX(v3h.X()); Pos.SetY(v3h.Y()); Pos.SetZ(v3h.Z()); Seq3d.Append(Pos); // point de contrainte courant myLinCont->Value(i)->D1(Ujp1,P3d,v1h,v2h); v3h=v1h^v2h; Pos.SetX(v3h.X()); Pos.SetY(v3h.Y()); Pos.SetZ(v3h.Z()); Seq3d.Append(Pos); } } } } //--------------------------------------------------------- // fonction : IsDone //--------------------------------------------------------- Standard_Boolean GeomPlate_BuildPlateSurface::IsDone() const { return myPlate.IsDone(); } //--------------------------------------------------------- // fonction : Surface //--------------------------------------------------------- Handle_GeomPlate_Surface GeomPlate_BuildPlateSurface::Surface() const { return myGeomPlateSurface ; } //--------------------------------------------------------- // fonction : SurfInit //--------------------------------------------------------- Handle_Geom_Surface GeomPlate_BuildPlateSurface::SurfInit() const { return mySurfInit ; } //--------------------------------------------------------- // fonction : Sense //--------------------------------------------------------- Handle_TColStd_HArray1OfInteger GeomPlate_BuildPlateSurface::Sense() const { Standard_Integer NTCurve = myLinCont->Length(); Handle(TColStd_HArray1OfInteger) Sens = new TColStd_HArray1OfInteger(1, NTCurve); for (Standard_Integer i=1; i<=NTCurve; i++) Sens->SetValue(i,mySense->Value(myInitOrder->Value(i))); return Sens; } //--------------------------------------------------------- // fonction : Curve2d //--------------------------------------------------------- Handle_TColGeom2d_HArray1OfCurve GeomPlate_BuildPlateSurface::Curves2d() const { Standard_Integer NTCurve = myLinCont->Length(); Handle(TColGeom2d_HArray1OfCurve) C2dfin = new TColGeom2d_HArray1OfCurve(1,NTCurve); for (Standard_Integer i=1; i<=NTCurve; i++) C2dfin->SetValue(i,myLinCont->Value(myInitOrder->Value(i))->Curve2dOnSurf()); return C2dfin; } //--------------------------------------------------------- //fonction : Order //--------------------------------------------------------- Handle_TColStd_HArray1OfInteger GeomPlate_BuildPlateSurface::Order() const { Handle_TColStd_HArray1OfInteger result= new TColStd_HArray1OfInteger(1,myLinCont->Length()); for (Standard_Integer i=1;i<=myLinCont->Length();i++) result->SetValue(myInitOrder->Value(i),i); return result; } //--------------------------------------------------------- // fonction : G0Error //--------------------------------------------------------- Standard_Real GeomPlate_BuildPlateSurface::G0Error() const { return myG0Error; } //--------------------------------------------------------- // fonction : G1Error //--------------------------------------------------------- Standard_Real GeomPlate_BuildPlateSurface::G1Error() const { return myG1Error; } //--------------------------------------------------------- // fonction : G2Error //--------------------------------------------------------- Standard_Real GeomPlate_BuildPlateSurface::G2Error() const { return myG2Error; } //======================================================================= //function : G0Error //purpose : //======================================================================= Standard_Real GeomPlate_BuildPlateSurface::G0Error( const Standard_Integer Index ) { Handle( TColStd_HArray1OfReal ) tdistance = new TColStd_HArray1OfReal( 1, myNbPtsOnCur ); Handle( TColStd_HArray1OfReal ) tangle = new TColStd_HArray1OfReal( 1, myNbPtsOnCur ); Handle( TColStd_HArray1OfReal ) tcurvature = new TColStd_HArray1OfReal( 1, myNbPtsOnCur ); EcartContraintesMil( Index, tdistance, tangle, tcurvature ); Standard_Real MaxDistance = 0.; for (Standard_Integer i = 1; i <= myNbPtsOnCur; i++) if (tdistance->Value(i) > MaxDistance) MaxDistance = tdistance->Value(i); return MaxDistance; } //======================================================================= //function : G1Error //purpose : //======================================================================= Standard_Real GeomPlate_BuildPlateSurface::G1Error( const Standard_Integer Index ) { Handle( TColStd_HArray1OfReal ) tdistance = new TColStd_HArray1OfReal( 1, myNbPtsOnCur ); Handle( TColStd_HArray1OfReal ) tangle = new TColStd_HArray1OfReal( 1, myNbPtsOnCur ); Handle( TColStd_HArray1OfReal ) tcurvature = new TColStd_HArray1OfReal( 1, myNbPtsOnCur ); EcartContraintesMil( Index, tdistance, tangle, tcurvature ); Standard_Real MaxAngle = 0.; for (Standard_Integer i = 1; i <= myNbPtsOnCur; i++) if (tangle->Value(i) > MaxAngle) MaxAngle = tangle->Value(i); return MaxAngle; } //======================================================================= //function : G2Error //purpose : //======================================================================= Standard_Real GeomPlate_BuildPlateSurface::G2Error( const Standard_Integer Index ) { Handle( TColStd_HArray1OfReal ) tdistance = new TColStd_HArray1OfReal( 1, myNbPtsOnCur ); Handle( TColStd_HArray1OfReal ) tangle = new TColStd_HArray1OfReal( 1, myNbPtsOnCur ); Handle( TColStd_HArray1OfReal ) tcurvature = new TColStd_HArray1OfReal( 1, myNbPtsOnCur ); EcartContraintesMil( Index, tdistance, tangle, tcurvature ); Standard_Real MaxCurvature = 0.; for (Standard_Integer i = 1; i <= myNbPtsOnCur; i++) if (tcurvature->Value(i) > MaxCurvature) MaxCurvature = tcurvature->Value(i); return MaxCurvature; } Handle(GeomPlate_CurveConstraint) GeomPlate_BuildPlateSurface::CurveConstraint( const Standard_Integer order) const { return myLinCont->Value(order); } Handle(GeomPlate_PointConstraint) GeomPlate_BuildPlateSurface::PointConstraint( const Standard_Integer order) const { return myPntCont->Value(order); } //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// // ========================================================= // M E T H O D E S P R I V E S // ========================================================= //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// //======================================================================= // function : CourbeJointive // purpose : Effectue un chainage des courbes pour pouvoir calculer // la surface d'init avec la methode du flux maxi. // Retourne vrai si il s'agit d'un contour ferme. //======================================================================= Standard_Boolean GeomPlate_BuildPlateSurface:: CourbeJointive(const Standard_Real tolerance) { Standard_Integer nbf = myLinCont->Length(); Standard_Real Ufinal1,Uinit1,Ufinal2,Uinit2; mySense = new TColStd_HArray1OfInteger(1,nbf,0); Standard_Boolean result=Standard_True; Standard_Integer j=1,i; gp_Pnt P1,P2; while (j <= (myNbBounds-1)) { Standard_Integer a=0; i=j+1; if (i > myNbBounds) { result = Standard_False; a=2; } while (a<1) { if (i > myNbBounds) { result = Standard_False; a=2; } else { Uinit1=myLinCont->Value(j)->FirstParameter(); Ufinal1=myLinCont->Value(j)->LastParameter(); Uinit2=myLinCont->Value(i)->FirstParameter(); Ufinal2=myLinCont->Value(i)->LastParameter(); if (mySense->Value(j)==1) Ufinal1=Uinit1; myLinCont->Value(j)->D0(Ufinal1,P1); myLinCont->Value(i)->D0(Uinit2,P2); if (P1.Distance(P2)Value(j+1); myLinCont->SetValue(j+1,myLinCont->Value(i)); myLinCont->SetValue(i,tampon); Standard_Integer Tmp=myInitOrder->Value(j+1); //Voir fonction TrierTab pour le fonctionnement de myInitOrder myInitOrder->SetValue(j+1,myInitOrder->Value(i)); myInitOrder->SetValue(i,Tmp); }; a=2; mySense->SetValue(j+1,0); } else { myLinCont->Value(i)->D0(Ufinal2,P2); if (P1.Distance(P2)Value(j+1); myLinCont->SetValue(j+1,myLinCont->Value(i)); myLinCont->SetValue(i,tampon); Standard_Integer Tmp=myInitOrder->Value(j+1); //Voir fonction TrierTab pour le fonctionnement de myInitOrder myInitOrder->SetValue(j+1,myInitOrder->Value(i)); myInitOrder->SetValue(i,Tmp); }; a=2; mySense->SetValue(j+1,1); } } } i++; } j++; } Uinit1=myLinCont->Value( myNbBounds )->FirstParameter(); Ufinal1=myLinCont->Value( myNbBounds )->LastParameter(); Uinit2=myLinCont->Value(1)->FirstParameter(); Ufinal2=myLinCont->Value(1)->LastParameter(); myLinCont->Value( myNbBounds )->D0(Ufinal1,P1); myLinCont->Value(1)->D0(Uinit2,P2); if ((mySense->Value( myNbBounds )==0) &&(P1.Distance(P2)Value( myNbBounds )->D0(Uinit1,P1); if ((mySense->Value( myNbBounds )==1) &&(P1.Distance(P2)Length(), NTPntCont = myPntCont->Length(); // Tableau des permutation pour conserver l'ordre initial voir TrierTab if (NTLinCont != 0) { myInitOrder = new TColStd_HArray1OfInteger(1,NTLinCont); for (Standard_Integer i = 1; i <= NTLinCont; i++) myInitOrder->SetValue( i, i ); } Standard_Boolean CourbeJoint=Standard_False; if (NTLinCont != 0 && ((CourbeJoint = CourbeJointive( myTol3d )) != 0) && IsOrderG1()) { nopt = 3; // Tableau contenant le nuage de point pour le calcul du plan Standard_Integer NbPoint = 20, Discr = NbPoint/4, pnum = 0; Handle( TColgp_HArray1OfPnt ) Pts = new TColgp_HArray1OfPnt( 1, (NbPoint+1)*NTLinCont+NTPntCont ); TColgp_SequenceOfVec Vecs, NewVecs; GeomPlate_SequenceOfAij Aset; Standard_Real Uinit, Ufinal, Uif; gp_Vec LastVec; Standard_Integer i ; for ( i = 1; i <= NTLinCont; i++) { Standard_Integer Order = myLinCont->Value(i)->Order(); NewVecs.Clear(); Uinit = myLinCont->Value(i)->FirstParameter(); Ufinal = myLinCont->Value(i)->LastParameter(); Uif = Ufinal - Uinit; if (mySense->Value(i) == 1) { Uinit = Ufinal; Uif = -Uif; } gp_Vec Vec1, Vec2, Normal; Standard_Boolean ToReverse = Standard_False; if (i > 1 && Order >= GeomAbs_G1) { gp_Pnt P; myLinCont->Value(i)->D1( Uinit, P, Vec1, Vec2 ); Normal = Vec1 ^ Vec2; if (LastVec.IsOpposite( Normal, AngTol )) ToReverse = Standard_True; } for (Standard_Integer j = 0; j <= NbPoint; j++) { // Nombre de point par courbe = 20 // Repartition lineaire Standard_Real Inter = j*Uif/(NbPoint); if (Order < GeomAbs_G1 || j % Discr != 0) myLinCont->Value(i)->D0( Uinit+Inter, Pts->ChangeValue(++pnum) ); else { myLinCont->Value(i)->D1( Uinit+Inter, Pts->ChangeValue(++pnum), Vec1, Vec2 ); Normal = Vec1 ^ Vec2; Normal.Normalize(); if (ToReverse) Normal.Reverse(); Standard_Boolean isNew = Standard_True; Standard_Integer k ; for ( k = 1; k <= Vecs.Length(); k++) if (Vecs(k).IsEqual( Normal, LinTol, AngTol )) { isNew = Standard_False; break; } if (isNew) for (k = 1; k <= NewVecs.Length(); k++) if (NewVecs(k).IsEqual( Normal, LinTol, AngTol )) { isNew = Standard_False; break; } if (isNew) NewVecs.Append( Normal ); } } if (Order >= GeomAbs_G1) { isHalfSpace = GeomPlate_BuildAveragePlane::HalfSpace( NewVecs, Vecs, Aset, LinTol, AngTol ); if (! isHalfSpace) break; LastVec = Normal; } } //for (i = 1; i <= NTLinCont; i++) if (isHalfSpace) { for (i = 1; i <= NTPntCont; i++) { Standard_Integer Order = myPntCont->Value(i)->Order(); NewVecs.Clear(); gp_Vec Vec1, Vec2, Normal; if (Order < GeomAbs_G1) myPntCont->Value(i)->D0( Pts->ChangeValue(++pnum) ); else { myPntCont->Value(i)->D1( Pts->ChangeValue(++pnum), Vec1, Vec2 ); Normal = Vec1 ^ Vec2; Normal.Normalize(); Standard_Boolean isNew = Standard_True; for (Standard_Integer k = 1; k <= Vecs.Length(); k++) if (Vecs(k).IsEqual( Normal, LinTol, AngTol )) { isNew = Standard_False; break; } if (isNew) { NewVecs.Append( Normal ); isHalfSpace = GeomPlate_BuildAveragePlane::HalfSpace( NewVecs, Vecs, Aset, LinTol, AngTol ); if (! isHalfSpace) { NewVecs(1).Reverse(); isHalfSpace = GeomPlate_BuildAveragePlane::HalfSpace( NewVecs, Vecs, Aset, LinTol, AngTol ); } if (! isHalfSpace) break; } } } //for (i = 1; i <= NTPntCont; i++) if (isHalfSpace) { Standard_Boolean NullExist = Standard_True; while (NullExist) { NullExist = Standard_False; for (i = 1; i <= Vecs.Length(); i++) if (Vecs(i).SquareMagnitude() == 0.) { NullExist = Standard_True; Vecs.Remove(i); break; } } GeomPlate_BuildAveragePlane BAP( Vecs, Pts ); Standard_Real u1,u2,v1,v2; BAP.MinMaxBox(u1,u2,v1,v2); // On agrandit le bazar pour les projections Standard_Real du = u2 - u1; Standard_Real dv = v2 - v1; u1 -= du; u2 += du; v1 -= dv; v2 += dv; mySurfInit = new Geom_RectangularTrimmedSurface( BAP.Plane(), u1,u2,v1,v2 ); } } //if (isHalfSpace) if (!isHalfSpace) { #ifdef PLATE_DEB cout<Value(i)->Length(); for (i=1;i<=NTLinCont;i++) { Standard_Integer NbPoint= (Standard_Integer )( NTPoint*(myLinCont->Value(i)->Length())/LenT); if (NbPoint<10) NbPoint=10; Npt+=NbPoint; } // Tableau contenant le nuage de point pour le calcul du plan Handle(TColgp_HArray1OfPnt) Pts = new TColgp_HArray1OfPnt(1,20*NTLinCont+NTPntCont); Standard_Integer NbPoint=20; Standard_Real Uinit,Ufinal, Uif; for ( i=1;i<=NTLinCont;i++) { Uinit=myLinCont->Value(i)->FirstParameter(); Ufinal=myLinCont->Value(i)->LastParameter(); Uif=Ufinal-Uinit; if (mySense->Value(i) == 1) { Uinit = Ufinal; Uif = -Uif; } for (Standard_Integer j=0; jValue(i)->D0(Uinit+Inter,P); Pts->SetValue(Np++,P); } } for (i=1;i<=NTPntCont;i++) { gp_Pnt P; myPntCont->Value(i)->D0(P); Pts->SetValue(Np++,P); } if (!CourbeJoint) myNbBounds = 0; GeomPlate_BuildAveragePlane BAP( Pts, NbPoint*myNbBounds, myTol3d/1000, popt, nopt ); if (!BAP.IsPlane()) Standard_Failure::Raise("the initial surface is not a plane."); Standard_Real u1,u2,v1,v2; BAP.MinMaxBox(u1,u2,v1,v2); // On agrandit le bazar pour les projections Standard_Real du = u2 - u1; Standard_Real dv = v2 - v1; u1 -= du; u2 += du; v1 -= dv; v2 += dv; mySurfInit = new Geom_RectangularTrimmedSurface(BAP.Plane(),u1,u2,v1,v2); } //if (nopt != 3) //Comparing metrics of curves and projected curves if (NTLinCont != 0 && myIsLinear) { Handle( Geom_Surface ) InitPlane = (Handle( Geom_RectangularTrimmedSurface )::DownCast(mySurfInit))->BasisSurface(); Standard_Real Ratio, R1 = 2., R2 = 0.6; //R1 = 3, R2 = 0.5;//R1 = 1.4, R2 = 0.8; //R1 = 5., R2 = 0.2; Handle( GeomAdaptor_HSurface ) hsur = new GeomAdaptor_HSurface( InitPlane ); Standard_Integer NbPoint = 20; // gp_Pnt P; // gp_Vec DerC, DerCproj, DU, DV; // gp_Pnt2d P2d; // gp_Vec2d DProj; for (Standard_Integer i = 1; i <= NTLinCont && myIsLinear; i++) { Standard_Real FirstPar = myLinCont->Value(i)->FirstParameter(); Standard_Real LastPar = myLinCont->Value(i)->LastParameter(); Standard_Real Uif = (LastPar - FirstPar)/(NbPoint); Handle( Adaptor3d_HCurve ) Curve = myLinCont->Value(i)->Curve3d(); ProjLib_CompProjectedCurve Projector( hsur, Curve, myTol3d, myTol3d ); Handle( ProjLib_HCompProjectedCurve ) ProjCurve = new ProjLib_HCompProjectedCurve(); ProjCurve->Set( Projector ); Adaptor3d_CurveOnSurface AProj(ProjCurve, hsur); gp_Pnt P; gp_Vec DerC, DerCproj; for (Standard_Integer j = 1 ; j < NbPoint && myIsLinear; j++) { Standard_Real Inter = FirstPar + j*Uif; Curve->D1(Inter, P, DerC); AProj.D1(Inter, P, DerCproj); Standard_Real A1 = DerC.Magnitude(); Standard_Real A2 = DerCproj.Magnitude(); if (A2 <= 1.e-20) { Ratio = 1.e20; } else { Ratio = A1 / A2; } if (Ratio > R1 || Ratio < R2) { myIsLinear = Standard_False; break; } } } #if PLATE_DEB if (! myIsLinear) cout <<"Metrics are too different :"<< Ratio<Bounds(u1,v1,u2,v2); GeomAdaptor_Surface Surf(mySurfInit); myTolU = Surf.UResolution(myTol3d); myTolV = Surf.VResolution(myTol3d); myProj.Initialize(Surf,u1,v1,u2,v2, myTolU,myTolV); //====================================================================== // Projection des courbes //====================================================================== Standard_Integer i; for (i = 1; i <= NTLinCont; i++) if (myLinCont->Value(i)->Curve2dOnSurf().IsNull()) myLinCont->ChangeValue(i)->SetCurve2dOnSurf(ProjectCurve(myLinCont->Value(i)->Curve3d())); //====================================================================== // Projection des points //====================================================================== for (i = 1; i<=NTPntCont; i++) { gp_Pnt P; myPntCont->Value(i)->D0(P); if (!myPntCont->Value(i)->HasPnt2dOnSurf()) myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P)); } //====================================================================== // Nbre de point par courbe //====================================================================== if ((NTLinCont !=0)&&(myNbPtsOnCur!=0)) CalculNbPtsInit(); //====================================================================== // Gestion des incompatibilites entre courbes //====================================================================== Handle(GeomPlate_HArray1OfSequenceOfReal) PntInter; Handle(GeomPlate_HArray1OfSequenceOfReal) PntG1G1; if (NTLinCont != 0) { PntInter = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont); PntG1G1 = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont); Intersect(PntInter, PntG1G1); } //==================================================================== // Discretisation des courbes //==================================================================== Discretise(PntInter, PntG1G1); //==================================================================== //Preparation des points de contrainte pour plate //==================================================================== LoadCurve( 0, 0 ); if (myPntCont->Length() != 0) LoadPoint( 0, 0 ); //==================================================================== //Resolution de la surface //==================================================================== myPlate.SolveTI(2, ComputeAnisotropie()); if (!myPlate.IsDone()) Standard_Failure::Raise("GeomPlate : abort calcul of Plate."); myGeomPlateSurface = new GeomPlate_Surface( mySurfInit, myPlate ); GeomPlate_MakeApprox App(myGeomPlateSurface, myTol3d, 1, 3, 15*myTol3d, -1, GeomAbs_C0); mySurfInit = App.Surface(); mySurfInitIsGive = Standard_True; myPlate.Init(); // Reset for (i=1;i<=NTLinCont;i++) { Handle( Geom2d_Curve ) NullCurve; NullCurve.Nullify(); myLinCont->ChangeValue(i)->SetCurve2dOnSurf( NullCurve); } } #if DRAW if (Affich) { char* name = new char[100]; sprintf(name,"surfinit_%d",++NbPlan); DrawTrSurf::Set(name, mySurfInit); } #endif } //--------------------------------------------------------- // fonction : Intersect //--------------------------------------------------------- // Recherche les intersections entre les courbe 2d // Si intersection et compatible( dans le cas G1-G1) // enleve le point sur une des deux courbes //--------------------------------------------------------- void GeomPlate_BuildPlateSurface:: Intersect(Handle(GeomPlate_HArray1OfSequenceOfReal)& PntInter, Handle(GeomPlate_HArray1OfSequenceOfReal)& PntG1G1) { Standard_Integer NTLinCont = myLinCont->Length(); Geom2dInt_GInter Intersection; Geom2dAdaptor_Curve Ci, Cj; IntRes2d_IntersectionPoint int2d; gp_Pnt P1,P2; gp_Pnt2d P2d; gp_Vec2d V2d; // if (!mySurfInitIsGive) for (Standard_Integer i=1;i<=NTLinCont;i++) { //Standard_Real NbPnt_i=myLinCont->Value(i)->NbPoints(); // Cherche l'intersection avec chaque courbe y compris elle meme. Ci.Load(myLinCont->Value(i)->Curve2dOnSurf()); for(Standard_Integer j=i; j<=NTLinCont; j++) { Cj.Load(myLinCont->Value(j)->Curve2dOnSurf()); if (i==j) Intersection.Perform(Ci, myTol2d*10, myTol2d*10); else Intersection.Perform(Ci, Cj, myTol2d*10, myTol2d*10); if (!Intersection.IsEmpty()) { // il existe une intersection Standard_Integer nbpt = Intersection.NbPoints(); // nombre de point d'intersection for (Standard_Integer k = 1; k <= nbpt; k++) { int2d = Intersection.Point(k); myLinCont->Value(i)->D0(int2d.ParamOnFirst(),P1); myLinCont->Value(j)->D0(int2d.ParamOnSecond(),P2); #if PLATE_DEB if (Affich> 1) { cout << " Intersection "<< k << " entre " << i << " &" << j << endl; cout <<" Distance = "<< P1.Distance(P2) << endl; } #endif if (P1.Distance( P2 ) < myTol3d) { // L'intersection 2d correspond a des points 3d proches // On note l'intervalle dans lequel il faudra enlever // un point pour eviter les doublons ce qui fait une // erreur dans plate. le point sur la courbe i est enleve // on conserve celui de la courbe j // la longueur de l'intervalle est une longueur 2d // correspondant en 3d a myTol3d Standard_Real tolint = Ci.Resolution(myTol3d); Ci.D1(int2d.ParamOnFirst(),P2d, V2d); Standard_Real aux = V2d.Magnitude(); if (aux > 1.e-7) { aux = myTol3d/aux; if (aux > 100*tolint) tolint*=100; else tolint = aux; } else tolint*=100; PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() - tolint); PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() + tolint); // Si G1-G1 if ( (myLinCont->Value(i)->Order()==1) &&(myLinCont->Value(j)->Order()==1)) { gp_Vec v11,v12,v13,v14,v15,v16,v21,v22,v23,v24,v25,v26; myLinCont->Value(i)->D2( int2d.ParamOnFirst(), P1, v11, v12, v13, v14, v15 ); myLinCont->Value(j)->D2( int2d.ParamOnSecond(), P2, v21, v22, v23, v24, v25 ); v16=v11^v12;v26=v21^v22; Standard_Real ant=v16.Angle(v26); if (ant>(PI/2)) ant= PI -ant; if ((Abs(v16*v15-v16*v25)>(myTol3d/1000)) ||(Abs(ant)>myTol3d/1000)) // Pas compatible ==> on enleve une zone en // contrainte G1 correspondant // a une tolerance 3D de 0.01 { Standard_Real coin; Standard_Real Tol= 100 * myTol3d; Standard_Real A1; gp_Pnt2d P1,P2; gp_Vec2d V1,V2; myLinCont->Value(i)->Curve2dOnSurf()->D1( int2d.ParamOnFirst(), P1, V1); myLinCont->Value(j)->Curve2dOnSurf()->D1( int2d.ParamOnSecond(), P2, V2); A1 = V1.Angle(V2); if (A1>(PI/2)) A1= PI - A1; if (Abs(Abs(A1)-PI)ChangeValue(i).Append(Par1); PntG1G1->ChangeValue(i).Append(Par2); coin = Cj.Resolution(Tol); Par1=int2d.ParamOnSecond()-coin; Par2=int2d.ParamOnSecond()+coin; // Stockage de l'intervalle pour la courbe j PntG1G1->ChangeValue(j).Append(Par1); PntG1G1->ChangeValue(j).Append(Par2); } } // Si G0-G1 if ((myLinCont->Value(i)->Order()==0 && myLinCont->Value(j)->Order()==1) || (myLinCont->Value(i)->Order()==1 && myLinCont->Value(j)->Order()==0)) { gp_Vec vec, vecU, vecV, N; if (myLinCont->Value(i)->Order() == 0) { Handle( Adaptor3d_HCurve ) theCurve = myLinCont->Value(i)->Curve3d(); theCurve->D1( int2d.ParamOnFirst(), P1, vec ); myLinCont->Value(j)->D1( int2d.ParamOnSecond(), P2, vecU, vecV ); } else { Handle( Adaptor3d_HCurve ) theCurve = myLinCont->Value(j)->Curve3d(); theCurve->D1( int2d.ParamOnSecond(), P2, vec ); myLinCont->Value(i)->D1( int2d.ParamOnFirst(), P1, vecU, vecV ); } N = vecU ^ vecV; Standard_Real Angle = vec.Angle( N ); Angle = Abs( PI/2-Angle ); if (Angle > myTolAng/10.) //????????? //if (Abs( scal ) > myTol3d/100) { // Pas compatible ==> on enleve une zone en // contrainte G0 et G1 correspondant // a une tolerance 3D de 0.01 Standard_Real coin; Standard_Real Tol= 100 * myTol3d; Standard_Real A1; gp_Pnt2d P1,P2; gp_Vec2d V1,V2; myLinCont->Value(i)->Curve2dOnSurf()->D1( int2d.ParamOnFirst(), P1, V1); myLinCont->Value(j)->Curve2dOnSurf()->D1( int2d.ParamOnSecond(), P2, V2); A1 = V1.Angle( V2 ); if (A1 > PI/2) A1= PI - A1; if (Abs(Abs(A1) - PI) < myTolAng) Tol = 100000 * myTol3d; #if PLATE_DEB if (Affich) cout <<"Angle entre Courbe "<Value(i)->Order() == 1) { coin = Ci.Resolution(Tol); coin *= Angle / myTolAng * 10.; #if PLATE_DEB cout<ChangeValue(i).Append( Par1 ); PntG1G1->ChangeValue(i).Append( Par2 ); } else { coin = Cj.Resolution(Tol); coin *= Angle / myTolAng * 10.; #if PLATE_DEB cout<ChangeValue(j).Append( Par1 ); PntG1G1->ChangeValue(j).Append( Par2 ); } } } } //if (P1.Distance( P2 ) < myTol3d) else { //L'intersection 2d correspond a des points 3d eloigne // On note l'intervalle dans lequel il faudra enlever // des points pour eviter les doublons ce qui fait une // erreur dans plate. le point sur la courbe i est enleve // on conserve celui de la courbe j. // la longueur de l'intervalle est une longueur 2d // correspondant a la distance des points en 3d a myTol3d Standard_Real tolint, Dist; Dist = P1.Distance(P2); tolint = Ci.Resolution(Dist); PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() - tolint); PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() + tolint); if (j!=i) { tolint = Cj.Resolution(Dist); PntInter->ChangeValue(j). Append( int2d.ParamOnSecond() - tolint); PntInter->ChangeValue(j). Append( int2d.ParamOnSecond() + tolint); } #ifdef PLATE_DEB cout<<"Attention: Deux points 3d ont la meme projection dist=" < 1) { Handle(Draw_Marker3D) mark = new (Draw_Marker3D)(P1, Draw_X, Draw_vert); char* name = new char[100]; sprintf(name,"mark_%d",++NbMark); Draw::Set(name, mark); } #endif } } } } } } //--------------------------------------------------------- // fonction : Discretise //--------------------------------------------------------- // Discretise les courbes selon le parametre de celle-ci // le tableau des sequences Parcont contiendera tous les // parametres des points sur les courbes // Le champ myPlateCont contiendera les parametre des points envoye a plate // il excluera les points en double et les zones imcompatibles. // La premiere partie correspond a la verfication de la compatibilite // et a la supprssion des points en double. //--------------------------------------------------------- void GeomPlate_BuildPlateSurface:: Discretise(const Handle(GeomPlate_HArray1OfSequenceOfReal)& PntInter, const Handle(GeomPlate_HArray1OfSequenceOfReal)& PntG1G1) { Standard_Integer NTLinCont = myLinCont->Length(); Standard_Boolean ACR; Handle(Geom2d_Curve) C2d; Geom2dAdaptor_Curve AC2d; // Handle(Adaptor_HCurve2d) HC2d; Handle(Law_Interpol) acrlaw = new (Law_Interpol) (); myPlateCont = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont); myParCont = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont); //=========================================================================== // Construction du tableau contenant les parametres des points de contrainte //=========================================================================== Standard_Real Uinit, Ufinal, Length2d=0, Inter; Standard_Real CurLength; Standard_Integer NbPnt_i, NbPtInter, NbPtG1G1; Handle(GeomPlate_CurveConstraint) LinCont; for (Standard_Integer i=1;i<=NTLinCont;i++) { LinCont = myLinCont->Value(i); Uinit=LinCont->FirstParameter(); Ufinal=LinCont->LastParameter(); // HC2d=LinCont->ProjectedCurve(); // if(HC2d.IsNull()) // ACR = (!HC2d.IsNull() || !C2d.IsNull()); C2d= LinCont->Curve2dOnSurf(); ACR = (!C2d.IsNull()); if (ACR) { // On Construit une loi proche de l'abscisse curviligne if(!C2d.IsNull()) AC2d.Load(C2d); // AC2d.Load(LinCont->Curve2dOnSurf()); Standard_Integer ii, Nbint = 20; Standard_Real U; TColgp_Array1OfPnt2d tabP2d(1, Nbint+1); tabP2d(1).SetY(Uinit); tabP2d(1).SetX(0.); tabP2d(Nbint+1).SetY(Ufinal); /* if (!HC2d.IsNull()) Length2d = GCPnts_AbscissaPoint::Length(HC2d->Curve2d(), Uinit, Ufinal); else*/ Length2d = GCPnts_AbscissaPoint::Length(AC2d, Uinit, Ufinal); tabP2d(Nbint+1).SetX(Length2d); for (ii = 2; ii<= Nbint; ii++) { U = Uinit + (Ufinal-Uinit)*((1-Cos((ii-1)*PI/(Nbint)))/2); tabP2d(ii).SetY(U); /* if (!HC2d.IsNull()) { Standard_Real L = GCPnts_AbscissaPoint::Length(HC2d->Curve2d(), Uinit, U); tabP2d(ii).SetX(L); } else*/ tabP2d(ii).SetX(GCPnts_AbscissaPoint::Length(AC2d, Uinit, U)); } acrlaw->Set(tabP2d); } NbPnt_i = (Standard_Integer )( LinCont->NbPoints()); NbPtInter= PntInter->Value(i).Length(); NbPtG1G1= PntG1G1->Value(i).Length(); #if PLATE_DEB if (Affich > 1) { cout << "Courbe : " << i << endl; cout << " NbPnt, NbPtInter, NbPtG1G1 :" << NbPnt_i << ", " << NbPtInter << ", " << NbPtG1G1 << endl; } #endif for (Standard_Integer j=1; j<=NbPnt_i; j++) { // repartition des points en cosinus selon l'ACR 2d // Afin d'eviter les points d'acumulation dans le 2d //Inter=Uinit+(Uif)*((-cos(PI*((j-1)/(NbPnt_i-1)))+1)/2); if (j==NbPnt_i) Inter=Ufinal;//pour parer au bug sur sun else if (ACR) { CurLength = Length2d*(1-Cos((j-1)*PI/(NbPnt_i-1)))/2; Inter = acrlaw->Value(CurLength); } else { Inter=Uinit+(Ufinal-Uinit)*((1-Cos((j-1)*PI/(NbPnt_i-1)))/2); } myParCont->ChangeValue(i).Append(Inter);// on ajoute le point if (NbPtInter!=0) { for(Standard_Integer l=1;l<=NbPtInter;l+=2) //on cherche si le point Inter est dans l'intervalle //PntInter[i] PntInter[i+1] //auquelle cas il ne faudrait pas le stocker (pb de doublons) { if ((Inter>PntInter->Value(i).Value(l)) &&(InterValue(i).Value(l+1))) { l=NbPtInter+2; // pour sortir de la boucle sans stocker le point } else { if (l+1>=NbPtInter) { // on a parcouru tout le tableau : Le point // n'appartient pas a un interval point commun if (NbPtG1G1!=0) // est qu'il existe un intervalle incompatible for(Standard_Integer k=1;k<=NbPtG1G1;k+=2) { if ((Inter>PntG1G1->Value(i).Value(k)) &&(InterValue(i).Value(k+1))) { k=NbPtG1G1+2; // pour sortir de la boucle // Ajouter les points de contrainte G0 gp_Pnt P3d,PP,Pdif; gp_Pnt2d P2d; AC2d.D0(Inter, P2d); LinCont->D0(Inter,P3d); mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP); Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1), -PP.Coord(2)+P3d.Coord(2), -PP.Coord(3)+P3d.Coord(3)); Plate_PinpointConstraint PC(P2d.XY(), Pdif.XYZ(),0,0); myPlate.Load(PC); } else // le point n'appartient pas a un interval G1 { if (k+1>=NbPtG1G1) { myPlateCont->ChangeValue(i).Append(Inter); // on ajoute le point } } } else { myPlateCont->ChangeValue(i).Append(Inter); // on ajoute le point } } } } } else { if (NbPtG1G1!=0) // est qu'il existe un intervalle incompatible for(Standard_Integer k=1;k<=NbPtG1G1;k+=2) { if ((Inter>PntG1G1->Value(i).Value(k)) &&(InterValue(i).Value(k+1))) { k=NbPtG1G1+2; // pour sortir de la boucle // Ajouter les points de contrainte G0 gp_Pnt P3d,PP,Pdif; gp_Pnt2d P2d; AC2d.D0(Inter, P2d); LinCont->D0(Inter,P3d); mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP); Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1), -PP.Coord(2)+P3d.Coord(2), -PP.Coord(3)+P3d.Coord(3)); Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0); myPlate.Load(PC); } else // le point n'appartient pas a un intervalle G1 { if (k+1>=NbPtG1G1) { myPlateCont->ChangeValue(i).Append(Inter); // on ajoute le point } } } else { if ( ( (!mySurfInitIsGive) &&(Geom2dAdaptor_Curve(LinCont->Curve2dOnSurf()).GetType()!=GeomAbs_Circle)) || ( (j>1) &&(jChangeValue(i).Append(Inter);// on ajoute le point } } } } } //--------------------------------------------------------- // fonction : CalculNbPtsInit //--------------------------------------------------------- // Calcul du nombre de points par courbe en fonction de la longueur // pour la premiere iteration //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::CalculNbPtsInit () { Standard_Real LenT=0; Standard_Integer NTLinCont=myLinCont->Length(); Standard_Integer NTPoint=(Standard_Integer )( myNbPtsOnCur*NTLinCont); Standard_Integer i ; for ( i=1;i<=NTLinCont;i++) LenT+=myLinCont->Value(i)->Length(); for ( i=1;i<=NTLinCont;i++) { Standard_Integer Cont=myLinCont->Value(i)->Order(); switch(Cont) { case 0 : // Cas G0 *1.2 myLinCont->ChangeValue(i)->SetNbPoints( Standard_Integer(1.2*NTPoint*(myLinCont->Value(i)->Length())/LenT)); break; case 1 : // Cas G1 *1 myLinCont->ChangeValue(i)->SetNbPoints( Standard_Integer(NTPoint*(myLinCont->Value(i)->Length())/LenT)); break; case 2 : // Cas G2 *0.7 myLinCont->ChangeValue(i)->SetNbPoints( Standard_Integer(0.7*NTPoint*(myLinCont->Value(i)->Length())/LenT)); break; } if (myLinCont->Value(i)->NbPoints()<3) myLinCont->ChangeValue(i)->SetNbPoints(3); } } //--------------------------------------------------------- // fonction : LoadCurve //--------------------------------------------------------- // A partir du tableau myParCont on charge tous les points notes dans plate //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::LoadCurve(const Standard_Integer NbBoucle, const Standard_Integer OrderMax ) { gp_Pnt P3d,Pdif,PP; gp_Pnt2d P2d; Standard_Integer NTLinCont=myLinCont->Length(), i, j; Standard_Integer Tang, Nt; for (i=1; i<=NTLinCont; i++){ Handle(GeomPlate_CurveConstraint) CC = myLinCont->Value(i); if (CC ->Order()!=-1) { Tang = Min(CC->Order(), OrderMax); Nt = myPlateCont->Value(i).Length(); if (Tang!=-1) for (j=1; j<=Nt; j++) {// Chargement des points G0 sur les frontieres CC ->D0(myPlateCont->Value(i).Value(j),P3d); if (!CC->ProjectedCurve().IsNull()) P2d = CC->ProjectedCurve()->Value(myPlateCont->Value(i).Value(j)); else { if (!CC->Curve2dOnSurf().IsNull()) P2d = CC->Curve2dOnSurf()->Value(myPlateCont->Value(i).Value(j)); else P2d = ProjectPoint(P3d); } mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP); Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1), -PP.Coord(2)+P3d.Coord(2), -PP.Coord(3)+P3d.Coord(3)); Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0); myPlate.Load(PC); // Chargement des points G1 if (Tang==1) { // ==1 gp_Vec V1,V2,V3,V4; CC->D1( myPlateCont->Value(i).Value(j), PP, V1, V2 ); mySurfInit->D1( P2d.Coord(1), P2d.Coord(2), PP, V3, V4 ); Plate_D1 D1final(V1.XYZ(),V2.XYZ()); Plate_D1 D1init(V3.XYZ(),V4.XYZ()); if (! myFree) { Plate_GtoCConstraint GCC( P2d.XY(), D1init, D1final ); myPlate.Load(GCC); } else if (NbBoucle == 1) { Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), D1init, D1final ); myPlate.Load( FreeGCC ); } else { gp_Vec DU, DV, Normal, DerPlateU, DerPlateV; Normal = V1^V2; //Normal.Normalize(); Standard_Real norm = Normal.Magnitude(); if (norm > 1.e-12) Normal /= norm; DerPlateU = myPrevPlate.EvaluateDerivative( P2d.XY(), 1, 0 ); DerPlateV = myPrevPlate.EvaluateDerivative( P2d.XY(), 0, 1 ); DU.SetLinearForm(-(V3 + DerPlateU).Dot(Normal), Normal, DerPlateU); DV.SetLinearForm(-(V4 + DerPlateV).Dot(Normal), Normal, DerPlateV); Plate_PinpointConstraint PinU( P2d.XY(), DU.XYZ(), 1, 0 ); Plate_PinpointConstraint PinV( P2d.XY(), DV.XYZ(), 0, 1 ); myPlate.Load( PinU ); myPlate.Load( PinV ); } } // Chargement des points G2 if (Tang==2) // ==2 { gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10; CC->D2(myPlateCont->Value(i).Value(j), PP,V1,V2,V5,V6,V7); mySurfInit->D2(P2d.Coord(1),P2d.Coord(2),PP,V3,V4,V8,V9,V10); Plate_D1 D1final(V1.XYZ(),V2.XYZ()); Plate_D1 D1init(V3.XYZ(),V4.XYZ()); Plate_D2 D2final (V5.XYZ(),V6.XYZ(),V7.XYZ()); Plate_D2 D2init (V8.XYZ(),V9.XYZ(),V10.XYZ()); // if (! myFree) // { Plate_GtoCConstraint GCC(P2d.XY(),D1init,D1final,D2init,D2final); myPlate.Load(GCC); // } // else // Good but too expansive // { // Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), // D1init, D1final, D2init, D2final ); // myPlate.Load( FreeGCC ); // } } } } } } //--------------------------------------------------------- //fonction : LoadPoint //--------------------------------------------------------- //void GeomPlate_BuildPlateSurface::LoadPoint(const Standard_Integer NbBoucle, void GeomPlate_BuildPlateSurface::LoadPoint(const Standard_Integer , const Standard_Integer OrderMax) { gp_Pnt P3d,Pdif,PP; gp_Pnt2d P2d; Standard_Integer NTPntCont=myPntCont->Length(); Standard_Integer Tang, i; // gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10; gp_Vec V1,V2,V3,V4; // Chargement des points de contraintes ponctuel for (i=1;i<=NTPntCont;i++) { myPntCont->Value(i)->D0(P3d); P2d = myPntCont->Value(i)->Pnt2dOnSurf(); mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP); Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1), -PP.Coord(2)+P3d.Coord(2), -PP.Coord(3)+P3d.Coord(3)); Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0); myPlate.Load(PC); Tang = Min(myPntCont->Value(i)->Order(), OrderMax); if (Tang==1) {// ==1 myPntCont->Value(i)->D1(PP,V1,V2); mySurfInit->D1(P2d.Coord(1),P2d.Coord(2),PP,V3,V4); Plate_D1 D1final(V1.XYZ(),V2.XYZ()); Plate_D1 D1init(V3.XYZ(),V4.XYZ()); if (! myFree) { Plate_GtoCConstraint GCC( P2d.XY(), D1init, D1final ); myPlate.Load( GCC ); } else { Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), D1init, D1final ); myPlate.Load( FreeGCC ); } } // Chargement des points G2 GeomPlate_PlateG0Criterion if (Tang==2) // ==2 { gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10; myPntCont->Value(i)->D2(PP,V1,V2,V5,V6,V7); // gp_Vec Tv2 = V1^V2; mySurfInit->D2(P2d.Coord(1),P2d.Coord(2),PP,V3,V4,V8,V9,V10); Plate_D1 D1final(V1.XYZ(),V2.XYZ()); Plate_D1 D1init(V3.XYZ(),V4.XYZ()); Plate_D2 D2final (V5.XYZ(),V6.XYZ(),V7.XYZ()); Plate_D2 D2init (V8.XYZ(),V9.XYZ(),V10.XYZ()); // if (! myFree) // { Plate_GtoCConstraint GCC( P2d.XY(), D1init, D1final, D2init, D2final ); myPlate.Load( GCC ); // } // else // Good but too expansive // { // Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), D1init, D1final, D2init//, D2final ); // myPlate.Load( FreeGCC ); // } } } } //--------------------------------------------------------- //fonction : VerifSurface //--------------------------------------------------------- Standard_Boolean GeomPlate_BuildPlateSurface:: VerifSurface(const Standard_Integer NbBoucle) { //====================================================================== // Calcul des erreurs //====================================================================== Standard_Integer NTLinCont=myLinCont->Length(); Standard_Boolean Result=Standard_True; // variable pour les calculs d erreur myG0Error=0,myG1Error =0, myG2Error=0; for (Standard_Integer i=1;i<=NTLinCont;i++) { Handle(GeomPlate_CurveConstraint) LinCont; LinCont = myLinCont->Value(i); if (LinCont->Order()!=-1) { Standard_Integer NbPts_i = myParCont->Value(i).Length(); if (NbPts_i<3) NbPts_i=4; Handle(TColStd_HArray1OfReal) tdist = new TColStd_HArray1OfReal(1,NbPts_i-1); Handle(TColStd_HArray1OfReal) tang = new TColStd_HArray1OfReal(1,NbPts_i-1); Handle(TColStd_HArray1OfReal) tcourb = new TColStd_HArray1OfReal(1,NbPts_i-1); EcartContraintesMil (i,tdist,tang,tcourb); Standard_Real diffDistMax=0,SdiffDist=0; Standard_Real diffAngMax=0,SdiffAng=0; Standard_Integer NdiffDist=0,NdiffAng=0; for (Standard_Integer j=1;jValue(j)>myG0Error) myG0Error=tdist->Value(j); if (tang->Value(j)>myG1Error) myG1Error=tang->Value(j); if (tcourb->Value(j)>myG2Error) myG2Error=tcourb->Value(j); Standard_Real U; if (myParCont->Value(i).Length()>3) U=(myParCont->Value(i).Value(j)+myParCont->Value(i).Value(j+1))/2; else U=LinCont->FirstParameter()+ ( LinCont->LastParameter()-LinCont->FirstParameter())*(j-1)/(NbPts_i-2); Standard_Real diffDist=tdist->Value(j)-LinCont->G0Criterion(U),diffAng; if (LinCont->Order()>0) diffAng=tang->Value(j)-LinCont->G1Criterion(U); else diffAng=0; // recherche de la variation de l'erreur maxi et calcul de la moyenne if (diffDist>0) { diffDist = diffDist/LinCont->G0Criterion(U); if (diffDist>diffDistMax) diffDistMax = diffDist; SdiffDist+=diffDist; NdiffDist++; #if DRAW if ((Affich) && (NbBoucle == myNbIter)) { gp_Pnt P; gp_Pnt2d P2d; LinCont->D0(U,P); Handle(Draw_Marker3D) mark = new (Draw_Marker3D)(P, Draw_X, Draw_orange); char* name = new char[100]; sprintf(name,"mark_%d",++NbMark); Draw::Set(name, mark); if (!LinCont->ProjectedCurve().IsNull()) P2d = LinCont->ProjectedCurve()->Value(U); else { if (!LinCont->Curve2dOnSurf().IsNull()) P2d = LinCont->Curve2dOnSurf()->Value(U); else P2d = ProjectPoint(P); } sprintf(name,"mark2d_%d",++NbMark); Handle(Draw_Marker2D) mark2d = new (Draw_Marker2D)(P2d, Draw_X, Draw_orange); Draw::Set(name, mark2d); } #endif } else if ((diffAng>0)&&(LinCont->Order()==1)) { diffAng=diffAng/myLinCont->Value(i)->G1Criterion(U); if (diffAng>diffAngMax) diffAngMax = diffAng; SdiffAng+=diffAng; NdiffAng++; #if DRAW if ((Affich) && (NbBoucle == myNbIter)) { gp_Pnt P; LinCont->D0(U,P); Handle(Draw_Marker3D) mark = new Draw_Marker3D(P, Draw_X, Draw_or); char* name = new char[100]; sprintf(name,"mark_%d",++NbMark); Draw::Set(name, mark); } #endif } } if (NdiffDist>0) {// au moins un point n'est pas acceptable en G0 Standard_Real Coef; if(LinCont->Order()== 0) Coef = 0.6*Log(diffDistMax+7.4); //7.4 correspond au calcul du coef mini = 1.2 donc e^1.2/0.6 else Coef = Log(diffDistMax+3.3); //3.3 correspond au calcul du coef mini = 1.2 donc e^1.2 if (Coef>3) Coef=3; //experimentalement apres le coef devient mauvais pour les cas L if ((NbBoucle>1)&&(diffDistMax>2)) { Coef=1.6; } if (LinCont->NbPoints()>=Floor(LinCont->NbPoints()*Coef)) Coef=2;// pour assurer une augmentation du nombre de points LinCont->SetNbPoints(Standard_Integer(LinCont->NbPoints() * Coef)); Result=Standard_False; } else if (NdiffAng>0) // au moins 1 points ne sont pas accepable en G1 { Standard_Real Coef=1.5; if ((LinCont->NbPoints()+1)>=Floor(LinCont->NbPoints()*Coef)) Coef=2; LinCont->SetNbPoints(Standard_Integer(LinCont->NbPoints()*Coef )) ; Result=Standard_False; } } } if (!Result) { if (myFree && NbBoucle == 1) myPrevPlate = myPlate; myPlate.Init(); } return Result; } //--------------------------------------------------------- //fonction : VerifPoint //--------------------------------------------------------- void GeomPlate_BuildPlateSurface:: VerifPoints (Standard_Real& Dist, Standard_Real& Ang, Standard_Real& Curv) const { Standard_Integer NTPntCont=myPntCont->Length(); gp_Pnt Pi, Pf; gp_Pnt2d P2d; gp_Vec v1i, v1f, v2i, v2f, v3i, v3f; Ang=0;Dist=0,Curv=0; Handle(GeomPlate_PointConstraint) PntCont; for (Standard_Integer i=1;i<=NTPntCont;i++) { PntCont = myPntCont->Value(i); switch (PntCont->Order()) { case 0 : P2d = PntCont->Pnt2dOnSurf(); PntCont->D0(Pi); myGeomPlateSurface->D0( P2d.Coord(1), P2d.Coord(2), Pf); Dist = Pf.Distance(Pi); break; case 1 : PntCont->D1(Pi, v1i, v2i); P2d = PntCont->Pnt2dOnSurf(); myGeomPlateSurface->D1( P2d.Coord(1), P2d.Coord(2), Pf, v1f, v2f); Dist = Pf.Distance(Pi); v3i = v1i^v2i; v3f=v1f^v2f; Ang=v3f.Angle(v3i); if (Ang>(PI/2)) Ang = PI -Ang; break; case 2 : Handle(Geom_Surface) Splate; Splate = Handle(Geom_Surface)::DownCast(myGeomPlateSurface); LocalAnalysis_SurfaceContinuity CG2; P2d = PntCont->Pnt2dOnSurf(); GeomLProp_SLProps Prop (Splate, P2d.Coord(1), P2d.Coord(2), 2, 0.001); CG2.ComputeAnalysis(Prop, PntCont->LPropSurf(),GeomAbs_G2); Dist=CG2.C0Value(); Ang=CG2.G1Angle(); Curv=CG2.G2CurvatureGap(); break; } } } Standard_Real GeomPlate_BuildPlateSurface::ComputeAnisotropie() const { if (myAnisotropie) { //Temporary return 1.0; } else return 1.0; } Standard_Boolean GeomPlate_BuildPlateSurface::IsOrderG1() const { Standard_Boolean result = Standard_True; for (Standard_Integer i = 1; i <= myLinCont->Length(); i++) if (myLinCont->Value(i)->Order() < 1) { result = Standard_False; break; } return result; }