1 // Created on: 1997-05-05
2 // Created by: Jerome LEMONIER
3 // Copyright (c) 1996-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
17 #include <Adaptor2d_HCurve2d.hxx>
18 #include <Adaptor3d_HCurve.hxx>
19 #include <Adaptor3d_CurveOnSurface.hxx>
20 #include <Approx_CurveOnSurface.hxx>
21 #include <Extrema_ExtPS.hxx>
22 #include <Extrema_POnSurf.hxx>
23 #include <GCPnts_AbscissaPoint.hxx>
24 #include <Geom2d_BezierCurve.hxx>
25 #include <Geom2d_BSplineCurve.hxx>
26 #include <Geom2d_Curve.hxx>
27 #include <Geom2dAdaptor_Curve.hxx>
28 #include <Geom2dInt_GInter.hxx>
29 #include <Geom_BSplineSurface.hxx>
30 #include <Geom_Plane.hxx>
31 #include <Geom_RectangularTrimmedSurface.hxx>
32 #include <Geom_Surface.hxx>
33 #include <GeomAdaptor.hxx>
34 #include <GeomAdaptor_HSurface.hxx>
35 #include <GeomAdaptor_Surface.hxx>
36 #include <GeomLProp_SLProps.hxx>
37 #include <GeomPlate_BuildAveragePlane.hxx>
38 #include <GeomPlate_BuildPlateSurface.hxx>
39 #include <GeomPlate_CurveConstraint.hxx>
40 #include <GeomPlate_HArray1OfSequenceOfReal.hxx>
41 #include <GeomPlate_MakeApprox.hxx>
42 #include <GeomPlate_PointConstraint.hxx>
43 #include <GeomPlate_SequenceOfAij.hxx>
44 #include <GeomPlate_Surface.hxx>
46 #include <gp_Pnt2d.hxx>
48 #include <gp_Vec2d.hxx>
49 #include <IntRes2d_IntersectionPoint.hxx>
50 #include <Law_Interpol.hxx>
51 #include <LocalAnalysis_SurfaceContinuity.hxx>
52 #include <Plate_D1.hxx>
53 #include <Plate_D2.hxx>
54 #include <Plate_FreeGtoCConstraint.hxx>
55 #include <Plate_GtoCConstraint.hxx>
56 #include <Plate_PinpointConstraint.hxx>
57 #include <Plate_Plate.hxx>
58 #include <Precision.hxx>
59 #include <ProjLib_CompProjectedCurve.hxx>
60 #include <ProjLib_HCompProjectedCurve.hxx>
61 #include <Standard_ConstructionError.hxx>
62 #include <Standard_RangeError.hxx>
63 #include <TColgp_Array1OfPnt2d.hxx>
64 #include <TColgp_HArray1OfPnt.hxx>
65 #include <TColgp_HArray2OfPnt.hxx>
66 #include <TColgp_SequenceOfVec.hxx>
67 #include <TColStd_HArray1OfReal.hxx>
68 #include <TColStd_SequenceOfInteger.hxx>
69 #include <Message_ProgressIndicator.hxx>
74 #include <DrawTrSurf.hxx>
75 #include <Draw_Marker3D.hxx>
76 #include <Draw_Marker2D.hxx>
79 // 1 : Display of Geometries and intermediary control
80 // 2 : Display of the number of constraints by curve + Intersection
81 // 3 : Dump of constraints in Plate
82 static Standard_Integer NbPlan = 0;
83 //static Standard_Integer NbCurv2d = 0;
84 static Standard_Integer NbMark = 0;
85 static Standard_Integer NbProj = 0;
89 #include <OSD_Chronometer.hxx>
90 #include <Geom_Surface.hxx>
91 static Standard_Integer Affich=0;
94 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
95 // =========================================================
96 // C O N S T R U C T O R S
97 // =========================================================
98 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
100 //---------------------------------------------------------
101 // Constructor compatible with the old version
102 //---------------------------------------------------------
103 GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface (
104 const Handle(TColStd_HArray1OfInteger)& NPoints,
105 const Handle(GeomPlate_HArray1OfHCurve)& TabCurve,
106 const Handle(TColStd_HArray1OfInteger)& Tang,
107 const Standard_Integer Degree,
108 const Standard_Integer NbIter,
109 const Standard_Real Tol2d,
110 const Standard_Real Tol3d,
111 const Standard_Real TolAng,
112 const Standard_Real ,
113 const Standard_Boolean Anisotropie
115 myAnisotropie(Anisotropie),
123 { Standard_Integer NTCurve=TabCurve->Length();// Number of linear constraints
124 myNbPtsOnCur = 0; // Different calculation of the number of points depending on the length
125 myLinCont = new GeomPlate_HSequenceOfCurveConstraint;
126 myPntCont = new GeomPlate_HSequenceOfPointConstraint;
128 throw Standard_ConstructionError("GeomPlate : Number of iteration must be >= 1");
131 throw Standard_ConstructionError("GeomPlate : the bounds Array is null");
132 if (Tang->Length()==0)
133 throw Standard_ConstructionError("GeomPlate : the constraints Array is null");
134 Standard_Integer nbp = 0;
136 for ( i=1;i<=NTCurve;i++)
137 { nbp+=NPoints->Value(i);
140 throw Standard_ConstructionError("GeomPlate : the resolution is impossible if the number of constraints points is 0");
142 throw Standard_ConstructionError("GeomPlate ; the degree resolution must be upper of 2");
143 // Filling fields passing from the old constructor to the new one
144 for(i=1;i<=NTCurve;i++)
145 { Handle(GeomPlate_CurveConstraint) Cont = new GeomPlate_CurveConstraint
146 ( TabCurve->Value(i),
149 myLinCont->Append(Cont);
151 mySurfInitIsGive=Standard_False;
153 myIsLinear = Standard_True;
154 myFree = Standard_False;
157 //------------------------------------------------------------------
158 // Constructor with initial surface and degree
159 //------------------------------------------------------------------
160 GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface (
161 const Handle(Geom_Surface)& Surf,
162 const Standard_Integer Degree,
163 const Standard_Integer NbPtsOnCur,
164 const Standard_Integer NbIter,
165 const Standard_Real Tol2d,
166 const Standard_Real Tol3d,
167 const Standard_Real TolAng,
168 const Standard_Real /*TolCurv*/,
169 const Standard_Boolean Anisotropie ) :
171 myAnisotropie(Anisotropie),
173 myNbPtsOnCur(NbPtsOnCur),
181 throw Standard_ConstructionError("GeomPlate : Number of iteration must be >= 1");
183 throw Standard_ConstructionError("GeomPlate : the degree must be above 2");
184 myLinCont = new GeomPlate_HSequenceOfCurveConstraint;
185 myPntCont = new GeomPlate_HSequenceOfPointConstraint;
186 mySurfInitIsGive=Standard_True;
188 myIsLinear = Standard_True;
189 myFree = Standard_False;
193 //---------------------------------------------------------
194 // Constructor with degree
195 //---------------------------------------------------------
196 GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface (
197 const Standard_Integer Degree,
198 const Standard_Integer NbPtsOnCur,
199 const Standard_Integer NbIter,
200 const Standard_Real Tol2d,
201 const Standard_Real Tol3d,
202 const Standard_Real TolAng,
203 const Standard_Real /*TolCurv*/,
204 const Standard_Boolean Anisotropie ) :
205 myAnisotropie(Anisotropie),
207 myNbPtsOnCur(NbPtsOnCur),
215 throw Standard_ConstructionError("GeomPlate : Number of iteration must be >= 1");
217 throw Standard_ConstructionError("GeomPlate : the degree resolution must be upper of 2");
218 myLinCont = new GeomPlate_HSequenceOfCurveConstraint;
219 myPntCont = new GeomPlate_HSequenceOfPointConstraint;
220 mySurfInitIsGive=Standard_False;
222 myIsLinear = Standard_True;
223 myFree = Standard_False;
226 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
227 // =========================================================
228 // P U B L I C M E T H O D S
229 // =========================================================
230 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
234 //-------------------------------------------------------------------------
235 // Function : TrierTab, internal Function, does not belong to class
236 //-------------------------------------------------------------------------
237 // Reorder the table of transformations
238 // After the call of CourbeJointive the order of curves is modified
239 // Ex : initial order of curves ==> A B C D E F
240 // In TabInit we note ==> 1 2 3 4 5 6
241 // after CourbeJointive ==> A E C B D F
242 // TabInit ==> 1 5 3 2 4 6
243 // after TrierTab the Table contains ==> 1 4 3 5 2 6
244 // It is also possible to access the 2nd curve by taking TabInit[2]
245 // i.e. the 4th from the table of classified curves
246 //-------------------------------------------------------------------------
247 static void TrierTab(Handle(TColStd_HArray1OfInteger)& Tab)
249 // Parse the table of transformations to find the initial order
250 Standard_Integer Nb=Tab->Length();
251 TColStd_Array1OfInteger TabTri(1,Nb);
252 for (Standard_Integer i=1;i<=Nb;i++)
253 TabTri.SetValue(Tab->Value(i),i);
254 Tab->ChangeArray1()=TabTri;
257 //---------------------------------------------------------
258 // Function : ProjectCurve
259 //---------------------------------------------------------
260 Handle(Geom2d_Curve) GeomPlate_BuildPlateSurface::ProjectCurve(const Handle(Adaptor3d_HCurve)& Curv)
262 // Project a curve on a plane
263 Handle(Geom2d_Curve) Curve2d ;
264 Handle(GeomAdaptor_HSurface) hsur = new GeomAdaptor_HSurface(mySurfInit);
267 ProjLib_CompProjectedCurve Projector(hsur, Curv, myTol3d/10, myTol3d/10);
269 Standard_Real UdebCheck, UfinCheck, ProjUdeb, ProjUfin;
270 UdebCheck = Curv->FirstParameter();
271 UfinCheck = Curv->LastParameter();
272 Projector.Bounds( 1, ProjUdeb, ProjUfin );
274 if (Projector.NbCurves() != 1 ||
275 Abs( UdebCheck -ProjUdeb ) > Precision::PConfusion() ||
276 Abs( UfinCheck -ProjUfin ) > Precision::PConfusion())
278 if (Projector.IsSinglePnt(1, P2d))
280 // solution in a point
281 TColgp_Array1OfPnt2d poles(1, 2);
283 Curve2d = new (Geom2d_BezierCurve) (poles);
287 Curve2d.Nullify(); // No continuous solution
289 std::cout << "BuildPlateSurace :: No continuous projection" << std::endl;
294 GeomAbs_Shape Continuity = GeomAbs_C1;
295 Standard_Integer MaxDegree = 10, MaxSeg;
296 Standard_Real Udeb, Ufin;
297 Handle(ProjLib_HCompProjectedCurve) HProjector =
298 new ProjLib_HCompProjectedCurve();
299 HProjector->Set(Projector);
300 Projector.Bounds(1, Udeb, Ufin);
302 MaxSeg = 20 + HProjector->NbIntervals(GeomAbs_C3);
303 Approx_CurveOnSurface appr(HProjector, hsur, Udeb, Ufin, myTol3d,
304 Continuity, MaxDegree, MaxSeg,
305 Standard_False, Standard_True);
307 Curve2d = appr.Curve2d();
312 sprintf(name,"proj_%d",++NbProj);
313 DrawTrSurf::Set(name, Curve2d);
318 //---------------------------------------------------------
319 // Function : ProjectedCurve
320 //---------------------------------------------------------
321 Handle(Adaptor2d_HCurve2d) GeomPlate_BuildPlateSurface::ProjectedCurve( Handle(Adaptor3d_HCurve)& Curv)
323 // Projection of a curve on the initial surface
325 Handle(GeomAdaptor_HSurface) hsur = new GeomAdaptor_HSurface(mySurfInit);
327 ProjLib_CompProjectedCurve Projector(hsur, Curv, myTolU/10, myTolV/10);
328 Handle(ProjLib_HCompProjectedCurve) HProjector =
329 new ProjLib_HCompProjectedCurve();
331 if (Projector.NbCurves() != 1) {
333 HProjector.Nullify(); // No continuous solution
335 std::cout << "BuildPlateSurace :: No continuous projection" << std::endl;
340 Standard_Real First1,Last1,First2,Last2;
341 First1=Curv->FirstParameter();
342 Last1 =Curv->LastParameter();
343 Projector.Bounds(1,First2,Last2);
345 if (Abs(First1-First2) <= Max(myTolU,myTolV) &&
346 Abs(Last1-Last2) <= Max(myTolU,myTolV))
349 HProjector->Set(Projector);
350 HProjector = Handle(ProjLib_HCompProjectedCurve)::
351 DownCast(HProjector->Trim(First2,Last2,Precision::PConfusion()));
355 HProjector.Nullify(); // No continuous solution
357 std::cout << "BuildPlateSurace :: No complete projection" << std::endl;
364 //---------------------------------------------------------
365 // Function : ProjectPoint
366 //---------------------------------------------------------
367 // Projects a point on the initial surface
368 //---------------------------------------------------------
369 gp_Pnt2d GeomPlate_BuildPlateSurface::ProjectPoint(const gp_Pnt &p3d)
372 Standard_Integer nearest = 1;
373 if( myProj.NbExt() > 1)
375 Standard_Real dist2mini = myProj.SquareDistance(1);
376 for (Standard_Integer i=2; i<=myProj.NbExt();i++)
378 if (myProj.SquareDistance(i) < dist2mini)
380 dist2mini = myProj.SquareDistance(i);
385 P=myProj.Point(nearest);
393 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
394 // =========================================================
395 // P U B L I C M E T H O D S
396 // =========================================================
397 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
398 //---------------------------------------------------------
400 //---------------------------------------------------------
401 // Initializes linear and point constraints
402 //---------------------------------------------------------
403 void GeomPlate_BuildPlateSurface::Init()
404 { myLinCont->Clear();
406 myPntCont = new GeomPlate_HSequenceOfPointConstraint;
407 myLinCont = new GeomPlate_HSequenceOfCurveConstraint;
410 //---------------------------------------------------------
411 // Function : LoadInitSurface
412 //---------------------------------------------------------
413 // Loads the initial surface
414 //---------------------------------------------------------
415 void GeomPlate_BuildPlateSurface::LoadInitSurface(const Handle(Geom_Surface)& Surf)
417 mySurfInitIsGive=Standard_True;
420 //---------------------------------------------------------
422 //---------------------------------------------------------
423 //---------------------------------------------------------
424 // Adds a linear constraint
425 //---------------------------------------------------------
426 void GeomPlate_BuildPlateSurface::
427 Add(const Handle(GeomPlate_CurveConstraint)& Cont)
428 { myLinCont->Append(Cont);
431 void GeomPlate_BuildPlateSurface::SetNbBounds( const Standard_Integer NbBounds )
433 myNbBounds = NbBounds;
436 //---------------------------------------------------------
438 //---------------------------------------------------------
439 //---------------------------------------------------------
440 // Adds a point constraint
441 //---------------------------------------------------------
442 void GeomPlate_BuildPlateSurface::
443 Add(const Handle(GeomPlate_PointConstraint)& Cont)
445 myPntCont->Append(Cont);
448 //---------------------------------------------------------
449 // Function : Perform
450 // Calculates the surface filled with loaded constraints
451 //---------------------------------------------------------
452 void GeomPlate_BuildPlateSurface::Perform(const Handle(Message_ProgressIndicator) & aProgress)
456 OSD_Chronometer Chrono;
462 myNbBounds = myLinCont->Length();
465 //=====================================================================
466 // Declaration of variables.
467 //=====================================================================
468 Standard_Integer NTLinCont = myLinCont->Length(),
469 NTPntCont = myPntCont->Length(), NbBoucle=0;
470 Standard_Boolean Fini=Standard_True;
471 if ((NTLinCont + NTPntCont) == 0)
474 std::cout << "WARNING : GeomPlate : The number of constraints is null." << std::endl;
480 //======================================================================
482 //======================================================================
483 if (!mySurfInitIsGive)
484 ComputeSurfInit(aProgress);
488 { // Table of transformations to preserve the initial order, see TrierTab
489 myInitOrder = new TColStd_HArray1OfInteger(1,NTLinCont);
490 for (Standard_Integer l=1;l<=NTLinCont;l++)
491 myInitOrder->SetValue(l,l);
492 if (!CourbeJointive(myTol3d))
493 {// throw Standard_Failure("Curves are not joined");
495 std::cout<<"WARNING : Courbes non jointives a "<<myTol3d<<" pres"<<std::endl;
498 TrierTab(myInitOrder); // Reorder the table of transformations
500 else if(NTLinCont > 0)//Patch
502 mySense = new TColStd_HArray1OfInteger( 1, NTLinCont, 0 );
503 myInitOrder = new TColStd_HArray1OfInteger( 1, NTLinCont, 1 );
507 if (mySurfInit.IsNull())
512 Standard_Real u1,v1,u2,v2;
513 mySurfInit->Bounds(u1,v1,u2,v2);
514 GeomAdaptor_Surface aSurfInit(mySurfInit);
515 myTolU = aSurfInit.UResolution(myTol3d);
516 myTolV = aSurfInit.VResolution(myTol3d);
517 myProj.Initialize(aSurfInit,u1,v1,u2,v2,
520 //======================================================================
521 // Projection of curves
522 //======================================================================
523 Standard_Boolean Ok = Standard_True;
524 for (Standard_Integer i = 1; i <= NTLinCont; i++)
525 if(myLinCont->Value(i)->Curve2dOnSurf().IsNull())
527 Handle( Geom2d_Curve ) Curve2d = ProjectCurve( myLinCont->Value(i)->Curve3d() );
528 if (Curve2d.IsNull())
533 myLinCont->ChangeValue(i)->SetCurve2dOnSurf( Curve2d );
537 GeomPlate_MakeApprox App(myGeomPlateSurface, myTol3d,
542 mySurfInit = App.Surface();
544 mySurfInit->Bounds(u1,v1,u2,v2);
545 GeomAdaptor_Surface Surf(mySurfInit);
546 myTolU = Surf.UResolution(myTol3d);
547 myTolV = Surf.VResolution(myTol3d);
548 myProj.Initialize(Surf,u1,v1,u2,v2,
552 for (Standard_Integer i = 1; i <= NTLinCont; i++)
554 Handle( Geom2d_Curve ) Curve2d = ProjectCurve( myLinCont->Value(i)->Curve3d() );
555 if (Curve2d.IsNull())
560 myLinCont->ChangeValue(i)->SetCurve2dOnSurf( Curve2d );
564 mySurfInit = myPlanarSurfInit;
566 mySurfInit->Bounds(u1,v1,u2,v2);
567 GeomAdaptor_Surface SurfNew(mySurfInit);
568 myTolU = SurfNew.UResolution(myTol3d);
569 myTolV = SurfNew.VResolution(myTol3d);
570 myProj.Initialize(SurfNew,u1,v1,u2,v2,
573 for (Standard_Integer i = 1; i <= NTLinCont; i++)
574 myLinCont->ChangeValue(i)->
575 SetCurve2dOnSurf(ProjectCurve( myLinCont->Value(i)->Curve3d() ) );
577 else { // Project the points
578 for (Standard_Integer i=1; i<=NTPntCont; i++) {
580 myPntCont->Value(i)->D0(P);
581 myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P));
586 //======================================================================
587 // Projection of points
588 //======================================================================
589 for (Standard_Integer i=1;i<=NTPntCont;i++) {
590 if (! myPntCont->Value(i)->HasPnt2dOnSurf()) {
592 myPntCont->Value(i)->D0(P);
593 myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P));
597 //======================================================================
598 // Number of points by curve
599 //======================================================================
600 if ((NTLinCont !=0)&&(myNbPtsOnCur!=0))
603 //======================================================================
604 // Management of incompatibilites between curves
605 //======================================================================
606 Handle(GeomPlate_HArray1OfSequenceOfReal) PntInter;
607 Handle(GeomPlate_HArray1OfSequenceOfReal) PntG1G1;
608 if (NTLinCont != 0) {
609 PntInter = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
610 PntG1G1 = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
611 Intersect(PntInter, PntG1G1);
614 //======================================================================
615 // Loop to obtain a better surface
616 //======================================================================
618 myFree = !myIsLinear;
623 if (Affich && NbBoucle) {
624 std::cout<<"Resultats boucle"<< NbBoucle << std::endl;
625 std::cout<<"DistMax="<<myG0Error<<std::endl;
627 std::cout<<"AngleMax="<<myG1Error<<std::endl;
629 std::cout<<"CourbMax="<<myG2Error<<std::endl;
634 { //====================================================================
635 // Calculate the total number of points and the maximum of points by curve
636 //====================================================================
637 Standard_Integer NPointMax=0;
638 for (Standard_Integer i=1;i<=NTLinCont;i++)
639 if ((myLinCont->Value(i)->NbPoints())>NPointMax)
640 NPointMax=(Standard_Integer )( myLinCont->Value(i)->NbPoints());
641 //====================================================================
642 // Discretization of curves
643 //====================================================================
644 Discretise(PntInter, PntG1G1);
645 //====================================================================
646 // Preparation of constraint points for plate
647 //====================================================================
648 LoadCurve( NbBoucle );
649 if ( myPntCont->Length() != 0)
650 LoadPoint( NbBoucle );
651 //====================================================================
652 // Construction of the surface
653 //====================================================================
655 myPlate.SolveTI(myDegree, ComputeAnisotropie(), aProgress);
657 if (!aProgress.IsNull() && aProgress->UserBreak())
662 if (!myPlate.IsDone())
665 std::cout << "WARNING : GeomPlate : calculation of Plate failed" << std::endl;
670 myGeomPlateSurface = new GeomPlate_Surface(mySurfInit,myPlate);
671 Standard_Real Umin,Umax,Vmin,Vmax;
672 myPlate.UVBox(Umin,Umax,Vmin,Vmax);
673 myGeomPlateSurface->SetBounds(Umin,Umax,Vmin,Vmax);
675 Fini = VerifSurface(NbBoucle);
676 if ((NbBoucle >= myNbIter)&&(!Fini))
679 std::cout<<"Warning: objective was not reached"<<std::endl;
681 Fini = Standard_True;
684 if ((NTPntCont != 0)&&(Fini))
685 { Standard_Real di,an,cu;
686 VerifPoints(di,an,cu);
690 { LoadPoint( NbBoucle );
691 //====================================================================
692 //Construction of the surface
693 //====================================================================
694 myPlate.SolveTI(myDegree, ComputeAnisotropie(), aProgress);
696 if (!aProgress.IsNull() && aProgress->UserBreak())
701 if (!myPlate.IsDone())
704 std::cout << "WARNING : GeomPlate : calculation of Plate failed" << std::endl;
709 myGeomPlateSurface = new GeomPlate_Surface(mySurfInit,myPlate);
710 Standard_Real Umin,Umax,Vmin,Vmax;
711 myPlate.UVBox(Umin,Umax,Vmin,Vmax);
712 myGeomPlateSurface->SetBounds(Umin,Umax,Vmin,Vmax);
713 Fini = Standard_True;
714 Standard_Real di,an,cu;
715 VerifPoints(di,an,cu);
717 } while (!Fini); // End loop for better surface
720 { std::cout<<"======== Global results ==========="<<std::endl;
721 std::cout<<"DistMax="<<myG0Error<<std::endl;
723 std::cout<<"AngleMax="<<myG1Error<<std::endl;
725 std::cout<<"CourbMax="<<myG2Error<<std::endl;
730 std::cout<<"*** END OF GEOMPLATE ***"<<std::endl;
731 std::cout<<"Time of calculation : "<<Tps<<std::endl;
732 std::cout<<"Number of loops : "<<NbBoucle<<std::endl;
736 //---------------------------------------------------------
737 // Function : EcartContraintesMIL
738 //---------------------------------------------------------
739 void GeomPlate_BuildPlateSurface::
740 EcartContraintesMil ( const Standard_Integer c,
741 Handle(TColStd_HArray1OfReal)& d,
742 Handle(TColStd_HArray1OfReal)& an,
743 Handle(TColStd_HArray1OfReal)& courb)
745 Standard_Integer NbPt=myParCont->Value(c).Length();
750 NbPt=myParCont->Value(c).Length();
751 gp_Vec v1i, v1f, v2i, v2f, v3i, v3f;
753 gp_Pnt2d P2d;Standard_Integer i;
754 Handle(GeomPlate_CurveConstraint) LinCont = myLinCont->Value(c);
755 switch (LinCont->Order())
757 for (i=1; i<NbPt; i++)
758 { U = ( myParCont->Value(c).Value(i) +
759 myParCont->Value(c).Value(i+1) )/2;
761 if (!LinCont->ProjectedCurve().IsNull())
762 P2d = LinCont->ProjectedCurve()->Value(U);
764 { if (!LinCont->Curve2dOnSurf().IsNull())
765 P2d = LinCont->Curve2dOnSurf()->Value(U);
767 P2d = ProjectPoint(Pi);
769 myGeomPlateSurface->D0( P2d.Coord(1), P2d.Coord(2), Pf);
772 d->ChangeValue(i) = Pf.Distance(Pi);
776 for (i=1; i<NbPt; i++)
777 { U = ( myParCont->Value(c).Value(i) +
778 myParCont->Value(c).Value(i+1) )/2;
779 LinCont->D1(U, Pi, v1i, v2i);
780 if (!LinCont->ProjectedCurve().IsNull())
781 P2d = LinCont->ProjectedCurve()->Value(U);
783 { if (!LinCont->Curve2dOnSurf().IsNull())
784 P2d = LinCont->Curve2dOnSurf()->Value(U);
786 P2d = ProjectPoint(Pi);
788 myGeomPlateSurface->D1( P2d.Coord(1), P2d.Coord(2), Pf, v1f, v2f);
789 d->ChangeValue(i) = Pf.Distance(Pi);
790 v3i = v1i^v2i; v3f=v1f^v2f;
791 Standard_Real angle=v3f.Angle(v3i);
793 an->ChangeValue(i) = M_PI -angle;
795 an->ChangeValue(i) = angle;
800 Handle(Geom_Surface) Splate (myGeomPlateSurface);
801 LocalAnalysis_SurfaceContinuity CG2;
802 for (i=1; i<NbPt; i++)
803 { U = ( myParCont->Value(c).Value(i) +
804 myParCont->Value(c).Value(i+1) )/2;
806 if (!LinCont->ProjectedCurve().IsNull())
807 P2d = LinCont->ProjectedCurve()->Value(U);
809 { if (!LinCont->Curve2dOnSurf().IsNull())
810 P2d = LinCont->Curve2dOnSurf()->Value(U);
812 P2d = ProjectPoint(Pi);
814 GeomLProp_SLProps Prop (Splate, P2d.Coord(1), P2d.Coord(2),
816 CG2.ComputeAnalysis(Prop, myLinCont->Value(c)->LPropSurf(U),
818 d->ChangeValue(i)=CG2.C0Value();
819 an->ChangeValue(i)=CG2.G1Angle();
820 courb->ChangeValue(i)=CG2.G2CurvatureGap();
828 //---------------------------------------------------------
829 // Function : Disc2dContour
830 //---------------------------------------------------------
831 void GeomPlate_BuildPlateSurface::
832 Disc2dContour ( const Standard_Integer /*nbp*/,
833 TColgp_SequenceOfXY& Seq2d)
836 if (Seq2d.Length()!=4)
837 std::cout<<"Number of points should be equal to 4 for Disc2dContour"<<std::endl;
842 // sampling in "cosine" + 3 points on each interval
843 Standard_Integer NTCurve = myLinCont->Length();
844 Standard_Integer NTPntCont = myPntCont->Length();
848 Standard_Real u1,v1,u2,v2;
851 mySurfInit->Bounds(u1,v1,u2,v2);
852 GeomAdaptor_Surface Surf(mySurfInit);
853 myProj.Initialize(Surf,u1,v1,u2,v2,
856 for( i=1; i<=NTPntCont; i++)
857 if (myPntCont->Value(i)->Order()!=-1)
858 { P2d = myPntCont->Value(i)->Pnt2dOnSurf();
859 UV.SetX(P2d.Coord(1));
860 UV.SetY(P2d.Coord(2));
863 for(i=1; i<=NTCurve; i++)
865 Handle(GeomPlate_CurveConstraint) LinCont = myLinCont->Value(i);
866 if (LinCont->Order()!=-1)
867 { Standard_Integer NbPt=myParCont->Value(i).Length();
868 // first point of constraint (j=0)
869 if (!LinCont->ProjectedCurve().IsNull())
870 P2d = LinCont->ProjectedCurve()->Value(myParCont->Value(i).Value(1));
873 { if (!LinCont->Curve2dOnSurf().IsNull())
874 P2d = LinCont->Curve2dOnSurf()->Value(myParCont->Value(i).Value(1));
878 LinCont->D0(myParCont->Value(i).Value(1),PP);
879 P2d = ProjectPoint(PP);
883 UV.SetX(P2d.Coord(1));
884 UV.SetY(P2d.Coord(2));
886 for (Standard_Integer j=2; j<NbPt; j++)
887 { Standard_Real Uj=myParCont->Value(i).Value(j),
888 Ujp1=myParCont->Value(i).Value(j+1);
889 if (!LinCont->ProjectedCurve().IsNull())
890 P2d = LinCont->ProjectedCurve()->Value((Ujp1+3*Uj)/4);
893 { if (!LinCont->Curve2dOnSurf().IsNull())
894 P2d = LinCont->Curve2dOnSurf()->Value((Ujp1+3*Uj)/4);
898 LinCont->D0((Ujp1+3*Uj)/4,PP);
899 P2d = ProjectPoint(PP);
902 UV.SetX(P2d.Coord(1));
903 UV.SetY(P2d.Coord(2));
905 // point 1/2 previous
906 if (!LinCont->ProjectedCurve().IsNull())
907 P2d = LinCont->ProjectedCurve()->Value((Ujp1+Uj)/2);
910 { if (!LinCont->Curve2dOnSurf().IsNull())
911 P2d = LinCont->Curve2dOnSurf()->Value((Ujp1+Uj)/2);
915 LinCont->D0((Ujp1+Uj)/2,PP);
916 P2d = ProjectPoint(PP);
920 UV.SetX(P2d.Coord(1));
921 UV.SetY(P2d.Coord(2));
923 // point 3/4 previous
924 if (!LinCont->ProjectedCurve().IsNull())
925 P2d = LinCont->ProjectedCurve()->Value((3*Ujp1+Uj)/4);
928 { if (!LinCont->Curve2dOnSurf().IsNull())
929 P2d = LinCont->Curve2dOnSurf()->Value((3*Ujp1+Uj)/4);
933 LinCont->D0((3*Ujp1+Uj)/4,PP);
934 P2d = ProjectPoint(PP);
938 UV.SetX(P2d.Coord(1));
939 UV.SetY(P2d.Coord(2));
941 // current constraint point
942 if (!LinCont->ProjectedCurve().IsNull())
943 P2d = LinCont->ProjectedCurve()->Value(Ujp1);
946 { if (!LinCont->Curve2dOnSurf().IsNull())
947 P2d = LinCont->Curve2dOnSurf()->Value(Ujp1);
951 LinCont->D0(Ujp1,PP);
952 P2d = ProjectPoint(PP);
956 UV.SetX(P2d.Coord(1));
957 UV.SetY(P2d.Coord(2));
964 //---------------------------------------------------------
965 // Function : Disc3dContour
966 //---------------------------------------------------------
967 void GeomPlate_BuildPlateSurface::
968 Disc3dContour (const Standard_Integer /*nbp*/,
969 const Standard_Integer iordre,
970 TColgp_SequenceOfXYZ& Seq3d)
973 if (Seq3d.Length()!=4)
974 std::cout<<"nbp should be equal to 4 for Disc3dContour"<<std::endl;
975 if (iordre!=0&&iordre!=1)
976 std::cout<<"incorrect order for Disc3dContour"<<std::endl;
980 // sampling in "cosine" + 3 points on each interval
981 Standard_Real u1,v1,u2,v2;
982 mySurfInit->Bounds(u1,v1,u2,v2);
983 GeomAdaptor_Surface Surf(mySurfInit);
984 myProj.Initialize(Surf,u1,v1,u2,v2,
985 Surf.UResolution(myTol3d),
986 Surf.VResolution(myTol3d));
987 Standard_Integer NTCurve = myLinCont->Length();
988 Standard_Integer NTPntCont = myPntCont->Length();
995 for( i=1; i<=NTPntCont; i++)
996 if (myPntCont->Value(i)->Order()!=-1)
998 { myPntCont->Value(i)->D0(P3d);
1005 myPntCont->Value(i)->D1(P3d,v1h,v2h);
1014 for(i=1; i<=NTCurve; i++)
1015 if (myLinCont->Value(i)->Order()!=-1)
1017 { Standard_Integer NbPt=myParCont->Value(i).Length();;
1018 // first constraint point (j=0)
1019 // Standard_Integer NbPt=myParCont->Length();
1021 myLinCont->Value(i)->D0(myParCont->Value(i).Value(1),P3d);
1028 myLinCont->Value(i)->D1(myParCont->Value(i).Value(1),P3d,v1h,v2h);
1036 for (Standard_Integer j=2; j<NbPt; j++)
1037 { Standard_Real Uj=myParCont->Value(i).Value(j),
1038 Ujp1=myParCont->Value(i).Value(j+1);
1040 // point 1/4 previous
1041 myLinCont->Value(i)->D0((Ujp1+3*Uj)/4,P3d);
1046 // point 1/2 previous
1047 myLinCont->Value(i)->D0((Ujp1+Uj)/2,P3d);
1052 // point 3/4 previous
1053 myLinCont->Value(i)->D0((3*Ujp1+Uj)/4,P3d);
1058 // current constraint point
1059 myLinCont->Value(i)->D0(Ujp1,P3d);
1066 // point 1/4 previous
1067 myLinCont->Value(i)->D1((Ujp1+3*Uj)/4,P3d,v1h,v2h);
1073 // point 1/2 previous
1074 myLinCont->Value(i)->D1((Ujp1+Uj)/2,P3d,v1h,v2h);
1080 // point 3/4 previous
1081 myLinCont->Value(i)->D1((3*Ujp1+Uj)/4,P3d,v1h,v2h);
1087 // current constraint point
1088 myLinCont->Value(i)->D1(Ujp1,P3d,v1h,v2h);
1101 //---------------------------------------------------------
1102 // Function : IsDone
1103 //---------------------------------------------------------
1104 Standard_Boolean GeomPlate_BuildPlateSurface::IsDone() const
1105 { return myPlate.IsDone();
1110 //---------------------------------------------------------
1111 // Function : Surface
1112 //---------------------------------------------------------
1113 Handle(GeomPlate_Surface) GeomPlate_BuildPlateSurface::Surface() const
1114 { return myGeomPlateSurface ;
1117 //---------------------------------------------------------
1118 // Function : SurfInit
1119 //---------------------------------------------------------
1120 Handle(Geom_Surface) GeomPlate_BuildPlateSurface::SurfInit() const
1121 { return mySurfInit ;
1125 //---------------------------------------------------------
1127 //---------------------------------------------------------
1128 Handle(TColStd_HArray1OfInteger) GeomPlate_BuildPlateSurface::Sense() const
1129 { Standard_Integer NTCurve = myLinCont->Length();
1130 Handle(TColStd_HArray1OfInteger) Sens = new TColStd_HArray1OfInteger(1,
1132 for (Standard_Integer i=1; i<=NTCurve; i++)
1133 Sens->SetValue(i,mySense->Value(myInitOrder->Value(i)));
1139 //---------------------------------------------------------
1140 // Function : Curve2d
1141 //---------------------------------------------------------
1142 Handle(TColGeom2d_HArray1OfCurve) GeomPlate_BuildPlateSurface::Curves2d() const
1143 { Standard_Integer NTCurve = myLinCont->Length();
1144 Handle(TColGeom2d_HArray1OfCurve) C2dfin =
1145 new TColGeom2d_HArray1OfCurve(1,NTCurve);
1147 for (Standard_Integer i=1; i<=NTCurve; i++)
1148 C2dfin->SetValue(i,myLinCont->Value(myInitOrder->Value(i))->Curve2dOnSurf());
1154 //---------------------------------------------------------
1156 //---------------------------------------------------------
1157 Handle(TColStd_HArray1OfInteger) GeomPlate_BuildPlateSurface::Order() const
1158 { Handle(TColStd_HArray1OfInteger) result=
1159 new TColStd_HArray1OfInteger(1,myLinCont->Length());
1160 for (Standard_Integer i=1;i<=myLinCont->Length();i++)
1161 result->SetValue(myInitOrder->Value(i),i);
1166 //---------------------------------------------------------
1167 // Function : G0Error
1168 //---------------------------------------------------------
1169 Standard_Real GeomPlate_BuildPlateSurface::G0Error() const
1173 //---------------------------------------------------------
1174 // Function : G1Error
1175 //---------------------------------------------------------
1176 Standard_Real GeomPlate_BuildPlateSurface::G1Error() const
1180 //---------------------------------------------------------
1181 // Function : G2Error
1182 //---------------------------------------------------------
1183 Standard_Real GeomPlate_BuildPlateSurface::G2Error() const
1187 //=======================================================================
1188 //function : G0Error
1190 //=======================================================================
1192 Standard_Real GeomPlate_BuildPlateSurface::G0Error( const Standard_Integer Index )
1194 Handle( TColStd_HArray1OfReal ) tdistance = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1195 Handle( TColStd_HArray1OfReal ) tangle = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1196 Handle( TColStd_HArray1OfReal ) tcurvature = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1197 EcartContraintesMil( Index, tdistance, tangle, tcurvature );
1198 Standard_Real MaxDistance = 0.;
1199 for (Standard_Integer i = 1; i <= myNbPtsOnCur; i++)
1200 if (tdistance->Value(i) > MaxDistance)
1201 MaxDistance = tdistance->Value(i);
1205 //=======================================================================
1206 //function : G1Error
1208 //=======================================================================
1210 Standard_Real GeomPlate_BuildPlateSurface::G1Error( const Standard_Integer Index )
1212 Handle( TColStd_HArray1OfReal ) tdistance = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1213 Handle( TColStd_HArray1OfReal ) tangle = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1214 Handle( TColStd_HArray1OfReal ) tcurvature = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1215 EcartContraintesMil( Index, tdistance, tangle, tcurvature );
1216 Standard_Real MaxAngle = 0.;
1217 for (Standard_Integer i = 1; i <= myNbPtsOnCur; i++)
1218 if (tangle->Value(i) > MaxAngle)
1219 MaxAngle = tangle->Value(i);
1223 //=======================================================================
1224 //function : G2Error
1226 //=======================================================================
1228 Standard_Real GeomPlate_BuildPlateSurface::G2Error( const Standard_Integer Index )
1230 Handle( TColStd_HArray1OfReal ) tdistance = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1231 Handle( TColStd_HArray1OfReal ) tangle = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1232 Handle( TColStd_HArray1OfReal ) tcurvature = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1233 EcartContraintesMil( Index, tdistance, tangle, tcurvature );
1234 Standard_Real MaxCurvature = 0.;
1235 for (Standard_Integer i = 1; i <= myNbPtsOnCur; i++)
1236 if (tcurvature->Value(i) > MaxCurvature)
1237 MaxCurvature = tcurvature->Value(i);
1238 return MaxCurvature;
1241 Handle(GeomPlate_CurveConstraint) GeomPlate_BuildPlateSurface::CurveConstraint( const Standard_Integer order) const
1243 return myLinCont->Value(order);
1245 Handle(GeomPlate_PointConstraint) GeomPlate_BuildPlateSurface::PointConstraint( const Standard_Integer order) const
1247 return myPntCont->Value(order);
1249 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
1250 // =========================================================
1251 // P R I V A T E M E T H O D S
1252 // =========================================================
1253 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
1255 //=======================================================================
1256 // Function : CourbeJointive
1257 // Purpose : Create a chain of curves to calculate the
1258 // initial surface with the method of max flow.
1259 // Return true if it is a closed contour.
1260 //=======================================================================
1262 Standard_Boolean GeomPlate_BuildPlateSurface::
1263 CourbeJointive(const Standard_Real tolerance)
1264 { Standard_Integer nbf = myLinCont->Length();
1265 Standard_Real Ufinal1,Uinit1,Ufinal2,Uinit2;
1266 mySense = new TColStd_HArray1OfInteger(1,nbf,0);
1267 Standard_Boolean result=Standard_True;
1268 Standard_Integer j=1,i;
1271 while (j <= (myNbBounds-1))
1273 Standard_Integer a=0;
1276 { result = Standard_False;
1280 { if (i > myNbBounds)
1281 { result = Standard_False;
1286 Uinit1=myLinCont->Value(j)->FirstParameter();
1287 Ufinal1=myLinCont->Value(j)->LastParameter();
1288 Uinit2=myLinCont->Value(i)->FirstParameter();
1289 Ufinal2=myLinCont->Value(i)->LastParameter();
1290 if (mySense->Value(j)==1)
1292 myLinCont->Value(j)->D0(Ufinal1,P1);
1293 myLinCont->Value(i)->D0(Uinit2,P2);
1294 if (P1.Distance(P2)<tolerance)
1296 Handle(GeomPlate_CurveConstraint) tampon = myLinCont->Value(j+1);
1297 myLinCont->SetValue(j+1,myLinCont->Value(i));
1298 myLinCont->SetValue(i,tampon);
1299 Standard_Integer Tmp=myInitOrder->Value(j+1);
1300 // See function TrierTab for the functioning of myInitOrder
1301 myInitOrder->SetValue(j+1,myInitOrder->Value(i));
1302 myInitOrder->SetValue(i,Tmp);
1307 mySense->SetValue(j+1,0);
1311 { myLinCont->Value(i)->D0(Ufinal2,P2);
1313 if (P1.Distance(P2)<tolerance)
1315 Handle(GeomPlate_CurveConstraint) tampon =
1316 myLinCont->Value(j+1);
1317 myLinCont->SetValue(j+1,myLinCont->Value(i));
1318 myLinCont->SetValue(i,tampon);
1319 Standard_Integer Tmp=myInitOrder->Value(j+1);
1320 // See function TrierTab for the functioning of myInitOrder
1321 myInitOrder->SetValue(j+1,myInitOrder->Value(i));
1322 myInitOrder->SetValue(i,Tmp);
1325 mySense->SetValue(j+1,1);
1333 Uinit1=myLinCont->Value( myNbBounds )->FirstParameter();
1334 Ufinal1=myLinCont->Value( myNbBounds )->LastParameter();
1335 Uinit2=myLinCont->Value(1)->FirstParameter();
1336 Ufinal2=myLinCont->Value(1)->LastParameter();
1337 myLinCont->Value( myNbBounds )->D0(Ufinal1,P1);
1338 myLinCont->Value(1)->D0(Uinit2,P2);
1339 if ((mySense->Value( myNbBounds )==0)
1340 &&(P1.Distance(P2)<tolerance))
1342 return ((Standard_True)&&(result));
1344 myLinCont->Value( myNbBounds )->D0(Uinit1,P1);
1345 if ((mySense->Value( myNbBounds )==1)
1346 &&(P1.Distance(P2)<tolerance))
1348 return ((Standard_True)&&(result));
1350 else return Standard_False;
1354 //-------------------------------------------------------------------------
1355 // Function : ComputeSurfInit
1356 //-------------------------------------------------------------------------
1357 // Calculate the initial surface either by the method of max flow or by
1358 // the method of the plane of inertia if the contour is not closed or if
1359 // there are point constraints.
1360 //-------------------------------------------------------------------------
1362 void GeomPlate_BuildPlateSurface::ComputeSurfInit(const Handle(Message_ProgressIndicator) & aProgress)
1364 Standard_Integer nopt=2, popt=2, Np=1;
1365 Standard_Boolean isHalfSpace = Standard_True;
1366 Standard_Real LinTol = 0.001, AngTol = 0.001; //AngTol = 0.0001; //LinTol = 0.0001
1368 // Option to calculate the initial plane
1369 Standard_Integer NTLinCont = myLinCont->Length(), NTPntCont = myPntCont->Length();
1370 // Table of transformation to preserve the initial order see TrierTab
1371 if (NTLinCont != 0) {
1372 myInitOrder = new TColStd_HArray1OfInteger(1,NTLinCont);
1373 for (Standard_Integer i = 1; i <= NTLinCont; i++)
1374 myInitOrder->SetValue( i, i );
1377 Standard_Boolean CourbeJoint = (NTLinCont != 0) && CourbeJointive (myTol3d);
1378 if (CourbeJoint && IsOrderG1())
1381 // Table contains the cloud of points for calculation of the plane
1382 Standard_Integer NbPoint = 20, Discr = NbPoint/4, pnum = 0;
1383 Handle( TColgp_HArray1OfPnt ) Pts = new TColgp_HArray1OfPnt( 1, (NbPoint+1)*NTLinCont+NTPntCont );
1384 TColgp_SequenceOfVec Vecs, NewVecs;
1385 GeomPlate_SequenceOfAij Aset;
1386 Standard_Real Uinit, Ufinal, Uif;
1388 Standard_Integer i ;
1389 for ( i = 1; i <= NTLinCont; i++)
1391 Standard_Integer Order = myLinCont->Value(i)->Order();
1395 Uinit = myLinCont->Value(i)->FirstParameter();
1396 Ufinal = myLinCont->Value(i)->LastParameter();
1397 Uif = Ufinal - Uinit;
1398 if (mySense->Value(i) == 1)
1404 gp_Vec Vec1, Vec2, Normal;
1405 Standard_Boolean ToReverse = Standard_False;
1406 if (i > 1 && Order >= GeomAbs_G1)
1409 myLinCont->Value(i)->D1( Uinit, P, Vec1, Vec2 );
1410 Normal = Vec1 ^ Vec2;
1411 if (LastVec.IsOpposite( Normal, AngTol ))
1412 ToReverse = Standard_True;
1415 for (Standard_Integer j = 0; j <= NbPoint; j++)
1416 { // Number of points per curve = 20
1417 // Linear distribution
1418 Standard_Real Inter = j*Uif/(NbPoint);
1419 if (Order < GeomAbs_G1 || j % Discr != 0)
1420 myLinCont->Value(i)->D0( Uinit+Inter, Pts->ChangeValue(++pnum) );
1423 myLinCont->Value(i)->D1( Uinit+Inter, Pts->ChangeValue(++pnum), Vec1, Vec2 );
1424 Normal = Vec1 ^ Vec2;
1428 Standard_Boolean isNew = Standard_True;
1429 Standard_Integer k ;
1430 for ( k = 1; k <= Vecs.Length(); k++)
1431 if (Vecs(k).IsEqual( Normal, LinTol, AngTol ))
1433 isNew = Standard_False;
1437 for (k = 1; k <= NewVecs.Length(); k++)
1438 if (NewVecs(k).IsEqual( Normal, LinTol, AngTol ))
1440 isNew = Standard_False;
1444 NewVecs.Append( Normal );
1447 if (Order >= GeomAbs_G1)
1449 isHalfSpace = GeomPlate_BuildAveragePlane::HalfSpace( NewVecs, Vecs, Aset, LinTol, AngTol );
1454 } //for (i = 1; i <= NTLinCont; i++)
1458 for (i = 1; i <= NTPntCont; i++)
1460 Standard_Integer Order = myPntCont->Value(i)->Order();
1463 gp_Vec Vec1, Vec2, Normal;
1464 if (Order < GeomAbs_G1)
1465 myPntCont->Value(i)->D0( Pts->ChangeValue(++pnum) );
1468 myPntCont->Value(i)->D1( Pts->ChangeValue(++pnum), Vec1, Vec2 );
1469 Normal = Vec1 ^ Vec2;
1471 Standard_Boolean isNew = Standard_True;
1472 for (Standard_Integer k = 1; k <= Vecs.Length(); k++)
1473 if (Vecs(k).IsEqual( Normal, LinTol, AngTol ))
1475 isNew = Standard_False;
1480 NewVecs.Append( Normal );
1481 isHalfSpace = GeomPlate_BuildAveragePlane::HalfSpace( NewVecs, Vecs, Aset, LinTol, AngTol );
1484 NewVecs(1).Reverse();
1485 isHalfSpace = GeomPlate_BuildAveragePlane::HalfSpace( NewVecs, Vecs, Aset, LinTol, AngTol );
1491 } //for (i = 1; i <= NTPntCont; i++)
1495 Standard_Boolean NullExist = Standard_True;
1498 NullExist = Standard_False;
1499 for (i = 1; i <= Vecs.Length(); i++)
1500 if (Vecs(i).SquareMagnitude() == 0.)
1502 NullExist = Standard_True;
1507 GeomPlate_BuildAveragePlane BAP( Vecs, Pts );
1508 Standard_Real u1,u2,v1,v2;
1509 BAP.MinMaxBox(u1,u2,v1,v2);
1510 // The space is greater for projections
1511 Standard_Real du = u2 - u1;
1512 Standard_Real dv = v2 - v1;
1513 u1 -= du; u2 += du; v1 -= dv; v2 += dv;
1514 mySurfInit = new Geom_RectangularTrimmedSurface( BAP.Plane(), u1,u2,v1,v2 );
1516 } //if (isHalfSpace)
1520 std::cout<<std::endl<<"Normals are not in half space"<<std::endl<<std::endl;
1522 myIsLinear = Standard_False;
1525 } //if (NTLinCont != 0 && (CourbeJoint = CourbeJointive( myTol3d )) && IsOrderG1())
1528 TrierTab( myInitOrder ); // Reorder the table of transformations
1532 if ( NTPntCont != 0)
1533 nopt = 1; //Calculate by the method of plane of inertia
1534 else if (!CourbeJoint || NTLinCont != myNbBounds)
1535 {// throw Standard_Failure("Curves are not joined");
1537 std::cout << "WARNING : Curves are non-adjacent with tolerance " << myTol3d << std::endl;
1542 Standard_Real LenT=0;
1543 Standard_Integer Npt=0;
1544 Standard_Integer NTPoint=20*NTLinCont;
1545 Standard_Integer i ;
1546 for ( i=1;i<=NTLinCont;i++)
1547 LenT+=myLinCont->Value(i)->Length();
1548 for (i=1;i<=NTLinCont;i++)
1549 { Standard_Integer NbPoint= (Standard_Integer )( NTPoint*(myLinCont->Value(i)->Length())/LenT);
1554 // Table containing a cloud of points for calculation of the plane
1555 Handle(TColgp_HArray1OfPnt) Pts = new TColgp_HArray1OfPnt(1,20*NTLinCont+NTPntCont);
1556 Standard_Integer NbPoint=20;
1557 Standard_Real Uinit,Ufinal, Uif;
1558 for ( i=1;i<=NTLinCont;i++)
1559 { Uinit=myLinCont->Value(i)->FirstParameter();
1560 Ufinal=myLinCont->Value(i)->LastParameter();
1562 if (mySense->Value(i) == 1)
1566 for (Standard_Integer j=0; j<NbPoint; j++)
1567 { // Number of points per curve = 20
1568 // Linear distribution
1569 Standard_Real Inter=j*Uif/(NbPoint);
1571 myLinCont->Value(i)->D0(Uinit+Inter,P);
1572 Pts->SetValue(Np++,P);
1575 for (i=1;i<=NTPntCont;i++)
1577 myPntCont->Value(i)->D0(P);
1578 Pts->SetValue(Np++,P);
1582 GeomPlate_BuildAveragePlane BAP( Pts, NbPoint*myNbBounds, myTol3d/1000, popt, nopt );
1586 std::cout << "WARNING : GeomPlate : the initial surface is not a plane." << std::endl;
1591 Standard_Real u1,u2,v1,v2;
1592 BAP.MinMaxBox(u1,u2,v1,v2);
1593 // The space is greater for projections
1594 Standard_Real du = u2 - u1;
1595 Standard_Real dv = v2 - v1;
1596 u1 -= du; u2 += du; v1 -= dv; v2 += dv;
1597 mySurfInit = new Geom_RectangularTrimmedSurface(BAP.Plane(),u1,u2,v1,v2);
1600 //Comparing metrics of curves and projected curves
1601 if (NTLinCont != 0 && myIsLinear)
1603 Handle( Geom_Surface ) InitPlane =
1604 (Handle( Geom_RectangularTrimmedSurface )::DownCast(mySurfInit))->BasisSurface();
1606 Standard_Real Ratio = 0., R1 = 2., R2 = 0.6; //R1 = 3, R2 = 0.5;//R1 = 1.4, R2 = 0.8; //R1 = 5., R2 = 0.2;
1607 Handle( GeomAdaptor_HSurface ) hsur =
1608 new GeomAdaptor_HSurface( InitPlane );
1609 Standard_Integer NbPoint = 20;
1611 // gp_Vec DerC, DerCproj, DU, DV;
1616 for (Standard_Integer i = 1; i <= NTLinCont && myIsLinear; i++)
1618 Standard_Real FirstPar = myLinCont->Value(i)->FirstParameter();
1619 Standard_Real LastPar = myLinCont->Value(i)->LastParameter();
1620 Standard_Real Uif = (LastPar - FirstPar)/(NbPoint);
1622 Handle( Adaptor3d_HCurve ) Curve = myLinCont->Value(i)->Curve3d();
1623 ProjLib_CompProjectedCurve Projector( hsur, Curve, myTol3d, myTol3d );
1624 Handle( ProjLib_HCompProjectedCurve ) ProjCurve =
1625 new ProjLib_HCompProjectedCurve();
1626 ProjCurve->Set( Projector );
1627 Adaptor3d_CurveOnSurface AProj(ProjCurve, hsur);
1630 gp_Vec DerC, DerCproj;
1631 for (Standard_Integer j = 1 ; j < NbPoint && myIsLinear; j++)
1633 Standard_Real Inter = FirstPar + j*Uif;
1634 Curve->D1(Inter, P, DerC);
1635 AProj.D1(Inter, P, DerCproj);
1637 Standard_Real A1 = DerC.Magnitude();
1638 Standard_Real A2 = DerCproj.Magnitude();
1645 if (Ratio > R1 || Ratio < R2)
1647 myIsLinear = Standard_False;
1654 std::cout <<"Metrics are too different :"<< Ratio<<std::endl;
1656 // myIsLinear = Standard_True; // !!
1657 } //comparing metrics of curves and projected curves
1662 myPlanarSurfInit = mySurfInit;
1666 sprintf(name,"planinit_%d",NbPlan+1);
1667 DrawTrSurf::Set(name, mySurfInit);
1670 Standard_Real u1,v1,u2,v2;
1671 mySurfInit->Bounds(u1,v1,u2,v2);
1672 GeomAdaptor_Surface Surf(mySurfInit);
1673 myTolU = Surf.UResolution(myTol3d);
1674 myTolV = Surf.VResolution(myTol3d);
1675 myProj.Initialize(Surf,u1,v1,u2,v2,
1678 //======================================================================
1679 // Projection of curves
1680 //======================================================================
1682 for (i = 1; i <= NTLinCont; i++)
1683 if (myLinCont->Value(i)->Curve2dOnSurf().IsNull())
1684 myLinCont->ChangeValue(i)->SetCurve2dOnSurf(ProjectCurve(myLinCont->Value(i)->Curve3d()));
1686 //======================================================================
1687 // Projection of points
1688 //======================================================================
1689 for (i = 1; i<=NTPntCont; i++)
1691 myPntCont->Value(i)->D0(P);
1692 if (!myPntCont->Value(i)->HasPnt2dOnSurf())
1693 myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P));
1696 //======================================================================
1697 // Number of points by curve
1698 //======================================================================
1699 if ((NTLinCont !=0)&&(myNbPtsOnCur!=0))
1702 //======================================================================
1703 // Management of incompatibilities between curves
1704 //======================================================================
1705 Handle(GeomPlate_HArray1OfSequenceOfReal) PntInter;
1706 Handle(GeomPlate_HArray1OfSequenceOfReal) PntG1G1;
1709 PntInter = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
1710 PntG1G1 = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
1711 Intersect(PntInter, PntG1G1);
1714 //====================================================================
1715 // Discretization of curves
1716 //====================================================================
1717 Discretise(PntInter, PntG1G1);
1718 //====================================================================
1719 //Preparation of points of constraint for plate
1720 //====================================================================
1722 if (myPntCont->Length() != 0)
1724 //====================================================================
1725 // Construction of the surface
1726 //====================================================================
1727 myPlate.SolveTI(2, ComputeAnisotropie(), aProgress);
1728 if (!aProgress.IsNull() && aProgress->UserBreak())
1733 if (!myPlate.IsDone())
1736 std::cout << "WARNING : GeomPlate : calculation of Plate failed" << std::endl;
1741 myGeomPlateSurface = new GeomPlate_Surface( mySurfInit, myPlate );
1743 GeomPlate_MakeApprox App(myGeomPlateSurface, myTol3d,
1747 mySurfInit = App.Surface();
1749 mySurfInitIsGive = Standard_True;
1750 myPlate.Init(); // Reset
1752 for (i=1;i<=NTLinCont;i++)
1754 Handle( Geom2d_Curve ) NullCurve;
1755 NullCurve.Nullify();
1756 myLinCont->ChangeValue(i)->SetCurve2dOnSurf( NullCurve);
1763 sprintf(name,"surfinit_%d",++NbPlan);
1764 DrawTrSurf::Set(name, mySurfInit);
1769 //---------------------------------------------------------
1770 // Function : Intersect
1771 //---------------------------------------------------------
1772 // Find intersections between 2d curves
1773 // If the intersection is compatible (in cases G1-G1)
1774 // remove the point on one of two curves
1775 //---------------------------------------------------------
1776 void GeomPlate_BuildPlateSurface::
1777 Intersect(Handle(GeomPlate_HArray1OfSequenceOfReal)& PntInter,
1778 Handle(GeomPlate_HArray1OfSequenceOfReal)& PntG1G1)
1780 Standard_Integer NTLinCont = myLinCont->Length();
1781 Geom2dInt_GInter Intersection;
1782 Geom2dAdaptor_Curve Ci, Cj;
1783 IntRes2d_IntersectionPoint int2d;
1787 // if (!mySurfInitIsGive)
1788 for (Standard_Integer i=1;i<=NTLinCont;i++)
1790 //Standard_Real NbPnt_i=myLinCont->Value(i)->NbPoints();
1791 // Find the intersection with each curve including the curve itself
1792 Ci.Load(myLinCont->Value(i)->Curve2dOnSurf());
1793 for(Standard_Integer j=i; j<=NTLinCont; j++)
1795 Cj.Load(myLinCont->Value(j)->Curve2dOnSurf());
1797 Intersection.Perform(Ci, myTol2d*10, myTol2d*10);
1799 Intersection.Perform(Ci, Cj, myTol2d*10, myTol2d*10);
1801 if (!Intersection.IsEmpty())
1802 { // there is one intersection
1803 Standard_Integer nbpt = Intersection.NbPoints();
1804 // number of points of intersection
1805 for (Standard_Integer k = 1; k <= nbpt; k++)
1806 { int2d = Intersection.Point(k);
1807 myLinCont->Value(i)->D0(int2d.ParamOnFirst(),P1);
1808 myLinCont->Value(j)->D0(int2d.ParamOnSecond(),P2);
1812 std::cout << " Intersection "<< k << " entre " << i
1813 << " &" << j << std::endl;
1814 std::cout <<" Distance = "<< P1.Distance(P2) << std::endl;
1817 if (P1.Distance( P2 ) < myTol3d)
1818 { // 2D intersection corresponds to close 3D points.
1819 // Note the interval, in which the point needs to be removed
1820 // to avoid duplications, which cause
1821 // error in plate. The point on curve i is removed;
1822 // the point on curve j is preserved;
1823 // the length of interval is a length 2d
1824 // corresponding in 3d to myTol3d
1825 Standard_Real tolint = Ci.Resolution(myTol3d);
1826 Ci.D1(int2d.ParamOnFirst(),P2d, V2d);
1827 Standard_Real aux = V2d.Magnitude();
1831 if (aux > 100*tolint) tolint*=100;
1836 PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() - tolint);
1837 PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() + tolint);
1839 if ( (myLinCont->Value(i)->Order()==1)
1840 &&(myLinCont->Value(j)->Order()==1))
1841 { gp_Vec v11,v12,v13,v14,v15,v16,v21,v22,v23,v24,v25,v26;
1842 myLinCont->Value(i)->D2( int2d.ParamOnFirst(), P1, v11, v12, v13, v14, v15 );
1843 myLinCont->Value(j)->D2( int2d.ParamOnSecond(), P2, v21, v22, v23, v24, v25 );
1844 v16=v11^v12;v26=v21^v22;
1845 Standard_Real ant=v16.Angle(v26);
1848 if ((Abs(v16*v15-v16*v25)>(myTol3d/1000))
1849 ||(Abs(ant)>myTol3d/1000))
1850 // Non-compatible ==> remove zone in constraint G1
1851 // corresponding to 3D tolerance of 0.01
1852 { Standard_Real coin;
1853 Standard_Real Tol= 100 * myTol3d;
1855 gp_Pnt2d P1temp,P2temp;
1857 myLinCont->Value(i)->Curve2dOnSurf()->D1( int2d.ParamOnFirst(), P1temp, V1);
1858 myLinCont->Value(j)->Curve2dOnSurf()->D1( int2d.ParamOnSecond(), P2temp, V2);
1862 if (Abs(Abs(A1)-M_PI)<myTolAng) Tol = 100000 * myTol3d;
1864 if (Affich) std::cout <<"Angle between curves "<<i<<","<<j
1865 <<" "<<Abs(Abs(A1)-M_PI)<<std::endl;
1868 coin = Ci.Resolution(Tol);
1869 Standard_Real Par1=int2d.ParamOnFirst()-coin,
1870 Par2=int2d.ParamOnFirst()+coin;
1871 // Storage of the interval for curve i
1872 PntG1G1->ChangeValue(i).Append(Par1);
1873 PntG1G1->ChangeValue(i).Append(Par2);
1874 coin = Cj.Resolution(Tol);
1875 Par1=int2d.ParamOnSecond()-coin;
1876 Par2=int2d.ParamOnSecond()+coin;
1877 // Storage of the interval for curve j
1878 PntG1G1->ChangeValue(j).Append(Par1);
1879 PntG1G1->ChangeValue(j).Append(Par2);
1883 if ((myLinCont->Value(i)->Order()==0 && myLinCont->Value(j)->Order()==1) ||
1884 (myLinCont->Value(i)->Order()==1 && myLinCont->Value(j)->Order()==0))
1886 gp_Vec vec, vecU, vecV, N;
1887 if (myLinCont->Value(i)->Order() == 0)
1889 Handle( Adaptor3d_HCurve ) theCurve = myLinCont->Value(i)->Curve3d();
1890 theCurve->D1( int2d.ParamOnFirst(), P1, vec );
1891 myLinCont->Value(j)->D1( int2d.ParamOnSecond(), P2, vecU, vecV );
1895 Handle( Adaptor3d_HCurve ) theCurve = myLinCont->Value(j)->Curve3d();
1896 theCurve->D1( int2d.ParamOnSecond(), P2, vec );
1897 myLinCont->Value(i)->D1( int2d.ParamOnFirst(), P1, vecU, vecV );
1900 Standard_Real Angle = vec.Angle( N );
1901 Angle = Abs( M_PI/2-Angle );
1902 if (Angle > myTolAng/10.) //????????? //if (Abs( scal ) > myTol3d/100)
1904 // Non-compatible ==> one removes zone in constraint G0 and G1
1905 // corresponding to 3D tolerance of 0.01
1907 Standard_Real Tol= 100 * myTol3d;
1909 gp_Pnt2d P1temp,P2temp;
1911 myLinCont->Value(i)->Curve2dOnSurf()->D1( int2d.ParamOnFirst(), P1temp, V1);
1912 myLinCont->Value(j)->Curve2dOnSurf()->D1( int2d.ParamOnSecond(), P2temp, V2);
1913 A1 = V1.Angle( V2 );
1916 if (Abs(Abs(A1) - M_PI) < myTolAng) Tol = 100000 * myTol3d;
1918 if (Affich) std::cout <<"Angle entre Courbe "<<i<<","<<j
1919 <<" "<<Abs(Abs(A1)-M_PI)<<std::endl;
1921 if (myLinCont->Value(i)->Order() == 1)
1923 coin = Ci.Resolution(Tol);
1924 coin *= Angle / myTolAng * 10.;
1926 std::cout<<std::endl<<"coin = "<<coin<<std::endl;
1928 Standard_Real Par1 = int2d.ParamOnFirst() - coin;
1929 Standard_Real Par2 = int2d.ParamOnFirst() + coin;
1930 // Storage of the interval for curve i
1931 PntG1G1->ChangeValue(i).Append( Par1 );
1932 PntG1G1->ChangeValue(i).Append( Par2 );
1936 coin = Cj.Resolution(Tol);
1937 coin *= Angle / myTolAng * 10.;
1939 std::cout<<std::endl<<"coin = "<<coin<<std::endl;
1941 Standard_Real Par1 = int2d.ParamOnSecond() - coin;
1942 Standard_Real Par2 = int2d.ParamOnSecond() + coin;
1943 // Storage of the interval for curve j
1944 PntG1G1->ChangeValue(j).Append( Par1 );
1945 PntG1G1->ChangeValue(j).Append( Par2 );
1949 } //if (P1.Distance( P2 ) < myTol3d)
1951 // 2D intersection corresponds to extended 3D points.
1952 // Note the interval where it is necessary to remove
1953 // the points to avoid duplications causing
1954 // error in plate. The point on curve i is removed,
1955 // the point on curve j is preserved.
1956 // The length of interval is 2D length
1957 // corresponding to the distance of points in 3D to myTol3d
1958 Standard_Real tolint, Dist;
1959 Dist = P1.Distance(P2);
1960 tolint = Ci.Resolution(Dist);
1961 PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() - tolint);
1962 PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() + tolint);
1965 tolint = Cj.Resolution(Dist);
1966 PntInter->ChangeValue(j).
1967 Append( int2d.ParamOnSecond() - tolint);
1968 PntInter->ChangeValue(j).
1969 Append( int2d.ParamOnSecond() + tolint);
1973 std::cout << "Attention: Two points 3d have the same projection dist = " << Dist << std::endl;
1978 Handle(Draw_Marker3D) mark = new (Draw_Marker3D)(P1, Draw_X, Draw_vert);
1980 sprintf(name,"mark_%d",++NbMark);
1981 Draw::Set(name, mark);
1991 //---------------------------------------------------------
1992 // Function : Discretize
1993 //---------------------------------------------------------
1994 // Discretize curves according to parameters
1995 // the table of sequences Parcont contains all
1996 // parameter of points on curves
1997 // Field myPlateCont contains parameter of points on a plate;
1998 // it excludes duplicate points and imcompatible zones.
1999 // The first part corresponds to verification of compatibility
2000 // and to removal of duplicate points.
2001 //---------------------------------------------------------
2002 void GeomPlate_BuildPlateSurface::
2003 Discretise(const Handle(GeomPlate_HArray1OfSequenceOfReal)& PntInter,
2004 const Handle(GeomPlate_HArray1OfSequenceOfReal)& PntG1G1)
2006 Standard_Integer NTLinCont = myLinCont->Length();
2007 Standard_Boolean ACR;
2008 Handle(Geom2d_Curve) C2d;
2009 Geom2dAdaptor_Curve AC2d;
2010 // Handle(Adaptor_HCurve2d) HC2d;
2011 Handle(Law_Interpol) acrlaw = new (Law_Interpol) ();
2012 myPlateCont = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
2013 myParCont = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
2016 //===========================================================================
2017 // Construction of the table containing parameters of constraint points
2018 //===========================================================================
2019 Standard_Real Uinit, Ufinal, Length2d=0, Inter;
2020 Standard_Real CurLength;
2021 Standard_Integer NbPnt_i, NbPtInter, NbPtG1G1;
2022 Handle(GeomPlate_CurveConstraint) LinCont;
2024 for (Standard_Integer i=1;i<=NTLinCont;i++) {
2025 LinCont = myLinCont->Value(i);
2026 Uinit=LinCont->FirstParameter();
2027 Ufinal=LinCont->LastParameter();
2028 // HC2d=LinCont->ProjectedCurve();
2029 // if(HC2d.IsNull())
2030 // ACR = (!HC2d.IsNull() || !C2d.IsNull());
2031 C2d= LinCont->Curve2dOnSurf();
2032 ACR = (!C2d.IsNull());
2034 // Construct a law close to curvilinear abscissa
2035 if(!C2d.IsNull()) AC2d.Load(C2d);
2036 // AC2d.Load(LinCont->Curve2dOnSurf());
2037 Standard_Integer ii, Nbint = 20;
2039 TColgp_Array1OfPnt2d tabP2d(1, Nbint+1);
2040 tabP2d(1).SetY(Uinit);
2042 tabP2d(Nbint+1).SetY(Ufinal);
2043 /* if (!HC2d.IsNull())
2045 Length2d = GCPnts_AbscissaPoint::Length(HC2d->Curve2d(), Uinit, Ufinal);
2047 Length2d = GCPnts_AbscissaPoint::Length(AC2d, Uinit, Ufinal);
2049 tabP2d(Nbint+1).SetX(Length2d);
2050 for (ii = 2; ii<= Nbint; ii++) {
2051 U = Uinit + (Ufinal-Uinit)*((1-Cos((ii-1)*M_PI/(Nbint)))/2);
2053 /* if (!HC2d.IsNull()) {
2054 Standard_Real L = GCPnts_AbscissaPoint::Length(HC2d->Curve2d(), Uinit, U);
2058 tabP2d(ii).SetX(GCPnts_AbscissaPoint::Length(AC2d, Uinit, U));
2060 acrlaw->Set(tabP2d);
2063 NbPnt_i = (Standard_Integer )( LinCont->NbPoints());
2064 NbPtInter= PntInter->Value(i).Length();
2065 NbPtG1G1= PntG1G1->Value(i).Length();
2069 std::cout << "Courbe : " << i << std::endl;
2070 std::cout << " NbPnt, NbPtInter, NbPtG1G1 :" << NbPnt_i << ", "
2071 << NbPtInter << ", " << NbPtG1G1 << std::endl;
2074 for (Standard_Integer j=1; j<=NbPnt_i; j++)
2076 // Distribution of points in cosine following ACR 2D
2077 // To avoid points of accumulation in 2D
2078 //Inter=Uinit+(Uif)*((-cos(M_PI*((j-1)/(NbPnt_i-1)))+1)/2);
2080 Inter=Ufinal;// to avoid bug on Sun
2082 CurLength = Length2d*(1-Cos((j-1)*M_PI/(NbPnt_i-1)))/2;
2083 Inter = acrlaw->Value(CurLength);
2086 Inter=Uinit+(Ufinal-Uinit)*((1-Cos((j-1)*M_PI/(NbPnt_i-1)))/2);
2088 myParCont->ChangeValue(i).Append(Inter);// add a point
2091 for(Standard_Integer l=1;l<=NbPtInter;l+=2)
2093 // check if the point Inter is in the interval
2094 // PntInter[i] PntInter[i+1]
2095 // in which case it is not necessary to store it (problem with duplicates)
2096 if ((Inter>PntInter->Value(i).Value(l))
2097 &&(Inter<PntInter->Value(i).Value(l+1)))
2100 // leave the loop without storing the point
2104 if (l+1>=NbPtInter) {
2105 // one has parsed the entire table : the point
2106 // does not belong to a common point interval
2109 // if there exists an incompatible interval
2110 for(Standard_Integer k=1;k<=NbPtG1G1;k+=2)
2112 if ((Inter>PntG1G1->Value(i).Value(k))
2113 &&(Inter<PntG1G1->Value(i).Value(k+1)))
2115 k=NbPtG1G1+2; // to leave the loop
2116 // Add points of constraint G0
2120 AC2d.D0(Inter, P2d);
2121 LinCont->D0(Inter,P3d);
2122 mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP);
2123 Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1),
2124 -PP.Coord(2)+P3d.Coord(2),
2125 -PP.Coord(3)+P3d.Coord(3));
2126 Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0);
2129 else // the point does not belong to interval G1
2133 myPlateCont->ChangeValue(i).Append(Inter);
2141 myPlateCont->ChangeValue(i).Append(Inter);
2150 if (NbPtG1G1!=0) // there exist an incompatible interval
2152 for(Standard_Integer k=1;k<=NbPtG1G1;k+=2)
2154 if ((Inter>PntG1G1->Value(i).Value(k))
2155 &&(Inter<PntG1G1->Value(i).Value(k+1)))
2157 k=NbPtG1G1+2; // to leave the loop
2158 // Add points of constraint G0
2162 AC2d.D0(Inter, P2d);
2163 LinCont->D0(Inter,P3d);
2164 mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP);
2165 Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1),
2166 -PP.Coord(2)+P3d.Coord(2),
2167 -PP.Coord(3)+P3d.Coord(3));
2168 Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0);
2172 else // the point does not belong to interval G1
2176 myPlateCont->ChangeValue(i).Append(Inter);
2184 if ( ( (!mySurfInitIsGive)
2185 &&(Geom2dAdaptor_Curve(LinCont->Curve2dOnSurf()).GetType()!=GeomAbs_Circle))
2186 || ( (j>1) &&(j<NbPnt_i))) // exclude extremeties
2187 myPlateCont->ChangeValue(i).Append(Inter);// add the point
2193 //---------------------------------------------------------
2194 // Function : CalculNbPtsInit
2195 //---------------------------------------------------------
2196 // Calculate the number of points by curve depending on the
2197 // length for the first iteration
2198 //---------------------------------------------------------
2199 void GeomPlate_BuildPlateSurface::CalculNbPtsInit ()
2201 Standard_Real LenT=0;
2202 Standard_Integer NTLinCont=myLinCont->Length();
2203 Standard_Integer NTPoint=(Standard_Integer )( myNbPtsOnCur*NTLinCont);
2204 Standard_Integer i ;
2206 for ( i=1;i<=NTLinCont;i++)
2207 LenT+=myLinCont->Value(i)->Length();
2208 for ( i=1;i<=NTLinCont;i++)
2209 { Standard_Integer Cont=myLinCont->Value(i)->Order();
2211 { case 0 : // Case G0 *1.2
2212 myLinCont->ChangeValue(i)->SetNbPoints(
2213 Standard_Integer(1.2*NTPoint*(myLinCont->Value(i)->Length())/LenT));
2215 case 1 : // Case G1 *1
2216 myLinCont->ChangeValue(i)->SetNbPoints(
2217 Standard_Integer(NTPoint*(myLinCont->Value(i)->Length())/LenT));
2219 case 2 : // Case G2 *0.7
2220 myLinCont->ChangeValue(i)->SetNbPoints(
2221 Standard_Integer(0.7*NTPoint*(myLinCont->Value(i)->Length())/LenT));
2224 if (myLinCont->Value(i)->NbPoints()<3)
2225 myLinCont->ChangeValue(i)->SetNbPoints(3);
2228 //---------------------------------------------------------
2229 // Function : LoadCurve
2230 //---------------------------------------------------------
2231 // Starting from table myParCont load all the points noted in plate
2232 //---------------------------------------------------------
2233 void GeomPlate_BuildPlateSurface::LoadCurve(const Standard_Integer NbBoucle,
2234 const Standard_Integer OrderMax )
2238 Standard_Integer NTLinCont=myLinCont->Length(), i, j;
2239 Standard_Integer Tang, Nt;
2242 for (i=1; i<=NTLinCont; i++){
2243 Handle(GeomPlate_CurveConstraint) CC = myLinCont->Value(i);
2244 if (CC ->Order()!=-1) {
2245 Tang = Min(CC->Order(), OrderMax);
2246 Nt = myPlateCont->Value(i).Length();
2248 for (j=1; j<=Nt; j++) {// Loading of points G0 on boundaries
2249 CC ->D0(myPlateCont->Value(i).Value(j),P3d);
2250 if (!CC->ProjectedCurve().IsNull())
2251 P2d = CC->ProjectedCurve()->Value(myPlateCont->Value(i).Value(j));
2254 if (!CC->Curve2dOnSurf().IsNull())
2255 P2d = CC->Curve2dOnSurf()->Value(myPlateCont->Value(i).Value(j));
2257 P2d = ProjectPoint(P3d);
2259 mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP);
2260 Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1),
2261 -PP.Coord(2)+P3d.Coord(2),
2262 -PP.Coord(3)+P3d.Coord(3));
2263 Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0);
2266 // Loading of points G1
2267 if (Tang==1) { // ==1
2269 CC->D1( myPlateCont->Value(i).Value(j), PP, V1, V2 );
2270 mySurfInit->D1( P2d.Coord(1), P2d.Coord(2), PP, V3, V4 );
2272 Plate_D1 D1final(V1.XYZ(),V2.XYZ());
2273 Plate_D1 D1init(V3.XYZ(),V4.XYZ());
2276 Plate_GtoCConstraint GCC( P2d.XY(), D1init, D1final );
2279 else if (NbBoucle == 1)
2281 Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), D1init, D1final );
2282 myPlate.Load( FreeGCC );
2285 gp_Vec DU, DV, Normal, DerPlateU, DerPlateV;
2288 //Normal.Normalize();
2289 Standard_Real norm = Normal.Magnitude();
2290 if (norm > 1.e-12) Normal /= norm;
2291 DerPlateU = myPrevPlate.EvaluateDerivative( P2d.XY(), 1, 0 );
2292 DerPlateV = myPrevPlate.EvaluateDerivative( P2d.XY(), 0, 1 );
2294 DU.SetLinearForm(-(V3 + DerPlateU).Dot(Normal), Normal, DerPlateU);
2295 DV.SetLinearForm(-(V4 + DerPlateV).Dot(Normal), Normal, DerPlateV);
2296 Plate_PinpointConstraint PinU( P2d.XY(), DU.XYZ(), 1, 0 );
2297 Plate_PinpointConstraint PinV( P2d.XY(), DV.XYZ(), 0, 1 );
2298 myPlate.Load( PinU );
2299 myPlate.Load( PinV );
2302 // Loading of points G2
2304 { gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10;
2305 CC->D2(myPlateCont->Value(i).Value(j),
2307 mySurfInit->D2(P2d.Coord(1),P2d.Coord(2),PP,V3,V4,V8,V9,V10);
2309 Plate_D1 D1final(V1.XYZ(),V2.XYZ());
2310 Plate_D1 D1init(V3.XYZ(),V4.XYZ());
2311 Plate_D2 D2final (V5.XYZ(),V6.XYZ(),V7.XYZ());
2312 Plate_D2 D2init (V8.XYZ(),V9.XYZ(),V10.XYZ());
2315 Plate_GtoCConstraint GCC(P2d.XY(),D1init,D1final,D2init,D2final);
2318 // else // Good but too expansive
2320 // Plate_FreeGtoCConstraint FreeGCC( P2d.XY(),
2321 // D1init, D1final, D2init, D2final );
2322 // myPlate.Load( FreeGCC );
2332 //---------------------------------------------------------
2333 // Function : LoadPoint
2334 //---------------------------------------------------------
2335 //void GeomPlate_BuildPlateSurface::LoadPoint(const Standard_Integer NbBoucle,
2336 void GeomPlate_BuildPlateSurface::LoadPoint(const Standard_Integer ,
2337 const Standard_Integer OrderMax)
2341 Standard_Integer NTPntCont=myPntCont->Length();
2342 Standard_Integer Tang, i;
2343 // gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10;
2345 // Loading of points of point constraints
2346 for (i=1;i<=NTPntCont;i++) {
2347 myPntCont->Value(i)->D0(P3d);
2348 P2d = myPntCont->Value(i)->Pnt2dOnSurf();
2349 mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP);
2350 Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1),
2351 -PP.Coord(2)+P3d.Coord(2),
2352 -PP.Coord(3)+P3d.Coord(3));
2353 Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0);
2355 Tang = Min(myPntCont->Value(i)->Order(), OrderMax);
2356 if (Tang==1) {// ==1
2358 myPntCont->Value(i)->D1(PP,V1,V2);
2359 mySurfInit->D1(P2d.Coord(1),P2d.Coord(2),PP,V3,V4);
2360 Plate_D1 D1final(V1.XYZ(),V2.XYZ());
2361 Plate_D1 D1init(V3.XYZ(),V4.XYZ());
2364 Plate_GtoCConstraint GCC( P2d.XY(), D1init, D1final );
2365 myPlate.Load( GCC );
2369 Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), D1init, D1final );
2370 myPlate.Load( FreeGCC );
2373 // Loading of points G2 GeomPlate_PlateG0Criterion
2375 { gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10;
2376 myPntCont->Value(i)->D2(PP,V1,V2,V5,V6,V7);
2377 // gp_Vec Tv2 = V1^V2;
2378 mySurfInit->D2(P2d.Coord(1),P2d.Coord(2),PP,V3,V4,V8,V9,V10);
2379 Plate_D1 D1final(V1.XYZ(),V2.XYZ());
2380 Plate_D1 D1init(V3.XYZ(),V4.XYZ());
2381 Plate_D2 D2final (V5.XYZ(),V6.XYZ(),V7.XYZ());
2382 Plate_D2 D2init (V8.XYZ(),V9.XYZ(),V10.XYZ());
2385 Plate_GtoCConstraint GCC( P2d.XY(), D1init, D1final, D2init, D2final );
2386 myPlate.Load( GCC );
2388 // else // Good but too expansive
2390 // Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), D1init, D1final, D2init//, D2final );
2391 // myPlate.Load( FreeGCC );
2397 //---------------------------------------------------------
2398 // Function : VerifSurface
2399 //---------------------------------------------------------
2400 Standard_Boolean GeomPlate_BuildPlateSurface::
2401 VerifSurface(const Standard_Integer NbBoucle)
2403 //======================================================================
2405 //======================================================================
2406 Standard_Integer NTLinCont=myLinCont->Length();
2407 Standard_Boolean Result=Standard_True;
2409 // variable for error calculation
2410 myG0Error=0,myG1Error =0, myG2Error=0;
2412 for (Standard_Integer i=1;i<=NTLinCont;i++) {
2413 Handle(GeomPlate_CurveConstraint) LinCont;
2414 LinCont = myLinCont->Value(i);
2415 if (LinCont->Order()!=-1) {
2416 Standard_Integer NbPts_i = myParCont->Value(i).Length();
2419 Handle(TColStd_HArray1OfReal) tdist =
2420 new TColStd_HArray1OfReal(1,NbPts_i-1);
2421 Handle(TColStd_HArray1OfReal) tang =
2422 new TColStd_HArray1OfReal(1,NbPts_i-1);
2423 Handle(TColStd_HArray1OfReal) tcourb =
2424 new TColStd_HArray1OfReal(1,NbPts_i-1);
2426 EcartContraintesMil (i,tdist,tang,tcourb);
2428 Standard_Real diffDistMax=0,SdiffDist=0;
2429 Standard_Real diffAngMax=0,SdiffAng=0;
2430 Standard_Integer NdiffDist=0,NdiffAng=0;
2433 for (Standard_Integer j=1;j<NbPts_i;j++)
2434 { if (tdist->Value(j)>myG0Error)
2435 myG0Error=tdist->Value(j);
2436 if (tang->Value(j)>myG1Error)
2437 myG1Error=tang->Value(j);
2438 if (tcourb->Value(j)>myG2Error)
2439 myG2Error=tcourb->Value(j);
2441 if (myParCont->Value(i).Length()>3)
2442 U=(myParCont->Value(i).Value(j)+myParCont->Value(i).Value(j+1))/2;
2444 U=LinCont->FirstParameter()+
2445 ( LinCont->LastParameter()-LinCont->FirstParameter())*(j-1)/(NbPts_i-2);
2446 Standard_Real diffDist=tdist->Value(j)-LinCont->G0Criterion(U),diffAng;
2447 if (LinCont->Order()>0)
2448 diffAng=tang->Value(j)-LinCont->G1Criterion(U);
2450 // find the maximum variation of error and calculate the average
2452 diffDist = diffDist/LinCont->G0Criterion(U);
2453 if (diffDist>diffDistMax)
2454 diffDistMax = diffDist;
2455 SdiffDist+=diffDist;
2458 if ((Affich) && (NbBoucle == myNbIter)) {
2462 Handle(Draw_Marker3D) mark =
2463 new (Draw_Marker3D)(P, Draw_X, Draw_orange);
2465 sprintf(name,"mark_%d",++NbMark);
2466 Draw::Set(name, mark);
2467 if (!LinCont->ProjectedCurve().IsNull())
2468 P2d = LinCont->ProjectedCurve()->Value(U);
2470 { if (!LinCont->Curve2dOnSurf().IsNull())
2471 P2d = LinCont->Curve2dOnSurf()->Value(U);
2473 P2d = ProjectPoint(P);
2475 sprintf(name,"mark2d_%d",++NbMark);
2476 Handle(Draw_Marker2D) mark2d =
2477 new (Draw_Marker2D)(P2d, Draw_X, Draw_orange);
2478 Draw::Set(name, mark2d);
2483 if ((diffAng>0)&&(LinCont->Order()==1)) {
2484 diffAng=diffAng/myLinCont->Value(i)->G1Criterion(U);
2485 if (diffAng>diffAngMax)
2486 diffAngMax = diffAng;
2490 if ((Affich) && (NbBoucle == myNbIter)) {
2493 Handle(Draw_Marker3D) mark =
2494 new Draw_Marker3D(P, Draw_X, Draw_or);
2496 sprintf(name,"mark_%d",++NbMark);
2497 Draw::Set(name, mark);
2503 if (NdiffDist>0) {// at least one point is not acceptable in G0
2505 if(LinCont->Order()== 0)
2506 Coef = 0.6*Log(diffDistMax+7.4);
2507 //7.4 corresponds to the calculation of min. coefficient = 1.2 is e^1.2/0.6
2509 Coef = Log(diffDistMax+3.3);
2510 //3.3 corresponds to calculation of min. coefficient = 1.2 donc e^1.2
2513 //experimentally after the coefficient becomes bad for L cases
2514 if ((NbBoucle>1)&&(diffDistMax>2))
2518 if (LinCont->NbPoints()>=Floor(LinCont->NbPoints()*Coef))
2519 Coef=2;// to provide increase of the number of points
2521 LinCont->SetNbPoints(Standard_Integer(LinCont->NbPoints() * Coef));
2522 Result=Standard_False;
2525 if (NdiffAng>0) // at least 1 point is not acceptable in G1
2526 { Standard_Real Coef=1.5;
2527 if ((LinCont->NbPoints()+1)>=Floor(LinCont->NbPoints()*Coef))
2530 LinCont->SetNbPoints(Standard_Integer(LinCont->NbPoints()*Coef )) ;
2531 Result=Standard_False;
2537 if (myFree && NbBoucle == 1)
2538 myPrevPlate = myPlate;
2546 //---------------------------------------------------------
2547 // Function : VerifPoint
2548 //---------------------------------------------------------
2549 void GeomPlate_BuildPlateSurface::
2550 VerifPoints (Standard_Real& Dist,
2552 Standard_Real& Curv) const
2553 { Standard_Integer NTPntCont=myPntCont->Length();
2556 gp_Vec v1i, v1f, v2i, v2f, v3i, v3f;
2557 Ang=0;Dist=0,Curv=0;
2558 Handle(GeomPlate_PointConstraint) PntCont;
2559 for (Standard_Integer i=1;i<=NTPntCont;i++) {
2560 PntCont = myPntCont->Value(i);
2561 switch (PntCont->Order())
2563 P2d = PntCont->Pnt2dOnSurf();
2565 myGeomPlateSurface->D0( P2d.Coord(1), P2d.Coord(2), Pf);
2566 Dist = Pf.Distance(Pi);
2569 PntCont->D1(Pi, v1i, v2i);
2570 P2d = PntCont->Pnt2dOnSurf();
2571 myGeomPlateSurface->D1( P2d.Coord(1), P2d.Coord(2), Pf, v1f, v2f);
2572 Dist = Pf.Distance(Pi);
2573 v3i = v1i^v2i; v3f=v1f^v2f;
2579 Handle(Geom_Surface) Splate (myGeomPlateSurface);
2580 LocalAnalysis_SurfaceContinuity CG2;
2581 P2d = PntCont->Pnt2dOnSurf();
2582 GeomLProp_SLProps Prop (Splate, P2d.Coord(1), P2d.Coord(2),
2584 CG2.ComputeAnalysis(Prop, PntCont->LPropSurf(),GeomAbs_G2);
2587 Curv=CG2.G2CurvatureGap();
2593 Standard_Real GeomPlate_BuildPlateSurface::ComputeAnisotropie() const
2603 Standard_Boolean GeomPlate_BuildPlateSurface::IsOrderG1() const
2605 Standard_Boolean result = Standard_True;
2606 for (Standard_Integer i = 1; i <= myLinCont->Length(); i++)
2607 if (myLinCont->Value(i)->Order() < 1)
2609 result = Standard_False;