Salome HOME
Copyright update 2020
[modules/geom.git] / src / ARCHIMEDE / Archimede_VolumeSection.cxx
index 9558be90f4c8a9aad9272459329ab5d9835653f8..5cf4f341b5176027b2c6897ccca7960d2abd515f 100644 (file)
@@ -1,33 +1,30 @@
-//  GEOM ARCHIMEDE : algorithm implementation
+// Copyright (C) 2007-2020  CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
 //
-//  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
-//  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS 
-// 
-//  This library is free software; you can redistribute it and/or 
-//  modify it under the terms of the GNU Lesser General Public 
-//  License as published by the Free Software Foundation; either 
-//  version 2.1 of the License. 
-// 
-//  This library is distributed in the hope that it will be useful, 
-//  but WITHOUT ANY WARRANTY; without even the implied warranty of 
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
-//  Lesser General Public License for more details. 
-// 
-//  You should have received a copy of the GNU Lesser General Public 
-//  License along with this library; if not, write to the Free Software 
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA 
-// 
-//  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org 
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
 //
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
+
+//  GEOM ARCHIMEDE : algorithm implementation
 //  File   : Archimede_VolumeSection.cxx
 //  Author : Nicolas REJNERI
 //  Module : GEOM
-//  $Header$
-
-#include <Standard_OStream.hxx>
-
+//
 #include "Archimede_VolumeSection.hxx"
 #include "utilities.h"
 
 #include <TopoDS_Face.hxx>
 #include <TopoDS_Shape.hxx>
 #include <math_Matrix.hxx>
-#include <math.h>
-#include <GC_MakePlane.hxx>
-#include <stdlib.h>
 #include <gp_Trsf.hxx>
 #include <gp_Dir.hxx>
 #include <gp_Ax1.hxx>
 #include <gp_Pnt.hxx>
-#include <gp_Pln.hxx>
 
 #include <GeomAPI_ProjectPointOnSurf.hxx>
 #include <Geom_RectangularTrimmedSurface.hxx>
@@ -95,29 +88,29 @@ void VolumeSection::CenterOfGravity()
       TopoDS_Face F = TopoDS::Face(ex.Current());
       Handle(Poly_Triangulation) Tr = BRep_Tool::Triangulation(F, L);
       if(Tr.IsNull())
-       MESSAGE("Error, null layer" )
+        MESSAGE("Error, null layer" )
       nbNodes = Tr->NbNodes();
       const TColgp_Array1OfPnt& Nodes = Tr->Nodes();
       
       // Calcul des dimensions de la boite englobante du solide
       
       for(i=1;i<=nbNodes;i++)
-       {
-         InitPoint = Nodes(i).Transformed(L.Transformation());
-         if(InitPoint.X() < Xmin)
-           Xmin = InitPoint.X();
-         if(InitPoint.X() > Xmax)
-           Xmax = InitPoint.X();
-         if(InitPoint.Y() < Ymin)
-           Ymin = InitPoint.Y();
-         if(InitPoint.Y() > Ymax)
-           Ymax = InitPoint.Y();
-         if(InitPoint.Z() < Zmin)
-           Zmin = InitPoint.Z();
-         if(InitPoint.Z() > Zmax)
-           Zmax = InitPoint.Z();
-         
-       }
+        {
+          InitPoint = Nodes(i).Transformed(L.Transformation());
+          if(InitPoint.X() < Xmin)
+            Xmin = InitPoint.X();
+          if(InitPoint.X() > Xmax)
+            Xmax = InitPoint.X();
+          if(InitPoint.Y() < Ymin)
+            Ymin = InitPoint.Y();
+          if(InitPoint.Y() > Ymax)
+            Ymax = InitPoint.Y();
+          if(InitPoint.Z() < Zmin)
+            Zmin = InitPoint.Z();
+          if(InitPoint.Z() > Zmax)
+            Zmax = InitPoint.Z();
+          
+        }
     }
   
   // Creation du point d'initialisation, c'est \80 dire le centre de gravit\89 
@@ -131,7 +124,7 @@ void VolumeSection::CenterOfGravity()
 Standard_Real VolumeSection::CalculateVolume(Standard_Real Elevation)
 {
   Standard_Integer i,noeud[3],flag[3];
-  Standard_Integer nbNodes;
+  //Standard_Integer nbNodes;
   TopExp_Explorer ex;
   TopLoc_Location L;
   Standard_Real z[3];
@@ -148,80 +141,80 @@ Standard_Real VolumeSection::CalculateVolume(Standard_Real Elevation)
       TopoDS_Face F = TopoDS::Face(ex.Current());
       Handle(Poly_Triangulation) Tr = BRep_Tool::Triangulation(F, L);
       if(Tr.IsNull())
-       MESSAGE("Error, null layer" )
+        MESSAGE("Error, null layer" )
       const Poly_Array1OfTriangle& triangles = Tr->Triangles();
       Standard_Integer nbTriangles = Tr->NbTriangles();
-      nbNodes = Tr->NbNodes();
+      //nbNodes = Tr->NbNodes();
       const TColgp_Array1OfPnt& Nodes = Tr->Nodes();
       
       // Calcul des volumes de chaque triangle, de chaque face 
       //en tenant compte des triangles coup\89s par le plan de section
       
       for (i=1;i<=nbTriangles;i++) 
-       {
-         Determinant=0;
-         //Gardons la meme orientation des noeuds
-         if (F.Orientation()  == TopAbs_REVERSED)
-           triangles(i).Get(noeud[0], noeud[2], noeud[1]);
-         else 
-           triangles(i).Get(noeud[0], noeud[1], noeud[2]);
-         
+        {
+          Determinant=0;
+          //Gardons la meme orientation des noeuds
+          if (F.Orientation()  == TopAbs_REVERSED)
+            triangles(i).Get(noeud[0], noeud[2], noeud[1]);
+          else 
+            triangles(i).Get(noeud[0], noeud[1], noeud[2]);
+          
           P[0] = Nodes(noeud[0]).Transformed(L.Transformation());
-         z[0] = P[0].Z();
-         P[1] = Nodes(noeud[1]).Transformed(L.Transformation());
-         z[1] = P[1].Z();
+          z[0] = P[0].Z();
+          P[1] = Nodes(noeud[1]).Transformed(L.Transformation());
+          z[1] = P[1].Z();
           P[2] = Nodes(noeud[2]).Transformed(L.Transformation());
           z[2] = P[2].Z();
 
-         // Determination des cas aux limites pour les triangles
-         Standard_Integer i,compteur=0;
+          // Determination des cas aux limites pour les triangles
+          Standard_Integer i,compteur=0;
 
-         for (i=0;i<=2;i++)
-           {
+          for (i=0;i<=2;i++)
+            {
               flag[i]=Standard_False;
-             if(z[i]>=Elevation)
-               {
-                 flag[i]=Standard_True;
-                 compteur++;
-               }
-           }
-         
-         switch(compteur)
-           {
-           case 0:
-             Determinant = ElementaryVolume(P[0],P[1],P[2]);
-             break;
-             
-           case 1:
-             for (i=0;i<=2;i++)
-               {
-                 if (flag[i]==Standard_True)
-                   {
-                     gp_Pnt Result1 = Intersection(P[i],P[(i+1)%3],Elevation);
-                     gp_Pnt Result2 = Intersection(P[i],P[(i+2)%3],Elevation);
-                     Determinant = ElementaryVolume(Result1,P[(i+1)%3],P[(i+2)%3])
-                       + ElementaryVolume(Result1,P[(i+2)%3],Result2);
-                   }
-               }
-             break;
-             
-           case 2:
-             for (i=0;i<=2;i++)
-               {
-                 if (flag[i]==Standard_False)
-                   {
-                     gp_Pnt Result1 = Intersection(P[i],P[(i+1)%3],Elevation);
-                     gp_Pnt Result2 = Intersection(P[i],P[(i+2)%3],Elevation);
-                     Determinant = ElementaryVolume(P[i],Result1,Result2);
-                   }
-               }
-             break;
-             
-           case 3:
-             break;
-           }
-         Volume += Determinant;
-       }
+              if(z[i]>=Elevation)
+                {
+                  flag[i]=Standard_True;
+                  compteur++;
+                }
+            }
+          
+          switch(compteur)
+            {
+            case 0:
+              Determinant = ElementaryVolume(P[0],P[1],P[2]);
+              break;
+              
+            case 1:
+              for (i=0;i<=2;i++)
+                {
+                  if (flag[i]==Standard_True)
+                    {
+                      gp_Pnt Result1 = Intersection(P[i],P[(i+1)%3],Elevation);
+                      gp_Pnt Result2 = Intersection(P[i],P[(i+2)%3],Elevation);
+                      Determinant = ElementaryVolume(Result1,P[(i+1)%3],P[(i+2)%3])
+                        + ElementaryVolume(Result1,P[(i+2)%3],Result2);
+                    }
+                }
+              break;
+              
+            case 2:
+              for (i=0;i<=2;i++)
+                {
+                  if (flag[i]==Standard_False)
+                    {
+                      gp_Pnt Result1 = Intersection(P[i],P[(i+1)%3],Elevation);
+                      gp_Pnt Result2 = Intersection(P[i],P[(i+2)%3],Elevation);
+                      Determinant = ElementaryVolume(P[i],Result1,Result2);
+                    }
+                }
+              break;
+              
+            case 3:
+              break;
+            }
+          Volume += Determinant;
+        }
     }
   
   return Volume;
@@ -265,30 +258,30 @@ Standard_Real VolumeSection::Archimede(Standard_Real Constante , Standard_Real E
   else
     {
       while((Bsup-Binf)>Epsilon)
-       { 
-         if((tempBinfVolume-Constante)*(tempCVolume-Constante)>0 && Abs(tempCVolume-Constante)>Epsilon)
-           {
-             Binf = c;
-             tempBinfVolume=tempCVolume;
-             
-             c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
-               /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
+        
+          if((tempBinfVolume-Constante)*(tempCVolume-Constante)>0 && Abs(tempCVolume-Constante)>Epsilon)
+            {
+              Binf = c;
+              tempBinfVolume=tempCVolume;
+              
+              c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
+                /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
  tempCVolume=CalculateVolume(c);
-           }
-         else if((tempBinfVolume-Constante)*(tempCVolume-Constante)<0 && Abs(tempCVolume-Constante)>Epsilon)
-           {
-             Bsup = c;
-             tempBsupVolume =tempCVolume;
+            }
+          else if((tempBinfVolume-Constante)*(tempCVolume-Constante)<0 && Abs(tempCVolume-Constante)>Epsilon)
+            {
+              Bsup = c;
+              tempBsupVolume =tempCVolume;
 
-             c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
-               /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
+              c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
+                /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
  tempCVolume=CalculateVolume(c);
-           }
-         else
-           {
-             goto endMethod;
-           }
-       }
+            }
+          else
+            {
+              goto endMethod;
+            }
+        }
       goto endMethod;
       
     }