]> SALOME platform Git repositories - modules/geom.git/blob - src/ARCHIMEDE/Archimede_VolumeSection.cxx
Salome HOME
NRI : 1_3_0 version.
[modules/geom.git] / 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 }