[occt.git] / src / GeomFill / GeomFill_Darboux.cxx

// Created on: 1998-03-04
// Created by: Roman BORISOV
// Copyright (c) 1998-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.


#include <Adaptor2d_HCurve2d.hxx>
#include <Adaptor3d_CurveOnSurface.hxx>
#include <Adaptor3d_CurveOnSurfacePtr.hxx>
#include <Adaptor3d_HCurve.hxx>
#include <Adaptor3d_HCurveOnSurface.hxx>
#include <Adaptor3d_HSurface.hxx>
#include <CSLib.hxx>
#include <Geom_UndefinedValue.hxx>
#include <GeomAdaptor_HSurface.hxx>
#include <GeomAdaptor_Surface.hxx>
#include <GeomFill_Darboux.hxx>
#include <GeomFill_TrihedronLaw.hxx>
#include <gp_Pnt2d.hxx>
#include <gp_Vec.hxx>
#include <gp_Vec2d.hxx>
#include <Standard_ConstructionError.hxx>
#include <Standard_OutOfRange.hxx>
#include <Standard_Type.hxx>
#include <TColgp_Array2OfVec.hxx>

IMPLEMENT_STANDARD_RTTIEXT(GeomFill_Darboux,GeomFill_TrihedronLaw)

//=======================================================================
//function : FDeriv
//purpose  : computes (F/|F|)'
//=======================================================================
static gp_Vec FDeriv(const gp_Vec& F, const gp_Vec& DF)
{
  Standard_Real Norma = F.Magnitude();
  gp_Vec Result = (DF - F*(F*DF)/(Norma*Norma))/Norma;
  return Result;
}

//=======================================================================
//function : DDeriv
//purpose  : computes (F/|F|)''
//=======================================================================
static gp_Vec DDeriv(const gp_Vec& F, const gp_Vec& DF, const gp_Vec& D2F)
{
  Standard_Real Norma = F.Magnitude();
  gp_Vec Result = (D2F - 2*DF*(F*DF)/(Norma*Norma))/Norma - 
     F*((DF.SquareMagnitude() + F*D2F 
        - 3*(F*DF)*(F*DF)/(Norma*Norma))/(Norma*Norma*Norma));
  return Result;
}

//=======================================================================
//function : NormalD0
//purpose  : computes Normal to Surface
//=======================================================================
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)
{
//  gp_Vec D1U,D1V,D2U,D2V,DUV;
  gp_Vec D1U,D1V;
  GeomAbs_Shape Cont = (Surf->Surface().UContinuity() < Surf->Surface().VContinuity()) ? 
    (Surf->Surface().UContinuity()) : (Surf->Surface().VContinuity());
  OrderU = OrderV = 0;
#ifdef CHECK  
  if (Cont == GeomAbs_C0) Geom_UndefinedValue::Raise();
#endif
  gp_Pnt P;
  Surf->D1(U, V, P, D1U, D1V);
  Standard_Real MagTol=0.000000001;
  CSLib_NormalStatus NStatus;
  CSLib::Normal(D1U, D1V, MagTol, NStatus, Normal);

  if (NStatus != CSLib_Defined) {
    if (Cont==GeomAbs_C0 || 
        Cont==GeomAbs_C1) {
      Geom_UndefinedValue::Raise();
      }
    Standard_Integer MaxOrder=3;
    TColgp_Array2OfVec DerNUV(0,MaxOrder,0,MaxOrder);
    TColgp_Array2OfVec DerSurf(0,MaxOrder+1,0,MaxOrder+1);
    Standard_Integer i,j;//OrderU,OrderV;
    Standard_Real Umin,Umax,Vmin,Vmax;
    Umin = Surf->Surface().FirstUParameter();
    Umax = Surf->Surface().LastUParameter();
    Vmin = Surf->Surface().FirstVParameter();
    Vmax = Surf->Surface().LastVParameter();
    for(i=1;i<=MaxOrder+1;i++){
      DerSurf.SetValue(i,0,Surf->DN(U,V,i,0));
    }
    
    for(i=0;i<=MaxOrder+1;i++)
      for(j=1;j<=MaxOrder+1;j++){
        DerSurf.SetValue(i,j,Surf->DN(U,V,i,j));
      }    
    
    for(i=0;i<=MaxOrder;i++)
      for(j=0;j<=MaxOrder;j++){
        DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurf));
      }
    
    CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,
                  NStatus,Normal,OrderU,OrderV);
    
    if (NStatus != CSLib_Defined) {
#ifdef OCCT_DEBUG
      cout << U << ", " << V<< endl;
      for(i=0;i<=MaxOrder;i++)
	for(j=0;j<=MaxOrder;j++){
	  cout <<"("<<i <<"," << j << ") = "<< DerSurf(i,j).X() <<","
	       << DerSurf(i,j).Y() <<"," << DerSurf(i,j).Z() << endl;
	}
#endif
      Geom_UndefinedValue::Raise();
    }
  }
}

//=======================================================================
//function : NormalD1
//purpose  : computes Normal to Surface and its first derivative
//=======================================================================
void NormalD1 (const Standard_Real U, const Standard_Real V, 
               const Handle(Adaptor3d_HSurface)& Surf, gp_Dir& Normal, 
               gp_Vec& D1UNormal, gp_Vec& D1VNormal)
{                                    
#ifdef CHECK  
  GeomAbs_Shape Cont = (Surf->Surface().UContinuity() < Surf->Surface().VContinuity()) ? 
    (Surf->Surface().UContinuity()) : (Surf->Surface().VContinuity());
  if (Cont==GeomAbs_C0 || 
      Cont==GeomAbs_C1) { 
    Geom_UndefinedDerivative::Raise();
  }
#endif
  gp_Vec d2u, d2v, d2uv;
  gp_Pnt P;
  Surf->D2(U, V, P, D1UNormal, D1VNormal,  d2u, d2v, d2uv);
  Standard_Real MagTol=0.000000001;
  CSLib_NormalStatus NStatus;
  CSLib::Normal (D1UNormal, D1VNormal, MagTol, NStatus, Normal);
  Standard_Integer MaxOrder;
  if (NStatus == CSLib_Defined) 
    MaxOrder=0;
  else 
    MaxOrder=3;
  Standard_Integer OrderU,OrderV;
  TColgp_Array2OfVec DerNUV(0,MaxOrder+1,0,MaxOrder+1);
  TColgp_Array2OfVec DerSurf(0,MaxOrder+2,0,MaxOrder+2);
  Standard_Integer i,j;
  Standard_Real Umin,Umax,Vmin,Vmax;
  Umin = Surf->Surface().FirstUParameter();
  Umax = Surf->Surface().LastUParameter();
  Vmin = Surf->Surface().FirstVParameter();
  Vmax = Surf->Surface().LastVParameter();
  
  DerSurf.SetValue(1, 0, D1UNormal);
    DerSurf.SetValue(0, 1, D1VNormal);
  DerSurf.SetValue(1, 1, d2uv);
  DerSurf.SetValue(2, 0, d2u);
  DerSurf.SetValue(0, 2, d2v);
  for(i=0;i<=MaxOrder+1;i++)
    for(j=i; j<=MaxOrder+2; j++ )
      if (i+j > 2) {
        DerSurf.SetValue(i,j,Surf->DN(U,V,i,j));
        if (i!=j) DerSurf.SetValue(j,i,Surf->DN(U,V,j,i));
      }
  
  for(i=0;i<=MaxOrder+1;i++)
    for(j=0;j<=MaxOrder+1;j++){
      DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurf));
      }
  
  CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,
                NStatus,Normal,OrderU,OrderV);
  if (NStatus != CSLib_Defined) Geom_UndefinedValue::Raise();

  D1UNormal = CSLib::DNNormal(1,0,DerNUV,OrderU,OrderV);
  D1VNormal = CSLib::DNNormal(0,1,DerNUV,OrderU,OrderV);
}

//=======================================================================
//function : NormalD2
//purpose  : computes Normal to Surface and its first and second derivatives
//=======================================================================
void NormalD2 (const Standard_Real U, const Standard_Real V, 
               const Handle(Adaptor3d_HSurface)& Surf, gp_Dir& Normal, 
               gp_Vec& D1UNormal, gp_Vec& D1VNormal,
               gp_Vec& D2UNormal, gp_Vec& D2VNormal, gp_Vec& D2UVNormal)
{
#ifdef CHECK  
  GeomAbs_Shape Cont = (Surf->Surface().UContinuity() < Surf->Surface().VContinuity()) ? 
    (Surf->Surface().UContinuity()) : (Surf->Surface().VContinuity());
     if (Cont == GeomAbs_C0 || Continuity == GeomAbs_C1 ||
         Cont == GeomAbs_C2) { Geom_UndefinedDerivative::Raise(); }
#endif
  gp_Vec d3u, d3uuv, d3uvv, d3v;
  gp_Pnt P;
  Surf->D3(U, V, P, D1UNormal, D1VNormal, 
           D2UNormal, D2VNormal, D2UVNormal,
           d3u,d3v, d3uuv, d3uvv);
  Standard_Real MagTol=0.000000001;
  CSLib_NormalStatus NStatus;
  CSLib::Normal (D1UNormal, D1VNormal, MagTol, NStatus, Normal);
  Standard_Integer MaxOrder;
  if (NStatus == CSLib_Defined) 
    MaxOrder=0;
  else 
    MaxOrder=3;
  Standard_Integer OrderU,OrderV;
  TColgp_Array2OfVec DerNUV(0,MaxOrder+2,0,MaxOrder+2);
  TColgp_Array2OfVec DerSurf(0,MaxOrder+3,0,MaxOrder+3);
  Standard_Integer i,j;
  
  Standard_Real Umin,Umax,Vmin,Vmax;
  Umin = Surf->Surface().FirstUParameter();
  Umax = Surf->Surface().LastUParameter();
  Vmin = Surf->Surface().FirstVParameter();
  Vmax = Surf->Surface().LastVParameter();
  
  DerSurf.SetValue(1, 0, D1UNormal);
  DerSurf.SetValue(0, 1, D1VNormal);
  DerSurf.SetValue(1, 1, D2UVNormal);
  DerSurf.SetValue(2, 0, D2UNormal);
  DerSurf.SetValue(0, 2, D2VNormal);
  DerSurf.SetValue(3, 0, d3u);
  DerSurf.SetValue(2, 1, d3uuv);
  DerSurf.SetValue(1, 2, d3uvv);
  DerSurf.SetValue(0, 3, d3v);
  for(i=0;i<=MaxOrder+2;i++)
    for(j=i; j<=MaxOrder+3; j++ )
      if (i+j > 3) {
        DerSurf.SetValue(i,j,Surf->DN(U,V,i,j));
        if (i!=j) DerSurf.SetValue(j,i,Surf->DN(U,V,j,i));
      }
  
  for(i=0;i<=MaxOrder+2;i++)
    for(j=0;j<=MaxOrder+2;j++){
      DerNUV.SetValue(i,j,CSLib::DNNUV(i,j,DerSurf));
    }
  
  CSLib::Normal(MaxOrder,DerNUV,MagTol,U,V,Umin,Umax,Vmin,Vmax,
                NStatus,Normal,OrderU,OrderV);
  if (NStatus != CSLib_Defined) Geom_UndefinedValue::Raise();
  
  D1UNormal = CSLib::DNNormal(1,0,DerNUV,OrderU,OrderV);
  D1VNormal = CSLib::DNNormal(0,1,DerNUV,OrderU,OrderV);
  D2UNormal = CSLib::DNNormal(2,0,DerNUV,OrderU,OrderV);
  D2VNormal = CSLib::DNNormal(0,2,DerNUV,OrderU,OrderV);
  D2UVNormal = CSLib::DNNormal(1,1,DerNUV,OrderU,OrderV);
}

GeomFill_Darboux::GeomFill_Darboux()
{
}

Handle(GeomFill_TrihedronLaw) GeomFill_Darboux::Copy() const
{
  Handle(GeomFill_Darboux) copy = new (GeomFill_Darboux)();
  if (!myCurve.IsNull()) copy->SetCurve(myCurve);
  return copy;
}

 Standard_Boolean GeomFill_Darboux::D0(const Standard_Real Param,gp_Vec& Tangent,gp_Vec& Normal,gp_Vec& BiNormal)
{
  gp_Pnt2d C2d;
  gp_Vec2d D2d;
  gp_Pnt S;
  gp_Vec dS_du, dS_dv;
  Handle(Adaptor2d_HCurve2d) myCurve2d = Adaptor3d_CurveOnSurfacePtr(&(myTrimmed->Curve()))->GetCurve();
  Handle(Adaptor3d_HSurface) mySupport = Adaptor3d_CurveOnSurfacePtr(&(myTrimmed->Curve()))->GetSurface();
  Standard_Integer OrderU, OrderV;
  myCurve2d->D1(Param, C2d, D2d);

//  Normal = dS_du.Crossed(dS_dv).Normalized();
  gp_Dir NormalDir;
  NormalD0(C2d.X(), C2d.Y(), mySupport, NormalDir, OrderU, OrderV);
  BiNormal.SetXYZ(NormalDir.XYZ()); 

  mySupport->D1(C2d.X(), C2d.Y(), S, dS_du, dS_dv);
//  if(D2d.Magnitude() <= Precision::Confusion())
//    DoTSingular(Param, Order);
  
  Tangent = D2d.X()*dS_du + D2d.Y()*dS_dv;
  Tangent.Normalize();

  Normal =  BiNormal;
  Normal ^= Tangent;
    

  return Standard_True;
}

 Standard_Boolean GeomFill_Darboux::D1(const Standard_Real Param,
				       gp_Vec& Tangent,gp_Vec& DTangent,
				       gp_Vec& Normal,gp_Vec& DNormal,
				       gp_Vec& BiNormal,gp_Vec& DBiNormal)
{
  gp_Pnt2d C2d;
  gp_Vec2d D2d, D2_2d;
  gp_Pnt S;
  gp_Vec dS_du, dS_dv, d2S_du, d2S_dv, d2S_duv, F, DF;
  Handle(Adaptor2d_HCurve2d) myCurve2d = Adaptor3d_CurveOnSurfacePtr(&(myTrimmed->Curve()))->GetCurve();
  Handle(Adaptor3d_HSurface) mySupport = Adaptor3d_CurveOnSurfacePtr(&(myTrimmed->Curve()))->GetSurface();
//  Standard_Integer Order;
  myCurve2d->D2(Param, C2d, D2d, D2_2d);
  mySupport->D2(C2d.X(), C2d.Y(), S, dS_du, dS_dv, 
		d2S_du, d2S_dv, d2S_duv);
//  if(D2d.Magnitude() <= Precision::Confusion())
//    DoTSingular(Param, Order);
  
  F = D2d.X()*dS_du + D2d.Y()*dS_dv;
  Tangent = F.Normalized();
  DF = D2_2d.X()*dS_du + D2_2d.Y()*dS_dv + 
       d2S_du*D2d.X()*D2d.X() + 2*d2S_duv*D2d.X()*D2d.Y() + 
       d2S_dv*D2d.Y()*D2d.Y();

  DTangent = FDeriv(F, DF);

  gp_Dir NormalDir;
  gp_Vec D1UNormal, D1VNormal;
  NormalD1(C2d.X(), C2d.Y(), mySupport, NormalDir, D1UNormal, D1VNormal);
  BiNormal.SetXYZ(NormalDir.XYZ());
  DBiNormal = D1UNormal*D2d.X() + D1VNormal*D2d.Y();
 
  Normal =  BiNormal;
  Normal ^= Tangent; 
  DNormal = BiNormal.Crossed(DTangent) + DBiNormal.Crossed(Tangent);
  
  return Standard_True;
}

 Standard_Boolean GeomFill_Darboux::D2(const Standard_Real Param,
				       gp_Vec& Tangent,gp_Vec& DTangent,gp_Vec& D2Tangent,
				       gp_Vec& Normal,gp_Vec& DNormal,gp_Vec& D2Normal,
				       gp_Vec& BiNormal,gp_Vec& DBiNormal,gp_Vec& D2BiNormal) 
{
  gp_Pnt2d C2d;
  gp_Vec2d D2d, D2_2d, D3_2d;
  gp_Pnt S;
  gp_Vec dS_du, dS_dv, d2S_du, d2S_dv, d2S_duv, 
         d3S_du, d3S_dv, d3S_duuv, d3S_duvv, F, DF, D2F;
  Handle(Adaptor2d_HCurve2d) myCurve2d = Adaptor3d_CurveOnSurfacePtr(&(myTrimmed->Curve()))->GetCurve();
  Handle(Adaptor3d_HSurface) mySupport = Adaptor3d_CurveOnSurfacePtr(&(myTrimmed->Curve()))->GetSurface();
//  Standard_Integer Order;
  myCurve2d->D3(Param, C2d, D2d, D2_2d, D3_2d);
  mySupport->D3(C2d.X(), C2d.Y(), S, dS_du, dS_dv, 
		d2S_du, d2S_dv, d2S_duv, d3S_du, d3S_dv, d3S_duuv, d3S_duvv);
//  if(D2d.Magnitude() <= Precision::Confusion())
//    DoTSingular(Param, Order);
  
  F = D2d.X()*dS_du + D2d.Y()*dS_dv;
  Tangent = F.Normalized();
  DF = D2_2d.X()*dS_du + D2_2d.Y()*dS_dv + 
       d2S_du*D2d.X()*D2d.X() + 2*d2S_duv*D2d.X()*D2d.Y() + 
       d2S_dv*D2d.Y()*D2d.Y();

  D2F = D3_2d.X()*dS_du + D3_2d.Y()*dS_dv + 
        D2_2d.X()*(D2d.X()*d2S_du + D2d.Y()*d2S_duv) +
	D2_2d.Y()*(D2d.X()*d2S_duv + D2d.Y()*d2S_dv) +
	D2d.X()*D2d.X()*(D2d.X()*d3S_du + D2d.Y()*d3S_duuv) +
	D2d.Y()*D2d.Y()*(D2d.X()*d3S_duvv + D2d.Y()*d3S_dv) +
	2*(
	   D2d.X()*D2_2d.X()*d2S_du + D2d.Y()*D2_2d.Y()*d2S_dv +
	   D2d.X()*D2d.Y()*(D2d.X()*d3S_duuv + D2d.Y()*d3S_duvv) +
	   d2S_duv*(D2_2d.X()*D2d.Y() + D2d.X()*D2_2d.Y())
	   );

  DTangent = FDeriv(F, DF);
  D2Tangent = DDeriv(F, DF, D2F);

  gp_Dir NormalDir;
  gp_Vec D1UNormal, D1VNormal, D2UNormal, D2VNormal, D2UVNormal;
  NormalD2(C2d.X(), C2d.Y(), mySupport, NormalDir, D1UNormal, D1VNormal, 
           D2UNormal, D2VNormal, D2UVNormal);
  BiNormal.SetXYZ(NormalDir.XYZ());
  DBiNormal = D1UNormal*D2d.X() + D1VNormal*D2d.Y();
  D2BiNormal = D1UNormal*D2_2d.X() + D1VNormal*D2_2d.Y() + D2UNormal*D2d.X()*D2d.X() +
    2*D2UVNormal*D2d.X()*D2d.Y() + D2VNormal*D2d.Y()*D2d.Y();

  Normal =  BiNormal;
  Normal ^= Tangent; 
  DNormal = BiNormal.Crossed(DTangent) + DBiNormal.Crossed(Tangent);
  D2Normal = BiNormal.Crossed(D2Tangent) + 2*DBiNormal.Crossed(DTangent) + D2BiNormal.Crossed(Tangent);

  return Standard_True;
}

 Standard_Integer GeomFill_Darboux::NbIntervals(const GeomAbs_Shape S) const
{
  return myCurve->NbIntervals(S);
}

 void GeomFill_Darboux::Intervals(TColStd_Array1OfReal& T,const GeomAbs_Shape S) const
{
  myCurve->Intervals(T, S);
}

 void GeomFill_Darboux::GetAverageLaw(gp_Vec& ATangent,gp_Vec& ANormal,gp_Vec& ABiNormal)
{
  Standard_Integer Num = 20; //order of digitalization
  gp_Vec T, N, BN;
  ATangent = gp_Vec(0, 0, 0);
  ANormal = gp_Vec(0, 0, 0);
  ABiNormal = gp_Vec(0, 0, 0);
  Standard_Real Step = (myTrimmed->LastParameter() - 
                        myTrimmed->FirstParameter()) / Num;
  Standard_Real Param;
  for (Standard_Integer i = 0; i <= Num; i++) {
    Param = myTrimmed->FirstParameter() + i*Step;
    if (Param > myTrimmed->LastParameter()) Param = myTrimmed->LastParameter();
    D0(Param, T, N, BN);
    ATangent += T;
    ANormal += N;
    ABiNormal += BN;
  }
  ATangent /= Num + 1;
  ANormal /= Num + 1;

  ATangent.Normalize();
  ABiNormal = ATangent.Crossed(ANormal).Normalized();
  ANormal = ABiNormal.Crossed(ATangent);  
}

 Standard_Boolean GeomFill_Darboux::IsConstant() const
{
  return (myCurve->GetType() == GeomAbs_Line);
}

 Standard_Boolean GeomFill_Darboux::IsOnlyBy3dCurve() const
{
  return Standard_False;
}
Commit	Line	Data
b311480e	1	// Created on: 1998-03-04
	2	// Created by: Roman BORISOV
	3	// Copyright (c) 1998-1999 Matra Datavision
973c2be1	4	// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e	5	//
973c2be1	6	// This file is part of Open CASCADE Technology software library.
b311480e	7	//
d5f74e42	8	// This library is free software; you can redistribute it and/or modify it under
d5f74e42	9	// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1	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.
b311480e	13	//
973c2be1	14	// Alternatively, this file may be used under the terms of Open CASCADE
973c2be1	15	// commercial license or contractual agreement.
7fd59977	16
42cf5bc1	17
7fd59977	18	#include <Adaptor2d_HCurve2d.hxx>
7fd59977	19	#include <Adaptor3d_CurveOnSurface.hxx>
42cf5bc1	20	#include <Adaptor3d_CurveOnSurfacePtr.hxx>
42cf5bc1	21	#include <Adaptor3d_HCurve.hxx>
7fd59977	22	#include <Adaptor3d_HCurveOnSurface.hxx>
7fd59977	23	#include <Adaptor3d_HSurface.hxx>
7fd59977	24	#include <CSLib.hxx>
7fd59977	25	#include <Geom_UndefinedValue.hxx>
	26	#include <GeomAdaptor_HSurface.hxx>
	27	#include <GeomAdaptor_Surface.hxx>
42cf5bc1	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>
7fd59977	37
92efcf78	38	IMPLEMENT_STANDARD_RTTIEXT(GeomFill_Darboux,GeomFill_TrihedronLaw)
92efcf78	39
7fd59977	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(FDF)/(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 - 2DF(FDF)/(NormaNorma))/Norma -
	59	F((DF.SquareMagnitude() + FD2F
	60	- 3(FDF)(FDF)/(NormaNorma))/(NormaNorma*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) Geom_UndefinedValue::Raise();
	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	Geom_UndefinedValue::Raise();
	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) {
0797d9d3	116	#ifdef OCCT_DEBUG
7fd59977	117	cout << U << ", " << V<< endl;
	118	for(i=0;i<=MaxOrder;i++)
	119	for(j=0;j<=MaxOrder;j++){
	120	cout <<"("<<i <<"," << j << ") = "<< DerSurf(i,j).X() <<","
	121	<< DerSurf(i,j).Y() <<"," << DerSurf(i,j).Z() << endl;
	122	}
	123	#endif
	124	Geom_UndefinedValue::Raise();
	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	Geom_UndefinedDerivative::Raise();
	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) Geom_UndefinedValue::Raise();
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) { Geom_UndefinedDerivative::Raise(); }
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) Geom_UndefinedValue::Raise();
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_duD2d.X()D2d.X() + 2d2S_duvD2d.X()*D2d.Y() +
325	d2S_dvD2d.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 = D1UNormalD2d.X() + D1VNormalD2d.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_duD2d.X()D2d.X() + 2d2S_duvD2d.X()*D2d.Y() +
365	d2S_dvD2d.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 = D1UNormalD2d.X() + D1VNormalD2d.Y();
387	D2BiNormal = D1UNormalD2_2d.X() + D1VNormalD2_2d.Y() + D2UNormalD2d.X()D2d.X() +
388	2D2UVNormalD2d.X()D2d.Y() + D2VNormalD2d.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	}