src/PARTITION/Partition_Loop3d.cxx

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