src/GeomFill/GeomFill_Darboux.cxx

   1 // Created on: 1998-03-04
   2 // Created by: Roman BORISOV
   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
  18 #include <Adaptor2d_HCurve2d.hxx>
  19 #include <Adaptor3d_CurveOnSurface.hxx>
  20 #include <Adaptor3d_CurveOnSurfacePtr.hxx>
  21 #include <Adaptor3d_HCurve.hxx>
  22 #include <Adaptor3d_HCurveOnSurface.hxx>
  23 #include <Adaptor3d_HSurface.hxx>
  24 #include <CSLib.hxx>
  25 #include <Geom_UndefinedValue.hxx>
  26 #include <GeomAdaptor_HSurface.hxx>
  27 #include <GeomAdaptor_Surface.hxx>
  28 #include <GeomFill_Darboux.hxx>
  29 #include <GeomFill_TrihedronLaw.hxx>
  30 #include <gp_Pnt2d.hxx>
  31 #include <gp_Vec.hxx>
  32 #include <gp_Vec2d.hxx>
  33 #include <Standard_ConstructionError.hxx>
  34 #include <Standard_OutOfRange.hxx>
  35 #include <Standard_Type.hxx>
  36 #include <TColgp_Array2OfVec.hxx>
  37
  38 IMPLEMENT_STANDARD_RTTIEXT(GeomFill_Darboux,GeomFill_TrihedronLaw)
  39
  40 //=======================================================================
  41 //function : FDeriv
  42 //purpose  : computes (F/|F|)'
  43 //=======================================================================
  44 static gp_Vec FDeriv(const gp_Vec& F, const gp_Vec& DF)
  45 {
  46   Standard_Real Norma = F.Magnitude();
  47   gp_Vec Result = (DF - F*(F*DF)/(Norma*Norma))/Norma;
  48   return Result;
  49 }
  50
  51 //=======================================================================
  52 //function : DDeriv
  53 //purpose  : computes (F/|F|)''
  54 //=======================================================================
  55 static gp_Vec DDeriv(const gp_Vec& F, const gp_Vec& DF, const gp_Vec& D2F)
  56 {
  57   Standard_Real Norma = F.Magnitude();
  58   gp_Vec Result = (D2F - 2*DF*(F*DF)/(Norma*Norma))/Norma -
  59      F*((DF.SquareMagnitude() + F*D2F
  60         - 3*(F*DF)*(F*DF)/(Norma*Norma))/(Norma*Norma*Norma));
  61   return Result;
  62 }
  63
  64 //=======================================================================
  65 //function : NormalD0
  66 //purpose  : computes Normal to Surface
  67 //=======================================================================
  68 static void NormalD0(const Standard_Real U, const Standard_Real V, const Handle(Adaptor3d_HSurface)& Surf, gp_Dir& Normal, Standard_Integer& OrderU, Standard_Integer& OrderV)
  69 {
  70 //  gp_Vec D1U,D1V,D2U,D2V,DUV;
  71   gp_Vec D1U,D1V;
  72   GeomAbs_Shape Cont = (Surf->Surface().UContinuity() < Surf->Surface().VContinuity()) ?
  73     (Surf->Surface().UContinuity()) : (Surf->Surface().VContinuity());
  74   OrderU = OrderV = 0;
  75 #ifdef CHECK
  76   if (Cont == GeomAbs_C0) throw Geom_UndefinedValue();
  77 #endif
  78   gp_Pnt P;
  79   Surf->D1(U, V, P, D1U, D1V);
  80   Standard_Real MagTol=0.000000001;
  81   CSLib_NormalStatus NStatus;
  82   CSLib::Normal(D1U, D1V, MagTol, NStatus, Normal);
  83
  84   if (NStatus != CSLib_Defined) {
  85     if (Cont==GeomAbs_C0 ||
  86         Cont==GeomAbs_C1) {
  87       throw Geom_UndefinedValue();
  88       }
  89     Standard_Integer MaxOrder=3;
  90     TColgp_Array2OfVec DerNUV(0,MaxOrder,0,MaxOrder);
  91     TColgp_Array2OfVec DerSurf(0,MaxOrder+1,0,MaxOrder+1);
  92     Standard_Integer i,j;//OrderU,OrderV;
  93     Standard_Real Umin,Umax,Vmin,Vmax;
  94     Umin = Surf->Surface().FirstUParameter();
  95     Umax = Surf->Surface().LastUParameter();
  96     Vmin = Surf->Surface().FirstVParameter();
  97     Vmax = Surf->Surface().LastVParameter();
  98     for(i=1;i<=MaxOrder+1;i++){
  99       DerSurf.SetValue(i,0,Surf->DN(U,V,i,0));
 100     }
 101
 102     for(i=0;i<=MaxOrder+1;i++)
 103       for(j=1;j<=MaxOrder+1;j++){
 104         DerSurf.SetValue(i,j,Surf->DN(U,V,i,j));
 105       }
 106
 107     for(i=0;i<=MaxOrder;i++)
 108       for(j=0;j<=MaxOrder;j++){
 109         DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurf));
 110       }
 111
 112     CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,
 113                   NStatus,Normal,OrderU,OrderV);
 114
 115     if (NStatus != CSLib_Defined) {
 116 #ifdef OCCT_DEBUG
 117       std::cout << U << ", " << V<< std::endl;
 118       for(i=0;i<=MaxOrder;i++)
 119         for(j=0;j<=MaxOrder;j++){
 120           std::cout <<"("<<i <<"," << j << ") = "<< DerSurf(i,j).X() <<","
 121                << DerSurf(i,j).Y() <<"," << DerSurf(i,j).Z() << std::endl;
 122         }
 123 #endif
 124       throw Geom_UndefinedValue();
 125     }
 126   }
 127 }
 128
 129 //=======================================================================
 130 //function : NormalD1
 131 //purpose  : computes Normal to Surface and its first derivative
 132 //=======================================================================
 133 void NormalD1 (const Standard_Real U, const Standard_Real V,
 134                const Handle(Adaptor3d_HSurface)& Surf, gp_Dir& Normal,
 135                gp_Vec& D1UNormal, gp_Vec& D1VNormal)
 136 {
 137 #ifdef CHECK
 138   GeomAbs_Shape Cont = (Surf->Surface().UContinuity() < Surf->Surface().VContinuity()) ?
 139     (Surf->Surface().UContinuity()) : (Surf->Surface().VContinuity());
 140   if (Cont==GeomAbs_C0 ||
 141       Cont==GeomAbs_C1) {
 142     throw Geom_UndefinedDerivative();
 143   }
 144 #endif
 145   gp_Vec d2u, d2v, d2uv;
 146   gp_Pnt P;
 147   Surf->D2(U, V, P, D1UNormal, D1VNormal,  d2u, d2v, d2uv);
 148   Standard_Real MagTol=0.000000001;
 149   CSLib_NormalStatus NStatus;
 150   CSLib::Normal (D1UNormal, D1VNormal, MagTol, NStatus, Normal);
 151   Standard_Integer MaxOrder;
 152   if (NStatus == CSLib_Defined)
 153     MaxOrder=0;
 154   else
 155     MaxOrder=3;
 156   Standard_Integer OrderU,OrderV;
 157   TColgp_Array2OfVec DerNUV(0,MaxOrder+1,0,MaxOrder+1);
 158   TColgp_Array2OfVec DerSurf(0,MaxOrder+2,0,MaxOrder+2);
 159   Standard_Integer i,j;
 160   Standard_Real Umin,Umax,Vmin,Vmax;
 161   Umin = Surf->Surface().FirstUParameter();
 162   Umax = Surf->Surface().LastUParameter();
 163   Vmin = Surf->Surface().FirstVParameter();
 164   Vmax = Surf->Surface().LastVParameter();
 165
 166   DerSurf.SetValue(1, 0, D1UNormal);
 167     DerSurf.SetValue(0, 1, D1VNormal);
 168   DerSurf.SetValue(1, 1, d2uv);
 169   DerSurf.SetValue(2, 0, d2u);
 170   DerSurf.SetValue(0, 2, d2v);
 171   for(i=0;i<=MaxOrder+1;i++)
 172     for(j=i; j<=MaxOrder+2; j++ )
 173       if (i+j > 2) {
 174         DerSurf.SetValue(i,j,Surf->DN(U,V,i,j));
 175         if (i!=j) DerSurf.SetValue(j,i,Surf->DN(U,V,j,i));
 176       }
 177
 178   for(i=0;i<=MaxOrder+1;i++)
 179     for(j=0;j<=MaxOrder+1;j++){
 180       DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurf));
 181       }
 182
 183   CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,
 184                 NStatus,Normal,OrderU,OrderV);
 185   if (NStatus != CSLib_Defined) throw Geom_UndefinedValue();
 186
 187   D1UNormal = CSLib::DNNormal(1,0,DerNUV,OrderU,OrderV);
 188   D1VNormal = CSLib::DNNormal(0,1,DerNUV,OrderU,OrderV);
 189 }
 190
 191 //=======================================================================
 192 //function : NormalD2
 193 //purpose  : computes Normal to Surface and its first and second derivatives
 194 //=======================================================================
 195 void NormalD2 (const Standard_Real U, const Standard_Real V,
 196                const Handle(Adaptor3d_HSurface)& Surf, gp_Dir& Normal,
 197                gp_Vec& D1UNormal, gp_Vec& D1VNormal,
 198                gp_Vec& D2UNormal, gp_Vec& D2VNormal, gp_Vec& D2UVNormal)
 199 {
 200 #ifdef CHECK
 201   GeomAbs_Shape Cont = (Surf->Surface().UContinuity() < Surf->Surface().VContinuity()) ?
 202     (Surf->Surface().UContinuity()) : (Surf->Surface().VContinuity());
 203      if (Cont == GeomAbs_C0 || Continuity == GeomAbs_C1 ||
 204          Cont == GeomAbs_C2) { throw Geom_UndefinedDerivative(); }
 205 #endif
 206   gp_Vec d3u, d3uuv, d3uvv, d3v;
 207   gp_Pnt P;
 208   Surf->D3(U, V, P, D1UNormal, D1VNormal,
 209            D2UNormal, D2VNormal, D2UVNormal,
 210            d3u,d3v, d3uuv, d3uvv);
 211   Standard_Real MagTol=0.000000001;
 212   CSLib_NormalStatus NStatus;
 213   CSLib::Normal (D1UNormal, D1VNormal, MagTol, NStatus, Normal);
 214   Standard_Integer MaxOrder;
 215   if (NStatus == CSLib_Defined)
 216     MaxOrder=0;
 217   else
 218     MaxOrder=3;
 219   Standard_Integer OrderU,OrderV;
 220   TColgp_Array2OfVec DerNUV(0,MaxOrder+2,0,MaxOrder+2);
 221   TColgp_Array2OfVec DerSurf(0,MaxOrder+3,0,MaxOrder+3);
 222   Standard_Integer i,j;
 223
 224   Standard_Real Umin,Umax,Vmin,Vmax;
 225   Umin = Surf->Surface().FirstUParameter();
 226   Umax = Surf->Surface().LastUParameter();
 227   Vmin = Surf->Surface().FirstVParameter();
 228   Vmax = Surf->Surface().LastVParameter();
 229
 230   DerSurf.SetValue(1, 0, D1UNormal);
 231   DerSurf.SetValue(0, 1, D1VNormal);
 232   DerSurf.SetValue(1, 1, D2UVNormal);
 233   DerSurf.SetValue(2, 0, D2UNormal);
 234   DerSurf.SetValue(0, 2, D2VNormal);
 235   DerSurf.SetValue(3, 0, d3u);
 236   DerSurf.SetValue(2, 1, d3uuv);
 237   DerSurf.SetValue(1, 2, d3uvv);
 238   DerSurf.SetValue(0, 3, d3v);
 239   for(i=0;i<=MaxOrder+2;i++)
 240     for(j=i; j<=MaxOrder+3; j++ )
 241       if (i+j > 3) {
 242         DerSurf.SetValue(i,j,Surf->DN(U,V,i,j));
 243         if (i!=j) DerSurf.SetValue(j,i,Surf->DN(U,V,j,i));
 244       }
 245
 246   for(i=0;i<=MaxOrder+2;i++)
 247     for(j=0;j<=MaxOrder+2;j++){
 248       DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurf));
 249     }
 250
 251   CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,
 252                 NStatus,Normal,OrderU,OrderV);
 253   if (NStatus != CSLib_Defined) throw Geom_UndefinedValue();
 254
 255   D1UNormal = CSLib::DNNormal(1,0,DerNUV,OrderU,OrderV);
 256   D1VNormal = CSLib::DNNormal(0,1,DerNUV,OrderU,OrderV);
 257   D2UNormal = CSLib::DNNormal(2,0,DerNUV,OrderU,OrderV);
 258   D2VNormal = CSLib::DNNormal(0,2,DerNUV,OrderU,OrderV);
 259   D2UVNormal = CSLib::DNNormal(1,1,DerNUV,OrderU,OrderV);
 260 }
 261
 262 GeomFill_Darboux::GeomFill_Darboux()
 263 {
 264 }
 265
 266 Handle(GeomFill_TrihedronLaw) GeomFill_Darboux::Copy() const
 267 {
 268   Handle(GeomFill_Darboux) copy = new (GeomFill_Darboux)();
 269   if (!myCurve.IsNull()) copy->SetCurve(myCurve);
 270   return copy;
 271 }
 272
 273  Standard_Boolean GeomFill_Darboux::D0(const Standard_Real Param,gp_Vec& Tangent,gp_Vec& Normal,gp_Vec& BiNormal)
 274 {
 275   gp_Pnt2d C2d;
 276   gp_Vec2d D2d;
 277   gp_Pnt S;
 278   gp_Vec dS_du, dS_dv;
 279   Handle(Adaptor2d_HCurve2d) myCurve2d = Adaptor3d_CurveOnSurfacePtr(&(myTrimmed->Curve()))->GetCurve();
 280   Handle(Adaptor3d_HSurface) mySupport = Adaptor3d_CurveOnSurfacePtr(&(myTrimmed->Curve()))->GetSurface();
 281   Standard_Integer OrderU, OrderV;
 282   myCurve2d->D1(Param, C2d, D2d);
 283
 284 //  Normal = dS_du.Crossed(dS_dv).Normalized();
 285   gp_Dir NormalDir;
 286   NormalD0(C2d.X(), C2d.Y(), mySupport, NormalDir, OrderU, OrderV);
 287   BiNormal.SetXYZ(NormalDir.XYZ());
 288
 289   mySupport->D1(C2d.X(), C2d.Y(), S, dS_du, dS_dv);
 290 //  if(D2d.Magnitude() <= Precision::Confusion())
 291 //    DoTSingular(Param, Order);
 292
 293   Tangent = D2d.X()*dS_du + D2d.Y()*dS_dv;
 294   Tangent.Normalize();
 295
 296   Normal =  BiNormal;
 297   Normal ^= Tangent;
 298
 299
 300   return Standard_True;
 301 }
 302
 303  Standard_Boolean GeomFill_Darboux::D1(const Standard_Real Param,
 304                                        gp_Vec& Tangent,gp_Vec& DTangent,
 305                                        gp_Vec& Normal,gp_Vec& DNormal,
 306                                        gp_Vec& BiNormal,gp_Vec& DBiNormal)
 307 {
 308   gp_Pnt2d C2d;
 309   gp_Vec2d D2d, D2_2d;
 310   gp_Pnt S;
 311   gp_Vec dS_du, dS_dv, d2S_du, d2S_dv, d2S_duv, F, DF;
 312   Handle(Adaptor2d_HCurve2d) myCurve2d = Adaptor3d_CurveOnSurfacePtr(&(myTrimmed->Curve()))->GetCurve();
 313   Handle(Adaptor3d_HSurface) mySupport = Adaptor3d_CurveOnSurfacePtr(&(myTrimmed->Curve()))->GetSurface();
 314 //  Standard_Integer Order;
 315   myCurve2d->D2(Param, C2d, D2d, D2_2d);
 316   mySupport->D2(C2d.X(), C2d.Y(), S, dS_du, dS_dv,
 317                 d2S_du, d2S_dv, d2S_duv);
 318 //  if(D2d.Magnitude() <= Precision::Confusion())
 319 //    DoTSingular(Param, Order);
 320
 321   F = D2d.X()*dS_du + D2d.Y()*dS_dv;
 322   Tangent = F.Normalized();
 323   DF = D2_2d.X()*dS_du + D2_2d.Y()*dS_dv +
 324        d2S_du*D2d.X()*D2d.X() + 2*d2S_duv*D2d.X()*D2d.Y() +
 325        d2S_dv*D2d.Y()*D2d.Y();
 326
 327   DTangent = FDeriv(F, DF);
 328
 329   gp_Dir NormalDir;
 330   gp_Vec D1UNormal, D1VNormal;
 331   NormalD1(C2d.X(), C2d.Y(), mySupport, NormalDir, D1UNormal, D1VNormal);
 332   BiNormal.SetXYZ(NormalDir.XYZ());
 333   DBiNormal = D1UNormal*D2d.X() + D1VNormal*D2d.Y();
 334
 335   Normal =  BiNormal;
 336   Normal ^= Tangent;
 337   DNormal = BiNormal.Crossed(DTangent) + DBiNormal.Crossed(Tangent);
 338
 339   return Standard_True;
 340 }
 341
 342  Standard_Boolean GeomFill_Darboux::D2(const Standard_Real Param,
 343                                        gp_Vec& Tangent,gp_Vec& DTangent,gp_Vec& D2Tangent,
 344                                        gp_Vec& Normal,gp_Vec& DNormal,gp_Vec& D2Normal,
 345                                        gp_Vec& BiNormal,gp_Vec& DBiNormal,gp_Vec& D2BiNormal)
 346 {
 347   gp_Pnt2d C2d;
 348   gp_Vec2d D2d, D2_2d, D3_2d;
 349   gp_Pnt S;
 350   gp_Vec dS_du, dS_dv, d2S_du, d2S_dv, d2S_duv,
 351          d3S_du, d3S_dv, d3S_duuv, d3S_duvv, F, DF, D2F;
 352   Handle(Adaptor2d_HCurve2d) myCurve2d = Adaptor3d_CurveOnSurfacePtr(&(myTrimmed->Curve()))->GetCurve();
 353   Handle(Adaptor3d_HSurface) mySupport = Adaptor3d_CurveOnSurfacePtr(&(myTrimmed->Curve()))->GetSurface();
 354 //  Standard_Integer Order;
 355   myCurve2d->D3(Param, C2d, D2d, D2_2d, D3_2d);
 356   mySupport->D3(C2d.X(), C2d.Y(), S, dS_du, dS_dv,
 357                 d2S_du, d2S_dv, d2S_duv, d3S_du, d3S_dv, d3S_duuv, d3S_duvv);
 358 //  if(D2d.Magnitude() <= Precision::Confusion())
 359 //    DoTSingular(Param, Order);
 360
 361   F = D2d.X()*dS_du + D2d.Y()*dS_dv;
 362   Tangent = F.Normalized();
 363   DF = D2_2d.X()*dS_du + D2_2d.Y()*dS_dv +
 364        d2S_du*D2d.X()*D2d.X() + 2*d2S_duv*D2d.X()*D2d.Y() +
 365        d2S_dv*D2d.Y()*D2d.Y();
 366
 367   D2F = D3_2d.X()*dS_du + D3_2d.Y()*dS_dv +
 368         D2_2d.X()*(D2d.X()*d2S_du + D2d.Y()*d2S_duv) +
 369         D2_2d.Y()*(D2d.X()*d2S_duv + D2d.Y()*d2S_dv) +
 370         D2d.X()*D2d.X()*(D2d.X()*d3S_du + D2d.Y()*d3S_duuv) +
 371         D2d.Y()*D2d.Y()*(D2d.X()*d3S_duvv + D2d.Y()*d3S_dv) +
 372         2*(
 373            D2d.X()*D2_2d.X()*d2S_du + D2d.Y()*D2_2d.Y()*d2S_dv +
 374            D2d.X()*D2d.Y()*(D2d.X()*d3S_duuv + D2d.Y()*d3S_duvv) +
 375            d2S_duv*(D2_2d.X()*D2d.Y() + D2d.X()*D2_2d.Y())
 376            );
 377
 378   DTangent = FDeriv(F, DF);
 379   D2Tangent = DDeriv(F, DF, D2F);
 380
 381   gp_Dir NormalDir;
 382   gp_Vec D1UNormal, D1VNormal, D2UNormal, D2VNormal, D2UVNormal;
 383   NormalD2(C2d.X(), C2d.Y(), mySupport, NormalDir, D1UNormal, D1VNormal,
 384            D2UNormal, D2VNormal, D2UVNormal);
 385   BiNormal.SetXYZ(NormalDir.XYZ());
 386   DBiNormal = D1UNormal*D2d.X() + D1VNormal*D2d.Y();
 387   D2BiNormal = D1UNormal*D2_2d.X() + D1VNormal*D2_2d.Y() + D2UNormal*D2d.X()*D2d.X() +
 388     2*D2UVNormal*D2d.X()*D2d.Y() + D2VNormal*D2d.Y()*D2d.Y();
 389
 390   Normal =  BiNormal;
 391   Normal ^= Tangent;
 392   DNormal = BiNormal.Crossed(DTangent) + DBiNormal.Crossed(Tangent);
 393   D2Normal = BiNormal.Crossed(D2Tangent) + 2*DBiNormal.Crossed(DTangent) + D2BiNormal.Crossed(Tangent);
 394
 395   return Standard_True;
 396 }
 397
 398  Standard_Integer GeomFill_Darboux::NbIntervals(const GeomAbs_Shape S) const
 399 {
 400   return myCurve->NbIntervals(S);
 401 }
 402
 403  void GeomFill_Darboux::Intervals(TColStd_Array1OfReal& T,const GeomAbs_Shape S) const
 404 {
 405   myCurve->Intervals(T, S);
 406 }
 407
 408  void GeomFill_Darboux::GetAverageLaw(gp_Vec& ATangent,gp_Vec& ANormal,gp_Vec& ABiNormal)
 409 {
 410   Standard_Integer Num = 20; //order of digitalization
 411   gp_Vec T, N, BN;
 412   ATangent = gp_Vec(0, 0, 0);
 413   ANormal = gp_Vec(0, 0, 0);
 414   ABiNormal = gp_Vec(0, 0, 0);
 415   Standard_Real Step = (myTrimmed->LastParameter() -
 416                         myTrimmed->FirstParameter()) / Num;
 417   Standard_Real Param;
 418   for (Standard_Integer i = 0; i <= Num; i++) {
 419     Param = myTrimmed->FirstParameter() + i*Step;
 420     if (Param > myTrimmed->LastParameter()) Param = myTrimmed->LastParameter();
 421     D0(Param, T, N, BN);
 422     ATangent += T;
 423     ANormal += N;
 424     ABiNormal += BN;
 425   }
 426   ATangent /= Num + 1;
 427   ANormal /= Num + 1;
 428
 429   ATangent.Normalize();
 430   ABiNormal = ATangent.Crossed(ANormal).Normalized();
 431   ANormal = ABiNormal.Crossed(ATangent);
 432 }
 433
 434  Standard_Boolean GeomFill_Darboux::IsConstant() const
 435 {
 436   return (myCurve->GetType() == GeomAbs_Line);
 437 }
 438
 439  Standard_Boolean GeomFill_Darboux::IsOnlyBy3dCurve() const
 440 {
 441   return Standard_False;
 442 }