src/ARCHIMEDE/Archimede_VolumeSection.cxx

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