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