src/PARTITION/Partition_Loop3d.cxx

   1 // Copyright (C) 2005  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   2 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   3 //
   4 // This library is free software; you can redistribute it and/or
   5 // modify it under the terms of the GNU Lesser General Public
   6 // License as published by the Free Software Foundation; either
   7 // version 2.1 of the License.
   8 //
   9 // This library is distributed in the hope that it will be useful
  10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 // Lesser General Public License for more details.
  13 //
  14 // You should have received a copy of the GNU Lesser General Public
  15 // License along with this library; if not, write to the Free Software
  16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 //
  20 //  GEOM PARTITION : partition algorithm
  21 //
  22 //  File   : Partition_Loop3d.cxx
  23 //  Module : GEOM
  24
  25 #include "Partition_Loop3d.ixx"
  26
  27 #include <TopExp_Explorer.hxx>
  28 #include <TopExp.hxx>
  29 #include <BRep_Builder.hxx>
  30 #include <TopTools_MapOfShape.hxx>
  31 #include <TopTools_ListIteratorOfListOfShape.hxx>
  32 #include <TopoDS_Shell.hxx>
  33 #include <TopoDS_Iterator.hxx>
  34 #include <TopoDS.hxx>
  35 #include <TopTools_MapIteratorOfMapOfShape.hxx>
  36 #include <gp_Vec.hxx>
  37 #include <gp_Pnt.hxx>
  38 #include <Geom2d_Curve.hxx>
  39 #include <BRep_Tool.hxx>
  40 #include <Geom_Surface.hxx>
  41 #include <gp_Pnt2d.hxx>
  42 #include <gp_Vec2d.hxx>
  43 #include <gp_Dir2d.hxx>
  44 #include <Geom_Curve.hxx>
  45
  46 using namespace std;
  47
  48 //=======================================================================
  49 //function : Partition_Loop3d
  50 //purpose  :
  51 //=======================================================================
  52
  53 Partition_Loop3d::Partition_Loop3d()
  54 {
  55 }
  56
  57 //=======================================================================
  58 //function : AddConstFaces
  59 //purpose  : Add faces of <S> as unique faces in the result.
  60 //=======================================================================
  61
  62 void Partition_Loop3d::AddConstFaces(const TopoDS_Shape& S)
  63 {
  64   TopExp_Explorer FaceExp(S, TopAbs_FACE);
  65   for (; FaceExp.More(); FaceExp.Next())
  66     myFaces.Append( FaceExp.Current() );
  67
  68   TopExp::MapShapesAndAncestors(S, TopAbs_EDGE, TopAbs_FACE, myEFMap);
  69 }
  70
  71 //=======================================================================
  72 //function : AddSectionFaces
  73 //purpose  : Add faces of <S> as double faces in the result.
  74 //=======================================================================
  75
  76 void Partition_Loop3d::AddSectionFaces(const TopoDS_Shape& S)
  77 {
  78   AddConstFaces( S );
  79   AddConstFaces( S.Reversed() );
  80 }
  81
  82 //=======================================================================
  83 //function : MakeShells
  84 //purpose  : Make and return shells.
  85 //           <AvoidFacesMap> can contain faces that must not be
  86 //           added to result shells.
  87 //=======================================================================
  88
  89 const TopTools_ListOfShape&
  90   Partition_Loop3d::MakeShells (const TopTools_MapOfOrientedShape& AvoidFacesMap)
  91 {
  92   myNewShells.Clear();
  93
  94   BRep_Builder Builder;
  95   TopTools_MapOfShape CheckedEdgesMap;
  96   TopTools_MapOfOrientedShape AddedFacesMap;
  97
  98   TopTools_ListIteratorOfListOfShape itF (myFaces);
  99   for (; itF.More(); itF.Next())
 100   {
 101     const TopoDS_Shape& FF = itF.Value();
 102     if (AvoidFacesMap.Contains( FF ) ||
 103         ! AddedFacesMap.Add( FF ) )
 104       continue;
 105
 106     // make a new shell
 107     TopoDS_Shell Shell;
 108     Builder.MakeShell(Shell);
 109     Builder.Add(Shell,FF);
 110
 111     // clear the maps from shapes added to previous Shell
 112     TopTools_MapIteratorOfMapOfShape itEM (CheckedEdgesMap);
 113     for (; itEM.More(); itEM.Next()) {
 114       TopTools_ListOfShape& FL = myEFMap.ChangeFromKey( itEM.Key());
 115       TopTools_ListIteratorOfListOfShape it (FL);
 116       while ( it.More()) {
 117         if (AddedFacesMap.Contains( it.Value()))
 118           FL.Remove( it );
 119         else
 120           it.Next();
 121       }
 122     }
 123     CheckedEdgesMap.Clear();
 124
 125
 126     // loop on faces added to Shell; add their neighbor faces to Shell and so on
 127     TopoDS_Iterator itAddedF (Shell);
 128     for (; itAddedF.More(); itAddedF.Next())
 129     {
 130       const TopoDS_Face& F = TopoDS::Face (itAddedF.Value());
 131
 132       // loop on edges of F; find a good neighbor face of F by E
 133       TopExp_Explorer EdgeExp(F, TopAbs_EDGE);
 134       for (; EdgeExp.More(); EdgeExp.Next())
 135       {
 136         const TopoDS_Edge& E = TopoDS::Edge( EdgeExp.Current());
 137         if (! CheckedEdgesMap.Add( E ))
 138           continue;
 139
 140         // candidate faces list
 141         const TopTools_ListOfShape& FL = myEFMap.ChangeFromKey(E);
 142         if (FL.IsEmpty())
 143           continue;
 144         // select one of neighbors
 145         TopoDS_Face SelF;
 146         if (FL.Extent() == 2) {
 147           if (! F.IsSame( FL.First() ))
 148             SelF = TopoDS::Face( FL.First() );
 149           else if (!F.IsSame( FL.Last() ))
 150             SelF = TopoDS::Face( FL.Last() );
 151         }
 152         else {
 153           // check if a face already added to Shell shares E
 154           TopTools_ListIteratorOfListOfShape it (FL);
 155           Standard_Boolean found = Standard_False;
 156           for (; !found && it.More(); it.Next())
 157             if (F != it.Value())
 158               found = AddedFacesMap.Contains( it.Value() );
 159           if (found)
 160             continue;
 161           // select basing on geometrical check
 162           Standard_Boolean GoodOri, inside;
 163           Standard_Real dot, MaxDot = -100;
 164           TopTools_ListOfShape TangFL; // tangent faces
 165           for ( it.Initialize( FL ) ; it.More(); it.Next()) {
 166             const TopoDS_Face& NeighborF = TopoDS::Face( it.Value());
 167             if (NeighborF.IsSame( F ))
 168               continue;
 169             inside = Partition_Loop3d::IsInside( E, F, NeighborF, 1, dot, GoodOri);
 170             if (!GoodOri)
 171               continue;
 172             if (!inside)
 173               dot = -dot - 3;
 174             if (dot < MaxDot)
 175               continue;
 176             if ( IsEqual( dot, MaxDot))
 177               TangFL.Append(SelF);
 178             else
 179               TangFL.Clear();
 180             MaxDot = dot;
 181             SelF = NeighborF;
 182           }
 183           if (!TangFL.IsEmpty()) {
 184             for (it.Initialize( TangFL ); it.More(); it.Next()) {
 185               const TopoDS_Face& NeighborF = TopoDS::Face( it.Value());
 186               if (Partition_Loop3d:: IsInside( E, SelF , NeighborF, 0, dot, GoodOri))
 187                 SelF = NeighborF;
 188             }
 189           }
 190         }
 191         if (!SelF.IsNull() &&
 192             AddedFacesMap.Add( SelF ) &&
 193             !AvoidFacesMap.Contains( SelF ))
 194           Builder.Add( Shell, SelF);
 195
 196       } // loop on edges of F
 197
 198     } // loop on the faces added to Shell
 199
 200     // Shell is complete
 201     myNewShells.Append( Shell );
 202
 203   } // loop on myFaces
 204
 205
 206   // prepare to the next call
 207   myFaces.Clear();
 208   myEFMap.Clear();
 209
 210   return myNewShells;
 211 }
 212
 213
 214
 215 //=======================================================================
 216 //function : Normal
 217 //purpose  :
 218 //=======================================================================
 219
 220 gp_Vec Partition_Loop3d::Normal(const TopoDS_Edge& E,
 221                                 const TopoDS_Face& F)
 222 {
 223   gp_Vec Norm, V1, V2;
 224   Standard_Real First, Last;
 225   gp_Pnt Ps;
 226
 227   Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface (E, F, First, Last);
 228   Handle(Geom_Surface) Sf = BRep_Tool::Surface(F);
 229
 230   gp_Pnt2d p = C2d->Value( 0.5*(First+Last) );
 231   Sf->D1(p.X(), p.Y(), Ps, V1, V2);
 232   Norm = V1.Crossed(V2);
 233
 234   if (F.Orientation() == TopAbs_REVERSED )
 235     Norm.Reverse();
 236
 237   return Norm;
 238 }
 239
 240 //=======================================================================
 241 //function : NextNormal
 242 //purpose  : find normal to F at point a little inside F near the middle of E
 243 //warning  : E must be properly oriented in F.
 244 //=======================================================================
 245
 246 static gp_Vec NextNormal(const TopoDS_Edge& E,
 247                          const TopoDS_Face& F)
 248 {
 249   Standard_Real First, Last;
 250
 251   Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface (E, F, First, Last);
 252   Handle(Geom_Surface) Sf = BRep_Tool::Surface(F);
 253
 254   gp_Pnt2d p;
 255   gp_Vec2d v;
 256   C2d->D1( 0.5*(First+Last), p, v);
 257   if (E.Orientation() != F.Orientation())
 258     v.Reverse();
 259   gp_Dir2d dir( -v.Y(), v.X() ); // dir inside F
 260
 261   Standard_Real duv = 1e-6; // this is not Ok and may give incorrect result if
 262   // resolutionUV of compared faces is very different. To have a good result,
 263   //it is necessary to get normal to faces at points equidistant from E in 3D
 264
 265   p.SetX( p.X() + dir.X()*duv );
 266   p.SetY( p.Y() + dir.Y()*duv );
 267
 268   gp_Pnt Ps;
 269   gp_Vec Norm, V1, V2, VV1, VV2;
 270   Sf->D1( p.X(), p.Y(), Ps, V1, V2);
 271   Norm = V1.Crossed(V2);
 272
 273   if (F.Orientation() == TopAbs_REVERSED )
 274     Norm.Reverse();
 275
 276   return Norm;
 277 }
 278
 279
 280 //=======================================================================
 281 //function : FindEinF
 282 //purpose  : find E in F
 283 //=======================================================================
 284
 285 static TopoDS_Edge FindEinF(const TopoDS_Edge& E,
 286                             const TopoDS_Face& F)
 287 {
 288   TopExp_Explorer expl (F, TopAbs_EDGE);
 289   for (; expl.More(); expl.Next())
 290     if( E.IsSame( expl.Current() ))
 291       return TopoDS::Edge(expl.Current());
 292   TopoDS_Edge nullE;
 293   return nullE;
 294 }
 295
 296 //=======================================================================
 297 //function : IsInside
 298 //purpose  : check if <F2> is inside <F1> by edge <E>.
 299 //           if <CountDot>, compute <Dot>: scalar production of
 300 //           normalized  vectors  pointing  inside  faces,  and
 301 //           check if faces are oriented well for sewing
 302 //=======================================================================
 303
 304 Standard_Boolean Partition_Loop3d::IsInside(const TopoDS_Edge& E,
 305                                             const TopoDS_Face& F1,
 306                                             const TopoDS_Face& F2,
 307                                             const Standard_Boolean CountDot,
 308                                             Standard_Real& Dot,
 309                                             Standard_Boolean& GoodOri)
 310 {
 311   Standard_Real f, l;
 312   gp_Pnt P;
 313   gp_Vec Vc1, Vc2, Vin1, Vin2, Nf1, Nf2;
 314   Handle(Geom_Curve) Curve = BRep_Tool::Curve(E,f,l);
 315   Curve->D1( 0.5*(f + l), P, Vc2);
 316   TopoDS_Edge E1, E2 = FindEinF (E, F2);
 317   if (E2.Orientation() == TopAbs_REVERSED ) Vc2.Reverse();
 318
 319   Nf1 = Normal(E,F1);
 320   Nf2 = Normal(E,F2);
 321
 322   Standard_Real sin =
 323     Nf1.CrossSquareMagnitude(Nf2) / Nf1.SquareMagnitude() / Nf2.SquareMagnitude();
 324   Standard_Boolean tangent = sin < 0.001;
 325
 326   Standard_Boolean inside = 0;
 327   if (tangent) {
 328     E1 = FindEinF (E, F1);
 329     gp_Vec NNf1 = NextNormal(E1,F1);
 330     gp_Vec NNf2 = NextNormal(E2,F2);
 331     Vin2 = NNf2.Crossed(Vc2);
 332     inside = Vin2 * NNf1 < 0;
 333   }
 334   else {
 335     Vin2 = Nf2.Crossed(Vc2);
 336     inside = Vin2 * Nf1 < 0;
 337   }
 338
 339   if (!CountDot) return inside;
 340
 341   if (tangent)
 342     Vin2 = Nf2.Crossed(Vc2);
 343   else
 344     E1 = FindEinF (E, F1);
 345
 346   Vc1 = Vc2;
 347   if (E1.Orientation() != E2.Orientation())
 348     Vc1.Reverse();
 349   Vin1 = Nf1.Crossed(Vc1);
 350
 351   if (tangent) {
 352     Standard_Real N1N2 = Nf1 * Nf2;
 353     GoodOri = (Vin2 * Vin1 < 0) ? N1N2 > 0 : N1N2 < 0;
 354   }
 355   else {
 356     Standard_Real V1N2 = Vin1 * Nf2;
 357     GoodOri = ( inside ? V1N2 <= 0 : V1N2 >= 0);
 358   }
 359
 360   Vin1.Normalize();
 361   Vin2.Normalize();
 362
 363   Dot = Vin2 * Vin1;
 364
 365   return inside;
 366 }
 367