0024750: Replace instantiations of TCollection generic classes by NCollection templat...
[occt.git] / src / GProp / GProp_GProps.cxx
1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15 #include <GProp_GProps.ixx>
16 #include <GProp.hxx>
17 #include <math_Jacobi.hxx>
18 #include <gp.hxx>
19 #include <gp.hxx>
20 #include <gp_XYZ.hxx>
21 #include <gp_Vec.hxx>
22
23
24 GProp_GProps::GProp_GProps () : g (gp::Origin()) , loc (gp::Origin()), dim (0.0) 
25
26   inertia = gp_Mat(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
27 }
28
29
30 GProp_GProps::GProp_GProps (const gp_Pnt& SystemLocation) : 
31        g (gp::Origin()), loc (SystemLocation), dim (0.0)
32
33   inertia = gp_Mat(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
34 }
35
36
37
38 void GProp_GProps::Add (const GProp_GProps& Item, const Standard_Real Density) {
39   if (Density <= gp::Resolution())  Standard_DomainError::Raise();
40   if (loc.Distance (Item.loc) <= gp::Resolution ()) {
41     gp_XYZ GXYZ = (Item.g.XYZ()).Multiplied (Item.dim * Density);
42     g.SetXYZ (g.XYZ().Multiplied (dim));
43     GXYZ.Add (g.XYZ());
44     dim = dim + Item.dim * Density;
45     if (Abs(dim) >= 1.e-20 ) {
46       GXYZ.Divide (dim);
47       g.SetXYZ (GXYZ);
48     }
49     else {
50       g.SetCoord(0., 0., 0.);
51     }
52     inertia = inertia + Item.inertia * Density;
53   }else{
54     gp_XYZ Itemloc = loc.XYZ() - Item.loc.XYZ();
55     gp_XYZ Itemg   = Item.loc.XYZ() + Item.g.XYZ();
56     gp_XYZ GXYZ    = Item.g.XYZ() - Itemloc;
57     GXYZ    = GXYZ.Multiplied (Item.dim * Density);
58     g.SetXYZ (g.XYZ().Multiplied (dim));
59     GXYZ.Add (g.XYZ());
60     dim = dim + Item.dim * Density;
61     if (Abs(dim) >= 1.e-20 ) {
62       GXYZ.Divide (dim);
63       g.SetXYZ (GXYZ);
64     }
65     else {
66       g.SetCoord(0., 0., 0.);
67     }
68     //We have to compute the inertia of the Item at the location point
69     //of the system using the Huyghens theorem
70     gp_Mat HMat;
71     gp_Mat ItemInertia = Item.inertia;
72     if (Item.g.XYZ().Modulus() > gp::Resolution()) {
73       //Computes the inertia of Item at its dim centre
74       GProp::HOperator (Itemg, Item.loc, Item.dim, HMat);
75       ItemInertia = ItemInertia - HMat;
76     }
77     //Computes the inertia of Item at the location point of the system
78     GProp::HOperator (Itemg, loc, Item.dim, HMat);
79     ItemInertia = ItemInertia + HMat;
80     inertia = inertia + ItemInertia * Density;
81   }
82 }
83
84 Standard_Real GProp_GProps::Mass () const { return dim;  }
85
86 gp_Pnt GProp_GProps::CentreOfMass () const 
87
88   return gp_Pnt(loc.XYZ() + g.XYZ()); 
89 }
90
91 gp_Mat GProp_GProps::MatrixOfInertia () const { 
92   gp_Mat HMat;
93   GProp::HOperator (g,gp::Origin(),dim, HMat);
94   return inertia - HMat;
95 }
96
97
98
99 void GProp_GProps::StaticMoments (Standard_Real& Ix, 
100                                   Standard_Real& Iy, 
101                                   Standard_Real& Iz) const {
102
103   gp_XYZ G = loc.XYZ() + g.XYZ();
104   Ix = G.X() * dim;
105   Iy = G.Y() * dim;
106   Iz = G.Z() * dim;
107 }
108
109
110
111
112 Standard_Real GProp_GProps::MomentOfInertia (const gp_Ax1& A) const {
113   // Moment of inertia / axis A
114   // 1] computes the math_Matrix of inertia / A.location()
115   // 2] applies this math_Matrix to A.Direction()
116   // 3] then computes the scalar product between this vector and 
117   //    A.Direction()
118
119       
120   if (loc.Distance (A.Location()) <= gp::Resolution()) {
121     return (A.Direction().XYZ()).Dot (
122            (A.Direction().XYZ()).Multiplied (inertia));
123   }
124   else {
125     gp_Mat HMat;
126     gp_Mat AxisInertia = MatrixOfInertia();
127     GProp::HOperator (gp_Pnt (loc.XYZ() + g.XYZ()), A.Location(), dim, HMat);
128     AxisInertia = AxisInertia + HMat;
129     return (A.Direction().XYZ()).Dot (
130            (A.Direction().XYZ()).Multiplied (AxisInertia));
131   }
132 }
133
134
135
136 Standard_Real GProp_GProps::RadiusOfGyration (const gp_Ax1& A) const {
137
138   return Sqrt (MomentOfInertia (A) / dim);
139 }
140
141
142
143 GProp_PrincipalProps GProp_GProps::PrincipalProperties () const {
144   
145   math_Matrix DiagMat (1, 3, 1, 3);
146   Standard_Integer i, j;
147   gp_Mat AxisInertia = MatrixOfInertia();
148   for (j = 1; j <= 3; j++) {
149     for (i = 1; i <= 3; i++) {
150        DiagMat (i, j) = AxisInertia.Value (i, j);
151     }
152   }
153   math_Jacobi J (DiagMat);
154   Standard_Real Ixx = J.Value (1);
155   Standard_Real Iyy = J.Value (2);
156   Standard_Real Izz = J.Value (3);
157   DiagMat  = J.Vectors ();      
158   gp_Vec Vxx (DiagMat (1,1), DiagMat (2, 1), DiagMat (3, 1));
159   gp_Vec Vyy (DiagMat (1,2), DiagMat (2, 2), DiagMat (3, 2));
160   gp_Vec Vzz (DiagMat (1,3), DiagMat (2, 3), DiagMat (3, 3));
161   //
162   // protection contre dim == 0.0e0 au cas ou on aurait rentre qu'un point
163   //
164   Standard_Real Rxx = 0.0e0 ;
165   Standard_Real Ryy = 0.0e0 ;
166   Standard_Real Rzz = 0.0e0 ;  
167   if (0.0e0 != dim) {
168     Rxx = Sqrt (Abs(Ixx / dim));
169     Ryy = Sqrt (Abs(Iyy / dim));
170     Rzz = Sqrt (Abs(Izz / dim));
171   }
172   return GProp_PrincipalProps (Ixx, Iyy, Izz, Rxx, Ryy, Rzz, Vxx, Vyy, Vzz,
173                              gp_Pnt(g.XYZ() + loc.XYZ()));
174 }
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192