0027309: Result of UnifySameDomain has InvalidMultiConnexity error
[occt.git] / src / GeomLib / GeomLib_CheckCurveOnSurface.cxx
CommitLineData
5adae760 1// Created by: Nikolai BUKHALOV
2// Copyright (c) 2015 OPEN CASCADE SAS
3//
4// This file is part of Open CASCADE Technology software library.
5//
6// This library is free software; you can redistribute it and/or modify it under
7// the terms of the GNU Lesser General Public License version 2.1 as published
8// by the Free Software Foundation, with special exception defined in the file
9// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10// distribution for complete text of the license and disclaimer of any warranty.
11//
12// Alternatively, this file may be used under the terms of Open CASCADE
13// commercial license or contractual agreement.
14
15#include <Adaptor2d_HCurve2d.hxx>
16#include <Adaptor3d_Curve.hxx>
17#include <Adaptor3d_CurveOnSurface.hxx>
18#include <Adaptor3d_HSurface.hxx>
19#include <Geom_BSplineCurve.hxx>
20#include <Geom_TrimmedCurve.hxx>
21#include <Geom2d_BSplineCurve.hxx>
22#include <Geom2d_TrimmedCurve.hxx>
23#include <Geom2dAdaptor_GHCurve.hxx>
24#include <GeomAdaptor_Curve.hxx>
25#include <GeomAdaptor_HSurface.hxx>
26#include <GeomLib_CheckCurveOnSurface.hxx>
27#include <gp_Pnt.hxx>
28#include <math_Matrix.hxx>
29#include <math_MultipleVarFunctionWithHessian.hxx>
30#include <math_NewtonMinimum.hxx>
31#include <math_PSO.hxx>
32#include <math_PSOParticlesPool.hxx>
33#include <OSD_Parallel.hxx>
34#include <Standard_ErrorHandler.hxx>
35#include <TColStd_Array1OfReal.hxx>
36
37class GeomLib_CheckCurveOnSurface_TargetFunc;
38
39static
40Standard_Boolean MinComputing(
41 GeomLib_CheckCurveOnSurface_TargetFunc& theFunction,
42 const Standard_Real theEpsilon, //1.0e-3
43 const Standard_Integer theNbParticles,
44 Standard_Real& theBestValue,
45 Standard_Real& theBestParameter);
46
47static Standard_Integer FillSubIntervals( const Handle(Geom_Curve)& theCurve3d,
48 const Handle(Geom2d_Curve)& theCurve2d,
49 const Standard_Real theFirst,
50 const Standard_Real theLast,
51 Standard_Integer &theNbParticles,
52 TColStd_Array1OfReal* const theSubIntervals = 0);
53
54//=======================================================================
55//class : GeomLib_CheckCurveOnSurface_TargetFunc
56//purpose : Target function (to be minimized)
57//=======================================================================
58class GeomLib_CheckCurveOnSurface_TargetFunc :
59 public math_MultipleVarFunctionWithHessian
60{
61 public:
62 GeomLib_CheckCurveOnSurface_TargetFunc( const Adaptor3d_Curve& theC3D,
63 const Adaptor3d_Curve& theAdCS,
64 const Standard_Real theFirst,
65 const Standard_Real theLast):
66 myCurve1(theC3D),
67 myCurve2(theAdCS),
68 myFirst(theFirst),
69 myLast(theLast)
70 {
71 }
72
73 //returns the number of parameters of the function
74 //(the function is one-dimension).
75 virtual Standard_Integer NbVariables() const {
76 return 1;
77 }
78
79 //returns value of the function when parameters are equal to theX
80 virtual Standard_Boolean Value(const math_Vector& theX,
81 Standard_Real& theFVal)
82 {
83 return Value(theX(1), theFVal);
84 }
85
86 //returns value of the one-dimension-function when parameter
87 //is equal to theX
88 Standard_Boolean Value( const Standard_Real theX,
89 Standard_Real& theFVal) const
90 {
91 try
92 {
93 OCC_CATCH_SIGNALS
94 if (!CheckParameter(theX))
95 return Standard_False;
96
97 const gp_Pnt aP1(myCurve1.Value(theX)),
98 aP2(myCurve2.Value(theX));
99
100 theFVal = -1.0*aP1.SquareDistance(aP2);
101 }
102 catch(Standard_Failure) {
103 return Standard_False;
104 }
105 //
106 return Standard_True;
107 }
108
109 //see analogical method for abstract owner class math_MultipleVarFunction
110 virtual Standard_Integer GetStateNumber()
111 {
112 return 0;
113 }
114
115 //returns the gradient of the function when parameters are
116 //equal to theX
117 virtual Standard_Boolean Gradient(const math_Vector& theX,
118 math_Vector& theGrad)
119 {
120 return Derive(theX(1), theGrad(1));
121 }
122
123 //returns 1st derivative of the the one-dimension-function when
124 //parameter is equal to theX
125 Standard_Boolean Derive(const Standard_Real theX, Standard_Real& theDeriv) const
126 {
127 try
128 {
129 OCC_CATCH_SIGNALS
130 if (!CheckParameter(theX))
131 {
132 return Standard_False;
133 }
134 //
135 gp_Pnt aP1, aP2;
136 gp_Vec aDC1, aDC2;
137 //
138 myCurve1.D1(theX, aP1, aDC1);
139 myCurve2.D1(theX, aP2, aDC2);
140
141 const gp_Vec aVec1(aP1, aP2), aVec2(aDC2-aDC1);
142 //
143 theDeriv = -2.0*aVec1.Dot(aVec2);
144 }
145 catch(Standard_Failure)
146 {
147 return Standard_False;
148 }
149
150 return Standard_True;
151 }
152
153 //returns value and gradient
154 virtual Standard_Boolean Values(const math_Vector& theX,
155 Standard_Real& theVal,
156 math_Vector& theGrad)
157 {
158 if (!Value(theX, theVal))
159 {
160 return Standard_False;
161 }
162 //
163 if (!Gradient(theX, theGrad)) {
164 return Standard_False;
165 }
166 //
167 return Standard_True;
168 }
169
170 //returns value, gradient and hessian
171 virtual Standard_Boolean Values(const math_Vector& theX,
172 Standard_Real& theVal,
173 math_Vector& theGrad,
174 math_Matrix& theHessian)
175 {
176 if (!Value(theX, theVal))
177 {
178 return Standard_False;
179 }
180 //
181 if (!Gradient(theX, theGrad))
182 {
183 return Standard_False;
184 }
185 //
186 theHessian(1,1) = theGrad(1);
187 //
188 return Standard_True;
189 }
190 //
191 Standard_Real FirstParameter() const
192 {
193 return myFirst;
194 }
195
196 //
197 Standard_Real LastParameter() const
198 {
199 return myLast;
200 }
201
202 private:
203 GeomLib_CheckCurveOnSurface_TargetFunc operator=(GeomLib_CheckCurveOnSurface_TargetFunc&);
204
205 //checks if the function can be computed when its parameter is
206 //equal to theParam
207 Standard_Boolean CheckParameter(const Standard_Real theParam) const
208 {
209 return ((myFirst <= theParam) && (theParam <= myLast));
210 }
211
212 const Adaptor3d_Curve& myCurve1;
213 const Adaptor3d_Curve& myCurve2;
214 const Standard_Real myFirst;
215 const Standard_Real myLast;
216};
217
218//=======================================================================
219//class : GeomLib_CheckCurveOnSurface_Local
220//purpose : Created for parallelization possibility only
221//=======================================================================
222class GeomLib_CheckCurveOnSurface_Local
223{
224public:
225 GeomLib_CheckCurveOnSurface_Local(
226 const Handle(Geom_Curve)& theCurve3D,
227 const Handle(Geom2d_Curve)& theCurve2D,
228 const Handle(Geom_Surface)& theSurface,
229 const TColStd_Array1OfReal& theIntervalsArr,
230 const Standard_Real theEpsilonRange,
231 const Standard_Integer theNbParticles):
232 myCurve3D(theCurve3D),
233 myCurve2D(theCurve2D),
234 mySurface(theSurface),
235 mySubIntervals(theIntervalsArr),
236 myEpsilonRange(theEpsilonRange),
237 myNbParticles(theNbParticles),
238 myArrOfDist(theIntervalsArr.Lower(), theIntervalsArr.Upper()-1),
239 myArrOfParam(theIntervalsArr.Lower(), theIntervalsArr.Upper()-1)
240 {
241 }
242
243 void operator()(const Standard_Integer& theIndex) const
244 {
245 //For every sub-interval (which is set by mySubIntervals array) this method
246 //computes optimal value of GeomLib_CheckCurveOnSurface_TargetFunc function.
247 //This optimal value will be put in corresponding (depending on theIndex - the
248 //identificator of the current interval in mySubIntervals array) cell of
249 //myArrOfDist and myArrOfParam arrays.
250 const GeomAdaptor_Curve anAC(myCurve3D);
251 const Handle(Adaptor2d_HCurve2d) anAd2dC = new Geom2dAdaptor_GHCurve(myCurve2D);
252 const Handle(Adaptor3d_HSurface) anAdS = new GeomAdaptor_HSurface(mySurface);
253
254 const Adaptor3d_CurveOnSurface anACS(anAd2dC, anAdS);
255
256 GeomLib_CheckCurveOnSurface_TargetFunc aFunc( anAC, anACS,
257 mySubIntervals.Value(theIndex),
258 mySubIntervals.Value(theIndex+1));
259
260 Standard_Real aMinDist = RealLast(), aPar = 0.0;
261 if(!MinComputing(aFunc, myEpsilonRange, myNbParticles, aMinDist, aPar))
262 {
263 myArrOfDist(theIndex) = RealLast();
264 myArrOfParam(theIndex) = aFunc.FirstParameter();
265 return;
266 }
267
268 myArrOfDist(theIndex) = aMinDist;
269 myArrOfParam(theIndex) = aPar;
270 }
271
272 //Returns optimal value (inverse of square of maximal distance)
273 void OptimalValues(Standard_Real& theMinimalValue, Standard_Real& theParameter) const
274 {
275 //This method looks for the minimal value of myArrOfDist.
276
277 const Standard_Integer aStartInd = myArrOfDist.Lower();
278 theMinimalValue = myArrOfDist(aStartInd);
279 theParameter = myArrOfParam(aStartInd);
280 for(Standard_Integer i = aStartInd + 1; i <= myArrOfDist.Upper(); i++)
281 {
282 if(myArrOfDist(i) < theMinimalValue)
283 {
284 theMinimalValue = myArrOfDist(i);
285 theParameter = myArrOfParam(i);
286 }
287 }
288 }
289
290private:
291 GeomLib_CheckCurveOnSurface_Local operator=(GeomLib_CheckCurveOnSurface_Local&);
292 const Handle(Geom_Curve)& myCurve3D;
293 const Handle(Geom2d_Curve)& myCurve2D;
294 const Handle(Geom_Surface)& mySurface;
295
296 const TColStd_Array1OfReal& mySubIntervals;
297 const Standard_Real myEpsilonRange;
298 const Standard_Integer myNbParticles;
299 mutable NCollection_Array1<Standard_Real> myArrOfDist;
300 mutable NCollection_Array1<Standard_Real> myArrOfParam;
301};
302
303//=======================================================================
304//function : GeomLib_CheckCurveOnSurface
305//purpose :
306//=======================================================================
307GeomLib_CheckCurveOnSurface::GeomLib_CheckCurveOnSurface()
308:
309 myFirst(0.),
310 myLast(0.),
311 myErrorStatus(0),
312 myMaxDistance(RealLast()),
6ca1c746 313 myMaxParameter(0.),
314 myTolRange(Precision::PConfusion())
5adae760 315{
316}
317
318//=======================================================================
319//function : GeomLib_CheckCurveOnSurface
320//purpose :
321//=======================================================================
322GeomLib_CheckCurveOnSurface::
323 GeomLib_CheckCurveOnSurface(const Handle(Geom_Curve)& theCurve,
324 const Handle(Geom_Surface)& theSurface,
325 const Standard_Real theFirst,
6ca1c746 326 const Standard_Real theLast,
327 const Standard_Real theTolRange):
5adae760 328 myCurve(theCurve),
329 mySurface(theSurface),
330 myFirst(theFirst),
331 myLast(theLast),
332 myErrorStatus(0),
333 myMaxDistance(RealLast()),
6ca1c746 334 myMaxParameter(0.),
335 myTolRange(theTolRange)
5adae760 336{
337}
338
339//=======================================================================
340//function : Init
341//purpose :
342//=======================================================================
343void GeomLib_CheckCurveOnSurface::Init()
344{
345 myCurve.Nullify();
346 mySurface.Nullify();
347 myFirst = 0.0;
348 myLast = 0.0;
349 myErrorStatus = 0;
350 myMaxDistance = RealLast();
351 myMaxParameter = 0.0;
6ca1c746 352 myTolRange = Precision::PConfusion();
5adae760 353}
354
355//=======================================================================
356//function : Init
357//purpose :
358//=======================================================================
359void GeomLib_CheckCurveOnSurface::Init( const Handle(Geom_Curve)& theCurve,
360 const Handle(Geom_Surface)& theSurface,
361 const Standard_Real theFirst,
6ca1c746 362 const Standard_Real theLast,
363 const Standard_Real theTolRange)
5adae760 364{
365 myCurve = theCurve;
366 mySurface = theSurface;
367 myFirst = theFirst;
368 myLast = theLast;
369 myErrorStatus = 0;
370 myMaxDistance = RealLast();
371 myMaxParameter = 0.0;
6ca1c746 372 myTolRange = theTolRange;
5adae760 373}
374
375//=======================================================================
376//function : Perform
377//purpose :
378//=======================================================================
5adae760 379#ifndef HAVE_TBB
380//After fixing bug # 26365, this fragment should be deleted
381//(together the text "#ifdef HAVE_TBB")
382
383void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
384 const Standard_Boolean)
385{
386 const Standard_Boolean isTheMTDisabled = Standard_True;
387#else
388void GeomLib_CheckCurveOnSurface::Perform(const Handle(Geom2d_Curve)& thePCurve,
389 const Standard_Boolean isTheMTDisabled)
390{
391#endif
392 if( myCurve.IsNull() ||
393 mySurface.IsNull() ||
394 thePCurve.IsNull())
395 {
396 myErrorStatus = 1;
397 return;
398 }
399
6ca1c746 400 if(((myCurve->FirstParameter() - myFirst) > myTolRange) ||
401 ((myCurve->LastParameter() - myLast) < -myTolRange) ||
402 ((thePCurve->FirstParameter() - myFirst) > myTolRange) ||
403 ((thePCurve->LastParameter() - myLast) < -myTolRange))
5adae760 404 {
405 myErrorStatus = 2;
406 return;
407 }
408
409 const Standard_Real anEpsilonRange = 1.e-3;
410
411 Standard_Integer aNbParticles = 3;
412
413 //Polynomial function with degree n has not more than n-1 maxima and
414 //minima (degree of 1st derivative is equal to n-1 => 1st derivative has
415 //no greater than n-1 roots). Consequently, this function has
416 //maximum n monotonicity intervals. That is a good idea to try to put
417 //at least one particle in every monotonicity interval. Therefore,
418 //number of particles should be equal to n.
419
420 const Standard_Integer aNbSubIntervals =
421 FillSubIntervals( myCurve, thePCurve,
422 myFirst, myLast, aNbParticles);
423
424 if(!aNbSubIntervals)
425 {
426 myErrorStatus = 3;
427 return;
428 }
429
430 try {
431 OCC_CATCH_SIGNALS
432
433 TColStd_Array1OfReal anIntervals(1, aNbSubIntervals+1);
434 FillSubIntervals(myCurve, thePCurve, myFirst, myLast, aNbParticles, &anIntervals);
435
436 GeomLib_CheckCurveOnSurface_Local aComp(myCurve, thePCurve,
437 mySurface, anIntervals, anEpsilonRange, aNbParticles);
438
439 OSD_Parallel::For(anIntervals.Lower(), anIntervals.Upper(), aComp, isTheMTDisabled);
440
441 aComp.OptimalValues(myMaxDistance, myMaxParameter);
442
443 myMaxDistance = sqrt(Abs(myMaxDistance));
444 }
445 catch (Standard_Failure) {
446 myErrorStatus = 3;
447 }
448}
449
450//=======================================================================
451// Function : FillSubIntervals
452// purpose : Divides [theFirst, theLast] interval on parts
453// in order to make searching-algorithm more precisely
454// (fills theSubIntervals array).
455// Returns number of subintervals.
456//=======================================================================
457Standard_Integer FillSubIntervals(const Handle(Geom_Curve)& theCurve3d,
458 const Handle(Geom2d_Curve)& theCurve2d,
459 const Standard_Real theFirst,
460 const Standard_Real theLast,
461 Standard_Integer &theNbParticles,
462 TColStd_Array1OfReal* const theSubIntervals)
463{
464 const Standard_Real anArrTempC[2] = {theFirst, theLast};
465 const TColStd_Array1OfReal anArrTemp(anArrTempC[0], 1, 2);
466
467 theNbParticles = 3;
468 Handle(Geom2d_BSplineCurve) aBS2DCurv;
469 Handle(Geom_BSplineCurve) aBS3DCurv;
470
471 //
472 if (theCurve3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
473 {
474 aBS3DCurv = Handle(Geom_BSplineCurve)::
475 DownCast(Handle(Geom_TrimmedCurve)::
476 DownCast(theCurve3d)->BasisCurve());
477 }
478 else
479 {
480 aBS3DCurv = Handle(Geom_BSplineCurve)::DownCast(theCurve3d);
481 }
482
483 if (theCurve2d->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)))
484 {
485 aBS2DCurv = Handle(Geom2d_BSplineCurve)::
486 DownCast(Handle(Geom2d_TrimmedCurve)::
487 DownCast(theCurve2d)->BasisCurve());
488 }
489 else
490 {
491 aBS2DCurv = Handle(Geom2d_BSplineCurve)::DownCast(theCurve2d);
492 }
493
494 const TColStd_Array1OfReal &anArrKnots3D = !aBS3DCurv.IsNull() ?
495 aBS3DCurv->Knots() :
496 anArrTemp;
497 const TColStd_Array1OfReal &anArrKnots2D = !aBS2DCurv.IsNull() ?
498 aBS2DCurv->Knots() :
499 anArrTemp;
500
501 Standard_Integer aNbSubIntervals = 1;
502
503 try
504 {
505 OCC_CATCH_SIGNALS
506 const Standard_Integer anIndMax3D = anArrKnots3D.Upper(),
507 anIndMax2D = anArrKnots2D.Upper();
508
509 Standard_Integer anIndex3D = anArrKnots3D.Lower(),
510 anIndex2D = anArrKnots2D.Lower();
511
512 if(theSubIntervals)
513 theSubIntervals->ChangeValue(aNbSubIntervals) = theFirst;
514
515 while((anIndex3D <= anIndMax3D) && (anIndex2D <= anIndMax2D))
516 {
517 const Standard_Real aVal3D = anArrKnots3D.Value(anIndex3D),
518 aVal2D = anArrKnots2D.Value(anIndex2D);
519 const Standard_Real aDelta = aVal3D - aVal2D;
520
521 if(aDelta < Precision::PConfusion())
522 {//aVal3D <= aVal2D
523 if((aVal3D > theFirst) && (aVal3D < theLast))
524 {
525 aNbSubIntervals++;
526
527 if(theSubIntervals)
528 theSubIntervals->ChangeValue(aNbSubIntervals) = aVal3D;
529 }
530
531 anIndex3D++;
532
533 if(-aDelta < Precision::PConfusion())
534 {//aVal3D == aVal2D
535 anIndex2D++;
536 }
537 }
538 else
539 {//aVal2D < aVal3D
540 if((aVal2D > theFirst) && (aVal2D < theLast))
541 {
542 aNbSubIntervals++;
543
544 if(theSubIntervals)
545 theSubIntervals->ChangeValue(aNbSubIntervals) = aVal2D;
546 }
547
548 anIndex2D++;
549 }
550 }
551
552 if(theSubIntervals)
553 theSubIntervals->ChangeValue(aNbSubIntervals+1) = theLast;
554
555 if(!aBS3DCurv.IsNull())
556 {
557 theNbParticles = Max(theNbParticles, aBS3DCurv->Degree());
558 }
559
560 if(!aBS2DCurv.IsNull())
561 {
562 theNbParticles = Max(theNbParticles, aBS2DCurv->Degree());
563 }
564 }
565 catch(Standard_Failure)
566 {
567#ifdef OCCT_DEBUG
568 cout << "ERROR! BRepLib_CheckCurveOnSurface.cxx, "
569 "FillSubIntervals(): Incorrect filling!" << endl;
570#endif
571
572 aNbSubIntervals = 0;
573 }
574
575 return aNbSubIntervals;
576}
577
578//=======================================================================
579//class : MinComputing
580//purpose : Performs computing minimal value
581//=======================================================================
582Standard_Boolean MinComputing (
583 GeomLib_CheckCurveOnSurface_TargetFunc& theFunction,
584 const Standard_Real theEpsilon, //1.0e-3
585 const Standard_Integer theNbParticles,
586 Standard_Real& theBestValue,
587 Standard_Real& theBestParameter)
588{
589 try
590 {
591 OCC_CATCH_SIGNALS
592
593 //They are used for finding a position of theNbParticles worst places
594 const Standard_Integer aNbControlPoints = 3*theNbParticles;
595 //
596 math_Vector aParInf(1, 1), aParSup(1, 1), anOutputParam(1, 1), aStepPar(1,1);
597 aParInf(1) = theFunction.FirstParameter();
598 aParSup(1) = theFunction.LastParameter();
599 theBestParameter = aParInf(1);
600 theBestValue = RealLast();
601
602 const Standard_Real aDeltaParam = aParSup(1) - aParInf(1);
603 if(aDeltaParam < Precision::PConfusion())
604 return Standard_False;
605
606 aStepPar(1) = theEpsilon*aDeltaParam;
607
608 math_PSOParticlesPool aParticles(theNbParticles, 1);
609
610 const Standard_Real aStep = aDeltaParam/(aNbControlPoints-1);
611 Standard_Integer aCount = 1;
612 for(Standard_Real aPrm = aParInf(1); aCount <= aNbControlPoints; aCount++,
613 aPrm = (aCount == aNbControlPoints)? aParSup(1) : aPrm+aStep)
614 {
615 Standard_Real aVal = RealLast();
616 theFunction.Value(aPrm, aVal);
617
618 PSO_Particle* aParticle = aParticles.GetWorstParticle();
619
620 if(aVal > aParticle->BestDistance)
621 continue;
622
623 aParticle->Position[0] = aPrm;
624 aParticle->BestPosition[0] = aPrm;
625 aParticle->Distance = aVal;
626 aParticle->BestDistance = aVal;
627 }
628
629 math_PSO aPSO(&theFunction, aParInf, aParSup, aStepPar);
630 aPSO.Perform(aParticles, theNbParticles, theBestValue, anOutputParam);
631
632 //Here, anOutputParam contains parameter, which is near to optimal.
633 //It needs to be more precise. Precision is made by math_NewtonMinimum.
634 math_NewtonMinimum anA(theFunction);
635 anA.Perform(theFunction, anOutputParam);
636
637 if(!anA.IsDone())
638 {
639#ifdef OCCT_DEBUG
640 cout << "BRepLib_CheckCurveOnSurface::Compute(): No solution found!" << endl;
641#endif
642 return Standard_False;
643 }
644
645 anA.Location(anOutputParam);
646 theBestParameter = anOutputParam(1);
647 theBestValue = anA.Minimum();
648 }
649 catch(Standard_Failure)
650 {
651#ifdef OCCT_DEBUG
652 cout << "BRepLib_CheckCurveOnSurface.cxx: Exception in MinComputing()!" << endl;
653#endif
654 return Standard_False;
655 }
656
657 return Standard_True;
658}