2 // File : Archimede_VolumeSection.cxx
3 // Created : Fri Feb 22 09:28:13 CET 2002
6 // Modified : Fri Feb 22 09:28:13 CET 2002
7 // Author : Nicolas REJNERI
10 // Copyright : Open CASCADE 2002
13 #include "Archimede_VolumeSection.hxx"
14 #include "utilities.h"
17 #include <BRepMesh_IncrementalMesh.hxx>
18 #include <TopExp_Explorer.hxx>
19 #include <TopLoc_Location.hxx>
20 #include <Poly_Triangulation.hxx>
21 #include <Poly_Array1OfTriangle.hxx>
22 #include <BRep_Tool.hxx>
24 #include <TopoDS_Face.hxx>
25 #include <TopoDS_Shape.hxx>
26 #include <math_Matrix.hxx>
28 #include <GC_MakePlane.hxx>
30 #include <gp_Trsf.hxx>
36 #include <GeomAPI_ProjectPointOnSurf.hxx>
37 #include <Geom_RectangularTrimmedSurface.hxx>
39 //-------------------------------------------------------------------------------------------------------
40 //----------------------------------- Methodes publiques -------------------------------------------------
41 //-------------------------------------------------------------------------------------------------------
43 // Maillage de la shape
44 VolumeSection::VolumeSection(TopoDS_Shape S , Standard_Real Precision):myShape(S),Tolerance(Precision)
46 // Maillage de la shape myShape
47 BRepMesh_IncrementalMesh(myShape,Tolerance);
50 TopoDS_Shape VolumeSection::GetShape()
55 void VolumeSection::SetPlane(Handle (Geom_Plane) P)
60 void VolumeSection::CenterOfGravity()
63 Standard_Integer nbNodes;
67 // Boucle sur les faces de la shape
76 for (ex.Init(myShape, TopAbs_FACE); ex.More(); ex.Next())
78 TopoDS_Face F = TopoDS::Face(ex.Current());
79 Handle(Poly_Triangulation) Tr = BRep_Tool::Triangulation(F, L);
81 MESSAGE("Error, null layer" )
82 nbNodes = Tr->NbNodes();
83 const TColgp_Array1OfPnt& Nodes = Tr->Nodes();
85 // Calcul des dimensions de la boite englobante du solide
87 for(i=1;i<=nbNodes;i++)
89 InitPoint = Nodes(i).Transformed(L.Transformation());
90 if(InitPoint.X() < Xmin)
92 if(InitPoint.X() > Xmax)
94 if(InitPoint.Y() < Ymin)
96 if(InitPoint.Y() > Ymax)
98 if(InitPoint.Z() < Zmin)
100 if(InitPoint.Z() > Zmax)
101 Zmax = InitPoint.Z();
106 // Creation du point d'initialisation, c'est à dire le centre de gravité
107 //géométrique de la boite englobante
109 InitPoint.SetX(0.5 * (Xmin + Xmax));
110 InitPoint.SetY(0.5 * (Ymin + Ymax));
114 Standard_Real VolumeSection::CalculateVolume(Standard_Real Elevation)
116 Standard_Integer i,noeud[3],flag[3];
117 Standard_Integer nbNodes;
121 Standard_Real Volume=0;
122 Standard_Real Determinant=0;
125 // Projection du point d'initialisation sur le plan de section
127 InitPoint.SetZ(Elevation);
129 for (ex.Init(myShape, TopAbs_FACE); ex.More(); ex.Next())
131 TopoDS_Face F = TopoDS::Face(ex.Current());
132 Handle(Poly_Triangulation) Tr = BRep_Tool::Triangulation(F, L);
134 MESSAGE("Error, null layer" )
135 const Poly_Array1OfTriangle& triangles = Tr->Triangles();
136 Standard_Integer nbTriangles = Tr->NbTriangles();
137 nbNodes = Tr->NbNodes();
138 const TColgp_Array1OfPnt& Nodes = Tr->Nodes();
140 // Calcul des volumes de chaque triangle, de chaque face
141 //en tenant compte des triangles coupés par le plan de section
143 for (i=1;i<=nbTriangles;i++)
146 //Gardons la meme orientation des noeuds
147 if (F.Orientation() == TopAbs_REVERSED)
148 triangles(i).Get(noeud[0], noeud[2], noeud[1]);
150 triangles(i).Get(noeud[0], noeud[1], noeud[2]);
152 P[0] = Nodes(noeud[0]).Transformed(L.Transformation());
154 P[1] = Nodes(noeud[1]).Transformed(L.Transformation());
156 P[2] = Nodes(noeud[2]).Transformed(L.Transformation());
159 // Determination des cas aux limites pour les triangles
160 Standard_Integer i,compteur=0;
164 flag[i]=Standard_False;
167 flag[i]=Standard_True;
175 Determinant = ElementaryVolume(P[0],P[1],P[2]);
181 if (flag[i]==Standard_True)
183 gp_Pnt Result1 = Intersection(P[i],P[(i+1)%3],Elevation);
184 gp_Pnt Result2 = Intersection(P[i],P[(i+2)%3],Elevation);
185 Determinant = ElementaryVolume(Result1,P[(i+1)%3],P[(i+2)%3])
186 + ElementaryVolume(Result1,P[(i+2)%3],Result2);
194 if (flag[i]==Standard_False)
196 gp_Pnt Result1 = Intersection(P[i],P[(i+1)%3],Elevation);
197 gp_Pnt Result2 = Intersection(P[i],P[(i+2)%3],Elevation);
198 Determinant = ElementaryVolume(P[i],Result1,Result2);
206 Volume += Determinant;
213 Standard_Real VolumeSection::Archimede(Standard_Real Constante , Standard_Real Epsilon)
215 // Resolution de l equation V(h) = Constante a l aide de l algorithme de dichotomie avec ponderation type
218 Standard_Real c,Binf,Bsup;
219 Standard_Real tempBsupVolume=0;
220 Standard_Real tempBinfVolume=0;
221 Standard_Real tempCVolume = 0;
227 MESSAGE("error, Bound + < Bound - in dichotomy")
230 tempBsupVolume = CalculateVolume(Bsup);
231 tempBinfVolume = CalculateVolume(Binf);
233 if (Constante>tempBsupVolume || Constante<tempBinfVolume)
235 MESSAGE("error, algorithm start Impossible. Wrong constant value" )
239 c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
240 /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
241 tempCVolume = CalculateVolume(c);
244 if(Abs(tempCVolume-Constante)<=Epsilon)
250 while((Bsup-Binf)>Epsilon)
252 if((tempBinfVolume-Constante)*(tempCVolume-Constante)>0 && Abs(tempCVolume-Constante)>Epsilon)
255 tempBinfVolume=tempCVolume;
257 c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
258 /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
259 tempCVolume=CalculateVolume(c);
261 else if((tempBinfVolume-Constante)*(tempCVolume-Constante)<0 && Abs(tempCVolume-Constante)>Epsilon)
264 tempBsupVolume =tempCVolume;
266 c = ((Binf*(tempBsupVolume-Constante))-(Bsup*(tempBinfVolume-Constante)))
267 /((tempBsupVolume-Constante)-(tempBinfVolume-Constante));
268 tempCVolume=CalculateVolume(c);
279 MESSAGE("La ligne de flottaison correspondant a la constante :"<<Constante<<" est a la cote Z = "<<c)
284 void VolumeSection::MakeRotation(gp_Dir PlaneDirection)
286 gp_Dir Zdirection(0.0,0.0,1.0);
287 Standard_Real VariationAngle = 0;
288 gp_Pnt RotationAxeLocation(0.0,0.0,0.0);
289 gp_Dir RotationAxeDirection(1.0,1.0,1.0);
290 gp_Ax1 RotationAxe(RotationAxeLocation,RotationAxeDirection);
291 gp_Trsf Transformation;
293 VariationAngle = Zdirection.Angle(PlaneDirection);
294 RotationAxe.SetDirection(PlaneDirection.Crossed(Zdirection));
295 Transformation.SetRotation(RotationAxe,VariationAngle);
296 TopLoc_Location L(Transformation);
298 myPlane->Transform(Transformation);
301 Handle (Geom_RectangularTrimmedSurface) VolumeSection::TrimSurf()
303 Standard_Real Umin,Umax,Vmin,Vmax;
304 gp_Pnt Pmin(Xmin,Ymin,Zmin);
305 GeomAPI_ProjectPointOnSurf Projection(Pmin,myPlane);
306 Projection.Parameters(1,Umin,Vmin);
307 gp_Pnt Pmax(Xmax,Ymax,Zmax);
308 GeomAPI_ProjectPointOnSurf Projection2(Pmax,myPlane);
309 Projection2.Parameters(1,Umax,Vmax);
310 Handle (Geom_RectangularTrimmedSurface) Plane = new Geom_RectangularTrimmedSurface(myPlane,Umin,Umax,Vmin,Vmax);
314 Handle (Geom_RectangularTrimmedSurface) VolumeSection::InvMakeRotation(gp_Dir PlaneDirection, Handle (Geom_RectangularTrimmedSurface) SurfTrim)
316 gp_Dir Zdirection(0.0,0.0,1.0);
317 Standard_Real VariationAngle = 0;
318 gp_Pnt RotationAxeLocation(0.0,0.0,0.0);
319 gp_Dir RotationAxeDirection(1.0,1.0,1.0);
320 gp_Ax1 RotationAxe(RotationAxeLocation,RotationAxeDirection);
321 gp_Trsf Transformation;
323 VariationAngle = Zdirection.Angle(PlaneDirection);
324 RotationAxe.SetDirection(PlaneDirection.Crossed(Zdirection));
325 Transformation.SetRotation(RotationAxe,-VariationAngle);
326 SurfTrim->Transform(Transformation);
327 TopLoc_Location L(Transformation);
333 Handle (Geom_RectangularTrimmedSurface) VolumeSection::AjustePlan(Handle (Geom_RectangularTrimmedSurface) SurfTrim, Standard_Real Cote, gp_Pnt PosPlan)
335 gp_Trsf Transformation;
336 gp_Pnt PosArchi(PosPlan.X(),PosPlan.Y(),Cote);
338 Transformation.SetTranslation(PosPlan,PosArchi);
339 SurfTrim->Transform(Transformation);
345 //-------------------------------------------------------------------------------------------------------
346 //----------------------------------- Methodes privees ---------------------------------------------------
347 //-------------------------------------------------------------------------------------------------------
350 //Fonction calculant l'intersection de la droite passant par les points P1 et P2
351 //avec le plan horizontal Z=Hauteur
352 gp_Pnt VolumeSection::Intersection(gp_Pnt P1,gp_Pnt P2,Standard_Real Hauteur)
354 Standard_Real constante;
357 constante = (Hauteur-P1.Z())/(P2.Z()-P1.Z());
358 Point.SetX(P1.X()*(1-constante) + constante*P2.X());
359 Point.SetY(P1.Y()*(1-constante) + constante*P2.Y());
365 //Fonction calculant le volume élémentaire de chaque tétraedre à partir de 3 points
366 Standard_Real VolumeSection::ElementaryVolume(gp_Pnt P1,gp_Pnt P2,gp_Pnt P3)
368 Standard_Real Determinant;
370 math_Matrix M(1,3,1,3);
372 M(1,1)=P1.X()-InitPoint.X();
373 M(1,2)=P2.X()-InitPoint.X();
374 M(1,3)=P3.X()-InitPoint.X();
375 M(2,1)=P1.Y()-InitPoint.Y();
376 M(2,2)=P2.Y()-InitPoint.Y();
377 M(2,3)=P3.Y()-InitPoint.Y();
378 M(3,1)=P1.Z()-InitPoint.Z();
379 M(3,2)=P2.Z()-InitPoint.Z();
380 M(3,3)=P3.Z()-InitPoint.Z();
382 Determinant = (1.0/6) * M.Determinant();
387 void VolumeSection::getZ( double& min, double& max)