43128e046c1c6a6d99a03b47008b9d3438336fac
[occt.git] / src / GProp / GProp_CGProps.gxx
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 <math.hxx>
16 #include <math_Vector.hxx>
17 #include <gp.hxx>
18 #include <gp_Vec.hxx>
19 #include <Standard_NotImplemented.hxx>
20
21 #include <TColStd_Array1OfReal.hxx>
22
23
24 GProp_CGProps::GProp_CGProps(){}
25
26 void GProp_CGProps::SetLocation(const gp_Pnt& CLocation)
27 {
28   loc = CLocation;
29 }
30
31 void GProp_CGProps::Perform (const Curve& C)
32 {
33
34   Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
35   dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
36
37   Standard_Real Lower    = Tool::FirstParameter  (C);
38   Standard_Real Upper    = Tool::LastParameter   (C);
39   Standard_Integer Order = Min(Tool::IntegrationOrder (C),
40                                math::GaussPointsMax());
41   
42   gp_Pnt P;    //value on the curve
43   gp_Vec V1;   //first derivative on the curve
44   Standard_Real ds;  //curvilign abscissae
45   Standard_Real ur, um, u;
46   Standard_Real x, y, z; 
47   Standard_Real xloc, yloc, zloc;
48
49   math_Vector GaussP (1, Order);
50   math_Vector GaussW (1, Order);
51   
52   //Recuperation des points de Gauss dans le fichier GaussPoints.
53   math::GaussPoints  (Order,GaussP);
54   math::GaussWeights (Order,GaussW);
55
56   // modified by NIZHNY-MKK  Thu Jun  9 12:13:21 2005.BEGIN
57   Standard_Integer nbIntervals = Tool::NbIntervals(C, GeomAbs_CN);
58   Standard_Boolean bHasIntervals = (nbIntervals > 1);
59   TColStd_Array1OfReal TI(1, nbIntervals + 1);
60
61   if(bHasIntervals) {
62     Tool::Intervals(C, TI, GeomAbs_CN);
63   }
64   else {
65     nbIntervals = 1;
66   }
67   Standard_Integer nIndex = 0;
68   Standard_Real UU1 = Min(Lower, Upper);
69   Standard_Real UU2 = Max(Lower, Upper);
70   
71   for(nIndex = 1; nIndex <= nbIntervals; nIndex++) {
72     if(bHasIntervals) {
73       Lower = Max(TI(nIndex), UU1);
74       Upper = Min(TI(nIndex+1), UU2);
75     }
76     else {
77       Lower = UU1;
78       Upper = UU2;
79     }
80
81     Standard_Real dimLocal, IxLocal, IyLocal, IzLocal, IxxLocal, IyyLocal, IzzLocal, IxyLocal, IxzLocal, IyzLocal;
82     dimLocal = IxLocal = IyLocal = IzLocal = IxxLocal = IyyLocal = IzzLocal = IxyLocal = IxzLocal = IyzLocal = 0.0;
83   // modified by NIZHNY-MKK  Thu Jun  9 12:13:32 2005.END
84
85     loc.Coord (xloc, yloc, zloc);
86
87     Standard_Integer i;
88
89     // Calcul des integrales aux points de gauss :
90     um = 0.5 * (Upper + Lower);
91     ur = 0.5 * (Upper - Lower);
92
93     for (i = 1; i <= Order; i++) {
94       u   = um + ur * GaussP (i);
95       Tool::D1 (C,u, P, V1); 
96       ds  = V1.Magnitude();
97       P.Coord (x, y, z);
98       x   -= xloc;
99       y   -= yloc;
100       z   -= zloc;
101       ds  *= GaussW (i);
102       dimLocal += ds; 
103       IxLocal  += x * ds;  
104       IyLocal  += y * ds;
105       IzLocal  += z * ds;
106       IxyLocal += x * y * ds;
107       IyzLocal += y * z * ds;
108       IxzLocal += x * z * ds;
109       x   *= x;      
110       y   *= y;      
111       z   *= z;      
112       IxxLocal += (y + z) * ds;
113       IyyLocal += (x + z) * ds;
114       IzzLocal += (x + y) * ds;
115     }
116     // modified by NIZHNY-MKK  Thu Jun  9 12:13:47 2005.BEGIN
117     dimLocal *= ur;
118     IxLocal  *= ur;
119     IyLocal  *= ur;
120     IzLocal  *= ur;
121     IxxLocal *= ur;
122     IyyLocal *= ur;
123     IzzLocal *= ur;
124     IxyLocal *= ur;
125     IxzLocal *= ur;
126     IyzLocal *= ur;
127
128     dim += dimLocal;
129     Ix += IxLocal;
130     Iy += IyLocal;
131     Iz += IzLocal;
132     Ixx += IxxLocal;
133     Iyy += IyyLocal;
134     Izz += IzzLocal;
135     Ixy += IxyLocal;
136     Ixz += IxzLocal;
137     Iyz += IyzLocal;
138   }
139   // modified by NIZHNY-MKK  Thu Jun  9 12:13:55 2005.END
140
141   inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
142                     gp_XYZ (-Ixy, Iyy, -Iyz),
143                     gp_XYZ (-Ixz, -Iyz, Izz));
144
145   if (Abs(dim) < gp::Resolution())
146     g = P;
147   else
148     g.SetCoord (Ix/dim, Iy/dim, Iz/dim);
149 }
150
151
152 GProp_CGProps::GProp_CGProps (const Curve& C, 
153                               const gp_Pnt&   CLocation)
154 {
155   SetLocation(CLocation);
156   Perform(C);
157 }