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