src/ARCHIMEDE/Archimede_VolumeSection.cxx

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