Back to index

salome-geom  6.5.0
Archimede_VolumeSection.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 //  GEOM ARCHIMEDE : algorithm implementation
00024 //  File   : Archimede_VolumeSection.cxx
00025 //  Author : Nicolas REJNERI
00026 //  Module : GEOM
00027 //  $Header: /home/server/cvs/GEOM/GEOM_SRC/src/ARCHIMEDE/Archimede_VolumeSection.cxx,v 1.10.2.1.6.2.14.1 2012-04-13 05:47:51 vsr Exp $
00028 //
00029 #include "Archimede_VolumeSection.hxx"
00030 #include "utilities.h"
00031 
00032 #include <BRepMesh_IncrementalMesh.hxx>
00033 #include <TopExp_Explorer.hxx>
00034 #include <TopLoc_Location.hxx>
00035 #include <Poly_Triangulation.hxx>
00036 #include <Poly_Array1OfTriangle.hxx>
00037 #include <BRep_Tool.hxx>
00038 #include <TopoDS.hxx>
00039 #include <TopoDS_Face.hxx>
00040 #include <TopoDS_Shape.hxx>
00041 #include <math_Matrix.hxx>
00042 #include <gp_Trsf.hxx>
00043 #include <gp_Dir.hxx>
00044 #include <gp_Ax1.hxx>
00045 #include <gp_Pnt.hxx>
00046 
00047 #include <GeomAPI_ProjectPointOnSurf.hxx>
00048 #include <Geom_RectangularTrimmedSurface.hxx>
00049 
00050 //------------------------------------------------------------------------------------------------------- 
00051 //----------------------------------- Methodes publiques -------------------------------------------------
00052 //------------------------------------------------------------------------------------------------------- 
00053 
00054 //  Maillage de la shape
00055 VolumeSection::VolumeSection(TopoDS_Shape S , Standard_Real Precision):myShape(S),Tolerance(Precision)
00056 {
00057   // Maillage de la shape myShape
00058   BRepMesh_IncrementalMesh(myShape,Tolerance);
00059 }
00060 
00061 TopoDS_Shape VolumeSection::GetShape()
00062 {
00063   return myShape;
00064 }
00065 
00066 void VolumeSection::SetPlane(Handle (Geom_Plane) P)
00067 {
00068   myPlane = P;
00069 }
00070 
00071 void VolumeSection::CenterOfGravity()
00072 {
00073   Standard_Integer i;
00074   Standard_Integer nbNodes;
00075   TopExp_Explorer ex;
00076   TopLoc_Location L;
00077   
00078   // Boucle sur les faces de la shape
00079   
00080   Xmin = 1000000000;
00081   Ymin = 1000000000;
00082   Zmin = 1000000000;
00083   Xmax = -1000000000;
00084   Ymax = -1000000000;
00085   Zmax = -1000000000;
00086   
00087   for (ex.Init(myShape, TopAbs_FACE); ex.More(); ex.Next()) 
00088     {
00089       TopoDS_Face F = TopoDS::Face(ex.Current());
00090       Handle(Poly_Triangulation) Tr = BRep_Tool::Triangulation(F, L);
00091       if(Tr.IsNull())
00092         MESSAGE("Error, null layer" )
00093       nbNodes = Tr->NbNodes();
00094       const TColgp_Array1OfPnt& Nodes = Tr->Nodes();
00095       
00096       // Calcul des dimensions de la boite englobante du solide
00097       
00098       for(i=1;i<=nbNodes;i++)
00099         {
00100           InitPoint = Nodes(i).Transformed(L.Transformation());
00101           if(InitPoint.X() < Xmin)
00102             Xmin = InitPoint.X();
00103           if(InitPoint.X() > Xmax)
00104             Xmax = InitPoint.X();
00105           if(InitPoint.Y() < Ymin)
00106             Ymin = InitPoint.Y();
00107           if(InitPoint.Y() > Ymax)
00108             Ymax = InitPoint.Y();
00109           if(InitPoint.Z() < Zmin)
00110             Zmin = InitPoint.Z();
00111           if(InitPoint.Z() > Zmax)
00112             Zmax = InitPoint.Z();
00113           
00114         }
00115     }
00116   
00117   // Creation du point d'initialisation, c'est € dire le centre de gravit‰ 
00118   //g‰om‰trique de la boite englobante
00119   
00120   InitPoint.SetX(0.5 * (Xmin + Xmax));
00121   InitPoint.SetY(0.5 * (Ymin + Ymax));
00122   InitPoint.SetZ(0);
00123 }
00124 
00125 Standard_Real VolumeSection::CalculateVolume(Standard_Real Elevation)
00126 {
00127   Standard_Integer i,noeud[3],flag[3];
00128   Standard_Integer nbNodes;
00129   TopExp_Explorer ex;
00130   TopLoc_Location L;
00131   Standard_Real z[3];
00132   Standard_Real Volume=0;
00133   Standard_Real Determinant=0;
00134   gp_Pnt P[3];
00135   
00136   // Projection du point d'initialisation sur le plan de section
00137   
00138   InitPoint.SetZ(Elevation);
00139 
00140   for (ex.Init(myShape, TopAbs_FACE); ex.More(); ex.Next()) 
00141     {
00142       TopoDS_Face F = TopoDS::Face(ex.Current());
00143       Handle(Poly_Triangulation) Tr = BRep_Tool::Triangulation(F, L);
00144       if(Tr.IsNull())
00145         MESSAGE("Error, null layer" )
00146       const Poly_Array1OfTriangle& triangles = Tr->Triangles();
00147       Standard_Integer nbTriangles = Tr->NbTriangles();
00148       nbNodes = Tr->NbNodes();
00149       const TColgp_Array1OfPnt& Nodes = Tr->Nodes();
00150       
00151       // Calcul des volumes de chaque triangle, de chaque face 
00152       //en tenant compte des triangles coup‰s par le plan de section
00153       
00154       for (i=1;i<=nbTriangles;i++) 
00155         {
00156           Determinant=0;
00157           //Gardons la meme orientation des noeuds
00158           if (F.Orientation()  == TopAbs_REVERSED)
00159             triangles(i).Get(noeud[0], noeud[2], noeud[1]);
00160           else 
00161             triangles(i).Get(noeud[0], noeud[1], noeud[2]);
00162           
00163           P[0] = Nodes(noeud[0]).Transformed(L.Transformation());
00164           z[0] = P[0].Z();
00165           P[1] = Nodes(noeud[1]).Transformed(L.Transformation());
00166           z[1] = P[1].Z();
00167           P[2] = Nodes(noeud[2]).Transformed(L.Transformation());
00168           z[2] = P[2].Z();
00169 
00170           // Determination des cas aux limites pour les triangles
00171           Standard_Integer i,compteur=0;
00172 
00173           for (i=0;i<=2;i++)
00174             {
00175               flag[i]=Standard_False;
00176               if(z[i]>=Elevation)
00177                 {
00178                   flag[i]=Standard_True;
00179                   compteur++;
00180                 }
00181             }
00182           
00183           switch(compteur)
00184             {
00185             case 0:
00186               Determinant = ElementaryVolume(P[0],P[1],P[2]);
00187               break;
00188               
00189             case 1:
00190               for (i=0;i<=2;i++)
00191                 {
00192                   if (flag[i]==Standard_True)
00193                     {
00194                       gp_Pnt Result1 = Intersection(P[i],P[(i+1)%3],Elevation);
00195                       gp_Pnt Result2 = Intersection(P[i],P[(i+2)%3],Elevation);
00196                       Determinant = ElementaryVolume(Result1,P[(i+1)%3],P[(i+2)%3])
00197                         + ElementaryVolume(Result1,P[(i+2)%3],Result2);
00198                     }
00199                 }
00200               break;
00201               
00202             case 2:
00203               for (i=0;i<=2;i++)
00204                 {
00205                   if (flag[i]==Standard_False)
00206                     {
00207                       gp_Pnt Result1 = Intersection(P[i],P[(i+1)%3],Elevation);
00208                       gp_Pnt Result2 = Intersection(P[i],P[(i+2)%3],Elevation);
00209                       Determinant = ElementaryVolume(P[i],Result1,Result2);
00210                     }
00211                 }
00212               break;
00213               
00214             case 3:
00215               break;
00216             }
00217           Volume += Determinant;
00218         }
00219     }
00220   
00221   return Volume;
00222 }
00223 
00224 Standard_Real VolumeSection::Archimede(Standard_Real Constante , Standard_Real Epsilon)
00225 {
00226   // Resolution de l equation V(h) = Constante a l aide de l algorithme de dichotomie avec ponderation type
00227   // Lagrange
00228   
00229   Standard_Real c,Binf,Bsup;
00230   Standard_Real tempBsupVolume=0;
00231   Standard_Real tempBinfVolume=0;
00232   Standard_Real tempCVolume = 0;
00233 
00234   Binf = Zmin;
00235   Bsup = Zmax;
00236   if(Binf>Bsup)
00237     {
00238       MESSAGE("error, Bound + < Bound - in dichotomy")
00239       return -1;
00240     }
00241   tempBsupVolume = CalculateVolume(Bsup);
00242   tempBinfVolume = CalculateVolume(Binf);
00243   
00244   if (Constante>tempBsupVolume || Constante<tempBinfVolume)
00245     {
00246       MESSAGE("error, algorithm start Impossible. Wrong constant value" )
00247       return -1;
00248     }
00249   
00250   c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
00251     /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
00252   tempCVolume = CalculateVolume(c);
00253   
00254   
00255   if(Abs(tempCVolume-Constante)<=Epsilon)
00256     {
00257       goto endMethod;
00258     }
00259   else
00260     {
00261       while((Bsup-Binf)>Epsilon)
00262         { 
00263           if((tempBinfVolume-Constante)*(tempCVolume-Constante)>0 && Abs(tempCVolume-Constante)>Epsilon)
00264             {
00265               Binf = c;
00266               tempBinfVolume=tempCVolume;
00267               
00268               c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
00269                 /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
00270  tempCVolume=CalculateVolume(c);
00271             }
00272           else if((tempBinfVolume-Constante)*(tempCVolume-Constante)<0 && Abs(tempCVolume-Constante)>Epsilon)
00273             {
00274               Bsup = c;
00275               tempBsupVolume =tempCVolume;
00276 
00277               c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
00278                 /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
00279  tempCVolume=CalculateVolume(c);
00280             }
00281           else
00282             {
00283               goto endMethod;
00284             }
00285         }
00286       goto endMethod;
00287       
00288     }
00289  endMethod:
00290   MESSAGE("La ligne de flottaison correspondant a la constante :"<<Constante<<" est a la cote Z = "<<c)
00291   
00292   return c;
00293 }
00294 
00295 void VolumeSection::MakeRotation(gp_Dir PlaneDirection)
00296 {
00297   gp_Dir Zdirection(0.0,0.0,1.0);
00298   Standard_Real VariationAngle = 0;
00299   gp_Pnt RotationAxeLocation(0.0,0.0,0.0);
00300   gp_Dir RotationAxeDirection(1.0,1.0,1.0);
00301   gp_Ax1 RotationAxe(RotationAxeLocation,RotationAxeDirection);
00302   gp_Trsf Transformation;
00303   
00304   VariationAngle = Zdirection.Angle(PlaneDirection);
00305   RotationAxe.SetDirection(PlaneDirection.Crossed(Zdirection));
00306   Transformation.SetRotation(RotationAxe,VariationAngle);
00307   TopLoc_Location L(Transformation);
00308   myShape.Move(L);
00309   myPlane->Transform(Transformation);
00310 }
00311 
00312 Handle (Geom_RectangularTrimmedSurface) VolumeSection::TrimSurf()
00313 {
00314   Standard_Real Umin,Umax,Vmin,Vmax;
00315   gp_Pnt Pmin(Xmin,Ymin,Zmin);
00316   GeomAPI_ProjectPointOnSurf Projection(Pmin,myPlane);
00317   Projection.Parameters(1,Umin,Vmin);
00318   gp_Pnt Pmax(Xmax,Ymax,Zmax);
00319   GeomAPI_ProjectPointOnSurf Projection2(Pmax,myPlane);
00320   Projection2.Parameters(1,Umax,Vmax);
00321   Handle (Geom_RectangularTrimmedSurface) Plane = new Geom_RectangularTrimmedSurface(myPlane,Umin,Umax,Vmin,Vmax);
00322   return Plane;
00323 }
00324 
00325 Handle (Geom_RectangularTrimmedSurface) VolumeSection::InvMakeRotation(gp_Dir PlaneDirection, Handle (Geom_RectangularTrimmedSurface) SurfTrim)
00326 {
00327   gp_Dir Zdirection(0.0,0.0,1.0);
00328   Standard_Real VariationAngle = 0;
00329   gp_Pnt RotationAxeLocation(0.0,0.0,0.0);
00330   gp_Dir RotationAxeDirection(1.0,1.0,1.0);
00331   gp_Ax1 RotationAxe(RotationAxeLocation,RotationAxeDirection);
00332   gp_Trsf Transformation;
00333 
00334   VariationAngle = Zdirection.Angle(PlaneDirection);
00335   RotationAxe.SetDirection(PlaneDirection.Crossed(Zdirection));
00336   Transformation.SetRotation(RotationAxe,-VariationAngle);
00337   SurfTrim->Transform(Transformation);
00338   TopLoc_Location L(Transformation);
00339   myShape.Move(L);
00340   
00341   return SurfTrim;
00342 }
00343 
00344 Handle (Geom_RectangularTrimmedSurface) VolumeSection::AjustePlan(Handle (Geom_RectangularTrimmedSurface) SurfTrim, Standard_Real Cote, gp_Pnt PosPlan)
00345 {
00346   gp_Trsf Transformation;
00347   gp_Pnt PosArchi(PosPlan.X(),PosPlan.Y(),Cote);
00348   
00349   Transformation.SetTranslation(PosPlan,PosArchi);
00350   SurfTrim->Transform(Transformation);
00351   
00352   return SurfTrim;
00353       
00354 }
00355 
00356 //------------------------------------------------------------------------------------------------------- 
00357 //----------------------------------- Methodes privees ---------------------------------------------------
00358 //------------------------------------------------------------------------------------------------------- 
00359 
00360 
00361  //Fonction calculant l'intersection de la droite passant par les points P1 et P2 
00362 //avec le plan horizontal Z=Hauteur
00363 gp_Pnt VolumeSection::Intersection(gp_Pnt P1,gp_Pnt P2,Standard_Real Hauteur)
00364 {
00365   Standard_Real constante;
00366   gp_Pnt Point;
00367 
00368   constante = (Hauteur-P1.Z())/(P2.Z()-P1.Z());
00369   Point.SetX(P1.X()*(1-constante) + constante*P2.X());
00370   Point.SetY(P1.Y()*(1-constante) + constante*P2.Y());
00371   Point.SetZ(Hauteur);
00372   
00373   return Point;
00374 }
00375 
00376 //Fonction calculant le volume ‰l‰mentaire de chaque t‰traedre € partir de 3 points
00377 Standard_Real VolumeSection::ElementaryVolume(gp_Pnt P1,gp_Pnt P2,gp_Pnt P3)
00378 {
00379   Standard_Real Determinant;
00380   
00381   math_Matrix M(1,3,1,3);
00382   
00383   M(1,1)=P1.X()-InitPoint.X();
00384   M(1,2)=P2.X()-InitPoint.X();
00385   M(1,3)=P3.X()-InitPoint.X();
00386   M(2,1)=P1.Y()-InitPoint.Y();
00387   M(2,2)=P2.Y()-InitPoint.Y();
00388   M(2,3)=P3.Y()-InitPoint.Y();
00389   M(3,1)=P1.Z()-InitPoint.Z();
00390   M(3,2)=P2.Z()-InitPoint.Z();
00391   M(3,3)=P3.Z()-InitPoint.Z();
00392   
00393   Determinant = (1.0/6) * M.Determinant();
00394 
00395   return Determinant;
00396 }
00397 
00398 void VolumeSection::getZ( double& min, double& max)
00399 {
00400   min = Zmin;
00401   max = Zmax;
00402 }
00403