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