Warning, /include/opencascade/ApproxInt_ImpPrmSvSurfaces.gxx is written in an unsupported language. File is not indexed.
0001 // Created on: 1993-03-17
0002 // Created by: Laurent BUCHARD
0003 // Copyright (c) 1993-1999 Matra Datavision
0004 // Copyright (c) 1999-2014 OPEN CASCADE SAS
0005 //
0006 // This file is part of Open CASCADE Technology software library.
0007 //
0008 // This library is free software; you can redistribute it and/or modify it under
0009 // the terms of the GNU Lesser General Public License version 2.1 as published
0010 // by the Free Software Foundation, with special exception defined in the file
0011 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0012 // distribution for complete text of the license and disclaimer of any warranty.
0013 //
0014 // Alternatively, this file may be used under the terms of Open CASCADE
0015 // commercial license or contractual agreement.
0016
0017 #include <IntSurf_PntOn2S.hxx>
0018 #include <TColStd_Array1OfReal.hxx>
0019 #include <math_FunctionSetRoot.hxx>
0020 #include <StdFail_NotDone.hxx>
0021 #include <GeomAbs_SurfaceType.hxx>
0022 #include <Precision.hxx>
0023
0024 //=======================================================================
0025 //function : IsSingular
0026 //purpose : Returns TRUE if vectors theDU || theDV or if at least one
0027 // of them has null-magnitude.
0028 // theSqLinTol is square of linear tolerance.
0029 // theAngTol is angular tolerance.
0030 //=======================================================================
0031 static Standard_Boolean IsSingular( const gp_Vec& theDU,
0032 const gp_Vec& theDV,
0033 const Standard_Real theSqLinTol,
0034 const Standard_Real theAngTol)
0035 {
0036 gp_Vec aDU(theDU), aDV(theDV);
0037
0038 const Standard_Real aSqMagnDU = aDU.SquareMagnitude(),
0039 aSqMagnDV = aDV.SquareMagnitude();
0040
0041 if(aSqMagnDU < theSqLinTol)
0042 return Standard_True;
0043
0044 aDU.Divide(sqrt(aSqMagnDU));
0045
0046 if(aSqMagnDV < theSqLinTol)
0047 return Standard_True;
0048
0049 aDV.Divide(sqrt(aSqMagnDV));
0050
0051 //Here aDU and aDV vectors have magnitude 1.0.
0052
0053 if(aDU.Crossed(aDV).SquareMagnitude() < theAngTol*theAngTol)
0054 return Standard_True;
0055
0056 return Standard_False;
0057 }
0058
0059 //=======================================================================
0060 //function : SingularProcessing
0061 //purpose : Computes 2D-representation (in UV-coordinates) of
0062 // theTg3D vector on the surface in case when
0063 // theDU.Crossed(theDV).Magnitude() == 0.0. Stores result in
0064 // theTg2D variable.
0065 // theDU and theDV are vectors of 1st derivative
0066 // (with respect to U and V variables correspondingly).
0067 // If theIsTo3DTgCompute == TRUE then theTg3D has not been
0068 // defined yet (it should be computed).
0069 // theLinTol is SQUARE of the tolerance.
0070 //
0071 //Algorithm:
0072 // Condition
0073 // Tg=theDU*theTg2D.X()+theDV*theTg2D.Y()
0074 // has to be satisfied strictly.
0075 // More over, vector Tg has to be NORMALIZED
0076 // (if theIsTo3DTgCompute == TRUE then new computed vector will
0077 // always have magnitude 1.0).
0078 //=======================================================================
0079 static Standard_Boolean SingularProcessing( const gp_Vec& theDU,
0080 const gp_Vec& theDV,
0081 const Standard_Boolean theIsTo3DTgCompute,
0082 const Standard_Real theLinTol,
0083 const Standard_Real theAngTol,
0084 gp_Vec& theTg3D,
0085 gp_Vec2d& theTg2D)
0086 {
0087 //Attention: @ \sin theAngTol \approx theAngTol @ (for cross-product)
0088
0089 //Really, vector theTg3D has to be normalized (if theIsTo3DTgCompute == FALSE).
0090 const Standard_Real aSQTan = theTg3D.SquareMagnitude();
0091
0092 const Standard_Real aSqMagnDU = theDU.SquareMagnitude(),
0093 aSqMagnDV = theDV.SquareMagnitude();
0094
0095 //There are some reasons of singularity
0096
0097 //1.
0098 if((aSqMagnDU < theLinTol) && (aSqMagnDV < theLinTol))
0099 {
0100 //For future, this case can be processed as same as in case of
0101 //osculating surfaces (expanding in Taylor series). Here,
0102 //we return only.
0103
0104 return Standard_False;
0105 }
0106
0107 //2.
0108 if(aSqMagnDU < theLinTol)
0109 {
0110 //In this case, theTg3D vector will be parallel with theDV.
0111 //Its true direction shall be precised later (the algorithm is
0112 //based on array of Walking-points).
0113
0114 if(theIsTo3DTgCompute)
0115 {
0116 //theTg3D will be normalized. Its magnitude is
0117 const Standard_Real aTgMagn = 1.0;
0118
0119 const Standard_Real aNorm = sqrt(aSqMagnDV);
0120 theTg3D = theDV.Divided(aNorm);
0121 theTg2D.SetCoord(0.0, aTgMagn/aNorm);
0122 }
0123 else
0124 {
0125 //theTg3D is already defined.
0126 //Here we check only, if this tangent is parallel to theDV.
0127
0128 if(theDV.Crossed(theTg3D).SquareMagnitude() <
0129 theAngTol*theAngTol*aSqMagnDV*aSQTan)
0130 {
0131 //theTg3D is parallel to theDV
0132
0133 //Use sign "+" if theTg3D and theDV are codirectional
0134 //and sign "-" if opposite
0135 const Standard_Real aDP = theTg3D.Dot(theDV);
0136 theTg2D.SetCoord(0.0, Sign(sqrt(aSQTan/aSqMagnDV), aDP));
0137 }
0138 else
0139 {
0140 //theTg3D is not parallel to theDV
0141 //It is abnormal
0142
0143 return Standard_False;
0144 }
0145 }
0146
0147 return Standard_True;
0148 }
0149
0150 //3.
0151 if(aSqMagnDV < theLinTol)
0152 {
0153 //In this case, theTg3D vector will be parallel with theDU.
0154 //Its true direction shall be precised later (the algorithm is
0155 //based on array of Walking-points).
0156
0157 if(theIsTo3DTgCompute)
0158 {
0159 //theTg3D will be normalized. Its magnitude is
0160 const Standard_Real aTgMagn = 1.0;
0161
0162 const Standard_Real aNorm = sqrt(aSqMagnDU);
0163 theTg3D = theDU.Divided(aNorm);
0164 theTg2D.SetCoord(aTgMagn/aNorm, 0.0);
0165 }
0166 else
0167 {
0168 //theTg3D is already defined.
0169 //Here we check only, if this tangent is parallel to theDU.
0170
0171 if(theDU.Crossed(theTg3D).SquareMagnitude() <
0172 theAngTol*theAngTol*aSqMagnDU*aSQTan)
0173 {
0174 //theTg3D is parallel to theDU
0175
0176 //Use sign "+" if theTg3D and theDU are codirectional
0177 //and sign "-" if opposite
0178 const Standard_Real aDP = theTg3D.Dot(theDU);
0179 theTg2D.SetCoord(Sign(sqrt(aSQTan/aSqMagnDU), aDP), 0.0);
0180 }
0181 else
0182 {
0183 //theTg3D is not parallel to theDU
0184 //It is abnormal
0185
0186 return Standard_False;
0187 }
0188 }
0189
0190 return Standard_True;
0191 }
0192
0193 //4. If aSqMagnDU > 0.0 && aSqMagnDV > 0.0 but theDV || theDU.
0194
0195 const Standard_Real aLenU = sqrt(aSqMagnDU),
0196 aLenV = sqrt(aSqMagnDV);
0197
0198 //aLenSum > 0.0 definitely
0199 const Standard_Real aLenSum = aLenU + aLenV;
0200
0201 if(theDV.Dot(theDU) > 0.0)
0202 {
0203 //Vectors theDV and theDU are codirectional.
0204
0205 if(theIsTo3DTgCompute)
0206 {
0207 theTg2D.SetCoord(1.0/aLenSum, 1.0/aLenSum);
0208 theTg3D = theDU*theTg2D.X() + theDV*theTg2D.Y();
0209 }
0210 else
0211 {
0212 //theTg3D is already defined.
0213 //Here we check only, if this tangent is parallel to theDU
0214 //(and theDV together).
0215
0216 if(theDU.Crossed(theTg3D).SquareMagnitude() <
0217 theAngTol*theAngTol*aSqMagnDU*aSQTan)
0218 {
0219 //theTg3D is parallel to theDU
0220
0221 const Standard_Real aDP = theTg3D.Dot(theDU);
0222 const Standard_Real aLenTg = Sign(sqrt(aSQTan), aDP);
0223 theTg2D.SetCoord(aLenTg/aLenSum, aLenTg/aLenSum);
0224 }
0225 else
0226 {
0227 //theTg3D is not parallel to theDU
0228 //It is abnormal
0229
0230 return Standard_False;
0231 }
0232 }
0233 }
0234 else
0235 {
0236 //Vectors theDV and theDU are opposite.
0237
0238 if(theIsTo3DTgCompute)
0239 {
0240 //Here we chose theDU as direction of theTg3D.
0241 //True direction shall be precised later (the algorithm is
0242 //based on array of Walking-points).
0243
0244 theTg2D.SetCoord(1.0/aLenSum, -1.0/aLenSum);
0245 theTg3D = theDU*theTg2D.X() + theDV*theTg2D.Y();
0246 }
0247 else
0248 {
0249 //theTg3D is already defined.
0250 //Here we check only, if this tangent is parallel to theDU
0251 //(and theDV together).
0252
0253 if(theDU.Crossed(theTg3D).SquareMagnitude() <
0254 theAngTol*theAngTol*aSqMagnDU*aSQTan)
0255 {
0256 //theTg3D is parallel to theDU
0257
0258 const Standard_Real aDP = theTg3D.Dot(theDU);
0259 const Standard_Real aLenTg = Sign(sqrt(aSQTan), aDP);
0260 theTg2D.SetCoord(aLenTg/aLenSum, -aLenTg/aLenSum);
0261 }
0262 else
0263 {
0264 //theTg3D is not parallel to theDU
0265 //It is abnormal
0266
0267 return Standard_False;
0268 }
0269 }
0270 }
0271
0272 return Standard_True;
0273 }
0274
0275 //=======================================================================
0276 //function : NonSingularProcessing
0277 //purpose : Computes 2D-representation (in UV-coordinates) of
0278 // theTg3D vector on the surface in case when
0279 // theDU.Crossed(theDV).Magnitude() > 0.0. Stores result in
0280 // theTg2D variable.
0281 // theDU and theDV are vectors of 1st derivative
0282 // (with respect to U and V variables correspondingly).
0283 // theLinTol is SQUARE of the tolerance.
0284 //
0285 //Algorithm:
0286 // Condition
0287 // Tg=theDU*theTg2D.X()+theDV*theTg2D.Y()
0288 // has to be satisfied strictly.
0289 // More over, vector Tg has always to be NORMALIZED.
0290 //=======================================================================
0291 static Standard_Boolean NonSingularProcessing(const gp_Vec& theDU,
0292 const gp_Vec& theDV,
0293 const gp_Vec& theTg3D,
0294 const Standard_Real theLinTol,
0295 const Standard_Real theAngTol,
0296 gp_Vec2d& theTg2D)
0297 {
0298 const gp_Vec aNormal = theDU.Crossed(theDV);
0299 const Standard_Real aSQMagn = aNormal.SquareMagnitude();
0300
0301 if(IsSingular(theDU, theDV, theLinTol, theAngTol))
0302 {
0303 gp_Vec aTg(theTg3D);
0304 return
0305 SingularProcessing(theDU, theDV, Standard_False,
0306 theLinTol, theAngTol, aTg, theTg2D);
0307 }
0308
0309 //If @\vec{T}=\vec{A}*U+\vec{B}*V@ then
0310
0311 // \left\{\begin{matrix}
0312 // \vec{A} \times \vec{T} = (\vec{A} \times \vec{B})*V
0313 // \vec{B} \times \vec{T} = (\vec{B} \times \vec{A})*U
0314 // \end{matrix}\right.
0315
0316 //From here, values of U and V can be found very easily
0317 //(if @\left \| \vec{A} \times \vec{B} \right \| > 0.0 @,
0318 //else it is singular case).
0319
0320 const gp_Vec aTgU(theTg3D.Crossed(theDU)), aTgV(theTg3D.Crossed(theDV));
0321 const Standard_Real aDeltaU = aTgV.SquareMagnitude()/aSQMagn;
0322 const Standard_Real aDeltaV = aTgU.SquareMagnitude()/aSQMagn;
0323
0324 theTg2D.SetCoord(Sign(sqrt(aDeltaU), aTgV.Dot(aNormal)), -Sign(sqrt(aDeltaV), aTgU.Dot(aNormal)));
0325
0326 return Standard_True;
0327 }
0328
0329 //--------------------------------------------------------------------------------
0330 ApproxInt_ImpPrmSvSurfaces::ApproxInt_ImpPrmSvSurfaces( const TheISurface& ISurf
0331 ,const ThePSurface& PSurf):
0332 MyIsTangent(Standard_False),
0333 MyHasBeenComputed(Standard_False),
0334 MyIsTangentbis(Standard_False),
0335 MyHasBeenComputedbis(Standard_False),
0336 MyImplicitFirst(Standard_True),
0337 MyZerImpFunc(PSurf,ISurf)
0338 {
0339 SetUseSolver(Standard_True);
0340 }
0341 //--------------------------------------------------------------------------------
0342 ApproxInt_ImpPrmSvSurfaces::ApproxInt_ImpPrmSvSurfaces( const ThePSurface& PSurf
0343 ,const TheISurface& ISurf):
0344 MyIsTangent(Standard_False),
0345 MyHasBeenComputed(Standard_False),
0346 MyIsTangentbis(Standard_False),
0347 MyHasBeenComputedbis(Standard_False),
0348 MyImplicitFirst(Standard_False),
0349 MyZerImpFunc(PSurf,ISurf)
0350 {
0351 SetUseSolver(Standard_True);
0352 }
0353 //--------------------------------------------------------------------------------
0354 void ApproxInt_ImpPrmSvSurfaces::Pnt(const Standard_Real u1,
0355 const Standard_Real v1,
0356 const Standard_Real u2,
0357 const Standard_Real v2,
0358 gp_Pnt& P) {
0359 gp_Pnt aP;
0360 gp_Vec aT;
0361 gp_Vec2d aTS1,aTS2;
0362 Standard_Real tu1=u1;
0363 Standard_Real tu2=u2;
0364 Standard_Real tv1=v1;
0365 Standard_Real tv2=v2;
0366 this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
0367 P=MyPnt;
0368 }
0369 //--------------------------------------------------------------------------------
0370 Standard_Boolean ApproxInt_ImpPrmSvSurfaces::Tangency(const Standard_Real u1,
0371 const Standard_Real v1,
0372 const Standard_Real u2,
0373 const Standard_Real v2,
0374 gp_Vec& T) {
0375 gp_Pnt aP;
0376 gp_Vec aT;
0377 gp_Vec2d aTS1,aTS2;
0378 Standard_Real tu1=u1;
0379 Standard_Real tu2=u2;
0380 Standard_Real tv1=v1;
0381 Standard_Real tv2=v2;
0382 Standard_Boolean t=this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
0383 T=MyTg;
0384 return(t);
0385 }
0386 //--------------------------------------------------------------------------------
0387 Standard_Boolean ApproxInt_ImpPrmSvSurfaces::TangencyOnSurf1(const Standard_Real u1,
0388 const Standard_Real v1,
0389 const Standard_Real u2,
0390 const Standard_Real v2,
0391 gp_Vec2d& T) {
0392 gp_Pnt aP;
0393 gp_Vec aT;
0394 gp_Vec2d aTS1,aTS2;
0395 Standard_Real tu1=u1;
0396 Standard_Real tu2=u2;
0397 Standard_Real tv1=v1;
0398 Standard_Real tv2=v2;
0399 Standard_Boolean t=this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
0400 T=MyTguv1;
0401 return(t);
0402 }
0403 //--------------------------------------------------------------------------------
0404 Standard_Boolean ApproxInt_ImpPrmSvSurfaces::TangencyOnSurf2(const Standard_Real u1,
0405 const Standard_Real v1,
0406 const Standard_Real u2,
0407 const Standard_Real v2,
0408 gp_Vec2d& T) {
0409 gp_Pnt aP;
0410 gp_Vec aT;
0411 gp_Vec2d aTS1,aTS2;
0412 Standard_Real tu1=u1;
0413 Standard_Real tu2=u2;
0414 Standard_Real tv1=v1;
0415 Standard_Real tv2=v2;
0416 Standard_Boolean t=this->Compute(tu1,tv1,tu2,tv2,aP,aT,aTS1,aTS2);
0417 T=MyTguv2;
0418 return(t);
0419 }
0420
0421 //=======================================================================
0422 //function : Compute
0423 //purpose : Computes point on curve, 3D and 2D-tangents of a curve and
0424 // parameters on the surfaces.
0425 //=======================================================================
0426 Standard_Boolean ApproxInt_ImpPrmSvSurfaces::Compute( Standard_Real& u1,
0427 Standard_Real& v1,
0428 Standard_Real& u2,
0429 Standard_Real& v2,
0430 gp_Pnt& P,
0431 gp_Vec& Tg,
0432 gp_Vec2d& Tguv1,
0433 gp_Vec2d& Tguv2)
0434 {
0435 const IntSurf_Quadric& aQSurf = MyZerImpFunc.ISurface();
0436 const ThePSurface& aPSurf = MyZerImpFunc.PSurface();
0437 gp_Vec2d& aQuadTg = MyImplicitFirst ? Tguv1 : Tguv2;
0438 gp_Vec2d& aPrmTg = MyImplicitFirst ? Tguv2 : Tguv1;
0439
0440 //for square
0441 const Standard_Real aNullValue = Precision::Approximation()*
0442 Precision::Approximation(),
0443 anAngTol = Precision::Angular();
0444
0445 Standard_Real tu1=u1;
0446 Standard_Real tu2=u2;
0447 Standard_Real tv1=v1;
0448 Standard_Real tv2=v2;
0449
0450 if(MyHasBeenComputed) {
0451 if( (MyParOnS1.X()==u1)&&(MyParOnS1.Y()==v1)
0452 &&(MyParOnS2.X()==u2)&&(MyParOnS2.Y()==v2)) {
0453 return(MyIsTangent);
0454 }
0455 else if(MyHasBeenComputedbis == Standard_False) {
0456 MyTgbis = MyTg;
0457 MyTguv1bis = MyTguv1;
0458 MyTguv2bis = MyTguv2;
0459 MyPntbis = MyPnt;
0460 MyParOnS1bis = MyParOnS1;
0461 MyParOnS2bis = MyParOnS2;
0462 MyIsTangentbis = MyIsTangent;
0463 MyHasBeenComputedbis = MyHasBeenComputed;
0464 }
0465 }
0466
0467 if(MyHasBeenComputedbis) {
0468 if( (MyParOnS1bis.X()==u1)&&(MyParOnS1bis.Y()==v1)
0469 &&(MyParOnS2bis.X()==u2)&&(MyParOnS2bis.Y()==v2)) {
0470
0471 gp_Vec TV(MyTg);
0472 gp_Vec2d TV1(MyTguv1);
0473 gp_Vec2d TV2(MyTguv2);
0474 gp_Pnt TP(MyPnt);
0475 gp_Pnt2d TP1(MyParOnS1);
0476 gp_Pnt2d TP2(MyParOnS2);
0477 Standard_Boolean TB=MyIsTangent;
0478
0479 MyTg = MyTgbis;
0480 MyTguv1 = MyTguv1bis;
0481 MyTguv2 = MyTguv2bis;
0482 MyPnt = MyPntbis;
0483 MyParOnS1 = MyParOnS1bis;
0484 MyParOnS2 = MyParOnS2bis;
0485 MyIsTangent = MyIsTangentbis;
0486
0487 MyTgbis = TV;
0488 MyTguv1bis = TV1;
0489 MyTguv2bis = TV2;
0490 MyPntbis = TP;
0491 MyParOnS1bis = TP1;
0492 MyParOnS2bis = TP2;
0493 MyIsTangentbis = TB;
0494
0495 return(MyIsTangent);
0496 }
0497 }
0498
0499 math_Vector X(1,2);
0500 math_Vector BornInf(1,2), BornSup(1,2), Tolerance(1,2);
0501 //--- ThePSurfaceTool::GetResolution(aPSurf,Tolerance(1),Tolerance(2));
0502 Tolerance(1) = 1.0e-8; Tolerance(2) = 1.0e-8;
0503 Standard_Real binfu,bsupu,binfv,bsupv;
0504 binfu = ThePSurfaceTool::FirstUParameter(aPSurf);
0505 binfv = ThePSurfaceTool::FirstVParameter(aPSurf);
0506 bsupu = ThePSurfaceTool::LastUParameter(aPSurf);
0507 bsupv = ThePSurfaceTool::LastVParameter(aPSurf);
0508 BornInf(1) = binfu; BornSup(1) = bsupu;
0509 BornInf(2) = binfv; BornSup(2) = bsupv;
0510 Standard_Real TranslationU = 0., TranslationV = 0.;
0511
0512 if (!FillInitialVectorOfSolution(u1, v1, u2, v2,
0513 binfu, bsupu, binfv, bsupv,
0514 X,
0515 TranslationU, TranslationV))
0516 {
0517 MyIsTangent=MyIsTangentbis=Standard_False;
0518 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
0519 return(Standard_False);
0520 }
0521
0522 Standard_Boolean aRsnldIsDone = Standard_False;
0523 Standard_Real PourTesterU = X(1);
0524 Standard_Real PourTesterV = X(2);
0525 if (GetUseSolver())
0526 {
0527 math_FunctionSetRoot Rsnld(MyZerImpFunc);
0528 Rsnld.SetTolerance(Tolerance);
0529 Rsnld.Perform(MyZerImpFunc, X, BornInf, BornSup);
0530 aRsnldIsDone = Rsnld.IsDone();
0531 if (aRsnldIsDone)
0532 Rsnld.Root(X);
0533 }
0534 if(aRsnldIsDone || !GetUseSolver())
0535 {
0536 MyHasBeenComputed = Standard_True;
0537
0538 Standard_Real DistAvantApresU = Abs(PourTesterU-X(1));
0539 Standard_Real DistAvantApresV = Abs(PourTesterV-X(2));
0540
0541 MyPnt = P = ThePSurfaceTool::Value(aPSurf, X(1), X(2));
0542
0543 if( (DistAvantApresV <= 0.001 ) &&
0544 (DistAvantApresU <= 0.001 ))
0545 {
0546 gp_Vec aD1uPrm,aD1vPrm;
0547 gp_Vec aD1uQuad,aD1vQuad;
0548
0549 if(MyImplicitFirst)
0550 {
0551 u2 = X(1)-TranslationU;
0552 v2 = X(2)-TranslationV;
0553
0554 if(aQSurf.TypeQuadric() != GeomAbs_Plane)
0555 {
0556 while(u1-tu1>M_PI) u1-=M_PI+M_PI;
0557 while(tu1-u1>M_PI) u1+=M_PI+M_PI;
0558 }
0559
0560 MyParOnS1.SetCoord(tu1,tv1);
0561 MyParOnS2.SetCoord(tu2,tv2);
0562
0563 gp_Pnt aP2;
0564
0565 ThePSurfaceTool::D1(aPSurf, X(1), X(2), P, aD1uPrm, aD1vPrm);
0566 aQSurf.D1(u1,v1, aP2, aD1uQuad, aD1vQuad);
0567
0568 //Middle-point of P-P2 segment
0569 P.BaryCenter(1.0, aP2, 1.0);
0570 }
0571 else
0572 {
0573 u1 = X(1)-TranslationU;
0574 v1 = X(2)-TranslationV;
0575 //aQSurf.Parameters(P, u2, v2);
0576 if(aQSurf.TypeQuadric() != GeomAbs_Plane)
0577 {
0578 while(u2-tu2>M_PI) u2-=M_PI+M_PI;
0579 while(tu2-u2>M_PI) u2+=M_PI+M_PI;
0580 }
0581
0582 MyParOnS1.SetCoord(tu1,tv1);
0583 MyParOnS2.SetCoord(tu2,tu2);
0584
0585 gp_Pnt aP2;
0586 ThePSurfaceTool::D1(aPSurf, X(1), X(2), P, aD1uPrm, aD1vPrm);
0587
0588 aQSurf.D1(u2, v2, aP2, aD1uQuad, aD1vQuad);
0589
0590 //Middle-point of P-P2 segment
0591 P.BaryCenter(1.0, aP2, 1.0);
0592 }
0593
0594 MyPnt = P;
0595
0596 //Normals to the surfaces
0597 gp_Vec aNormalPrm(aD1uPrm.Crossed(aD1vPrm)),
0598 aNormalImp(aQSurf.Normale(MyPnt));
0599
0600 const Standard_Real aSQMagnPrm = aNormalPrm.SquareMagnitude(),
0601 aSQMagnImp = aNormalImp.SquareMagnitude();
0602
0603 Standard_Boolean isPrmSingular = Standard_False,
0604 isImpSingular = Standard_False;
0605
0606 if(IsSingular(aD1uPrm, aD1vPrm, aNullValue, anAngTol))
0607 {
0608 isPrmSingular = Standard_True;
0609 if(!SingularProcessing(aD1uPrm, aD1vPrm, Standard_True,
0610 aNullValue, anAngTol, Tg, aPrmTg))
0611 {
0612 MyIsTangent = Standard_False;
0613 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
0614 return Standard_False;
0615 }
0616
0617 MyTg = Tg;
0618 }
0619 else
0620 {
0621 aNormalPrm.Divide(sqrt(aSQMagnPrm));
0622 }
0623
0624 //Analogically for implicit surface
0625 if(aSQMagnImp < aNullValue)
0626 {
0627 isImpSingular = Standard_True;
0628
0629 if(!SingularProcessing(aD1uQuad, aD1vQuad, !isPrmSingular,
0630 aNullValue, anAngTol, Tg, aQuadTg))
0631 {
0632 MyIsTangent = Standard_False;
0633 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
0634 return Standard_False;
0635 }
0636
0637 MyTg = Tg;
0638 }
0639 else
0640 {
0641 aNormalImp.Divide(sqrt(aSQMagnImp));
0642 }
0643
0644 if(isImpSingular && isPrmSingular)
0645 {
0646 //All is OK. All abnormal cases were processed above.
0647
0648 MyTguv1 = Tguv1;
0649 MyTguv2 = Tguv2;
0650
0651 MyIsTangent=Standard_True;
0652 return MyIsTangent;
0653 }
0654 else if(!(isImpSingular || isPrmSingular))
0655 {
0656 //Processing pure non-singular case
0657 //(3D- and 2D-tangents are still not defined)
0658
0659 //Ask to pay attention to the fact that here
0660 //aNormalImp and aNormalPrm are normalized.
0661 //Therefore, @ \left \| \vec{Tg} \right \| = 0.0 @
0662 //if and only if (aNormalImp || aNormalPrm).
0663 Tg = aNormalImp.Crossed(aNormalPrm);
0664 }
0665
0666 const Standard_Real aSQMagnTg = Tg.SquareMagnitude();
0667
0668 if(aSQMagnTg < aNullValue)
0669 {
0670 MyIsTangent = Standard_False;
0671 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
0672 return Standard_False;
0673 }
0674
0675 //Normalize Tg vector
0676 Tg.Divide(sqrt(aSQMagnTg));
0677 MyTg = Tg;
0678
0679 if(!isPrmSingular)
0680 {
0681 //If isPrmSingular==TRUE then aPrmTg has already been computed.
0682
0683 if(!NonSingularProcessing(aD1uPrm, aD1vPrm, Tg, aNullValue, anAngTol, aPrmTg))
0684 {
0685 MyIsTangent = Standard_False;
0686 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
0687 return Standard_False;
0688 }
0689 }
0690
0691 if(!isImpSingular)
0692 {
0693 //If isImpSingular==TRUE then aQuadTg has already been computed.
0694
0695 if(!NonSingularProcessing(aD1uQuad, aD1vQuad, Tg, aNullValue, anAngTol, aQuadTg))
0696 {
0697 MyIsTangent = Standard_False;
0698 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
0699 return Standard_False;
0700 }
0701 }
0702
0703 MyTguv1 = Tguv1;
0704 MyTguv2 = Tguv2;
0705
0706 MyIsTangent=Standard_True;
0707
0708 #ifdef OCCT_DEBUG
0709 //cout << "+++++++++++++++++ ApproxInt_ImpPrmSvSurfaces::Compute(...) ++++++++++" << endl;
0710 //printf( "P2d_1(%+10.20f, %+10.20f); P2d_2(%+10.20f, %+10.20f)\n"
0711 // "P(%+10.20f, %+10.20f, %+10.20f);\n"
0712 // "Tg = {%+10.20f, %+10.20f, %+10.20f};\n"
0713 // "Tguv1 = {%+10.20f, %+10.20f};\n"
0714 // "Tguv2 = {%+10.20f, %+10.20f}\n",
0715 // u1, v1, u2, v2,
0716 // P.X(), P.Y(), P.Z(),
0717 // Tg.X(), Tg.Y(), Tg.Z(),
0718 // Tguv1.X(), Tguv1.Y(), Tguv2.X(), Tguv2.Y());
0719 //cout << "-----------------------------------------------------------------------" << endl;
0720 #endif
0721
0722 return Standard_True;
0723 }
0724 else {
0725 //-- cout<<" ApproxInt_ImpImpSvSurfaces.gxx : Distance apres recadrage Trop Grande "<<endl;
0726
0727 MyIsTangent=Standard_False;
0728 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
0729 return Standard_False;
0730 }
0731 }
0732 else {
0733 MyIsTangent = Standard_False;
0734 MyHasBeenComputed = MyHasBeenComputedbis = Standard_False;
0735 return Standard_False;
0736 }
0737 }
0738
0739 //=======================================================================
0740 //function : SeekPoint
0741 //purpose : Computes point on curve and
0742 // parameters on the surfaces.
0743 //=======================================================================
0744 Standard_Boolean ApproxInt_ImpPrmSvSurfaces::SeekPoint(const Standard_Real u1,
0745 const Standard_Real v1,
0746 const Standard_Real u2,
0747 const Standard_Real v2,
0748 IntSurf_PntOn2S& Point) {
0749 gp_Pnt aP;
0750 gp_Vec aT;
0751 gp_Vec2d aTS1,aTS2;
0752 const IntSurf_Quadric& aQSurf = MyZerImpFunc.ISurface();
0753 const ThePSurface& aPSurf = MyZerImpFunc.PSurface();
0754
0755 math_Vector X(1,2);
0756 math_Vector BornInf(1,2), BornSup(1,2), Tolerance(1,2);
0757 //--- ThePSurfaceTool::GetResolution(aPSurf,Tolerance(1),Tolerance(2));
0758 Tolerance(1) = 1.0e-8; Tolerance(2) = 1.0e-8;
0759 Standard_Real binfu,bsupu,binfv,bsupv;
0760 binfu = ThePSurfaceTool::FirstUParameter(aPSurf);
0761 binfv = ThePSurfaceTool::FirstVParameter(aPSurf);
0762 bsupu = ThePSurfaceTool::LastUParameter(aPSurf);
0763 bsupv = ThePSurfaceTool::LastVParameter(aPSurf);
0764 BornInf(1) = binfu; BornSup(1) = bsupu;
0765 BornInf(2) = binfv; BornSup(2) = bsupv;
0766 Standard_Real TranslationU = 0., TranslationV = 0.;
0767
0768 if (!FillInitialVectorOfSolution(u1, v1, u2, v2,
0769 binfu, bsupu, binfv, bsupv,
0770 X,
0771 TranslationU, TranslationV))
0772 return Standard_False;
0773
0774 Standard_Real NewU1, NewV1, NewU2, NewV2;
0775
0776 math_FunctionSetRoot Rsnld(MyZerImpFunc);
0777 Rsnld.SetTolerance(Tolerance);
0778 Rsnld.Perform(MyZerImpFunc,X,BornInf,BornSup);
0779 if(Rsnld.IsDone()) {
0780 MyHasBeenComputed = Standard_True;
0781 Rsnld.Root(X);
0782
0783 MyPnt = ThePSurfaceTool::Value(aPSurf, X(1), X(2));
0784
0785 if(MyImplicitFirst)
0786 {
0787 NewU2 = X(1)-TranslationU;
0788 NewV2 = X(2)-TranslationV;
0789
0790 aQSurf.Parameters(MyPnt, NewU1, NewV1);
0791 //adjust U
0792 if (aQSurf.TypeQuadric() != GeomAbs_Plane)
0793 {
0794 Standard_Real sign = (NewU1 > u1)? -1 : 1;
0795 while (Abs(u1 - NewU1) > M_PI)
0796 NewU1 += sign*(M_PI+M_PI);
0797 }
0798 }
0799 else
0800 {
0801 NewU1 = X(1)-TranslationU;
0802 NewV1 = X(2)-TranslationV;
0803
0804 aQSurf.Parameters(MyPnt, NewU2, NewV2);
0805 //adjust U
0806 if (aQSurf.TypeQuadric() != GeomAbs_Plane)
0807 {
0808 Standard_Real sign = (NewU2 > u2)? -1 : 1;
0809 while (Abs(u2 - NewU2) > M_PI)
0810 NewU2 += sign*(M_PI+M_PI);
0811 }
0812 }
0813 }
0814 else
0815 return Standard_False;
0816
0817 Point.SetValue(MyPnt, NewU1, NewV1, NewU2, NewV2);
0818 return Standard_True;
0819 }
0820 //--------------------------------------------------------------------------------
0821
0822 Standard_Boolean
0823 ApproxInt_ImpPrmSvSurfaces::FillInitialVectorOfSolution(const Standard_Real u1,
0824 const Standard_Real v1,
0825 const Standard_Real u2,
0826 const Standard_Real v2,
0827 const Standard_Real binfu,
0828 const Standard_Real bsupu,
0829 const Standard_Real binfv,
0830 const Standard_Real bsupv,
0831 math_Vector& X,
0832 Standard_Real& TranslationU,
0833 Standard_Real& TranslationV)
0834 {
0835 const ThePSurface& aPSurf = MyZerImpFunc.PSurface();
0836
0837 math_Vector F(1,1);
0838
0839 TranslationU = 0.0;
0840 TranslationV = 0.0;
0841
0842 if(MyImplicitFirst) {
0843 if(u2<binfu-0.0000000001) {
0844 if(ThePSurfaceTool::IsUPeriodic(aPSurf)) {
0845 Standard_Real d = ThePSurfaceTool::UPeriod(aPSurf);
0846 do { TranslationU+=d; } while(u2+TranslationU < binfu);
0847 }
0848 else
0849 return(Standard_False);
0850 }
0851 else if(u2>bsupu+0.0000000001) {
0852 if(ThePSurfaceTool::IsUPeriodic(aPSurf)) {
0853 Standard_Real d = ThePSurfaceTool::UPeriod(aPSurf);
0854 do { TranslationU-=d; } while(u2+TranslationU > bsupu);
0855 }
0856 else
0857 return(Standard_False);
0858 }
0859 if(v2<binfv-0.0000000001) {
0860 if(ThePSurfaceTool::IsVPeriodic(aPSurf)) {
0861 Standard_Real d = ThePSurfaceTool::VPeriod(aPSurf);
0862 do { TranslationV+=d; } while(v2+TranslationV < binfv);
0863 }
0864 else
0865 return(Standard_False);
0866 }
0867 else if(v2>bsupv+0.0000000001) {
0868 if(ThePSurfaceTool::IsVPeriodic(aPSurf)) {
0869 Standard_Real d = ThePSurfaceTool::VPeriod(aPSurf);
0870 do { TranslationV-=d; } while(v2+TranslationV > bsupv);
0871 }
0872 else
0873 return(Standard_False);
0874 }
0875 X(1) = u2+TranslationU;
0876 X(2) = v2+TranslationV;
0877 }
0878 else {
0879 if(u1<binfu-0.0000000001) {
0880 if(ThePSurfaceTool::IsUPeriodic(aPSurf)) {
0881 Standard_Real d = ThePSurfaceTool::UPeriod(aPSurf);
0882 do { TranslationU+=d; } while(u1+TranslationU < binfu);
0883 }
0884 else
0885 return(Standard_False);
0886 }
0887 else if(u1>bsupu+0.0000000001) {
0888 if(ThePSurfaceTool::IsUPeriodic(aPSurf)) {
0889 Standard_Real d = ThePSurfaceTool::UPeriod(aPSurf);
0890 do { TranslationU-=d; } while(u1+TranslationU > bsupu);
0891 }
0892 else
0893 return(Standard_False);
0894 }
0895 if(v1<binfv-0.0000000001) {
0896 if(ThePSurfaceTool::IsVPeriodic(aPSurf)) {
0897 Standard_Real d = ThePSurfaceTool::VPeriod(aPSurf);
0898 do { TranslationV+=d; } while(v1+TranslationV < binfv);
0899 }
0900 else
0901 return(Standard_False);
0902 }
0903 else if(v1>bsupv+0.0000000001) {
0904 if(ThePSurfaceTool::IsVPeriodic(aPSurf)) {
0905 Standard_Real d = ThePSurfaceTool::VPeriod(aPSurf);
0906 do { TranslationV-=d; } while(v1+TranslationV > bsupv);
0907 }
0908 else
0909 return(Standard_False);
0910 }
0911 X(1) = u1+TranslationU;
0912 X(2) = v1+TranslationV;
0913 }
0914
0915 //----------------------------------------------------
0916 //Make a small step from boundaries in order to avoid
0917 //finding "outboundaried" solution (Rsnld -> NotDone).
0918 if (GetUseSolver())
0919 {
0920 Standard_Real du = Max(Precision::Confusion(), ThePSurfaceTool::UResolution(aPSurf, Precision::Confusion()));
0921 Standard_Real dv = Max(Precision::Confusion(), ThePSurfaceTool::VResolution(aPSurf, Precision::Confusion()));
0922 if (X(1) - 0.0000000001 <= binfu) X(1) = X(1) + du;
0923 if (X(1) + 0.0000000001 >= bsupu) X(1) = X(1) - du;
0924 if (X(2) - 0.0000000001 <= binfv) X(2) = X(2) + dv;
0925 if (X(2) + 0.0000000001 >= bsupv) X(2) = X(2) - dv;
0926 }
0927
0928 return Standard_True;
0929 }
