0023404: Create SquareConfusion function in Precision package for speed and convenience
[occt.git] / src / LocOpe / LocOpe_SplitDrafts.cxx
1 // Created on: 1996-10-02
2 // Created by: Jacques GOUSSARD
3 // Copyright (c) 1996-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
5 //
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
10 //
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13 //
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
20
21
22
23 #include <LocOpe_SplitDrafts.ixx>
24
25 #include <TopExp_Explorer.hxx>
26 #include <TopTools_DataMapIteratorOfDataMapOfShapeListOfShape.hxx>
27
28 #include <BRep_Builder.hxx>
29 #include <BRep_Tool.hxx>
30
31 #include <BRepTools_Substitution.hxx>
32
33 #include <LocOpe_WiresOnShape.hxx>
34 #include <LocOpe_Spliter.hxx>
35 #include <LocOpe_SplitShape.hxx>
36 #include <LocOpe_FindEdges.hxx>
37 #include <LocOpe_BuildShape.hxx>
38
39
40 #include <TopoDS_Vertex.hxx>
41 #include <TopoDS_Edge.hxx>
42
43 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
44 #include <TopTools_DataMapIteratorOfDataMapOfShapeShape.hxx>
45
46 #include <IntAna_QuadQuadGeo.hxx>
47 #include <IntCurveSurface_HInter.hxx>
48 #include <IntCurveSurface_IntersectionPoint.hxx>
49 #include <IntCurveSurface_IntersectionSegment.hxx>
50 #include <GeomInt_IntSS.hxx>
51 #include <Extrema_ExtPC.hxx>
52 #include <GeomAdaptor_HCurve.hxx>
53 #include <GeomAdaptor_HSurface.hxx>
54 #include <GeomAdaptor_Curve.hxx>
55 #include <GeomAdaptor_Surface.hxx>
56
57 #include <Geom_Surface.hxx>
58 #include <Geom_RectangularTrimmedSurface.hxx>
59 #include <Geom_Plane.hxx>
60 #include <Geom_Curve.hxx>
61 #include <Geom_TrimmedCurve.hxx>
62 #include <Geom_Line.hxx>
63
64 #include <Geom2d_Curve.hxx>
65 #include <Geom2d_Line.hxx>
66 #include <gp_Pnt2d.hxx>
67 #include <gp_Vec2d.hxx>
68
69 #include <GeomFill_Pipe.hxx>
70
71 #include <GProp_GProps.hxx>
72
73 #include <Standard_ConstructionError.hxx>
74
75 #include <TopoDS.hxx>
76 #include <TopExp.hxx>
77 #include <Precision.hxx>
78 #include <BRepGProp.hxx>
79 #include <gp.hxx>
80
81 static Standard_Boolean NewPlane(const TopoDS_Face&,
82                                  const gp_Dir&,
83                                  const gp_Pln&,
84                                  const Standard_Real,
85                                  gp_Pln&,
86                                  gp_Ax1&,
87                                  const Standard_Boolean);
88
89 static void MakeFace(TopoDS_Face&,
90                      TopTools_ListOfShape&);
91
92 static TopoDS_Edge  NewEdge(const TopoDS_Edge&,
93                             const TopoDS_Face&,
94                             const Handle(Geom_Surface)&,
95                             const TopoDS_Vertex&,
96                             const TopoDS_Vertex&);
97
98
99 static Standard_Boolean Contains(const TopTools_ListOfShape&,
100                                  const TopoDS_Shape&);
101
102
103 //=======================================================================
104 //function : Init
105 //purpose  : 
106 //=======================================================================
107
108 void LocOpe_SplitDrafts::Init(const TopoDS_Shape& S)
109 {
110   myShape = S;
111   myResult.Nullify();
112   myMap.Clear();
113 }
114
115 //=======================================================================
116 //function : Perform
117 //purpose  : 
118 //=======================================================================
119
120 void LocOpe_SplitDrafts::Perform(const TopoDS_Face& F,
121                                  const TopoDS_Wire& W,
122                                  const gp_Dir& Extr,
123                                  const gp_Pln& NPl,
124                                  const Standard_Real Angle)
125 {
126   Perform(F,W,Extr,NPl,Angle,Extr,NPl,Angle,Standard_True,Standard_False);
127 }
128
129
130 //=======================================================================
131 //function : Perform
132 //purpose  : 
133 //=======================================================================
134
135 void LocOpe_SplitDrafts::Perform(const TopoDS_Face& F,
136                                  const TopoDS_Wire& W,
137                                  const gp_Dir& Extrg,
138                                  const gp_Pln& NPlg,
139                                  const Standard_Real Angleg,
140                                  const gp_Dir& Extrd,
141                                  const gp_Pln& NPld,
142                                  const Standard_Real Angled,
143                                  const Standard_Boolean ModLeft,
144                                  const Standard_Boolean ModRight)
145
146 {
147   Standard_Integer j ;
148
149   myResult.Nullify();
150   myMap.Clear();
151   if (myShape.IsNull() || F.IsNull() || W.IsNull()) {
152     Standard_NullObject::Raise();
153   }    
154
155   if (!ModLeft && !ModRight) {
156     Standard_ConstructionError::Raise();
157   }
158
159 #ifdef DEB
160   TopAbs_Orientation OriF;
161 #else
162   TopAbs_Orientation OriF = TopAbs_FORWARD;
163 #endif
164   Standard_Boolean FinS = Standard_False;
165   TopExp_Explorer exp,exp2;
166   for (exp.Init(myShape,TopAbs_FACE); exp.More(); exp.Next()) {
167     const TopoDS_Shape& fac = exp.Current();
168     TopTools_ListOfShape thelist;
169     myMap.Bind(fac, thelist);
170     if (fac.IsSame(F)) {
171       OriF = fac.Orientation();
172       FinS = Standard_True;
173     }
174   }
175
176   if (!FinS) {
177     cout << "LocOpe_SplitDrafts:!Fins Standard_ConstructionError::Raise()" << endl;
178     Standard_ConstructionError::Raise();
179   }    
180
181   gp_Pln NewPlg,NewPld;
182   gp_Ax1 NormalFg,NormalFd;
183   TopoDS_Shape aLocalFace = F.Oriented(OriF);
184
185   if (!NewPlane(TopoDS::Face(aLocalFace),
186                 Extrg,NPlg,Angleg,NewPlg,NormalFg,ModLeft) ||
187       !NewPlane(TopoDS::Face(aLocalFace),
188                 Extrd,NPld,Angled,NewPld,NormalFd,ModRight)) {
189     //  if (!NewPlane(TopoDS::Face(F.Oriented(OriF)),
190     //          Extrg,NPlg,Angleg,NewPlg,NormalFg,ModLeft) ||
191 //      !NewPlane(TopoDS::Face(F.Oriented(OriF)),
192     //          Extrd,NPld,Angled,NewPld,NormalFd,ModRight)) {
193     return;
194   }
195
196
197   TopTools_ListIteratorOfListOfShape itl;
198   BRep_Builder B;
199
200   Handle(Geom_Surface) NewSg = new Geom_Plane(NewPlg);
201   Handle(Geom_Surface) NewSd = new Geom_Plane(NewPld);
202   Handle(Geom_Line) theLinePipe = new Geom_Line(NormalFg); // ou NormalFd
203   GeomInt_IntSS i2s(NewSg,NewSd,Precision::Confusion());
204
205   TopTools_MapOfShape theMap;
206   Handle(GeomAdaptor_HCurve) HAC = new GeomAdaptor_HCurve;
207   GeomAdaptor_Curve AC;
208   Handle(GeomAdaptor_HSurface) HAS = new GeomAdaptor_HSurface;
209   GeomAdaptor_Surface AS;
210   IntCurveSurface_HInter intcs;
211
212   TopoDS_Wire theW = W;
213   if (i2s.IsDone() && i2s.NbLines() > 0) {
214     // on split le wire" << endl;
215
216     GeomFill_Pipe thePipe;
217     thePipe.GenerateParticularCase(Standard_True);
218     thePipe.Init(theLinePipe,i2s.Line(1));
219     thePipe.Perform(Standard_True);
220
221     Handle(Geom_Surface) Spl = thePipe.Surface();
222     AS.Load(Spl);
223     HAS->Set(AS);
224     
225     LocOpe_SplitShape splw(W);
226
227     for (exp.Init(W,TopAbs_EDGE); exp.More(); exp.Next()) {
228       const TopoDS_Edge& edg = TopoDS::Edge(exp.Current());
229       if (theMap.Add(edg)) {
230         TopLoc_Location Loc;
231         Standard_Real f,l;
232         Handle(Geom_Curve) C = BRep_Tool::Curve(edg,f,l);
233         AC.Load(C);
234         HAC->Set(AC);
235         intcs.Perform(HAC,HAS);
236         if (!intcs.IsDone()) {
237           continue; // voir ce qu`on peut faire de mieux
238         }
239
240         if (intcs.NbSegments() >= 2) {
241           continue; // Not yet implemented...and probably never"
242         }
243
244         if (intcs.NbSegments() == 1) {
245           const IntCurveSurface_IntersectionPoint& P1 = 
246             intcs.Segment(1).FirstPoint();
247           const IntCurveSurface_IntersectionPoint& P2 = 
248             intcs.Segment(1).SecondPoint();
249           const gp_Pnt& pf = P1.Pnt();
250           const gp_Pnt& pl = P2.Pnt();
251           TopoDS_Vertex Vf,Vl;
252           TopExp::Vertices(edg,Vf,Vl);
253           gp_Pnt Pf = BRep_Tool::Pnt(Vf);
254           gp_Pnt Pl = BRep_Tool::Pnt(Vl);
255           Standard_Real Tolf = BRep_Tool::Tolerance(Vf);
256           Standard_Real Toll = BRep_Tool::Tolerance(Vl);
257           Tolf *= Tolf;
258           Toll *= Toll;
259
260           Standard_Real dff = pf.SquareDistance(Pf);
261           Standard_Real dfl = pf.SquareDistance(Pl);
262           Standard_Real dlf = pl.SquareDistance(Pf);
263           Standard_Real dll = pl.SquareDistance(Pl);
264
265           if ((dff <= Tolf && dll <= Toll) ||
266               (dlf <= Tolf && dfl <= Toll)) {
267             continue;
268           }
269           else {
270             // on segmente edg en pf et pl
271             TopoDS_Vertex Vnewf,Vnewl;
272             B.MakeVertex(Vnewf,pf,Precision::Confusion());
273             B.MakeVertex(Vnewl,pl,Precision::Confusion());
274             if (P1.W() >= f && P1.W() <= l &&
275                 P2.W() >= f && P2.W() <= l) {
276               splw.Add(Vnewf,P1.W(),edg);
277               splw.Add(Vnewl,P2.W(),edg);
278             }
279             else {
280               continue;
281             }
282           }
283         }
284         else if (intcs.NbPoints() != 0) {
285           TopoDS_Vertex Vf,Vl;
286           TopExp::Vertices(edg,Vf,Vl);
287           gp_Pnt Pf = BRep_Tool::Pnt(Vf);
288           gp_Pnt Pl = BRep_Tool::Pnt(Vl);
289           Standard_Real Tolf = BRep_Tool::Tolerance(Vf);
290           Standard_Real Toll = BRep_Tool::Tolerance(Vl);
291           Tolf *= Tolf;
292           Toll *= Toll;
293
294           for (Standard_Integer i = 1; i <= intcs.NbPoints(); i++) {
295             const IntCurveSurface_IntersectionPoint& Pi = intcs.Point(i);
296             const gp_Pnt& pi = Pi.Pnt();
297             Standard_Real dif = pi.SquareDistance(Pf);
298             Standard_Real dil = pi.SquareDistance(Pl);
299             if (dif <= Tolf) {
300             }
301             else if (dil <= Toll) {
302             }
303             else {
304               if (Pi.W() >= f && Pi.W() <= l) {
305                 TopoDS_Vertex Vnew;
306                 B.MakeVertex(Vnew,pi,Precision::Confusion());
307                 splw.Add(Vnew,Pi.W(),edg);
308               }
309             }
310           }
311         }
312       }
313     }
314
315     const TopTools_ListOfShape& lres = splw.DescendantShapes(W);
316     if (lres.Extent() != 1) {
317       return;
318     }
319
320     if (!W.IsSame(lres.First())) {
321       theW.Nullify();
322       theW = TopoDS::Wire(lres.First());
323     }
324
325     for (exp.ReInit(); exp.More(); exp.Next()) {
326       if (!myMap.IsBound(exp.Current())) {
327         TopTools_ListOfShape thelist1;
328         myMap.Bind(exp.Current(), thelist1);
329         for (itl.Initialize(splw.DescendantShapes(exp.Current())); 
330              itl.More(); itl.Next()) {
331           myMap(exp.Current()).Append(itl.Value());
332         }
333         for (exp2.Init(exp.Current(),TopAbs_VERTEX);exp2.More();exp2.Next()) {
334           if (!myMap.IsBound(exp2.Current())) {
335             TopTools_ListOfShape thelist2;
336             myMap.Bind(exp2.Current(), thelist2);
337             myMap(exp2.Current()).Append(exp2.Current());
338           }
339         }
340       }
341     }
342   }
343   else {
344     for (exp.Init(W,TopAbs_EDGE); exp.More(); exp.Next()) {
345       if (!myMap.IsBound(exp.Current())) {
346         TopTools_ListOfShape thelist3;
347         myMap.Bind(exp.Current(), thelist3);
348         myMap(exp.Current()).Append(exp.Current());
349         for (exp2.Init(exp.Current(),TopAbs_VERTEX);exp2.More();exp2.Next()) {
350           if (!myMap.IsBound(exp2.Current())) {
351             TopTools_ListOfShape thelist4;
352             myMap.Bind(exp2.Current(), thelist4);
353             myMap(exp2.Current()).Append(exp2.Current());
354           }
355         }
356       }
357     }
358   }
359
360   // On split la face par le wire
361   
362   Handle(LocOpe_WiresOnShape) WonS = new LocOpe_WiresOnShape(myShape);
363   LocOpe_Spliter Spls(myShape);
364   WonS->Bind(theW,F);
365
366 // JAG Le code suivant marchera apres integration de thick0
367 //  LocOpe_FindEdges fined(W,F);
368 //  for (fined.InitIterator(); fined.More(); fined.Next()) {
369 //    WonS->Bind(fined.EdgeFrom(),fined.EdgeTo());
370 //  }
371
372   Spls.Perform(WonS);
373   if (!Spls.IsDone()) {
374     return;
375   }
376
377   TopoDS_Shape Res = Spls.ResultingShape();
378   const TopTools_ListOfShape& theLeft = Spls.DirectLeft();
379
380   // Descendants
381   for (exp.Init(myShape,TopAbs_FACE); exp.More(); exp.Next()) {
382     const TopoDS_Shape& fac = exp.Current();
383     for (itl.Initialize(Spls.DescendantShapes(fac)); itl.More(); itl.Next()) {
384       myMap(fac).Append(itl.Value());
385     }
386   }
387
388   TopTools_DataMapOfShapeShape MapW;
389   for (exp.Init(theW,TopAbs_EDGE); exp.More(); exp.Next()) {
390     if (!MapW.IsBound(exp.Current())) {
391       MapW.Bind(exp.Current(),TopoDS_Shape());
392       for (exp2.Init(exp.Current(),TopAbs_VERTEX); exp2.More(); exp2.Next()) {
393         if (!MapW.IsBound(exp2.Current())) {
394           MapW.Bind(exp2.Current(),TopoDS_Shape());
395         }
396
397       }
398     }
399   }
400   
401
402
403   TopTools_IndexedDataMapOfShapeListOfShape theMapEF;
404   TopExp::MapShapesAndAncestors(Res,TopAbs_EDGE,TopAbs_FACE,theMapEF);
405
406   // On stocke les geometries potentiellement generees par les edges
407   TopTools_IndexedDataMapOfShapeShape MapEV; // genere
408   TopTools_DataMapOfShapeListOfShape MapSg,MapSd; // image a gauche et a droite
409
410   Standard_Integer Nbedges,index;
411   for (itl.Initialize(myMap(F)); itl.More(); itl.Next()) {
412     const TopoDS_Shape& fac = TopoDS::Face(itl.Value());
413     for (exp.Init(fac,TopAbs_EDGE); exp.More(); exp.Next()) {
414       const TopoDS_Shape& edg = exp.Current();
415       if (MapEV.FindIndex(edg) != 0) {
416         continue;
417       }
418       if (MapW.IsBound(edg)) { // edge du wire initial
419         TopLoc_Location Loc;
420         Standard_Real f,l;
421         Handle(Geom_Curve) C = BRep_Tool::Curve(TopoDS::Edge(edg),Loc,f,l);
422         if (C->DynamicType() == STANDARD_TYPE(Geom_TrimmedCurve)) {
423           C = Handle(Geom_TrimmedCurve)::DownCast(C)->BasisCurve();
424         }
425         C = Handle(Geom_Curve)::
426           DownCast(C->Transformed(Loc.Transformation()));
427
428         GeomFill_Pipe thePipe;
429         thePipe.GenerateParticularCase(Standard_True);
430         thePipe.Init(theLinePipe,C);
431         thePipe.Perform(Standard_True);
432
433         Handle(Geom_Surface) thePS = thePipe.Surface();
434         if (thePS->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
435           thePS = Handle(Geom_RectangularTrimmedSurface)::DownCast(thePS)
436             ->BasisSurface();
437         }
438
439         TopoDS_Face NewFace;
440         B.MakeFace(NewFace,thePS,Precision::Confusion());
441         MapEV.Add(edg,NewFace);
442       }
443       else { // on recupere la face.
444         index = theMapEF.FindIndex(edg);
445         if (theMapEF(index).Extent() != 2) {
446           return; // NotDone
447         }
448         TopoDS_Face theFace;
449         if (theMapEF(index).First().IsSame(fac)) {
450           MapEV.Add(edg,theMapEF(index).Last());
451         }
452         else {
453           MapEV.Add(edg,theMapEF(index).First());
454         }
455       }
456     }
457   }
458
459
460   TopTools_DataMapOfShapeShape MapSonS;
461
462   Nbedges = MapEV.Extent();
463   for (index = 1; index <= Nbedges; index++) {
464     for (exp.Init(MapEV.FindKey(index),TopAbs_VERTEX); 
465          exp.More(); exp.Next()) {
466       const TopoDS_Vertex& vtx = TopoDS::Vertex(exp.Current());
467       if (MapEV.FindIndex(vtx)!= 0) {
468         continue;
469       }
470
471       // Localisation du vertex :
472       //    - entre 2 edges d`origine : on recupere l`edge qui n`est 
473       //                                pas dans F
474       //    - entre 2 edges du wire   : droite
475       //    - mixte                   : intersection de surfaces
476       for ( j = 1; j<=Nbedges; j++) {
477         if (j == index) {
478           continue;
479         }
480         for (exp2.Init(MapEV.FindKey(j),TopAbs_VERTEX);
481              exp2.More(); exp2.Next()) {
482           const TopoDS_Shape& vtx2 = exp2.Current();
483           if (vtx2.IsSame(vtx)) {
484             break;
485           }
486         }
487         if (exp2.More()) {
488           break;
489         }
490       }
491       Standard_Integer Choice = 0;
492       const TopoDS_Shape& edg1 = MapEV.FindKey(index);
493       TopoDS_Shape edg2;
494       if (j <= Nbedges) {
495         edg2 = MapEV.FindKey(j);
496       }
497       else {
498         edg2 = edg1;
499       }
500       if (MapW.IsBound(edg1)) {
501         if (j>Nbedges) { // doit correspondre a edge ferme
502           Choice = 2; // droite
503         }
504         else if (MapW.IsBound(MapEV.FindKey(j))) {
505           Choice = 2; // droite
506         }
507         else {
508           Choice = 3; // mixte
509         }
510       }
511       else {
512         if (j>Nbedges) { // doit correspondre a edge ferme
513           Choice = 1; // edge a retrouver
514         }
515         else if (!MapW.IsBound(MapEV.FindKey(j))) {
516           Choice = 1; // edge a retrouver
517         }
518         else {
519           Choice = 3; // mixte
520         }
521       }
522       Handle(Geom_Curve) Newc;
523       Handle(Geom2d_Curve) newCs1,newCs2;
524       Standard_Real knownp=0;
525       TopoDS_Edge Ebind;
526       switch (Choice) {
527       case 1:
528         {
529           for (exp2.Init(Res,TopAbs_EDGE); exp2.More(); exp2.Next()) {
530             if (exp2.Current().IsSame(edg1) || exp2.Current().IsSame(edg2)) {
531               continue;
532             }
533 //          for (TopExp_Explorer exp3(exp2.Current().Oriented(TopAbs_FORWARD),
534             TopExp_Explorer exp3(exp2.Current().Oriented(TopAbs_FORWARD),
535                                       TopAbs_VERTEX) ;
536             for ( ; exp3.More(); exp3.Next())  {
537               if (exp3.Current().IsSame(vtx)) {
538                 break;
539               }
540             }
541             if (exp3.More()) {
542               break;
543             }
544           }
545           if (exp2.More()) {
546             Standard_Real f,l;
547             TopLoc_Location Loc;
548             Newc = BRep_Tool::Curve(TopoDS::Edge(exp2.Current()),Loc,f,l);
549             Newc = Handle(Geom_Curve)::DownCast
550               (Newc->Transformed(Loc.Transformation()));
551             Ebind = TopoDS::Edge(exp2.Current());
552             knownp = BRep_Tool::Parameter(vtx,Ebind);
553           }
554           else { // droite ??? il vaudrait mieux sortir
555             return;
556
557 //          gp_Lin theLine(NormalFg);
558 //          theLine.Translate(NormalF.Location(),BRep_Tool::Pnt(vtx));
559 //          Newc = new Geom_Line(theLine);
560 //          knownp = 0.;
561           }
562         }
563         break;
564       case 2:
565         {
566           gp_Lin theLine(NormalFg);
567           theLine.Translate(NormalFg.Location(),BRep_Tool::Pnt(vtx));
568           Newc = new Geom_Line(theLine);
569           knownp = 0.;
570         }
571         break;
572       case 3:
573         {
574           const TopoDS_Face& F1 = TopoDS::Face(MapEV.FindFromKey(edg1));
575           const TopoDS_Face& F2 = TopoDS::Face(MapEV.FindFromKey(edg2));
576           Handle(Geom_Surface) S1 = BRep_Tool::Surface(F1);
577           Handle(Geom_Surface) S2 = BRep_Tool::Surface(F2);
578           Standard_Boolean AppS1 = Standard_False;
579           Standard_Boolean AppS2 = Standard_False;
580           if (S1->DynamicType() != STANDARD_TYPE(Geom_Plane)) {
581             AppS1 = Standard_True;
582           }
583           if (S2->DynamicType() != STANDARD_TYPE(Geom_Plane)) {
584             AppS2 = Standard_True;
585           }
586           i2s.Perform(S1,S2,Precision::Confusion(),Standard_True,AppS1,AppS2);
587           if (!i2s.IsDone() || i2s.NbLines() <= 0) {
588             return;
589           }
590         
591           Standard_Real pmin=0, Dist2, Dist2Min, Glob2Min = RealLast();
592           GeomAdaptor_Curve TheCurve;
593
594           Standard_Integer i,imin,k;
595           gp_Pnt pv = BRep_Tool::Pnt(vtx);
596           imin = 0;
597           for (i=1; i<= i2s.NbLines(); i++) {
598             TheCurve.Load(i2s.Line(i));
599             Extrema_ExtPC myExtPC(pv,TheCurve);
600
601             if (myExtPC.IsDone()) {
602               gp_Pnt p1b,p2b;
603               Standard_Real thepmin = TheCurve.FirstParameter();
604               myExtPC.TrimmedSquareDistances(Dist2Min,Dist2,p1b,p2b);
605               if (Dist2 < Dist2Min) {
606                 thepmin = TheCurve.LastParameter();
607               }
608               for (k=1; k<=myExtPC.NbExt(); k++) {
609                 Dist2 = myExtPC.SquareDistance(k);
610                 if (Dist2 < Dist2Min) {
611                   Dist2Min = Dist2;
612                   thepmin = myExtPC.Point(k).Parameter();
613                 }
614               }
615               
616               if (Dist2Min  < Glob2Min) {
617                 Glob2Min = Dist2Min;
618                 pmin = thepmin;
619                 imin = i;
620               }
621             }
622           }
623           if (imin == 0) {
624             return;
625           }
626           
627           Newc = i2s.Line(imin);
628           knownp = pmin;
629           if (AppS1) {
630             newCs1 = i2s.LineOnS1(imin);
631           }
632           if (AppS2) {
633             newCs2 = i2s.LineOnS2(imin);
634           }
635         }
636         break;
637       }
638
639
640       // Determination des vertex par intersection sur Plg ou/et Pld
641
642       AC.Load(Newc);
643       HAC->Set(AC);
644       Standard_Integer nbfois = 2;
645       TopoDS_Vertex vtx1,vtx2;
646       Standard_Real p1=0,p2=0;
647       Standard_Boolean IsLeft=Standard_False;
648       if (Choice == 1) { 
649         // edge retrouve : on ne fait qu`une seule intersection
650         // il faut utiliser Plg ou Pld
651
652         Standard_Integer indedgf = theMapEF.FindIndex(edg1);
653         for (itl.Initialize(theMapEF(indedgf)); itl.More(); itl.Next()) {
654           if (Contains(myMap(F),itl.Value())) {
655             if (Contains(theLeft,itl.Value())) {
656               AS.Load(NewSg);
657               IsLeft = Standard_True;
658             }
659             else {
660               AS.Load(NewSd);
661               IsLeft = Standard_False;
662             }
663             
664             nbfois = 1;
665             vtx2 = vtx;
666             p2 = knownp;
667             break;
668           }
669         }
670         if (!itl.More()) {
671           cout << "LocOpe_SplitDrafts: betite probleme "<< endl;
672           return;
673         }
674
675       }
676       else {
677         AS.Load(NewSg);
678       }
679
680       for (Standard_Integer it = 1; it<=nbfois; it++) {
681         if (it == 2) {
682           AS.Load(NewSd);
683         }
684         HAS->Set(AS);
685
686         intcs.Perform(HAC,HAS);
687         if (!intcs.IsDone()) {
688           return; // voir ce qu`on peut faire de mieux
689         }
690         Standard_Integer imin = 1;
691         Standard_Real delta = Abs(knownp - intcs.Point(1).W());
692         for (Standard_Integer i = 2;  i<= intcs.NbPoints(); i++) {
693           Standard_Real newdelta =  Abs(knownp - intcs.Point(i).W());
694           if (newdelta < delta) {
695             imin = i;
696             delta = newdelta;
697           }
698         }
699         if (it == 1) {
700           B.MakeVertex(vtx1,intcs.Point(imin).Pnt(),Precision::Confusion());
701           p1 = intcs.Point(imin).W();
702           knownp = p1;
703         }
704         else {
705           B.MakeVertex(vtx2,intcs.Point(imin).Pnt(),Precision::Confusion());
706           p2 = intcs.Point(imin).W();
707         }
708       }
709       if (Abs(p1-p2) > Precision::PConfusion()) {
710         TopoDS_Edge NewEdge;
711         B.MakeEdge(NewEdge,Newc,Precision::Confusion());
712         if (p1 < p2) {
713           B.Add(NewEdge,vtx1.Oriented(TopAbs_FORWARD));
714           B.Add(NewEdge,vtx2.Oriented(TopAbs_REVERSED));
715         }
716         else {
717           B.Add(NewEdge,vtx1.Oriented(TopAbs_REVERSED));
718           B.Add(NewEdge,vtx2.Oriented(TopAbs_FORWARD));
719         }
720         B.UpdateVertex(vtx1,p1,NewEdge,Precision::Confusion());
721         B.UpdateVertex(vtx2,p2,NewEdge,Precision::Confusion());
722         if (!newCs1.IsNull()) {
723           B.UpdateEdge(NewEdge,newCs1,
724                        TopoDS::Face(MapEV.FindFromKey(edg1)),
725                        Precision::Confusion());
726         }
727
728         if (!newCs2.IsNull()) {
729           B.UpdateEdge(NewEdge,newCs2,
730                        TopoDS::Face(MapEV.FindFromKey(edg2)),
731                        Precision::Confusion());
732         }
733
734       
735         MapEV.Add(vtx,NewEdge);
736
737         if (Choice == 1) {
738           TopoDS_Shape aLocalEdge = Ebind.EmptyCopied();
739           TopoDS_Edge NE = TopoDS::Edge(aLocalEdge);
740 //        TopoDS_Edge NE = TopoDS::Edge(Ebind.EmptyCopied());
741           for (exp2.Init(Ebind,TopAbs_VERTEX); exp2.More(); exp2.Next()) {
742             const TopoDS_Vertex& thevtx = TopoDS::Vertex(exp2.Current());
743             if (thevtx.IsSame(vtx)) {
744               B.Add(NE,vtx1.Oriented(thevtx.Orientation()));
745               B.UpdateVertex(vtx1,p1,NE,Precision::Confusion());
746             }
747             else {
748               B.Add(NE,thevtx);
749               Standard_Real theprm = BRep_Tool::Parameter(thevtx,Ebind);
750               B.UpdateVertex(thevtx,theprm,NE,BRep_Tool::Tolerance(thevtx));
751             }
752           }
753           MapSonS.Bind(Ebind,NE.Oriented(TopAbs_FORWARD));
754           if (IsLeft) {
755             TopTools_ListOfShape thelist5;
756             MapSg.Bind(vtx, thelist5);
757             MapSg(vtx).Append(vtx1);
758           }
759           else {
760             TopTools_ListOfShape thelist6;
761             MapSd.Bind(vtx, thelist6);
762             MapSd(vtx).Append(vtx1);
763           }
764         }
765         else {
766           TopTools_ListOfShape thelist7, thelist8;
767           MapSg.Bind(vtx, thelist7);
768           MapSd.Bind(vtx, thelist8);
769           MapSg(vtx).Append(vtx1);
770           MapSd(vtx).Append(vtx2);
771         }
772       }
773       else {
774         MapEV.Add(vtx,vtx2); // on peut avoir vtx2 = vtx si choix == 1
775         if (Choice == 1) {
776           if (IsLeft) {
777             TopTools_ListOfShape thelist9;
778             MapSg.Bind(vtx, thelist9);
779             MapSg(vtx).Append(vtx);
780           }
781           else {
782             TopTools_ListOfShape thelist10;
783             MapSd.Bind(vtx, thelist10);
784             MapSd(vtx).Append(vtx);
785           }
786         }
787         else {
788           TopTools_ListOfShape thelist11, thelist12;
789           MapSg.Bind(vtx, thelist11);
790           MapSd.Bind(vtx, thelist12);
791           MapSg(vtx).Append(vtx2);
792           MapSd(vtx).Append(vtx2);
793         }
794       }
795     }
796   }
797
798
799   theMap.Clear();
800   for (exp.Init(theW,TopAbs_EDGE); exp.More(); exp.Next()) {
801     const TopoDS_Edge& edg = TopoDS::Edge(exp.Current());
802     if (!theMap.Add(edg)) { // precaution sans doute inutile...
803       continue;
804     }
805     Standard_Integer indedg = MapEV.FindIndex(edg);
806     TopoDS_Face& GenF = TopoDS::Face(MapEV(indedg));
807     TopTools_ListOfShape thelist13, thelist14;
808     MapSg.Bind(edg, thelist13);  // genere a gauche
809     MapSd.Bind(edg, thelist14);  // genere a droite
810     TopoDS_Vertex Vf,Vl;
811     TopoDS_Shape aLocalEdge = edg.Oriented(TopAbs_FORWARD);
812     TopExp::Vertices(TopoDS::Edge(aLocalEdge),Vf,Vl);
813 //    TopExp::Vertices(TopoDS::Edge(edg.Oriented(TopAbs_FORWARD)),Vf,Vl);
814     TopoDS_Shape Gvf = MapEV.FindFromKey(Vf); 
815     TopoDS_Shape Gvl = MapEV.FindFromKey(Vl); 
816
817 /* Le code suivant est OK. On essaie de l`ameliorer
818
819     if (Gvf.ShapeType() == TopAbs_VERTEX &&
820         Gvl.ShapeType() == TopAbs_VERTEX) {
821       // en fait on doit pouvoir avoir 1 face a 2 cotes...
822       if (Gvf.IsSame(Vf)) {
823         MapW(edg) = edg;
824         MapSg(edg).Append(edg.Oriented(TopAbs_FORWARD));
825         MapSd(edg).Append(edg.Oriented(TopAbs_FORWARD));
826       }
827       else {
828         TopoDS_Edge NewEdg = NewEdge(edg,
829                                      GenF,NewSg,
830                                      TopoDS::Vertex(Gvf),
831                                      TopoDS::Vertex(Gvl));
832         if (NewEdg.IsNull()) {
833           return;
834         }
835         MapW(edg) = NewEdg;
836         MapSg(edg).Append(NewEdg);
837         MapSd(edg).Append(NewEdg);
838       }
839     }
840     else if (Gvf.ShapeType() == TopAbs_VERTEX  ||
841              Gvl.ShapeType() == TopAbs_VERTEX) {      // face triangulaire
842       TopoDS_Vertex Vfd,Vld,Vfg,Vlg;
843       if (Gvf.ShapeType() == TopAbs_VERTEX) {
844         Vfg = TopoDS::Vertex(Gvf);
845         Vfd = Vfg;
846         Vlg = TopoDS::Vertex(MapSg(Vl).First());
847         Vld = TopoDS::Vertex(MapSd(Vl).First());
848       }
849       else {
850         Vlg = TopoDS::Vertex(Gvl);
851         Vld = Vlg;
852         Vfg = TopoDS::Vertex(MapSg(Vf).First());
853         Vfd = TopoDS::Vertex(MapSd(Vf).First());
854       }
855
856       TopoDS_Edge NewEdgg = NewEdge(edg,GenF,NewSg,Vfg,Vlg);
857       if (NewEdgg.IsNull()) {
858         return;
859       }
860
861       TopoDS_Edge NewEdgd = NewEdge(edg,GenF,NewSd,Vfd,Vld);
862       if (NewEdgg.IsNull()) {
863         return;
864       }
865       MapSg(edg).Append(NewEdgg);
866       MapSd(edg).Append(NewEdgd);
867
868       TopTools_ListOfShape theedges;
869       theedges.Append(NewEdgg);
870       theedges.Append(NewEdgd);
871       if (Gvf.ShapeType() == TopAbs_EDGE) {
872         theedges.Append(Gvf);
873       }
874       else {//if (Gvl.ShapeType() == TopAbs_EDGE) {
875         theedges.Append(Gvl);
876       }
877       MakeFace(GenF,theedges);
878       MapW(edg) = GenF;
879     }
880     else {
881       // une face a 4 cotes
882       TopoDS_Vertex Vfd,Vld,Vfg,Vlg;
883
884       Vfg = TopoDS::Vertex(MapSg(Vf).First());
885       Vfd = TopoDS::Vertex(MapSd(Vf).First());
886       Vlg = TopoDS::Vertex(MapSg(Vl).First());
887       Vld = TopoDS::Vertex(MapSd(Vl).First());
888       
889       TopoDS_Vertex VVf1,VVl1,VVf2,VVl2;
890       TopExp::Vertices(TopoDS::Edge(Gvf.Oriented(TopAbs_FORWARD)),VVf1,VVl1);
891       TopExp::Vertices(TopoDS::Edge(Gvl.Oriented(TopAbs_FORWARD)),VVf2,VVl2);
892
893       TopoDS_Edge NewEdgg = NewEdge(edg,GenF,NewSg,Vfg,Vlg);
894       if (NewEdgg.IsNull()) {
895         return;
896       }
897       
898       TopoDS_Edge NewEdgd = NewEdge(edg,GenF,NewSd,Vfd,Vld);
899       if (NewEdgd.IsNull()) {
900         return;
901       }
902
903       if ((VVf1.IsSame(Vfg) && VVf2.IsSame(Vlg)) ||
904           (VVf1.IsSame(Vfd) && VVf2.IsSame(Vld))) {
905         // 4 cotes
906         MapSg(edg).Append(NewEdgg);
907         MapSd(edg).Append(NewEdgd);
908         
909         TopTools_ListOfShape theedges;
910         theedges.Append(NewEdgg);
911         theedges.Append(NewEdgd);
912         theedges.Append(Gvf);
913         theedges.Append(Gvl);
914         
915         MakeFace(GenF,theedges);
916         MapW(edg) = GenF;
917       }
918       else {
919 #ifdef DEB
920         cout << "Pb d'analyse" << endl;
921 #endif
922         return;
923       }
924     }
925 */
926     // nouveau code
927
928     TopoDS_Vertex Vfd,Vld,Vfg,Vlg;
929     if (Gvf.ShapeType() == TopAbs_VERTEX) {
930       Vfg = TopoDS::Vertex(Gvf);
931       Vfd = Vfg;
932     }
933     else {
934       Vfg = TopoDS::Vertex(MapSg(Vf).First());
935       Vfd = TopoDS::Vertex(MapSd(Vf).First());
936     }
937     if (Gvl.ShapeType() == TopAbs_VERTEX) {
938       Vlg = TopoDS::Vertex(Gvl);
939       Vld = Vlg;
940     }
941     else {
942       Vlg = TopoDS::Vertex(MapSg(Vl).First());
943       Vld = TopoDS::Vertex(MapSd(Vl).First());
944     }
945
946     TopoDS_Edge NewEdgg = NewEdge(edg,GenF,NewSg,Vfg,Vlg);
947     if (NewEdgg.IsNull()) {
948       return;
949     }
950     
951     TopoDS_Edge NewEdgd = NewEdge(edg,GenF,NewSd,Vfd,Vld);
952     if (NewEdgg.IsNull()) {
953       return;
954     }
955
956     Standard_Boolean isedg = Standard_False;
957     if (Gvf.ShapeType() == TopAbs_VERTEX &&
958         Gvl.ShapeType() == TopAbs_VERTEX) {
959       // edg ou face a 2 cotes
960
961       // Comparaison NewEdgg et NewEdgd
962       Standard_Real fg,lg,fd,ld;
963       Handle(Geom_Curve) Cg = BRep_Tool::Curve(NewEdgg,fg,lg);
964       Handle(Geom_Curve) Cd = BRep_Tool::Curve(NewEdgd,fd,ld);
965       Standard_Real prmg = (fg+lg)/2.;
966       Standard_Real prmd = (fd+ld)/2.;
967       gp_Pnt pg = Cg->Value(prmg);
968       gp_Pnt pd = Cd->Value(prmd);
969       Standard_Real Tol = Max(BRep_Tool::Tolerance(NewEdgg),
970                               BRep_Tool::Tolerance(NewEdgg));
971       if (pg.SquareDistance(pd) <= Tol*Tol) {
972         isedg = Standard_True;
973         // raffinement pour essayer de partager l`edge de depart...
974         Standard_Boolean modified = Standard_True;
975         if (Gvf.IsSame(Vf) && Gvl.IsSame(Vl)) {
976           // Comparaison avec l`edge de depart
977           Cd = BRep_Tool::Curve(edg,fd,ld);
978           prmd = (fd+ld)/2.;
979           pd = Cd->Value(prmd);
980           Tol = Max(BRep_Tool::Tolerance(NewEdgg),
981                     BRep_Tool::Tolerance(edg));
982           if (pg.SquareDistance(pd) <= Tol*Tol) {
983             modified = Standard_False;
984           }
985         }
986
987         if (!modified) {
988           MapW(edg) = edg;
989           MapSg(edg).Append(edg);
990           MapSd(edg).Append(edg);
991         }
992         else {
993           MapW(edg) = NewEdgg;
994           MapSg(edg).Append(NewEdgg);
995           MapSd(edg).Append(NewEdgg);
996         }
997       }
998     }
999
1000     if (!isedg) {
1001       // face a 2 ou 3 ou 4 cotes
1002       MapSg(edg).Append(NewEdgg);
1003       MapSd(edg).Append(NewEdgd);
1004       
1005       TopTools_ListOfShape theedges;
1006       theedges.Append(NewEdgg);
1007       theedges.Append(NewEdgd);
1008       if (Gvf.ShapeType() == TopAbs_EDGE) {
1009         theedges.Append(Gvf);
1010       }
1011       if (Gvl.ShapeType() == TopAbs_EDGE) {
1012         theedges.Append(Gvl);
1013       }
1014       MakeFace(GenF,theedges);
1015       MapW(edg) = GenF;
1016     }
1017
1018   }
1019
1020
1021   TopTools_MapOfShape mapedgadded;
1022   TopTools_ListOfShape thefaces;
1023
1024   for (itl.Initialize(myMap(F)); itl.More(); itl.Next()) {
1025     const TopoDS_Face& fac = TopoDS::Face(itl.Value());
1026     theMap.Clear();
1027     TopoDS_Face DrftFace; // elle est FORWARD
1028     Standard_Boolean IsLeft;
1029     if (Contains(theLeft,fac)) {
1030       B.MakeFace(DrftFace,NewSg,BRep_Tool::Tolerance(fac));
1031       IsLeft = Standard_True;
1032     }
1033     else {
1034       B.MakeFace(DrftFace,NewSd,BRep_Tool::Tolerance(fac));
1035       IsLeft = Standard_False;
1036     }
1037     
1038     TopExp_Explorer exp3;
1039     for (exp3.Init(fac.Oriented(TopAbs_FORWARD),TopAbs_WIRE);
1040          exp3.More(); exp3.Next()) {
1041       const TopoDS_Shape& wir = exp3.Current();
1042       TopoDS_Wire NewWireOnF;
1043       B.MakeWire(NewWireOnF);
1044       for (exp.Init(wir.Oriented(TopAbs_FORWARD),TopAbs_EDGE);
1045            exp.More(); exp.Next()) {
1046         const TopoDS_Edge& edg = TopoDS::Edge(exp.Current());
1047         if (!theMap.Add(edg)) { // precaution sans doute inutile...
1048           continue;
1049         }
1050         if (MapW.IsBound(edg)) { // edge du wire d`origine
1051           TopTools_ListIteratorOfListOfShape itld;
1052           TopAbs_Orientation ored = edg.Orientation();
1053           if (IsLeft) {
1054             itld.Initialize(MapSg(edg));
1055           }
1056           else {
1057             itld.Initialize(MapSd(edg));
1058           }
1059           for (; itld.More(); itld.Next()) {
1060             if (itld.Value().Orientation() == TopAbs_REVERSED) {
1061               ored = TopAbs::Reverse(ored);
1062             }
1063             TopoDS_Shape aLocalEdge = itld.Value().Oriented(ored);
1064             B.Add(NewWireOnF,TopoDS::Edge(aLocalEdge));
1065 //          B.Add(NewWireOnF,TopoDS::Edge(itld.Value().Oriented(ored)));
1066           }
1067         }
1068         else {
1069           Handle(Geom_Surface) NewS;
1070           if (IsLeft) {
1071             NewS = NewSg;
1072           }
1073           else {
1074             NewS = NewSd;
1075           }
1076           Standard_Integer indedg = MapEV.FindIndex(edg);
1077           const TopoDS_Face& GenF = TopoDS::Face(MapEV(indedg));
1078           TopoDS_Vertex Vf,Vl;
1079           TopoDS_Shape aLocalEdge = edg.Oriented(TopAbs_FORWARD);
1080           TopExp::Vertices(TopoDS::Edge(aLocalEdge),Vf,Vl);
1081 //        TopExp::Vertices(TopoDS::Edge(edg.Oriented(TopAbs_FORWARD)),Vf,Vl);
1082           TopoDS_Shape Gvf = MapEV.FindFromKey(Vf); 
1083           TopoDS_Shape Gvl = MapEV.FindFromKey(Vl); 
1084           if (Gvf.ShapeType() == TopAbs_VERTEX &&
1085               Gvl.ShapeType() == TopAbs_VERTEX) {
1086             if (!Gvf.IsSame(Vf) || !Gvl.IsSame(Vl)) {
1087               TopoDS_Edge NewEdg = NewEdge(edg,GenF,NewS,
1088                                            TopoDS::Vertex(Gvf),
1089                                            TopoDS::Vertex(Gvl));
1090               if (NewEdg.IsNull()) {
1091                 return;
1092               }
1093
1094               MapSonS.Bind(edg,NewEdg);
1095
1096               if (NewEdg.Orientation() == TopAbs_REVERSED) {
1097                 NewEdg.Orientation(TopAbs::Reverse(edg.Orientation()));
1098               }
1099               else {
1100                 NewEdg.Orientation(edg.Orientation());
1101               }         
1102               B.Add(NewWireOnF,NewEdg);
1103             }
1104             else { // Frozen???
1105               B.Add(NewWireOnF,edg);
1106             }
1107           }
1108           else {
1109             TopoDS_Vertex Vff,Vll;
1110             if (Gvf.ShapeType() == TopAbs_VERTEX) {
1111               Vff = TopoDS::Vertex(Gvf);
1112             }
1113             else {
1114               if (IsLeft) {
1115                 Vff = TopoDS::Vertex(MapSg(Vf).First());
1116               }
1117               else {
1118                 Vff = TopoDS::Vertex(MapSd(Vf).First());
1119               }
1120             }
1121             if (Gvl.ShapeType() == TopAbs_VERTEX) {
1122               Vll = TopoDS::Vertex(Gvl);
1123             }
1124             else {
1125               if (IsLeft) {
1126                 Vll = TopoDS::Vertex(MapSg(Vl).First());
1127               }
1128               else {
1129                 Vll = TopoDS::Vertex(MapSd(Vl).First());
1130               }
1131             }
1132             
1133             TopoDS_Edge NewEdg = NewEdge(edg,GenF,NewS,Vff,Vll);
1134             if (NewEdg.IsNull()) {
1135               return;
1136             }
1137
1138             if (!MapW.IsBound(Vf) && !MapW.IsBound(Vl)) {
1139               MapSonS.Bind(edg,NewEdg);
1140             }
1141 //          else if (MapW.IsBound(Vf) && MapW.IsBound(Vl)) {
1142               
1143
1144 //          }
1145             else {
1146               if (MapW.IsBound(Vf)) {
1147                 if (Gvf.ShapeType() != TopAbs_EDGE || 
1148                     mapedgadded.Contains(Gvf)) {
1149                   MapSonS.Bind(edg,NewEdg);
1150                 }
1151                 else {
1152                   TopoDS_Wire NewWir;
1153                   B.MakeWire(NewWir);
1154                   B.Add(NewWir,NewEdg);
1155
1156                   TopoDS_Vertex Vf2,Vl2;
1157                   TopExp::Vertices(TopoDS::Edge(Gvf),Vf2,Vl2);
1158
1159                   //TopAbs_Orientation ornw = NewEdg.Orientation();
1160
1161                   // ici bug orientation : voir tspdrft6
1162
1163 //                if ((ornw == TopAbs_FORWARD && Vl2.IsSame(Vff)) ||
1164 //                    (ornw == TopAbs_REVERSED && Vl2.IsSame(Vll))) {
1165                   if (Vl2.IsSame(Vff)) {
1166                     B.Add(NewWir,Gvf.Oriented(TopAbs_FORWARD));
1167                   }
1168                   else {
1169                     B.Add(NewWir,Gvf.Oriented(TopAbs_REVERSED));
1170                   }
1171                   mapedgadded.Add(Gvf);
1172                   MapSonS.Bind(edg,NewWir); // NewWire est FORWARD
1173                 }
1174               }
1175               else {
1176                 if (Gvl.ShapeType() != TopAbs_EDGE ||
1177                     mapedgadded.Contains(Gvl)) {
1178                   MapSonS.Bind(edg,NewEdg);
1179                 }
1180                 else {
1181                   TopoDS_Wire NewWir;
1182                   B.MakeWire(NewWir);
1183                   B.Add(NewWir,NewEdg);
1184
1185                   TopoDS_Vertex Vf2,Vl2;
1186                   TopExp::Vertices(TopoDS::Edge(Gvl),Vf2,Vl2);
1187
1188                   //TopAbs_Orientation ornw = NewEdg.Orientation();
1189
1190                   // ici bug orientation : voir tspdrft6
1191
1192 //                if ((ornw == TopAbs_FORWARD && Vl2.IsSame(Vff)) ||
1193 //                    (ornw == TopAbs_REVERSED && Vl2.IsSame(Vll))) {
1194                   if (Vf2.IsSame(Vll)) {
1195                     B.Add(NewWir,Gvl.Oriented(TopAbs_FORWARD));
1196                   }
1197                   else {
1198                     B.Add(NewWir,Gvl.Oriented(TopAbs_REVERSED));
1199                   }
1200                   mapedgadded.Add(Gvl);
1201                   MapSonS.Bind(edg,NewWir); // NewWire est FORWARD
1202                 }
1203               }
1204             }
1205             if (NewEdg.Orientation() == TopAbs_REVERSED) {
1206               NewEdg.Orientation(TopAbs::Reverse(edg.Orientation()));
1207             }
1208             else {
1209               NewEdg.Orientation(edg.Orientation());
1210             }           
1211             B.Add(NewWireOnF,NewEdg);
1212           }
1213         }
1214       }
1215       B.Add(DrftFace,NewWireOnF.Oriented(wir.Orientation()));
1216     }
1217     thefaces.Append(DrftFace);
1218   }
1219
1220   BRepTools_Substitution theSubs;
1221   TopTools_DataMapIteratorOfDataMapOfShapeShape itdmss;
1222   for (itdmss.Initialize(MapSonS);
1223        itdmss.More(); itdmss.Next()) {
1224     TopTools_ListOfShape lsubs;
1225     for (exp.Init(itdmss.Value(),TopAbs_EDGE); exp.More(); exp.Next()) {
1226       lsubs.Append(exp.Current());
1227     }
1228     theSubs.Substitute(itdmss.Key(),lsubs);
1229   }
1230
1231   // on reconstruit les faces
1232   for (exp.Init(Res,TopAbs_FACE); exp.More(); exp.Next()) {
1233     if (Contains(myMap(F),exp.Current())) {
1234       continue;
1235     }
1236     theSubs.Build(exp.Current());
1237   }
1238
1239   // Stockage des descendants
1240 //  for (TopTools_DataMapIteratorOfDataMapOfShapeListOfShape itdmsls(myMap);
1241   TopTools_DataMapIteratorOfDataMapOfShapeListOfShape itdmsls(myMap) ;
1242   for ( ; itdmsls.More(); itdmsls.Next()) {
1243     if (itdmsls.Key().ShapeType() == TopAbs_EDGE) {
1244       TopTools_ListOfShape thedesc;
1245       theMap.Clear();
1246       for (itl.Initialize(itdmsls.Value());itl.More(); itl.Next()) {
1247         if (theMap.Add(MapW(itl.Value()))) {
1248           thedesc.Append(MapW(itl.Value()));
1249         }
1250       }
1251       myMap(itdmsls.Key()) = thedesc;
1252     }
1253     else if (itdmsls.Key().IsSame(F)) {
1254       myMap(F).Clear();
1255       for (itl.Initialize(thefaces); itl.More(); itl.Next()) {
1256         myMap(F).Append(itl.Value());
1257       }
1258     }
1259     else {
1260       TopTools_ListOfShape thedesc;
1261       theMap.Clear();
1262       for (itl.Initialize(itdmsls.Value());itl.More(); itl.Next()) {
1263         if (theSubs.IsCopied(itl.Value())) {
1264           if (theSubs.Copy(itl.Value()).Extent() != 1) {
1265 #ifdef DEB
1266             cout << "Invalid number of descendant" << endl;
1267 #endif
1268             return;
1269           }
1270           else {
1271             if (theMap.Add(theSubs.Copy(itl.Value()).First())) {
1272               thedesc.Append(theSubs.Copy(itl.Value()).First());
1273             }
1274           }
1275         }
1276         else if (theMap.Add(itl.Value())) {
1277           thedesc.Append(itl.Value());
1278         }
1279       }
1280       myMap(itdmsls.Key()) = thedesc;
1281     }
1282   }
1283
1284   theMap.Clear();
1285   thefaces.Clear();
1286   for (itdmsls.Initialize(myMap);itdmsls.More(); itdmsls.Next()) {
1287     for (itl.Initialize(itdmsls.Value());itl.More(); itl.Next()) {
1288       if (itl.Value().ShapeType() == TopAbs_FACE &&
1289           theMap.Add(itl.Value())) {
1290         thefaces.Append(itl.Value());
1291       }
1292     }
1293   }
1294   LocOpe_BuildShape BS(thefaces);
1295   myResult = BS.Shape();
1296 }
1297
1298
1299 //=======================================================================
1300 //function : Shape
1301 //purpose  : 
1302 //=======================================================================
1303
1304 const TopoDS_Shape& LocOpe_SplitDrafts::Shape () const
1305 {
1306   if (myResult.IsNull()) {
1307     StdFail_NotDone::Raise();
1308   }
1309   return myResult;
1310 }
1311
1312 //=======================================================================
1313 //function : ShapesFromShape
1314 //purpose  : 
1315 //=======================================================================
1316
1317 const TopTools_ListOfShape& LocOpe_SplitDrafts::ShapesFromShape 
1318    (const TopoDS_Shape& S) const
1319 {
1320   if (myResult.IsNull()) {
1321     StdFail_NotDone::Raise();
1322   }
1323   return myMap(S);
1324 }
1325
1326
1327 //=======================================================================
1328 //function : NewPlane
1329 //purpose  : 
1330 //=======================================================================
1331
1332 static Standard_Boolean NewPlane(const TopoDS_Face& F,
1333                                  const gp_Dir& Extr,
1334                                  const gp_Pln& Neutr,
1335                                  const Standard_Real Ang,
1336                                  gp_Pln& Newpl,
1337                                  gp_Ax1& NormalF,
1338                                  const Standard_Boolean Modify)
1339 {
1340
1341
1342   // Determination du nouveau plan incline
1343   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
1344   if (S->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
1345     S = Handle(Geom_RectangularTrimmedSurface)::DownCast(S)->BasisSurface();
1346   }
1347   Handle(Geom_Plane) P = Handle(Geom_Plane)::DownCast(S);
1348   if (P.IsNull()) {
1349     return Standard_False;
1350   }
1351
1352   gp_Pln Plorig = P->Pln();
1353   if (!Modify) {
1354     Newpl = Plorig;
1355     NormalF = Newpl.Axis();
1356     if (Newpl.Direct() && F.Orientation() == TopAbs_REVERSED ||
1357         !Newpl.Direct() && F.Orientation() == TopAbs_FORWARD) {
1358       NormalF.Reverse();
1359     }
1360     return Standard_True;
1361   }
1362
1363   gp_Ax1 Axe;
1364   Standard_Real Theta;
1365
1366   IntAna_QuadQuadGeo i2pl(Plorig,Neutr,
1367                           Precision::Angular(),Precision::Confusion());
1368   
1369   if (i2pl.IsDone() && i2pl.TypeInter() == IntAna_Line) {
1370     gp_Lin LinInters = i2pl.Line(1);
1371     gp_Dir nx = LinInters.Direction();
1372     NormalF = Plorig.Axis();
1373     gp_Dir ny = NormalF.Direction().Crossed(nx);
1374     Standard_Real a = Extr.Dot(nx);
1375     if (Abs(a) <=1-Precision::Angular()) { 
1376       Standard_Real b = Extr.Dot(ny);
1377       Standard_Real c = Extr.Dot(NormalF.Direction());
1378       Standard_Boolean direct(Plorig.Direct());
1379       TopAbs_Orientation Oris = F.Orientation();
1380       if ((direct && Oris == TopAbs_REVERSED) ||
1381           (!direct && Oris == TopAbs_FORWARD)) {
1382         b = -b;
1383         c = -c;
1384         NormalF.Reverse();
1385       }
1386       Standard_Real denom = Sqrt(1-a*a);
1387       Standard_Real Sina = Sin(Ang);
1388       if (denom>Abs(Sina)) {
1389         Standard_Real phi = ATan2(b/denom,c/denom);
1390         Standard_Real theta0 = ACos(Sina/denom); 
1391         Theta = theta0 - phi;
1392         if (Cos(Theta) <0.) {
1393           Theta = -theta0 -phi;
1394         }
1395         Axe = LinInters.Position();
1396         Newpl = Plorig.Rotated(Axe,Theta);
1397         return Standard_True;
1398       }
1399     }
1400   }
1401   cout << "fin newplane return standard_false" << endl;
1402   return Standard_False;
1403 }
1404
1405
1406 //=======================================================================
1407 //function : MakeFace
1408 //purpose  : 
1409 //=======================================================================
1410
1411 static void MakeFace(TopoDS_Face& F,
1412                      TopTools_ListOfShape& ledg)
1413 {
1414
1415   // ledg est une liste d'edge
1416
1417   BRep_Builder B;
1418
1419   // Verification de l`existence des p-curves. Celles qui manquent
1420   // correspondent necessairement a des isos (et meme des iso u).
1421
1422   Standard_Real f,l;
1423 //  for (TopTools_ListIteratorOfListOfShape itl(ledg); 
1424   TopTools_ListIteratorOfListOfShape itl(ledg) ;
1425   for ( ; itl.More(); itl.Next()) {
1426     TopoDS_Edge& edg = TopoDS::Edge(itl.Value());
1427     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(edg,F,f,l);
1428     if (C2d.IsNull()) {
1429       BRep_Tool::Range(edg,f,l);
1430       TopoDS_Vertex V1,V2;
1431       TopExp::Vertices(edg,V1,V2);
1432         TopTools_ListIteratorOfListOfShape itl2;
1433       for (itl2.Initialize(ledg); 
1434            itl2.More(); itl2.Next()) {
1435         const TopoDS_Edge& edg2 = TopoDS::Edge(itl2.Value());
1436         if (edg2.IsSame(edg)) {
1437           continue;
1438         }
1439         TopoDS_Vertex Vp1,Vp2;
1440         TopExp::Vertices(edg2,Vp1,Vp2);
1441         if (Vp1.IsSame(V1) || Vp2.IsSame(V1) ||
1442             Vp1.IsSame(V2) || Vp2.IsSame(V2)) {
1443           Standard_Real f2,l2;
1444           Handle(Geom2d_Curve) C22d = BRep_Tool::CurveOnSurface(edg2,F,f2,l2);
1445           if (!C22d.IsNull()) {
1446             gp_Pnt2d pt2d;
1447             if (Vp1.IsSame(V1)) {
1448               pt2d = C22d->Value(f2);
1449               pt2d.SetY(pt2d.Y()-f);
1450             }
1451             else if (Vp2.IsSame(V1)) {
1452               pt2d = C22d->Value(l2);
1453               pt2d.SetY(pt2d.Y()-f);
1454             }
1455             else if (Vp1.IsSame(V2)) {
1456               pt2d = C22d->Value(f2);
1457               pt2d.SetY(pt2d.Y()-l);
1458             }
1459             else if (Vp2.IsSame(V2)) {
1460               pt2d = C22d->Value(l2);
1461               pt2d.SetY(pt2d.Y()-l);
1462             }
1463             C2d = new Geom2d_Line(pt2d,gp::DY2d());
1464             B.UpdateEdge(edg,C2d,F,BRep_Tool::Tolerance(edg));
1465             break;
1466           }
1467         }
1468       }
1469       if (C2d.IsNull()) {
1470         cout << "Ca merde violemment" << endl;
1471       }
1472     }
1473   }
1474
1475   TopTools_ListOfShape lwires;
1476   Standard_Boolean alldone = ledg.IsEmpty();
1477   while (!alldone) {
1478     TopoDS_Wire Wnew;
1479     B.MakeWire(Wnew);
1480     TopoDS_Shape aLocalShape = ledg.First();
1481     const TopoDS_Edge& edg = TopoDS::Edge(aLocalShape);
1482 //    const TopoDS_Edge& edg = TopoDS::Edge(ledg.First());
1483     TopoDS_Vertex VFirst,VLast;
1484     if (edg.Orientation() == TopAbs_FORWARD) {
1485       TopExp::Vertices(edg,VFirst,VLast);
1486     }
1487     else {
1488       TopExp::Vertices(edg,VLast,VFirst);
1489     }
1490     B.Add(Wnew,edg);
1491     ledg.RemoveFirst();
1492     // on suppose VFirst et VLast non nuls
1493     Standard_Boolean wdone = (ledg.IsEmpty() || VFirst.IsSame(VLast));
1494     while (!wdone) {
1495       TopoDS_Vertex VF,VL;
1496 #ifdef DEB
1497       TopAbs_Orientation oredg;
1498 #else
1499       TopAbs_Orientation oredg = TopAbs_FORWARD;
1500 #endif
1501       for (itl.Initialize(ledg); itl.More(); itl.Next()) {
1502         const TopoDS_Edge& edg2 = TopoDS::Edge(itl.Value());
1503         TopoDS_Shape aLocalShape  = edg2.Oriented(TopAbs_FORWARD);
1504         TopExp::Vertices(TopoDS::Edge(aLocalShape),VF,VL);
1505 //      TopExp::Vertices(TopoDS::Edge(edg2.Oriented(TopAbs_FORWARD)),VF,VL);
1506         if (VF.IsSame(VLast)) {
1507           VLast = VL;
1508           oredg = TopAbs_FORWARD;
1509           break;
1510         }
1511         else if (VL.IsSame(VFirst)) {
1512           VFirst = VF;
1513           oredg = TopAbs_FORWARD;
1514           break;
1515         }
1516         else if (VF.IsSame(VFirst)) {
1517           VFirst = VL;
1518           oredg = TopAbs_REVERSED;
1519           break;
1520         }
1521         else if (VL.IsSame(VLast)) {
1522           VLast = VF;
1523           oredg = TopAbs_REVERSED;
1524           break;
1525         }
1526
1527       }
1528       if (!itl.More()) {
1529         wdone = Standard_True;
1530       }
1531       else {
1532         TopoDS_Shape aLocalShape = itl.Value().Oriented(oredg);
1533         B.Add(Wnew,TopoDS::Edge(aLocalShape));
1534 //      B.Add(Wnew,TopoDS::Edge(itl.Value().Oriented(oredg)));
1535         ledg.Remove(itl);
1536         wdone = (ledg.IsEmpty() || VFirst.IsSame(VLast));
1537       }
1538     }
1539     lwires.Append(Wnew);
1540     alldone = ledg.IsEmpty();
1541   }
1542
1543
1544
1545
1546   F.Orientation(TopAbs_FORWARD);
1547   for (itl.Initialize(lwires); itl.More(); itl.Next()) {
1548     TopoDS_Shape aLocalShape = F.EmptyCopied();
1549     TopoDS_Face NewFace = TopoDS::Face(aLocalShape);
1550 //    TopoDS_Face NewFace = TopoDS::Face(F.EmptyCopied());
1551     B.Add(NewFace,itl.Value());
1552     GProp_GProps GP;
1553     BRepGProp::SurfaceProperties(NewFace,GP);
1554     if (GP.Mass() < 0) {
1555       itl.Value().Reverse();
1556     }
1557   }
1558   if (lwires.Extent() == 1) {
1559     B.Add(F,lwires.First());
1560   }
1561   else {
1562     cout << "Not yet implemented : nbwire >= 2" << endl;
1563   }
1564
1565 }
1566
1567
1568 //=======================================================================
1569 //function : Contains
1570 //purpose  : 
1571 //=======================================================================
1572
1573 static Standard_Boolean Contains(const TopTools_ListOfShape& ll,
1574                                  const TopoDS_Shape& s)
1575 {
1576   TopTools_ListIteratorOfListOfShape itl;
1577   for (itl.Initialize(ll); itl.More(); itl.Next()) {
1578     if (itl.Value().IsSame(s)) {
1579       return Standard_True;
1580     }
1581   }
1582   return Standard_False;
1583 }
1584
1585
1586
1587 //=======================================================================
1588 //function : Contains
1589 //purpose  : 
1590 //=======================================================================
1591
1592 static TopoDS_Edge  NewEdge(const TopoDS_Edge& edg,
1593                             const TopoDS_Face& F,
1594                             const Handle(Geom_Surface)& NewS,
1595                             const TopoDS_Vertex& V1,
1596                             const TopoDS_Vertex& V2)
1597 {
1598   TopoDS_Edge NewEdg;
1599   Handle(Geom_Surface) S1 = BRep_Tool::Surface(F);
1600   Standard_Boolean AppS1 = Standard_False;
1601   if (S1->DynamicType() != STANDARD_TYPE(Geom_Plane)) {
1602     AppS1 = Standard_True;
1603   }
1604
1605
1606   GeomInt_IntSS i2s(S1,NewS,Precision::Confusion(),Standard_True,AppS1);
1607   if (!i2s.IsDone() || i2s.NbLines() <= 0) {
1608     return NewEdg;
1609   }
1610
1611   BRep_Builder B;
1612 //  Standard_Real pmin, Dist, DistMin;
1613   Standard_Real Dist2, Dist2Min;
1614   Standard_Real prmf=0,prml=0;
1615   GeomAdaptor_Curve TheCurve;
1616         
1617   Standard_Integer i,imin,k;
1618   gp_Pnt pvf = BRep_Tool::Pnt(V1);
1619   gp_Pnt pvl = BRep_Tool::Pnt(V2);
1620   imin = 0;
1621   for (i=1; i<= i2s.NbLines(); i++) {
1622     TheCurve.Load(i2s.Line(i));
1623     Extrema_ExtPC myExtPC(pvf,TheCurve);
1624     
1625     if (myExtPC.IsDone()) {
1626       gp_Pnt p1b,p2b;
1627       Standard_Real thepmin = TheCurve.FirstParameter();
1628       myExtPC.TrimmedSquareDistances(Dist2Min,Dist2,p1b,p2b);
1629       if (Dist2 < Dist2Min && !TheCurve.IsPeriodic()) {
1630         Dist2Min = Dist2;
1631         thepmin = TheCurve.LastParameter();
1632       }
1633       for (k=1; k<=myExtPC.NbExt(); k++) {
1634         Dist2 = myExtPC.SquareDistance(k);
1635         if (Dist2 < Dist2Min) {
1636           Dist2Min = Dist2;
1637           thepmin = myExtPC.Point(k).Parameter();
1638         }
1639       }
1640       
1641       if (Dist2Min  <= Precision::SquareConfusion()) {
1642         prmf = thepmin;
1643         myExtPC.Perform(pvl);
1644         if (myExtPC.IsDone()) {
1645           thepmin = TheCurve.LastParameter();
1646           myExtPC.TrimmedSquareDistances(Dist2,Dist2Min,p1b,p2b);
1647           if (Dist2 < Dist2Min && !TheCurve.IsClosed()) {
1648             Dist2Min = Dist2;
1649             thepmin = TheCurve.FirstParameter();
1650           }
1651           for (k=1; k<=myExtPC.NbExt(); k++) {
1652             Dist2 = myExtPC.SquareDistance(k);
1653             if (Dist2 < Dist2Min) {
1654               Dist2Min = Dist2;
1655               thepmin = myExtPC.Point(k).Parameter();
1656             }
1657           }
1658           
1659       if (Dist2Min  <= Precision::SquareConfusion()) {
1660             prml = thepmin;
1661             break;
1662           }
1663         }
1664       }
1665     }
1666   }
1667
1668   if (i <= i2s.NbLines()) {
1669     Standard_Boolean rev = Standard_False;
1670     TopoDS_Vertex Vf = V1;
1671     TopoDS_Vertex Vl = V2;
1672     Handle(Geom_Curve) Cimg = i2s.Line(i);
1673     Handle(Geom2d_Curve) Cimg2d;
1674     if (AppS1) {
1675       Cimg2d = i2s.LineOnS1(i);
1676     }
1677
1678     if (Cimg->IsPeriodic()) {
1679
1680       Standard_Real period = Cimg->Period();
1681       Standard_Real imf = Cimg->FirstParameter();
1682       Standard_Real iml = Cimg->LastParameter();
1683
1684       Standard_Real f,l;
1685       BRep_Tool::Range(edg,f,l);
1686       Standard_Real delt = l-f;
1687       Standard_Real delt1 = Abs(prml-prmf);
1688       Standard_Real delt2 = Abs(period-delt1);
1689       
1690       if (delt1 == 0 || delt2 == 0) {
1691 //      prmf = 0;
1692 //      prml = period;
1693         prmf = imf;
1694         prml = iml;
1695       }
1696       else {
1697         if (Abs(delt1-delt) > Abs(delt2-delt)) {
1698           // le bon ecart est delt2...
1699           if (prml > prmf) {
1700             if (prml < iml) {
1701               prmf += period;
1702             }
1703             else {
1704               prml -= period;
1705             }
1706           }
1707           else {
1708             if (prmf < iml) {
1709               prml += period;
1710             }
1711             else {
1712               prmf -= period;
1713             }
1714           }
1715         }
1716         else if (Abs(delt1-delt) < Abs(delt2-delt)) {
1717           if (prmf >= iml && prml >= iml) {
1718             prmf -= period;
1719             prml -= period;
1720           }
1721           else if (prmf <= imf && prml <= imf) {
1722             prmf += period;
1723             prml += period;
1724           }
1725         }
1726         else { // egalite; on priveligie l'ordre f,l
1727           if (prmf > prml) {
1728             prmf -= period;
1729           }
1730           if (prmf >= iml && prml >= iml) {
1731             prmf -= period;
1732             prml -= period;
1733           }
1734           else if (prmf <= imf && prml <= imf) {
1735             prmf += period;
1736             prml += period;
1737           }
1738         }
1739       }
1740 #ifdef DEB
1741       Standard_Real ptol = Precision::PConfusion();
1742       if (prmf < imf - ptol || prmf > iml + ptol ||
1743           prml < imf - ptol || prml > iml + ptol) {
1744         cout << "Ca ne va pas aller" << endl;
1745       }
1746 #endif
1747
1748
1749     }
1750
1751     if (S1->IsUPeriodic()) {
1752
1753       Standard_Real speriod = S1->UPeriod();
1754 //      Standard_Real f,l;
1755       gp_Pnt2d pf,pl;
1756       pf = Cimg2d->Value(prmf);
1757       pl = Cimg2d->Value(prml);
1758
1759       Standard_Real Uf = pf.X();
1760       Standard_Real Ul = pl.X();
1761       Standard_Real ptra = 0.0;
1762
1763       Standard_Real Ustart = Min(Uf,Ul);
1764       while (Ustart < -Precision::PConfusion()) {
1765         Ustart += speriod;
1766         ptra += speriod;
1767       }
1768       while (Ustart > speriod - Precision::PConfusion()) {
1769         Ustart -= speriod;
1770         ptra -= speriod;
1771       }
1772       if (ptra != 0) {
1773         Cimg2d->Translate(gp_Vec2d(ptra,0.));
1774       }
1775
1776
1777     }
1778     if (prmf < prml) {
1779       Vf.Orientation(TopAbs_FORWARD);
1780       Vl.Orientation(TopAbs_REVERSED);
1781     }
1782     else {
1783       Vf.Orientation(TopAbs_REVERSED);
1784       Vl.Orientation(TopAbs_FORWARD);
1785       rev = Standard_True;;
1786     }
1787
1788     B.MakeEdge(NewEdg,Cimg,Precision::Confusion());
1789
1790     B.Add(NewEdg,Vf);
1791     B.Add(NewEdg,Vl);
1792     B.UpdateVertex(Vf,prmf,NewEdg,Precision::Confusion());
1793     B.UpdateVertex(Vl,prml,NewEdg,Precision::Confusion());
1794     if (AppS1) {
1795       B.UpdateEdge(NewEdg,Cimg2d,F,Precision::Confusion());
1796     }
1797
1798
1799     if (rev) {
1800       NewEdg.Orientation(TopAbs_REVERSED);
1801     }
1802   }
1803   return NewEdg;
1804 }