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
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
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.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 #include <GeomPlate_BuildPlateSurface.hxx>
18
19 #include <Adaptor2d_Curve2d.hxx>
20 #include <Adaptor3d_Curve.hxx>
21 #include <Adaptor3d_CurveOnSurface.hxx>
22 #include <Approx_CurveOnSurface.hxx>
23 #include <Extrema_ExtPS.hxx>
24 #include <Extrema_POnSurf.hxx>
25 #include <GCPnts_AbscissaPoint.hxx>
26 #include <Geom2d_BezierCurve.hxx>
27 #include <Geom2d_BSplineCurve.hxx>
28 #include <Geom2d_Curve.hxx>
29 #include <Geom2dAdaptor_Curve.hxx>
30 #include <Geom2dInt_GInter.hxx>
31 #include <Geom_BSplineSurface.hxx>
32 #include <Geom_Plane.hxx>
33 #include <Geom_RectangularTrimmedSurface.hxx>
34 #include <Geom_Surface.hxx>
35 #include <GeomAdaptor.hxx>
36 #include <GeomAdaptor_Surface.hxx>
37 #include <GeomLProp_SLProps.hxx>
38 #include <GeomPlate_BuildAveragePlane.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>
45 #include <gp_Pnt.hxx>
46 #include <gp_Pnt2d.hxx>
47 #include <gp_Vec.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_ProgressScope.hxx>
70
71 #include <stdio.h>
72
73 #ifdef DRAW
74 #include <DrawTrSurf.hxx>
75 #include <Draw_Marker3D.hxx>
76 #include <Draw_Marker2D.hxx>
77 #include <Draw.hxx>
78 // 0 : No display
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;
86 #endif
87
88 #ifdef OCCT_DEBUG
89 #include <OSD_Chronometer.hxx>
90 #include <Geom_Surface.hxx>
91 static Standard_Integer Affich=0;
92 #endif
93
94 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
95 // =========================================================
96 // C O N S T R U C T O R S
97 // =========================================================
98 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
99
100 //---------------------------------------------------------
101 // Constructor compatible with the old version
102 //---------------------------------------------------------
GeomPlate_BuildPlateSurface(const Handle (TColStd_HArray1OfInteger)& NPoints,const Handle (GeomPlate_HArray1OfHCurve)& TabCurve,const Handle (TColStd_HArray1OfInteger)& Tang,const Standard_Integer Degree,const Standard_Integer NbIter,const Standard_Real Tol2d,const Standard_Real Tol3d,const Standard_Real TolAng,const Standard_Real,const Standard_Boolean Anisotropie)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
114 ) :
115 myAnisotropie(Anisotropie),
116 myDegree(Degree),
117 myNbIter(NbIter),
118 myProj(),
119 myTol2d(Tol2d),
120 myTol3d(Tol3d),
121 myTolAng(TolAng),
122 myNbBounds(0)
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;
127 if (myNbIter<1)
128 throw Standard_ConstructionError("GeomPlate : Number of iteration must be >= 1");
129
130 if (NTCurve==0)
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;
135 Standard_Integer i ;
136 for ( i=1;i<=NTCurve;i++)
137 { nbp+=NPoints->Value(i);
138 }
139 if (nbp==0)
140 throw Standard_ConstructionError("GeomPlate : the resolution is impossible if the number of constraints points is 0");
141 if (myDegree<2)
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),
147 Tang->Value(i),
148 NPoints->Value(i));
149 myLinCont->Append(Cont);
150 }
151 mySurfInitIsGive=Standard_False;
152
153 myIsLinear = Standard_True;
154 myFree = Standard_False;
155 }
156
157 //------------------------------------------------------------------
158 // Constructor with initial surface and degree
159 //------------------------------------------------------------------
GeomPlate_BuildPlateSurface(const Handle (Geom_Surface)& Surf,const Standard_Integer Degree,const Standard_Integer NbPtsOnCur,const Standard_Integer NbIter,const Standard_Real Tol2d,const Standard_Real Tol3d,const Standard_Real TolAng,const Standard_Real,const Standard_Boolean Anisotropie)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 ) :
170 mySurfInit(Surf),
171 myAnisotropie(Anisotropie),
172 myDegree(Degree),
173 myNbPtsOnCur(NbPtsOnCur),
174 myNbIter(NbIter),
175 myProj(),
176 myTol2d(Tol2d),
177 myTol3d(Tol3d),
178 myTolAng(TolAng),
179 myNbBounds(0)
180 { if (myNbIter<1)
181 throw Standard_ConstructionError("GeomPlate : Number of iteration must be >= 1");
182 if (myDegree<2)
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;
187
188 myIsLinear = Standard_True;
189 myFree = Standard_False;
190 }
191
192
193 //---------------------------------------------------------
194 // Constructor with degree
195 //---------------------------------------------------------
GeomPlate_BuildPlateSurface(const Standard_Integer Degree,const Standard_Integer NbPtsOnCur,const Standard_Integer NbIter,const Standard_Real Tol2d,const Standard_Real Tol3d,const Standard_Real TolAng,const Standard_Real,const Standard_Boolean Anisotropie)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),
206 myDegree(Degree),
207 myNbPtsOnCur(NbPtsOnCur),
208 myNbIter(NbIter),
209 myProj(),
210 myTol2d(Tol2d),
211 myTol3d(Tol3d),
212 myTolAng(TolAng),
213 myNbBounds(0)
214 { if (myNbIter<1)
215 throw Standard_ConstructionError("GeomPlate : Number of iteration must be >= 1");
216 if (myDegree<2)
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;
221
222 myIsLinear = Standard_True;
223 myFree = Standard_False;
224 }
225
226 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
227 // =========================================================
228 // P U B L I C M E T H O D S
229 // =========================================================
230 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
231
232
233
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 //-------------------------------------------------------------------------
TrierTab(Handle (TColStd_HArray1OfInteger)& Tab)247 static void TrierTab(Handle(TColStd_HArray1OfInteger)& Tab)
248 {
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;
255 }
256
257 //---------------------------------------------------------
258 // Function : ProjectCurve
259 //---------------------------------------------------------
Handle(Geom2d_Curve)260 Handle(Geom2d_Curve) GeomPlate_BuildPlateSurface::ProjectCurve(const Handle(Adaptor3d_Curve)& Curv)
261 {
262 // Project a curve on a plane
263 Handle(Geom2d_Curve) Curve2d ;
264 Handle(GeomAdaptor_Surface) hsur = new GeomAdaptor_Surface(mySurfInit);
265 gp_Pnt2d P2d;
266
267 Handle(ProjLib_HCompProjectedCurve) HProjector = new ProjLib_HCompProjectedCurve (hsur, Curv, myTol3d/10, myTol3d/10);
268
269 Standard_Real UdebCheck, UfinCheck, ProjUdeb, ProjUfin;
270 UdebCheck = Curv->FirstParameter();
271 UfinCheck = Curv->LastParameter();
272 HProjector->Bounds( 1, ProjUdeb, ProjUfin );
273
274 if (HProjector->NbCurves() != 1 ||
275 Abs( UdebCheck -ProjUdeb ) > Precision::PConfusion() ||
276 Abs( UfinCheck -ProjUfin ) > Precision::PConfusion())
277 {
278 if (HProjector->IsSinglePnt(1, P2d))
279 {
280 // solution in a point
281 TColgp_Array1OfPnt2d poles(1, 2);
282 poles.Init(P2d);
283 Curve2d = new (Geom2d_BezierCurve) (poles);
284 }
285 else
286 {
287 Curve2d.Nullify(); // No continuous solution
288 #ifdef OCCT_DEBUG
289 std::cout << "BuildPlateSurace :: No continuous projection" << std::endl;
290 #endif
291 }
292 }
293 else {
294 GeomAbs_Shape Continuity = GeomAbs_C1;
295 Standard_Integer MaxDegree = 10, MaxSeg;
296 Standard_Real Udeb, Ufin;
297 HProjector->Bounds(1, Udeb, Ufin);
298
299 MaxSeg = 20 + HProjector->NbIntervals(GeomAbs_C3);
300 Approx_CurveOnSurface appr(HProjector, hsur, Udeb, Ufin, myTol3d);
301 appr.Perform(MaxSeg, MaxDegree, Continuity, Standard_False, Standard_True);
302
303 Curve2d = appr.Curve2d();
304 }
305 #ifdef DRAW
306 if (Affich) {
307 char name[256];
308 sprintf(name,"proj_%d",++NbProj);
309 DrawTrSurf::Set(name, Curve2d);
310 }
311 #endif
312 return Curve2d;
313 }
314 //---------------------------------------------------------
315 // Function : ProjectedCurve
316 //---------------------------------------------------------
Handle(Adaptor2d_Curve2d)317 Handle(Adaptor2d_Curve2d) GeomPlate_BuildPlateSurface::ProjectedCurve( Handle(Adaptor3d_Curve)& Curv)
318 {
319 // Projection of a curve on the initial surface
320
321 Handle(GeomAdaptor_Surface) hsur = new GeomAdaptor_Surface(mySurfInit);
322
323 Handle(ProjLib_HCompProjectedCurve) HProjector = new ProjLib_HCompProjectedCurve (hsur, Curv, myTolU/10, myTolV/10);
324 if (HProjector->NbCurves() != 1)
325 {
326 HProjector.Nullify(); // No continuous solution
327 #ifdef OCCT_DEBUG
328 std::cout << "BuildPlateSurace :: No continuous projection" << std::endl;
329 #endif
330 }
331 else
332 {
333 Standard_Real First1,Last1,First2,Last2;
334 First1=Curv->FirstParameter();
335 Last1 =Curv->LastParameter();
336 HProjector->Bounds(1,First2,Last2);
337
338 if (Abs(First1-First2) <= Max(myTolU,myTolV) &&
339 Abs(Last1-Last2) <= Max(myTolU,myTolV))
340 {
341 HProjector = Handle(ProjLib_HCompProjectedCurve)::DownCast(HProjector->Trim(First2,Last2,Precision::PConfusion()));
342 }
343 else
344 {
345 HProjector.Nullify(); // No continuous solution
346 #ifdef OCCT_DEBUG
347 std::cout << "BuildPlateSurace :: No complete projection" << std::endl;
348 #endif
349 }
350 }
351 return HProjector;
352 }
353
354 //---------------------------------------------------------
355 // Function : ProjectPoint
356 //---------------------------------------------------------
357 // Projects a point on the initial surface
358 //---------------------------------------------------------
ProjectPoint(const gp_Pnt & p3d)359 gp_Pnt2d GeomPlate_BuildPlateSurface::ProjectPoint(const gp_Pnt &p3d)
360 { Extrema_POnSurf P;
361 myProj.Perform(p3d);
362 Standard_Integer nearest = 1;
363 if( myProj.NbExt() > 1)
364 {
365 Standard_Real dist2mini = myProj.SquareDistance(1);
366 for (Standard_Integer i=2; i<=myProj.NbExt();i++)
367 {
368 if (myProj.SquareDistance(i) < dist2mini)
369 {
370 dist2mini = myProj.SquareDistance(i);
371 nearest = i;
372 }
373 }
374 }
375 P=myProj.Point(nearest);
376 Standard_Real u,v;
377 P.Parameter(u,v);
378 gp_Pnt2d p2d;
379 p2d.SetCoord(u,v);
380
381 return p2d;
382 }
383 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
384 // =========================================================
385 // P U B L I C M E T H O D S
386 // =========================================================
387 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
388 //---------------------------------------------------------
389 // Function : Init
390 //---------------------------------------------------------
391 // Initializes linear and point constraints
392 //---------------------------------------------------------
Init()393 void GeomPlate_BuildPlateSurface::Init()
394 { myLinCont->Clear();
395 myPntCont->Clear();
396 myPntCont = new GeomPlate_HSequenceOfPointConstraint;
397 myLinCont = new GeomPlate_HSequenceOfCurveConstraint;
398 }
399
400 //---------------------------------------------------------
401 // Function : LoadInitSurface
402 //---------------------------------------------------------
403 // Loads the initial surface
404 //---------------------------------------------------------
LoadInitSurface(const Handle (Geom_Surface)& Surf)405 void GeomPlate_BuildPlateSurface::LoadInitSurface(const Handle(Geom_Surface)& Surf)
406 { mySurfInit = Surf;
407 mySurfInitIsGive=Standard_True;
408 }
409
410 //---------------------------------------------------------
411 // Function : Add
412 //---------------------------------------------------------
413 //---------------------------------------------------------
414 // Adds a linear constraint
415 //---------------------------------------------------------
416 void GeomPlate_BuildPlateSurface::
Add(const Handle (GeomPlate_CurveConstraint)& Cont)417 Add(const Handle(GeomPlate_CurveConstraint)& Cont)
418 { myLinCont->Append(Cont);
419 }
420
SetNbBounds(const Standard_Integer NbBounds)421 void GeomPlate_BuildPlateSurface::SetNbBounds( const Standard_Integer NbBounds )
422 {
423 myNbBounds = NbBounds;
424 }
425
426 //---------------------------------------------------------
427 // Function : Add
428 //---------------------------------------------------------
429 //---------------------------------------------------------
430 // Adds a point constraint
431 //---------------------------------------------------------
432 void GeomPlate_BuildPlateSurface::
Add(const Handle (GeomPlate_PointConstraint)& Cont)433 Add(const Handle(GeomPlate_PointConstraint)& Cont)
434 {
435 myPntCont->Append(Cont);
436 }
437
438 //---------------------------------------------------------
439 // Function : Perform
440 // Calculates the surface filled with loaded constraints
441 //---------------------------------------------------------
Perform(const Message_ProgressRange & theProgress)442 void GeomPlate_BuildPlateSurface::Perform(const Message_ProgressRange& theProgress)
443 {
444 #ifdef OCCT_DEBUG
445 // Timing
446 OSD_Chronometer Chrono;
447 Chrono.Reset();
448 Chrono.Start();
449 #endif
450
451 if (myNbBounds == 0)
452 myNbBounds = myLinCont->Length();
453
454 myPlate.Init();
455 //=====================================================================
456 // Declaration of variables.
457 //=====================================================================
458 Standard_Integer NTLinCont = myLinCont->Length(),
459 NTPntCont = myPntCont->Length(), NbBoucle=0;
460 Standard_Boolean Fini=Standard_True;
461 if ((NTLinCont + NTPntCont) == 0)
462 {
463 #ifdef OCCT_DEBUG
464 std::cout << "WARNING : GeomPlate : The number of constraints is null." << std::endl;
465 #endif
466
467 return;
468 }
469
470 //======================================================================
471 // Initial Surface
472 //======================================================================
473 Message_ProgressScope aPS(theProgress, "Calculating the surface filled", 100, Standard_True);
474 if (!mySurfInitIsGive)
475 {
476 ComputeSurfInit (aPS.Next(10));
477 if (aPS.UserBreak())
478 return;
479 }
480
481 else {
482 if (NTLinCont>=2)
483 { // Table of transformations to preserve the initial order, see TrierTab
484 myInitOrder = new TColStd_HArray1OfInteger(1,NTLinCont);
485 for (Standard_Integer l=1;l<=NTLinCont;l++)
486 myInitOrder->SetValue(l,l);
487 if (!CourbeJointive(myTol3d))
488 {// throw Standard_Failure("Curves are not joined");
489 #ifdef OCCT_DEBUG
490 std::cout<<"WARNING : Courbes non jointives a "<<myTol3d<<" pres"<<std::endl;
491 #endif
492 }
493 TrierTab(myInitOrder); // Reorder the table of transformations
494 }
495 else if(NTLinCont > 0)//Patch
496 {
497 mySense = new TColStd_HArray1OfInteger( 1, NTLinCont, 0 );
498 myInitOrder = new TColStd_HArray1OfInteger( 1, NTLinCont, 1 );
499 }
500 }
501
502 if (mySurfInit.IsNull())
503 {
504 return;
505 }
506
507 Standard_Real u1,v1,u2,v2;
508 mySurfInit->Bounds(u1,v1,u2,v2);
509 GeomAdaptor_Surface aSurfInit(mySurfInit);
510 myTolU = aSurfInit.UResolution(myTol3d);
511 myTolV = aSurfInit.VResolution(myTol3d);
512 myProj.Initialize(aSurfInit,u1,v1,u2,v2,
513 myTolU,myTolV);
514
515 //======================================================================
516 // Projection of curves
517 //======================================================================
518 Standard_Boolean Ok = Standard_True;
519 for (Standard_Integer i = 1; i <= NTLinCont; i++)
520 if(myLinCont->Value(i)->Curve2dOnSurf().IsNull())
521 {
522 Handle( Geom2d_Curve ) Curve2d = ProjectCurve( myLinCont->Value(i)->Curve3d() );
523 if (Curve2d.IsNull())
524 {
525 Ok = Standard_False;
526 break;
527 }
528 myLinCont->ChangeValue(i)->SetCurve2dOnSurf( Curve2d );
529 }
530 if (!Ok)
531 {
532 GeomPlate_MakeApprox App(myGeomPlateSurface, myTol3d,
533 1, 3,
534 15*myTol3d,
535 -1, GeomAbs_C0,
536 1.3);
537 mySurfInit = App.Surface();
538
539 mySurfInit->Bounds(u1,v1,u2,v2);
540 GeomAdaptor_Surface Surf(mySurfInit);
541 myTolU = Surf.UResolution(myTol3d);
542 myTolV = Surf.VResolution(myTol3d);
543 myProj.Initialize(Surf,u1,v1,u2,v2,
544 myTolU,myTolV);
545
546 Ok = Standard_True;
547 for (Standard_Integer i = 1; i <= NTLinCont; i++)
548 {
549 Handle( Geom2d_Curve ) Curve2d = ProjectCurve( myLinCont->Value(i)->Curve3d() );
550 if (Curve2d.IsNull())
551 {
552 Ok = Standard_False;
553 break;
554 }
555 myLinCont->ChangeValue(i)->SetCurve2dOnSurf( Curve2d );
556 }
557 if (!Ok)
558 {
559 mySurfInit = myPlanarSurfInit;
560
561 mySurfInit->Bounds(u1,v1,u2,v2);
562 GeomAdaptor_Surface SurfNew(mySurfInit);
563 myTolU = SurfNew.UResolution(myTol3d);
564 myTolV = SurfNew.VResolution(myTol3d);
565 myProj.Initialize(SurfNew,u1,v1,u2,v2,
566 myTolU,myTolV);
567
568 for (Standard_Integer i = 1; i <= NTLinCont; i++)
569 myLinCont->ChangeValue(i)->
570 SetCurve2dOnSurf(ProjectCurve( myLinCont->Value(i)->Curve3d() ) );
571 }
572 else { // Project the points
573 for (Standard_Integer i=1; i<=NTPntCont; i++) {
574 gp_Pnt P;
575 myPntCont->Value(i)->D0(P);
576 myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P));
577 }
578 }
579 }
580
581 //======================================================================
582 // Projection of points
583 //======================================================================
584 for (Standard_Integer i=1;i<=NTPntCont;i++) {
585 if (! myPntCont->Value(i)->HasPnt2dOnSurf()) {
586 gp_Pnt P;
587 myPntCont->Value(i)->D0(P);
588 myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P));
589 }
590 }
591
592 //======================================================================
593 // Number of points by curve
594 //======================================================================
595 if ((NTLinCont !=0)&&(myNbPtsOnCur!=0))
596 CalculNbPtsInit();
597
598 //======================================================================
599 // Management of incompatibilites between curves
600 //======================================================================
601 Handle(GeomPlate_HArray1OfSequenceOfReal) PntInter;
602 Handle(GeomPlate_HArray1OfSequenceOfReal) PntG1G1;
603 if (NTLinCont != 0) {
604 PntInter = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
605 PntG1G1 = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
606 Intersect(PntInter, PntG1G1);
607 }
608
609 //======================================================================
610 // Loop to obtain a better surface
611 //======================================================================
612
613 myFree = !myIsLinear;
614
615 do
616 {
617 #ifdef OCCT_DEBUG
618 if (Affich && NbBoucle) {
619 std::cout<<"Resultats boucle"<< NbBoucle << std::endl;
620 std::cout<<"DistMax="<<myG0Error<<std::endl;
621 if (myG1Error!=0)
622 std::cout<<"AngleMax="<<myG1Error<<std::endl;
623 if (myG2Error!=0)
624 std::cout<<"CourbMax="<<myG2Error<<std::endl;
625 }
626 #endif
627 NbBoucle++;
628 if (NTLinCont!=0)
629 { //====================================================================
630 // Calculate the total number of points and the maximum of points by curve
631 //====================================================================
632 Standard_Integer NPointMax=0;
633 for (Standard_Integer i=1;i<=NTLinCont;i++)
634 if ((myLinCont->Value(i)->NbPoints())>NPointMax)
635 NPointMax=(Standard_Integer )( myLinCont->Value(i)->NbPoints());
636 //====================================================================
637 // Discretization of curves
638 //====================================================================
639 Discretise(PntInter, PntG1G1);
640 //====================================================================
641 // Preparation of constraint points for plate
642 //====================================================================
643 LoadCurve( NbBoucle );
644 if ( myPntCont->Length() != 0)
645 LoadPoint( NbBoucle );
646 //====================================================================
647 // Construction of the surface
648 //====================================================================
649
650 myPlate.SolveTI(myDegree, ComputeAnisotropie(), aPS.Next(90));
651
652 if (aPS.UserBreak())
653 {
654 return;
655 }
656
657 if (!myPlate.IsDone())
658 {
659 #ifdef OCCT_DEBUG
660 std::cout << "WARNING : GeomPlate : calculation of Plate failed" << std::endl;
661 #endif
662 return;
663 }
664
665 myGeomPlateSurface = new GeomPlate_Surface(mySurfInit,myPlate);
666 Standard_Real Umin,Umax,Vmin,Vmax;
667 myPlate.UVBox(Umin,Umax,Vmin,Vmax);
668 myGeomPlateSurface->SetBounds(Umin,Umax,Vmin,Vmax);
669
670 Fini = VerifSurface(NbBoucle);
671 if ((NbBoucle >= myNbIter)&&(!Fini))
672 {
673 #ifdef OCCT_DEBUG
674 std::cout<<"Warning: objective was not reached"<<std::endl;
675 #endif
676 Fini = Standard_True;
677 }
678
679 if ((NTPntCont != 0)&&(Fini))
680 { Standard_Real di,an,cu;
681 VerifPoints(di,an,cu);
682 }
683 }
684 else
685 { LoadPoint( NbBoucle );
686 //====================================================================
687 //Construction of the surface
688 //====================================================================
689 myPlate.SolveTI(myDegree, ComputeAnisotropie(), aPS.Next(90));
690
691 if (aPS.UserBreak())
692 {
693 return;
694 }
695
696 if (!myPlate.IsDone())
697 {
698 #ifdef OCCT_DEBUG
699 std::cout << "WARNING : GeomPlate : calculation of Plate failed" << std::endl;
700 #endif
701 return;
702 }
703
704 myGeomPlateSurface = new GeomPlate_Surface(mySurfInit,myPlate);
705 Standard_Real Umin,Umax,Vmin,Vmax;
706 myPlate.UVBox(Umin,Umax,Vmin,Vmax);
707 myGeomPlateSurface->SetBounds(Umin,Umax,Vmin,Vmax);
708 Fini = Standard_True;
709 Standard_Real di,an,cu;
710 VerifPoints(di,an,cu);
711 }
712 } while (!Fini); // End loop for better surface
713 #ifdef OCCT_DEBUG
714 if (NTLinCont != 0)
715 { std::cout<<"======== Global results ==========="<<std::endl;
716 std::cout<<"DistMax="<<myG0Error<<std::endl;
717 if (myG1Error!=0)
718 std::cout<<"AngleMax="<<myG1Error<<std::endl;
719 if (myG2Error!=0)
720 std::cout<<"CourbMax="<<myG2Error<<std::endl;
721 }
722 Chrono.Stop();
723 Standard_Real Tps;
724 Chrono.Show(Tps);
725 std::cout<<"*** END OF GEOMPLATE ***"<<std::endl;
726 std::cout<<"Time of calculation : "<<Tps<<std::endl;
727 std::cout<<"Number of loops : "<<NbBoucle<<std::endl;
728 #endif
729 }
730
731 //---------------------------------------------------------
732 // Function : EcartContraintesMIL
733 //---------------------------------------------------------
734 void GeomPlate_BuildPlateSurface::
EcartContraintesMil(const Standard_Integer c,Handle (TColStd_HArray1OfReal)& d,Handle (TColStd_HArray1OfReal)& an,Handle (TColStd_HArray1OfReal)& courb)735 EcartContraintesMil ( const Standard_Integer c,
736 Handle(TColStd_HArray1OfReal)& d,
737 Handle(TColStd_HArray1OfReal)& an,
738 Handle(TColStd_HArray1OfReal)& courb)
739 {
740 Standard_Integer NbPt=myParCont->Value(c).Length();
741 Standard_Real U;
742 if (NbPt<3)
743 NbPt=4;
744 else
745 NbPt=myParCont->Value(c).Length();
746 gp_Vec v1i, v1f, v2i, v2f, v3i, v3f;
747 gp_Pnt Pi, Pf;
748 gp_Pnt2d P2d;Standard_Integer i;
749 Handle(GeomPlate_CurveConstraint) LinCont = myLinCont->Value(c);
750 switch (LinCont->Order())
751 { case 0 :
752 for (i=1; i<NbPt; i++)
753 { U = ( myParCont->Value(c).Value(i) +
754 myParCont->Value(c).Value(i+1) )/2;
755 LinCont->D0(U, Pi);
756 if (!LinCont->ProjectedCurve().IsNull())
757 P2d = LinCont->ProjectedCurve()->Value(U);
758 else
759 { if (!LinCont->Curve2dOnSurf().IsNull())
760 P2d = LinCont->Curve2dOnSurf()->Value(U);
761 else
762 P2d = ProjectPoint(Pi);
763 }
764 myGeomPlateSurface->D0( P2d.Coord(1), P2d.Coord(2), Pf);
765 an->Init(0);
766 courb->Init(0);
767 d->ChangeValue(i) = Pf.Distance(Pi);
768 }
769 break;
770 case 1 :
771 for (i=1; i<NbPt; i++)
772 { U = ( myParCont->Value(c).Value(i) +
773 myParCont->Value(c).Value(i+1) )/2;
774 LinCont->D1(U, Pi, v1i, v2i);
775 if (!LinCont->ProjectedCurve().IsNull())
776 P2d = LinCont->ProjectedCurve()->Value(U);
777 else
778 { if (!LinCont->Curve2dOnSurf().IsNull())
779 P2d = LinCont->Curve2dOnSurf()->Value(U);
780 else
781 P2d = ProjectPoint(Pi);
782 }
783 myGeomPlateSurface->D1( P2d.Coord(1), P2d.Coord(2), Pf, v1f, v2f);
784 d->ChangeValue(i) = Pf.Distance(Pi);
785 v3i = v1i^v2i; v3f=v1f^v2f;
786 Standard_Real angle=v3f.Angle(v3i);
787 if (angle>(M_PI/2))
788 an->ChangeValue(i) = M_PI -angle;
789 else
790 an->ChangeValue(i) = angle;
791 courb->Init(0);
792 }
793 break;
794 case 2 :
795 Handle(Geom_Surface) Splate (myGeomPlateSurface);
796 LocalAnalysis_SurfaceContinuity CG2;
797 for (i=1; i<NbPt; i++)
798 { U = ( myParCont->Value(c).Value(i) +
799 myParCont->Value(c).Value(i+1) )/2;
800 LinCont->D0(U, Pi);
801 if (!LinCont->ProjectedCurve().IsNull())
802 P2d = LinCont->ProjectedCurve()->Value(U);
803 else
804 { if (!LinCont->Curve2dOnSurf().IsNull())
805 P2d = LinCont->Curve2dOnSurf()->Value(U);
806 else
807 P2d = ProjectPoint(Pi);
808 }
809 GeomLProp_SLProps Prop (Splate, P2d.Coord(1), P2d.Coord(2),
810 2, 0.001);
811 CG2.ComputeAnalysis(Prop, myLinCont->Value(c)->LPropSurf(U),
812 GeomAbs_G2);
813 d->ChangeValue(i)=CG2.C0Value();
814 an->ChangeValue(i)=CG2.G1Angle();
815 courb->ChangeValue(i)=CG2.G2CurvatureGap();
816 }
817 break;
818 }
819 }
820
821
822
823 //---------------------------------------------------------
824 // Function : Disc2dContour
825 //---------------------------------------------------------
826 void GeomPlate_BuildPlateSurface::
Disc2dContour(const Standard_Integer,TColgp_SequenceOfXY & Seq2d)827 Disc2dContour ( const Standard_Integer /*nbp*/,
828 TColgp_SequenceOfXY& Seq2d)
829 {
830 #ifdef OCCT_DEBUG
831 if (Seq2d.Length()!=4)
832 std::cout<<"Number of points should be equal to 4 for Disc2dContour"<<std::endl;
833 #endif
834 // initialization
835 Seq2d.Clear();
836
837 // sampling in "cosine" + 3 points on each interval
838 Standard_Integer NTCurve = myLinCont->Length();
839 Standard_Integer NTPntCont = myPntCont->Length();
840 gp_Pnt2d P2d;
841 gp_XY UV;
842 gp_Pnt PP;
843 Standard_Real u1,v1,u2,v2;
844 Standard_Integer i ;
845
846 mySurfInit->Bounds(u1,v1,u2,v2);
847 GeomAdaptor_Surface Surf(mySurfInit);
848 myProj.Initialize(Surf,u1,v1,u2,v2,
849 myTolU,myTolV);
850
851 for( i=1; i<=NTPntCont; i++)
852 if (myPntCont->Value(i)->Order()!=-1)
853 { P2d = myPntCont->Value(i)->Pnt2dOnSurf();
854 UV.SetX(P2d.Coord(1));
855 UV.SetY(P2d.Coord(2));
856 Seq2d.Append(UV);
857 }
858 for(i=1; i<=NTCurve; i++)
859 {
860 Handle(GeomPlate_CurveConstraint) LinCont = myLinCont->Value(i);
861 if (LinCont->Order()!=-1)
862 { Standard_Integer NbPt=myParCont->Value(i).Length();
863 // first point of constraint (j=0)
864 if (!LinCont->ProjectedCurve().IsNull())
865 P2d = LinCont->ProjectedCurve()->Value(myParCont->Value(i).Value(1));
866
867 else
868 { if (!LinCont->Curve2dOnSurf().IsNull())
869 P2d = LinCont->Curve2dOnSurf()->Value(myParCont->Value(i).Value(1));
870
871 else
872 {
873 LinCont->D0(myParCont->Value(i).Value(1),PP);
874 P2d = ProjectPoint(PP);
875 }
876 }
877
878 UV.SetX(P2d.Coord(1));
879 UV.SetY(P2d.Coord(2));
880 Seq2d.Append(UV);
881 for (Standard_Integer j=2; j<NbPt; j++)
882 { Standard_Real Uj=myParCont->Value(i).Value(j),
883 Ujp1=myParCont->Value(i).Value(j+1);
884 if (!LinCont->ProjectedCurve().IsNull())
885 P2d = LinCont->ProjectedCurve()->Value((Ujp1+3*Uj)/4);
886
887 else
888 { if (!LinCont->Curve2dOnSurf().IsNull())
889 P2d = LinCont->Curve2dOnSurf()->Value((Ujp1+3*Uj)/4);
890
891 else
892 {
893 LinCont->D0((Ujp1+3*Uj)/4,PP);
894 P2d = ProjectPoint(PP);
895 }
896 }
897 UV.SetX(P2d.Coord(1));
898 UV.SetY(P2d.Coord(2));
899 Seq2d.Append(UV);
900 // point 1/2 previous
901 if (!LinCont->ProjectedCurve().IsNull())
902 P2d = LinCont->ProjectedCurve()->Value((Ujp1+Uj)/2);
903
904 else
905 { if (!LinCont->Curve2dOnSurf().IsNull())
906 P2d = LinCont->Curve2dOnSurf()->Value((Ujp1+Uj)/2);
907
908 else
909 {
910 LinCont->D0((Ujp1+Uj)/2,PP);
911 P2d = ProjectPoint(PP);
912 }
913 }
914
915 UV.SetX(P2d.Coord(1));
916 UV.SetY(P2d.Coord(2));
917 Seq2d.Append(UV);
918 // point 3/4 previous
919 if (!LinCont->ProjectedCurve().IsNull())
920 P2d = LinCont->ProjectedCurve()->Value((3*Ujp1+Uj)/4);
921
922 else
923 { if (!LinCont->Curve2dOnSurf().IsNull())
924 P2d = LinCont->Curve2dOnSurf()->Value((3*Ujp1+Uj)/4);
925
926 else
927 {
928 LinCont->D0((3*Ujp1+Uj)/4,PP);
929 P2d = ProjectPoint(PP);
930 }
931 }
932
933 UV.SetX(P2d.Coord(1));
934 UV.SetY(P2d.Coord(2));
935 Seq2d.Append(UV);
936 // current constraint point
937 if (!LinCont->ProjectedCurve().IsNull())
938 P2d = LinCont->ProjectedCurve()->Value(Ujp1);
939
940 else
941 { if (!LinCont->Curve2dOnSurf().IsNull())
942 P2d = LinCont->Curve2dOnSurf()->Value(Ujp1);
943
944 else
945 {
946 LinCont->D0(Ujp1,PP);
947 P2d = ProjectPoint(PP);
948 }
949 }
950
951 UV.SetX(P2d.Coord(1));
952 UV.SetY(P2d.Coord(2));
953 Seq2d.Append(UV);
954 }
955 }
956 }
957 }
958
959 //---------------------------------------------------------
960 // Function : Disc3dContour
961 //---------------------------------------------------------
962 void GeomPlate_BuildPlateSurface::
Disc3dContour(const Standard_Integer,const Standard_Integer iordre,TColgp_SequenceOfXYZ & Seq3d)963 Disc3dContour (const Standard_Integer /*nbp*/,
964 const Standard_Integer iordre,
965 TColgp_SequenceOfXYZ& Seq3d)
966 {
967 #ifdef OCCT_DEBUG
968 if (Seq3d.Length()!=4)
969 std::cout<<"nbp should be equal to 4 for Disc3dContour"<<std::endl;
970 if (iordre!=0&&iordre!=1)
971 std::cout<<"incorrect order for Disc3dContour"<<std::endl;
972 #endif
973 // initialization
974 Seq3d.Clear();
975 // sampling in "cosine" + 3 points on each interval
976 Standard_Real u1,v1,u2,v2;
977 mySurfInit->Bounds(u1,v1,u2,v2);
978 GeomAdaptor_Surface Surf(mySurfInit);
979 myProj.Initialize(Surf,u1,v1,u2,v2,
980 Surf.UResolution(myTol3d),
981 Surf.VResolution(myTol3d));
982 Standard_Integer NTCurve = myLinCont->Length();
983 Standard_Integer NTPntCont = myPntCont->Length();
984 // gp_Pnt2d P2d;
985 gp_Pnt P3d;
986 gp_Vec v1h,v2h,v3h;
987 gp_XYZ Pos;
988 Standard_Integer i ;
989
990 for( i=1; i<=NTPntCont; i++)
991 if (myPntCont->Value(i)->Order()!=-1)
992 { if (iordre==0)
993 { myPntCont->Value(i)->D0(P3d);
994 Pos.SetX(P3d.X());
995 Pos.SetY(P3d.Y());
996 Pos.SetZ(P3d.Z());
997 Seq3d.Append(Pos);
998 }
999 else {
1000 myPntCont->Value(i)->D1(P3d,v1h,v2h);
1001 v3h=v1h^v2h;
1002 Pos.SetX(v3h.X());
1003 Pos.SetY(v3h.Y());
1004 Pos.SetZ(v3h.Z());
1005 Seq3d.Append(Pos);
1006 }
1007 }
1008
1009 for(i=1; i<=NTCurve; i++)
1010 if (myLinCont->Value(i)->Order()!=-1)
1011
1012 { Standard_Integer NbPt=myParCont->Value(i).Length();
1013 // first constraint point (j=0)
1014 // Standard_Integer NbPt=myParCont->Length();
1015 if (iordre==0) {
1016 myLinCont->Value(i)->D0(myParCont->Value(i).Value(1),P3d);
1017 Pos.SetX(P3d.X());
1018 Pos.SetY(P3d.Y());
1019 Pos.SetZ(P3d.Z());
1020 Seq3d.Append(Pos);
1021 }
1022 else {
1023 myLinCont->Value(i)->D1(myParCont->Value(i).Value(1),P3d,v1h,v2h);
1024 v3h=v1h^v2h;
1025 Pos.SetX(v3h.X());
1026 Pos.SetY(v3h.Y());
1027 Pos.SetZ(v3h.Z());
1028 Seq3d.Append(Pos);
1029 }
1030
1031 for (Standard_Integer j=2; j<NbPt; j++)
1032 { Standard_Real Uj=myParCont->Value(i).Value(j),
1033 Ujp1=myParCont->Value(i).Value(j+1);
1034 if (iordre==0) {
1035 // point 1/4 previous
1036 myLinCont->Value(i)->D0((Ujp1+3*Uj)/4,P3d);
1037 Pos.SetX(P3d.X());
1038 Pos.SetY(P3d.Y());
1039 Pos.SetZ(P3d.Z());
1040 Seq3d.Append(Pos);
1041 // point 1/2 previous
1042 myLinCont->Value(i)->D0((Ujp1+Uj)/2,P3d);
1043 Pos.SetX(P3d.X());
1044 Pos.SetY(P3d.Y());
1045 Pos.SetZ(P3d.Z());
1046 Seq3d.Append(Pos);
1047 // point 3/4 previous
1048 myLinCont->Value(i)->D0((3*Ujp1+Uj)/4,P3d);
1049 Pos.SetX(P3d.X());
1050 Pos.SetY(P3d.Y());
1051 Pos.SetZ(P3d.Z());
1052 Seq3d.Append(Pos);
1053 // current constraint point
1054 myLinCont->Value(i)->D0(Ujp1,P3d);
1055 Pos.SetX(P3d.X());
1056 Pos.SetY(P3d.Y());
1057 Pos.SetZ(P3d.Z());
1058 Seq3d.Append(Pos);
1059 }
1060 else {
1061 // point 1/4 previous
1062 myLinCont->Value(i)->D1((Ujp1+3*Uj)/4,P3d,v1h,v2h);
1063 v3h=v1h^v2h;
1064 Pos.SetX(v3h.X());
1065 Pos.SetY(v3h.Y());
1066 Pos.SetZ(v3h.Z());
1067 Seq3d.Append(Pos);
1068 // point 1/2 previous
1069 myLinCont->Value(i)->D1((Ujp1+Uj)/2,P3d,v1h,v2h);
1070 v3h=v1h^v2h;
1071 Pos.SetX(v3h.X());
1072 Pos.SetY(v3h.Y());
1073 Pos.SetZ(v3h.Z());
1074 Seq3d.Append(Pos);
1075 // point 3/4 previous
1076 myLinCont->Value(i)->D1((3*Ujp1+Uj)/4,P3d,v1h,v2h);
1077 v3h=v1h^v2h;
1078 Pos.SetX(v3h.X());
1079 Pos.SetY(v3h.Y());
1080 Pos.SetZ(v3h.Z());
1081 Seq3d.Append(Pos);
1082 // current constraint point
1083 myLinCont->Value(i)->D1(Ujp1,P3d,v1h,v2h);
1084 v3h=v1h^v2h;
1085 Pos.SetX(v3h.X());
1086 Pos.SetY(v3h.Y());
1087 Pos.SetZ(v3h.Z());
1088 Seq3d.Append(Pos);
1089 }
1090 }
1091 }
1092
1093 }
1094
1095
1096 //---------------------------------------------------------
1097 // Function : IsDone
1098 //---------------------------------------------------------
IsDone() const1099 Standard_Boolean GeomPlate_BuildPlateSurface::IsDone() const
1100 { return myPlate.IsDone();
1101 }
1102
1103
1104
1105 //---------------------------------------------------------
1106 // Function : Surface
1107 //---------------------------------------------------------
Handle(GeomPlate_Surface)1108 Handle(GeomPlate_Surface) GeomPlate_BuildPlateSurface::Surface() const
1109 { return myGeomPlateSurface ;
1110 }
1111
1112 //---------------------------------------------------------
1113 // Function : SurfInit
1114 //---------------------------------------------------------
Handle(Geom_Surface)1115 Handle(Geom_Surface) GeomPlate_BuildPlateSurface::SurfInit() const
1116 { return mySurfInit ;
1117 }
1118
1119
1120 //---------------------------------------------------------
1121 // Function : Sense
1122 //---------------------------------------------------------
Handle(TColStd_HArray1OfInteger)1123 Handle(TColStd_HArray1OfInteger) GeomPlate_BuildPlateSurface::Sense() const
1124 { Standard_Integer NTCurve = myLinCont->Length();
1125 Handle(TColStd_HArray1OfInteger) Sens = new TColStd_HArray1OfInteger(1,
1126 NTCurve);
1127 for (Standard_Integer i=1; i<=NTCurve; i++)
1128 Sens->SetValue(i,mySense->Value(myInitOrder->Value(i)));
1129 return Sens;
1130 }
1131
1132
1133
1134 //---------------------------------------------------------
1135 // Function : Curve2d
1136 //---------------------------------------------------------
Handle(TColGeom2d_HArray1OfCurve)1137 Handle(TColGeom2d_HArray1OfCurve) GeomPlate_BuildPlateSurface::Curves2d() const
1138 { Standard_Integer NTCurve = myLinCont->Length();
1139 Handle(TColGeom2d_HArray1OfCurve) C2dfin =
1140 new TColGeom2d_HArray1OfCurve(1,NTCurve);
1141
1142 for (Standard_Integer i=1; i<=NTCurve; i++)
1143 C2dfin->SetValue(i,myLinCont->Value(myInitOrder->Value(i))->Curve2dOnSurf());
1144 return C2dfin;
1145
1146 }
1147
1148
1149 //---------------------------------------------------------
1150 // Function : Order
1151 //---------------------------------------------------------
Handle(TColStd_HArray1OfInteger)1152 Handle(TColStd_HArray1OfInteger) GeomPlate_BuildPlateSurface::Order() const
1153 { Handle(TColStd_HArray1OfInteger) result=
1154 new TColStd_HArray1OfInteger(1,myLinCont->Length());
1155 for (Standard_Integer i=1;i<=myLinCont->Length();i++)
1156 result->SetValue(myInitOrder->Value(i),i);
1157 return result;
1158 }
1159
1160
1161 //---------------------------------------------------------
1162 // Function : G0Error
1163 //---------------------------------------------------------
G0Error() const1164 Standard_Real GeomPlate_BuildPlateSurface::G0Error() const
1165 { return myG0Error;
1166 }
1167
1168 //---------------------------------------------------------
1169 // Function : G1Error
1170 //---------------------------------------------------------
G1Error() const1171 Standard_Real GeomPlate_BuildPlateSurface::G1Error() const
1172 { return myG1Error;
1173 }
1174
1175 //---------------------------------------------------------
1176 // Function : G2Error
1177 //---------------------------------------------------------
G2Error() const1178 Standard_Real GeomPlate_BuildPlateSurface::G2Error() const
1179 { return myG2Error;
1180 }
1181
1182 //=======================================================================
1183 //function : G0Error
1184 //purpose :
1185 //=======================================================================
1186
G0Error(const Standard_Integer Index)1187 Standard_Real GeomPlate_BuildPlateSurface::G0Error( const Standard_Integer Index )
1188 {
1189 Handle( TColStd_HArray1OfReal ) tdistance = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1190 Handle( TColStd_HArray1OfReal ) tangle = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1191 Handle( TColStd_HArray1OfReal ) tcurvature = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1192 EcartContraintesMil( Index, tdistance, tangle, tcurvature );
1193 Standard_Real MaxDistance = 0.;
1194 for (Standard_Integer i = 1; i <= myNbPtsOnCur; i++)
1195 if (tdistance->Value(i) > MaxDistance)
1196 MaxDistance = tdistance->Value(i);
1197 return MaxDistance;
1198 }
1199
1200 //=======================================================================
1201 //function : G1Error
1202 //purpose :
1203 //=======================================================================
1204
G1Error(const Standard_Integer Index)1205 Standard_Real GeomPlate_BuildPlateSurface::G1Error( const Standard_Integer Index )
1206 {
1207 Handle( TColStd_HArray1OfReal ) tdistance = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1208 Handle( TColStd_HArray1OfReal ) tangle = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1209 Handle( TColStd_HArray1OfReal ) tcurvature = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1210 EcartContraintesMil( Index, tdistance, tangle, tcurvature );
1211 Standard_Real MaxAngle = 0.;
1212 for (Standard_Integer i = 1; i <= myNbPtsOnCur; i++)
1213 if (tangle->Value(i) > MaxAngle)
1214 MaxAngle = tangle->Value(i);
1215 return MaxAngle;
1216 }
1217
1218 //=======================================================================
1219 //function : G2Error
1220 //purpose :
1221 //=======================================================================
1222
G2Error(const Standard_Integer Index)1223 Standard_Real GeomPlate_BuildPlateSurface::G2Error( const Standard_Integer Index )
1224 {
1225 Handle( TColStd_HArray1OfReal ) tdistance = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1226 Handle( TColStd_HArray1OfReal ) tangle = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1227 Handle( TColStd_HArray1OfReal ) tcurvature = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1228 EcartContraintesMil( Index, tdistance, tangle, tcurvature );
1229 Standard_Real MaxCurvature = 0.;
1230 for (Standard_Integer i = 1; i <= myNbPtsOnCur; i++)
1231 if (tcurvature->Value(i) > MaxCurvature)
1232 MaxCurvature = tcurvature->Value(i);
1233 return MaxCurvature;
1234 }
1235
Handle(GeomPlate_CurveConstraint)1236 Handle(GeomPlate_CurveConstraint) GeomPlate_BuildPlateSurface::CurveConstraint( const Standard_Integer order) const
1237 {
1238 return myLinCont->Value(order);
1239 }
Handle(GeomPlate_PointConstraint)1240 Handle(GeomPlate_PointConstraint) GeomPlate_BuildPlateSurface::PointConstraint( const Standard_Integer order) const
1241 {
1242 return myPntCont->Value(order);
1243 }
1244 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
1245 // =========================================================
1246 // P R I V A T E M E T H O D S
1247 // =========================================================
1248 //\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
1249
1250 //=======================================================================
1251 // Function : CourbeJointive
1252 // Purpose : Create a chain of curves to calculate the
1253 // initial surface with the method of max flow.
1254 // Return true if it is a closed contour.
1255 //=======================================================================
1256
1257 Standard_Boolean GeomPlate_BuildPlateSurface::
CourbeJointive(const Standard_Real tolerance)1258 CourbeJointive(const Standard_Real tolerance)
1259 { Standard_Integer nbf = myLinCont->Length();
1260 Standard_Real Ufinal1,Uinit1,Ufinal2,Uinit2;
1261 mySense = new TColStd_HArray1OfInteger(1,nbf,0);
1262 Standard_Boolean result=Standard_True;
1263 Standard_Integer j=1,i;
1264 gp_Pnt P1,P2;
1265
1266 while (j <= (myNbBounds-1))
1267 {
1268 Standard_Integer a=0;
1269 i=j+1;
1270 if (i > myNbBounds)
1271 { result = Standard_False;
1272 a=2;
1273 }
1274 while (a<1)
1275 { if (i > myNbBounds)
1276 { result = Standard_False;
1277 a=2;
1278 }
1279 else
1280 {
1281 Uinit1=myLinCont->Value(j)->FirstParameter();
1282 Ufinal1=myLinCont->Value(j)->LastParameter();
1283 Uinit2=myLinCont->Value(i)->FirstParameter();
1284 Ufinal2=myLinCont->Value(i)->LastParameter();
1285 if (mySense->Value(j)==1)
1286 Ufinal1=Uinit1;
1287 myLinCont->Value(j)->D0(Ufinal1,P1);
1288 myLinCont->Value(i)->D0(Uinit2,P2);
1289 if (P1.Distance(P2)<tolerance)
1290 { if (i!=j+1) {
1291 Handle(GeomPlate_CurveConstraint) tampon = myLinCont->Value(j+1);
1292 myLinCont->SetValue(j+1,myLinCont->Value(i));
1293 myLinCont->SetValue(i,tampon);
1294 Standard_Integer Tmp=myInitOrder->Value(j+1);
1295 // See function TrierTab for the functioning of myInitOrder
1296 myInitOrder->SetValue(j+1,myInitOrder->Value(i));
1297 myInitOrder->SetValue(i,Tmp);
1298
1299
1300 };
1301 a=2;
1302 mySense->SetValue(j+1,0);
1303
1304 }
1305 else
1306 { myLinCont->Value(i)->D0(Ufinal2,P2);
1307
1308 if (P1.Distance(P2)<tolerance)
1309 { if (i!=j+1) {
1310 Handle(GeomPlate_CurveConstraint) tampon =
1311 myLinCont->Value(j+1);
1312 myLinCont->SetValue(j+1,myLinCont->Value(i));
1313 myLinCont->SetValue(i,tampon);
1314 Standard_Integer Tmp=myInitOrder->Value(j+1);
1315 // See function TrierTab for the functioning of myInitOrder
1316 myInitOrder->SetValue(j+1,myInitOrder->Value(i));
1317 myInitOrder->SetValue(i,Tmp);
1318 };
1319 a=2;
1320 mySense->SetValue(j+1,1);
1321 }
1322 }
1323 }
1324 i++;
1325 }
1326 j++;
1327 }
1328 Uinit1=myLinCont->Value( myNbBounds )->FirstParameter();
1329 Ufinal1=myLinCont->Value( myNbBounds )->LastParameter();
1330 Uinit2=myLinCont->Value(1)->FirstParameter();
1331 Ufinal2=myLinCont->Value(1)->LastParameter();
1332 myLinCont->Value( myNbBounds )->D0(Ufinal1,P1);
1333 myLinCont->Value(1)->D0(Uinit2,P2);
1334 if ((mySense->Value( myNbBounds )==0)
1335 &&(P1.Distance(P2)<tolerance))
1336 {
1337 return ((Standard_True)&&(result));
1338 }
1339 myLinCont->Value( myNbBounds )->D0(Uinit1,P1);
1340 if ((mySense->Value( myNbBounds )==1)
1341 &&(P1.Distance(P2)<tolerance))
1342 {
1343 return ((Standard_True)&&(result));
1344 }
1345 else return Standard_False;
1346 }
1347
1348
1349 //-------------------------------------------------------------------------
1350 // Function : ComputeSurfInit
1351 //-------------------------------------------------------------------------
1352 // Calculate the initial surface either by the method of max flow or by
1353 // the method of the plane of inertia if the contour is not closed or if
1354 // there are point constraints.
1355 //-------------------------------------------------------------------------
1356
ComputeSurfInit(const Message_ProgressRange & theProgress)1357 void GeomPlate_BuildPlateSurface::ComputeSurfInit(const Message_ProgressRange& theProgress)
1358 {
1359 Standard_Integer nopt=2, popt=2, Np=1;
1360 Standard_Boolean isHalfSpace = Standard_True;
1361 Standard_Real LinTol = 0.001, AngTol = 0.001; //AngTol = 0.0001; //LinTol = 0.0001
1362
1363 // Option to calculate the initial plane
1364 Standard_Integer NTLinCont = myLinCont->Length(), NTPntCont = myPntCont->Length();
1365 // Table of transformation to preserve the initial order see TrierTab
1366 if (NTLinCont != 0) {
1367 myInitOrder = new TColStd_HArray1OfInteger(1,NTLinCont);
1368 for (Standard_Integer i = 1; i <= NTLinCont; i++)
1369 myInitOrder->SetValue( i, i );
1370 }
1371
1372 Standard_Boolean CourbeJoint = (NTLinCont != 0) && CourbeJointive (myTol3d);
1373 if (CourbeJoint && IsOrderG1())
1374 {
1375 nopt = 3;
1376 // Table contains the cloud of points for calculation of the plane
1377 Standard_Integer NbPoint = 20, Discr = NbPoint/4, pnum = 0;
1378 Handle( TColgp_HArray1OfPnt ) Pts = new TColgp_HArray1OfPnt( 1, (NbPoint+1)*NTLinCont+NTPntCont );
1379 TColgp_SequenceOfVec Vecs, NewVecs;
1380 GeomPlate_SequenceOfAij Aset;
1381 Standard_Real Uinit, Ufinal, Uif;
1382 gp_Vec LastVec;
1383 Standard_Integer i ;
1384 for ( i = 1; i <= NTLinCont; i++)
1385 {
1386 Standard_Integer Order = myLinCont->Value(i)->Order();
1387
1388 NewVecs.Clear();
1389
1390 Uinit = myLinCont->Value(i)->FirstParameter();
1391 Ufinal = myLinCont->Value(i)->LastParameter();
1392 Uif = Ufinal - Uinit;
1393 if (mySense->Value(i) == 1)
1394 {
1395 Uinit = Ufinal;
1396 Uif = -Uif;
1397 }
1398
1399 gp_Vec Vec1, Vec2, Normal;
1400 Standard_Boolean ToReverse = Standard_False;
1401 if (i > 1 && Order >= GeomAbs_G1)
1402 {
1403 gp_Pnt P;
1404 myLinCont->Value(i)->D1( Uinit, P, Vec1, Vec2 );
1405 Normal = Vec1 ^ Vec2;
1406 if (LastVec.IsOpposite( Normal, AngTol ))
1407 ToReverse = Standard_True;
1408 }
1409
1410 for (Standard_Integer j = 0; j <= NbPoint; j++)
1411 { // Number of points per curve = 20
1412 // Linear distribution
1413 Standard_Real Inter = j*Uif/(NbPoint);
1414 if (Order < GeomAbs_G1 || j % Discr != 0)
1415 myLinCont->Value(i)->D0( Uinit+Inter, Pts->ChangeValue(++pnum) );
1416 else
1417 {
1418 myLinCont->Value(i)->D1( Uinit+Inter, Pts->ChangeValue(++pnum), Vec1, Vec2 );
1419 Normal = Vec1 ^ Vec2;
1420 Normal.Normalize();
1421 if (ToReverse)
1422 Normal.Reverse();
1423 Standard_Boolean isNew = Standard_True;
1424 Standard_Integer k ;
1425 for ( k = 1; k <= Vecs.Length(); k++)
1426 if (Vecs(k).IsEqual( Normal, LinTol, AngTol ))
1427 {
1428 isNew = Standard_False;
1429 break;
1430 }
1431 if (isNew)
1432 for (k = 1; k <= NewVecs.Length(); k++)
1433 if (NewVecs(k).IsEqual( Normal, LinTol, AngTol ))
1434 {
1435 isNew = Standard_False;
1436 break;
1437 }
1438 if (isNew)
1439 NewVecs.Append( Normal );
1440 }
1441 }
1442 if (Order >= GeomAbs_G1)
1443 {
1444 isHalfSpace = GeomPlate_BuildAveragePlane::HalfSpace( NewVecs, Vecs, Aset, LinTol, AngTol );
1445 if (! isHalfSpace)
1446 break;
1447 LastVec = Normal;
1448 }
1449 } //for (i = 1; i <= NTLinCont; i++)
1450
1451 if (isHalfSpace)
1452 {
1453 for (i = 1; i <= NTPntCont; i++)
1454 {
1455 Standard_Integer Order = myPntCont->Value(i)->Order();
1456
1457 NewVecs.Clear();
1458 gp_Vec Vec1, Vec2, Normal;
1459 if (Order < GeomAbs_G1)
1460 myPntCont->Value(i)->D0( Pts->ChangeValue(++pnum) );
1461 else
1462 {
1463 myPntCont->Value(i)->D1( Pts->ChangeValue(++pnum), Vec1, Vec2 );
1464 Normal = Vec1 ^ Vec2;
1465 Normal.Normalize();
1466 Standard_Boolean isNew = Standard_True;
1467 for (Standard_Integer k = 1; k <= Vecs.Length(); k++)
1468 if (Vecs(k).IsEqual( Normal, LinTol, AngTol ))
1469 {
1470 isNew = Standard_False;
1471 break;
1472 }
1473 if (isNew)
1474 {
1475 NewVecs.Append( Normal );
1476 isHalfSpace = GeomPlate_BuildAveragePlane::HalfSpace( NewVecs, Vecs, Aset, LinTol, AngTol );
1477 if (! isHalfSpace)
1478 {
1479 NewVecs(1).Reverse();
1480 isHalfSpace = GeomPlate_BuildAveragePlane::HalfSpace( NewVecs, Vecs, Aset, LinTol, AngTol );
1481 }
1482 if (! isHalfSpace)
1483 break;
1484 }
1485 }
1486 } //for (i = 1; i <= NTPntCont; i++)
1487
1488 if (isHalfSpace)
1489 {
1490 Standard_Boolean NullExist = Standard_True;
1491 while (NullExist)
1492 {
1493 NullExist = Standard_False;
1494 for (i = 1; i <= Vecs.Length(); i++)
1495 if (Vecs(i).SquareMagnitude() == 0.)
1496 {
1497 NullExist = Standard_True;
1498 Vecs.Remove(i);
1499 break;
1500 }
1501 }
1502 GeomPlate_BuildAveragePlane BAP( Vecs, Pts );
1503 Standard_Real u1,u2,v1,v2;
1504 BAP.MinMaxBox(u1,u2,v1,v2);
1505 // The space is greater for projections
1506 Standard_Real du = u2 - u1;
1507 Standard_Real dv = v2 - v1;
1508 u1 -= du; u2 += du; v1 -= dv; v2 += dv;
1509 mySurfInit = new Geom_RectangularTrimmedSurface( BAP.Plane(), u1,u2,v1,v2 );
1510 }
1511 } //if (isHalfSpace)
1512 if (!isHalfSpace)
1513 {
1514 #ifdef OCCT_DEBUG
1515 std::cout<<std::endl<<"Normals are not in half space"<<std::endl<<std::endl;
1516 #endif
1517 myIsLinear = Standard_False;
1518 nopt = 2;
1519 }
1520 } //if (NTLinCont != 0 && (CourbeJoint = CourbeJointive( myTol3d )) && IsOrderG1())
1521
1522 if (NTLinCont != 0)
1523 TrierTab( myInitOrder ); // Reorder the table of transformations
1524
1525 if (nopt != 3)
1526 {
1527 if ( NTPntCont != 0)
1528 nopt = 1; //Calculate by the method of plane of inertia
1529 else if (!CourbeJoint || NTLinCont != myNbBounds)
1530 {// throw Standard_Failure("Curves are not joined");
1531 #ifdef OCCT_DEBUG
1532 std::cout << "WARNING : Curves are non-adjacent with tolerance " << myTol3d << std::endl;
1533 #endif
1534 nopt = 1;
1535 }
1536
1537 Standard_Real LenT=0;
1538 Standard_Integer Npt=0;
1539 Standard_Integer NTPoint=20*NTLinCont;
1540 Standard_Integer i ;
1541 for ( i=1;i<=NTLinCont;i++)
1542 LenT+=myLinCont->Value(i)->Length();
1543 for (i=1;i<=NTLinCont;i++)
1544 { Standard_Integer NbPoint= (Standard_Integer )( NTPoint*(myLinCont->Value(i)->Length())/LenT);
1545 if (NbPoint<10)
1546 NbPoint=10;
1547
1548 (void )Npt; // unused but set for debug
1549 Npt += NbPoint;
1550 }
1551 // Table containing a cloud of points for calculation of the plane
1552 Handle(TColgp_HArray1OfPnt) Pts = new TColgp_HArray1OfPnt(1,20*NTLinCont+NTPntCont);
1553 Standard_Integer NbPoint=20;
1554 Standard_Real Uinit,Ufinal, Uif;
1555 for ( i=1;i<=NTLinCont;i++)
1556 { Uinit=myLinCont->Value(i)->FirstParameter();
1557 Ufinal=myLinCont->Value(i)->LastParameter();
1558 Uif=Ufinal-Uinit;
1559 if (mySense->Value(i) == 1)
1560 { Uinit = Ufinal;
1561 Uif = -Uif;
1562 }
1563 for (Standard_Integer j=0; j<NbPoint; j++)
1564 { // Number of points per curve = 20
1565 // Linear distribution
1566 Standard_Real Inter=j*Uif/(NbPoint);
1567 gp_Pnt P;
1568 myLinCont->Value(i)->D0(Uinit+Inter,P);
1569 Pts->SetValue(Np++,P);
1570 }
1571 }
1572 for (i=1;i<=NTPntCont;i++)
1573 { gp_Pnt P;
1574 myPntCont->Value(i)->D0(P);
1575 Pts->SetValue(Np++,P);
1576 }
1577 if (!CourbeJoint)
1578 myNbBounds = 0;
1579 GeomPlate_BuildAveragePlane BAP( Pts, NbPoint*myNbBounds, myTol3d/1000, popt, nopt );
1580 if (!BAP.IsPlane())
1581 {
1582 #ifdef OCCT_DEBUG
1583 std::cout << "WARNING : GeomPlate : the initial surface is not a plane." << std::endl;
1584 #endif
1585
1586 return;
1587 }
1588 Standard_Real u1,u2,v1,v2;
1589 BAP.MinMaxBox(u1,u2,v1,v2);
1590 // The space is greater for projections
1591 Standard_Real du = u2 - u1;
1592 Standard_Real dv = v2 - v1;
1593 u1 -= du; u2 += du; v1 -= dv; v2 += dv;
1594 mySurfInit = new Geom_RectangularTrimmedSurface(BAP.Plane(),u1,u2,v1,v2);
1595 } //if (nopt != 3)
1596
1597 //Comparing metrics of curves and projected curves
1598 if (NTLinCont != 0 && myIsLinear)
1599 {
1600 Handle( Geom_Surface ) InitPlane =
1601 (Handle( Geom_RectangularTrimmedSurface )::DownCast(mySurfInit))->BasisSurface();
1602
1603 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;
1604 Handle( GeomAdaptor_Surface ) hsur =
1605 new GeomAdaptor_Surface( InitPlane );
1606 Standard_Integer NbPoint = 20;
1607 // gp_Pnt P;
1608 // gp_Vec DerC, DerCproj, DU, DV;
1609 // gp_Pnt2d P2d;
1610 // gp_Vec2d DProj;
1611
1612
1613 for (Standard_Integer i = 1; i <= NTLinCont && myIsLinear; i++)
1614 {
1615 Standard_Real FirstPar = myLinCont->Value(i)->FirstParameter();
1616 Standard_Real LastPar = myLinCont->Value(i)->LastParameter();
1617 Standard_Real Uif = (LastPar - FirstPar)/(NbPoint);
1618
1619 Handle( Adaptor3d_Curve ) Curve = myLinCont->Value(i)->Curve3d();
1620 Handle( ProjLib_HCompProjectedCurve ) ProjCurve = new ProjLib_HCompProjectedCurve (hsur, Curve, myTol3d, myTol3d);
1621 Adaptor3d_CurveOnSurface AProj(ProjCurve, hsur);
1622
1623 gp_Pnt P;
1624 gp_Vec DerC, DerCproj;
1625 for (Standard_Integer j = 1 ; j < NbPoint && myIsLinear; j++)
1626 {
1627 Standard_Real Inter = FirstPar + j*Uif;
1628 Curve->D1(Inter, P, DerC);
1629 AProj.D1(Inter, P, DerCproj);
1630
1631 Standard_Real A1 = DerC.Magnitude();
1632 Standard_Real A2 = DerCproj.Magnitude();
1633 if (A2 <= 1.e-20) {
1634 Ratio = 1.e20;
1635 }
1636 else {
1637 Ratio = A1 / A2;
1638 }
1639 if (Ratio > R1 || Ratio < R2)
1640 {
1641 myIsLinear = Standard_False;
1642 break;
1643 }
1644 }
1645 }
1646 #ifdef OCCT_DEBUG
1647 if (! myIsLinear)
1648 std::cout <<"Metrics are too different :"<< Ratio<<std::endl;
1649 #endif
1650 // myIsLinear = Standard_True; // !!
1651 } //comparing metrics of curves and projected curves
1652
1653
1654 if (! myIsLinear)
1655 {
1656 myPlanarSurfInit = mySurfInit;
1657 #ifdef DRAW
1658 if (Affich) {
1659 char name[256];
1660 sprintf(name,"planinit_%d",NbPlan+1);
1661 DrawTrSurf::Set(name, mySurfInit);
1662 }
1663 #endif
1664 Standard_Real u1,v1,u2,v2;
1665 mySurfInit->Bounds(u1,v1,u2,v2);
1666 GeomAdaptor_Surface Surf(mySurfInit);
1667 myTolU = Surf.UResolution(myTol3d);
1668 myTolV = Surf.VResolution(myTol3d);
1669 myProj.Initialize(Surf,u1,v1,u2,v2,
1670 myTolU,myTolV);
1671
1672 //======================================================================
1673 // Projection of curves
1674 //======================================================================
1675 Standard_Integer i;
1676 for (i = 1; i <= NTLinCont; i++)
1677 if (myLinCont->Value(i)->Curve2dOnSurf().IsNull())
1678 myLinCont->ChangeValue(i)->SetCurve2dOnSurf(ProjectCurve(myLinCont->Value(i)->Curve3d()));
1679
1680 //======================================================================
1681 // Projection of points
1682 //======================================================================
1683 for (i = 1; i<=NTPntCont; i++)
1684 { gp_Pnt P;
1685 myPntCont->Value(i)->D0(P);
1686 if (!myPntCont->Value(i)->HasPnt2dOnSurf())
1687 myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P));
1688 }
1689
1690 //======================================================================
1691 // Number of points by curve
1692 //======================================================================
1693 if ((NTLinCont !=0)&&(myNbPtsOnCur!=0))
1694 CalculNbPtsInit();
1695
1696 //======================================================================
1697 // Management of incompatibilities between curves
1698 //======================================================================
1699 Handle(GeomPlate_HArray1OfSequenceOfReal) PntInter;
1700 Handle(GeomPlate_HArray1OfSequenceOfReal) PntG1G1;
1701 if (NTLinCont != 0)
1702 {
1703 PntInter = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
1704 PntG1G1 = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
1705 Intersect(PntInter, PntG1G1);
1706 }
1707
1708 //====================================================================
1709 // Discretization of curves
1710 //====================================================================
1711 Discretise(PntInter, PntG1G1);
1712 //====================================================================
1713 //Preparation of points of constraint for plate
1714 //====================================================================
1715 LoadCurve( 0, 0 );
1716 if (myPntCont->Length() != 0)
1717 LoadPoint( 0, 0 );
1718 //====================================================================
1719 // Construction of the surface
1720 //====================================================================
1721 Message_ProgressScope aPS(theProgress, "ComputeSurfInit", 1);
1722 myPlate.SolveTI(2, ComputeAnisotropie(), aPS.Next());
1723 if (theProgress.UserBreak())
1724 {
1725 return;
1726 }
1727
1728 if (!myPlate.IsDone())
1729 {
1730 #ifdef OCCT_DEBUG
1731 std::cout << "WARNING : GeomPlate : calculation of Plate failed" << std::endl;
1732 #endif
1733 return;
1734 }
1735
1736 myGeomPlateSurface = new GeomPlate_Surface( mySurfInit, myPlate );
1737
1738 GeomPlate_MakeApprox App(myGeomPlateSurface, myTol3d,
1739 1, 3,
1740 15*myTol3d,
1741 -1, GeomAbs_C0);
1742 mySurfInit = App.Surface();
1743
1744 mySurfInitIsGive = Standard_True;
1745 myPlate.Init(); // Reset
1746
1747 for (i=1;i<=NTLinCont;i++)
1748 {
1749 Handle( Geom2d_Curve ) NullCurve;
1750 NullCurve.Nullify();
1751 myLinCont->ChangeValue(i)->SetCurve2dOnSurf( NullCurve);
1752 }
1753 }
1754
1755 #ifdef DRAW
1756 if (Affich) {
1757 char name[256];
1758 sprintf(name,"surfinit_%d",++NbPlan);
1759 DrawTrSurf::Set(name, mySurfInit);
1760 }
1761 #endif
1762 }
1763
1764 //---------------------------------------------------------
1765 // Function : Intersect
1766 //---------------------------------------------------------
1767 // Find intersections between 2d curves
1768 // If the intersection is compatible (in cases G1-G1)
1769 // remove the point on one of two curves
1770 //---------------------------------------------------------
1771 void GeomPlate_BuildPlateSurface::
Intersect(Handle (GeomPlate_HArray1OfSequenceOfReal)& PntInter,Handle (GeomPlate_HArray1OfSequenceOfReal)& PntG1G1)1772 Intersect(Handle(GeomPlate_HArray1OfSequenceOfReal)& PntInter,
1773 Handle(GeomPlate_HArray1OfSequenceOfReal)& PntG1G1)
1774 {
1775 Standard_Integer NTLinCont = myLinCont->Length();
1776 Geom2dInt_GInter Intersection;
1777 Geom2dAdaptor_Curve Ci, Cj;
1778 IntRes2d_IntersectionPoint int2d;
1779 gp_Pnt P1,P2;
1780 gp_Pnt2d P2d;
1781 gp_Vec2d V2d;
1782 // if (!mySurfInitIsGive)
1783 for (Standard_Integer i=1;i<=NTLinCont;i++)
1784 {
1785 //Standard_Real NbPnt_i=myLinCont->Value(i)->NbPoints();
1786 // Find the intersection with each curve including the curve itself
1787 Ci.Load(myLinCont->Value(i)->Curve2dOnSurf());
1788 for(Standard_Integer j=i; j<=NTLinCont; j++)
1789 {
1790 Cj.Load(myLinCont->Value(j)->Curve2dOnSurf());
1791 if (i==j)
1792 Intersection.Perform(Ci, myTol2d*10, myTol2d*10);
1793 else
1794 Intersection.Perform(Ci, Cj, myTol2d*10, myTol2d*10);
1795
1796 if (!Intersection.IsEmpty())
1797 { // there is one intersection
1798 Standard_Integer nbpt = Intersection.NbPoints();
1799 // number of points of intersection
1800 for (Standard_Integer k = 1; k <= nbpt; k++)
1801 { int2d = Intersection.Point(k);
1802 myLinCont->Value(i)->D0(int2d.ParamOnFirst(),P1);
1803 myLinCont->Value(j)->D0(int2d.ParamOnSecond(),P2);
1804 #ifdef OCCT_DEBUG
1805 if (Affich> 1)
1806 {
1807 std::cout << " Intersection "<< k << " entre " << i
1808 << " &" << j << std::endl;
1809 std::cout <<" Distance = "<< P1.Distance(P2) << std::endl;
1810 }
1811 #endif
1812 if (P1.Distance( P2 ) < myTol3d)
1813 { // 2D intersection corresponds to close 3D points.
1814 // Note the interval, in which the point needs to be removed
1815 // to avoid duplications, which cause
1816 // error in plate. The point on curve i is removed;
1817 // the point on curve j is preserved;
1818 // the length of interval is a length 2d
1819 // corresponding in 3d to myTol3d
1820 Standard_Real tolint = Ci.Resolution(myTol3d);
1821 Ci.D1(int2d.ParamOnFirst(),P2d, V2d);
1822 Standard_Real aux = V2d.Magnitude();
1823 if (aux > 1.e-7)
1824 {
1825 aux = myTol3d/aux;
1826 if (aux > 100*tolint) tolint*=100;
1827 else tolint = aux;
1828 }
1829 else tolint*=100;
1830
1831 PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() - tolint);
1832 PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() + tolint);
1833 // If G1-G1
1834 if ( (myLinCont->Value(i)->Order()==1)
1835 &&(myLinCont->Value(j)->Order()==1))
1836 { gp_Vec v11,v12,v13,v14,v15,v16,v21,v22,v23,v24,v25,v26;
1837 myLinCont->Value(i)->D2( int2d.ParamOnFirst(), P1, v11, v12, v13, v14, v15 );
1838 myLinCont->Value(j)->D2( int2d.ParamOnSecond(), P2, v21, v22, v23, v24, v25 );
1839 v16=v11^v12;v26=v21^v22;
1840 Standard_Real ant=v16.Angle(v26);
1841 if (ant>(M_PI/2))
1842 ant= M_PI -ant;
1843 if ((Abs(v16*v15-v16*v25)>(myTol3d/1000))
1844 ||(Abs(ant)>myTol3d/1000))
1845 // Non-compatible ==> remove zone in constraint G1
1846 // corresponding to 3D tolerance of 0.01
1847 { Standard_Real coin;
1848 Standard_Real Tol= 100 * myTol3d;
1849 Standard_Real A1;
1850 gp_Pnt2d P1temp,P2temp;
1851 gp_Vec2d V1,V2;
1852 myLinCont->Value(i)->Curve2dOnSurf()->D1( int2d.ParamOnFirst(), P1temp, V1);
1853 myLinCont->Value(j)->Curve2dOnSurf()->D1( int2d.ParamOnSecond(), P2temp, V2);
1854 A1 = V1.Angle(V2);
1855 if (A1>(M_PI/2))
1856 A1= M_PI - A1;
1857 if (Abs(Abs(A1)-M_PI)<myTolAng) Tol = 100000 * myTol3d;
1858 #ifdef OCCT_DEBUG
1859 if (Affich) std::cout <<"Angle between curves "<<i<<","<<j
1860 <<" "<<Abs(Abs(A1)-M_PI)<<std::endl;
1861 #endif
1862
1863 coin = Ci.Resolution(Tol);
1864 Standard_Real Par1=int2d.ParamOnFirst()-coin,
1865 Par2=int2d.ParamOnFirst()+coin;
1866 // Storage of the interval for curve i
1867 PntG1G1->ChangeValue(i).Append(Par1);
1868 PntG1G1->ChangeValue(i).Append(Par2);
1869 coin = Cj.Resolution(Tol);
1870 Par1=int2d.ParamOnSecond()-coin;
1871 Par2=int2d.ParamOnSecond()+coin;
1872 // Storage of the interval for curve j
1873 PntG1G1->ChangeValue(j).Append(Par1);
1874 PntG1G1->ChangeValue(j).Append(Par2);
1875 }
1876 }
1877 // If G0-G1
1878 if ((myLinCont->Value(i)->Order()==0 && myLinCont->Value(j)->Order()==1) ||
1879 (myLinCont->Value(i)->Order()==1 && myLinCont->Value(j)->Order()==0))
1880 {
1881 gp_Vec vec, vecU, vecV, N;
1882 if (myLinCont->Value(i)->Order() == 0)
1883 {
1884 Handle( Adaptor3d_Curve ) theCurve = myLinCont->Value(i)->Curve3d();
1885 theCurve->D1( int2d.ParamOnFirst(), P1, vec );
1886 myLinCont->Value(j)->D1( int2d.ParamOnSecond(), P2, vecU, vecV );
1887 }
1888 else
1889 {
1890 Handle( Adaptor3d_Curve ) theCurve = myLinCont->Value(j)->Curve3d();
1891 theCurve->D1( int2d.ParamOnSecond(), P2, vec );
1892 myLinCont->Value(i)->D1( int2d.ParamOnFirst(), P1, vecU, vecV );
1893 }
1894 N = vecU ^ vecV;
1895 Standard_Real Angle = vec.Angle( N );
1896 Angle = Abs( M_PI/2-Angle );
1897 if (Angle > myTolAng/10.) //????????? //if (Abs( scal ) > myTol3d/100)
1898 {
1899 // Non-compatible ==> one removes zone in constraint G0 and G1
1900 // corresponding to 3D tolerance of 0.01
1901 Standard_Real coin;
1902 Standard_Real Tol= 100 * myTol3d;
1903 Standard_Real A1;
1904 gp_Pnt2d P1temp,P2temp;
1905 gp_Vec2d V1,V2;
1906 myLinCont->Value(i)->Curve2dOnSurf()->D1( int2d.ParamOnFirst(), P1temp, V1);
1907 myLinCont->Value(j)->Curve2dOnSurf()->D1( int2d.ParamOnSecond(), P2temp, V2);
1908 A1 = V1.Angle( V2 );
1909 if (A1 > M_PI/2)
1910 A1= M_PI - A1;
1911 if (Abs(Abs(A1) - M_PI) < myTolAng) Tol = 100000 * myTol3d;
1912 #ifdef OCCT_DEBUG
1913 if (Affich) std::cout <<"Angle entre Courbe "<<i<<","<<j
1914 <<" "<<Abs(Abs(A1)-M_PI)<<std::endl;
1915 #endif
1916 if (myLinCont->Value(i)->Order() == 1)
1917 {
1918 coin = Ci.Resolution(Tol);
1919 coin *= Angle / myTolAng * 10.;
1920 #ifdef OCCT_DEBUG
1921 std::cout<<std::endl<<"coin = "<<coin<<std::endl;
1922 #endif
1923 Standard_Real Par1 = int2d.ParamOnFirst() - coin;
1924 Standard_Real Par2 = int2d.ParamOnFirst() + coin;
1925 // Storage of the interval for curve i
1926 PntG1G1->ChangeValue(i).Append( Par1 );
1927 PntG1G1->ChangeValue(i).Append( Par2 );
1928 }
1929 else
1930 {
1931 coin = Cj.Resolution(Tol);
1932 coin *= Angle / myTolAng * 10.;
1933 #ifdef OCCT_DEBUG
1934 std::cout<<std::endl<<"coin = "<<coin<<std::endl;
1935 #endif
1936 Standard_Real Par1 = int2d.ParamOnSecond() - coin;
1937 Standard_Real Par2 = int2d.ParamOnSecond() + coin;
1938 // Storage of the interval for curve j
1939 PntG1G1->ChangeValue(j).Append( Par1 );
1940 PntG1G1->ChangeValue(j).Append( Par2 );
1941 }
1942 }
1943 }
1944 } //if (P1.Distance( P2 ) < myTol3d)
1945 else {
1946 // 2D intersection corresponds to extended 3D points.
1947 // Note the interval where it is necessary to remove
1948 // the points to avoid duplications causing
1949 // error in plate. The point on curve i is removed,
1950 // the point on curve j is preserved.
1951 // The length of interval is 2D length
1952 // corresponding to the distance of points in 3D to myTol3d
1953 Standard_Real tolint, Dist;
1954 Dist = P1.Distance(P2);
1955 tolint = Ci.Resolution(Dist);
1956 PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() - tolint);
1957 PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() + tolint);
1958 if (j!=i)
1959 {
1960 tolint = Cj.Resolution(Dist);
1961 PntInter->ChangeValue(j).
1962 Append( int2d.ParamOnSecond() - tolint);
1963 PntInter->ChangeValue(j).
1964 Append( int2d.ParamOnSecond() + tolint);
1965 }
1966
1967 #ifdef OCCT_DEBUG
1968 std::cout << "Attention: Two points 3d have the same projection dist = " << Dist << std::endl;
1969 #endif
1970 #ifdef DRAW
1971 if (Affich > 1)
1972 {
1973 Handle(Draw_Marker3D) mark = new (Draw_Marker3D)(P1, Draw_X, Draw_vert);
1974 char name[256];
1975 sprintf(name,"mark_%d",++NbMark);
1976 Draw::Set(name, mark);
1977 }
1978 #endif
1979 }
1980 }
1981 }
1982 }
1983 }
1984 }
1985
1986 //---------------------------------------------------------
1987 // Function : Discretize
1988 //---------------------------------------------------------
1989 // Discretize curves according to parameters
1990 // the table of sequences Parcont contains all
1991 // parameter of points on curves
1992 // Field myPlateCont contains parameter of points on a plate;
1993 // it excludes duplicate points and imcompatible zones.
1994 // The first part corresponds to verification of compatibility
1995 // and to removal of duplicate points.
1996 //---------------------------------------------------------
1997 void GeomPlate_BuildPlateSurface::
Discretise(const Handle (GeomPlate_HArray1OfSequenceOfReal)& PntInter,const Handle (GeomPlate_HArray1OfSequenceOfReal)& PntG1G1)1998 Discretise(const Handle(GeomPlate_HArray1OfSequenceOfReal)& PntInter,
1999 const Handle(GeomPlate_HArray1OfSequenceOfReal)& PntG1G1)
2000 {
2001 Standard_Integer NTLinCont = myLinCont->Length();
2002 Standard_Boolean ACR;
2003 Handle(Geom2d_Curve) C2d;
2004 Geom2dAdaptor_Curve AC2d;
2005 // Handle(Adaptor_HCurve2d) HC2d;
2006 Handle(Law_Interpol) acrlaw = new (Law_Interpol) ();
2007 myPlateCont = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
2008 myParCont = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
2009
2010
2011 //===========================================================================
2012 // Construction of the table containing parameters of constraint points
2013 //===========================================================================
2014 Standard_Real Uinit, Ufinal, Length2d=0, Inter;
2015 Standard_Real CurLength;
2016 Standard_Integer NbPnt_i, NbPtInter, NbPtG1G1;
2017 Handle(GeomPlate_CurveConstraint) LinCont;
2018
2019 for (Standard_Integer i=1;i<=NTLinCont;i++) {
2020 LinCont = myLinCont->Value(i);
2021 Uinit=LinCont->FirstParameter();
2022 Ufinal=LinCont->LastParameter();
2023 // HC2d=LinCont->ProjectedCurve();
2024 // if(HC2d.IsNull())
2025 // ACR = (!HC2d.IsNull() || !C2d.IsNull());
2026 C2d= LinCont->Curve2dOnSurf();
2027 ACR = (!C2d.IsNull());
2028 if (ACR) {
2029 // Construct a law close to curvilinear abscissa
2030 if(!C2d.IsNull()) AC2d.Load(C2d);
2031 // AC2d.Load(LinCont->Curve2dOnSurf());
2032 Standard_Integer ii, Nbint = 20;
2033 Standard_Real U;
2034 TColgp_Array1OfPnt2d tabP2d(1, Nbint+1);
2035 tabP2d(1).SetY(Uinit);
2036 tabP2d(1).SetX(0.);
2037 tabP2d(Nbint+1).SetY(Ufinal);
2038 /* if (!HC2d.IsNull())
2039
2040 Length2d = GCPnts_AbscissaPoint::Length(HC2d->Curve2d(), Uinit, Ufinal);
2041 else*/
2042 Length2d = GCPnts_AbscissaPoint::Length(AC2d, Uinit, Ufinal);
2043
2044 tabP2d(Nbint+1).SetX(Length2d);
2045 for (ii = 2; ii<= Nbint; ii++) {
2046 U = Uinit + (Ufinal-Uinit)*((1-Cos((ii-1)*M_PI/(Nbint)))/2);
2047 tabP2d(ii).SetY(U);
2048 /* if (!HC2d.IsNull()) {
2049 Standard_Real L = GCPnts_AbscissaPoint::Length(HC2d->Curve2d(), Uinit, U);
2050 tabP2d(ii).SetX(L);
2051 }
2052 else*/
2053 tabP2d(ii).SetX(GCPnts_AbscissaPoint::Length(AC2d, Uinit, U));
2054 }
2055 acrlaw->Set(tabP2d);
2056 }
2057
2058 NbPnt_i = (Standard_Integer )( LinCont->NbPoints());
2059 NbPtInter= PntInter->Value(i).Length();
2060 NbPtG1G1= PntG1G1->Value(i).Length();
2061
2062 #ifdef OCCT_DEBUG
2063 if (Affich > 1) {
2064 std::cout << "Courbe : " << i << std::endl;
2065 std::cout << " NbPnt, NbPtInter, NbPtG1G1 :" << NbPnt_i << ", "
2066 << NbPtInter << ", " << NbPtG1G1 << std::endl;
2067 }
2068 #endif
2069 for (Standard_Integer j=1; j<=NbPnt_i; j++)
2070 {
2071 // Distribution of points in cosine following ACR 2D
2072 // To avoid points of accumulation in 2D
2073 //Inter=Uinit+(Uif)*((-cos(M_PI*((j-1)/(NbPnt_i-1)))+1)/2);
2074 if (j==NbPnt_i)
2075 Inter=Ufinal;// to avoid bug on Sun
2076 else if (ACR) {
2077 CurLength = Length2d*(1-Cos((j-1)*M_PI/(NbPnt_i-1)))/2;
2078 Inter = acrlaw->Value(CurLength);
2079 }
2080 else {
2081 Inter=Uinit+(Ufinal-Uinit)*((1-Cos((j-1)*M_PI/(NbPnt_i-1)))/2);
2082 }
2083 myParCont->ChangeValue(i).Append(Inter);// add a point
2084 if (NbPtInter!=0)
2085 {
2086 for(Standard_Integer l=1;l<=NbPtInter;l+=2)
2087 {
2088 // check if the point Inter is in the interval
2089 // PntInter[i] PntInter[i+1]
2090 // in which case it is not necessary to store it (problem with duplicates)
2091 if ((Inter>PntInter->Value(i).Value(l))
2092 &&(Inter<PntInter->Value(i).Value(l+1)))
2093 {
2094 l=NbPtInter+2;
2095 // leave the loop without storing the point
2096 }
2097 else
2098 {
2099 if (l+1>=NbPtInter) {
2100 // one has parsed the entire table : the point
2101 // does not belong to a common point interval
2102 if (NbPtG1G1!=0)
2103 {
2104 // if there exists an incompatible interval
2105 for(Standard_Integer k=1;k<=NbPtG1G1;k+=2)
2106 {
2107 if ((Inter>PntG1G1->Value(i).Value(k))
2108 &&(Inter<PntG1G1->Value(i).Value(k+1)))
2109 {
2110 k=NbPtG1G1+2; // to leave the loop
2111 // Add points of constraint G0
2112 gp_Pnt P3d,PP,Pdif;
2113 gp_Pnt2d P2d;
2114
2115 AC2d.D0(Inter, P2d);
2116 LinCont->D0(Inter,P3d);
2117 mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP);
2118 Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1),
2119 -PP.Coord(2)+P3d.Coord(2),
2120 -PP.Coord(3)+P3d.Coord(3));
2121 Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0);
2122 myPlate.Load(PC);
2123 }
2124 else // the point does not belong to interval G1
2125 {
2126 if (k+1>=NbPtG1G1)
2127 {
2128 myPlateCont->ChangeValue(i).Append(Inter);
2129 // add the point
2130 }
2131 }
2132 }
2133 }
2134 else
2135 {
2136 myPlateCont->ChangeValue(i).Append(Inter);
2137 // add the point
2138 }
2139 }
2140 }
2141 }
2142 }
2143 else
2144 {
2145 if (NbPtG1G1!=0) // there exist an incompatible interval
2146 {
2147 for(Standard_Integer k=1;k<=NbPtG1G1;k+=2)
2148 {
2149 if ((Inter>PntG1G1->Value(i).Value(k))
2150 &&(Inter<PntG1G1->Value(i).Value(k+1)))
2151 {
2152 k=NbPtG1G1+2; // to leave the loop
2153 // Add points of constraint G0
2154 gp_Pnt P3d,PP,Pdif;
2155 gp_Pnt2d P2d;
2156
2157 AC2d.D0(Inter, P2d);
2158 LinCont->D0(Inter,P3d);
2159 mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP);
2160 Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1),
2161 -PP.Coord(2)+P3d.Coord(2),
2162 -PP.Coord(3)+P3d.Coord(3));
2163 Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0);
2164 myPlate.Load(PC);
2165
2166 }
2167 else // the point does not belong to interval G1
2168 {
2169 if (k+1>=NbPtG1G1)
2170 {
2171 myPlateCont->ChangeValue(i).Append(Inter);
2172 // add the point
2173 }
2174 }
2175 }
2176 }
2177 else
2178 {
2179 if ( ( (!mySurfInitIsGive)
2180 &&(Geom2dAdaptor_Curve(LinCont->Curve2dOnSurf()).GetType()!=GeomAbs_Circle))
2181 || ( (j>1) &&(j<NbPnt_i))) // exclude extremeties
2182 myPlateCont->ChangeValue(i).Append(Inter);// add the point
2183 }
2184 }
2185 }
2186 }
2187 }
2188 //---------------------------------------------------------
2189 // Function : CalculNbPtsInit
2190 //---------------------------------------------------------
2191 // Calculate the number of points by curve depending on the
2192 // length for the first iteration
2193 //---------------------------------------------------------
CalculNbPtsInit()2194 void GeomPlate_BuildPlateSurface::CalculNbPtsInit ()
2195 {
2196 Standard_Real LenT=0;
2197 Standard_Integer NTLinCont=myLinCont->Length();
2198 Standard_Integer NTPoint=(Standard_Integer )( myNbPtsOnCur*NTLinCont);
2199 Standard_Integer i ;
2200
2201 for ( i=1;i<=NTLinCont;i++)
2202 LenT+=myLinCont->Value(i)->Length();
2203 for ( i=1;i<=NTLinCont;i++)
2204 { Standard_Integer Cont=myLinCont->Value(i)->Order();
2205 switch(Cont)
2206 { case 0 : // Case G0 *1.2
2207 myLinCont->ChangeValue(i)->SetNbPoints(
2208 Standard_Integer(1.2*NTPoint*(myLinCont->Value(i)->Length())/LenT));
2209 break;
2210 case 1 : // Case G1 *1
2211 myLinCont->ChangeValue(i)->SetNbPoints(
2212 Standard_Integer(NTPoint*(myLinCont->Value(i)->Length())/LenT));
2213 break;
2214 case 2 : // Case G2 *0.7
2215 myLinCont->ChangeValue(i)->SetNbPoints(
2216 Standard_Integer(0.7*NTPoint*(myLinCont->Value(i)->Length())/LenT));
2217 break;
2218 }
2219 if (myLinCont->Value(i)->NbPoints()<3)
2220 myLinCont->ChangeValue(i)->SetNbPoints(3);
2221 }
2222 }
2223 //---------------------------------------------------------
2224 // Function : LoadCurve
2225 //---------------------------------------------------------
2226 // Starting from table myParCont load all the points noted in plate
2227 //---------------------------------------------------------
LoadCurve(const Standard_Integer NbBoucle,const Standard_Integer OrderMax)2228 void GeomPlate_BuildPlateSurface::LoadCurve(const Standard_Integer NbBoucle,
2229 const Standard_Integer OrderMax )
2230 {
2231 gp_Pnt P3d,Pdif,PP;
2232 gp_Pnt2d P2d;
2233 Standard_Integer NTLinCont=myLinCont->Length(), i, j;
2234 Standard_Integer Tang, Nt;
2235
2236
2237 for (i=1; i<=NTLinCont; i++){
2238 Handle(GeomPlate_CurveConstraint) CC = myLinCont->Value(i);
2239 if (CC ->Order()!=-1) {
2240 Tang = Min(CC->Order(), OrderMax);
2241 Nt = myPlateCont->Value(i).Length();
2242 if (Tang!=-1)
2243 for (j=1; j<=Nt; j++) {// Loading of points G0 on boundaries
2244 CC ->D0(myPlateCont->Value(i).Value(j),P3d);
2245 if (!CC->ProjectedCurve().IsNull())
2246 P2d = CC->ProjectedCurve()->Value(myPlateCont->Value(i).Value(j));
2247
2248 else {
2249 if (!CC->Curve2dOnSurf().IsNull())
2250 P2d = CC->Curve2dOnSurf()->Value(myPlateCont->Value(i).Value(j));
2251 else
2252 P2d = ProjectPoint(P3d);
2253 }
2254 mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP);
2255 Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1),
2256 -PP.Coord(2)+P3d.Coord(2),
2257 -PP.Coord(3)+P3d.Coord(3));
2258 Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0);
2259 myPlate.Load(PC);
2260
2261 // Loading of points G1
2262 if (Tang==1) { // ==1
2263 gp_Vec V1,V2,V3,V4;
2264 CC->D1( myPlateCont->Value(i).Value(j), PP, V1, V2 );
2265 mySurfInit->D1( P2d.Coord(1), P2d.Coord(2), PP, V3, V4 );
2266
2267 Plate_D1 D1final(V1.XYZ(),V2.XYZ());
2268 Plate_D1 D1init(V3.XYZ(),V4.XYZ());
2269 if (! myFree)
2270 {
2271 Plate_GtoCConstraint GCC( P2d.XY(), D1init, D1final );
2272 myPlate.Load(GCC);
2273 }
2274 else if (NbBoucle == 1)
2275 {
2276 Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), D1init, D1final );
2277 myPlate.Load( FreeGCC );
2278 }
2279 else {
2280 gp_Vec DU, DV, Normal, DerPlateU, DerPlateV;
2281
2282 Normal = V1^V2;
2283 //Normal.Normalize();
2284 Standard_Real norm = Normal.Magnitude();
2285 if (norm > 1.e-12) Normal /= norm;
2286 DerPlateU = myPrevPlate.EvaluateDerivative( P2d.XY(), 1, 0 );
2287 DerPlateV = myPrevPlate.EvaluateDerivative( P2d.XY(), 0, 1 );
2288
2289 DU.SetLinearForm(-(V3 + DerPlateU).Dot(Normal), Normal, DerPlateU);
2290 DV.SetLinearForm(-(V4 + DerPlateV).Dot(Normal), Normal, DerPlateV);
2291 Plate_PinpointConstraint PinU( P2d.XY(), DU.XYZ(), 1, 0 );
2292 Plate_PinpointConstraint PinV( P2d.XY(), DV.XYZ(), 0, 1 );
2293 myPlate.Load( PinU );
2294 myPlate.Load( PinV );
2295 }
2296 }
2297 // Loading of points G2
2298 if (Tang==2) // ==2
2299 { gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10;
2300 CC->D2(myPlateCont->Value(i).Value(j),
2301 PP,V1,V2,V5,V6,V7);
2302 mySurfInit->D2(P2d.Coord(1),P2d.Coord(2),PP,V3,V4,V8,V9,V10);
2303
2304 Plate_D1 D1final(V1.XYZ(),V2.XYZ());
2305 Plate_D1 D1init(V3.XYZ(),V4.XYZ());
2306 Plate_D2 D2final (V5.XYZ(),V6.XYZ(),V7.XYZ());
2307 Plate_D2 D2init (V8.XYZ(),V9.XYZ(),V10.XYZ());
2308 // if (! myFree)
2309 // {
2310 Plate_GtoCConstraint GCC(P2d.XY(),D1init,D1final,D2init,D2final);
2311 myPlate.Load(GCC);
2312 // }
2313 // else // Good but too expansive
2314 // {
2315 // Plate_FreeGtoCConstraint FreeGCC( P2d.XY(),
2316 // D1init, D1final, D2init, D2final );
2317 // myPlate.Load( FreeGCC );
2318 // }
2319
2320 }
2321 }
2322 }
2323 }
2324 }
2325
2326
2327 //---------------------------------------------------------
2328 // Function : LoadPoint
2329 //---------------------------------------------------------
2330 //void GeomPlate_BuildPlateSurface::LoadPoint(const Standard_Integer NbBoucle,
LoadPoint(const Standard_Integer,const Standard_Integer OrderMax)2331 void GeomPlate_BuildPlateSurface::LoadPoint(const Standard_Integer ,
2332 const Standard_Integer OrderMax)
2333 {
2334 gp_Pnt P3d,Pdif,PP;
2335 gp_Pnt2d P2d;
2336 Standard_Integer NTPntCont=myPntCont->Length();
2337 Standard_Integer Tang, i;
2338 // gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10;
2339
2340 // Loading of points of point constraints
2341 for (i=1;i<=NTPntCont;i++) {
2342 myPntCont->Value(i)->D0(P3d);
2343 P2d = myPntCont->Value(i)->Pnt2dOnSurf();
2344 mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP);
2345 Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1),
2346 -PP.Coord(2)+P3d.Coord(2),
2347 -PP.Coord(3)+P3d.Coord(3));
2348 Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0);
2349 myPlate.Load(PC);
2350 Tang = Min(myPntCont->Value(i)->Order(), OrderMax);
2351 if (Tang==1) {// ==1
2352 gp_Vec V1,V2,V3,V4;
2353 myPntCont->Value(i)->D1(PP,V1,V2);
2354 mySurfInit->D1(P2d.Coord(1),P2d.Coord(2),PP,V3,V4);
2355 Plate_D1 D1final(V1.XYZ(),V2.XYZ());
2356 Plate_D1 D1init(V3.XYZ(),V4.XYZ());
2357 if (! myFree)
2358 {
2359 Plate_GtoCConstraint GCC( P2d.XY(), D1init, D1final );
2360 myPlate.Load( GCC );
2361 }
2362 else
2363 {
2364 Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), D1init, D1final );
2365 myPlate.Load( FreeGCC );
2366 }
2367 }
2368 // Loading of points G2 GeomPlate_PlateG0Criterion
2369 if (Tang==2) // ==2
2370 { gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10;
2371 myPntCont->Value(i)->D2(PP,V1,V2,V5,V6,V7);
2372 // gp_Vec Tv2 = V1^V2;
2373 mySurfInit->D2(P2d.Coord(1),P2d.Coord(2),PP,V3,V4,V8,V9,V10);
2374 Plate_D1 D1final(V1.XYZ(),V2.XYZ());
2375 Plate_D1 D1init(V3.XYZ(),V4.XYZ());
2376 Plate_D2 D2final (V5.XYZ(),V6.XYZ(),V7.XYZ());
2377 Plate_D2 D2init (V8.XYZ(),V9.XYZ(),V10.XYZ());
2378 // if (! myFree)
2379 // {
2380 Plate_GtoCConstraint GCC( P2d.XY(), D1init, D1final, D2init, D2final );
2381 myPlate.Load( GCC );
2382 // }
2383 // else // Good but too expansive
2384 // {
2385 // Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), D1init, D1final, D2init//, D2final );
2386 // myPlate.Load( FreeGCC );
2387 // }
2388 }
2389 }
2390
2391 }
2392 //---------------------------------------------------------
2393 // Function : VerifSurface
2394 //---------------------------------------------------------
2395 Standard_Boolean GeomPlate_BuildPlateSurface::
VerifSurface(const Standard_Integer NbBoucle)2396 VerifSurface(const Standard_Integer NbBoucle)
2397 {
2398 //======================================================================
2399 // Calculate errors
2400 //======================================================================
2401 Standard_Integer NTLinCont=myLinCont->Length();
2402 Standard_Boolean Result=Standard_True;
2403
2404 // variable for error calculation
2405 myG0Error=0,myG1Error =0, myG2Error=0;
2406
2407 for (Standard_Integer i=1;i<=NTLinCont;i++) {
2408 Handle(GeomPlate_CurveConstraint) LinCont;
2409 LinCont = myLinCont->Value(i);
2410 if (LinCont->Order()!=-1) {
2411 Standard_Integer NbPts_i = myParCont->Value(i).Length();
2412 if (NbPts_i<3)
2413 NbPts_i=4;
2414 Handle(TColStd_HArray1OfReal) tdist =
2415 new TColStd_HArray1OfReal(1,NbPts_i-1);
2416 Handle(TColStd_HArray1OfReal) tang =
2417 new TColStd_HArray1OfReal(1,NbPts_i-1);
2418 Handle(TColStd_HArray1OfReal) tcourb =
2419 new TColStd_HArray1OfReal(1,NbPts_i-1);
2420
2421 EcartContraintesMil (i,tdist,tang,tcourb);
2422
2423 Standard_Real diffDistMax=0, diffAngMax=0;
2424 //Standard_Real SdiffDist=0, SdiffAng=0;
2425 Standard_Integer NdiffDist=0,NdiffAng=0;
2426
2427
2428 for (Standard_Integer j=1;j<NbPts_i;j++)
2429 { if (tdist->Value(j)>myG0Error)
2430 myG0Error=tdist->Value(j);
2431 if (tang->Value(j)>myG1Error)
2432 myG1Error=tang->Value(j);
2433 if (tcourb->Value(j)>myG2Error)
2434 myG2Error=tcourb->Value(j);
2435 Standard_Real U;
2436 if (myParCont->Value(i).Length()>3)
2437 U=(myParCont->Value(i).Value(j)+myParCont->Value(i).Value(j+1))/2;
2438 else
2439 U=LinCont->FirstParameter()+
2440 ( LinCont->LastParameter()-LinCont->FirstParameter())*(j-1)/(NbPts_i-2);
2441 Standard_Real diffDist=tdist->Value(j)-LinCont->G0Criterion(U),diffAng;
2442 if (LinCont->Order()>0)
2443 diffAng=tang->Value(j)-LinCont->G1Criterion(U);
2444 else diffAng=0;
2445 // find the maximum variation of error and calculate the average
2446 if (diffDist>0) {
2447 diffDist = diffDist/LinCont->G0Criterion(U);
2448 if (diffDist>diffDistMax)
2449 diffDistMax = diffDist;
2450 //SdiffDist+=diffDist;
2451 NdiffDist++;
2452 #ifdef DRAW
2453 if ((Affich) && (NbBoucle == myNbIter)) {
2454 gp_Pnt P;
2455 gp_Pnt2d P2d;
2456 LinCont->D0(U,P);
2457 Handle(Draw_Marker3D) mark =
2458 new (Draw_Marker3D)(P, Draw_X, Draw_orange);
2459 char name[256];
2460 sprintf(name,"mark_%d",++NbMark);
2461 Draw::Set(name, mark);
2462 if (!LinCont->ProjectedCurve().IsNull())
2463 P2d = LinCont->ProjectedCurve()->Value(U);
2464 else
2465 { if (!LinCont->Curve2dOnSurf().IsNull())
2466 P2d = LinCont->Curve2dOnSurf()->Value(U);
2467 else
2468 P2d = ProjectPoint(P);
2469 }
2470 sprintf(name,"mark2d_%d",++NbMark);
2471 Handle(Draw_Marker2D) mark2d =
2472 new (Draw_Marker2D)(P2d, Draw_X, Draw_orange);
2473 Draw::Set(name, mark2d);
2474 }
2475 #endif
2476 }
2477 else
2478 if ((diffAng>0)&&(LinCont->Order()==1)) {
2479 diffAng=diffAng/myLinCont->Value(i)->G1Criterion(U);
2480 if (diffAng>diffAngMax)
2481 diffAngMax = diffAng;
2482 //SdiffAng+=diffAng;
2483 NdiffAng++;
2484 #ifdef DRAW
2485 if ((Affich) && (NbBoucle == myNbIter)) {
2486 gp_Pnt P;
2487 LinCont->D0(U,P);
2488 Handle(Draw_Marker3D) mark =
2489 new Draw_Marker3D(P, Draw_X, Draw_or);
2490 char name[256];
2491 sprintf(name,"mark_%d",++NbMark);
2492 Draw::Set(name, mark);
2493 }
2494 #endif
2495 }
2496 }
2497
2498 if (NdiffDist>0) {// at least one point is not acceptable in G0
2499 Standard_Real Coef;
2500 if(LinCont->Order()== 0)
2501 Coef = 0.6*Log(diffDistMax+7.4);
2502 //7.4 corresponds to the calculation of min. coefficient = 1.2 is e^1.2/0.6
2503 else
2504 Coef = Log(diffDistMax+3.3);
2505 //3.3 corresponds to calculation of min. coefficient = 1.2 donc e^1.2
2506 if (Coef>3)
2507 Coef=3;
2508 //experimentally after the coefficient becomes bad for L cases
2509 if ((NbBoucle>1)&&(diffDistMax>2))
2510 { Coef=1.6;
2511 }
2512
2513 if (LinCont->NbPoints()>=Floor(LinCont->NbPoints()*Coef))
2514 Coef=2;// to provide increase of the number of points
2515
2516 LinCont->SetNbPoints(Standard_Integer(LinCont->NbPoints() * Coef));
2517 Result=Standard_False;
2518 }
2519 else
2520 if (NdiffAng>0) // at least 1 point is not acceptable in G1
2521 { Standard_Real Coef=1.5;
2522 if ((LinCont->NbPoints()+1)>=Floor(LinCont->NbPoints()*Coef))
2523 Coef=2;
2524
2525 LinCont->SetNbPoints(Standard_Integer(LinCont->NbPoints()*Coef )) ;
2526 Result=Standard_False;
2527 }
2528 }
2529 }
2530 if (!Result)
2531 {
2532 if (myFree && NbBoucle == 1)
2533 myPrevPlate = myPlate;
2534 myPlate.Init();
2535 }
2536 return Result;
2537 }
2538
2539
2540
2541 //---------------------------------------------------------
2542 // Function : VerifPoint
2543 //---------------------------------------------------------
2544 void GeomPlate_BuildPlateSurface::
VerifPoints(Standard_Real & Dist,Standard_Real & Ang,Standard_Real & Curv) const2545 VerifPoints (Standard_Real& Dist,
2546 Standard_Real& Ang,
2547 Standard_Real& Curv) const
2548 { Standard_Integer NTPntCont=myPntCont->Length();
2549 gp_Pnt Pi, Pf;
2550 gp_Pnt2d P2d;
2551 gp_Vec v1i, v1f, v2i, v2f, v3i, v3f;
2552 Ang=0;Dist=0,Curv=0;
2553 Handle(GeomPlate_PointConstraint) PntCont;
2554 for (Standard_Integer i=1;i<=NTPntCont;i++) {
2555 PntCont = myPntCont->Value(i);
2556 switch (PntCont->Order())
2557 { case 0 :
2558 P2d = PntCont->Pnt2dOnSurf();
2559 PntCont->D0(Pi);
2560 myGeomPlateSurface->D0( P2d.Coord(1), P2d.Coord(2), Pf);
2561 Dist = Pf.Distance(Pi);
2562 break;
2563 case 1 :
2564 PntCont->D1(Pi, v1i, v2i);
2565 P2d = PntCont->Pnt2dOnSurf();
2566 myGeomPlateSurface->D1( P2d.Coord(1), P2d.Coord(2), Pf, v1f, v2f);
2567 Dist = Pf.Distance(Pi);
2568 v3i = v1i^v2i; v3f=v1f^v2f;
2569 Ang=v3f.Angle(v3i);
2570 if (Ang>(M_PI/2))
2571 Ang = M_PI -Ang;
2572 break;
2573 case 2 :
2574 Handle(Geom_Surface) Splate (myGeomPlateSurface);
2575 LocalAnalysis_SurfaceContinuity CG2;
2576 P2d = PntCont->Pnt2dOnSurf();
2577 GeomLProp_SLProps Prop (Splate, P2d.Coord(1), P2d.Coord(2),
2578 2, 0.001);
2579 CG2.ComputeAnalysis(Prop, PntCont->LPropSurf(),GeomAbs_G2);
2580 Dist=CG2.C0Value();
2581 Ang=CG2.G1Angle();
2582 Curv=CG2.G2CurvatureGap();
2583 break;
2584 }
2585 }
2586 }
2587
ComputeAnisotropie() const2588 Standard_Real GeomPlate_BuildPlateSurface::ComputeAnisotropie() const
2589 {
2590 if (myAnisotropie)
2591 {
2592 //Temporary
2593 return 1.0;
2594 }
2595 else return 1.0;
2596 }
2597
IsOrderG1() const2598 Standard_Boolean GeomPlate_BuildPlateSurface::IsOrderG1() const
2599 {
2600 Standard_Boolean result = Standard_True;
2601 for (Standard_Integer i = 1; i <= myLinCont->Length(); i++)
2602 if (myLinCont->Value(i)->Order() < 1)
2603 {
2604 result = Standard_False;
2605 break;
2606 }
2607 return result;
2608 }
2609