// Created on: 1995-04-26
// Created by: Flore Lantheaume
// Copyright (c) 1995-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 <Adaptor3d_TopolTool.hxx>
#include <BRepBlend_Chamfer.hxx>
#include <BRepBlend_ChamfInv.hxx>
#include <BRepBlend_ConstThroat.hxx>
#include <BRepBlend_ConstThroatInv.hxx>
#include <BRepBlend_ConstThroatWithPenetration.hxx>
#include <BRepBlend_ConstThroatWithPenetrationInv.hxx>
#include <BRepBlend_ChAsym.hxx>
#include <BRepBlend_ChAsymInv.hxx>
#include <BRepBlend_Line.hxx>
#include <BRepBlend_Walking.hxx>
#include <BRepTools.hxx>
#include <ChFi3d.hxx>
#include <ChFi3d_Builder_0.hxx>
#include <ChFi3d_ChBuilder.hxx>
#include <ChFiDS_ChamfSpine.hxx>
#include <ChFiDS_CircSection.hxx>
#include <ChFiDS_SurfData.hxx>
#include <NCollection_Sequence.hxx>
#include <NCollection_HSequence.hxx>
#include <ChFiDS_ElSpine.hxx>
#include <ChFiDS_Regul.hxx>
#include <NCollection_List.hxx>
#include <ChFiDS_Stripe.hxx>
#include <NCollection_Array1.hxx>
#include <NCollection_HArray1.hxx>
#include <ChFiDS_Spine.hxx>
#include <ElSLib.hxx>
#include <Extrema_GenLocateExtPS.hxx>
#include <Standard_ConstructionError.hxx>
#include <Standard_DomainError.hxx>
#include <Standard_NotImplemented.hxx>
#include <TopoDS_Edge.hxx>
#include <TopoDS_Face.hxx>
#include <TopoDS_Shape.hxx>
#include <TopoDS_Vertex.hxx>
#include <TopOpeBRepBuild_HBuilder.hxx>
#include <TopOpeBRepDS_HDataStructure.hxx>
#include <memory>

#ifdef OCCT_DEBUG
extern bool ChFi3d_GettraceCHRON();
#endif

//=======================================================================
// function : SearchCommonFaces
// purpose  : search the 2 common faces <F1> and <F2> of the edge <E>
//           uses the EFMap and take the 2 first good faces
//=======================================================================

void SearchCommonFaces(const ChFiDS_Map&  EFMap,
                       const TopoDS_Edge& E,
                       TopoDS_Face&       F1,
                       TopoDS_Face&       F2)
{
  TopoDS_Face                              Fc;
  NCollection_List<TopoDS_Shape>::Iterator It;

  F1.Nullify();
  F2.Nullify();
  for (It.Initialize(EFMap(E)); It.More(); It.Next())
  {
    Fc = TopoDS::Face(It.Value());
    if (F1.IsNull())
      F1 = Fc;
    else if (!Fc.IsSame(F1))
    {
      F2 = Fc;
      break;
    }
  }

  if (!F1.IsNull() && F2.IsNull() && BRepTools::IsReallyClosed(E, F1))
    F2 = F1;
}

//=======================================================================
// function : ExtentSpinesOnCommonFace
// purpose  : Extend spines of two chamfers by distance dis1,dis2
//           on their common face
//           Two guide lines Spine1 and Spine2 cross in V
//           isfirst(i) = False if Spine(i) is oriented to V (i = 1,2)
//=======================================================================

void ExtentSpineOnCommonFace(occ::handle<ChFiDS_Spine>& Spine1,
                             occ::handle<ChFiDS_Spine>& Spine2,
                             const TopoDS_Vertex&       V,
                             const double               dis1,
                             const double               dis2,
                             const bool                 isfirst1,
                             const bool                 isfirst2)
{
  double tolesp = 1.e-7;

  // alpha, the opening angle between two
  // tangents of two guidelines in V is found
  double tga1, tga2;
  double d1plus = 0., d2plus = 0.;

  gp_Pnt tmp;
  gp_Vec tg1, tg2;
  Spine1->D1(Spine1->Absc(V), tmp, tg1);
  Spine2->D1(Spine2->Absc(V), tmp, tg2);
  tg1.Normalize();
  tg2.Normalize();
  if (isfirst1)
    tg1.Reverse();
  if (isfirst2)
    tg2.Reverse();

  double cosalpha, sinalpha;
  cosalpha = tg1.Dot(tg2);
  sinalpha = sqrt(1 - cosalpha * cosalpha);

  // a1+a2 = alpha
  double temp;
  temp = cosalpha + dis2 / dis1;
  if (std::abs(temp) > tolesp)
  {
    tga1   = sinalpha / temp;
    d1plus = dis1 / tga1;
  }
  temp = cosalpha + dis1 / dis2;
  if (std::abs(temp) > tolesp)
  {
    tga2   = sinalpha / temp;
    d2plus = dis2 / tga2;
  }

  // extension by the calculated distance
  if (d1plus > 0.)
  {
    d1plus *= 3.;
    if (isfirst1)
    {
      Spine1->SetFirstParameter(-d1plus);
      Spine1->SetFirstTgt(0.);
    }
    else
    {
      double param = Spine1->LastParameter(Spine1->NbEdges());
      Spine1->SetLastParameter(d1plus + param);
      Spine1->SetLastTgt(param);
    }
  }
  if (d2plus > 0.)
  {
    d2plus *= 1.5;
    if (isfirst2)
    {
      Spine2->SetFirstParameter(-d2plus);
      Spine2->SetFirstTgt(0.);
    }
    else
    {
      double param = Spine2->LastParameter(Spine2->NbEdges());
      Spine2->SetLastParameter(d2plus + param);
      Spine2->SetLastTgt(param);
    }
  }
}

//=================================================================================================

ChFi3d_ChBuilder::ChFi3d_ChBuilder(const TopoDS_Shape& S, const double Ta)
    : ChFi3d_Builder(S, Ta)
{
  myMode = ChFiDS_ClassicChamfer;
}

//=======================================================================
// function : Add
// purpose  : create a new stripe with a spine containing the edge <E>
//           add on the spine the tangential edges to <E>
//=======================================================================

void ChFi3d_ChBuilder::Add(const TopoDS_Edge& E)
{
  TopoDS_Face dummy;

  if (!Contains(E) && myEFMap.Contains(E))
  {
    occ::handle<ChFiDS_Stripe> Stripe    = new ChFiDS_Stripe();
    occ::handle<ChFiDS_Spine>& Sp        = Stripe->ChangeSpine();
    Sp                                   = new ChFiDS_ChamfSpine(tolesp);
    occ::handle<ChFiDS_ChamfSpine> Spine = occ::down_cast<ChFiDS_ChamfSpine>(Sp);

    TopoDS_Edge E_wnt = E;
    E_wnt.Orientation(TopAbs_FORWARD);
    Spine->SetEdges(E_wnt);
    if (PerformElement(Spine, -1, dummy))
    {
      PerformExtremity(Spine);
      Spine->Load();
      myListStripe.Append(Stripe);
    }
  }
}

//=======================================================================
// function : Add
// purpose  : create a new stripe with a ChamfSpine containing the edge <E>
//           with the distance <Dis>   on the face <F>
//           . Add the tangential edges continuous to E in the spine
//
//=======================================================================

void ChFi3d_ChBuilder::Add(const double Dis, const TopoDS_Edge& E)
{
  if (!Contains(E) && myEFMap.Contains(E))
  {

    TopoDS_Face dummy;

    TopoDS_Edge E_wnt = E;
    E_wnt.Orientation(TopAbs_FORWARD);

    occ::handle<ChFiDS_Stripe> Stripe    = new ChFiDS_Stripe();
    occ::handle<ChFiDS_Spine>& Sp        = Stripe->ChangeSpine();
    Sp                                   = new ChFiDS_ChamfSpine(tolesp);
    occ::handle<ChFiDS_ChamfSpine> Spine = occ::down_cast<ChFiDS_ChamfSpine>(Sp);

    Spine->SetMode(myMode);

    Spine->SetEdges(E_wnt);
    if (PerformElement(Spine, -1, dummy))
    {
      Spine->Load();
      myListStripe.Append(Stripe);

      Spine->SetDist(Dis);

      PerformExtremity(Spine);
    }
  }
}

//=======================================================================
// function : SetDist
// purpose  : set the distances <Dis>of the contour of
//           index IC
//=======================================================================

void ChFi3d_ChBuilder::SetDist(const double Dis, const int IC, const TopoDS_Face& F)
{

  if (IC <= NbElements())
  {
    occ::handle<ChFiDS_ChamfSpine> csp = occ::down_cast<ChFiDS_ChamfSpine>(Value(IC));

    // Search the first edge which has a common face equal to F
    TopoDS_Face F1, F2, FirstF1, FirstF2;
    // TopAbs_Orientation Or1,Or2;
    // int Choix, ChoixConge;
    BRepAdaptor_Surface Sb1, Sb2;
    int                 i     = 1;
    bool                Found = false;
    while ((i <= csp->NbEdges()) && (!Found))
    {
      SearchCommonFaces(myEFMap, csp->Edges(i), F1, F2);
      if (i == 1)
      {
        FirstF1 = F1;
        FirstF2 = F2;
      }
      Found = (F1.IsSame(F) || F2.IsSame(F));
      i++;
    }

    if (Found)
    {
      if (F2.IsSame(F))
      {
        F2 = F1;
        F1 = F;
      }
      csp->SetDist(Dis);
    }
    else
      throw Standard_DomainError("the face is not common to any of edges of the contour");
  }
}

//=================================================================================================

void ChFi3d_ChBuilder::GetDist(const int IC, double& Dis) const
{
  occ::handle<ChFiDS_ChamfSpine> chsp = occ::down_cast<ChFiDS_ChamfSpine>(Value(IC));
  chsp->GetDist(Dis);
}

//=======================================================================
// function : Add
// purpose  : create a new stripe with a ChAsymSpine containing the edge <E>
//           with the distance <Dis> and the angle Angle  on the face <F>
//           . Add the tangential edges continuous to E in the spine
//
//=======================================================================

void ChFi3d_ChBuilder::Add(const double       Dis1,
                           const double       Dis2,
                           const TopoDS_Edge& E,
                           const TopoDS_Face& F)
{
  if (!Contains(E) && myEFMap.Contains(E))
  {

    TopoDS_Edge E_wnt = E;
    E_wnt.Orientation(TopAbs_FORWARD);

    occ::handle<ChFiDS_Stripe> Stripe    = new ChFiDS_Stripe();
    occ::handle<ChFiDS_Spine>& Sp        = Stripe->ChangeSpine();
    Sp                                   = new ChFiDS_ChamfSpine(tolesp);
    occ::handle<ChFiDS_ChamfSpine> Spine = occ::down_cast<ChFiDS_ChamfSpine>(Sp);

    Spine->SetMode(myMode);
    double Offset = -1;
    if (myMode == ChFiDS_ConstThroatWithPenetrationChamfer)
    {
      Offset = std::min(Dis1, Dis2);
    }

    Spine->SetEdges(E_wnt);
    if (PerformElement(Spine, Offset, F))
    {
      Spine->Load();
      myListStripe.Append(Stripe);

      Spine->SetDists(Dis1, Dis2);

      PerformExtremity(Spine);
    }
  }
}

//=======================================================================
// function : SetDists
// purpose  : set the distances <Dis1> and <Dis2> of the contour of
//           index IC
//=======================================================================

void ChFi3d_ChBuilder::SetDists(const double       Dis1,
                                const double       Dis2,
                                const int          IC,
                                const TopoDS_Face& F)
{

  if (IC <= NbElements())
  {
    occ::handle<ChFiDS_ChamfSpine> csp = occ::down_cast<ChFiDS_ChamfSpine>(Value(IC));

    // Search the first edge which has a common face equal to F
    TopoDS_Face         F1, F2, FirstF1, FirstF2;
    TopAbs_Orientation  Or1, Or2;
    int                 Choix, ChoixConge;
    BRepAdaptor_Surface Sb1, Sb2;
    int                 i     = 1;
    bool                Found = false;
    while ((i <= csp->NbEdges()) && (!Found))
    {
      SearchCommonFaces(myEFMap, csp->Edges(i), F1, F2);
      if (i == 1)
      {
        FirstF1 = F1;
        FirstF2 = F2;
      }
      Found = (F1.IsSame(F) || F2.IsSame(F));
      i++;
    }

    if (Found)
    {
      if (F2.IsSame(F))
      {
        F2 = F1;
        F1 = F;
      }
      Sb1.Initialize(F1);
      Sb2.Initialize(F2);
      Choix = ChFi3d::ConcaveSide(Sb1, Sb2, csp->Edges(i - 1), Or1, Or2);
      Sb1.Initialize(FirstF1);
      Sb2.Initialize(FirstF2);
      ChoixConge = ChFi3d::ConcaveSide(Sb1, Sb2, csp->Edges(1), Or1, Or2);
      if (ChoixConge % 2 != Choix % 2)
        csp->SetDists(Dis2, Dis1);
      else
        csp->SetDists(Dis1, Dis2);
    }
    else
      throw Standard_DomainError("the face is not common to any of edges of the contour");
  }
}

//=================================================================================================

void ChFi3d_ChBuilder::Dists(const int IC, double& dis1, double& dis2) const
{
  occ::handle<ChFiDS_ChamfSpine> chsp = occ::down_cast<ChFiDS_ChamfSpine>(Value(IC));
  double                         temp1, temp2;
  chsp->Dists(temp1, temp2);
  dis1 = temp1;
  dis2 = temp2;
}

//=======================================================================
// function : Add
// purpose  : create a new stripe with a ChAsymSpine containing the edge <E>
//           with the distance <Dis> and the angle Angle  on the face <F>
//           . Add the tangential edges continuous to E in the spine
//
//=======================================================================

void ChFi3d_ChBuilder::AddDA(const double       Dis1,
                             const double       Angle,
                             const TopoDS_Edge& E,
                             const TopoDS_Face& F)
{
  if (!Contains(E) && myEFMap.Contains(E))
  {

    TopoDS_Edge E_wnt = E;
    E_wnt.Orientation(TopAbs_FORWARD);

    occ::handle<ChFiDS_Stripe> Stripe    = new ChFiDS_Stripe();
    occ::handle<ChFiDS_Spine>& Sp        = Stripe->ChangeSpine();
    Sp                                   = new ChFiDS_ChamfSpine(tolesp);
    occ::handle<ChFiDS_ChamfSpine> Spine = occ::down_cast<ChFiDS_ChamfSpine>(Sp);

    Spine->SetEdges(E_wnt);
    if (PerformElement(Spine, -1, F))
    {
      Spine->Load();
      myListStripe.Append(Stripe);

      Spine->SetDistAngle(Dis1, Angle);

      PerformExtremity(Spine);
    }
  }
}

//=======================================================================
// function : SetDistAngle
// purpose  : set the distances <Dis> and <Angle> of the contour of
//           index IC
//=======================================================================

void ChFi3d_ChBuilder::SetDistAngle(const double       Dis,
                                    const double       Angle,
                                    const int          IC,
                                    const TopoDS_Face& F)
{

  if (IC <= NbElements())
  {
    occ::handle<ChFiDS_ChamfSpine> csp = occ::down_cast<ChFiDS_ChamfSpine>(Value(IC));

    // Search the first edge which has a common face equal to F
    TopoDS_Face         F1, F2, FirstF1, FirstF2;
    BRepAdaptor_Surface Sb1, Sb2;
    int                 i     = 1;
    bool                Found = false;

    while ((i <= csp->NbEdges()) && (!Found))
    {
      SearchCommonFaces(myEFMap, csp->Edges(i), F1, F2);
      if (i == 1)
      {
        FirstF1 = F1;
        FirstF2 = F2;
      }
      Found = (F1.IsSame(F) || F2.IsSame(F));
      i++;
    }

    if (Found)
    {
      if (F2.IsSame(F))
      {
        F2 = F1;
        F1 = F;
      }
      Sb1.Initialize(F1);
      Sb2.Initialize(F2);
      Sb1.Initialize(FirstF1);
      Sb2.Initialize(FirstF2);
      csp->SetDistAngle(Dis, Angle);
    }
    else
      throw Standard_DomainError("the face is not common to any edges of the contour");
  }
}

//=================================================================================================

void ChFi3d_ChBuilder::GetDistAngle(const int IC, double& Dis, double& Angle) const
{
  occ::handle<ChFiDS_ChamfSpine> chsp = occ::down_cast<ChFiDS_ChamfSpine>(Value(IC));

  chsp->GetDistAngle(Dis, Angle);
}

//=======================================================================
// function : SetMode
// purpose  : set the mode of chamfer
//=======================================================================

void ChFi3d_ChBuilder::SetMode(const ChFiDS_ChamfMode theMode)
{
  myMode = theMode;
}

//=================================================================================================

ChFiDS_ChamfMethod ChFi3d_ChBuilder::IsChamfer(const int IC) const
{
  occ::handle<ChFiDS_ChamfSpine> chsp = occ::down_cast<ChFiDS_ChamfSpine>(Value(IC));

  ChFiDS_ChamfMethod ret = chsp->IsChamfer();

  return ret;
}

//=================================================================================================

ChFiDS_ChamfMode ChFi3d_ChBuilder::Mode() const
{
  return myMode;
}

//=================================================================================================

void ChFi3d_ChBuilder::ResetContour(const int IC)
{
  if (IC <= NbElements())
  {
    occ::handle<ChFiDS_ChamfSpine> chsp = occ::down_cast<ChFiDS_ChamfSpine>(Value(IC));
    chsp->Reset(true);
  }
}

//--------------------------------AJOUT----------------------------------
//=================================================================================================

void ChFi3d_ChBuilder::Simulate(const int IC)
{
#ifdef OCCT_DEBUG
  if (ChFi3d_GettraceCHRON())
  {
    simul.Reset();
    elspine.Reset();
    chemine.Reset();
    simul.Start();
  }
#endif
  NCollection_List<occ::handle<ChFiDS_Stripe>>::Iterator itel;
  int                                                    i = 1;
  for (itel.Initialize(myListStripe); itel.More(); itel.Next(), i++)
  {
    if (i == IC)
    {
      PerformSetOfSurf(itel.ChangeValue(), true);
      break;
    }
  }
#ifdef OCCT_DEBUG
  if (ChFi3d_GettraceCHRON())
  {
    simul.Stop();
    std::cout << "Total simulation time : ";
    simul.Show();
    std::cout << "Spine construction time : ";
    elspine.Show();
    std::cout << "and progression time : ";
    chemine.Show();
  }
#endif
}

//---------------------------AJOUT---------------------------------------
//=================================================================================================

int ChFi3d_ChBuilder::NbSurf(const int IC) const
{
  NCollection_List<occ::handle<ChFiDS_Stripe>>::Iterator itel;
  int                                                    i = 1;
  for (itel.Initialize(myListStripe); itel.More(); itel.Next(), i++)
  {
    if (i == IC)
    {
      return itel.Value()->SetOfSurfData()->Length();
    }
  }
  return 0;
}

//--------------------------AJOUT---------------------------------------
//=================================================================================================

occ::handle<NCollection_HArray1<ChFiDS_CircSection>> ChFi3d_ChBuilder::Sect(const int IC,
                                                                            const int IS) const
{
  NCollection_List<occ::handle<ChFiDS_Stripe>>::Iterator itel;
  int                                                    i = 1;
  occ::handle<NCollection_HArray1<ChFiDS_CircSection>>   res;
  for (itel.Initialize(myListStripe); itel.More(); itel.Next(), i++)
  {
    if (i == IC)
    {
      occ::handle<Standard_Transient> bid = itel.Value()->SetOfSurfData()->Value(IS)->Simul();
      res = occ::down_cast<NCollection_HArray1<ChFiDS_CircSection>>(bid);
      return res;
    }
  }
  return occ::handle<NCollection_HArray1<ChFiDS_CircSection>>();
}

//-------------------MODIFS---------------------------------------------
//=======================================================================
// function : SimulKPart
// purpose  : Stores simulating sections in simul
//=======================================================================

void ChFi3d_ChBuilder::SimulKPart(const occ::handle<ChFiDS_SurfData>& SD) const
{
  TopOpeBRepDS_DataStructure& DStr = myDS->ChangeDS();
  occ::handle<Geom_Surface>   S    = DStr.Surface(SD->Surf()).Surface();
  gp_Pnt2d                    p1f =
    SD->InterferenceOnS1().PCurveOnSurf()->Value(SD->InterferenceOnS1().FirstParameter());
  gp_Pnt2d p1l =
    SD->InterferenceOnS1().PCurveOnSurf()->Value(SD->InterferenceOnS1().LastParameter());
  gp_Pnt2d p2f =
    SD->InterferenceOnS2().PCurveOnSurf()->Value(SD->InterferenceOnS2().FirstParameter());
  gp_Pnt2d p2l =
    SD->InterferenceOnS2().PCurveOnSurf()->Value(SD->InterferenceOnS2().LastParameter());
  GeomAdaptor_Surface                                  AS(S);
  occ::handle<NCollection_HArray1<ChFiDS_CircSection>> sec;
  double                                               u1, v1, u2, v2;
  GeomAbs_SurfaceType                                  typ = AS.GetType();
  switch (typ)
  {
    case GeomAbs_Plane: {
      v1                       = p1f.Y();
      v2                       = p2f.Y();
      u1                       = std::max(p1f.X(), p2f.X());
      u2                       = std::min(p1l.X(), p2l.X());
      sec                      = new NCollection_HArray1<ChFiDS_CircSection>(1, 2);
      gp_Pln              Pl   = AS.Plane();
      ChFiDS_CircSection& sec1 = sec->ChangeValue(1);
      ChFiDS_CircSection& sec2 = sec->ChangeValue(2);
      sec1.Set(ElSLib::PlaneUIso(Pl.Position(), u1), v1, v2);
      sec2.Set(ElSLib::PlaneUIso(Pl.Position(), u2), v1, v2);
    }
    break;
    case GeomAbs_Cone: {
      v1          = p1f.Y();
      v2          = p2f.Y();
      u1          = std::max(p1f.X(), p2f.X());
      u2          = std::min(p1l.X(), p2l.X());
      double  ang = (u2 - u1);
      gp_Cone Co  = AS.Cone();
      double  rad = Co.RefRadius(), sang = Co.SemiAngle();
      int     n = (int)(36. * ang / M_PI + 1);
      if (n < 2)
        n = 2;
      sec = new NCollection_HArray1<ChFiDS_CircSection>(1, n);
      for (int i = 1; i <= n; i++)
      {
        ChFiDS_CircSection& isec = sec->ChangeValue(i);
        double              u    = u1 + (i - 1) * (u2 - u1) / (n - 1);
        isec.Set(ElSLib::ConeUIso(Co.Position(), rad, sang, u), v1, v2);
      }
    }
    break;
    default:
      break;
  }
  SD->SetSimul(sec);
}

//------------------------MODIFS---------------------------------------
//=================================================================================================

bool ChFi3d_ChBuilder::SimulSurf(occ::handle<ChFiDS_SurfData>&           Data,
                                 const occ::handle<ChFiDS_ElSpine>&      HGuide,
                                 const occ::handle<ChFiDS_Spine>&        Spine,
                                 const int                               Choix,
                                 const occ::handle<BRepAdaptor_Surface>& S1,
                                 const occ::handle<Adaptor3d_TopolTool>& I1,
                                 const occ::handle<BRepAdaptor_Surface>& S2,
                                 const occ::handle<Adaptor3d_TopolTool>& I2,
                                 const double                            TolGuide,
                                 double&                                 First,
                                 double&                                 Last,
                                 const bool                              Inside,
                                 const bool                              Appro,
                                 const bool                              Forward,
                                 const bool                              RecOnS1,
                                 const bool                              RecOnS2,
                                 const math_Vector&                      Soldep,
                                 int&                                    intf,
                                 int&                                    intl)

{
  occ::handle<ChFiDS_ChamfSpine> chsp = occ::down_cast<ChFiDS_ChamfSpine>(Spine);

  if (chsp.IsNull())
    throw Standard_ConstructionError("SimulSurf : this is not the spine of a chamfer");

  double radius;
  // Flexible parameters!
  double la = HGuide->LastParameter(), fi = HGuide->FirstParameter();
  double longueur    = la - fi;
  double MaxStep     = longueur * 0.05;
  double radiusspine = 0, locfleche, w;
  gp_Pnt Pbid;
  gp_Vec d1, d2;
  // As ElSpine is parameterized by a curvilinear quasi-abscissa,
  // the min radius is estimated as 1/D2 max;
  // for(int i = 0; i <= 20; i++){
  int i;
  for (i = 0; i <= 20; i++)
  {
    w = fi + i * MaxStep;
    HGuide->D2(w, Pbid, d1, d2);
    double temp = d2.SquareMagnitude();
    if (temp > radiusspine)
      radiusspine = temp;
  }

  occ::handle<BRepBlend_Line> lin;
  double                      PFirst = First;
  if (intf)
    First = chsp->FirstParameter(1);
  if (intl)
    Last = chsp->LastParameter(chsp->NbEdges());

  occ::handle<ChFiDS_ElSpine> OffsetHGuide;

  if (chsp->IsChamfer() == ChFiDS_Sym)
  {
    double dis;
    chsp->GetDist(dis);
    radius    = std::max(dis, radiusspine);
    locfleche = radius * 1.e-2; // graphic criterion

    std::unique_ptr<BlendFunc_GenChamfer>  pFunc;
    std::unique_ptr<BlendFunc_GenChamfInv> pFInv;
    if (chsp->Mode() == ChFiDS_ClassicChamfer)
    {
      pFunc.reset(new BRepBlend_Chamfer(S1, S2, HGuide));
      pFInv.reset(new BRepBlend_ChamfInv(S1, S2, HGuide));
    }
    else
    {
      pFunc.reset(new BRepBlend_ConstThroat(S1, S2, HGuide));
      pFInv.reset(new BRepBlend_ConstThroatInv(S1, S2, HGuide));
    }
    pFunc->Set(dis, dis, Choix);
    pFInv->Set(dis, dis, Choix);

    done = SimulData(Data,
                     HGuide,
                     OffsetHGuide,
                     lin,
                     S1,
                     I1,
                     S2,
                     I2,
                     *pFunc,
                     *pFInv,
                     PFirst,
                     MaxStep,
                     locfleche,
                     TolGuide,
                     First,
                     Last,
                     Inside,
                     Appro,
                     Forward,
                     Soldep,
                     4,
                     RecOnS1,
                     RecOnS2);

    if (!done)
      return false;
    occ::handle<NCollection_HArray1<ChFiDS_CircSection>> sec;
    gp_Pnt2d                                             pf1, pl1, pf2, pl2;

    int nbp = lin->NbPoints();
    sec     = new NCollection_HArray1<ChFiDS_CircSection>(1, nbp);
    for (i = 1; i <= nbp; i++)
    {
      ChFiDS_CircSection& isec = sec->ChangeValue(i);
      double              u1, v1, u2, v2, ww, p1, p2;
      gp_Lin              line;
      const Blend_Point&  p = lin->Point(i);
      p.ParametersOnS1(u1, v1);
      p.ParametersOnS2(u2, v2);
      ww = p.Parameter();
      pFunc->Section(ww, u1, v1, u2, v2, p1, p2, line);
      isec.Set(line, p1, p2);
      if (i == 1)
      {
        pf1.SetCoord(u1, v1);
        pf2.SetCoord(u2, v2);
      }
      if (i == nbp)
      {
        pl1.SetCoord(u1, v1);
        pl2.SetCoord(u2, v2);
      }
    }

    Data->SetSimul(sec);
    Data->Set2dPoints(pf1, pl1, pf2, pl2);
    ChFi3d_FilCommonPoint(lin->StartPointOnFirst(),
                          lin->TransitionOnS1(),
                          true,
                          Data->ChangeVertexFirstOnS1(),
                          tolapp3d);
    ChFi3d_FilCommonPoint(lin->EndPointOnFirst(),
                          lin->TransitionOnS1(),
                          false,
                          Data->ChangeVertexLastOnS1(),
                          tolapp3d);
    ChFi3d_FilCommonPoint(lin->StartPointOnSecond(),
                          lin->TransitionOnS2(),
                          true,
                          Data->ChangeVertexFirstOnS2(),
                          tolapp3d);
    ChFi3d_FilCommonPoint(lin->EndPointOnSecond(),
                          lin->TransitionOnS2(),
                          false,
                          Data->ChangeVertexLastOnS2(),
                          tolapp3d);

    bool reverse = (!Forward || Inside);
    if (intf && reverse)
    {
      bool                      ok  = false;
      const ChFiDS_CommonPoint& cp1 = Data->VertexFirstOnS1();
      if (cp1.IsOnArc())
      {
        TopoDS_Face F1 = S1->Face();
        TopoDS_Face bid;
        intf = !SearchFace(Spine, cp1, F1, bid);
        ok   = intf != 0;
      }
      const ChFiDS_CommonPoint& cp2 = Data->VertexFirstOnS2();
      if (cp2.IsOnArc() && !ok)
      {
        TopoDS_Face F2 = S2->Face();
        TopoDS_Face bid;
        intf = !SearchFace(Spine, cp2, F2, bid);
      }
    }
    if (intl)
    {
      bool                      ok  = false;
      const ChFiDS_CommonPoint& cp1 = Data->VertexLastOnS1();
      if (cp1.IsOnArc())
      {
        TopoDS_Face F1 = S1->Face();
        TopoDS_Face bid;
        intl = !SearchFace(Spine, cp1, F1, bid);
        ok   = intl != 0;
      }
      const ChFiDS_CommonPoint& cp2 = Data->VertexLastOnS2();
      if (cp2.IsOnArc() && !ok)
      {
        TopoDS_Face F2 = S2->Face();
        TopoDS_Face bid;
        intl = !SearchFace(Spine, cp2, F2, bid);
      }
    }
  }
  else if (chsp->IsChamfer() == ChFiDS_TwoDist)
  {
    double dis1, dis2;
    chsp->Dists(dis1, dis2);
    radius    = std::max(dis1, dis2);
    radius    = std::max(radius, radiusspine);
    locfleche = radius * 1.e-2; // graphic criterion

    std::unique_ptr<BlendFunc_GenChamfer>  pFunc;
    std::unique_ptr<BlendFunc_GenChamfInv> pFInv;
    if (chsp->Mode() == ChFiDS_ClassicChamfer)
    {
      pFunc.reset(new BRepBlend_Chamfer(S1, S2, HGuide));
      pFInv.reset(new BRepBlend_ChamfInv(S1, S2, HGuide));
      pFunc->Set(dis1, dis2, Choix);
      pFInv->Set(dis1, dis2, Choix);
    }
    else
    {
      NCollection_List<occ::handle<ChFiDS_ElSpine>>& ll        = Spine->ChangeElSpines();
      NCollection_List<occ::handle<ChFiDS_ElSpine>>& ll_offset = Spine->ChangeOffsetElSpines();
      NCollection_List<occ::handle<ChFiDS_ElSpine>>::Iterator ILES(ll), ILES_offset(ll_offset);
      for (; ILES.More(); ILES.Next(), ILES_offset.Next())
      {
        const occ::handle<ChFiDS_ElSpine>& aHElSpine = ILES.Value();
        if (aHElSpine == HGuide)
          OffsetHGuide = ILES_offset.Value();
      }

      if (OffsetHGuide.IsNull())
      {
        std::cout << std::endl << "Construction of offset guide failed!" << std::endl;
        // exception
      }
      pFunc.reset(new BRepBlend_ConstThroatWithPenetration(S1, S2, OffsetHGuide));
      pFInv.reset(new BRepBlend_ConstThroatWithPenetrationInv(S1, S2, OffsetHGuide));
      double Throat = std::max(dis1, dis2);
      pFunc->Set(Throat, Throat, Choix);
      pFInv->Set(Throat, Throat, Choix);
    }

    done = SimulData(Data,
                     HGuide,
                     OffsetHGuide,
                     lin,
                     S1,
                     I1,
                     S2,
                     I2,
                     *pFunc,
                     *pFInv,
                     PFirst,
                     MaxStep,
                     locfleche,
                     TolGuide,
                     First,
                     Last,
                     Inside,
                     Appro,
                     Forward,
                     Soldep,
                     4,
                     RecOnS1,
                     RecOnS2);

    if (!done)
      return false;
    occ::handle<NCollection_HArray1<ChFiDS_CircSection>> sec;
    gp_Pnt2d                                             pf1, pl1, pf2, pl2;

    int nbp = lin->NbPoints();
    sec     = new NCollection_HArray1<ChFiDS_CircSection>(1, nbp);
    for (i = 1; i <= nbp; i++)
    {
      ChFiDS_CircSection& isec = sec->ChangeValue(i);
      double              u1, v1, u2, v2, ww, p1, p2;
      gp_Lin              line;
      const Blend_Point&  p = lin->Point(i);
      p.ParametersOnS1(u1, v1);
      p.ParametersOnS2(u2, v2);
      ww = p.Parameter();
      pFunc->Section(ww, u1, v1, u2, v2, p1, p2, line);
      isec.Set(line, p1, p2);
      if (i == 1)
      {
        pf1.SetCoord(u1, v1);
        pf2.SetCoord(u2, v2);
      }
      if (i == nbp)
      {
        pl1.SetCoord(u1, v1);
        pl2.SetCoord(u2, v2);
      }
    }

    Data->SetSimul(sec);
    Data->Set2dPoints(pf1, pl1, pf2, pl2);
    ChFi3d_FilCommonPoint(lin->StartPointOnFirst(),
                          lin->TransitionOnS1(),
                          true,
                          Data->ChangeVertexFirstOnS1(),
                          tolapp3d);
    ChFi3d_FilCommonPoint(lin->EndPointOnFirst(),
                          lin->TransitionOnS1(),
                          false,
                          Data->ChangeVertexLastOnS1(),
                          tolapp3d);
    ChFi3d_FilCommonPoint(lin->StartPointOnSecond(),
                          lin->TransitionOnS2(),
                          true,
                          Data->ChangeVertexFirstOnS2(),
                          tolapp3d);
    ChFi3d_FilCommonPoint(lin->EndPointOnSecond(),
                          lin->TransitionOnS2(),
                          false,
                          Data->ChangeVertexLastOnS2(),
                          tolapp3d);

    bool reverse = (!Forward || Inside);
    if (intf && reverse)
    {
      bool                      ok  = false;
      const ChFiDS_CommonPoint& cp1 = Data->VertexFirstOnS1();
      if (cp1.IsOnArc())
      {
        TopoDS_Face F1 = S1->Face();
        TopoDS_Face bid;
        intf = !SearchFace(Spine, cp1, F1, bid);
        ok   = intf != 0;
      }
      const ChFiDS_CommonPoint& cp2 = Data->VertexFirstOnS2();
      if (cp2.IsOnArc() && !ok)
      {
        TopoDS_Face F2 = S2->Face();
        TopoDS_Face bid;
        intf = !SearchFace(Spine, cp2, F2, bid);
      }
    }
    if (intl)
    {
      bool                      ok  = false;
      const ChFiDS_CommonPoint& cp1 = Data->VertexLastOnS1();
      if (cp1.IsOnArc())
      {
        TopoDS_Face F1 = S1->Face();
        TopoDS_Face bid;
        intl = !SearchFace(Spine, cp1, F1, bid);
        ok   = intl != 0;
      }
      const ChFiDS_CommonPoint& cp2 = Data->VertexLastOnS2();
      if (cp2.IsOnArc() && !ok)
      {
        TopoDS_Face F2 = S2->Face();
        TopoDS_Face bid;
        intl = !SearchFace(Spine, cp2, F2, bid);
      }
    }
  }
  else
  { // distance and angle
    double dis, angle;
    chsp->GetDistAngle(dis, angle);
    radius    = std::max(dis, dis * tan(angle));
    radius    = std::max(radius, radiusspine);
    locfleche = radius * 1.e-2; // graphic criterion

    int Ch = Choix;

    BRepBlend_ChAsym    Func(S1, S2, HGuide);
    BRepBlend_ChAsymInv FInv(S1, S2, HGuide);

    Func.Set(dis, angle, Ch);
    FInv.Set(dis, angle, Ch);

    done = SimulData(Data,
                     HGuide,
                     OffsetHGuide,
                     lin,
                     S1,
                     I1,
                     S2,
                     I2,
                     Func,
                     FInv,
                     PFirst,
                     MaxStep,
                     locfleche,
                     TolGuide,
                     First,
                     Last,
                     Inside,
                     Appro,
                     Forward,
                     Soldep,
                     4,
                     RecOnS1,
                     RecOnS2);

    if (!done)
      return false;
    occ::handle<NCollection_HArray1<ChFiDS_CircSection>> sec;
    gp_Pnt2d                                             pf1, pl1, pf2, pl2;

    int nbp = lin->NbPoints();
    sec     = new NCollection_HArray1<ChFiDS_CircSection>(1, nbp);
    for (i = 1; i <= nbp; i++)
    {
      ChFiDS_CircSection& isec = sec->ChangeValue(i);
      double              u1, v1, u2, v2, ww, p1, p2;
      gp_Lin              line;
      const Blend_Point&  p = lin->Point(i);
      p.ParametersOnS1(u1, v1);
      p.ParametersOnS2(u2, v2);
      ww = p.Parameter();
      Func.Section(ww, u1, v1, u2, v2, p1, p2, line);
      isec.Set(line, p1, p2);
      if (i == 1)
      {
        pf1.SetCoord(u1, v1);
        pf2.SetCoord(u2, v2);
      }
      if (i == nbp)
      {
        pl1.SetCoord(u1, v1);
        pl2.SetCoord(u2, v2);
      }
    }

    Data->SetSimul(sec);
    Data->Set2dPoints(pf1, pl1, pf2, pl2);
    ChFi3d_FilCommonPoint(lin->StartPointOnFirst(),
                          lin->TransitionOnS1(),
                          true,
                          Data->ChangeVertexFirstOnS1(),
                          tolapp3d);
    ChFi3d_FilCommonPoint(lin->EndPointOnFirst(),
                          lin->TransitionOnS1(),
                          false,
                          Data->ChangeVertexLastOnS1(),
                          tolapp3d);
    ChFi3d_FilCommonPoint(lin->StartPointOnSecond(),
                          lin->TransitionOnS2(),
                          true,
                          Data->ChangeVertexFirstOnS2(),
                          tolapp3d);
    ChFi3d_FilCommonPoint(lin->EndPointOnSecond(),
                          lin->TransitionOnS2(),
                          false,
                          Data->ChangeVertexLastOnS2(),
                          tolapp3d);

    bool reverse = (!Forward || Inside);
    if (intf && reverse)
    {
      bool                      ok  = false;
      const ChFiDS_CommonPoint& cp1 = Data->VertexFirstOnS1();
      if (cp1.IsOnArc())
      {
        TopoDS_Face F1 = S1->Face();
        TopoDS_Face bid;
        intf = !SearchFace(Spine, cp1, F1, bid);
        ok   = intf != 0;
      }
      const ChFiDS_CommonPoint& cp2 = Data->VertexFirstOnS2();
      if (cp2.IsOnArc() && !ok)
      {
        TopoDS_Face F2 = S2->Face();
        TopoDS_Face bid;
        intf = !SearchFace(Spine, cp2, F2, bid);
      }
    }

    if (intl)
    {
      bool                      ok  = false;
      const ChFiDS_CommonPoint& cp1 = Data->VertexLastOnS1();
      if (cp1.IsOnArc())
      {
        TopoDS_Face F1 = S1->Face();
        TopoDS_Face bid;
        intl = !SearchFace(Spine, cp1, F1, bid);
        ok   = intl != 0;
      }
      const ChFiDS_CommonPoint& cp2 = Data->VertexLastOnS2();
      if (cp2.IsOnArc() && !ok)
      {
        TopoDS_Face F2 = S2->Face();
        TopoDS_Face bid;
        intl = !SearchFace(Spine, cp2, F2, bid);
      }
    }
  } // distance and angle
  return true;
}

void ChFi3d_ChBuilder::SimulSurf(occ::handle<ChFiDS_SurfData>&,
                                 const occ::handle<ChFiDS_ElSpine>&,
                                 const occ::handle<ChFiDS_Spine>&,
                                 const int,
                                 const occ::handle<BRepAdaptor_Surface>&,
                                 const occ::handle<Adaptor3d_TopolTool>&,
                                 const occ::handle<BRepAdaptor_Curve2d>&,
                                 const occ::handle<BRepAdaptor_Surface>&,
                                 const occ::handle<BRepAdaptor_Curve2d>&,
                                 bool&,
                                 const occ::handle<BRepAdaptor_Surface>&,
                                 const occ::handle<Adaptor3d_TopolTool>&,
                                 const TopAbs_Orientation,
                                 const double,
                                 const double,
                                 double&,
                                 double&,
                                 const bool,
                                 const bool,
                                 const bool,
                                 const bool,
                                 const bool,
                                 const bool,
                                 const math_Vector&)
{
  throw Standard_Failure("SimulSurf Not Implemented");
}

void ChFi3d_ChBuilder::SimulSurf(occ::handle<ChFiDS_SurfData>&,
                                 const occ::handle<ChFiDS_ElSpine>&,
                                 const occ::handle<ChFiDS_Spine>&,
                                 const int,
                                 const occ::handle<BRepAdaptor_Surface>&,
                                 const occ::handle<Adaptor3d_TopolTool>&,
                                 const TopAbs_Orientation,
                                 const occ::handle<BRepAdaptor_Surface>&,
                                 const occ::handle<Adaptor3d_TopolTool>&,
                                 const occ::handle<BRepAdaptor_Curve2d>&,
                                 const occ::handle<BRepAdaptor_Surface>&,
                                 const occ::handle<BRepAdaptor_Curve2d>&,
                                 bool&,
                                 const double,
                                 const double,
                                 double&,
                                 double&,
                                 const bool,
                                 const bool,
                                 const bool,
                                 const bool,
                                 const bool,
                                 const bool,
                                 const math_Vector&)
{
  throw Standard_Failure("SimulSurf Not Implemented");
}

void ChFi3d_ChBuilder::SimulSurf(occ::handle<ChFiDS_SurfData>&,
                                 const occ::handle<ChFiDS_ElSpine>&,
                                 const occ::handle<ChFiDS_Spine>&,
                                 const int,
                                 const occ::handle<BRepAdaptor_Surface>&,
                                 const occ::handle<Adaptor3d_TopolTool>&,
                                 const occ::handle<BRepAdaptor_Curve2d>&,
                                 const occ::handle<BRepAdaptor_Surface>&,
                                 const occ::handle<BRepAdaptor_Curve2d>&,
                                 bool&,
                                 const TopAbs_Orientation,
                                 const occ::handle<BRepAdaptor_Surface>&,
                                 const occ::handle<Adaptor3d_TopolTool>&,
                                 const occ::handle<BRepAdaptor_Curve2d>&,
                                 const occ::handle<BRepAdaptor_Surface>&,
                                 const occ::handle<BRepAdaptor_Curve2d>&,
                                 bool&,
                                 const TopAbs_Orientation,
                                 const double,
                                 const double,
                                 double&,
                                 double&,
                                 const bool,
                                 const bool,
                                 const bool,
                                 const bool,
                                 const bool,
                                 const bool,
                                 const bool,
                                 const math_Vector&)
{
  throw Standard_Failure("SimulSurf Not Implemented");
}

//------------------------MODIFS---------------------------------------
//=======================================================================
// function : PerformFirstSection
// purpose  : to implement the first section if there is no KPart
//=======================================================================

bool ChFi3d_ChBuilder::PerformFirstSection(const occ::handle<ChFiDS_Spine>&        Spine,
                                           const occ::handle<ChFiDS_ElSpine>&      HGuide,
                                           const int                               Choix,
                                           occ::handle<BRepAdaptor_Surface>&       S1,
                                           occ::handle<BRepAdaptor_Surface>&       S2,
                                           const occ::handle<Adaptor3d_TopolTool>& I1,
                                           const occ::handle<Adaptor3d_TopolTool>& I2,
                                           const double                            Par,
                                           math_Vector&                            SolDep,
                                           TopAbs_State&                           Pos1,
                                           TopAbs_State&                           Pos2) const
{
  occ::handle<ChFiDS_ChamfSpine> chsp = occ::down_cast<ChFiDS_ChamfSpine>(Spine);

  if (chsp.IsNull())
    throw Standard_ConstructionError("PerformSurf : this is not the spine of a chamfer");

  double TolGuide = HGuide->Resolution(tolapp3d);

  if (chsp->IsChamfer() == ChFiDS_Sym)
  {
    double dis;
    chsp->GetDist(dis);

    std::unique_ptr<BlendFunc_GenChamfer> pFunc;
    if (chsp->Mode() == ChFiDS_ClassicChamfer)
    {
      pFunc.reset(new BRepBlend_Chamfer(S1, S2, HGuide));
    }
    else
    {
      pFunc.reset(new BRepBlend_ConstThroat(S1, S2, HGuide));
    }
    pFunc->Set(dis, dis, Choix);
    BRepBlend_Walking TheWalk(S1, S2, I1, I2, HGuide);

    // calculate an approximate starting solution
    gp_Vec TgF, TgL, tmp1, tmp2, d1gui;
    gp_Pnt pt1, pt2, ptgui;
    gp_XYZ temp;

    HGuide->D1(Par, ptgui, d1gui);
    //  ptgui = S1->Value(SolDep(1),SolDep(2));

    pFunc->Set(Par);
    pFunc->Tangent(SolDep(1), SolDep(2), SolDep(3), SolDep(4), TgF, TgL, tmp1, tmp2);

    bool   rev1 = false;
    bool   rev2 = false;
    double sign = (TgF.Crossed(d1gui)).Dot(TgL);

    if (Choix % 2 == 1)
      rev1 = true;
    else
      rev2 = true;

    if (sign < 0.)
    {
      rev1 = !rev1;
      rev2 = !rev2;
    }

    if (rev1)
      TgF.Reverse();
    if (rev2)
      TgL.Reverse();

    temp = (TgF.XYZ()).Multiplied(dis);
    pt1.SetXYZ((ptgui.XYZ()).Added(temp));
    temp = (TgL.XYZ()).Multiplied(dis);
    pt2.SetXYZ((ptgui.XYZ()).Added(temp));

    double tol = tolesp * 1.e2;
    //    double u,v;
    Extrema_GenLocateExtPS proj1(*S1, tol, tol);
    proj1.Perform(pt1, SolDep(1), SolDep(2));
    Extrema_GenLocateExtPS proj2(*S2, tol, tol);
    proj2.Perform(pt2, SolDep(3), SolDep(4));

    if (proj1.IsDone())
    {
      (proj1.Point()).Parameter(SolDep(1), SolDep(2));
    }
    if (proj2.IsDone())
    {
      (proj2.Point()).Parameter(SolDep(3), SolDep(4));
    }

    return TheWalk.PerformFirstSection(*pFunc, Par, SolDep, tolapp3d, TolGuide, Pos1, Pos2);
  }
  else if (chsp->IsChamfer() == ChFiDS_TwoDist)
  {
    double dis1, dis2;
    chsp->Dists(dis1, dis2);

    std::unique_ptr<BlendFunc_GenChamfer> pFunc;
    if (chsp->Mode() == ChFiDS_ClassicChamfer)
    {
      pFunc.reset(new BRepBlend_Chamfer(S1, S2, HGuide));
      pFunc->Set(dis1, dis2, Choix);
    }
    else
    {
      occ::handle<ChFiDS_ElSpine>                    OffsetHGuide;
      NCollection_List<occ::handle<ChFiDS_ElSpine>>& ll        = Spine->ChangeElSpines();
      NCollection_List<occ::handle<ChFiDS_ElSpine>>& ll_offset = Spine->ChangeOffsetElSpines();
      NCollection_List<occ::handle<ChFiDS_ElSpine>>::Iterator ILES(ll), ILES_offset(ll_offset);
      for (; ILES.More(); ILES.Next(), ILES_offset.Next())
      {
        const occ::handle<ChFiDS_ElSpine>& aHElSpine = ILES.Value();
        if (aHElSpine == HGuide)
          OffsetHGuide = ILES_offset.Value();
      }

      if (OffsetHGuide.IsNull())
      {
        std::cout << std::endl << "Construction of offset guide failed!" << std::endl;
        // exception
      }
      pFunc.reset(new BRepBlend_ConstThroatWithPenetration(S1, S2, OffsetHGuide));
      double Throat = std::max(dis1, dis2);
      pFunc->Set(Throat, Throat, Choix); // dis2?
    }
    BRepBlend_Walking TheWalk(S1, S2, I1, I2, HGuide);

    // calculate an approximate starting solution
    gp_Vec TgF, TgL, tmp1, tmp2, d1gui;
    gp_Pnt pt1, pt2, ptgui;
    gp_XYZ temp;

    HGuide->D1(Par, ptgui, d1gui);
    //  ptgui = S1->Value(SolDep(1),SolDep(2));

    pFunc->Set(Par);
    pFunc->Tangent(SolDep(1), SolDep(2), SolDep(3), SolDep(4), TgF, TgL, tmp1, tmp2);

    bool   rev1 = false;
    bool   rev2 = false;
    double sign = (TgF.Crossed(d1gui)).Dot(TgL);

    if (Choix % 2 == 1)
      rev1 = true;
    else
      rev2 = true;

    if (sign < 0.)
    {
      rev1 = !rev1;
      rev2 = !rev2;
    }

    if (rev1)
      TgF.Reverse();
    if (rev2)
      TgL.Reverse();

    double aDist1 = dis1, aDist2 = dis2;
    if (chsp->Mode() == ChFiDS_ConstThroatWithPenetrationChamfer)
    {
      /*
      double Alpha = TgF.Angle(TgL);
      double SinAlpha = std::sin(Alpha);
      double CosAlpha = std::cos(Alpha);
      double TanAlpha = std::tan(Alpha);
      double dis1dis1 = dis1*dis1, dis2dis2 = dis2*dis2;
      aDist2 = sqrt(dis1dis1 - dis2dis2) - dis2/TanAlpha;
      double CosBeta = sqrt(1-dis2dis2/dis1dis1)*CosAlpha + dis2/dis1*SinAlpha;
      double FullDist1 = dis1/CosBeta;
      aDist1 = FullDist1 - dis2/SinAlpha;
      */
    }

    temp = (TgF.XYZ()).Multiplied(aDist1);
    pt1.SetXYZ((ptgui.XYZ()).Added(temp));
    temp = (TgL.XYZ()).Multiplied(aDist2);
    pt2.SetXYZ((ptgui.XYZ()).Added(temp));

    double tol = tolesp * 1.e2;
    //    double u,v;

    Extrema_GenLocateExtPS proj1(*S1, tol, tol);
    proj1.Perform(pt1, SolDep(1), SolDep(2));
    Extrema_GenLocateExtPS proj2(*S2, tol, tol);
    proj2.Perform(pt2, SolDep(3), SolDep(4));

    if (proj1.IsDone())
    {
      (proj1.Point()).Parameter(SolDep(1), SolDep(2));
    }
    if (proj2.IsDone())
    {
      (proj2.Point()).Parameter(SolDep(3), SolDep(4));
    }

    return TheWalk.PerformFirstSection(*pFunc, Par, SolDep, tolapp3d, TolGuide, Pos1, Pos2);
  }
  else
  { // distance and angle
    double dis1, angle;
    chsp->GetDistAngle(dis1, angle);

    int Ch = Choix;

    BRepBlend_ChAsym Func(S1, S2, HGuide);
    Func.Set(dis1, angle, Ch);
    BRepBlend_Walking TheWalk(S1, S2, I1, I2, HGuide);

    // calculate an approximate starting solution
    gp_Vec TgF, TgL, tmp1, tmp2, d1gui;
    gp_Pnt pt1, pt2, ptgui;
    gp_XYZ temp;

    HGuide->D1(Par, ptgui, d1gui);
    //  ptgui = S1->Value(SolDep(1),SolDep(2));

    Func.Set(Par);
    Func.Tangent(SolDep(1), SolDep(2), SolDep(3), SolDep(4), TgF, TgL, tmp1, tmp2);

    bool   rev1 = false;
    bool   rev2 = false;
    double sign = (TgF.Crossed(d1gui)).Dot(TgL);

    if (Ch % 2 == 1)
      rev1 = true;
    else
      rev2 = true;

    if (sign < 0.)
    {
      rev1 = !rev1;
      rev2 = !rev2;
    }

    if (rev1)
      TgF.Reverse();
    if (rev2)
      TgL.Reverse();

    temp = (TgF.XYZ()).Multiplied(dis1);
    pt1.SetXYZ((ptgui.XYZ()).Added(temp));

    double dis2, tmpcos, tmpsin;
    tmpcos = TgF.Dot(TgL);
    tmpsin = sqrt(1. - tmpcos * tmpcos);

    dis2 = dis1 / (tmpcos + tmpsin / tan(angle));

    temp = (TgL.XYZ()).Multiplied(dis2);
    pt2.SetXYZ((ptgui.XYZ()).Added(temp));

    double tol = tolesp * 1.e2;
    //      double u,v;
    Extrema_GenLocateExtPS proj1(*S1, tol, tol);
    proj1.Perform(pt1, SolDep(1), SolDep(2));
    Extrema_GenLocateExtPS proj2(*S2, tol, tol);
    proj2.Perform(pt2, SolDep(3), SolDep(4));
    if (proj1.IsDone())
    {
      (proj1.Point()).Parameter(SolDep(1), SolDep(2));
    }
    if (proj2.IsDone())
    {
      (proj2.Point()).Parameter(SolDep(3), SolDep(4));
    }

    return TheWalk.PerformFirstSection(Func, Par, SolDep, tolapp3d, TolGuide, Pos1, Pos2);
  } // distance and angle
}

//=================================================================================================

bool ChFi3d_ChBuilder::PerformSurf(NCollection_Sequence<occ::handle<ChFiDS_SurfData>>& SeqData,
                                   const occ::handle<ChFiDS_ElSpine>&                  HGuide,
                                   const occ::handle<ChFiDS_Spine>&                    Spine,
                                   const int                                           Choix,
                                   const occ::handle<BRepAdaptor_Surface>&             S1,
                                   const occ::handle<Adaptor3d_TopolTool>&             I1,
                                   const occ::handle<BRepAdaptor_Surface>&             S2,
                                   const occ::handle<Adaptor3d_TopolTool>&             I2,
                                   const double                                        MaxStep,
                                   const double                                        Fleche,
                                   const double                                        TolGuide,
                                   double&                                             First,
                                   double&                                             Last,
                                   const bool                                          Inside,
                                   const bool                                          Appro,
                                   const bool                                          Forward,
                                   const bool                                          RecOnS1,
                                   const bool                                          RecOnS2,
                                   const math_Vector&                                  Soldep,
                                   int&                                                intf,
                                   int&                                                intl)

{
  occ::handle<ChFiDS_SurfData>   Data = SeqData(1);
  occ::handle<ChFiDS_ChamfSpine> chsp = occ::down_cast<ChFiDS_ChamfSpine>(Spine);

  if (chsp.IsNull())
    throw Standard_ConstructionError("PerformSurf : this is not the spine of a chamfer");

  bool                        gd1, gd2, gf1, gf2;
  occ::handle<BRepBlend_Line> lin;
  TopAbs_Orientation          Or     = S1->Face().Orientation();
  double                      PFirst = First;
  if (intf)
    First = chsp->FirstParameter(1);
  if (intl)
    Last = chsp->LastParameter(chsp->NbEdges());

  if (chsp->IsChamfer() == ChFiDS_Sym)
  {

    std::unique_ptr<BlendFunc_GenChamfer>  pFunc;
    std::unique_ptr<BlendFunc_GenChamfInv> pFInv;
    if (chsp->Mode() == ChFiDS_ClassicChamfer)
    {
      pFunc.reset(new BRepBlend_Chamfer(S1, S2, HGuide));
      pFInv.reset(new BRepBlend_ChamfInv(S1, S2, HGuide));
    }
    else
    {
      pFunc.reset(new BRepBlend_ConstThroat(S1, S2, HGuide));
      pFInv.reset(new BRepBlend_ConstThroatInv(S1, S2, HGuide));
    }
    double dis;
    chsp->GetDist(dis);
    pFunc->Set(dis, dis, Choix);
    pFInv->Set(dis, dis, Choix);

    done = ComputeData(Data,
                       HGuide,
                       Spine,
                       lin,
                       S1,
                       I1,
                       S2,
                       I2,
                       *pFunc,
                       *pFInv,
                       PFirst,
                       MaxStep,
                       Fleche,
                       TolGuide,
                       First,
                       Last,
                       Inside,
                       Appro,
                       Forward,
                       Soldep,
                       intf,
                       intl,
                       gd1,
                       gd2,
                       gf1,
                       gf2,
                       RecOnS1,
                       RecOnS2);
    if (!done)
      return false; // ratrappage possible PMN 14/05/1998
    done = CompleteData(Data, *pFunc, lin, S1, S2, Or, gd1, gd2, gf1, gf2);
    if (!done)
      throw Standard_Failure("PerformSurf : Fail of approximation!");
  }
  else if (chsp->IsChamfer() == ChFiDS_TwoDist)
  {
    double d1, d2;
    chsp->Dists(d1, d2);

    std::unique_ptr<BlendFunc_GenChamfer>  pFunc;
    std::unique_ptr<BlendFunc_GenChamfInv> pFInv;
    if (chsp->Mode() == ChFiDS_ClassicChamfer)
    {
      pFunc.reset(new BRepBlend_Chamfer(S1, S2, HGuide));
      pFInv.reset(new BRepBlend_ChamfInv(S1, S2, HGuide));
      pFunc->Set(d1, d2, Choix);
      pFInv->Set(d1, d2, Choix);
    }
    else
    {
      occ::handle<ChFiDS_ElSpine>                    OffsetHGuide;
      NCollection_List<occ::handle<ChFiDS_ElSpine>>& ll        = Spine->ChangeElSpines();
      NCollection_List<occ::handle<ChFiDS_ElSpine>>& ll_offset = Spine->ChangeOffsetElSpines();
      NCollection_List<occ::handle<ChFiDS_ElSpine>>::Iterator ILES(ll), ILES_offset(ll_offset);
      for (; ILES.More(); ILES.Next(), ILES_offset.Next())
      {
        const occ::handle<ChFiDS_ElSpine>& aHElSpine = ILES.Value();
        if (aHElSpine == HGuide)
          OffsetHGuide = ILES_offset.Value();
      }

      if (OffsetHGuide.IsNull())
      {
        std::cout << std::endl << "Construction of offset guide failed!" << std::endl;
        // exception
      }
      pFunc.reset(new BRepBlend_ConstThroatWithPenetration(S1, S2, OffsetHGuide));
      pFInv.reset(new BRepBlend_ConstThroatWithPenetrationInv(S1, S2, OffsetHGuide));
      double Throat = std::max(d1, d2);
      pFunc->Set(Throat, Throat, Choix);
      pFInv->Set(Throat, Throat, Choix);
    }

    done = ComputeData(Data,
                       HGuide,
                       Spine,
                       lin,
                       S1,
                       I1,
                       S2,
                       I2,
                       *pFunc,
                       *pFInv,
                       PFirst,
                       MaxStep,
                       Fleche,
                       TolGuide,
                       First,
                       Last,
                       Inside,
                       Appro,
                       Forward,
                       Soldep,
                       intf,
                       intl,
                       gd1,
                       gd2,
                       gf1,
                       gf2,
                       RecOnS1,
                       RecOnS2);
    if (!done)
      return false; // ratrappage possible PMN 14/05/1998
    done = CompleteData(Data, *pFunc, lin, S1, S2, Or, gd1, gd2, gf1, gf2);
    if (!done)
      throw Standard_Failure("PerformSurf : Fail of approximation!");
  }
  else
  { // distance and angle
    double d1, angle;
    chsp->GetDistAngle(d1, angle);

    int Ch = Choix;

    BRepBlend_ChAsym    Func(S1, S2, HGuide);
    BRepBlend_ChAsymInv FInv(S1, S2, HGuide);
    Func.Set(d1, angle, Ch);
    FInv.Set(d1, angle, Ch);

    done = ComputeData(Data,
                       HGuide,
                       Spine,
                       lin,
                       S1,
                       I1,
                       S2,
                       I2,
                       Func,
                       FInv,
                       PFirst,
                       MaxStep,
                       Fleche,
                       TolGuide,
                       First,
                       Last,
                       Inside,
                       Appro,
                       Forward,
                       Soldep,
                       intf,
                       intl,
                       gd1,
                       gd2,
                       gf1,
                       gf2,
                       RecOnS1,
                       RecOnS2);

    if (!done)
      return false; // ratrappage possible PMN 14/05/1998
    done = CompleteData(Data, Func, lin, S1, S2, Or, gd1, gd2, gf1, gf2);
    if (!done)
      throw Standard_Failure("PerformSurf : Fail of approximation!");
  }

  return true;
}

void ChFi3d_ChBuilder::PerformSurf(NCollection_Sequence<occ::handle<ChFiDS_SurfData>>&,
                                   const occ::handle<ChFiDS_ElSpine>&,
                                   const occ::handle<ChFiDS_Spine>&,
                                   const int,
                                   const occ::handle<BRepAdaptor_Surface>&,
                                   const occ::handle<Adaptor3d_TopolTool>&,
                                   const occ::handle<BRepAdaptor_Curve2d>&,
                                   const occ::handle<BRepAdaptor_Surface>&,
                                   const occ::handle<BRepAdaptor_Curve2d>&,
                                   bool&,
                                   const occ::handle<BRepAdaptor_Surface>&,
                                   const occ::handle<Adaptor3d_TopolTool>&,
                                   const TopAbs_Orientation,
                                   const double,
                                   const double,
                                   const double,
                                   double&,
                                   double&,
                                   const bool,
                                   const bool,
                                   const bool,
                                   const bool,
                                   const bool,
                                   const bool,
                                   const math_Vector&)
{
  throw Standard_Failure("PerformSurf Not Implemented");
}

void ChFi3d_ChBuilder::PerformSurf(NCollection_Sequence<occ::handle<ChFiDS_SurfData>>&,
                                   const occ::handle<ChFiDS_ElSpine>&,
                                   const occ::handle<ChFiDS_Spine>&,
                                   const int,
                                   const occ::handle<BRepAdaptor_Surface>&,
                                   const occ::handle<Adaptor3d_TopolTool>&,
                                   const TopAbs_Orientation,
                                   const occ::handle<BRepAdaptor_Surface>&,
                                   const occ::handle<Adaptor3d_TopolTool>&,
                                   const occ::handle<BRepAdaptor_Curve2d>&,
                                   const occ::handle<BRepAdaptor_Surface>&,
                                   const occ::handle<BRepAdaptor_Curve2d>&,
                                   bool&,
                                   const double,
                                   const double,
                                   const double,
                                   double&,
                                   double&,
                                   const bool,
                                   const bool,
                                   const bool,
                                   const bool,
                                   const bool,
                                   const bool,
                                   const math_Vector&)
{
  throw Standard_Failure("PerformSurf Not Implemented");
}

void ChFi3d_ChBuilder::PerformSurf(NCollection_Sequence<occ::handle<ChFiDS_SurfData>>&,
                                   const occ::handle<ChFiDS_ElSpine>&,
                                   const occ::handle<ChFiDS_Spine>&,
                                   const int,
                                   const occ::handle<BRepAdaptor_Surface>&,
                                   const occ::handle<Adaptor3d_TopolTool>&,
                                   const occ::handle<BRepAdaptor_Curve2d>&,
                                   const occ::handle<BRepAdaptor_Surface>&,
                                   const occ::handle<BRepAdaptor_Curve2d>&,
                                   bool&,
                                   const TopAbs_Orientation,
                                   const occ::handle<BRepAdaptor_Surface>&,
                                   const occ::handle<Adaptor3d_TopolTool>&,
                                   const occ::handle<BRepAdaptor_Curve2d>&,
                                   const occ::handle<BRepAdaptor_Surface>&,
                                   const occ::handle<BRepAdaptor_Curve2d>&,
                                   bool&,
                                   const TopAbs_Orientation,
                                   const double,
                                   const double,
                                   const double,
                                   double&,
                                   double&,
                                   const bool,
                                   const bool,
                                   const bool,
                                   const bool,
                                   const bool,
                                   const bool,
                                   const bool,
                                   const math_Vector&)
{
  throw Standard_Failure("PerformSurf Not Implemented");
}

//=======================================================================
// function : ExtentOneCorner
// purpose  : extends the spine of the stripe S on the side of the vertex V
// PMN : 28/11/97 : Reproduces the code of fillets, and it seems to work better...
//=======================================================================

void ChFi3d_ChBuilder::ExtentOneCorner(const TopoDS_Vertex& V, const occ::handle<ChFiDS_Stripe>& S)
{
  int                       Sens  = 0;
  double                    Coeff = 0.5;
  occ::handle<ChFiDS_Spine> Spine = S->Spine();
  ChFi3d_IndexOfSurfData(V, S, Sens);
  if (Spine->IsTangencyExtremity((Sens == 1)))
    return; // No extension on queue
  double dU = Spine->LastParameter(Spine->NbEdges());
  if (Sens == 1)
  {
    Spine->SetFirstParameter(-dU * Coeff);
    Spine->SetFirstTgt(0.);
  }
  else
  {
    Spine->SetLastParameter(dU * (1. + Coeff));
    Spine->SetLastTgt(dU);
  }
}

//=======================================================================
// function : ExtentTwoCorner
// purpose  : extends the spines of the stripes contained in the list LS,
//           on the side of the vertex V
//=======================================================================

void ChFi3d_ChBuilder::ExtentTwoCorner(const TopoDS_Vertex&                                V,
                                       const NCollection_List<occ::handle<ChFiDS_Stripe>>& LS)
{
  int                                                    Sens = 0;
  NCollection_List<occ::handle<ChFiDS_Stripe>>::Iterator itel(LS);
  bool                                                   FF = true;
  bool                                                   isfirst[2];
  int                                                    Iedge[2];
  Iedge[0] = 1;
  Iedge[1] = 1;
  occ::handle<ChFiDS_Stripe> Stripe[2];
  occ::handle<ChFiDS_Spine>  Spine[2];

  int i = 0;
  for (; itel.More(); itel.Next(), i++)
  {
    ChFi3d_IndexOfSurfData(V, itel.Value(), Sens);
    if (!FF)
      if (Stripe[1] == itel.Value())
        Sens = -Sens;

    Stripe[i]  = itel.Value();
    isfirst[i] = (Sens == 1);
    Spine[i]   = Stripe[i]->Spine();
    if (!isfirst[i])
      Iedge[i] = Spine[i]->NbEdges();
    FF = false;
  }

  occ::handle<ChFiDS_ChamfSpine> chsp[2];
  double                         d[4], dis[2] = {0.0, 0.0};
  int                            j;
  TopoDS_Face                    F[4];
  double                         tmpang, tmd;

  for (i = 0, j = 0; i < 2; i++, j += 2)
  {
    chsp[i] = occ::down_cast<ChFiDS_ChamfSpine>(Spine[i]);
    ConexFaces(Spine[i], Iedge[i], F[j], F[j + 1]);

    if (chsp[i]->IsChamfer() == ChFiDS_Sym)
    {
      chsp[i]->GetDist(d[j]);
      d[j + 1] = d[j];
    }
    else if (chsp[i]->IsChamfer() == ChFiDS_TwoDist)
    {
      chsp[i]->Dists(d[j], d[j + 1]);
    }
    else
    {
      chsp[i]->GetDistAngle(tmd, tmpang);
      // an approximate calculation of distance 2 is done
      d[j]     = tmd;
      d[j + 1] = tmd * tan(tmpang);
    }
  }

  bool notfound = true;
  i             = 0;
  while (notfound && (i < 2))
  {
    j = 0;
    while (notfound && (j < 2))
    {
      if (F[i].IsSame(F[j + 2]))
      {
        dis[0] = d[i];
        //	dOnArc[0] = d[(i+1)%2];

        dis[1] = d[j + 2];
        //	dOnArc[1] = d[(j+1)%2 + 2];
        notfound = false;
      }
      j++;
    }
    i++;
  }
  // ExtentTwoCorner

  ChFiDS_State State[2];

  for (i = 0; i < 2; i++)
  {
    if (isfirst[i])
      State[i] = Spine[i]->FirstStatus();
    else
      State[i] = Spine[i]->LastStatus();
  }

  if (State[0] == ChFiDS_AllSame)
  {
    /*
       // The greatest intersection of the chamfer is found (on the incident edge)
        // with the face at end
        i = 0;
        j = 1;
        if(dOnArc[j] > dOnArc[i]) {
          int temp = i;
          i = j;
          j = temp;
        }
        ExtentOneCorner( V, Stripe[i] ); */

    // it is necessary that two chamfers touch the face at end
    for (j = 0; j < 2; j++)
      ExtentOneCorner(V, Stripe[j]);
  }
  else if ((State[0] == ChFiDS_OnSame) && (State[1] == ChFiDS_OnSame))
  {

    ExtentSpineOnCommonFace(Spine[0], Spine[1], V, dis[0], dis[1], isfirst[0], isfirst[1]);
  }
}

//=================================================================================================

void ChFi3d_ChBuilder::ExtentThreeCorner(const TopoDS_Vertex&                                V,
                                         const NCollection_List<occ::handle<ChFiDS_Stripe>>& LS)
{
  int                                          Sens = 0;
  NCollection_List<occ::handle<ChFiDS_Stripe>> check;
  bool                                         isfirst[3];
  int                                          Iedge[3];
  Iedge[0] = 1;
  Iedge[1] = 1;
  Iedge[2] = 1;
  occ::handle<ChFiDS_Spine> Spine[3];

  int i = 0;
  for (NCollection_List<occ::handle<ChFiDS_Stripe>>::Iterator itel(LS); itel.More();
       itel.Next(), i++)
  {
    occ::handle<ChFiDS_Stripe> Stripe = itel.Value();
    ChFi3d_IndexOfSurfData(V, Stripe, Sens);
    for (NCollection_List<occ::handle<ChFiDS_Stripe>>::Iterator ich(check); ich.More(); ich.Next())
    {
      if (Stripe == ich.Value())
      {
        Sens = -Sens;
        break;
      }
    }

    isfirst[i] = (Sens == 1);
    Spine[i]   = Stripe->Spine();
    if (!isfirst[i])
      Iedge[i] = Spine[i]->NbEdges();

    check.Append(Stripe);
  }

  double      d[3][2], tmd, tmpangle;
  int         j;
  TopoDS_Face F[3][2];

  occ::handle<ChFiDS_ChamfSpine> chsp[3];

  for (i = 0; i < 3; i++)
  {
    chsp[i] = occ::down_cast<ChFiDS_ChamfSpine>(Spine[i]);
    ConexFaces(Spine[i], Iedge[i], F[i][0], F[i][1]);

    if (chsp[i]->IsChamfer() == ChFiDS_Sym)
    {
      chsp[i]->GetDist(d[i][0]);
      d[i][1] = d[i][0];
    }
    else if (chsp[i]->IsChamfer() == ChFiDS_TwoDist)
    {
      chsp[i]->Dists(d[i][0], d[i][1]);
    }
    else
    {
      chsp[i]->GetDistAngle(tmd, tmpangle);
      // an approximate calculation of distance 2 is done

      d[i][0] = tmd;
      d[i][1] = tmd * tan(tmpangle);
    }
  }

  // dis[i][j] distance from chamfer i on the common face with
  // chamfer j
  double dis[3][3];

  for (i = 0; i < 3; i++)
  {
    //    for (int ii=0; ii<3; ii++) {
    //      j = (i+ii)%3;
    j             = (i + 1) % 3;
    bool notfound = true;
    int  k, l;
    k = 0;
    while (notfound && (k < 2))
    {
      l = 0;
      while (notfound && (l < 2))
      {
        if (F[i][k].IsSame(F[j][l]))
        {
          dis[i][j] = d[i][k];
          dis[j][i] = d[j][l];
          notfound  = false;
        }
        l++;
      }
      k++;
    }
    //    }
  }

  // ExtentThreeCorner
  for (i = 0; i < 3; i++)
  {
    j = (i + 1) % 3;
    ExtentSpineOnCommonFace(Spine[i], Spine[j], V, dis[i][j], dis[j][i], isfirst[i], isfirst[j]);
  }
}

//=================================================================================================

void ChFi3d_ChBuilder::SetRegul()

{
  NCollection_List<ChFiDS_Regul>::Iterator it;
  NCollection_List<TopoDS_Shape>::Iterator itc;
  NCollection_List<TopoDS_Shape>::Iterator its1;
  NCollection_List<TopoDS_Shape>::Iterator its2;
  BRepAdaptor_Surface                      S;
  BRepAdaptor_Curve2d                      PC;
  double                                   u, v, t;
  gp_Pnt                                   p;
  gp_Vec                                   n1, n2, du, dv;
  BRep_Builder                             B;
  double                                   Seuil  = M_PI / 360.;
  double                                   Seuil2 = Seuil * Seuil;
  for (it.Initialize(myRegul); it.More(); it.Next())
  {
    const ChFiDS_Regul& reg = it.Value();
    itc.Initialize(myCoup->NewEdges(reg.Curve()));
    if (itc.More())
    {
      TopoDS_Edge E = TopoDS::Edge(itc.Value());
      if (reg.IsSurface1() && reg.IsSurface2())
      {
        its1.Initialize(myCoup->NewFaces(reg.S1()));
        its2.Initialize(myCoup->NewFaces(reg.S2()));
        if (its1.More() && its2.More())
        {
          TopoDS_Face F1 = TopoDS::Face(its1.Value());
          TopoDS_Face F2 = TopoDS::Face(its2.Value());
          S.Initialize(F1, false);
          PC.Initialize(E, F1);
          t = 0.5 * (PC.FirstParameter() + PC.LastParameter());
          PC.Value(t).Coord(u, v);
          S.D1(u, v, p, du, dv);
          n1 = du.Crossed(dv);

          S.Initialize(F2, false);
          PC.Initialize(E, F2);
          PC.Value(t).Coord(u, v);
          S.D1(u, v, p, du, dv);
          n2 = du.Crossed(dv);

          if (n1.SquareMagnitude() > 1.e-14 && n2.SquareMagnitude() > 1.e-14)
          {
            n1.Normalize();
            n2.Normalize();
            double sina2 = n1.Crossed(n2).SquareMagnitude();
            if (sina2 < Seuil2)
            {
              GeomAbs_Shape cont = ChFi3d_evalconti(E, F1, F2);
              B.Continuity(E, F1, F2, cont);
            }
          }
        }
      }
    }
  }
}

//=======================================================================
// function : ConexFaces
// purpose  : F1, F2 are connected to edge so that F1 corresponds to distance
//=======================================================================

void ChFi3d_ChBuilder::ConexFaces(const occ::handle<ChFiDS_Spine>& Spine,
                                  const int                        IEdge,
                                  TopoDS_Face&                     F1,
                                  TopoDS_Face&                     F2) const
{
  BRepAdaptor_Surface Sb1, Sb2;
  TopAbs_Orientation  tmp1, tmp2;
  int                 RC, Choix;
  TopoDS_Face         f1, f2, ff1, ff2;

  // calculate the reference orientation
  //  ChFi3d_Builder::StripeOrientations is private
  SearchCommonFaces(myEFMap, Spine->Edges(1), ff1, ff2);
  ff1.Orientation(TopAbs_FORWARD);
  Sb1.Initialize(ff1);
  ff2.Orientation(TopAbs_FORWARD);
  Sb2.Initialize(ff2);
  RC = ChFi3d::ConcaveSide(Sb1, Sb2, Spine->Edges(1), tmp1, tmp2);

  // calculate the connected faces
  SearchCommonFaces(myEFMap, Spine->Edges(IEdge), f1, f2);
  Sb1.Initialize(f1);
  Sb2.Initialize(f2);
  Choix = ChFi3d::ConcaveSide(Sb1, Sb2, Spine->Edges(IEdge), tmp1, tmp2);

  if (RC % 2 != Choix % 2)
  {
    F1 = f2;
    F2 = f1;
  }
  else
  {
    F1 = f1;
    F2 = f2;
  }
}
