src/BRepBlend/BRepBlend_AppFuncRoot.cxx

   1 // Created on: 1998-05-12
   2 // Created by: Philippe NOUAILLE
   3 // Copyright (c) 1998-1999 Matra Datavision
   4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
   5 //
   6 // The content of this file is subject to the Open CASCADE Technology Public
   7 // License Version 6.5 (the "License"). You may not use the content of this file
   8 // except in compliance with the License. Please obtain a copy of the License
   9 // at http://www.opencascade.org and read it completely before using this file.
  10 //
  11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  13 //
  14 // The Original Code and all software distributed under the License is
  15 // distributed on an "AS IS" basis, without warranty of any kind, and the
  16 // Initial Developer hereby disclaims all such warranties, including without
  17 // limitation, any warranties of merchantability, fitness for a particular
  18 // purpose or non-infringement. Please see the License for the specific terms
  19 // and conditions governing the rights and limitations under the License.
  20
  21
  22
  23 #include <BRepBlend_AppFuncRoot.ixx>
  24 #include <Blend_AppFunction.hxx>
  25
  26 #include <Blend_Point.hxx>
  27 #include <BRepBlend_Line.hxx>
  28
  29 #include <math_FunctionSetRoot.hxx>
  30
  31 #include <TColgp_HArray1OfPnt.hxx>
  32 #include <TColgp_HArray1OfPnt2d.hxx>
  33 #include <TColgp_HArray1OfVec.hxx>
  34 #include <TColgp_HArray1OfVec2d.hxx>
  35
  36 #include <TColStd_HArray1OfReal.hxx>
  37 #include <TColStd_HArray1OfInteger.hxx>
  38
  39 BRepBlend_AppFuncRoot::BRepBlend_AppFuncRoot(Handle(BRepBlend_Line)& Line,
  40                                              Blend_AppFunction&      Func,
  41                                              const Standard_Real     Tol3d,
  42                                              const Standard_Real     Tol2d)
  43 :myLine(Line),
  44  myFunc(&Func),
  45  myTolerance(1,Func.NbVariables()),
  46  X1(1,Func.NbVariables()),
  47  X2(1,Func.NbVariables()),
  48  XInit(1,Func.NbVariables()),
  49  Sol(1,Func.NbVariables())
  50 {
  51   Standard_Integer NbPoles, NbKnots, Degree, NbPoles2d;
  52   Standard_Integer ii;
  53
  54   //  Tolerances
  55   Func.GetTolerance(myTolerance, Tol3d);
  56   Standard_Integer dim = Func.NbVariables();
  57   for (ii=1; ii<= dim; ii++) {
  58     if (myTolerance(ii)>Tol2d) { myTolerance(ii) = Tol2d;}
  59   }
  60
  61   //  Tables
  62   Func.GetShape( NbPoles, NbKnots, Degree, NbPoles2d);
  63
  64   // Calculation of BaryCentre (rationnal case).
  65   if (Func.IsRational()) {
  66     Standard_Real Xmax =-1.e100, Xmin = 1.e100,
  67     Ymax =-1.e100, Ymin = 1.e100,
  68     Zmax =-1.e100, Zmin = 1.e100;
  69     Blend_Point P;
  70     for (ii=1; ii<=myLine->NbPoints(); ii++) {
  71       P = myLine->Point(ii);
  72       Xmax = Max ( Max(P.PointOnS1().X(), P.PointOnS2().X()), Xmax);
  73       Xmin = Min ( Min(P.PointOnS1().X(), P.PointOnS2().X()), Xmin);
  74       Ymax = Max ( Max(P.PointOnS1().Y(), P.PointOnS2().Y()), Ymax);
  75       Ymin = Min ( Min(P.PointOnS1().Y(), P.PointOnS2().Y()), Ymin);
  76       Zmax = Max ( Max(P.PointOnS1().Z(), P.PointOnS2().Z()), Zmax);
  77       Zmin = Min ( Min(P.PointOnS1().Z(), P.PointOnS2().Z()), Zmin);
  78
  79       myBary.SetCoord((Xmax+Xmin)/2, (Ymax+Ymin)/2, (Zmax+Zmin)/2);
  80     }
  81   }
  82   else {myBary.SetCoord(0,0,0);}
  83 }
  84
  85 //================================================================================
  86 // Function: D0
  87 // Purpose : Calculation of section for v = Param, if calculation fails
  88 //           Standard_False is raised.
  89 //================================================================================
  90 Standard_Boolean BRepBlend_AppFuncRoot::D0(const Standard_Real Param,
  91                                            const Standard_Real First,
  92                                            const Standard_Real Last,TColgp_Array1OfPnt& Poles,
  93                                            TColgp_Array1OfPnt2d& Poles2d,TColStd_Array1OfReal& Weigths)
  94 {
  95   Standard_Boolean Ok=Standard_True;
  96   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
  97   Ok = SearchPoint( *Func, Param, myPnt);
  98
  99   if (Ok) (*Func).Section(myPnt,
 100                           Poles,
 101                           Poles2d,
 102                           Weigths);
 103   return Ok;
 104 }
 105
 106 //================================================================================
 107 // Function: D1
 108 // Purpose : Calculation of the partial derivative of the section corresponding to v
 109 //           for v = Param, if the calculation fails Standard_False is raised.
 110 //================================================================================
 111 Standard_Boolean BRepBlend_AppFuncRoot::D1(const Standard_Real Param,
 112                                            const Standard_Real First,
 113                                            const Standard_Real Last,
 114                                            TColgp_Array1OfPnt& Poles,
 115                                            TColgp_Array1OfVec& DPoles,
 116                                            TColgp_Array1OfPnt2d& Poles2d,
 117                                            TColgp_Array1OfVec2d& DPoles2d,
 118                                            TColStd_Array1OfReal& Weigths,
 119                                            TColStd_Array1OfReal& DWeigths)
 120 {
 121   Standard_Boolean Ok=Standard_True;
 122   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 123
 124   Ok = SearchPoint( *Func, Param, myPnt);
 125
 126   if (Ok) {
 127     Ok = (*Func).Section(myPnt,
 128                          Poles, DPoles,
 129                          Poles2d, DPoles2d,
 130                          Weigths, DWeigths);
 131   }
 132
 133   return Ok;
 134 }
 135
 136 //===========================================================================
 137 // Function: D2
 138 // Purpose : Calculation of the derivative and second partial of the
 139 //           section corresponding to v.
 140 //           For v = Param, if the calculation fails Standard_False is raised.
 141 //===========================================================================
 142 Standard_Boolean BRepBlend_AppFuncRoot::D2(const Standard_Real Param,
 143                                            const Standard_Real First,
 144                                            const Standard_Real Last,
 145                                            TColgp_Array1OfPnt& Poles,
 146                                            TColgp_Array1OfVec& DPoles,
 147                                            TColgp_Array1OfVec& D2Poles,
 148                                            TColgp_Array1OfPnt2d& Poles2d,
 149                                            TColgp_Array1OfVec2d& DPoles2d,
 150                                            TColgp_Array1OfVec2d& D2Poles2d,
 151                                            TColStd_Array1OfReal& Weigths,
 152                                            TColStd_Array1OfReal& DWeigths,
 153                                            TColStd_Array1OfReal& D2Weigths)
 154 {
 155   Standard_Boolean Ok=Standard_True;
 156   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 157
 158   Ok = SearchPoint( *Func, Param, myPnt);
 159   if (Ok) {
 160     Ok = (*Func).Section(myPnt,
 161                          Poles, DPoles, D2Poles,
 162                          Poles2d, DPoles2d, D2Poles2d,
 163                          Weigths, DWeigths, D2Weigths);
 164   }
 165   return Ok;
 166 }
 167
 168 Standard_Integer BRepBlend_AppFuncRoot::Nb2dCurves() const
 169 {
 170   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 171   Standard_Integer i,j,k,nbpol2d;
 172   (*Func).GetShape(i,j,k,nbpol2d);
 173   return nbpol2d;
 174 }
 175
 176 void BRepBlend_AppFuncRoot::SectionShape(Standard_Integer& NbPoles,
 177                                          Standard_Integer& NbKnots,
 178                                          Standard_Integer& Degree) const
 179 {
 180   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 181   Standard_Integer ii;
 182   (*Func).GetShape( NbPoles, NbKnots, Degree, ii);
 183 }
 184
 185 void BRepBlend_AppFuncRoot::Knots(TColStd_Array1OfReal& TKnots) const
 186 {
 187   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 188   Func->Knots(TKnots);
 189 }
 190
 191 void BRepBlend_AppFuncRoot::Mults(TColStd_Array1OfInteger& TMults) const
 192 {
 193   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 194   Func->Mults(TMults);
 195 }
 196
 197 Standard_Boolean BRepBlend_AppFuncRoot::IsRational() const
 198 {
 199   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 200   return  (*Func).IsRational();
 201 }
 202
 203 Standard_Integer BRepBlend_AppFuncRoot::NbIntervals(const GeomAbs_Shape S) const
 204 {
 205   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 206   return  Func->NbIntervals(S);
 207 }
 208
 209 void BRepBlend_AppFuncRoot::Intervals(TColStd_Array1OfReal& T,const GeomAbs_Shape S) const
 210 {
 211   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 212   Func->Intervals(T, S);
 213 }
 214
 215 void BRepBlend_AppFuncRoot::SetInterval(const Standard_Real First,const Standard_Real Last)
 216 {
 217   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 218   Func->Set(First, Last);
 219 }
 220
 221 void BRepBlend_AppFuncRoot::Resolution(const Standard_Integer Index,
 222                                        const Standard_Real Tol,
 223                                        Standard_Real& TolU,
 224                                        Standard_Real& TolV) const
 225 {
 226   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 227   Func->Resolution(Index,Tol,TolU,TolV);
 228 }
 229
 230 void BRepBlend_AppFuncRoot::GetTolerance(const Standard_Real BoundTol,
 231                                          const Standard_Real SurfTol,
 232                                          const Standard_Real AngleTol,
 233                                          TColStd_Array1OfReal& Tol3d) const
 234 {
 235   Standard_Integer ii;
 236   math_Vector V3d(1, Tol3d.Length()), V1d(1, Tol3d.Length());
 237   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 238
 239   Func->GetTolerance(BoundTol, SurfTol, AngleTol, V3d, V1d);
 240   for (ii=1; ii<=Tol3d.Length(); ii++) Tol3d(ii) = V3d(ii);
 241 }
 242
 243 void BRepBlend_AppFuncRoot::SetTolerance(const Standard_Real Tol3d,
 244                                          const Standard_Real Tol2d)
 245 {
 246   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 247   Standard_Integer ii, dim = Func->NbVariables();
 248   Func->GetTolerance(myTolerance, Tol3d);
 249   for (ii=1; ii<=dim; ii++) {
 250     if (myTolerance(ii)>Tol2d) { myTolerance(ii) = Tol2d;}
 251   }
 252 }
 253
 254 gp_Pnt BRepBlend_AppFuncRoot::BarycentreOfSurf() const
 255 {
 256   return myBary;
 257 }
 258
 259 Standard_Real BRepBlend_AppFuncRoot::MaximalSection() const
 260 {
 261   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 262   return Func->GetSectionSize();
 263 }
 264
 265 void BRepBlend_AppFuncRoot::GetMinimalWeight(TColStd_Array1OfReal& Weigths) const
 266 {
 267   Blend_AppFunction* Func = (Blend_AppFunction*)myFunc;
 268   Func->GetMinimalWeight(Weigths);
 269 }
 270
 271
 272 //================================================================================
 273 //
 274 // Function : SearchPoint
 275 //
 276 // Purpose : Find point solution with parameter Param (on 2 Surfaces)
 277 //
 278 // Algorithm :
 279 //     1) Approximative solution is found from already calculated Points
 280 //     2) Convergence is done by a method of type Newton
 281 //
 282 // Possible causes of fails :
 283 //        - Singularity on surfaces.
 284 //        - no information oin the "line" resulting from processing.
 285 //
 286 //================================================================================
 287
 288 Standard_Boolean BRepBlend_AppFuncRoot::SearchPoint(Blend_AppFunction& Func,
 289                                                     const Standard_Real Param,
 290                                                     Blend_Point& Pnt)
 291 {
 292   Standard_Boolean Trouve;
 293   Standard_Integer dim = Func.NbVariables();
 294   // (1) Find a point of init
 295   Standard_Integer I1=1, I2=myLine->NbPoints(), Index;
 296   Standard_Real t1, t2;
 297
 298   //  (1.a) It is checked if it is inside
 299   if (Param < myLine->Point(I1).Parameter()) {return Standard_False;}
 300   if (Param > myLine->Point(I2).Parameter()) {return Standard_False;}
 301
 302   //  (1.b) Find the interval
 303   Trouve = SearchLocation(Param, I1, I2, Index);
 304
 305   //  (1.c) If the point is already calculated it is returned
 306   if (Trouve) {
 307     Pnt = myLine->Point(Index);
 308     Vec(XInit,Pnt);
 309   }
 310   else {
 311     //  (1.d) Intialisation by linear interpolation
 312     Pnt = myLine->Point(Index);
 313     Vec(X1,Pnt);
 314     t1 = Pnt.Parameter();
 315
 316     Pnt = myLine->Point(Index+1);
 317     Vec(X2,Pnt);
 318     t2 = Pnt.Parameter();
 319
 320     Standard_Real Parammt1 = (Param-t1) / (t2-t1);
 321     Standard_Real t2mParam = (t2-Param) / (t2-t1);
 322     for(Standard_Integer i = 1; i <= dim; i++){
 323       XInit(i) = X2(i) * Parammt1 + X1(i) * t2mParam;
 324     }
 325   }
 326
 327   // (2) Calculation of the solution ------------------------
 328   Func.Set(Param);
 329   Func.GetBounds(X1, X2);
 330   math_FunctionSetRoot rsnld(Func, myTolerance, 30);
 331
 332   rsnld.Perform(Func, XInit, X1, X2);
 333
 334   if (!rsnld.IsDone()) {
 335 # ifdef DEB
 336     cout << "AppFunc : RNLD Not done en t = " <<  Param  << endl;
 337 # endif
 338     return Standard_False;
 339   }
 340   rsnld.Root(Sol);
 341
 342   // (3) Storage of the point
 343   Point(Func,Param,Sol,Pnt);
 344
 345   // (4) Insertion of the point if the calculation seems long.
 346   if ((!Trouve)&&(rsnld.NbIterations()>3)) {
 347 #ifdef DEB
 348     cout << "Evaluation in t = " <<  Param << "given" << endl;
 349     rsnld.Dump(cout);
 350 #endif
 351     myLine->InsertBefore(Index+1, Pnt);
 352   }
 353   return Standard_True;
 354 }
 355
 356
 357 //=============================================================================
 358 //
 359 // Function : SearchLocation
 360 //
 361 // Purpose : Binary search of the line of the parametric interval containing
 362 //           Param in the list of calculated points (myline)
 363 //           if the point of parameter Param is already stored in the list
 364 //           True is raised and ParamIndex corresponds to line of Point.
 365 //           Complexity of this algorithm is log(n)/log(2)
 366 //================================================================================
 367 Standard_Boolean BRepBlend_AppFuncRoot::SearchLocation(const Standard_Real Param,
 368                                                        const Standard_Integer FirstIndex,
 369                                                        const Standard_Integer LastIndex,
 370                                                        Standard_Integer& ParamIndex) const
 371 {
 372   Standard_Integer Ideb = FirstIndex, Ifin =  LastIndex, Idemi;
 373   Standard_Real Valeur;
 374
 375   Valeur = myLine->Point(Ideb).Parameter();
 376   if (Param == Valeur) {
 377     ParamIndex = Ideb;
 378     return Standard_True;
 379   }
 380
 381   Valeur = myLine->Point(Ifin).Parameter();
 382   if (Param == Valeur) {
 383     ParamIndex = Ifin;
 384     return Standard_True;
 385   }
 386
 387   while ( Ideb+1 != Ifin) {
 388     Idemi = (Ideb+Ifin)/2;
 389     Valeur = myLine->Point(Idemi).Parameter();
 390     if (Valeur < Param) {Ideb = Idemi;}
 391     else {
 392       if ( Valeur > Param) { Ifin = Idemi;}
 393       else                 { ParamIndex = Idemi;
 394                              return Standard_True;}
 395     }
 396   }
 397
 398   ParamIndex = Ideb;
 399   return Standard_False;
 400 }
 401