// Created on: 1996-10-02
// Created by: Jacques GOUSSARD
// Copyright (c) 1996-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 <LocOpe_SplitDrafts.hxx>

#include <BRep_Builder.hxx>
#include <BRep_Tool.hxx>
#include <BRepGProp.hxx>
#include <BRepTools_Substitution.hxx>
#include <Extrema_ExtPC.hxx>
#include <Geom2d_Curve.hxx>
#include <Geom2d_Line.hxx>
#include <Geom_Curve.hxx>
#include <Geom_Line.hxx>
#include <Geom_Plane.hxx>
#include <Geom_RectangularTrimmedSurface.hxx>
#include <Geom_Surface.hxx>
#include <Geom_TrimmedCurve.hxx>
#include <GeomAdaptor_Curve.hxx>
#include <GeomFill_Pipe.hxx>
#include <GeomInt_IntSS.hxx>
#include <gp.hxx>
#include <gp_Dir.hxx>
#include <gp_Pln.hxx>
#include <gp_Pnt2d.hxx>
#include <gp_Vec2d.hxx>
#include <GProp_GProps.hxx>
#include <IntAna_QuadQuadGeo.hxx>
#include <IntCurveSurface_HInter.hxx>
#include <IntCurveSurface_IntersectionPoint.hxx>
#include <IntCurveSurface_IntersectionSegment.hxx>
#include <LocOpe_BuildShape.hxx>
#include <LocOpe_Spliter.hxx>
#include <LocOpe_SplitShape.hxx>
#include <LocOpe_WiresOnShape.hxx>
#include <Precision.hxx>
#include <Standard_ConstructionError.hxx>
#include <Standard_NullObject.hxx>
#include <StdFail_NotDone.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Edge.hxx>
#include <TopoDS_Face.hxx>
#include <TopoDS_Shape.hxx>
#include <TopoDS_Vertex.hxx>
#include <TopoDS_Wire.hxx>
#include <TopTools_ShapeMapHasher.hxx>
#include <NCollection_DataMap.hxx>
#include <NCollection_List.hxx>
#include <NCollection_IndexedDataMap.hxx>

static bool NewPlane(const TopoDS_Face&,
                     const gp_Dir&,
                     const gp_Pln&,
                     const double,
                     gp_Pln&,
                     gp_Ax1&,
                     const bool);

static void MakeFace(TopoDS_Face&, NCollection_List<TopoDS_Shape>&);

static TopoDS_Edge NewEdge(const TopoDS_Edge&,
                           const TopoDS_Face&,
                           const occ::handle<Geom_Surface>&,
                           const TopoDS_Vertex&,
                           const TopoDS_Vertex&);

static bool Contains(const NCollection_List<TopoDS_Shape>&, const TopoDS_Shape&);

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

void LocOpe_SplitDrafts::Init(const TopoDS_Shape& S)
{
  myShape = S;
  myResult.Nullify();
  myMap.Clear();
}

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

void LocOpe_SplitDrafts::Perform(const TopoDS_Face& F,
                                 const TopoDS_Wire& W,
                                 const gp_Dir&      Extr,
                                 const gp_Pln&      NPl,
                                 const double       Angle)
{
  Perform(F, W, Extr, NPl, Angle, Extr, NPl, Angle, true, false);
}

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

void LocOpe_SplitDrafts::Perform(const TopoDS_Face& F,
                                 const TopoDS_Wire& W,
                                 const gp_Dir&      Extrg,
                                 const gp_Pln&      NPlg,
                                 const double       Angleg,
                                 const gp_Dir&      Extrd,
                                 const gp_Pln&      NPld,
                                 const double       Angled,
                                 const bool         ModLeft,
                                 const bool         ModRight)

{
  int j;

  myResult.Nullify();
  myMap.Clear();
  if (myShape.IsNull() || F.IsNull() || W.IsNull())
  {
    throw Standard_NullObject();
  }

  if (!ModLeft && !ModRight)
  {
    throw Standard_ConstructionError();
  }

  TopAbs_Orientation OriF = TopAbs_FORWARD;

  bool            FinS = false;
  TopExp_Explorer exp, exp2;
  for (exp.Init(myShape, TopAbs_FACE); exp.More(); exp.Next())
  {
    const TopoDS_Shape&            fac = exp.Current();
    NCollection_List<TopoDS_Shape> thelist;
    myMap.Bind(fac, thelist);
    if (fac.IsSame(F))
    {
      OriF = fac.Orientation();
      FinS = true;
    }
  }

  if (!FinS)
  {
    std::cout << "LocOpe_SplitDrafts:!Fins throw Standard_ConstructionError()" << std::endl;
    throw Standard_ConstructionError();
  }

  gp_Pln       NewPlg, NewPld;
  gp_Ax1       NormalFg, NormalFd;
  TopoDS_Shape aLocalFace = F.Oriented(OriF);

  if (!NewPlane(TopoDS::Face(aLocalFace), Extrg, NPlg, Angleg, NewPlg, NormalFg, ModLeft)
      || !NewPlane(TopoDS::Face(aLocalFace), Extrd, NPld, Angled, NewPld, NormalFd, ModRight))
  {
    //  if (!NewPlane(TopoDS::Face(F.Oriented(OriF)),
    //		Extrg,NPlg,Angleg,NewPlg,NormalFg,ModLeft) ||
    //      !NewPlane(TopoDS::Face(F.Oriented(OriF)),
    //		Extrd,NPld,Angled,NewPld,NormalFd,ModRight)) {
    return;
  }

  NCollection_List<TopoDS_Shape>::Iterator itl;
  BRep_Builder                             B;

  occ::handle<Geom_Surface> NewSg       = new Geom_Plane(NewPlg);
  occ::handle<Geom_Surface> NewSd       = new Geom_Plane(NewPld);
  occ::handle<Geom_Line>    theLinePipe = new Geom_Line(NormalFg); // ou NormalFd
  GeomInt_IntSS             i2s(NewSg, NewSd, Precision::Confusion());

  NCollection_Map<TopoDS_Shape, TopTools_ShapeMapHasher> theMap;
  occ::handle<GeomAdaptor_Curve>                         HAC = new GeomAdaptor_Curve;
  occ::handle<GeomAdaptor_Surface>                       HAS = new GeomAdaptor_Surface;
  IntCurveSurface_HInter                                 intcs;

  TopoDS_Wire theW = W;
  if (i2s.IsDone() && i2s.NbLines() > 0)
  {
    // on split le wire" << std::endl;

    GeomFill_Pipe thePipe;
    thePipe.GenerateParticularCase(true);
    thePipe.Init(theLinePipe, i2s.Line(1));
    thePipe.Perform(true);
    if (!thePipe.IsDone())
      throw Standard_ConstructionError("GeomFill_Pipe : Cannot make a surface");

    occ::handle<Geom_Surface> Spl = thePipe.Surface();
    HAS->Load(Spl);

    LocOpe_SplitShape splw(W);

    for (exp.Init(W, TopAbs_EDGE); exp.More(); exp.Next())
    {
      const TopoDS_Edge& edg = TopoDS::Edge(exp.Current());
      if (theMap.Add(edg))
      {
        TopLoc_Location         Loc;
        double                  f, l;
        occ::handle<Geom_Curve> C = BRep_Tool::Curve(edg, f, l);
        HAC->Load(C);
        intcs.Perform(HAC, HAS);
        if (!intcs.IsDone())
        {
          continue; // voir ce qu`on peut faire de mieux
        }

        if (intcs.NbSegments() >= 2)
        {
          continue; // Not yet implemented...and probably never"
        }

        if (intcs.NbSegments() == 1)
        {
          const IntCurveSurface_IntersectionPoint& P1 = intcs.Segment(1).FirstPoint();
          const IntCurveSurface_IntersectionPoint& P2 = intcs.Segment(1).SecondPoint();
          const gp_Pnt&                            pf = P1.Pnt();
          const gp_Pnt&                            pl = P2.Pnt();
          TopoDS_Vertex                            Vf, Vl;
          TopExp::Vertices(edg, Vf, Vl);
          gp_Pnt Pf   = BRep_Tool::Pnt(Vf);
          gp_Pnt Pl   = BRep_Tool::Pnt(Vl);
          double Tolf = BRep_Tool::Tolerance(Vf);
          double Toll = BRep_Tool::Tolerance(Vl);
          Tolf *= Tolf;
          Toll *= Toll;

          double dff = pf.SquareDistance(Pf);
          double dfl = pf.SquareDistance(Pl);
          double dlf = pl.SquareDistance(Pf);
          double dll = pl.SquareDistance(Pl);

          if ((dff <= Tolf && dll <= Toll) || (dlf <= Tolf && dfl <= Toll))
          {
            continue;
          }
          else
          {
            // on segmente edg en pf et pl
            TopoDS_Vertex Vnewf, Vnewl;
            B.MakeVertex(Vnewf, pf, Precision::Confusion());
            B.MakeVertex(Vnewl, pl, Precision::Confusion());
            if (P1.W() >= f && P1.W() <= l && P2.W() >= f && P2.W() <= l)
            {
              splw.Add(Vnewf, P1.W(), edg);
              splw.Add(Vnewl, P2.W(), edg);
            }
            else
            {
              continue;
            }
          }
        }
        else if (intcs.NbPoints() != 0)
        {
          TopoDS_Vertex Vf, Vl;
          TopExp::Vertices(edg, Vf, Vl);
          gp_Pnt Pf   = BRep_Tool::Pnt(Vf);
          gp_Pnt Pl   = BRep_Tool::Pnt(Vl);
          double Tolf = BRep_Tool::Tolerance(Vf);
          double Toll = BRep_Tool::Tolerance(Vl);
          Tolf *= Tolf;
          Toll *= Toll;

          for (int i = 1; i <= intcs.NbPoints(); i++)
          {
            const IntCurveSurface_IntersectionPoint& Pi  = intcs.Point(i);
            const gp_Pnt&                            pi  = Pi.Pnt();
            double                                   dif = pi.SquareDistance(Pf);
            double                                   dil = pi.SquareDistance(Pl);
            if (dif <= Tolf)
            {
            }
            else if (dil <= Toll)
            {
            }
            else
            {
              if (Pi.W() >= f && Pi.W() <= l)
              {
                TopoDS_Vertex Vnew;
                B.MakeVertex(Vnew, pi, Precision::Confusion());
                splw.Add(Vnew, Pi.W(), edg);
              }
            }
          }
        }
      }
    }

    const NCollection_List<TopoDS_Shape>& lres = splw.DescendantShapes(W);
    if (lres.Extent() != 1)
    {
      return;
    }

    if (!W.IsSame(lres.First()))
    {
      theW.Nullify();
      theW = TopoDS::Wire(lres.First());
    }

    for (exp.ReInit(); exp.More(); exp.Next())
    {
      if (!myMap.IsBound(exp.Current()))
      {
        NCollection_List<TopoDS_Shape> thelist1;
        myMap.Bind(exp.Current(), thelist1);
        for (itl.Initialize(splw.DescendantShapes(exp.Current())); itl.More(); itl.Next())
        {
          myMap(exp.Current()).Append(itl.Value());
        }
        for (exp2.Init(exp.Current(), TopAbs_VERTEX); exp2.More(); exp2.Next())
        {
          if (!myMap.IsBound(exp2.Current()))
          {
            NCollection_List<TopoDS_Shape> thelist2;
            myMap.Bind(exp2.Current(), thelist2);
            myMap(exp2.Current()).Append(exp2.Current());
          }
        }
      }
    }
  }
  else
  {
    for (exp.Init(W, TopAbs_EDGE); exp.More(); exp.Next())
    {
      if (!myMap.IsBound(exp.Current()))
      {
        NCollection_List<TopoDS_Shape> thelist3;
        myMap.Bind(exp.Current(), thelist3);
        myMap(exp.Current()).Append(exp.Current());
        for (exp2.Init(exp.Current(), TopAbs_VERTEX); exp2.More(); exp2.Next())
        {
          if (!myMap.IsBound(exp2.Current()))
          {
            NCollection_List<TopoDS_Shape> thelist4;
            myMap.Bind(exp2.Current(), thelist4);
            myMap(exp2.Current()).Append(exp2.Current());
          }
        }
      }
    }
  }

  // On split la face par le wire

  occ::handle<LocOpe_WiresOnShape> WonS = new LocOpe_WiresOnShape(myShape);
  LocOpe_Spliter                   Spls(myShape);
  WonS->Bind(theW, F);

  // JAG Le code suivant marchera apres integration de thick0
  //  LocOpe_FindEdges fined(W,F);
  //  for (fined.InitIterator(); fined.More(); fined.Next()) {
  //    WonS->Bind(fined.EdgeFrom(),fined.EdgeTo());
  //  }

  Spls.Perform(WonS);
  if (!Spls.IsDone())
  {
    return;
  }

  TopoDS_Shape                          Res     = Spls.ResultingShape();
  const NCollection_List<TopoDS_Shape>& theLeft = Spls.DirectLeft();

  // Descendants
  for (exp.Init(myShape, TopAbs_FACE); exp.More(); exp.Next())
  {
    const TopoDS_Shape& fac = exp.Current();
    for (itl.Initialize(Spls.DescendantShapes(fac)); itl.More(); itl.Next())
    {
      myMap(fac).Append(itl.Value());
    }
  }

  NCollection_DataMap<TopoDS_Shape, TopoDS_Shape, TopTools_ShapeMapHasher> MapW;
  for (exp.Init(theW, TopAbs_EDGE); exp.More(); exp.Next())
  {
    if (!MapW.IsBound(exp.Current()))
    {
      MapW.Bind(exp.Current(), TopoDS_Shape());
      for (exp2.Init(exp.Current(), TopAbs_VERTEX); exp2.More(); exp2.Next())
      {
        if (!MapW.IsBound(exp2.Current()))
        {
          MapW.Bind(exp2.Current(), TopoDS_Shape());
        }
      }
    }
  }

  NCollection_IndexedDataMap<TopoDS_Shape, NCollection_List<TopoDS_Shape>, TopTools_ShapeMapHasher>
    theMapEF;
  TopExp::MapShapesAndAncestors(Res, TopAbs_EDGE, TopAbs_FACE, theMapEF);

  // On stocke les geometries potentiellement generees par les edges
  NCollection_IndexedDataMap<TopoDS_Shape, TopoDS_Shape, TopTools_ShapeMapHasher> MapEV; // genere
  NCollection_DataMap<TopoDS_Shape, NCollection_List<TopoDS_Shape>, TopTools_ShapeMapHasher> MapSg,
    MapSd; // image a gauche et a droite

  int Nbedges, index;
  for (itl.Initialize(myMap(F)); itl.More(); itl.Next())
  {
    const TopoDS_Shape& fac = TopoDS::Face(itl.Value());
    for (exp.Init(fac, TopAbs_EDGE); exp.More(); exp.Next())
    {
      const TopoDS_Shape& edg = exp.Current();
      if (MapEV.FindIndex(edg) != 0)
      {
        continue;
      }
      if (MapW.IsBound(edg))
      { // edge du wire initial
        TopLoc_Location         Loc;
        double                  f, l;
        occ::handle<Geom_Curve> C = BRep_Tool::Curve(TopoDS::Edge(edg), Loc, f, l);
        if (C->DynamicType() == STANDARD_TYPE(Geom_TrimmedCurve))
        {
          C = occ::down_cast<Geom_TrimmedCurve>(C)->BasisCurve();
        }
        C = occ::down_cast<Geom_Curve>(C->Transformed(Loc.Transformation()));

        GeomFill_Pipe thePipe;
        thePipe.GenerateParticularCase(true);
        thePipe.Init(theLinePipe, C);
        thePipe.Perform(true);
        if (!thePipe.IsDone())
          throw Standard_ConstructionError("GeomFill_Pipe : Cannot make a surface");

        occ::handle<Geom_Surface> thePS = thePipe.Surface();
        if (thePS->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface))
        {
          thePS = occ::down_cast<Geom_RectangularTrimmedSurface>(thePS)->BasisSurface();
        }

        TopoDS_Face NewFace;
        B.MakeFace(NewFace, thePS, Precision::Confusion());
        MapEV.Add(edg, NewFace);
      }
      else
      { // on recupere la face.
        index = theMapEF.FindIndex(edg);
        if (theMapEF(index).Extent() != 2)
        {
          return; // NotDone
        }
        TopoDS_Face theFace;
        if (theMapEF(index).First().IsSame(fac))
        {
          MapEV.Add(edg, theMapEF(index).Last());
        }
        else
        {
          MapEV.Add(edg, theMapEF(index).First());
        }
      }
    }
  }

  NCollection_DataMap<TopoDS_Shape, TopoDS_Shape, TopTools_ShapeMapHasher> MapSonS;

  Nbedges = MapEV.Extent();
  for (index = 1; index <= Nbedges; index++)
  {
    for (exp.Init(MapEV.FindKey(index), TopAbs_VERTEX); exp.More(); exp.Next())
    {
      const TopoDS_Vertex& vtx = TopoDS::Vertex(exp.Current());
      if (MapEV.FindIndex(vtx) != 0)
      {
        continue;
      }

      // Localisation du vertex :
      //    - entre 2 edges d`origine : on recupere l`edge qui n`est
      //                                pas dans F
      //    - entre 2 edges du wire   : droite
      //    - mixte                   : intersection de surfaces
      for (j = 1; j <= Nbedges; j++)
      {
        if (j == index)
        {
          continue;
        }
        for (exp2.Init(MapEV.FindKey(j), TopAbs_VERTEX); exp2.More(); exp2.Next())
        {
          const TopoDS_Shape& vtx2 = exp2.Current();
          if (vtx2.IsSame(vtx))
          {
            break;
          }
        }
        if (exp2.More())
        {
          break;
        }
      }
      int                 Choice = 0;
      const TopoDS_Shape& edg1   = MapEV.FindKey(index);
      TopoDS_Shape        edg2;
      if (j <= Nbedges)
      {
        edg2 = MapEV.FindKey(j);
      }
      else
      {
        edg2 = edg1;
      }
      if (MapW.IsBound(edg1))
      {
        if (j > Nbedges)
        {             // doit correspondre a edge ferme
          Choice = 2; // droite
        }
        else if (MapW.IsBound(MapEV.FindKey(j)))
        {
          Choice = 2; // droite
        }
        else
        {
          Choice = 3; // mixte
        }
      }
      else
      {
        if (j > Nbedges)
        {             // doit correspondre a edge ferme
          Choice = 1; // edge a retrouver
        }
        else if (!MapW.IsBound(MapEV.FindKey(j)))
        {
          Choice = 1; // edge a retrouver
        }
        else
        {
          Choice = 3; // mixte
        }
      }
      occ::handle<Geom_Curve>   Newc;
      occ::handle<Geom2d_Curve> newCs1, newCs2;
      double                    knownp = 0;
      TopoDS_Edge               Ebind;
      switch (Choice)
      {
        case 1: {
          for (exp2.Init(Res, TopAbs_EDGE); exp2.More(); exp2.Next())
          {
            if (exp2.Current().IsSame(edg1) || exp2.Current().IsSame(edg2))
            {
              continue;
            }
            //	    for (TopExp_Explorer exp3(exp2.Current().Oriented(TopAbs_FORWARD),
            TopExp_Explorer exp3(exp2.Current().Oriented(TopAbs_FORWARD), TopAbs_VERTEX);
            for (; exp3.More(); exp3.Next())
            {
              if (exp3.Current().IsSame(vtx))
              {
                break;
              }
            }
            if (exp3.More())
            {
              break;
            }
          }
          if (exp2.More())
          {
            double          f, l;
            TopLoc_Location Loc;
            Newc   = BRep_Tool::Curve(TopoDS::Edge(exp2.Current()), Loc, f, l);
            Newc   = occ::down_cast<Geom_Curve>(Newc->Transformed(Loc.Transformation()));
            Ebind  = TopoDS::Edge(exp2.Current());
            knownp = BRep_Tool::Parameter(vtx, Ebind);
          }
          else
          { // droite ??? il vaudrait mieux sortir
            return;

            //	    gp_Lin theLine(NormalFg);
            //	    theLine.Translate(NormalF.Location(),BRep_Tool::Pnt(vtx));
            //	    Newc = new Geom_Line(theLine);
            //	    knownp = 0.;
          }
        }
        break;
        case 2: {
          gp_Lin theLine(NormalFg);
          theLine.Translate(NormalFg.Location(), BRep_Tool::Pnt(vtx));
          Newc   = new Geom_Line(theLine);
          knownp = 0.;
        }
        break;
        case 3: {
          const TopoDS_Face&        F1    = TopoDS::Face(MapEV.FindFromKey(edg1));
          const TopoDS_Face&        F2    = TopoDS::Face(MapEV.FindFromKey(edg2));
          occ::handle<Geom_Surface> S1    = BRep_Tool::Surface(F1);
          occ::handle<Geom_Surface> S2    = BRep_Tool::Surface(F2);
          bool                      AppS1 = false;
          bool                      AppS2 = false;
          if (S1->DynamicType() != STANDARD_TYPE(Geom_Plane))
          {
            AppS1 = true;
          }
          if (S2->DynamicType() != STANDARD_TYPE(Geom_Plane))
          {
            AppS2 = true;
          }
          i2s.Perform(S1, S2, Precision::Confusion(), true, AppS1, AppS2);
          if (!i2s.IsDone() || i2s.NbLines() <= 0)
          {
            return;
          }

          double            pmin = 0, Dist2, Dist2Min, Glob2Min = RealLast();
          GeomAdaptor_Curve TheCurve;

          int    i, imin, k;
          gp_Pnt pv = BRep_Tool::Pnt(vtx);
          imin      = 0;
          for (i = 1; i <= i2s.NbLines(); i++)
          {
            TheCurve.Load(i2s.Line(i));
            Extrema_ExtPC myExtPC(pv, TheCurve);

            if (myExtPC.IsDone())
            {
              gp_Pnt p1b, p2b;
              double thepmin = TheCurve.FirstParameter();
              myExtPC.TrimmedSquareDistances(Dist2Min, Dist2, p1b, p2b);
              if (Dist2 < Dist2Min)
              {
                thepmin = TheCurve.LastParameter();
              }
              for (k = 1; k <= myExtPC.NbExt(); k++)
              {
                Dist2 = myExtPC.SquareDistance(k);
                if (Dist2 < Dist2Min)
                {
                  Dist2Min = Dist2;
                  thepmin  = myExtPC.Point(k).Parameter();
                }
              }

              if (Dist2Min < Glob2Min)
              {
                Glob2Min = Dist2Min;
                pmin     = thepmin;
                imin     = i;
              }
            }
          }
          if (imin == 0)
          {
            return;
          }

          Newc   = i2s.Line(imin);
          knownp = pmin;
          if (AppS1)
          {
            newCs1 = i2s.LineOnS1(imin);
          }
          if (AppS2)
          {
            newCs2 = i2s.LineOnS2(imin);
          }
        }
        break;
      }

      // Determination des vertex par intersection sur Plg ou/et Pld

      HAC->Load(Newc);
      int           nbfois = 2;
      TopoDS_Vertex vtx1, vtx2;
      double        p1 = 0, p2 = 0;
      bool          IsLeft = false;
      if (Choice == 1)
      {
        // edge retrouve : on ne fait qu`une seule intersection
        // il faut utiliser Plg ou Pld

        int indedgf = theMapEF.FindIndex(edg1);
        for (itl.Initialize(theMapEF(indedgf)); itl.More(); itl.Next())
        {
          if (Contains(myMap(F), itl.Value()))
          {
            if (Contains(theLeft, itl.Value()))
            {
              HAS->Load(NewSg);
              IsLeft = true;
            }
            else
            {
              HAS->Load(NewSd);
              IsLeft = false;
            }

            nbfois = 1;
            vtx2   = vtx;
            p2     = knownp;
            break;
          }
        }
        if (!itl.More())
        {
          std::cout << "LocOpe_SplitDrafts: betite probleme " << std::endl;
          return;
        }
      }
      else
      {
        HAS->Load(NewSg);
      }

      for (int it = 1; it <= nbfois; it++)
      {
        if (it == 2)
        {
          HAS->Load(NewSd);
        }

        intcs.Perform(HAC, HAS);
        if (!intcs.IsDone())
        {
          return; // voir ce qu`on peut faire de mieux
        }
        int    imin  = 1;
        double delta = std::abs(knownp - intcs.Point(1).W());
        for (int i = 2; i <= intcs.NbPoints(); i++)
        {
          double newdelta = std::abs(knownp - intcs.Point(i).W());
          if (newdelta < delta)
          {
            imin  = i;
            delta = newdelta;
          }
        }
        if (it == 1)
        {
          B.MakeVertex(vtx1, intcs.Point(imin).Pnt(), Precision::Confusion());
          p1     = intcs.Point(imin).W();
          knownp = p1;
        }
        else
        {
          B.MakeVertex(vtx2, intcs.Point(imin).Pnt(), Precision::Confusion());
          p2 = intcs.Point(imin).W();
        }
      }
      if (std::abs(p1 - p2) > Precision::PConfusion())
      {
        TopoDS_Edge NewEdge;
        B.MakeEdge(NewEdge, Newc, Precision::Confusion());
        if (p1 < p2)
        {
          B.Add(NewEdge, vtx1.Oriented(TopAbs_FORWARD));
          B.Add(NewEdge, vtx2.Oriented(TopAbs_REVERSED));
        }
        else
        {
          B.Add(NewEdge, vtx1.Oriented(TopAbs_REVERSED));
          B.Add(NewEdge, vtx2.Oriented(TopAbs_FORWARD));
        }
        B.UpdateVertex(vtx1, p1, NewEdge, Precision::Confusion());
        B.UpdateVertex(vtx2, p2, NewEdge, Precision::Confusion());
        if (!newCs1.IsNull())
        {
          B.UpdateEdge(NewEdge,
                       newCs1,
                       TopoDS::Face(MapEV.FindFromKey(edg1)),
                       Precision::Confusion());
        }

        if (!newCs2.IsNull())
        {
          B.UpdateEdge(NewEdge,
                       newCs2,
                       TopoDS::Face(MapEV.FindFromKey(edg2)),
                       Precision::Confusion());
        }

        MapEV.Add(vtx, NewEdge);

        if (Choice == 1)
        {
          TopoDS_Shape aLocalEdge = Ebind.EmptyCopied();
          TopoDS_Edge  NE         = TopoDS::Edge(aLocalEdge);
          //	  TopoDS_Edge NE = TopoDS::Edge(Ebind.EmptyCopied());
          for (exp2.Init(Ebind, TopAbs_VERTEX); exp2.More(); exp2.Next())
          {
            const TopoDS_Vertex& thevtx = TopoDS::Vertex(exp2.Current());
            if (thevtx.IsSame(vtx))
            {
              B.Add(NE, vtx1.Oriented(thevtx.Orientation()));
              B.UpdateVertex(vtx1, p1, NE, Precision::Confusion());
            }
            else
            {
              B.Add(NE, thevtx);
              double theprm = BRep_Tool::Parameter(thevtx, Ebind);
              B.UpdateVertex(thevtx, theprm, NE, BRep_Tool::Tolerance(thevtx));
            }
          }
          MapSonS.Bind(Ebind, NE.Oriented(TopAbs_FORWARD));
          if (IsLeft)
          {
            NCollection_List<TopoDS_Shape> thelist5;
            MapSg.Bind(vtx, thelist5);
            MapSg(vtx).Append(vtx1);
          }
          else
          {
            NCollection_List<TopoDS_Shape> thelist6;
            MapSd.Bind(vtx, thelist6);
            MapSd(vtx).Append(vtx1);
          }
        }
        else
        {
          NCollection_List<TopoDS_Shape> thelist7, thelist8;
          MapSg.Bind(vtx, thelist7);
          MapSd.Bind(vtx, thelist8);
          MapSg(vtx).Append(vtx1);
          MapSd(vtx).Append(vtx2);
        }
      }
      else
      {
        MapEV.Add(vtx, vtx2); // on peut avoir vtx2 = vtx si choix == 1
        if (Choice == 1)
        {
          if (IsLeft)
          {
            NCollection_List<TopoDS_Shape> thelist9;
            MapSg.Bind(vtx, thelist9);
            MapSg(vtx).Append(vtx);
          }
          else
          {
            NCollection_List<TopoDS_Shape> thelist10;
            MapSd.Bind(vtx, thelist10);
            MapSd(vtx).Append(vtx);
          }
        }
        else
        {
          NCollection_List<TopoDS_Shape> thelist11, thelist12;
          MapSg.Bind(vtx, thelist11);
          MapSd.Bind(vtx, thelist12);
          MapSg(vtx).Append(vtx2);
          MapSd(vtx).Append(vtx2);
        }
      }
    }
  }

  theMap.Clear();
  for (exp.Init(theW, TopAbs_EDGE); exp.More(); exp.Next())
  {
    const TopoDS_Edge& edg = TopoDS::Edge(exp.Current());
    if (!theMap.Add(edg))
    { // precaution sans doute inutile...
      continue;
    }
    int                            indedg = MapEV.FindIndex(edg);
    TopoDS_Face&                   GenF   = TopoDS::Face(MapEV(indedg));
    NCollection_List<TopoDS_Shape> thelist13, thelist14;
    MapSg.Bind(edg, thelist13); // genere a gauche
    MapSd.Bind(edg, thelist14); // genere a droite
    TopoDS_Vertex Vf, Vl;
    TopoDS_Shape  aLocalEdge = edg.Oriented(TopAbs_FORWARD);
    TopExp::Vertices(TopoDS::Edge(aLocalEdge), Vf, Vl);
    //    TopExp::Vertices(TopoDS::Edge(edg.Oriented(TopAbs_FORWARD)),Vf,Vl);
    TopoDS_Shape Gvf = MapEV.FindFromKey(Vf);
    TopoDS_Shape Gvl = MapEV.FindFromKey(Vl);

    /* Le code suivant est OK. On essaie de l`ameliorer

        if (Gvf.ShapeType() == TopAbs_VERTEX &&
        Gvl.ShapeType() == TopAbs_VERTEX) {
          // en fait on doit pouvoir avoir 1 face a 2 cotes...
          if (Gvf.IsSame(Vf)) {
        MapW(edg) = edg;
        MapSg(edg).Append(edg.Oriented(TopAbs_FORWARD));
        MapSd(edg).Append(edg.Oriented(TopAbs_FORWARD));
          }
          else {
        TopoDS_Edge NewEdg = NewEdge(edg,
                         GenF,NewSg,
                         TopoDS::Vertex(Gvf),
                         TopoDS::Vertex(Gvl));
        if (NewEdg.IsNull()) {
          return;
        }
        MapW(edg) = NewEdg;
        MapSg(edg).Append(NewEdg);
        MapSd(edg).Append(NewEdg);
          }
        }
        else if (Gvf.ShapeType() == TopAbs_VERTEX  ||
             Gvl.ShapeType() == TopAbs_VERTEX) {      // face triangulaire
          TopoDS_Vertex Vfd,Vld,Vfg,Vlg;
          if (Gvf.ShapeType() == TopAbs_VERTEX) {
        Vfg = TopoDS::Vertex(Gvf);
        Vfd = Vfg;
        Vlg = TopoDS::Vertex(MapSg(Vl).First());
        Vld = TopoDS::Vertex(MapSd(Vl).First());
          }
          else {
        Vlg = TopoDS::Vertex(Gvl);
        Vld = Vlg;
        Vfg = TopoDS::Vertex(MapSg(Vf).First());
        Vfd = TopoDS::Vertex(MapSd(Vf).First());
          }

          TopoDS_Edge NewEdgg = NewEdge(edg,GenF,NewSg,Vfg,Vlg);
          if (NewEdgg.IsNull()) {
        return;
          }

          TopoDS_Edge NewEdgd = NewEdge(edg,GenF,NewSd,Vfd,Vld);
          if (NewEdgg.IsNull()) {
        return;
          }
          MapSg(edg).Append(NewEdgg);
          MapSd(edg).Append(NewEdgd);

          NCollection_List<TopoDS_Shape> theedges;
          theedges.Append(NewEdgg);
          theedges.Append(NewEdgd);
          if (Gvf.ShapeType() == TopAbs_EDGE) {
        theedges.Append(Gvf);
          }
          else {//if (Gvl.ShapeType() == TopAbs_EDGE) {
        theedges.Append(Gvl);
          }
          MakeFace(GenF,theedges);
          MapW(edg) = GenF;
        }
        else {
          // une face a 4 cotes
          TopoDS_Vertex Vfd,Vld,Vfg,Vlg;

          Vfg = TopoDS::Vertex(MapSg(Vf).First());
          Vfd = TopoDS::Vertex(MapSd(Vf).First());
          Vlg = TopoDS::Vertex(MapSg(Vl).First());
          Vld = TopoDS::Vertex(MapSd(Vl).First());

          TopoDS_Vertex VVf1,VVl1,VVf2,VVl2;
          TopExp::Vertices(TopoDS::Edge(Gvf.Oriented(TopAbs_FORWARD)),VVf1,VVl1);
          TopExp::Vertices(TopoDS::Edge(Gvl.Oriented(TopAbs_FORWARD)),VVf2,VVl2);

          TopoDS_Edge NewEdgg = NewEdge(edg,GenF,NewSg,Vfg,Vlg);
          if (NewEdgg.IsNull()) {
        return;
          }

          TopoDS_Edge NewEdgd = NewEdge(edg,GenF,NewSd,Vfd,Vld);
          if (NewEdgd.IsNull()) {
        return;
          }

          if ((VVf1.IsSame(Vfg) && VVf2.IsSame(Vlg)) ||
          (VVf1.IsSame(Vfd) && VVf2.IsSame(Vld))) {
        // 4 cotes
        MapSg(edg).Append(NewEdgg);
        MapSd(edg).Append(NewEdgd);

        NCollection_List<TopoDS_Shape> theedges;
        theedges.Append(NewEdgg);
        theedges.Append(NewEdgd);
        theedges.Append(Gvf);
        theedges.Append(Gvl);

        MakeFace(GenF,theedges);
        MapW(edg) = GenF;
          }
          else {
    #ifdef OCCT_DEBUG
        std::cout << "Pb d'analyse" << std::endl;
    #endif
        return;
          }
        }
    */
    // nouveau code

    TopoDS_Vertex Vfd, Vld, Vfg, Vlg;
    if (Gvf.ShapeType() == TopAbs_VERTEX)
    {
      Vfg = TopoDS::Vertex(Gvf);
      Vfd = Vfg;
    }
    else
    {
      Vfg = TopoDS::Vertex(MapSg(Vf).First());
      Vfd = TopoDS::Vertex(MapSd(Vf).First());
    }
    if (Gvl.ShapeType() == TopAbs_VERTEX)
    {
      Vlg = TopoDS::Vertex(Gvl);
      Vld = Vlg;
    }
    else
    {
      Vlg = TopoDS::Vertex(MapSg(Vl).First());
      Vld = TopoDS::Vertex(MapSd(Vl).First());
    }

    TopoDS_Edge NewEdgg = NewEdge(edg, GenF, NewSg, Vfg, Vlg);
    if (NewEdgg.IsNull())
    {
      return;
    }

    TopoDS_Edge NewEdgd = NewEdge(edg, GenF, NewSd, Vfd, Vld);
    if (NewEdgg.IsNull())
    {
      return;
    }

    bool isedg = false;
    if (Gvf.ShapeType() == TopAbs_VERTEX && Gvl.ShapeType() == TopAbs_VERTEX)
    {
      // edg ou face a 2 cotes

      // Comparaison NewEdgg et NewEdgd
      double                  fg, lg, fd, ld;
      occ::handle<Geom_Curve> Cg   = BRep_Tool::Curve(NewEdgg, fg, lg);
      occ::handle<Geom_Curve> Cd   = BRep_Tool::Curve(NewEdgd, fd, ld);
      double                  prmg = (fg + lg) / 2.;
      double                  prmd = (fd + ld) / 2.;
      gp_Pnt                  pg   = Cg->Value(prmg);
      gp_Pnt                  pd   = Cd->Value(prmd);
      double Tol = std::max(BRep_Tool::Tolerance(NewEdgg), BRep_Tool::Tolerance(NewEdgg));
      if (pg.SquareDistance(pd) <= Tol * Tol)
      {
        isedg = true;
        // raffinement pour essayer de partager l`edge de depart...
        bool modified = true;
        if (Gvf.IsSame(Vf) && Gvl.IsSame(Vl))
        {
          // Comparaison avec l`edge de depart
          Cd   = BRep_Tool::Curve(edg, fd, ld);
          prmd = (fd + ld) / 2.;
          pd   = Cd->Value(prmd);
          Tol  = std::max(BRep_Tool::Tolerance(NewEdgg), BRep_Tool::Tolerance(edg));
          if (pg.SquareDistance(pd) <= Tol * Tol)
          {
            modified = false;
          }
        }

        if (!modified)
        {
          MapW(edg) = edg;
          MapSg(edg).Append(edg);
          MapSd(edg).Append(edg);
        }
        else
        {
          MapW(edg) = NewEdgg;
          MapSg(edg).Append(NewEdgg);
          MapSd(edg).Append(NewEdgg);
        }
      }
    }

    if (!isedg)
    {
      // face a 2 ou 3 ou 4 cotes
      MapSg(edg).Append(NewEdgg);
      MapSd(edg).Append(NewEdgd);

      NCollection_List<TopoDS_Shape> theedges;
      theedges.Append(NewEdgg);
      theedges.Append(NewEdgd);
      if (Gvf.ShapeType() == TopAbs_EDGE)
      {
        theedges.Append(Gvf);
      }
      if (Gvl.ShapeType() == TopAbs_EDGE)
      {
        theedges.Append(Gvl);
      }
      MakeFace(GenF, theedges);
      MapW(edg) = GenF;
    }
  }

  NCollection_Map<TopoDS_Shape, TopTools_ShapeMapHasher> mapedgadded;
  NCollection_List<TopoDS_Shape>                         thefaces;

  for (itl.Initialize(myMap(F)); itl.More(); itl.Next())
  {
    const TopoDS_Face& fac = TopoDS::Face(itl.Value());
    theMap.Clear();
    TopoDS_Face DrftFace; // elle est FORWARD
    bool        IsLeft;
    if (Contains(theLeft, fac))
    {
      B.MakeFace(DrftFace, NewSg, BRep_Tool::Tolerance(fac));
      IsLeft = true;
    }
    else
    {
      B.MakeFace(DrftFace, NewSd, BRep_Tool::Tolerance(fac));
      IsLeft = false;
    }

    TopExp_Explorer exp3;
    for (exp3.Init(fac.Oriented(TopAbs_FORWARD), TopAbs_WIRE); exp3.More(); exp3.Next())
    {
      const TopoDS_Shape& wir = exp3.Current();
      TopoDS_Wire         NewWireOnF;
      B.MakeWire(NewWireOnF);
      for (exp.Init(wir.Oriented(TopAbs_FORWARD), TopAbs_EDGE); exp.More(); exp.Next())
      {
        const TopoDS_Edge& edg = TopoDS::Edge(exp.Current());
        if (!theMap.Add(edg))
        { // precaution sans doute inutile...
          continue;
        }
        if (MapW.IsBound(edg))
        { // edge du wire d`origine
          NCollection_List<TopoDS_Shape>::Iterator itld;
          TopAbs_Orientation                       ored = edg.Orientation();
          if (IsLeft)
          {
            itld.Initialize(MapSg(edg));
          }
          else
          {
            itld.Initialize(MapSd(edg));
          }
          for (; itld.More(); itld.Next())
          {
            if (itld.Value().Orientation() == TopAbs_REVERSED)
            {
              ored = TopAbs::Reverse(ored);
            }
            TopoDS_Shape aLocalEdge = itld.Value().Oriented(ored);
            B.Add(NewWireOnF, TopoDS::Edge(aLocalEdge));
            //	    B.Add(NewWireOnF,TopoDS::Edge(itld.Value().Oriented(ored)));
          }
        }
        else
        {
          occ::handle<Geom_Surface> NewS;
          if (IsLeft)
          {
            NewS = NewSg;
          }
          else
          {
            NewS = NewSd;
          }
          int                indedg = MapEV.FindIndex(edg);
          const TopoDS_Face& GenF   = TopoDS::Face(MapEV(indedg));
          TopoDS_Vertex      Vf, Vl;
          TopoDS_Shape       aLocalEdge = edg.Oriented(TopAbs_FORWARD);
          TopExp::Vertices(TopoDS::Edge(aLocalEdge), Vf, Vl);
          //	  TopExp::Vertices(TopoDS::Edge(edg.Oriented(TopAbs_FORWARD)),Vf,Vl);
          TopoDS_Shape Gvf = MapEV.FindFromKey(Vf);
          TopoDS_Shape Gvl = MapEV.FindFromKey(Vl);
          if (Gvf.ShapeType() == TopAbs_VERTEX && Gvl.ShapeType() == TopAbs_VERTEX)
          {
            if (!Gvf.IsSame(Vf) || !Gvl.IsSame(Vl))
            {
              TopoDS_Edge NewEdg =
                NewEdge(edg, GenF, NewS, TopoDS::Vertex(Gvf), TopoDS::Vertex(Gvl));
              if (NewEdg.IsNull())
              {
                return;
              }

              MapSonS.Bind(edg, NewEdg);

              if (NewEdg.Orientation() == TopAbs_REVERSED)
              {
                NewEdg.Orientation(TopAbs::Reverse(edg.Orientation()));
              }
              else
              {
                NewEdg.Orientation(edg.Orientation());
              }
              B.Add(NewWireOnF, NewEdg);
            }
            else
            { // Frozen???
              B.Add(NewWireOnF, edg);
            }
          }
          else
          {
            TopoDS_Vertex Vff, Vll;
            if (Gvf.ShapeType() == TopAbs_VERTEX)
            {
              Vff = TopoDS::Vertex(Gvf);
            }
            else
            {
              if (IsLeft)
              {
                Vff = TopoDS::Vertex(MapSg(Vf).First());
              }
              else
              {
                Vff = TopoDS::Vertex(MapSd(Vf).First());
              }
            }
            if (Gvl.ShapeType() == TopAbs_VERTEX)
            {
              Vll = TopoDS::Vertex(Gvl);
            }
            else
            {
              if (IsLeft)
              {
                Vll = TopoDS::Vertex(MapSg(Vl).First());
              }
              else
              {
                Vll = TopoDS::Vertex(MapSd(Vl).First());
              }
            }

            TopoDS_Edge NewEdg = NewEdge(edg, GenF, NewS, Vff, Vll);
            if (NewEdg.IsNull())
            {
              return;
            }

            if (!MapW.IsBound(Vf) && !MapW.IsBound(Vl))
            {
              MapSonS.Bind(edg, NewEdg);
            }
            //	    else if (MapW.IsBound(Vf) && MapW.IsBound(Vl)) {

            //	    }
            else
            {
              if (MapW.IsBound(Vf))
              {
                if (Gvf.ShapeType() != TopAbs_EDGE || mapedgadded.Contains(Gvf))
                {
                  MapSonS.Bind(edg, NewEdg);
                }
                else
                {
                  TopoDS_Wire NewWir;
                  B.MakeWire(NewWir);
                  B.Add(NewWir, NewEdg);

                  TopoDS_Vertex Vf2, Vl2;
                  TopExp::Vertices(TopoDS::Edge(Gvf), Vf2, Vl2);

                  // TopAbs_Orientation ornw = NewEdg.Orientation();

                  // ici bug orientation : voir tspdrft6

                  //		  if ((ornw == TopAbs_FORWARD && Vl2.IsSame(Vff)) ||
                  //		      (ornw == TopAbs_REVERSED && Vl2.IsSame(Vll))) {
                  if (Vl2.IsSame(Vff))
                  {
                    B.Add(NewWir, Gvf.Oriented(TopAbs_FORWARD));
                  }
                  else
                  {
                    B.Add(NewWir, Gvf.Oriented(TopAbs_REVERSED));
                  }
                  mapedgadded.Add(Gvf);
                  MapSonS.Bind(edg, NewWir); // NewWire est FORWARD
                }
              }
              else
              {
                if (Gvl.ShapeType() != TopAbs_EDGE || mapedgadded.Contains(Gvl))
                {
                  MapSonS.Bind(edg, NewEdg);
                }
                else
                {
                  TopoDS_Wire NewWir;
                  B.MakeWire(NewWir);
                  B.Add(NewWir, NewEdg);

                  TopoDS_Vertex Vf2, Vl2;
                  TopExp::Vertices(TopoDS::Edge(Gvl), Vf2, Vl2);

                  // TopAbs_Orientation ornw = NewEdg.Orientation();

                  // ici bug orientation : voir tspdrft6

                  //		  if ((ornw == TopAbs_FORWARD && Vl2.IsSame(Vff)) ||
                  //		      (ornw == TopAbs_REVERSED && Vl2.IsSame(Vll))) {
                  if (Vf2.IsSame(Vll))
                  {
                    B.Add(NewWir, Gvl.Oriented(TopAbs_FORWARD));
                  }
                  else
                  {
                    B.Add(NewWir, Gvl.Oriented(TopAbs_REVERSED));
                  }
                  mapedgadded.Add(Gvl);
                  MapSonS.Bind(edg, NewWir); // NewWire est FORWARD
                }
              }
            }
            if (NewEdg.Orientation() == TopAbs_REVERSED)
            {
              NewEdg.Orientation(TopAbs::Reverse(edg.Orientation()));
            }
            else
            {
              NewEdg.Orientation(edg.Orientation());
            }
            B.Add(NewWireOnF, NewEdg);
          }
        }
      }
      B.Add(DrftFace, NewWireOnF.Oriented(wir.Orientation()));
    }
    thefaces.Append(DrftFace);
  }

  BRepTools_Substitution                                                             theSubs;
  NCollection_DataMap<TopoDS_Shape, TopoDS_Shape, TopTools_ShapeMapHasher>::Iterator itdmss;
  for (itdmss.Initialize(MapSonS); itdmss.More(); itdmss.Next())
  {
    NCollection_List<TopoDS_Shape> lsubs;
    for (exp.Init(itdmss.Value(), TopAbs_EDGE); exp.More(); exp.Next())
    {
      lsubs.Append(exp.Current());
    }
    theSubs.Substitute(itdmss.Key(), lsubs);
  }

  // on reconstruit les faces
  for (exp.Init(Res, TopAbs_FACE); exp.More(); exp.Next())
  {
    if (Contains(myMap(F), exp.Current()))
    {
      continue;
    }
    theSubs.Build(exp.Current());
  }

  // Stockage des descendants
  //  for (NCollection_DataMap<TopoDS_Shape, NCollection_List<TopoDS_Shape>,
  //  TopTools_ShapeMapHasher>::Iterator itdmsls(myMap);
  NCollection_DataMap<TopoDS_Shape, NCollection_List<TopoDS_Shape>, TopTools_ShapeMapHasher>::
    Iterator itdmsls(myMap);
  for (; itdmsls.More(); itdmsls.Next())
  {
    if (itdmsls.Key().ShapeType() == TopAbs_EDGE)
    {
      NCollection_List<TopoDS_Shape> thedesc;
      theMap.Clear();
      for (itl.Initialize(itdmsls.Value()); itl.More(); itl.Next())
      {
        if (theMap.Add(MapW(itl.Value())))
        {
          thedesc.Append(MapW(itl.Value()));
        }
      }
      myMap(itdmsls.Key()) = thedesc;
    }
    else if (itdmsls.Key().IsSame(F))
    {
      myMap(F).Clear();
      for (itl.Initialize(thefaces); itl.More(); itl.Next())
      {
        myMap(F).Append(itl.Value());
      }
    }
    else
    {
      NCollection_List<TopoDS_Shape> thedesc;
      theMap.Clear();
      for (itl.Initialize(itdmsls.Value()); itl.More(); itl.Next())
      {
        if (theSubs.IsCopied(itl.Value()))
        {
          if (theSubs.Copy(itl.Value()).Extent() != 1)
          {
#ifdef OCCT_DEBUG
            std::cout << "Invalid number of descendant" << std::endl;
#endif
            return;
          }
          else
          {
            if (theMap.Add(theSubs.Copy(itl.Value()).First()))
            {
              thedesc.Append(theSubs.Copy(itl.Value()).First());
            }
          }
        }
        else if (theMap.Add(itl.Value()))
        {
          thedesc.Append(itl.Value());
        }
      }
      myMap(itdmsls.Key()) = thedesc;
    }
  }

  theMap.Clear();
  thefaces.Clear();
  for (itdmsls.Initialize(myMap); itdmsls.More(); itdmsls.Next())
  {
    for (itl.Initialize(itdmsls.Value()); itl.More(); itl.Next())
    {
      if (itl.Value().ShapeType() == TopAbs_FACE && theMap.Add(itl.Value()))
      {
        thefaces.Append(itl.Value());
      }
    }
  }
  LocOpe_BuildShape BS(thefaces);
  myResult = BS.Shape();
}

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

const TopoDS_Shape& LocOpe_SplitDrafts::Shape() const
{
  if (myResult.IsNull())
  {
    throw StdFail_NotDone();
  }
  return myResult;
}

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

const NCollection_List<TopoDS_Shape>& LocOpe_SplitDrafts::ShapesFromShape(
  const TopoDS_Shape& S) const
{
  if (myResult.IsNull())
  {
    throw StdFail_NotDone();
  }
  return myMap(S);
}

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

static bool NewPlane(const TopoDS_Face& F,
                     const gp_Dir&      Extr,
                     const gp_Pln&      Neutr,
                     const double       Ang,
                     gp_Pln&            Newpl,
                     gp_Ax1&            NormalF,
                     const bool         Modify)
{

  // Determination du nouveau plan incline
  occ::handle<Geom_Surface> S = BRep_Tool::Surface(F);
  if (S->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface))
  {
    S = occ::down_cast<Geom_RectangularTrimmedSurface>(S)->BasisSurface();
  }
  occ::handle<Geom_Plane> P = occ::down_cast<Geom_Plane>(S);
  if (P.IsNull())
  {
    return false;
  }

  gp_Pln Plorig = P->Pln();
  if (!Modify)
  {
    Newpl   = Plorig;
    NormalF = Newpl.Axis();
    if ((Newpl.Direct() && F.Orientation() == TopAbs_REVERSED)
        || (!Newpl.Direct() && F.Orientation() == TopAbs_FORWARD))
    {
      NormalF.Reverse();
    }
    return true;
  }

  gp_Ax1 Axe;
  double Theta;

  IntAna_QuadQuadGeo i2pl(Plorig, Neutr, Precision::Angular(), Precision::Confusion());

  if (i2pl.IsDone() && i2pl.TypeInter() == IntAna_Line)
  {
    gp_Lin LinInters = i2pl.Line(1);
    gp_Dir nx        = LinInters.Direction();
    NormalF          = Plorig.Axis();
    gp_Dir ny        = NormalF.Direction().Crossed(nx);
    double a         = Extr.Dot(nx);
    if (std::abs(a) <= 1 - Precision::Angular())
    {
      double             b = Extr.Dot(ny);
      double             c = Extr.Dot(NormalF.Direction());
      bool               direct(Plorig.Direct());
      TopAbs_Orientation Oris = F.Orientation();
      if ((direct && Oris == TopAbs_REVERSED) || (!direct && Oris == TopAbs_FORWARD))
      {
        b = -b;
        c = -c;
        NormalF.Reverse();
      }
      double denom = std::sqrt(1 - a * a);
      double Sina  = std::sin(Ang);
      if (denom > std::abs(Sina))
      {
        double phi    = std::atan2(b / denom, c / denom);
        double theta0 = std::acos(Sina / denom);
        Theta         = theta0 - phi;
        if (std::cos(Theta) < 0.)
        {
          Theta = -theta0 - phi;
        }
        Axe   = LinInters.Position();
        Newpl = Plorig.Rotated(Axe, Theta);
        return true;
      }
    }
  }
  std::cout << "fin newplane return standard_false" << std::endl;
  return false;
}

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

static void MakeFace(TopoDS_Face& F, NCollection_List<TopoDS_Shape>& ledg)
{

  // ledg est une liste d'edge

  BRep_Builder B;

  // Verification de l`existence des p-curves. Celles qui manquent
  // correspondent necessairement a des isos (et meme des iso u).

  double f, l;
  //  for (NCollection_List<TopoDS_Shape>::Iterator itl(ledg);
  NCollection_List<TopoDS_Shape>::Iterator itl(ledg);
  for (; itl.More(); itl.Next())
  {
    const TopoDS_Edge&        edg = TopoDS::Edge(itl.Value());
    occ::handle<Geom2d_Curve> C2d = BRep_Tool::CurveOnSurface(edg, F, f, l);
    if (C2d.IsNull())
    {
      BRep_Tool::Range(edg, f, l);
      TopoDS_Vertex V1, V2;
      TopExp::Vertices(edg, V1, V2);
      NCollection_List<TopoDS_Shape>::Iterator itl2;
      for (itl2.Initialize(ledg); itl2.More(); itl2.Next())
      {
        const TopoDS_Edge& edg2 = TopoDS::Edge(itl2.Value());
        if (edg2.IsSame(edg))
        {
          continue;
        }
        TopoDS_Vertex Vp1, Vp2;
        TopExp::Vertices(edg2, Vp1, Vp2);
        if (Vp1.IsSame(V1) || Vp2.IsSame(V1) || Vp1.IsSame(V2) || Vp2.IsSame(V2))
        {
          double                    f2, l2;
          occ::handle<Geom2d_Curve> C22d = BRep_Tool::CurveOnSurface(edg2, F, f2, l2);
          if (!C22d.IsNull())
          {
            gp_Pnt2d pt2d;
            if (Vp1.IsSame(V1))
            {
              pt2d = C22d->Value(f2);
              pt2d.SetY(pt2d.Y() - f);
            }
            else if (Vp2.IsSame(V1))
            {
              pt2d = C22d->Value(l2);
              pt2d.SetY(pt2d.Y() - f);
            }
            else if (Vp1.IsSame(V2))
            {
              pt2d = C22d->Value(f2);
              pt2d.SetY(pt2d.Y() - l);
            }
            else if (Vp2.IsSame(V2))
            {
              pt2d = C22d->Value(l2);
              pt2d.SetY(pt2d.Y() - l);
            }
            C2d = new Geom2d_Line(pt2d, gp::DY2d());
            B.UpdateEdge(edg, C2d, F, BRep_Tool::Tolerance(edg));
            break;
          }
        }
      }
      if (C2d.IsNull())
      {
        std::cout << "Ca merde violemment" << std::endl;
      }
    }
  }

  NCollection_List<TopoDS_Shape> lwires;
  bool                           alldone = ledg.IsEmpty();
  while (!alldone)
  {
    TopoDS_Wire Wnew;
    B.MakeWire(Wnew);
    TopoDS_Shape       aLocalShapeCur = ledg.First();
    const TopoDS_Edge& edg            = TopoDS::Edge(aLocalShapeCur);
    //    const TopoDS_Edge& edg = TopoDS::Edge(ledg.First());
    TopoDS_Vertex VFirst, VLast;
    if (edg.Orientation() == TopAbs_FORWARD)
    {
      TopExp::Vertices(edg, VFirst, VLast);
    }
    else
    {
      TopExp::Vertices(edg, VLast, VFirst);
    }
    B.Add(Wnew, edg);
    ledg.RemoveFirst();
    // on suppose VFirst et VLast non nuls
    bool wdone = (ledg.IsEmpty() || VFirst.IsSame(VLast));
    while (!wdone)
    {
      TopoDS_Vertex VF, VL;

      TopAbs_Orientation oredg = TopAbs_FORWARD;

      for (itl.Initialize(ledg); itl.More(); itl.Next())
      {
        const TopoDS_Edge& edg2        = TopoDS::Edge(itl.Value());
        TopoDS_Shape       aLocalShape = edg2.Oriented(TopAbs_FORWARD);
        TopExp::Vertices(TopoDS::Edge(aLocalShape), VF, VL);
        //	TopExp::Vertices(TopoDS::Edge(edg2.Oriented(TopAbs_FORWARD)),VF,VL);
        if (VF.IsSame(VLast))
        {
          VLast = VL;
          oredg = TopAbs_FORWARD;
          break;
        }
        else if (VL.IsSame(VFirst))
        {
          VFirst = VF;
          oredg  = TopAbs_FORWARD;
          break;
        }
        else if (VF.IsSame(VFirst))
        {
          VFirst = VL;
          oredg  = TopAbs_REVERSED;
          break;
        }
        else if (VL.IsSame(VLast))
        {
          VLast = VF;
          oredg = TopAbs_REVERSED;
          break;
        }
      }
      if (!itl.More())
      {
        wdone = true;
      }
      else
      {
        TopoDS_Shape aLocalShape = itl.Value().Oriented(oredg);
        B.Add(Wnew, TopoDS::Edge(aLocalShape));
        //	B.Add(Wnew,TopoDS::Edge(itl.Value().Oriented(oredg)));
        ledg.Remove(itl);
        wdone = (ledg.IsEmpty() || VFirst.IsSame(VLast));
      }
    }
    lwires.Append(Wnew);
    alldone = ledg.IsEmpty();
  }

  F.Orientation(TopAbs_FORWARD);
  for (itl.Initialize(lwires); itl.More(); itl.Next())
  {
    TopoDS_Shape aLocalShape = F.EmptyCopied();
    TopoDS_Face  NewFace     = TopoDS::Face(aLocalShape);
    //    TopoDS_Face NewFace = TopoDS::Face(F.EmptyCopied());
    B.Add(NewFace, itl.Value());
    GProp_GProps GP;
    BRepGProp::SurfaceProperties(NewFace, GP);
    if (GP.Mass() < 0)
    {
      itl.ChangeValue().Reverse();
    }
  }
  if (lwires.Extent() == 1)
  {
    B.Add(F, lwires.First());
  }
  else
  {
    std::cout << "Not yet implemented : nbwire >= 2" << std::endl;
  }
}

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

static bool Contains(const NCollection_List<TopoDS_Shape>& ll, const TopoDS_Shape& s)
{
  NCollection_List<TopoDS_Shape>::Iterator itl;
  for (itl.Initialize(ll); itl.More(); itl.Next())
  {
    if (itl.Value().IsSame(s))
    {
      return true;
    }
  }
  return false;
}

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

static TopoDS_Edge NewEdge(const TopoDS_Edge&               edg,
                           const TopoDS_Face&               F,
                           const occ::handle<Geom_Surface>& NewS,
                           const TopoDS_Vertex&             V1,
                           const TopoDS_Vertex&             V2)
{
  TopoDS_Edge               NewEdg;
  occ::handle<Geom_Surface> S1    = BRep_Tool::Surface(F);
  bool                      AppS1 = false;
  if (S1->DynamicType() != STANDARD_TYPE(Geom_Plane))
  {
    AppS1 = true;
  }

  GeomInt_IntSS i2s(S1, NewS, Precision::Confusion(), true, AppS1);
  if (!i2s.IsDone() || i2s.NbLines() <= 0)
  {
    return NewEdg;
  }

  BRep_Builder B;
  //  double pmin, Dist, DistMin;
  double            Dist2, Dist2Min;
  double            prmf = 0, prml = 0;
  GeomAdaptor_Curve TheCurve;

  int    i, k;
  gp_Pnt pvf = BRep_Tool::Pnt(V1);
  gp_Pnt pvl = BRep_Tool::Pnt(V2);
  for (i = 1; i <= i2s.NbLines(); i++)
  {
    TheCurve.Load(i2s.Line(i));
    Extrema_ExtPC myExtPC(pvf, TheCurve);

    if (myExtPC.IsDone())
    {
      gp_Pnt p1b, p2b;
      double thepmin = TheCurve.FirstParameter();
      myExtPC.TrimmedSquareDistances(Dist2Min, Dist2, p1b, p2b);
      if (Dist2 < Dist2Min && !TheCurve.IsPeriodic())
      {
        Dist2Min = Dist2;
        thepmin  = TheCurve.LastParameter();
      }
      for (k = 1; k <= myExtPC.NbExt(); k++)
      {
        Dist2 = myExtPC.SquareDistance(k);
        if (Dist2 < Dist2Min)
        {
          Dist2Min = Dist2;
          thepmin  = myExtPC.Point(k).Parameter();
        }
      }

      if (Dist2Min <= Precision::SquareConfusion())
      {
        prmf = thepmin;
        myExtPC.Perform(pvl);
        if (myExtPC.IsDone())
        {
          thepmin = TheCurve.LastParameter();
          myExtPC.TrimmedSquareDistances(Dist2, Dist2Min, p1b, p2b);
          if (Dist2 < Dist2Min && !TheCurve.IsClosed())
          {
            Dist2Min = Dist2;
            thepmin  = TheCurve.FirstParameter();
          }
          for (k = 1; k <= myExtPC.NbExt(); k++)
          {
            Dist2 = myExtPC.SquareDistance(k);
            if (Dist2 < Dist2Min)
            {
              Dist2Min = Dist2;
              thepmin  = myExtPC.Point(k).Parameter();
            }
          }

          if (Dist2Min <= Precision::SquareConfusion())
          {
            prml = thepmin;
            break;
          }
        }
      }
    }
  }

  if (i <= i2s.NbLines())
  {
    bool                           rev  = false;
    TopoDS_Vertex                  Vf   = V1;
    TopoDS_Vertex                  Vl   = V2;
    const occ::handle<Geom_Curve>& Cimg = i2s.Line(i);
    occ::handle<Geom2d_Curve>      Cimg2d;
    if (AppS1)
    {
      Cimg2d = i2s.LineOnS1(i);
    }

    if (Cimg->IsPeriodic())
    {

      double period = Cimg->Period();
      double imf    = Cimg->FirstParameter();
      double iml    = Cimg->LastParameter();

      double f, l;
      BRep_Tool::Range(edg, f, l);
      double delt  = l - f;
      double delt1 = std::abs(prml - prmf);
      double delt2 = std::abs(period - delt1);

      if (delt1 == 0 || delt2 == 0)
      {
        //	prmf = 0;
        //	prml = period;
        prmf = imf;
        prml = iml;
      }
      else
      {
        if (std::abs(delt1 - delt) > std::abs(delt2 - delt))
        {
          // le bon ecart est delt2...
          if (prml > prmf)
          {
            if (prml < iml)
            {
              prmf += period;
            }
            else
            {
              prml -= period;
            }
          }
          else
          {
            if (prmf < iml)
            {
              prml += period;
            }
            else
            {
              prmf -= period;
            }
          }
        }
        else if (std::abs(delt1 - delt) < std::abs(delt2 - delt))
        {
          if (prmf >= iml && prml >= iml)
          {
            prmf -= period;
            prml -= period;
          }
          else if (prmf <= imf && prml <= imf)
          {
            prmf += period;
            prml += period;
          }
        }
        else
        { // egalite; on priveligie l'ordre f,l
          if (prmf > prml)
          {
            prmf -= period;
          }
          if (prmf >= iml && prml >= iml)
          {
            prmf -= period;
            prml -= period;
          }
          else if (prmf <= imf && prml <= imf)
          {
            prmf += period;
            prml += period;
          }
        }
      }
#ifdef OCCT_DEBUG
      double ptol = Precision::PConfusion();
      if (prmf < imf - ptol || prmf > iml + ptol || prml < imf - ptol || prml > iml + ptol)
      {
        std::cout << "Ca ne va pas aller" << std::endl;
      }
#endif
    }

    if (S1->IsUPeriodic())
    {

      double speriod = S1->UPeriod();
      //      double f,l;
      gp_Pnt2d pf, pl;
      pf = Cimg2d->Value(prmf);
      pl = Cimg2d->Value(prml);

      double Uf   = pf.X();
      double Ul   = pl.X();
      double ptra = 0.0;

      double Ustart = std::min(Uf, Ul);
      while (Ustart < -Precision::PConfusion())
      {
        Ustart += speriod;
        ptra += speriod;
      }
      while (Ustart > speriod - Precision::PConfusion())
      {
        Ustart -= speriod;
        ptra -= speriod;
      }
      if (ptra != 0)
      {
        Cimg2d->Translate(gp_Vec2d(ptra, 0.));
      }
    }
    if (prmf < prml)
    {
      Vf.Orientation(TopAbs_FORWARD);
      Vl.Orientation(TopAbs_REVERSED);
    }
    else
    {
      Vf.Orientation(TopAbs_REVERSED);
      Vl.Orientation(TopAbs_FORWARD);
      rev = true;
    }

    B.MakeEdge(NewEdg, Cimg, Precision::Confusion());

    B.Add(NewEdg, Vf);
    B.Add(NewEdg, Vl);
    B.UpdateVertex(Vf, prmf, NewEdg, Precision::Confusion());
    B.UpdateVertex(Vl, prml, NewEdg, Precision::Confusion());
    if (AppS1)
    {
      B.UpdateEdge(NewEdg, Cimg2d, F, Precision::Confusion());
    }

    if (rev)
    {
      NewEdg.Orientation(TopAbs_REVERSED);
    }
  }
  return NewEdg;
}
