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