Logo Search packages:      
Sourcecode: opencascade version File versions  Download package

GeomFill_SectionPlacement.cxx

// File:    GeomFill_SectionPlacement.cxx
// Created: Mon Dec 15 17:09:47 1997
// Author:  Philippe MANGIN
//          <pmn@sgi29>


#include <GeomFill_SectionPlacement.ixx>

#include <GeomLib.hxx>
#include <Geom_Plane.hxx>
#include <GeomLProp_CLProps.hxx>
#include <GeomAbs_CurveType.hxx> 
#include <GeomAdaptor_HSurface.hxx>

#include <gp_Ax3.hxx>
#include <gp_Ax2.hxx>
#include <gp_Pnt.hxx>
#include <gp_Dir.hxx>
#include <gp_Trsf.hxx>
#include <TColgp_HArray1OfPnt.hxx>
#include <Bnd_Box.hxx>
#include <BndLib_Add3dCurve.hxx>

#include <Precision.hxx>
#include <gp.hxx>
#include <Extrema_ExtCC.hxx>
#include <Extrema_POnCurv.hxx>
#include <IntCurveSurface_HInter.hxx>
#include <IntCurveSurface_IntersectionPoint.hxx>

#include <Geom_BSplineCurve.hxx>
#include <Geom_TrimmedCurve.hxx>
#include <Geom_Line.hxx>
#include <Geom_Conic.hxx>
#include <Geom_BezierCurve.hxx>

#include <Geom_CartesianPoint.hxx>


//===============================================================
// Function :
// Purpose :
//===============================================================
static void Tangente(const Adaptor3d_Curve& Path,
                 const Standard_Real Param,
                 gp_Pnt& P,
                 gp_Vec& Tang)
{
  Path.D1(Param, P, Tang);
  Standard_Real Norm = Tang.Magnitude();

  for (Standard_Integer ii=2; (ii<12) && (Norm < Precision::Confusion()); 
       ii++) {
    Tang =  Path.DN(Param, ii);
    Norm = Tang.Magnitude();
  }

  if (Norm > 100*gp::Resolution()) Tang /= Norm; 
}
                   

static Standard_Real Penalite(const Standard_Real angle,
                        const Standard_Real dist)
{
  Standard_Real penal;
  
  if (dist < 1)
    penal = Sqrt(dist);
  else if (dist <2)
    penal = Pow(dist, 2);
  else
    penal = dist + 2;

  if (angle > 1.e-3) {
    penal += 1./angle -2./PI;
  }
  else {
    penal += 1.e3;
  }

  return penal;
}

static Standard_Real EvalAngle(const gp_Vec& V1, 
                         const gp_Vec& V2)
{
 Standard_Real angle;
 angle = V1.Angle(V2);
 if (angle > PI/2) angle = PI - angle;
 return angle;
}

static Standard_Integer NbSamples(const Handle(Geom_Curve)& aCurve)
{
  Standard_Real nbs = 100.; //on default

  Handle(Geom_Curve) theCurve = aCurve;
  if (aCurve->IsInstance(STANDARD_TYPE(Geom_TrimmedCurve)))
    theCurve = (Handle(Geom_TrimmedCurve)::DownCast(aCurve))->BasisCurve();

  if (theCurve->IsInstance(STANDARD_TYPE(Geom_Line)))
    nbs = 1;
  else if (theCurve->IsKind(STANDARD_TYPE(Geom_Conic)))
    nbs = 4;
  else if (theCurve->IsInstance(STANDARD_TYPE(Geom_BezierCurve)))
    {
      Handle(Geom_BezierCurve) BC = Handle(Geom_BezierCurve)::DownCast(theCurve);
      nbs = 3 + BC->NbPoles();
    }
  else if (theCurve->IsInstance(STANDARD_TYPE(Geom_BSplineCurve)))
    {
      Handle(Geom_BSplineCurve) BC = Handle(Geom_BSplineCurve)::DownCast(theCurve);
      nbs  = BC->NbKnots();
      nbs *= BC->Degree();
      Standard_Real ratio =
      (aCurve->LastParameter() - aCurve->FirstParameter())/(BC->LastParameter() - BC->FirstParameter());
      nbs *= ratio;
      if(nbs < 4.0)
      nbs = 4;
    }

  if (nbs > 300.)
    nbs = 300;

  return ((Standard_Integer)nbs);
}

//===============================================================
// Function :DistMini
// Purpose : Examine un extrema pour updater <Dist> & <Param>
//===============================================================
static void DistMini(const Extrema_ExtPC& Ext,
                 const Adaptor3d_Curve& C,
                 Standard_Real& Dist,
                 Standard_Real& Param)
{
  Standard_Real dist1, dist2;
  Standard_Integer ii;
  gp_Pnt P1, P2;
  Dist = RealLast();

  Ext.TrimmedDistances(dist1, dist2, P1, P2);
  if ( (dist1<Dist) || (dist2<Dist) ) {
    if (dist1 < dist2) {
      Dist = dist1;
      Param  = C.FirstParameter();
    }
    else {
      Dist = dist2;
      Param  = C.LastParameter();
    }
  }

  if (Ext.IsDone())
    {
      for (ii=1; ii<= Ext.NbExt(); ii++) {
      if (Ext.Value(ii) < Dist) {
        Dist = Ext.Value(ii);
        Param  =  Ext.Point(ii).Parameter();
      }
      }
    }
}
                 

//===============================================================
// Function : Constructeur
// Purpose :
//===============================================================
GeomFill_SectionPlacement::
GeomFill_SectionPlacement(const Handle(GeomFill_LocationLaw)& L,
                    const Handle(Geom_Geometry)& Section) :
                    myLaw(L), /* myAdpSection(Section),  mySection(Section), */
                    Dist(RealLast()), AngleMax(0.)
{

  done   = Standard_False;
  isplan = Standard_False;
  myIsPoint = Standard_False;

  if (Section->IsInstance(STANDARD_TYPE(Geom_CartesianPoint)))
    {
      myIsPoint = Standard_True;
      myPoint   = Handle(Geom_CartesianPoint)::DownCast(Section)->Pnt();
      isplan    = Standard_True;
    }
  else
    {
      Handle(Geom_Curve) CurveSection = Handle(Geom_Curve)::DownCast(Section);
      myAdpSection.Load(CurveSection);
      mySection = CurveSection;
    }

  Standard_Integer i, j, NbPoles=0;
  Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;

  // Boite d'encombrement de la section pour en deduire le gabarit
  Bnd_Box box;
  if (myIsPoint)
    box.Add(myPoint);
  else
    BndLib_Add3dCurve::Add(myAdpSection, 1.e-4, box);
  box.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);

  Standard_Real DX = aXmax - aXmin ;
  Standard_Real DY = aYmax - aYmin ;
  Standard_Real DZ = aZmax - aZmin ;
  Gabarit = Sqrt( DX*DX+ DY*DY + DZ*DZ )/2. ;

  Gabarit += Precision::Confusion(); // Cas des toute petite
  

  // Initialisation de TheAxe pour les cas singulier
  if (!myIsPoint)
    {
      gp_Pnt P;
      gp_Vec V;
      Tangente(myAdpSection, 
             (myAdpSection.FirstParameter()+myAdpSection.LastParameter())/2,
             P, V);
      TheAxe.SetLocation(P);
      TheAxe.SetDirection(V);
 
      // y a t'il un Plan moyen ?
      GeomAbs_CurveType TheType = myAdpSection.GetType();
      switch  (TheType) {
      case GeomAbs_Circle:
      {
        isplan = Standard_True;
        TheAxe =  myAdpSection.Circle().Axis();
        break;
      }
      case GeomAbs_Ellipse:
      {
        isplan = Standard_True;
        TheAxe =  myAdpSection.Ellipse().Axis();
        break;
      }
      case GeomAbs_Hyperbola:
      {
        isplan = Standard_True;
        TheAxe =  myAdpSection.Hyperbola().Axis();
        break;
      }
      case GeomAbs_Parabola:
      {
        isplan = Standard_True;
        TheAxe =  myAdpSection.Parabola().Axis();
        break;
      }
      case GeomAbs_Line:
      {
        NbPoles = 0; // Pas de Plan !!
        break;
      }
      case GeomAbs_BezierCurve:
      case GeomAbs_BSplineCurve:
      {
        NbPoles = myAdpSection.NbPoles();
        break;
      }
      default:
      NbPoles = 21;
      }
      
      
      if (!isplan && NbPoles>2)
      {
        // Calcul d'un plan moyen.
        Handle(TColgp_HArray1OfPnt) Pnts;
        Standard_Real first = myAdpSection.FirstParameter();
        Standard_Real last =  myAdpSection.LastParameter();
        Standard_Real t, delta;
        if (myAdpSection.GetType() == GeomAbs_BSplineCurve)
          {
            Handle(Geom_BSplineCurve) BC =
            Handle(Geom_BSplineCurve)::DownCast(myAdpSection.Curve());
            Standard_Integer I1, I2, I3, I4;
            BC->LocateU( first, Precision::Confusion(), I1, I2 );
            BC->LocateU( last,  Precision::Confusion(), I3, I4 );
            Standard_Integer NbKnots = I3 - I2 + 1;
            
            Standard_Integer NbLocalPnts = 10;
            Standard_Integer NbPnts = (NbKnots-1) * NbLocalPnts;
            if (I1 != I2)
            NbPnts += NbLocalPnts;
            if (I3 != I4)
            NbPnts += NbLocalPnts;
            if (!myAdpSection.IsClosed())
            NbPnts++;
            Pnts = new TColgp_HArray1OfPnt(1, NbPnts);
            Standard_Integer nb = 1;
            if (I1 != I2)
            {
              Standard_Real locallast = (BC->Knot(I2) < last)? BC->Knot(I2) : last;
              delta = (locallast - first) / NbLocalPnts;
              for (j = 0; j < NbLocalPnts; j++)
                {
                  t = first + j*delta;
                  Pnts->SetValue( nb++, myAdpSection.Value(t) );
                }
            }
            for (i = I2; i < I3; i++)
            {
              t = BC->Knot(i);
              delta = (BC->Knot(i+1) - t) / NbLocalPnts;
              for (j = 0; j < NbLocalPnts; j++)
                {
                  t += delta;
                  Pnts->SetValue( nb++, myAdpSection.Value(t) );
                }
            }
            if (I3 != I4 && first < BC->Knot(I3))
            {
              t = BC->Knot(I3);
              delta = (last - t) / NbLocalPnts;
              for (j = 0; j < NbLocalPnts; j++)
                {
                  t += delta;
                  Pnts->SetValue( nb++, myAdpSection.Value(t) );
                }
            }
            if (!myAdpSection.IsClosed())
            Pnts->SetValue( nb, myAdpSection.Value(last) );
          }
        else // other type
          {
            Standard_Integer NbPnts = NbPoles-1;
            if (!myAdpSection.IsClosed())
            NbPnts++;
            Pnts = new TColgp_HArray1OfPnt(1, NbPnts);
            delta = (last - first) / (NbPoles-1);
            for (i = 0; i < NbPoles-1; i++)
            {
              t = first + i*delta;
              Pnts->SetValue( i+1, myAdpSection.Value(t) );
            }
            if (!myAdpSection.IsClosed())
            Pnts->SetValue( NbPnts, myAdpSection.Value(last) );
          }
        
        Standard_Boolean issing;
        gp_Ax2 axe;
        GeomLib::AxeOfInertia(Pnts->Array1(), axe, issing, Precision::Confusion());
        if (!issing) {
          isplan = Standard_True;
          TheAxe.SetLocation(axe.Location());
          TheAxe.SetDirection(axe.Direction());
        }
      }
      
      
      myExt.Initialize(myAdpSection, 
                   myAdpSection.FirstParameter(), 
                   myAdpSection.LastParameter(), 
                   Precision::Confusion());
    }
}

//===============================================================
// Function :SetLocation
// Purpose :
//===============================================================
void GeomFill_SectionPlacement::
00365 SetLocation(const Handle(GeomFill_LocationLaw)& L) 
{
  myLaw = L;
}

//===============================================================
// Function : Perform
// Purpose : Le plus simple
//===============================================================
void GeomFill_SectionPlacement::Perform(const Standard_Real Tol) 
{
  Handle(Adaptor3d_HCurve) Path;
  Path = myLaw->GetCurve();
  Perform(Path, Tol);
}

//===============================================================
// Function :Perform
// Purpose : Recherche automatique
//===============================================================
void GeomFill_SectionPlacement::Perform(const Handle(Adaptor3d_HCurve)& Path,
                              const Standard_Real Tol) 
{
  Standard_Real IntTol = 1.e-5;
  Standard_Real DistCenter = Precision::Infinite();

  if (myIsPoint)
    {
      Extrema_ExtPC Projector(myPoint, Path->Curve(), Precision::Confusion());
      DistMini( Projector, Path->Curve(), Dist, PathParam );
      AngleMax = PI/2;
    }
  else
    {
      PathParam = Path->FirstParameter();
      SecParam =  myAdpSection.FirstParameter();
      
      Standard_Real distaux, taux, alpha;
      gp_Pnt PonPath, PonSec, P;
      gp_Vec VRef, dp1;
      VRef.SetXYZ(TheAxe.Direction().XYZ());
      
      Tangente( Path->Curve(), PathParam, PonPath, dp1);
      PonSec = myAdpSection.Value(SecParam);
      Dist =  PonPath.Distance(PonSec);
      if (Dist > Tol) { // On Cherche un meilleur point sur la section
      myExt.Perform(PonPath);  
      if ( myExt.IsDone() ) { 
        DistMini(myExt, myAdpSection, Dist, SecParam);
        PonSec = myAdpSection.Value(SecParam);
      }
      }
      AngleMax = EvalAngle(VRef, dp1);
      if (isplan) AngleMax = PI/2 - AngleMax;
      
      Standard_Boolean Trouve = Standard_False; 
      Standard_Integer ii;
      
      if (isplan) {
      // (1.1) Distances Point-Plan
      Standard_Real DistPlan;
      gp_Vec V1 (PonPath, TheAxe.Location());
      DistPlan = Abs(V1.Dot(VRef));
      if (DistPlan <= IntTol)
        DistCenter = V1.Magnitude();
      
      gp_Pnt Plast = Path->Value( Path->LastParameter() );
      V1.SetXYZ( TheAxe.Location().XYZ() - Plast.XYZ() );
      DistPlan = Abs(V1.Dot(VRef));
      if (DistPlan <= IntTol)
        {
          Standard_Real aDist = V1.Magnitude();
          if (aDist < DistCenter)
            {
            DistCenter = aDist;
            PonPath = Plast;
            PathParam = Path->LastParameter();
            }
        }
      
      // (1.2) Intersection Plan-courbe
      gp_Ax3 axe (TheAxe.Location(), TheAxe.Direction());
      Handle(Geom_Plane) plan = new (Geom_Plane)(axe);
      Handle(GeomAdaptor_HSurface) adplan = 
        new (GeomAdaptor_HSurface)(plan);
      IntCurveSurface_HInter Intersector;
      Intersector.Perform(Path, adplan);
      if (Intersector.IsDone())
        {
          Standard_Real w;
          gp_Pnt P;
          Standard_Real aDist;
          for (ii=1; ii<=Intersector.NbPoints(); ii++)
            {
            w = Intersector.Point(ii).W();
            P = Path->Value( w );
            aDist = P.Distance( TheAxe.Location() );
            if (aDist < DistCenter)
              {
                DistCenter = aDist;
                PonPath = P;
                PathParam = w;
              }
            }
        }
      if (!Intersector.IsDone() || Intersector.NbPoints() == 0)
        {
          Standard_Integer NbPnts = NbSamples( mySection );
          TColgp_Array1OfPnt Pnts( 1, NbPnts+1 );
          Standard_Real delta = (mySection->LastParameter()-mySection->FirstParameter())/NbPnts;
          for (ii = 0; ii <= NbPnts; ii++)
            Pnts(ii+1) = mySection->Value( mySection->FirstParameter() + ii*delta );
          
          gp_Pnt BaryCenter;
          gp_Dir Xdir, Ydir;
          Standard_Real Xgap, Ygap, Zgap;
          GeomLib::Inertia( Pnts, BaryCenter, Xdir, Ydir, Xgap, Ygap, Zgap );
          
          gp_Pnt Pfirst = Path->Value( Path->FirstParameter() );
          if (Pfirst.Distance(BaryCenter) < Plast.Distance(BaryCenter))
            PathParam = Path->FirstParameter();
          else
            {
            PathParam = Path->LastParameter();
            Tangente( Path->Curve(), PathParam, PonPath, dp1);
            PonSec = myAdpSection.Value(SecParam);
            Dist =  PonPath.Distance(PonSec);
            if (Dist > Tol) { // On Cherche un meilleur point sur la section
              myExt.Perform(PonPath);  
              if ( myExt.IsDone() ) { 
                DistMini(myExt, myAdpSection, Dist, SecParam);
                PonSec = myAdpSection.Value(SecParam);
              }
            }
            AngleMax = EvalAngle(VRef, dp1);
            AngleMax = PI/2 - AngleMax;
            }
        }
    
      /*
         // (1.1) Distances Point-Plan
         Standard_Real DistPlan;
         gp_Vec V1 (PonPath, TheAxe.Location());
         DistPlan = Abs(V1.Dot(VRef));
         
         // On examine l'autre extremite
         gp_Pnt P;
         Tangente(Path->Curve(), Path->LastParameter(), P, dp1);
         V1.SetXYZ(TheAxe.Location().XYZ()-P.XYZ());
         if  (Abs(V1.Dot(VRef)) <= DistPlan ) { // On prend l'autre extremite
         alpha =  PI/2 - EvalAngle(VRef, dp1);
         distaux =  PonPath.Distance(PonSec);
         if (distaux > Tol) {
         myExt.Perform(P);  
         if ( myExt.IsDone() ) 
         DistMini(myExt, myAdpSection, distaux, taux);
         }
         else 
         taux = SecParam;
         
         if (Choix(distaux, alpha)) {
         Dist = distaux;
         AngleMax = alpha;
         PonPath = P;
         PathParam = Path->LastParameter();
         }
         }
         
         // (1.2) Intersection Plan-courbe
         gp_Ax3 axe (TheAxe.Location(), TheAxe.Direction());
         Handle(Geom_Plane) plan = new (Geom_Plane)(axe);
         Handle(GeomAdaptor_HSurface) adplan = 
         new (GeomAdaptor_HSurface)(plan);
         IntCurveSurface_HInter Intersector;
         Intersector.Perform(Path, adplan);
         if (Intersector.IsDone()) {
         Standard_Real w;
         gp_Vec V;
         for (ii=1; ii<=Intersector.NbPoints(); ii++) {
         w = Intersector.Point(ii).W();
         //(1.3) test d'angle
         Tangente( Path->Curve(), w, P, V);
         alpha = PI/2 - EvalAngle(V, VRef);
         //(1.4) Test de distance Point-Courbe
         myExt.Perform(P);  
         if ( myExt.IsDone() ) {
         DistMini(myExt, myAdpSection, distaux, taux);
         if (Choix(distaux, alpha)) {
         Dist = distaux;
         SecParam  = taux;
         AngleMax = alpha;
         PonPath = P;
         PathParam = w;
         PonSec = myAdpSection.Value(SecParam);
         }
         }
         else {
         distaux = P.Distance(PonSec);
         if (Choix(distaux, alpha)) {
         Dist = distaux;
         AngleMax = alpha;
         PonPath = P;
         PathParam = w;
         }
         }
         }
         }
    */
#if DEB
      if (Intersector.NbPoints() == 0) {
        Intersector.Dump();
      } 
#endif
      } 
      
      // Cas General
      if (! isplan) {
      // (2.1) Distance avec les extremites ...
      myExt.Perform(PonPath);  
      if ( myExt.IsDone() ) {
        DistMini(myExt, myAdpSection, distaux, taux);
        if (distaux < Dist) {
          Dist = distaux;
          SecParam  = taux;
        } 
      }
      Trouve = (Dist <= Tol);
      if (!Trouve) {
        Tangente( Path->Curve(), Path->LastParameter(), P, dp1);
        alpha = EvalAngle(VRef, dp1);
        myExt.Perform(P);     
        if ( myExt.IsDone() ) {
          if ( myExt.IsDone() ) {
            DistMini(myExt, myAdpSection, distaux, taux);
            if (Choix(distaux, alpha)) {
            Dist = distaux;
            SecParam  = taux;
            AngleMax = alpha;
            PonPath = P;
            PathParam = Path->LastParameter();
            }
          }
        }
        Trouve = (Dist <= Tol); 
      }
      
      // (2.2) Distance courbe-courbe
      if (!Trouve) {
        Extrema_ExtCC Ext(Path->Curve(), myAdpSection, 
                      Path->FirstParameter(), Path->LastParameter(),
                      myAdpSection.FirstParameter(), 
                      myAdpSection.LastParameter(),
                      Path->Resolution(Tol/100), 
                      myAdpSection.Resolution(Tol/100));
        if (Ext.IsDone()) {
          Extrema_POnCurv P1, P2;
          for (ii=1; ii<=Ext.NbExt(); ii++) {
            distaux = Ext.Value(ii);
            Ext.Points(ii, P1, P2);
            Tangente(Path->Curve(), P1.Parameter(), P, dp1);
            alpha =  EvalAngle(VRef, dp1);
            if (Choix(distaux, alpha)) {
            Trouve = Standard_True;
            Dist = distaux;
            PathParam = P1.Parameter();
            SecParam = P2.Parameter();
            PonSec =  P2.Value();
            PonPath = P;
            AngleMax = alpha;
            }
          }
        }
        if (!Trouve){ 
          // Si l'on a toujours rien, on essai une distance point/path
          // c'est la derniere chance.
          Extrema_ExtPC PExt;
          PExt.Initialize(Path->Curve(), 
                      Path->FirstParameter(),  
                      Path->LastParameter(),
                      Precision::Confusion());
          PExt.Perform(PonSec);
          if ( PExt.IsDone() ) {
            // modified for OCC13595  Tue Oct 17 15:00:08 2006.BEGIN
            //          DistMini(PExt, myAdpSection, distaux, taux);
            DistMini(PExt, Path->Curve(), distaux, taux);
            // modified for OCC13595  Tue Oct 17 15:00:11 2006.END
            Tangente(Path->Curve(), taux, P, dp1);
            alpha = EvalAngle(VRef, dp1);
            if (Choix(distaux, alpha)) {
            Dist = distaux;
            PonPath = P;
            AngleMax = alpha;
            PathParam = taux;
            }
          }
        }
      }
      }
    }

  done = Standard_True;
}


//===============================================================
// Function : Perform
// Purpose : Calcul le placement pour un parametre donne.
//===============================================================
void GeomFill_SectionPlacement::Perform(const Standard_Real Param,
                              const Standard_Real Tol) 
{
  done = Standard_True;
  Handle(Adaptor3d_HCurve) Path;
  Path = myLaw->GetCurve();

  PathParam = Param;
  if (myIsPoint)
    {
      gp_Pnt PonPath = Path->Value( PathParam );
      Dist = PonPath.Distance(myPoint);
      AngleMax = PI/2;
    }
  else
    {
      SecParam =  myAdpSection.FirstParameter();
      
      //  Standard_Real distaux, taux, alpha;
      //  gp_Pnt PonPath, PonSec, P;
      gp_Pnt PonPath, PonSec;
      gp_Vec VRef, dp1;
      VRef.SetXYZ(TheAxe.Direction().XYZ()); 
      
      Tangente( Path->Curve(), PathParam, PonPath, dp1);
      PonSec = myAdpSection.Value(SecParam);
      Dist =  PonPath.Distance(PonSec);
      if (Dist > Tol) { // On Cherche un meilleur point sur la section
      myExt.Perform(PonPath);  
      if ( myExt.IsDone() ) { 
        DistMini(myExt, myAdpSection, Dist, SecParam);
        PonSec = myAdpSection.Value(SecParam);
      }
      }
      AngleMax = EvalAngle(VRef, dp1);
      if (isplan) AngleMax = PI/2 - AngleMax;
    }

  done = Standard_True;
}

//===============================================================
// Function :IsDone
// Purpose :
//===============================================================
 Standard_Boolean GeomFill_SectionPlacement::IsDone() const
{
  return done; 
}

//===============================================================
// Function : ParameterOnPath
// Purpose :
//===============================================================
 Standard_Real GeomFill_SectionPlacement::ParameterOnPath() const
{
  return PathParam;
}

//===============================================================
// Function : ParameterOnSection
// Purpose :
//===============================================================
 Standard_Real GeomFill_SectionPlacement::ParameterOnSection() const
{
 return SecParam;
}

//===============================================================
// Function : Distance
// Purpose :
//===============================================================
 Standard_Real GeomFill_SectionPlacement::Distance() const
{
  return Dist;
}

//===============================================================
// Function : Angle
// Purpose :
//===============================================================
 Standard_Real GeomFill_SectionPlacement::Angle() const
{
  return AngleMax;
}

//===============================================================
// Function : Transformation
// Purpose :
//===============================================================
 gp_Trsf GeomFill_SectionPlacement::Transformation(
         const Standard_Boolean WithTranslation,
       const Standard_Boolean WithCorrection) const
{
  gp_Vec V;
  gp_Mat M;
  gp_Dir DT, DN, D;
// modified by NIZHNY-MKK  Fri Oct 17 15:27:07 2003
  gp_Pnt P(0., 0., 0.), PSection(0., 0., 0.);

  // Calcul des reperes
  myLaw->D0(PathParam, M, V);

  P.SetXYZ(V.XYZ());
  D.SetXYZ (M.Column(3));
  DN.SetXYZ(M.Column(1));
  gp_Ax3 Paxe(P, D, DN);


  if (WithTranslation || WithCorrection) {
    if (myIsPoint)
      PSection = myPoint;
    else
      PSection = mySection->Value(SecParam);   
  }
  // modified by NIZHNY-MKK  Wed Oct  8 15:03:19 2003.BEGIN
  gp_Trsf Rot;

  if (WithCorrection && !myIsPoint) { 
    Standard_Real angle;
    
    if (!isplan)
      Standard_Failure::Raise("Illegal usage: can't rotate non-planar profile");
    
    gp_Dir ProfileNormal = TheAxe.Direction();
    gp_Dir SpineStartDir = Paxe.Direction();

    if (! ProfileNormal.IsParallel( SpineStartDir, Precision::Angular() )) {
      gp_Dir DirAxeOfRotation = ProfileNormal ^ SpineStartDir;
      angle = ProfileNormal.AngleWithRef( SpineStartDir, DirAxeOfRotation );
      gp_Ax1 AxeOfRotation( TheAxe.Location(), DirAxeOfRotation );
      Rot.SetRotation( AxeOfRotation, angle );
    }
    PSection.Transform(Rot);
  }
  // modified by NIZHNY-MKK  Wed Oct  8 15:03:21 2003.END

  if (WithTranslation) { 
    P.ChangeCoord().SetLinearForm(-1,PSection.XYZ(), 
                          V.XYZ());
  }
  else {
    P.SetCoord(0., 0., 0.);
  }

  gp_Ax3 Saxe(P, gp::DZ(), gp::DX());

  // Transfo
  gp_Trsf Tf;
  Tf.SetTransformation(Saxe, Paxe);

  if (WithCorrection) {
    // modified by NIZHNY-MKK  Fri Oct 17 15:26:36 2003.BEGIN
    //     Standard_Real angle;
    //     gp_Trsf Rot;
    
    //     if (!isplan)
    //       Standard_Failure::Raise("Illegal usage: can't rotate non-planar profile");
    
    //     gp_Dir ProfileNormal = TheAxe.Direction();
    //     gp_Dir SpineStartDir = Paxe.Direction();
    //     if (! ProfileNormal.IsParallel( SpineStartDir, Precision::Angular() ))
    //       {
    //      gp_Dir DirAxeOfRotation = ProfileNormal ^ SpineStartDir;
    //      angle = ProfileNormal.AngleWithRef( SpineStartDir, DirAxeOfRotation );
    //      gp_Ax1 AxeOfRotation( TheAxe.Location(), DirAxeOfRotation );
    //      Rot.SetRotation( AxeOfRotation, angle );
    //       }
    // modified by NIZHNY-MKK  Fri Oct 17 15:26:42 2003.END

    Tf *= Rot;
  }

  return Tf;
}

//===============================================================
// Function : Section
// Purpose :
//===============================================================
 Handle(Geom_Curve) GeomFill_SectionPlacement::
 Section(const Standard_Boolean WithTranslation) const
{
  Handle(Geom_Curve) TheSection =  
    Handle(Geom_Curve)::DownCast(mySection->Copy());
  TheSection->Transform(Transformation(WithTranslation, Standard_False));
  return TheSection;
}


//===============================================================
// Function :
// Purpose :
//===============================================================
 Handle(Geom_Curve) GeomFill_SectionPlacement::
 ModifiedSection(const Standard_Boolean WithTranslation) const
{
  Handle(Geom_Curve) TheSection =  
    Handle(Geom_Curve)::DownCast(mySection->Copy());
  TheSection->Transform(Transformation(WithTranslation, Standard_True));
  return TheSection;
}

//===============================================================
// Function :SectionAxis
// Purpose :
//===============================================================
 void GeomFill_SectionPlacement::SectionAxis(const gp_Mat& M,
                                   gp_Vec& T,
                                   gp_Vec& N, 
                                   gp_Vec& BN) const
{
  Standard_Real Eps = 1.e-10;
  gp_Dir D;
  gp_Vec PathNormal;
  GeomLProp_CLProps CP(mySection, SecParam, 2, Eps);
  if (CP.IsTangentDefined()) {
    CP.Tangent(D);
    T.SetXYZ(D.XYZ());
    T.Normalize();
    if (CP.Curvature() > Eps) {
      CP.Normal(D);
      N.SetXYZ(D.XYZ());
    }
    else { 
      // Cas ambigu, on essai de recuperer la normal a la trajectoire
      // Reste le probleme des points d'inflexions, qui n'est pas 
      // bien traiter par LProp (pas de recuperation de la normal) : 
      // A voir...
      PathNormal.SetXYZ(M.Column(1));
      PathNormal.Normalize();
      BN = T ^ PathNormal;
      if (BN.Magnitude() > Eps) {
      BN.Normalize();
      //N = BN ^ T;
      }
      N = BN ^ T;
    }
  }
  else { // Cas indefinie, on prend le triedre 
        // complet sur la trajectoire
    T.SetXYZ(M.Column(3));
    N.SetXYZ(M.Column(2));
  }
  BN = T ^ N;
}


//===============================================================
// Function :Choix
// Purpose : Decide si le couple (dist, angle) est "meilleur" que
// couple courrant.
//===============================================================
 Standard_Boolean GeomFill_SectionPlacement::Choix(const Standard_Real dist,
                                       const Standard_Real  angle) const
{
  Standard_Real evoldist, evolangle;
  evoldist = dist - Dist;
  evolangle = angle - AngleMax;
  // (1) Si la gain en distance est > que le gabarit, on prend
  if (evoldist < - Gabarit) 
    return Standard_True;

 //  (2) si l'ecart en distance est de l'ordre du gabarit
  if (Abs(evoldist) < Gabarit) {
//  (2.1) si le gain en angle est important on garde
    if  (evolangle > 0.5)
       return Standard_True;
 //  (2.2) si la variation d'angle est moderee on evalue
 //  une fonction de penalite
    if (Penalite(angle, dist/Gabarit) < Penalite(AngleMax, Dist/Gabarit) )
      return Standard_True;
  }

  return Standard_False;
}









Generated by  Doxygen 1.6.0   Back to index