src/ARCHIMEDE/Archimede_VolumeSection.cxx

   1 // Copyright (C) 2007-2012  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
  23 //  GEOM ARCHIMEDE : algorithm implementation
  24 //  File   : Archimede_VolumeSection.cxx
  25 //  Author : Nicolas REJNERI
  26 //  Module : GEOM
  27 //  $Header$
  28 //
  29 #include "Archimede_VolumeSection.hxx"
  30 #include "utilities.h"
  31
  32 #include <BRepMesh_IncrementalMesh.hxx>
  33 #include <TopExp_Explorer.hxx>
  34 #include <TopLoc_Location.hxx>
  35 #include <Poly_Triangulation.hxx>
  36 #include <Poly_Array1OfTriangle.hxx>
  37 #include <BRep_Tool.hxx>
  38 #include <TopoDS.hxx>
  39 #include <TopoDS_Face.hxx>
  40 #include <TopoDS_Shape.hxx>
  41 #include <math_Matrix.hxx>
  42 #include <gp_Trsf.hxx>
  43 #include <gp_Dir.hxx>
  44 #include <gp_Ax1.hxx>
  45 #include <gp_Pnt.hxx>
  46
  47 #include <GeomAPI_ProjectPointOnSurf.hxx>
  48 #include <Geom_RectangularTrimmedSurface.hxx>
  49
  50 //-------------------------------------------------------------------------------------------------------
  51 //----------------------------------- Methodes publiques -------------------------------------------------
  52 //-------------------------------------------------------------------------------------------------------
  53
  54 //  Maillage de la shape
  55 VolumeSection::VolumeSection(TopoDS_Shape S , Standard_Real Precision):myShape(S),Tolerance(Precision)
  56 {
  57   // Maillage de la shape myShape
  58   BRepMesh_IncrementalMesh(myShape,Tolerance);
  59 }
  60
  61 TopoDS_Shape VolumeSection::GetShape()
  62 {
  63   return myShape;
  64 }
  65
  66 void VolumeSection::SetPlane(Handle (Geom_Plane) P)
  67 {
  68   myPlane = P;
  69 }
  70
  71 void VolumeSection::CenterOfGravity()
  72 {
  73   Standard_Integer i;
  74   Standard_Integer nbNodes;
  75   TopExp_Explorer ex;
  76   TopLoc_Location L;
  77
  78   // Boucle sur les faces de la shape
  79
  80   Xmin = 1000000000;
  81   Ymin = 1000000000;
  82   Zmin = 1000000000;
  83   Xmax = -1000000000;
  84   Ymax = -1000000000;
  85   Zmax = -1000000000;
  86
  87   for (ex.Init(myShape, TopAbs_FACE); ex.More(); ex.Next())
  88     {
  89       TopoDS_Face F = TopoDS::Face(ex.Current());
  90       Handle(Poly_Triangulation) Tr = BRep_Tool::Triangulation(F, L);
  91       if(Tr.IsNull())
  92         MESSAGE("Error, null layer" )
  93       nbNodes = Tr->NbNodes();
  94       const TColgp_Array1OfPnt& Nodes = Tr->Nodes();
  95
  96       // Calcul des dimensions de la boite englobante du solide
  97
  98       for(i=1;i<=nbNodes;i++)
  99         {
 100           InitPoint = Nodes(i).Transformed(L.Transformation());
 101           if(InitPoint.X() < Xmin)
 102             Xmin = InitPoint.X();
 103           if(InitPoint.X() > Xmax)
 104             Xmax = InitPoint.X();
 105           if(InitPoint.Y() < Ymin)
 106             Ymin = InitPoint.Y();
 107           if(InitPoint.Y() > Ymax)
 108             Ymax = InitPoint.Y();
 109           if(InitPoint.Z() < Zmin)
 110             Zmin = InitPoint.Z();
 111           if(InitPoint.Z() > Zmax)
 112             Zmax = InitPoint.Z();
 113
 114         }
 115     }
 116
 117   // Creation du point d'initialisation, c'est \80 dire le centre de gravit\89
 118   //g\89om\89trique de la boite englobante
 119
 120   InitPoint.SetX(0.5 * (Xmin + Xmax));
 121   InitPoint.SetY(0.5 * (Ymin + Ymax));
 122   InitPoint.SetZ(0);
 123 }
 124
 125 Standard_Real VolumeSection::CalculateVolume(Standard_Real Elevation)
 126 {
 127   Standard_Integer i,noeud[3],flag[3];
 128   Standard_Integer nbNodes;
 129   TopExp_Explorer ex;
 130   TopLoc_Location L;
 131   Standard_Real z[3];
 132   Standard_Real Volume=0;
 133   Standard_Real Determinant=0;
 134   gp_Pnt P[3];
 135
 136   // Projection du point d'initialisation sur le plan de section
 137
 138   InitPoint.SetZ(Elevation);
 139
 140   for (ex.Init(myShape, TopAbs_FACE); ex.More(); ex.Next())
 141     {
 142       TopoDS_Face F = TopoDS::Face(ex.Current());
 143       Handle(Poly_Triangulation) Tr = BRep_Tool::Triangulation(F, L);
 144       if(Tr.IsNull())
 145         MESSAGE("Error, null layer" )
 146       const Poly_Array1OfTriangle& triangles = Tr->Triangles();
 147       Standard_Integer nbTriangles = Tr->NbTriangles();
 148       nbNodes = Tr->NbNodes();
 149       const TColgp_Array1OfPnt& Nodes = Tr->Nodes();
 150
 151       // Calcul des volumes de chaque triangle, de chaque face
 152       //en tenant compte des triangles coup\89s par le plan de section
 153
 154       for (i=1;i<=nbTriangles;i++)
 155         {
 156           Determinant=0;
 157           //Gardons la meme orientation des noeuds
 158           if (F.Orientation()  == TopAbs_REVERSED)
 159             triangles(i).Get(noeud[0], noeud[2], noeud[1]);
 160           else
 161             triangles(i).Get(noeud[0], noeud[1], noeud[2]);
 162
 163           P[0] = Nodes(noeud[0]).Transformed(L.Transformation());
 164           z[0] = P[0].Z();
 165           P[1] = Nodes(noeud[1]).Transformed(L.Transformation());
 166           z[1] = P[1].Z();
 167           P[2] = Nodes(noeud[2]).Transformed(L.Transformation());
 168           z[2] = P[2].Z();
 169
 170           // Determination des cas aux limites pour les triangles
 171           Standard_Integer i,compteur=0;
 172
 173           for (i=0;i<=2;i++)
 174             {
 175               flag[i]=Standard_False;
 176               if(z[i]>=Elevation)
 177                 {
 178                   flag[i]=Standard_True;
 179                   compteur++;
 180                 }
 181             }
 182
 183           switch(compteur)
 184             {
 185             case 0:
 186               Determinant = ElementaryVolume(P[0],P[1],P[2]);
 187               break;
 188
 189             case 1:
 190               for (i=0;i<=2;i++)
 191                 {
 192                   if (flag[i]==Standard_True)
 193                     {
 194                       gp_Pnt Result1 = Intersection(P[i],P[(i+1)%3],Elevation);
 195                       gp_Pnt Result2 = Intersection(P[i],P[(i+2)%3],Elevation);
 196                       Determinant = ElementaryVolume(Result1,P[(i+1)%3],P[(i+2)%3])
 197                         + ElementaryVolume(Result1,P[(i+2)%3],Result2);
 198                     }
 199                 }
 200               break;
 201
 202             case 2:
 203               for (i=0;i<=2;i++)
 204                 {
 205                   if (flag[i]==Standard_False)
 206                     {
 207                       gp_Pnt Result1 = Intersection(P[i],P[(i+1)%3],Elevation);
 208                       gp_Pnt Result2 = Intersection(P[i],P[(i+2)%3],Elevation);
 209                       Determinant = ElementaryVolume(P[i],Result1,Result2);
 210                     }
 211                 }
 212               break;
 213
 214             case 3:
 215               break;
 216             }
 217           Volume += Determinant;
 218         }
 219     }
 220
 221   return Volume;
 222 }
 223
 224 Standard_Real VolumeSection::Archimede(Standard_Real Constante , Standard_Real Epsilon)
 225 {
 226   // Resolution de l equation V(h) = Constante a l aide de l algorithme de dichotomie avec ponderation type
 227   // Lagrange
 228
 229   Standard_Real c,Binf,Bsup;
 230   Standard_Real tempBsupVolume=0;
 231   Standard_Real tempBinfVolume=0;
 232   Standard_Real tempCVolume = 0;
 233
 234   Binf = Zmin;
 235   Bsup = Zmax;
 236   if(Binf>Bsup)
 237     {
 238       MESSAGE("error, Bound + < Bound - in dichotomy")
 239       return -1;
 240     }
 241   tempBsupVolume = CalculateVolume(Bsup);
 242   tempBinfVolume = CalculateVolume(Binf);
 243
 244   if (Constante>tempBsupVolume || Constante<tempBinfVolume)
 245     {
 246       MESSAGE("error, algorithm start Impossible. Wrong constant value" )
 247       return -1;
 248     }
 249
 250   c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
 251     /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
 252   tempCVolume = CalculateVolume(c);
 253
 254
 255   if(Abs(tempCVolume-Constante)<=Epsilon)
 256     {
 257       goto endMethod;
 258     }
 259   else
 260     {
 261       while((Bsup-Binf)>Epsilon)
 262         {
 263           if((tempBinfVolume-Constante)*(tempCVolume-Constante)>0 && Abs(tempCVolume-Constante)>Epsilon)
 264             {
 265               Binf = c;
 266               tempBinfVolume=tempCVolume;
 267
 268               c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
 269                 /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
 270  tempCVolume=CalculateVolume(c);
 271             }
 272           else if((tempBinfVolume-Constante)*(tempCVolume-Constante)<0 && Abs(tempCVolume-Constante)>Epsilon)
 273             {
 274               Bsup = c;
 275               tempBsupVolume =tempCVolume;
 276
 277               c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
 278                 /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
 279  tempCVolume=CalculateVolume(c);
 280             }
 281           else
 282             {
 283               goto endMethod;
 284             }
 285         }
 286       goto endMethod;
 287
 288     }
 289  endMethod:
 290   MESSAGE("La ligne de flottaison correspondant a la constante :"<<Constante<<" est a la cote Z = "<<c)
 291
 292   return c;
 293 }
 294
 295 void VolumeSection::MakeRotation(gp_Dir PlaneDirection)
 296 {
 297   gp_Dir Zdirection(0.0,0.0,1.0);
 298   Standard_Real VariationAngle = 0;
 299   gp_Pnt RotationAxeLocation(0.0,0.0,0.0);
 300   gp_Dir RotationAxeDirection(1.0,1.0,1.0);
 301   gp_Ax1 RotationAxe(RotationAxeLocation,RotationAxeDirection);
 302   gp_Trsf Transformation;
 303
 304   VariationAngle = Zdirection.Angle(PlaneDirection);
 305   RotationAxe.SetDirection(PlaneDirection.Crossed(Zdirection));
 306   Transformation.SetRotation(RotationAxe,VariationAngle);
 307   TopLoc_Location L(Transformation);
 308   myShape.Move(L);
 309   myPlane->Transform(Transformation);
 310 }
 311
 312 Handle (Geom_RectangularTrimmedSurface) VolumeSection::TrimSurf()
 313 {
 314   Standard_Real Umin,Umax,Vmin,Vmax;
 315   gp_Pnt Pmin(Xmin,Ymin,Zmin);
 316   GeomAPI_ProjectPointOnSurf Projection(Pmin,myPlane);
 317   Projection.Parameters(1,Umin,Vmin);
 318   gp_Pnt Pmax(Xmax,Ymax,Zmax);
 319   GeomAPI_ProjectPointOnSurf Projection2(Pmax,myPlane);
 320   Projection2.Parameters(1,Umax,Vmax);
 321   Handle (Geom_RectangularTrimmedSurface) Plane = new Geom_RectangularTrimmedSurface(myPlane,Umin,Umax,Vmin,Vmax);
 322   return Plane;
 323 }
 324
 325 Handle (Geom_RectangularTrimmedSurface) VolumeSection::InvMakeRotation(gp_Dir PlaneDirection, Handle (Geom_RectangularTrimmedSurface) SurfTrim)
 326 {
 327   gp_Dir Zdirection(0.0,0.0,1.0);
 328   Standard_Real VariationAngle = 0;
 329   gp_Pnt RotationAxeLocation(0.0,0.0,0.0);
 330   gp_Dir RotationAxeDirection(1.0,1.0,1.0);
 331   gp_Ax1 RotationAxe(RotationAxeLocation,RotationAxeDirection);
 332   gp_Trsf Transformation;
 333
 334   VariationAngle = Zdirection.Angle(PlaneDirection);
 335   RotationAxe.SetDirection(PlaneDirection.Crossed(Zdirection));
 336   Transformation.SetRotation(RotationAxe,-VariationAngle);
 337   SurfTrim->Transform(Transformation);
 338   TopLoc_Location L(Transformation);
 339   myShape.Move(L);
 340
 341   return SurfTrim;
 342 }
 343
 344 Handle (Geom_RectangularTrimmedSurface) VolumeSection::AjustePlan(Handle (Geom_RectangularTrimmedSurface) SurfTrim, Standard_Real Cote, gp_Pnt PosPlan)
 345 {
 346   gp_Trsf Transformation;
 347   gp_Pnt PosArchi(PosPlan.X(),PosPlan.Y(),Cote);
 348
 349   Transformation.SetTranslation(PosPlan,PosArchi);
 350   SurfTrim->Transform(Transformation);
 351
 352   return SurfTrim;
 353
 354 }
 355
 356 //-------------------------------------------------------------------------------------------------------
 357 //----------------------------------- Methodes privees ---------------------------------------------------
 358 //-------------------------------------------------------------------------------------------------------
 359
 360
 361  //Fonction calculant l'intersection de la droite passant par les points P1 et P2
 362 //avec le plan horizontal Z=Hauteur
 363 gp_Pnt VolumeSection::Intersection(gp_Pnt P1,gp_Pnt P2,Standard_Real Hauteur)
 364 {
 365   Standard_Real constante;
 366   gp_Pnt Point;
 367
 368   constante = (Hauteur-P1.Z())/(P2.Z()-P1.Z());
 369   Point.SetX(P1.X()*(1-constante) + constante*P2.X());
 370   Point.SetY(P1.Y()*(1-constante) + constante*P2.Y());
 371   Point.SetZ(Hauteur);
 372
 373   return Point;
 374 }
 375
 376 //Fonction calculant le volume \89l\89mentaire de chaque t\89traedre \80 partir de 3 points
 377 Standard_Real VolumeSection::ElementaryVolume(gp_Pnt P1,gp_Pnt P2,gp_Pnt P3)
 378 {
 379   Standard_Real Determinant;
 380
 381   math_Matrix M(1,3,1,3);
 382
 383   M(1,1)=P1.X()-InitPoint.X();
 384   M(1,2)=P2.X()-InitPoint.X();
 385   M(1,3)=P3.X()-InitPoint.X();
 386   M(2,1)=P1.Y()-InitPoint.Y();
 387   M(2,2)=P2.Y()-InitPoint.Y();
 388   M(2,3)=P3.Y()-InitPoint.Y();
 389   M(3,1)=P1.Z()-InitPoint.Z();
 390   M(3,2)=P2.Z()-InitPoint.Z();
 391   M(3,3)=P3.Z()-InitPoint.Z();
 392
 393   Determinant = (1.0/6) * M.Determinant();
 394
 395   return Determinant;
 396 }
 397
 398 void VolumeSection::getZ( double& min, double& max)
 399 {
 400   min = Zmin;
 401   max = Zmax;
 402 }
 403