// Created on: 1997-05-05 // Created by: Jerome LEMONIER // Copyright (c) 1996-1999 Matra Datavision // Copyright (c) 1999-2014 OPEN CASCADE SAS // // This file is part of Open CASCADE Technology software library. // // This library is free software; you can redistribute it and/or modify it under // the terms of the GNU Lesser General Public License version 2.1 as published // by the Free Software Foundation, with special exception defined in the file // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT // distribution for complete text of the license and disclaimer of any warranty. // // Alternatively, this file may be used under the terms of Open CASCADE // commercial license or contractual agreement. #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #ifdef OCCT_DEBUG #include static int Affich = 0; #endif //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// // ========================================================= // C O N S T R U C T O R S // ========================================================= //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// //--------------------------------------------------------- // Constructor compatible with the old version //--------------------------------------------------------- GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface( const occ::handle>& NPoints, const occ::handle>>& TabCurve, const occ::handle>& Tang, const int Degree, const int NbIter, const double Tol2d, const double Tol3d, const double TolAng, const double, const bool Anisotropie) : myAnisotropie(Anisotropie), myDegree(Degree), myNbIter(NbIter), myTol2d(Tol2d), myTol3d(Tol3d), myTolAng(TolAng), myNbBounds(0) { int NTCurve = TabCurve->Length(); // Number of linear constraints myNbPtsOnCur = 0; // Different calculation of the number of points depending on the length myLinCont = new NCollection_HSequence>; myPntCont = new NCollection_HSequence>; if (myNbIter < 1) { throw Standard_ConstructionError("GeomPlate : Number of iteration must be >= 1"); } if (NTCurve == 0) { throw Standard_ConstructionError("GeomPlate : the bounds Array is null"); } if (Tang->Length() == 0) { throw Standard_ConstructionError("GeomPlate : the constraints Array is null"); } int nbp = 0; int i; for (i = 1; i <= NTCurve; i++) { nbp += NPoints->Value(i); } if (nbp == 0) { throw Standard_ConstructionError( "GeomPlate : the resolution is impossible if the number of constraints points is 0"); } if (myDegree < 2) { throw Standard_ConstructionError("GeomPlate ; the degree resolution must be upper of 2"); } // Filling fields passing from the old constructor to the new one for (i = 1; i <= NTCurve; i++) { occ::handle Cont = new GeomPlate_CurveConstraint(TabCurve->Value(i), Tang->Value(i), NPoints->Value(i)); myLinCont->Append(Cont); } mySurfInitIsGive = false; myIsLinear = true; myFree = false; } //------------------------------------------------------------------ // Constructor with initial surface and degree //------------------------------------------------------------------ GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface(const occ::handle& Surf, const int Degree, const int NbPtsOnCur, const int NbIter, const double Tol2d, const double Tol3d, const double TolAng, const double /*TolCurv*/, const bool Anisotropie) : mySurfInit(Surf), myAnisotropie(Anisotropie), myDegree(Degree), myNbPtsOnCur(NbPtsOnCur), myNbIter(NbIter), myTol2d(Tol2d), myTol3d(Tol3d), myTolAng(TolAng), myNbBounds(0) { if (myNbIter < 1) { throw Standard_ConstructionError("GeomPlate : Number of iteration must be >= 1"); } if (myDegree < 2) { throw Standard_ConstructionError("GeomPlate : the degree must be above 2"); } myLinCont = new NCollection_HSequence>; myPntCont = new NCollection_HSequence>; mySurfInitIsGive = true; myIsLinear = true; myFree = false; } //--------------------------------------------------------- // Constructor with degree //--------------------------------------------------------- GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface(const int Degree, const int NbPtsOnCur, const int NbIter, const double Tol2d, const double Tol3d, const double TolAng, const double /*TolCurv*/, const bool Anisotropie) : myAnisotropie(Anisotropie), myDegree(Degree), myNbPtsOnCur(NbPtsOnCur), myNbIter(NbIter), myTol2d(Tol2d), myTol3d(Tol3d), myTolAng(TolAng), myNbBounds(0) { if (myNbIter < 1) { throw Standard_ConstructionError("GeomPlate : Number of iteration must be >= 1"); } if (myDegree < 2) { throw Standard_ConstructionError("GeomPlate : the degree resolution must be upper of 2"); } myLinCont = new NCollection_HSequence>; myPntCont = new NCollection_HSequence>; mySurfInitIsGive = false; myIsLinear = true; myFree = false; } //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// // ========================================================= // P U B L I C M E T H O D S // ========================================================= //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// //================================================================================================= // Reorder the table of transformations // After the call of CourbeJointive the order of curves is modified // Ex : initial order of curves ==> A B C D E F // In TabInit we note ==> 1 2 3 4 5 6 // after CourbeJointive ==> A E C B D F // TabInit ==> 1 5 3 2 4 6 // after TrierTab the Table contains ==> 1 4 3 5 2 6 // It is also possible to access the 2nd curve by taking TabInit[2] // i.e. the 4th from the table of classified curves //------------------------------------------------------------------------- static void TrierTab(occ::handle>& Tab) { // Parse the table of transformations to find the initial order int Nb = Tab->Length(); NCollection_Array1 TabTri(1, Nb); for (int i = 1; i <= Nb; i++) { TabTri.SetValue(Tab->Value(i), i); } Tab->ChangeArray1() = TabTri; } //--------------------------------------------------------- // Function : ProjectCurve //--------------------------------------------------------- occ::handle GeomPlate_BuildPlateSurface::ProjectCurve( const occ::handle& Curv) { // Project a curve on a plane occ::handle Curve2d; occ::handle hsur = new GeomAdaptor_Surface(mySurfInit); gp_Pnt2d P2d; occ::handle HProjector = new ProjLib_HCompProjectedCurve(hsur, Curv, myTol3d / 10, myTol3d / 10); double UdebCheck, UfinCheck, ProjUdeb, ProjUfin; UdebCheck = Curv->FirstParameter(); UfinCheck = Curv->LastParameter(); HProjector->Bounds(1, ProjUdeb, ProjUfin); if (HProjector->NbCurves() != 1 || std::abs(UdebCheck - ProjUdeb) > Precision::PConfusion() || std::abs(UfinCheck - ProjUfin) > Precision::PConfusion()) { if (HProjector->IsSinglePnt(1, P2d)) { // solution in a point NCollection_Array1 poles(1, 2); poles.Init(P2d); Curve2d = new (Geom2d_BezierCurve)(poles); } else { Curve2d.Nullify(); // No continuous solution #ifdef OCCT_DEBUG std::cout << "BuildPlateSurace :: No continuous projection" << std::endl; #endif } } else { GeomAbs_Shape Continuity = GeomAbs_C1; int MaxDegree = 10, MaxSeg; double Udeb, Ufin; HProjector->Bounds(1, Udeb, Ufin); MaxSeg = 20 + HProjector->NbIntervals(GeomAbs_C3); Approx_CurveOnSurface appr(HProjector, hsur, Udeb, Ufin, myTol3d); appr.Perform(MaxSeg, MaxDegree, Continuity, false, true); Curve2d = appr.Curve2d(); } return Curve2d; } //--------------------------------------------------------- // Function : ProjectedCurve //--------------------------------------------------------- occ::handle GeomPlate_BuildPlateSurface::ProjectedCurve( occ::handle& Curv) { // Projection of a curve on the initial surface occ::handle hsur = new GeomAdaptor_Surface(mySurfInit); occ::handle HProjector = new ProjLib_HCompProjectedCurve(hsur, Curv, myTolU / 10, myTolV / 10); if (HProjector->NbCurves() != 1) { HProjector.Nullify(); // No continuous solution #ifdef OCCT_DEBUG std::cout << "BuildPlateSurace :: No continuous projection" << std::endl; #endif } else { double First1, Last1, First2, Last2; First1 = Curv->FirstParameter(); Last1 = Curv->LastParameter(); HProjector->Bounds(1, First2, Last2); if (std::abs(First1 - First2) <= std::max(myTolU, myTolV) && std::abs(Last1 - Last2) <= std::max(myTolU, myTolV)) { HProjector = occ::down_cast( HProjector->Trim(First2, Last2, Precision::PConfusion())); } else { HProjector.Nullify(); // No continuous solution #ifdef OCCT_DEBUG std::cout << "BuildPlateSurace :: No complete projection" << std::endl; #endif } } return HProjector; } //--------------------------------------------------------- // Function : ProjectPoint //--------------------------------------------------------- // Projects a point on the initial surface //--------------------------------------------------------- gp_Pnt2d GeomPlate_BuildPlateSurface::ProjectPoint(const gp_Pnt& p3d) { Extrema_POnSurf P; myProj.Perform(p3d); int nearest = 1; if (myProj.NbExt() > 1) { double dist2mini = myProj.SquareDistance(1); for (int i = 2; i <= myProj.NbExt(); i++) { if (myProj.SquareDistance(i) < dist2mini) { dist2mini = myProj.SquareDistance(i); nearest = i; } } } P = myProj.Point(nearest); double u, v; P.Parameter(u, v); gp_Pnt2d p2d; p2d.SetCoord(u, v); return p2d; } //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// // ========================================================= // P U B L I C M E T H O D S // ========================================================= //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// //--------------------------------------------------------- // Function : Init //--------------------------------------------------------- // Initializes linear and point constraints //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::Init() { myLinCont->Clear(); myPntCont->Clear(); myPntCont = new NCollection_HSequence>; myLinCont = new NCollection_HSequence>; } //--------------------------------------------------------- // Function : LoadInitSurface //--------------------------------------------------------- // Loads the initial surface //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::LoadInitSurface(const occ::handle& Surf) { mySurfInit = Surf; mySurfInitIsGive = true; } //--------------------------------------------------------- // Function : Add //--------------------------------------------------------- //--------------------------------------------------------- // Adds a linear constraint //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::Add(const occ::handle& Cont) { myLinCont->Append(Cont); } void GeomPlate_BuildPlateSurface::SetNbBounds(const int NbBounds) { myNbBounds = NbBounds; } //--------------------------------------------------------- // Function : Add //--------------------------------------------------------- //--------------------------------------------------------- // Adds a point constraint //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::Add(const occ::handle& Cont) { myPntCont->Append(Cont); } //--------------------------------------------------------- // Function : Perform // Calculates the surface filled with loaded constraints //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::Perform(const Message_ProgressRange& theProgress) { myGeomPlateSurface.Nullify(); #ifdef OCCT_DEBUG // Timing OSD_Chronometer Chrono; Chrono.Reset(); Chrono.Start(); #endif if (myNbBounds == 0) { myNbBounds = myLinCont->Length(); } myPlate.Init(); //===================================================================== // Declaration of variables. //===================================================================== int NTLinCont = myLinCont->Length(), NTPntCont = myPntCont->Length(), NbBoucle = 0; bool Fini = true; if ((NTLinCont + NTPntCont) == 0) { #ifdef OCCT_DEBUG std::cout << "WARNING : GeomPlate : The number of constraints is null." << std::endl; #endif return; } //====================================================================== // Initial Surface //====================================================================== Message_ProgressScope aPS(theProgress, "Calculating the surface filled", 100, true); if (!mySurfInitIsGive) { ComputeSurfInit(aPS.Next(10)); if (aPS.UserBreak()) { return; } } else { if (NTLinCont >= 2) { // Table of transformations to preserve the initial order, see TrierTab myInitOrder = new NCollection_HArray1(1, NTLinCont); for (int l = 1; l <= NTLinCont; l++) { myInitOrder->SetValue(l, l); } if (!CourbeJointive(myTol3d)) { // throw Standard_Failure("Curves are not joined"); #ifdef OCCT_DEBUG std::cout << "WARNING : Courbes non jointives a " << myTol3d << " pres" << std::endl; #endif } TrierTab(myInitOrder); // Reorder the table of transformations } else if (NTLinCont > 0) // Patch { mySense = new NCollection_HArray1(1, NTLinCont, 0); myInitOrder = new NCollection_HArray1(1, NTLinCont, 1); } } if (mySurfInit.IsNull()) { return; } double u1, v1, u2, v2; mySurfInit->Bounds(u1, v1, u2, v2); GeomAdaptor_Surface aSurfInit(mySurfInit); myTolU = aSurfInit.UResolution(myTol3d); myTolV = aSurfInit.VResolution(myTol3d); myProj.Initialize(aSurfInit, u1, v1, u2, v2, myTolU, myTolV); //====================================================================== // Projection of curves //====================================================================== bool Ok = true; for (int i = 1; i <= NTLinCont; i++) { if (myLinCont->Value(i)->Curve2dOnSurf().IsNull()) { occ::handle Curve2d = ProjectCurve(myLinCont->Value(i)->Curve3d()); if (Curve2d.IsNull()) { Ok = 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 = true; for (int i = 1; i <= NTLinCont; i++) { occ::handle Curve2d = ProjectCurve(myLinCont->Value(i)->Curve3d()); if (Curve2d.IsNull()) { Ok = false; break; } myLinCont->ChangeValue(i)->SetCurve2dOnSurf(Curve2d); } if (!Ok) { mySurfInit = myPlanarSurfInit; mySurfInit->Bounds(u1, v1, u2, v2); GeomAdaptor_Surface SurfNew(mySurfInit); myTolU = SurfNew.UResolution(myTol3d); myTolV = SurfNew.VResolution(myTol3d); myProj.Initialize(SurfNew, u1, v1, u2, v2, myTolU, myTolV); for (int i = 1; i <= NTLinCont; i++) { myLinCont->ChangeValue(i)->SetCurve2dOnSurf(ProjectCurve(myLinCont->Value(i)->Curve3d())); } } else { // Project the points for (int i = 1; i <= NTPntCont; i++) { gp_Pnt P; myPntCont->Value(i)->D0(P); myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P)); } } } //====================================================================== // Projection of points //====================================================================== for (int i = 1; i <= NTPntCont; i++) { if (!myPntCont->Value(i)->HasPnt2dOnSurf()) { gp_Pnt P; myPntCont->Value(i)->D0(P); myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P)); } } //====================================================================== // Number of points by curve //====================================================================== if ((NTLinCont != 0) && (myNbPtsOnCur != 0)) { CalculNbPtsInit(); } //====================================================================== // Management of incompatibilites between curves //====================================================================== occ::handle>> PntInter; occ::handle>> PntG1G1; if (NTLinCont != 0) { PntInter = new NCollection_HArray1>(1, NTLinCont); PntG1G1 = new NCollection_HArray1>(1, NTLinCont); Intersect(PntInter, PntG1G1); } //====================================================================== // Loop to obtain a better surface //====================================================================== myFree = !myIsLinear; do { #ifdef OCCT_DEBUG if (Affich && NbBoucle) { std::cout << "Resultats boucle" << NbBoucle << std::endl; std::cout << "DistMax=" << myG0Error << std::endl; if (myG1Error != 0) std::cout << "AngleMax=" << myG1Error << std::endl; if (myG2Error != 0) std::cout << "CourbMax=" << myG2Error << std::endl; } #endif NbBoucle++; if (NTLinCont != 0) { //==================================================================== // Calculate the total number of points and the maximum of points by curve //==================================================================== int NPointMax = 0; for (int i = 1; i <= NTLinCont; i++) { if ((myLinCont->Value(i)->NbPoints()) > NPointMax) { NPointMax = (int)(myLinCont->Value(i)->NbPoints()); } } //==================================================================== // Discretization of curves //==================================================================== Discretise(PntInter, PntG1G1); //==================================================================== // Preparation of constraint points for plate //==================================================================== LoadCurve(NbBoucle); if (myPntCont->Length() != 0) { LoadPoint(NbBoucle); } //==================================================================== // Construction of the surface //==================================================================== myPlate.SolveTI(myDegree, ComputeAnisotropie(), aPS.Next(90)); if (aPS.UserBreak()) { return; } if (!myPlate.IsDone()) { #ifdef OCCT_DEBUG std::cout << "WARNING : GeomPlate : calculation of Plate failed" << std::endl; #endif return; } myGeomPlateSurface = new GeomPlate_Surface(mySurfInit, myPlate); double Umin, Umax, Vmin, Vmax; myPlate.UVBox(Umin, Umax, Vmin, Vmax); myGeomPlateSurface->SetBounds(Umin, Umax, Vmin, Vmax); Fini = VerifSurface(NbBoucle); if ((NbBoucle >= myNbIter) && (!Fini)) { #ifdef OCCT_DEBUG std::cout << "Warning: objective was not reached" << std::endl; #endif Fini = true; } if ((NTPntCont != 0) && (Fini)) { double di, an, cu; VerifPoints(di, an, cu); } } else { LoadPoint(NbBoucle); //==================================================================== // Construction of the surface //==================================================================== myPlate.SolveTI(myDegree, ComputeAnisotropie(), aPS.Next(90)); if (aPS.UserBreak()) { return; } if (!myPlate.IsDone()) { #ifdef OCCT_DEBUG std::cout << "WARNING : GeomPlate : calculation of Plate failed" << std::endl; #endif return; } myGeomPlateSurface = new GeomPlate_Surface(mySurfInit, myPlate); double Umin, Umax, Vmin, Vmax; myPlate.UVBox(Umin, Umax, Vmin, Vmax); myGeomPlateSurface->SetBounds(Umin, Umax, Vmin, Vmax); Fini = true; double di, an, cu; VerifPoints(di, an, cu); } } while (!Fini); // End loop for better surface #ifdef OCCT_DEBUG if (NTLinCont != 0) { std::cout << "======== Global results ===========" << std::endl; std::cout << "DistMax=" << myG0Error << std::endl; if (myG1Error != 0) std::cout << "AngleMax=" << myG1Error << std::endl; if (myG2Error != 0) std::cout << "CourbMax=" << myG2Error << std::endl; } Chrono.Stop(); double Tps; Chrono.Show(Tps); std::cout << "*** END OF GEOMPLATE ***" << std::endl; std::cout << "Time of calculation : " << Tps << std::endl; std::cout << "Number of loops : " << NbBoucle << std::endl; #endif } //--------------------------------------------------------- // Function : EcartContraintesMIL //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::EcartContraintesMil( const int c, occ::handle>& d, occ::handle>& an, occ::handle>& courb) { int NbPt = myParCont->Value(c).Length(); double 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; int i; occ::handle LinCont = myLinCont->Value(c); switch (LinCont->Order()) { case 0: for (i = 1; i < NbPt; i++) { U = (myParCont->Value(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; i < NbPt; i++) { U = (myParCont->Value(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; double angle = v3f.Angle(v3i); if (angle > (M_PI / 2)) { an->ChangeValue(i) = M_PI - angle; } else { an->ChangeValue(i) = angle; } courb->Init(0); } break; case 2: occ::handle Splate(myGeomPlateSurface); LocalAnalysis_SurfaceContinuity CG2; for (i = 1; i < NbPt; i++) { U = (myParCont->Value(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; } } //--------------------------------------------------------- // Function : Disc2dContour //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::Disc2dContour(const int /*nbp*/, NCollection_Sequence& Seq2d) { #ifdef OCCT_DEBUG if (Seq2d.Length() != 4) std::cout << "Number of points should be equal to 4 for Disc2dContour" << std::endl; #endif // initialization Seq2d.Clear(); // sampling in "cosine" + 3 points on each interval int NTCurve = myLinCont->Length(); int NTPntCont = myPntCont->Length(); gp_Pnt2d P2d; gp_XY UV; gp_Pnt PP; double u1, v1, u2, v2; int 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++) { occ::handle LinCont = myLinCont->Value(i); if (LinCont->Order() != -1) { int NbPt = myParCont->Value(i).Length(); // first point of constraint (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 (int j = 2; j < NbPt; j++) { double Uj = myParCont->Value(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 1/2 previous 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 previous 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); // current constraint point 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); } } } } //--------------------------------------------------------- // Function : Disc3dContour //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::Disc3dContour(const int /*nbp*/, const int iordre, NCollection_Sequence& Seq3d) { #ifdef OCCT_DEBUG if (Seq3d.Length() != 4) std::cout << "nbp should be equal to 4 for Disc3dContour" << std::endl; if (iordre != 0 && iordre != 1) std::cout << "incorrect order for Disc3dContour" << std::endl; #endif // initialization Seq3d.Clear(); // sampling in "cosine" + 3 points on each interval double u1, v1, u2, v2; mySurfInit->Bounds(u1, v1, u2, v2); GeomAdaptor_Surface Surf(mySurfInit); myProj.Initialize(Surf, u1, v1, u2, v2, Surf.UResolution(myTol3d), Surf.VResolution(myTol3d)); int NTCurve = myLinCont->Length(); int NTPntCont = myPntCont->Length(); // gp_Pnt2d P2d; gp_Pnt P3d; gp_Vec v1h, v2h, v3h; gp_XYZ Pos; int 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) { int NbPt = myParCont->Value(i).Length(); // first constraint point (j=0) // int 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 (int j = 2; j < NbPt; j++) { double Uj = myParCont->Value(i).Value(j), Ujp1 = myParCont->Value(i).Value(j + 1); if (iordre == 0) { // point 1/4 previous 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 1/2 previous 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 previous myLinCont->Value(i)->D0((3 * Ujp1 + Uj) / 4, P3d); Pos.SetX(P3d.X()); Pos.SetY(P3d.Y()); Pos.SetZ(P3d.Z()); Seq3d.Append(Pos); // current constraint point 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 previous 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 1/2 previous 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 previous 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); // current constraint point 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); } } } } } //--------------------------------------------------------- // Function : IsDone //--------------------------------------------------------- bool GeomPlate_BuildPlateSurface::IsDone() const { return myPlate.IsDone(); } //--------------------------------------------------------- // Function : Surface //--------------------------------------------------------- occ::handle GeomPlate_BuildPlateSurface::Surface() const { return myGeomPlateSurface; } //--------------------------------------------------------- // Function : SurfInit //--------------------------------------------------------- occ::handle GeomPlate_BuildPlateSurface::SurfInit() const { return mySurfInit; } //--------------------------------------------------------- // Function : Sense //--------------------------------------------------------- occ::handle> GeomPlate_BuildPlateSurface::Sense() const { int NTCurve = myLinCont->Length(); occ::handle> Sens = new NCollection_HArray1(1, NTCurve); for (int i = 1; i <= NTCurve; i++) { Sens->SetValue(i, mySense->Value(myInitOrder->Value(i))); } return Sens; } //--------------------------------------------------------- // Function : Curve2d //--------------------------------------------------------- occ::handle>> GeomPlate_BuildPlateSurface::Curves2d() const { int NTCurve = myLinCont->Length(); occ::handle>> C2dfin = new NCollection_HArray1>(1, NTCurve); for (int i = 1; i <= NTCurve; i++) { C2dfin->SetValue(i, myLinCont->Value(myInitOrder->Value(i))->Curve2dOnSurf()); } return C2dfin; } //--------------------------------------------------------- // Function : Order //--------------------------------------------------------- occ::handle> GeomPlate_BuildPlateSurface::Order() const { occ::handle> result = new NCollection_HArray1(1, myLinCont->Length()); for (int i = 1; i <= myLinCont->Length(); i++) { result->SetValue(myInitOrder->Value(i), i); } return result; } //--------------------------------------------------------- // Function : G0Error //--------------------------------------------------------- double GeomPlate_BuildPlateSurface::G0Error() const { return myG0Error; } //--------------------------------------------------------- // Function : G1Error //--------------------------------------------------------- double GeomPlate_BuildPlateSurface::G1Error() const { return myG1Error; } //--------------------------------------------------------- // Function : G2Error //--------------------------------------------------------- double GeomPlate_BuildPlateSurface::G2Error() const { return myG2Error; } //================================================================================================= double GeomPlate_BuildPlateSurface::G0Error(const int Index) { occ::handle> tdistance = new NCollection_HArray1(1, myNbPtsOnCur); occ::handle> tangle = new NCollection_HArray1(1, myNbPtsOnCur); occ::handle> tcurvature = new NCollection_HArray1(1, myNbPtsOnCur); EcartContraintesMil(Index, tdistance, tangle, tcurvature); double MaxDistance = 0.; for (int i = 1; i <= myNbPtsOnCur; i++) { if (tdistance->Value(i) > MaxDistance) { MaxDistance = tdistance->Value(i); } } return MaxDistance; } //================================================================================================= double GeomPlate_BuildPlateSurface::G1Error(const int Index) { occ::handle> tdistance = new NCollection_HArray1(1, myNbPtsOnCur); occ::handle> tangle = new NCollection_HArray1(1, myNbPtsOnCur); occ::handle> tcurvature = new NCollection_HArray1(1, myNbPtsOnCur); EcartContraintesMil(Index, tdistance, tangle, tcurvature); double MaxAngle = 0.; for (int i = 1; i <= myNbPtsOnCur; i++) { if (tangle->Value(i) > MaxAngle) { MaxAngle = tangle->Value(i); } } return MaxAngle; } //================================================================================================= double GeomPlate_BuildPlateSurface::G2Error(const int Index) { occ::handle> tdistance = new NCollection_HArray1(1, myNbPtsOnCur); occ::handle> tangle = new NCollection_HArray1(1, myNbPtsOnCur); occ::handle> tcurvature = new NCollection_HArray1(1, myNbPtsOnCur); EcartContraintesMil(Index, tdistance, tangle, tcurvature); double MaxCurvature = 0.; for (int i = 1; i <= myNbPtsOnCur; i++) { if (tcurvature->Value(i) > MaxCurvature) { MaxCurvature = tcurvature->Value(i); } } return MaxCurvature; } occ::handle GeomPlate_BuildPlateSurface::CurveConstraint( const int order) const { return myLinCont->Value(order); } occ::handle GeomPlate_BuildPlateSurface::PointConstraint( const int order) const { return myPntCont->Value(order); } //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// // ========================================================= // P R I V A T E M E T H O D S // ========================================================= //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\// //======================================================================= // Function : CourbeJointive // Purpose : Create a chain of curves to calculate the // initial surface with the method of max flow. // Return true if it is a closed contour. //======================================================================= bool GeomPlate_BuildPlateSurface::CourbeJointive(const double tolerance) { int nbf = myLinCont->Length(); double Ufinal1, Uinit1, Ufinal2, Uinit2; mySense = new NCollection_HArray1(1, nbf, 0); bool result = true; int j = 1, i; gp_Pnt P1, P2; while (j <= (myNbBounds - 1)) { int a = 0; i = j + 1; if (i > myNbBounds) { result = false; a = 2; } while (a < 1) { if (i > myNbBounds) { result = 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) < tolerance) { if (i != j + 1) { occ::handle tampon = myLinCont->Value(j + 1); myLinCont->SetValue(j + 1, myLinCont->Value(i)); myLinCont->SetValue(i, tampon); int Tmp = myInitOrder->Value(j + 1); // See function TrierTab for the functioning of 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) < tolerance) { if (i != j + 1) { occ::handle tampon = myLinCont->Value(j + 1); myLinCont->SetValue(j + 1, myLinCont->Value(i)); myLinCont->SetValue(i, tampon); int Tmp = myInitOrder->Value(j + 1); // See function TrierTab for the functioning of 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) < tolerance)) { return ((result)); } myLinCont->Value(myNbBounds)->D0(Uinit1, P1); if ((mySense->Value(myNbBounds) == 1) && (P1.Distance(P2) < tolerance)) { return ((result)); } else { return false; } } //================================================================================================= // Calculate the initial surface either by the method of max flow or by // the method of the plane of inertia if the contour is not closed or if // there are point constraints. //------------------------------------------------------------------------- void GeomPlate_BuildPlateSurface::ComputeSurfInit(const Message_ProgressRange& theProgress) { int nopt = 2, popt = 2, Np = 1; bool isHalfSpace = true; double LinTol = 0.001, AngTol = 0.001; // AngTol = 0.0001; //LinTol = 0.0001 // Option to calculate the initial plane int NTLinCont = myLinCont->Length(), NTPntCont = myPntCont->Length(); // Table of transformation to preserve the initial order see TrierTab if (NTLinCont != 0) { myInitOrder = new NCollection_HArray1(1, NTLinCont); for (int i = 1; i <= NTLinCont; i++) { myInitOrder->SetValue(i, i); } } bool CourbeJoint = (NTLinCont != 0) && CourbeJointive(myTol3d); if (CourbeJoint && IsOrderG1()) { nopt = 3; // Table contains the cloud of points for calculation of the plane int NbPoint = 20, Discr = NbPoint / 4, pnum = 0; occ::handle> Pts = new NCollection_HArray1(1, (NbPoint + 1) * NTLinCont + NTPntCont); NCollection_Sequence Vecs, NewVecs; NCollection_Sequence Aset; double Uinit, Ufinal, Uif; gp_Vec LastVec; int i; for (i = 1; i <= NTLinCont; i++) { int 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; bool ToReverse = 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 = true; } } for (int j = 0; j <= NbPoint; j++) { // Number of points per curve = 20 // Linear distribution double 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(); } bool isNew = true; int k; for (k = 1; k <= Vecs.Length(); k++) { if (Vecs(k).IsEqual(Normal, LinTol, AngTol)) { isNew = false; break; } } if (isNew) { for (k = 1; k <= NewVecs.Length(); k++) { if (NewVecs(k).IsEqual(Normal, LinTol, AngTol)) { isNew = 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++) { int 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(); bool isNew = true; for (int k = 1; k <= Vecs.Length(); k++) { if (Vecs(k).IsEqual(Normal, LinTol, AngTol)) { isNew = 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) { bool NullExist = true; while (NullExist) { NullExist = false; for (i = 1; i <= Vecs.Length(); i++) { if (Vecs(i).SquareMagnitude() == 0.) { NullExist = true; Vecs.Remove(i); break; } } } GeomPlate_BuildAveragePlane BAP(Vecs, Pts); double u1, u2, v1, v2; BAP.MinMaxBox(u1, u2, v1, v2); // The space is greater for projections double du = u2 - u1; double 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 OCCT_DEBUG std::cout << std::endl << "Normals are not in half space" << std::endl << std::endl; #endif myIsLinear = false; nopt = 2; } } // if (NTLinCont != 0 && (CourbeJoint = CourbeJointive( myTol3d )) && IsOrderG1()) if (NTLinCont != 0) { TrierTab(myInitOrder); // Reorder the table of transformations } if (nopt != 3) { if (NTPntCont != 0) { nopt = 1; // Calculate by the method of plane of inertia } else if (!CourbeJoint || NTLinCont != myNbBounds) { // throw Standard_Failure("Curves are not joined"); #ifdef OCCT_DEBUG std::cout << "WARNING : Curves are non-adjacent with tolerance " << myTol3d << std::endl; #endif nopt = 1; } double LenT = 0; int Npt = 0; int NTPoint = 20 * NTLinCont; int i; for (i = 1; i <= NTLinCont; i++) { LenT += myLinCont->Value(i)->Length(); } for (i = 1; i <= NTLinCont; i++) { int NbPoint = (int)(NTPoint * (myLinCont->Value(i)->Length()) / LenT); if (NbPoint < 10) { NbPoint = 10; } (void)Npt; // unused but set for debug Npt += NbPoint; } // Table containing a cloud of points for calculation of the plane occ::handle> Pts = new NCollection_HArray1(1, 20 * NTLinCont + NTPntCont); int NbPoint = 20; double 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 (int j = 0; j < NbPoint; j++) { // Number of points per curve = 20 // Linear distribution double Inter = j * Uif / (NbPoint); gp_Pnt P; myLinCont->Value(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()) { #ifdef OCCT_DEBUG std::cout << "WARNING : GeomPlate : the initial surface is not a plane." << std::endl; #endif return; } double u1, u2, v1, v2; BAP.MinMaxBox(u1, u2, v1, v2); // The space is greater for projections double du = u2 - u1; double 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) { occ::handle InitPlane = (occ::down_cast(mySurfInit))->BasisSurface(); // clang-format off double Ratio = 0., R1 = 2., R2 = 0.6; //R1 = 3, R2 = 0.5;//R1 = 1.4, R2 = 0.8; //R1 = 5., R2 = 0.2; // clang-format on occ::handle hsur = new GeomAdaptor_Surface(InitPlane); int NbPoint = 20; // gp_Pnt P; // gp_Vec DerC, DerCproj, DU, DV; // gp_Pnt2d P2d; // gp_Vec2d DProj; for (int i = 1; i <= NTLinCont && myIsLinear; i++) { double FirstPar = myLinCont->Value(i)->FirstParameter(); double LastPar = myLinCont->Value(i)->LastParameter(); double Uif = (LastPar - FirstPar) / (NbPoint); occ::handle Curve = myLinCont->Value(i)->Curve3d(); occ::handle ProjCurve = new ProjLib_HCompProjectedCurve(hsur, Curve, myTol3d, myTol3d); Adaptor3d_CurveOnSurface AProj(ProjCurve, hsur); gp_Pnt P; gp_Vec DerC, DerCproj; for (int j = 1; j < NbPoint && myIsLinear; j++) { double Inter = FirstPar + j * Uif; Curve->D1(Inter, P, DerC); AProj.D1(Inter, P, DerCproj); double A1 = DerC.Magnitude(); double A2 = DerCproj.Magnitude(); if (A2 <= 1.e-20) { Ratio = 1.e20; } else { Ratio = A1 / A2; } if (Ratio > R1 || Ratio < R2) { myIsLinear = false; break; } } } #ifdef OCCT_DEBUG if (!myIsLinear) std::cout << "Metrics are too different :" << Ratio << std::endl; #endif // myIsLinear = true; // !! } // comparing metrics of curves and projected curves if (!myIsLinear) { myPlanarSurfInit = mySurfInit; double 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 of curves //====================================================================== int i; for (i = 1; i <= NTLinCont; i++) { if (myLinCont->Value(i)->Curve2dOnSurf().IsNull()) { myLinCont->ChangeValue(i)->SetCurve2dOnSurf(ProjectCurve(myLinCont->Value(i)->Curve3d())); } } //====================================================================== // Projection of 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)); } } //====================================================================== // Number of points by curve //====================================================================== if ((NTLinCont != 0) && (myNbPtsOnCur != 0)) { CalculNbPtsInit(); } //====================================================================== // Management of incompatibilities between curves //====================================================================== occ::handle>> PntInter; occ::handle>> PntG1G1; if (NTLinCont != 0) { PntInter = new NCollection_HArray1>(1, NTLinCont); PntG1G1 = new NCollection_HArray1>(1, NTLinCont); Intersect(PntInter, PntG1G1); } //==================================================================== // Discretization of curves //==================================================================== Discretise(PntInter, PntG1G1); //==================================================================== // Preparation of points of constraint for plate //==================================================================== LoadCurve(0, 0); if (myPntCont->Length() != 0) { LoadPoint(0, 0); } //==================================================================== // Construction of the surface //==================================================================== Message_ProgressScope aPS(theProgress, "ComputeSurfInit", 1); myPlate.SolveTI(2, ComputeAnisotropie(), aPS.Next()); if (theProgress.UserBreak()) { return; } if (!myPlate.IsDone()) { #ifdef OCCT_DEBUG std::cout << "WARNING : GeomPlate : calculation of Plate failed" << std::endl; #endif return; } myGeomPlateSurface = new GeomPlate_Surface(mySurfInit, myPlate); GeomPlate_MakeApprox App(myGeomPlateSurface, myTol3d, 1, 3, 15 * myTol3d, -1, GeomAbs_C0); mySurfInit = App.Surface(); mySurfInitIsGive = true; myPlate.Init(); // Reset for (i = 1; i <= NTLinCont; i++) { occ::handle NullCurve; NullCurve.Nullify(); myLinCont->ChangeValue(i)->SetCurve2dOnSurf(NullCurve); } } } //--------------------------------------------------------- // Function : Intersect //--------------------------------------------------------- // Find intersections between 2d curves // If the intersection is compatible (in cases G1-G1) // remove the point on one of two curves //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::Intersect( occ::handle>>& PntInter, occ::handle>>& PntG1G1) { int 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 (int i = 1; i <= NTLinCont; i++) { // double NbPnt_i=myLinCont->Value(i)->NbPoints(); // Find the intersection with each curve including the curve itself Ci.Load(myLinCont->Value(i)->Curve2dOnSurf()); for (int 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()) { // there is one intersection int nbpt = Intersection.NbPoints(); // number of points of intersection for (int k = 1; k <= nbpt; k++) { int2d = Intersection.Point(k); myLinCont->Value(i)->D0(int2d.ParamOnFirst(), P1); myLinCont->Value(j)->D0(int2d.ParamOnSecond(), P2); #ifdef OCCT_DEBUG if (Affich > 1) { std::cout << " Intersection " << k << " entre " << i << " &" << j << std::endl; std::cout << " Distance = " << P1.Distance(P2) << std::endl; } #endif if (P1.Distance(P2) < myTol3d) { // 2D intersection corresponds to close 3D points. // Note the interval, in which the point needs to be removed // to avoid duplications, which cause // error in plate. The point on curve i is removed; // the point on curve j is preserved; // the length of interval is a length 2d // corresponding in 3d to myTol3d double tolint = Ci.Resolution(myTol3d); Ci.D1(int2d.ParamOnFirst(), P2d, V2d); double 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); // If 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; double ant = v16.Angle(v26); if (ant > (M_PI / 2)) { ant = M_PI - ant; } if ((std::abs(v16 * v15 - v16 * v25) > (myTol3d / 1000)) || (std::abs(ant) > myTol3d / 1000)) // Non-compatible ==> remove zone in constraint G1 // corresponding to 3D tolerance of 0.01 { double coin; double Tol = 100 * myTol3d; double A1; gp_Pnt2d P1temp, P2temp; gp_Vec2d V1, V2; myLinCont->Value(i)->Curve2dOnSurf()->D1(int2d.ParamOnFirst(), P1temp, V1); myLinCont->Value(j)->Curve2dOnSurf()->D1(int2d.ParamOnSecond(), P2temp, V2); A1 = V1.Angle(V2); if (A1 > (M_PI / 2)) { A1 = M_PI - A1; } if (std::abs(std::abs(A1) - M_PI) < myTolAng) { Tol = 100000 * myTol3d; } #ifdef OCCT_DEBUG if (Affich) std::cout << "Angle between curves " << i << "," << j << " " << std::abs(std::abs(A1) - M_PI) << std::endl; #endif coin = Ci.Resolution(Tol); double Par1 = int2d.ParamOnFirst() - coin, Par2 = int2d.ParamOnFirst() + coin; // Storage of the interval for curve i PntG1G1->ChangeValue(i).Append(Par1); PntG1G1->ChangeValue(i).Append(Par2); coin = Cj.Resolution(Tol); Par1 = int2d.ParamOnSecond() - coin; Par2 = int2d.ParamOnSecond() + coin; // Storage of the interval for curve j PntG1G1->ChangeValue(j).Append(Par1); PntG1G1->ChangeValue(j).Append(Par2); } } // If 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) { occ::handle theCurve = myLinCont->Value(i)->Curve3d(); theCurve->D1(int2d.ParamOnFirst(), P1, vec); myLinCont->Value(j)->D1(int2d.ParamOnSecond(), P2, vecU, vecV); } else { occ::handle theCurve = myLinCont->Value(j)->Curve3d(); theCurve->D1(int2d.ParamOnSecond(), P2, vec); myLinCont->Value(i)->D1(int2d.ParamOnFirst(), P1, vecU, vecV); } N = vecU ^ vecV; double Angle = vec.Angle(N); Angle = std::abs(M_PI / 2 - Angle); if (Angle > myTolAng / 10.) //????????? //if (std::abs( scal ) > myTol3d/100) { // Non-compatible ==> one removes zone in constraint G0 and G1 // corresponding to 3D tolerance of 0.01 double coin; double Tol = 100 * myTol3d; double A1; gp_Pnt2d P1temp, P2temp; gp_Vec2d V1, V2; myLinCont->Value(i)->Curve2dOnSurf()->D1(int2d.ParamOnFirst(), P1temp, V1); myLinCont->Value(j)->Curve2dOnSurf()->D1(int2d.ParamOnSecond(), P2temp, V2); A1 = V1.Angle(V2); if (A1 > M_PI / 2) { A1 = M_PI - A1; } if (std::abs(std::abs(A1) - M_PI) < myTolAng) { Tol = 100000 * myTol3d; } #ifdef OCCT_DEBUG if (Affich) std::cout << "Angle entre Courbe " << i << "," << j << " " << std::abs(std::abs(A1) - M_PI) << std::endl; #endif if (myLinCont->Value(i)->Order() == 1) { coin = Ci.Resolution(Tol); coin *= Angle / myTolAng * 10.; #ifdef OCCT_DEBUG std::cout << std::endl << "coin = " << coin << std::endl; #endif double Par1 = int2d.ParamOnFirst() - coin; double Par2 = int2d.ParamOnFirst() + coin; // Storage of the interval for curve i PntG1G1->ChangeValue(i).Append(Par1); PntG1G1->ChangeValue(i).Append(Par2); } else { coin = Cj.Resolution(Tol); coin *= Angle / myTolAng * 10.; #ifdef OCCT_DEBUG std::cout << std::endl << "coin = " << coin << std::endl; #endif double Par1 = int2d.ParamOnSecond() - coin; double Par2 = int2d.ParamOnSecond() + coin; // Storage of the interval for curve j PntG1G1->ChangeValue(j).Append(Par1); PntG1G1->ChangeValue(j).Append(Par2); } } } } // if (P1.Distance( P2 ) < myTol3d) else { // 2D intersection corresponds to extended 3D points. // Note the interval where it is necessary to remove // the points to avoid duplications causing // error in plate. The point on curve i is removed, // the point on curve j is preserved. // The length of interval is 2D length // corresponding to the distance of points in 3D to myTol3d double 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 OCCT_DEBUG std::cout << "Attention: Two points 3d have the same projection dist = " << Dist << std::endl; #endif } } } } } } //--------------------------------------------------------- // Function : Discretize //--------------------------------------------------------- // Discretize curves according to parameters // the table of sequences Parcont contains all // parameter of points on curves // Field myPlateCont contains parameter of points on a plate; // it excludes duplicate points and incompatible zones. // The first part corresponds to verification of compatibility // and to removal of duplicate points. //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::Discretise( const occ::handle>>& PntInter, const occ::handle>>& PntG1G1) { int NTLinCont = myLinCont->Length(); bool ACR; occ::handle C2d; Geom2dAdaptor_Curve AC2d; // occ::handle HC2d; occ::handle acrlaw = new (Law_Interpol)(); myPlateCont = new NCollection_HArray1>(1, NTLinCont); myParCont = new NCollection_HArray1>(1, NTLinCont); //=========================================================================== // Construction of the table containing parameters of constraint points //=========================================================================== double Uinit, Ufinal, Length2d = 0, Inter; double CurLength; int NbPnt_i, NbPtInter, NbPtG1G1; occ::handle LinCont; for (int 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) { // Construct a law close to curvilinear abscissa if (!C2d.IsNull()) { AC2d.Load(C2d); } // AC2d.Load(LinCont->Curve2dOnSurf()); int ii, Nbint = 20; double U; NCollection_Array1 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 - std::cos((ii - 1) * M_PI / (Nbint))) / 2); tabP2d(ii).SetY(U); /* if (!HC2d.IsNull()) { double 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 = (int)(LinCont->NbPoints()); NbPtInter = PntInter->Value(i).Length(); NbPtG1G1 = PntG1G1->Value(i).Length(); #ifdef OCCT_DEBUG if (Affich > 1) { std::cout << "Courbe : " << i << std::endl; std::cout << " NbPnt, NbPtInter, NbPtG1G1 :" << NbPnt_i << ", " << NbPtInter << ", " << NbPtG1G1 << std::endl; } #endif for (int j = 1; j <= NbPnt_i; j++) { // Distribution of points in cosine following ACR 2D // To avoid points of accumulation in 2D // Inter=Uinit+(Uif)*((-cos(M_PI*((j-1)/(NbPnt_i-1)))+1)/2); if (j == NbPnt_i) { Inter = Ufinal; // to avoid bug on Sun } else if (ACR) { CurLength = Length2d * (1 - std::cos((j - 1) * M_PI / (NbPnt_i - 1))) / 2; Inter = acrlaw->Value(CurLength); } else { Inter = Uinit + (Ufinal - Uinit) * ((1 - std::cos((j - 1) * M_PI / (NbPnt_i - 1))) / 2); } myParCont->ChangeValue(i).Append(Inter); // add a point if (NbPtInter != 0) { for (int l = 1; l <= NbPtInter; l += 2) { // check if the point Inter is in the interval // PntInter[i] PntInter[i+1] // in which case it is not necessary to store it (problem with duplicates) if ((Inter > PntInter->Value(i).Value(l)) && (Inter < PntInter->Value(i).Value(l + 1))) { l = NbPtInter + 2; // leave the loop without storing the point } else { if (l + 1 >= NbPtInter) { // one has parsed the entire table : the point // does not belong to a common point interval if (NbPtG1G1 != 0) { // if there exists an incompatible interval for (int k = 1; k <= NbPtG1G1; k += 2) { if ((Inter > PntG1G1->Value(i).Value(k)) && (Inter < PntG1G1->Value(i).Value(k + 1))) { k = NbPtG1G1 + 2; // to leave the loop // Add points of constraint 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 // the point does not belong to interval G1 { if (k + 1 >= NbPtG1G1) { myPlateCont->ChangeValue(i).Append(Inter); // add the point } } } } else { myPlateCont->ChangeValue(i).Append(Inter); // add the point } } } } } else { if (NbPtG1G1 != 0) // there exist an incompatible interval { for (int k = 1; k <= NbPtG1G1; k += 2) { if ((Inter > PntG1G1->Value(i).Value(k)) && (Inter < PntG1G1->Value(i).Value(k + 1))) { k = NbPtG1G1 + 2; // to leave the loop // Add points of constraint 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 // the point does not belong to interval G1 { if (k + 1 >= NbPtG1G1) { myPlateCont->ChangeValue(i).Append(Inter); // add the point } } } } else { if (((!mySurfInitIsGive) && (Geom2dAdaptor_Curve(LinCont->Curve2dOnSurf()).GetType() != GeomAbs_Circle)) || ((j > 1) && (j < NbPnt_i))) { // exclude extremities myPlateCont->ChangeValue(i).Append(Inter); // add the point } } } } } } //--------------------------------------------------------- // Function : CalculNbPtsInit //--------------------------------------------------------- // Calculate the number of points by curve depending on the // length for the first iteration //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::CalculNbPtsInit() { double LenT = 0; int NTLinCont = myLinCont->Length(); int NTPoint = (int)(myNbPtsOnCur * NTLinCont); int i; for (i = 1; i <= NTLinCont; i++) { LenT += myLinCont->Value(i)->Length(); } for (i = 1; i <= NTLinCont; i++) { int Cont = myLinCont->Value(i)->Order(); switch (Cont) { case 0: // Case G0 *1.2 myLinCont->ChangeValue(i)->SetNbPoints( int(1.2 * NTPoint * (myLinCont->Value(i)->Length()) / LenT)); break; case 1: // Case G1 *1 myLinCont->ChangeValue(i)->SetNbPoints( int(NTPoint * (myLinCont->Value(i)->Length()) / LenT)); break; case 2: // Case G2 *0.7 myLinCont->ChangeValue(i)->SetNbPoints( int(0.7 * NTPoint * (myLinCont->Value(i)->Length()) / LenT)); break; } if (myLinCont->Value(i)->NbPoints() < 3) { myLinCont->ChangeValue(i)->SetNbPoints(3); } } } //--------------------------------------------------------- // Function : LoadCurve //--------------------------------------------------------- // Starting from table myParCont load all the points noted in plate //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::LoadCurve(const int NbBoucle, const int OrderMax) { gp_Pnt P3d, Pdif, PP; gp_Pnt2d P2d; int NTLinCont = myLinCont->Length(), i, j; int Tang, Nt; for (i = 1; i <= NTLinCont; i++) { occ::handle CC = myLinCont->Value(i); if (CC->Order() != -1) { Tang = std::min(CC->Order(), OrderMax); Nt = myPlateCont->Value(i).Length(); if (Tang != -1) { for (j = 1; j <= Nt; j++) { // Loading of points G0 on boundaries 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); // Loading of 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(); double 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); } } // Loading of 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 ); // } } } } } } } //--------------------------------------------------------- // Function : LoadPoint //--------------------------------------------------------- // void GeomPlate_BuildPlateSurface::LoadPoint(const int NbBoucle, void GeomPlate_BuildPlateSurface::LoadPoint(const int, const int OrderMax) { gp_Pnt P3d, Pdif, PP; gp_Pnt2d P2d; int NTPntCont = myPntCont->Length(); int Tang, i; // gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10; // Loading of points of point constraints 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 = std::min(myPntCont->Value(i)->Order(), OrderMax); if (Tang == 1) { // ==1 gp_Vec V1, V2, V3, V4; 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); } } // Loading of 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 ); // } } } } //--------------------------------------------------------- // Function : VerifSurface //--------------------------------------------------------- bool GeomPlate_BuildPlateSurface::VerifSurface(const int NbBoucle) { //====================================================================== // Calculate errors //====================================================================== int NTLinCont = myLinCont->Length(); bool Result = true; // variable for error calculation myG0Error = 0, myG1Error = 0, myG2Error = 0; for (int i = 1; i <= NTLinCont; i++) { occ::handle LinCont; LinCont = myLinCont->Value(i); if (LinCont->Order() != -1) { int NbPts_i = myParCont->Value(i).Length(); if (NbPts_i < 3) { NbPts_i = 4; } occ::handle> tdist = new NCollection_HArray1(1, NbPts_i - 1); occ::handle> tang = new NCollection_HArray1(1, NbPts_i - 1); occ::handle> tcourb = new NCollection_HArray1(1, NbPts_i - 1); EcartContraintesMil(i, tdist, tang, tcourb); double diffDistMax = 0, diffAngMax = 0; // double SdiffDist=0, SdiffAng=0; int NdiffDist = 0, NdiffAng = 0; for (int j = 1; j < NbPts_i; j++) { if (tdist->Value(j) > myG0Error) { myG0Error = tdist->Value(j); } if (tang->Value(j) > myG1Error) { myG1Error = tang->Value(j); } if (tcourb->Value(j) > myG2Error) { myG2Error = tcourb->Value(j); } double 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); } double diffDist = tdist->Value(j) - LinCont->G0Criterion(U), diffAng; if (LinCont->Order() > 0) { diffAng = tang->Value(j) - LinCont->G1Criterion(U); } else { diffAng = 0; } // find the maximum variation of error and calculate the average if (diffDist > 0) { diffDist = diffDist / LinCont->G0Criterion(U); if (diffDist > diffDistMax) { diffDistMax = diffDist; } // SdiffDist+=diffDist; NdiffDist++; } else if ((diffAng > 0) && (LinCont->Order() == 1)) { diffAng = diffAng / myLinCont->Value(i)->G1Criterion(U); if (diffAng > diffAngMax) { diffAngMax = diffAng; } // SdiffAng+=diffAng; NdiffAng++; } } if (NdiffDist > 0) { // at least one point is not acceptable in G0 double Coef; if (LinCont->Order() == 0) { Coef = 0.6 * std::log(diffDistMax + 7.4); // 7.4 corresponds to the calculation of min. coefficient = 1.2 is e^1.2/0.6 } else { Coef = std::log(diffDistMax + 3.3); } // 3.3 corresponds to calculation of min. coefficient = 1.2 donc e^1.2 if (Coef > 3) { Coef = 3; } // experimentally after the coefficient becomes bad for L cases if ((NbBoucle > 1) && (diffDistMax > 2)) { Coef = 1.6; } if (LinCont->NbPoints() >= std::floor(LinCont->NbPoints() * Coef)) { Coef = 2; // to provide increase of the number of points } LinCont->SetNbPoints(int(LinCont->NbPoints() * Coef)); Result = false; } else if (NdiffAng > 0) // at least 1 point is not acceptable in G1 { double Coef = 1.5; if ((LinCont->NbPoints() + 1) >= std::floor(LinCont->NbPoints() * Coef)) { Coef = 2; } LinCont->SetNbPoints(int(LinCont->NbPoints() * Coef)); Result = false; } } } if (!Result) { if (myFree && NbBoucle == 1) { myPrevPlate = myPlate; } myPlate.Init(); } return Result; } //--------------------------------------------------------- // Function : VerifPoint //--------------------------------------------------------- void GeomPlate_BuildPlateSurface::VerifPoints(double& Dist, double& Ang, double& Curv) const { int NTPntCont = myPntCont->Length(); gp_Pnt Pi, Pf; gp_Pnt2d P2d; gp_Vec v1i, v1f, v2i, v2f, v3i, v3f; Ang = 0; Dist = 0, Curv = 0; occ::handle PntCont; for (int 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 > (M_PI / 2)) { Ang = M_PI - Ang; } break; case 2: occ::handle Splate(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; } } } double GeomPlate_BuildPlateSurface::ComputeAnisotropie() const { if (myAnisotropie) { // Temporary return 1.0; } else { return 1.0; } } bool GeomPlate_BuildPlateSurface::IsOrderG1() const { bool result = true; for (int i = 1; i <= myLinCont->Length(); i++) { if (myLinCont->Value(i)->Order() < 1) { result = false; break; } } return result; }