-// GEOM ARCHIMEDE : algorithm implementation
+// Copyright (C) 2007-2023 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>
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 = Tr->Node(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
- //g\89om\89trique de la boite englobante
+ // Creation du point d'initialisation, c'est e dire le centre de gravite
+ // geometrique de la boite englobante
InitPoint.SetX(0.5 * (Xmin + Xmax));
InitPoint.SetY(0.5 * (Ymin + Ymax));
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];
Standard_Real Volume=0;
Standard_Real Determinant=0;
gp_Pnt P[3];
-
+
// Projection du point d'initialisation sur le plan de section
-
+
InitPoint.SetZ(Elevation);
for (ex.Init(myShape, TopAbs_FACE); ex.More(); ex.Next())
TopoDS_Face F = TopoDS::Face(ex.Current());
Handle(Poly_Triangulation) Tr = BRep_Tool::Triangulation(F, L);
if(Tr.IsNull())
- MESSAGE("Error, null layer" )
- const Poly_Array1OfTriangle& triangles = Tr->Triangles();
+ MESSAGE("Error, null layer" );
Standard_Integer nbTriangles = Tr->NbTriangles();
- 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]);
-
- 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();
- P[2] = Nodes(noeud[2]).Transformed(L.Transformation());
+ // en tenant compte des triangles coupes par le plan de section
+
+ for (i=1;i<=nbTriangles;i++)
+ {
+ Determinant=0;
+ //Gardons la meme orientation des noeuds
+ if (F.Orientation() == TopAbs_REVERSED)
+ Tr->Triangle(i).Get(noeud[0], noeud[2], noeud[1]);
+ else
+ Tr->Triangle(i).Get(noeud[0], noeud[1], noeud[2]);
+
+ P[0] = Tr->Node(noeud[0]).Transformed(L.Transformation());
+ z[0] = P[0].Z();
+ P[1] = Tr->Node(noeud[1]).Transformed(L.Transformation());
+ z[1] = P[1].Z();
+ P[2] = Tr->Node(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;
Bsup = Zmax;
if(Binf>Bsup)
{
- MESSAGE("error, Bound + < Bound - in dichotomy")
+ MESSAGE("error, Bound + < Bound - in dichotomy");
return -1;
}
tempBsupVolume = CalculateVolume(Bsup);
if (Constante>tempBsupVolume || Constante<tempBinfVolume)
{
- MESSAGE("error, algorithm start Impossible. Wrong constant value" )
+ MESSAGE("error, algorithm start Impossible. Wrong constant value" );
return -1;
}
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;
-
- c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
- /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
+ }
+ 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));
tempCVolume=CalculateVolume(c);
- }
- else
- {
- goto endMethod;
- }
- }
+ }
+ else
+ {
+ goto endMethod;
+ }
+ }
goto endMethod;
}
endMethod:
- MESSAGE("La ligne de flottaison correspondant a la constante :"<<Constante<<" est a la cote Z = "<<c)
+ MESSAGE("La ligne de flottaison correspondant a la constante :"<<Constante<<" est a la cote Z = "<<c);
return c;
}
return Point;
}
-//Fonction calculant le volume \89l\89mentaire de chaque t\89traedre \80 partir de 3 points
+// Fonction calculant le volume elementaire de chaque tetraedre e partir de 3 points
Standard_Real VolumeSection::ElementaryVolume(gp_Pnt P1,gp_Pnt P2,gp_Pnt P3)
{
Standard_Real Determinant;
