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