src/ARCHIMEDE/Archimede_VolumeSection.cxx

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