projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0024171: Eliminate CLang compiler warning -Wreorder
[occt.git] / src / GProp / GProp_SGProps.gxx
CommitLineData
b311480e 1// Copyright (c) 1995-1999 Matra Datavision
2// Copyright (c) 1999-2012 OPEN CASCADE SAS
3//
4// The content of this file is subject to the Open CASCADE Technology Public
5// License Version 6.5 (the "License"). You may not use the content of this file
6// except in compliance with the License. Please obtain a copy of the License
7// at http://www.opencascade.org and read it completely before using this file.
8//
9// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
10// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
11//
12// The Original Code and all software distributed under the License is
13// distributed on an "AS IS" basis, without warranty of any kind, and the
14// Initial Developer hereby disclaims all such warranties, including without
15// limitation, any warranties of merchantability, fitness for a particular
16// purpose or non-infringement. Please see the License for the specific terms
17// and conditions governing the rights and limitations under the License.
18
7fd59977 19#include <Standard_NotImplemented.hxx>
20#include <math_Vector.hxx>
21#include <math.hxx>
22#include <gp_Pnt2d.hxx>
23#include <gp_Vec2d.hxx>
24#include <gp_Pnt.hxx>
25#include <gp_Vec.hxx>
26
27#include <TColStd_Array1OfReal.hxx>
28#include <Precision.hxx>
29
30class HMath_Vector{
31 math_Vector *pvec;
32 void operator=(const math_Vector&){}
33 public:
34 HMath_Vector(){ pvec = 0;}
35 HMath_Vector(math_Vector* pv){ pvec = pv;}
36 ~HMath_Vector(){ if(pvec != 0) delete pvec;}
37 void operator=(math_Vector* pv){ if(pvec != pv && pvec != 0) delete pvec; pvec = pv;}
38 Standard_Real& operator()(Standard_Integer i){ return (*pvec).operator()(i);}
39 const Standard_Real& operator()(Standard_Integer i) const{ return (*pvec).operator()(i);}
40 const math_Vector* operator->() const{ return pvec;}
41 math_Vector* operator->(){ return pvec;}
42 math_Vector* Vector(){ return pvec;}
43 math_Vector* Init(Standard_Real v, Standard_Integer i = 0, Standard_Integer iEnd = 0){
44 if(pvec == 0) return pvec;
45 if(iEnd - i == 0) pvec->Init(v);
46 else { Standard_Integer End = (iEnd <= pvec->Upper()) ? iEnd : pvec->Upper();
47 for(; i <= End; i++) pvec->operator()(i) = v; }
48 return pvec;
49 }
50};
51
52static Standard_Real EPS_PARAM = 1.e-12;
53static Standard_Real EPS_DIM = 1.e-20;
54static Standard_Real ERROR_ALGEBR_RATIO = 2.0/3.0;
55
56static Standard_Integer GPM = 61;
57static Standard_Integer SUBS_POWER = 32;
58static Standard_Integer SM = 1953;
59
60static math_Vector LGaussP0(1,GPM);
61static math_Vector LGaussW0(1,GPM);
62static math_Vector LGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
63static math_Vector LGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
64
65static math_Vector* LGaussP[] = {&LGaussP0,&LGaussP1};
66static math_Vector* LGaussW[] = {&LGaussW0,&LGaussW1};
67
68static HMath_Vector L1 = new math_Vector(1,SM,0.0);
69static HMath_Vector L2 = new math_Vector(1,SM,0.0);
70static HMath_Vector DimL = new math_Vector(1,SM,0.0);
71static HMath_Vector ErrL = new math_Vector(1,SM,0.0);
72static HMath_Vector ErrUL = new math_Vector(1,SM,0.0);
73static HMath_Vector IxL = new math_Vector(1,SM,0.0);
74static HMath_Vector IyL = new math_Vector(1,SM,0.0);
75static HMath_Vector IzL = new math_Vector(1,SM,0.0);
76static HMath_Vector IxxL = new math_Vector(1,SM,0.0);
77static HMath_Vector IyyL = new math_Vector(1,SM,0.0);
78static HMath_Vector IzzL = new math_Vector(1,SM,0.0);
79static HMath_Vector IxyL = new math_Vector(1,SM,0.0);
80static HMath_Vector IxzL = new math_Vector(1,SM,0.0);
81static HMath_Vector IyzL = new math_Vector(1,SM,0.0);
82
83static math_Vector UGaussP0(1,GPM);
84static math_Vector UGaussW0(1,GPM);
85static math_Vector UGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
86static math_Vector UGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
87
88static math_Vector* UGaussP[] = {&UGaussP0,&UGaussP1};
89static math_Vector* UGaussW[] = {&UGaussW0,&UGaussW1};
90
91static HMath_Vector U1 = new math_Vector(1,SM,0.0);
92static HMath_Vector U2 = new math_Vector(1,SM,0.0);
93static HMath_Vector DimU = new math_Vector(1,SM,0.0);
94static HMath_Vector ErrU = new math_Vector(1,SM,0.0);
95static HMath_Vector IxU = new math_Vector(1,SM,0.0);
96static HMath_Vector IyU = new math_Vector(1,SM,0.0);
97static HMath_Vector IzU = new math_Vector(1,SM,0.0);
98static HMath_Vector IxxU = new math_Vector(1,SM,0.0);
99static HMath_Vector IyyU = new math_Vector(1,SM,0.0);
100static HMath_Vector IzzU = new math_Vector(1,SM,0.0);
101static HMath_Vector IxyU = new math_Vector(1,SM,0.0);
102static HMath_Vector IxzU = new math_Vector(1,SM,0.0);
103static HMath_Vector IyzU = new math_Vector(1,SM,0.0);
104
105static Standard_Integer FillIntervalBounds(Standard_Real A,
106 Standard_Real B,
107 const TColStd_Array1OfReal& Knots,
108 HMath_Vector& VA,
109 HMath_Vector& VB)
110{
111 Standard_Integer i = 1, iEnd = Knots.Upper(), j = 1, k = 1;
112 VA(j++) = A;
113 for(; i <= iEnd; i++){
114 Standard_Real kn = Knots(i);
115 if(A < kn)
116 if(kn < B) VA(j++) = VB(k++) = kn; else break;
117 }
118 VB(k) = B;
119 return k;
120}
121
122static inline Standard_Integer MaxSubs(Standard_Integer n, Standard_Integer coeff = SUBS_POWER){
123// return n = IntegerLast()/coeff < n? IntegerLast(): n*coeff + 1;
124 return Min((n * coeff + 1),SM);
125}
126
127static Standard_Integer LFillIntervalBounds(Standard_Real A,
128 Standard_Real B,
129 const TColStd_Array1OfReal& Knots,
130 const Standard_Integer NumSubs)
131{
132 Standard_Integer iEnd = Knots.Upper(), jEnd = L1->Upper();
133 if(iEnd - 1 > jEnd){
134 iEnd = MaxSubs(iEnd-1,NumSubs);
135 L1 = new math_Vector(1,iEnd);
136 L2 = new math_Vector(1,iEnd);
137 DimL = new math_Vector(1,iEnd);
138 ErrL = new math_Vector(1,iEnd,0.0);
139 ErrUL = new math_Vector(1,iEnd,0.0);
140 IxL = new math_Vector(1,iEnd);
141 IyL = new math_Vector(1,iEnd);
142 IzL = new math_Vector(1,iEnd);
143 IxxL = new math_Vector(1,iEnd);
144 IyyL = new math_Vector(1,iEnd);
145 IzzL = new math_Vector(1,iEnd);
146 IxyL = new math_Vector(1,iEnd);
147 IxzL = new math_Vector(1,iEnd);
148 IyzL = new math_Vector(1,iEnd);
149 }
150 return FillIntervalBounds(A, B, Knots, L1, L2);
151}
152
153static Standard_Integer UFillIntervalBounds(Standard_Real A,
154 Standard_Real B,
155 const TColStd_Array1OfReal& Knots,
156 const Standard_Integer NumSubs)
157{
158 Standard_Integer iEnd = Knots.Upper(), jEnd = U1->Upper();
159 if(iEnd - 1 > jEnd){
160 iEnd = MaxSubs(iEnd-1,NumSubs);
161 U1 = new math_Vector(1,iEnd);
162 U2 = new math_Vector(1,iEnd);
163 DimU = new math_Vector(1,iEnd);
164 ErrU = new math_Vector(1,iEnd,0.0);
165 IxU = new math_Vector(1,iEnd);
166 IyU = new math_Vector(1,iEnd);
167 IzU = new math_Vector(1,iEnd);
168 IxxU = new math_Vector(1,iEnd);
169 IyyU = new math_Vector(1,iEnd);
170 IzzU = new math_Vector(1,iEnd);
171 IxyU = new math_Vector(1,iEnd);
172 IxzU = new math_Vector(1,iEnd);
173 IyzU = new math_Vector(1,iEnd);
174 }
175 return FillIntervalBounds(A, B, Knots, U1, U2);
176}
177
178static Standard_Real CCompute(Face& S,
179 Domain& D,
180 const gp_Pnt& loc,
181 Standard_Real& Dim,
182 gp_Pnt& g,
183 gp_Mat& inertia,
184 const Standard_Real EpsDim,
185 const Standard_Boolean isErrorCalculation,
186 const Standard_Boolean isVerifyComputation)
187{
188 Standard_Boolean isNaturalRestriction = S.NaturalRestriction();
189
190 Standard_Integer NumSubs = SUBS_POWER;
191
192 Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
193 Dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
194 Standard_Real x, y, z;
195 //boundary curve parametrization
196 Standard_Real l1, l2, lm, lr, l;
197 //Face parametrization in U and V direction
198 Standard_Real BV1, BV2, v;
199 Standard_Real BU1, BU2, u1, u2, um, ur, u;
200 S.Bounds (BU1, BU2, BV1, BV2); u1 = BU1;
201 //location point used to compute the inertia
202 Standard_Real xloc, yloc, zloc;
203 loc.Coord (xloc, yloc, zloc); // use member of parent class
204 //Jacobien (x, y, z) -> (u, v) = ||n||
205 Standard_Real ds;
206 //On the Face
207 gp_Pnt Ps;
208 gp_Vec VNor;
209 //On the boundary curve u-v
210 gp_Pnt2d Puv;
211 gp_Vec2d Vuv;
212 Standard_Real Dul; // Dul = Du / Dl
213 Standard_Real CDim[2], CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz;
214 Standard_Real LocDim[2], LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy, LocIxz, LocIyz;
215
216 Standard_Real ErrorU, ErrorL, ErrorLMax = 0.0, Eps=0.0, EpsL=0.0, EpsU=0.0;
217
218 Standard_Integer iD = 0, NbLSubs, iLS, iLSubEnd, iGL, iGLEnd, NbLGaussP[2], LRange[2], iL, kL, kLEnd, IL, JL;
219 Standard_Integer i, NbUSubs, iUS, iUSubEnd, iGU, iGUEnd, NbUGaussP[2], URange[2], iU, kU, kUEnd, IU, JU;
220 Standard_Integer UMaxSubs, LMaxSubs;
221 iGLEnd = isErrorCalculation? 2: 1;
222 for(i = 0; i < 2; i++) {
223 LocDim[i] = 0.0;
224 CDim[i] = 0.0;
225 }
226
227 NbUGaussP[0] = S.SIntOrder(EpsDim);
228 NbUGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbUGaussP[0]));
229 math::GaussPoints(NbUGaussP[0],UGaussP0); math::GaussWeights(NbUGaussP[0],UGaussW0);
230 math::GaussPoints(NbUGaussP[1],UGaussP1); math::GaussWeights(NbUGaussP[1],UGaussW1);
231
232 NbUSubs = S.SUIntSubs();
233 TColStd_Array1OfReal UKnots(1,NbUSubs+1);
234 S.UKnots(UKnots);
235
236
237 while (isNaturalRestriction || D.More()) {
238 if(isNaturalRestriction){
239 NbLGaussP[0] = Min(2*NbUGaussP[0],math::GaussPointsMax());
240 }else{
241 S.Load(D.Value()); ++iD;
242 NbLGaussP[0] = S.LIntOrder(EpsDim);
243 }
244
245
246 NbLGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbLGaussP[0]));
247 math::GaussPoints(NbLGaussP[0],LGaussP0); math::GaussWeights(NbLGaussP[0],LGaussW0);
248 math::GaussPoints(NbLGaussP[1],LGaussP1); math::GaussWeights(NbLGaussP[1],LGaussW1);
249
250 NbLSubs = isNaturalRestriction? S.SVIntSubs(): S.LIntSubs();
251
252 TColStd_Array1OfReal LKnots(1,NbLSubs+1);
253 if(isNaturalRestriction){
254 S.VKnots(LKnots);
255 l1 = BV1; l2 = BV2;
256 }else{
257 S.LKnots(LKnots);
258 l1 = S.FirstParameter(); l2 = S.LastParameter();
259 }
260 ErrorL = 0.0;
261 kLEnd = 1; JL = 0;
262 //OCC503(apo): if(Abs(l2-l1) < EPS_PARAM) continue;
263 if(Abs(l2-l1) > EPS_PARAM) {
264 iLSubEnd = LFillIntervalBounds(l1, l2, LKnots, NumSubs);
265 LMaxSubs = MaxSubs(iLSubEnd);
266 if(LMaxSubs > DimL.Vector()->Upper()) LMaxSubs = DimL.Vector()->Upper();
267 DimL.Init(0.0,1,LMaxSubs); ErrL.Init(0.0,1,LMaxSubs); ErrUL.Init(0.0,1,LMaxSubs);
268 do{// while: L
269 if(++JL > iLSubEnd){
270 LRange[0] = IL = ErrL->Max(); LRange[1] = JL;
271 L1(JL) = (L1(IL) + L2(IL))/2.0; L2(JL) = L2(IL); L2(IL) = L1(JL);
272 }else LRange[0] = IL = JL;
273 if(JL == LMaxSubs || Abs(L2(JL) - L1(JL)) < EPS_PARAM)
274 if(kLEnd == 1){
275 DimL(JL) = ErrL(JL) = IxL(JL) = IyL(JL) = IzL(JL) =
276 IxxL(JL) = IyyL(JL) = IzzL(JL) = IxyL(JL) = IxzL(JL) = IyzL(JL) = 0.0;
277 }else{
278 JL--;
279 EpsL = ErrorL; Eps = EpsL/0.9;
280 break;
281 }
282 else
283 for(kL=0; kL < kLEnd; kL++){
284 iLS = LRange[kL];
285 lm = 0.5*(L2(iLS) + L1(iLS));
286 lr = 0.5*(L2(iLS) - L1(iLS));
287 CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
288 for(iGL=0; iGL < iGLEnd; iGL++){//
289 CDim[iGL] = 0.0;
290 for(iL=1; iL<=NbLGaussP[iGL]; iL++){
291 l = lm + lr*(*LGaussP[iGL])(iL);
292 if(isNaturalRestriction){
293 v = l; u2 = BU2; Dul = (*LGaussW[iGL])(iL);
294 }else{
295 S.D12d (l, Puv, Vuv);
296 Dul = Vuv.Y()*(*LGaussW[iGL])(iL); // Dul = Du / Dl
297 if(Abs(Dul) < EPS_PARAM) continue;
298 v = Puv.Y(); u2 = Puv.X();
299 //Check on cause out off bounds of value current parameter
300 if(v < BV1) v = BV1; else if(v > BV2) v = BV2;
301 if(u2 < BU1) u2 = BU1; else if(u2 > BU2) u2 = BU2;
302 }
303 ErrUL(iLS) = 0.0;
304 kUEnd = 1; JU = 0;
305 if(Abs(u2-u1) < EPS_PARAM) continue;
306 iUSubEnd = UFillIntervalBounds(u1, u2, UKnots, NumSubs);
307 UMaxSubs = MaxSubs(iUSubEnd);
308 if(UMaxSubs > DimU.Vector()->Upper()) UMaxSubs = DimU.Vector()->Upper();
309 DimU.Init(0.0,1,UMaxSubs); ErrU.Init(0.0,1,UMaxSubs); ErrorU = 0.0;
310 do{//while: U
311 if(++JU > iUSubEnd){
312 URange[0] = IU = ErrU->Max(); URange[1] = JU;
313 U1(JU) = (U1(IU)+U2(IU))/2.0; U2(JU) = U2(IU); U2(IU) = U1(JU);
314 }else URange[0] = IU = JU;
315 if(JU == UMaxSubs || Abs(U2(JU) - U1(JU)) < EPS_PARAM)
316 if(kUEnd == 1){
317 DimU(JU) = ErrU(JU) = IxU(JU) = IyU(JU) = IzU(JU) =
318 IxxU(JU) = IyyU(JU) = IzzU(JU) = IxyU(JU) = IxzU(JU) = IyzU(JU) = 0.0;
319 }else{
320 JU--;
321 EpsU = ErrorU; Eps = EpsU*Abs((u2-u1)*Dul)/0.1; EpsL = 0.9*Eps;
322 break;
323 }
324 else
325 for(kU=0; kU < kUEnd; kU++){
326 iUS = URange[kU];
327 um = 0.5*(U2(iUS) + U1(iUS));
328 ur = 0.5*(U2(iUS) - U1(iUS));
329 LocIx = LocIy = LocIz = LocIxx = LocIyy = LocIzz = LocIxy = LocIxz = LocIyz = 0.0;
330 iGUEnd = iGLEnd - iGL;
331 for(iGU=0; iGU < iGUEnd; iGU++){//
332 LocDim[iGU] = 0.0;
333 for(iU=1; iU<=NbUGaussP[iGU]; iU++){
334 u = um + ur*(*UGaussP[iGU])(iU);
335 S.Normal(u, v, Ps, VNor);
336 ds = VNor.Magnitude(); //Jacobien(x,y,z) -> (u,v)=||n||
337 ds *= (*UGaussW[iGU])(iU);
338 LocDim[iGU] += ds;
339 if(iGU > 0) continue;
340 Ps.Coord(x, y, z);
341 x -= xloc; y -= yloc; z -= zloc;
342 LocIx += x*ds; LocIy += y*ds; LocIz += z*ds;
343 LocIxy += x*y*ds; LocIyz += y*z*ds; LocIxz += x*z*ds;
344 x *= x; y *= y; z *= z;
345 LocIxx += (y+z)*ds; LocIyy += (x+z)*ds; LocIzz += (x+y)*ds;
346 }//for: iU
347 }//for: iGU
348 DimU(iUS) = LocDim[0]*ur;
349 if(iGL > 0) continue;
350 ErrU(iUS) = Abs(LocDim[1]-LocDim[0])*ur;
351 IxU(iUS) = LocIx*ur; IyU(iUS) = LocIy*ur; IzU(iUS) = LocIz*ur;
352 IxxU(iUS) = LocIxx*ur; IyyU(iUS) = LocIyy*ur; IzzU(iUS) = LocIzz*ur;
353 IxyU(iUS) = LocIxy*ur; IxzU(iUS) = LocIxz*ur; IyzU(iUS) = LocIyz*ur;
354 }//for: kU (iUS)
355 if(JU == iUSubEnd) kUEnd = 2;
356 if(kUEnd == 2) ErrorU = ErrU(ErrU->Max());
357 }while((ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1);
358 for(i=1; i<=JU; i++) CDim[iGL] += DimU(i)*Dul;
359 if(iGL > 0) continue;
360 ErrUL(iLS) = ErrorU*Abs((u2-u1)*Dul);
361 for(i=1; i<=JU; i++){
362 CIx += IxU(i)*Dul; CIy += IyU(i)*Dul; CIz += IzU(i)*Dul;
363 CIxx += IxxU(i)*Dul; CIyy += IyyU(i)*Dul; CIzz += IzzU(i)*Dul;
364 CIxy += IxyU(i)*Dul; CIxz += IxzU(i)*Dul; CIyz += IyzU(i)*Dul;
365 }
366 }//for: iL
367 }//for: iGL
368 DimL(iLS) = CDim[0]*lr;
369 if(iGLEnd == 2) ErrL(iLS) = Abs(CDim[1]-CDim[0])*lr + ErrUL(iLS);
370 IxL(iLS) = CIx*lr; IyL(iLS) = CIy*lr; IzL(iLS) = CIz*lr;
371 IxxL(iLS) = CIxx*lr; IyyL(iLS) = CIyy*lr; IzzL(iLS) = CIzz*lr;
372 IxyL(iLS) = CIxy*lr; IxzL(iLS) = CIxz*lr; IyzL(iLS) = CIyz*lr;
373 }//for: (kL)iLS
374 // Calculate/correct epsilon of computation by current value of Dim
375 //That is need for not spend time for
376 if(JL == iLSubEnd){
377 kLEnd = 2;
378 Standard_Real DDim = 0.0;
379 for(i=1; i<=JL; i++) DDim += DimL(i);
380 DDim = Abs(DDim*EpsDim);
381 if(DDim > Eps) {
382 Eps = DDim; EpsL = 0.9*Eps;
383 }
384 }
385 if(kLEnd == 2) ErrorL = ErrL(ErrL->Max());
386 }while((ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1);
387 for(i=1; i<=JL; i++){
388 Dim += DimL(i);
389 Ix += IxL(i); Iy += IyL(i); Iz += IzL(i);
390 Ixx += IxxL(i); Iyy += IyyL(i); Izz += IzzL(i);
391 Ixy += IxyL(i); Ixz += IxzL(i); Iyz += IyzL(i);
392 }
393 ErrorLMax = Max(ErrorLMax, ErrorL);
394 }
395 if(isNaturalRestriction) break;
396 D.Next();
397 }
398 if(Abs(Dim) >= EPS_DIM){
399 Ix /= Dim; Iy /= Dim; Iz /= Dim;
400 g.SetCoord (Ix, Iy, Iz);
401 }else{
402 Dim =0.0;
403 g.SetCoord (0., 0.,0.);
404 }
405 inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
406 gp_XYZ (-Ixy, Iyy, -Iyz),
407 gp_XYZ (-Ixz, -Iyz, Izz));
408
409 if(iGLEnd == 2) Eps = Dim != 0.0? ErrorLMax/Abs(Dim): 0.0;
410 else Eps = EpsDim;
411 return Eps;
412}
413
414static Standard_Real Compute(Face& S, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia,
415 Standard_Real EpsDim)
416{
417 Standard_Boolean isErrorCalculation = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
418 Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
419 EpsDim = Abs(EpsDim);
420 Domain D;
421 return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
422}
423
424static Standard_Real Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia,
425 Standard_Real EpsDim)
426{
427 Standard_Boolean isErrorCalculation = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
428 Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
429 EpsDim = Abs(EpsDim);
430 return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
431}
432
433static void Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& dim, gp_Pnt& g, gp_Mat& inertia){
434 Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
435 dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
436
437 Standard_Real x, y, z;
438 Standard_Integer NbCGaussgp_Pnts = 0;
439
440 Standard_Real l1, l2, lm, lr, l; //boundary curve parametrization
441 Standard_Real v1, v2, vm, vr, v; //Face parametrization in v direction
442 Standard_Real u1, u2, um, ur, u;
443 Standard_Real ds; //Jacobien (x, y, z) -> (u, v) = ||n||
444
445 gp_Pnt P; //On the Face
446 gp_Vec VNor;
447
448 gp_Pnt2d Puv; //On the boundary curve u-v
449 gp_Vec2d Vuv;
450 Standard_Real Dul; // Dul = Du / Dl
451 Standard_Real CArea, CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz;
452 Standard_Real LocArea, LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy,
453 LocIxz, LocIyz;
454
455
456 S.Bounds (u1, u2, v1, v2);
457
458 Standard_Integer NbUGaussgp_Pnts = Min(S.UIntegrationOrder (),
459 math::GaussPointsMax());
460 Standard_Integer NbVGaussgp_Pnts = Min(S.VIntegrationOrder (),
461 math::GaussPointsMax());
462
463 Standard_Integer NbGaussgp_Pnts = Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts);
464
465 //Number of Gauss points for the integration
466 //on the Face
467 math_Vector GaussSPV (1, NbGaussgp_Pnts);
468 math_Vector GaussSWV (1, NbGaussgp_Pnts);
469 math::GaussPoints (NbGaussgp_Pnts,GaussSPV);
470 math::GaussWeights (NbGaussgp_Pnts,GaussSWV);
471
472
473 //location point used to compute the inertia
474 Standard_Real xloc, yloc, zloc;
475 loc.Coord (xloc, yloc, zloc);
476
477 while (D.More()) {
478
479 S.Load(D.Value());
480 NbCGaussgp_Pnts = Min(S.IntegrationOrder (), math::GaussPointsMax());
481
482 math_Vector GaussCP (1, NbCGaussgp_Pnts);
483 math_Vector GaussCW (1, NbCGaussgp_Pnts);
484 math::GaussPoints (NbCGaussgp_Pnts,GaussCP);
485 math::GaussWeights (NbCGaussgp_Pnts,GaussCW);
486
487 CArea = 0.0;
488 CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
489 l1 = S.FirstParameter ();
490 l2 = S.LastParameter ();
491 lm = 0.5 * (l2 + l1);
492 lr = 0.5 * (l2 - l1);
493
494 Puv = S.Value2d (lm);
495 vm = Puv.Y();
496 Puv = S.Value2d (lr);
497 vr = Puv.Y();
498
499 for (Standard_Integer i = 1; i <= NbCGaussgp_Pnts; i++) {
500 l = lm + lr * GaussCP (i);
501 S.D12d(l, Puv, Vuv);
502 v = Puv.Y();
503 u2 = Puv.X();
504 Dul = Vuv.Y();
505 Dul *= GaussCW (i);
506 um = 0.5 * (u2 + u1);
507 ur = 0.5 * (u2 - u1);
508 LocArea = LocIx = LocIy = LocIz = LocIxx = LocIyy = LocIzz =
509 LocIxy = LocIxz = LocIyz = 0.0;
510 for (Standard_Integer j = 1; j <= NbGaussgp_Pnts; j++) {
511 u = um + ur * GaussSPV (j);
512 S.Normal (u, v, P, VNor);
513 ds = VNor.Magnitude(); //normal.Magnitude
514 ds = ds * Dul * GaussSWV (j);
515 LocArea += ds;
516 P.Coord (x, y, z);
517 x -= xloc;
518 y -= yloc;
519 z -= zloc;
520 LocIx += x * ds;
521 LocIy += y * ds;
522 LocIz += z * ds;
523 LocIxy += x * y * ds;
524 LocIyz += y * z * ds;
525 LocIxz += x * z * ds;
526 x *= x;
527 y *= y;
528 z *= z;
529 LocIxx += (y + z) * ds;
530 LocIyy += (x + z) * ds;
531 LocIzz += (x + y) * ds;
532 }
533 CArea += LocArea * ur;
534 CIx += LocIx * ur;
535 CIy += LocIy * ur;
536 CIz += LocIz * ur;
537 CIxx += LocIxx * ur;
538 CIyy += LocIyy * ur;
539 CIzz += LocIzz * ur;
540 CIxy += LocIxy * ur;
541 CIxz += LocIxz * ur;
542 CIyz += LocIyz * ur;
543 }
544 dim += CArea * lr;
545 Ix += CIx * lr;
546 Iy += CIy * lr;
547 Iz += CIz * lr;
548 Ixx += CIxx * lr;
549 Iyy += CIyy * lr;
550 Izz += CIzz * lr;
551 Ixy += CIxy * lr;
552 Ixz += CIxz * lr;
553 Iyz += CIyz * lr;
554 D.Next();
555 }
556 if (Abs(dim) >= EPS_DIM) {
557 Ix /= dim;
558 Iy /= dim;
559 Iz /= dim;
560 g.SetCoord (Ix, Iy, Iz);
561 }
562 else {
563 dim =0.;
564 g.SetCoord (0., 0.,0.);
565 }
566 inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
567 gp_XYZ (-Ixy, Iyy, -Iyz),
568 gp_XYZ (-Ixz, -Iyz, Izz));
569}
570
571
572
573static void Compute(const Face& S, const gp_Pnt& loc, Standard_Real& dim, gp_Pnt& g, gp_Mat& inertia){
574 Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
575 dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
576
577 Standard_Real LowerU, UpperU, LowerV, UpperV;
578 S.Bounds (LowerU, UpperU, LowerV, UpperV);
579 Standard_Integer UOrder = Min(S.UIntegrationOrder (),
580 math::GaussPointsMax());
581 Standard_Integer VOrder = Min(S.VIntegrationOrder (),
582 math::GaussPointsMax());
583 gp_Pnt P;
584 gp_Vec VNor;
585 Standard_Real dsi, ds;
586 Standard_Real ur, um, u, vr, vm, v;
587 Standard_Real x, y, z;
588 Standard_Real Ixi, Iyi, Izi, Ixxi, Iyyi, Izzi, Ixyi, Ixzi, Iyzi;
589 Standard_Real xloc, yloc, zloc;
590 loc.Coord (xloc, yloc, zloc);
591
592 Standard_Integer i, j;
593 math_Vector GaussPU (1, UOrder); //gauss points and weights
594 math_Vector GaussWU (1, UOrder);
595 math_Vector GaussPV (1, VOrder);
596 math_Vector GaussWV (1, VOrder);
597
598 //Recuperation des points de Gauss dans le fichier GaussPoints.
599 math::GaussPoints (UOrder,GaussPU);
600 math::GaussWeights (UOrder,GaussWU);
601 math::GaussPoints (VOrder,GaussPV);
602 math::GaussWeights (VOrder,GaussWV);
603
604 // Calcul des integrales aux points de gauss :
605 um = 0.5 * (UpperU + LowerU);
606 vm = 0.5 * (UpperV + LowerV);
607 ur = 0.5 * (UpperU - LowerU);
608 vr = 0.5 * (UpperV - LowerV);
609
610 for (j = 1; j <= VOrder; j++) {
611 v = vm + vr * GaussPV (j);
612 dsi = Ixi = Iyi = Izi = Ixxi = Iyyi = Izzi = Ixyi = Ixzi = Iyzi = 0.0;
613
614 for (i = 1; i <= UOrder; i++) {
615 u = um + ur * GaussPU (i);
616 S.Normal (u, v, P, VNor);
617 ds = VNor.Magnitude() * GaussWU (i);
618 P.Coord (x, y, z);
619 x -= xloc;
620 y -= yloc;
621 z -= zloc;
622 dsi += ds;
623 Ixi += x * ds;
624 Iyi += y * ds;
625 Izi += z * ds;
626 Ixyi += x * y * ds;
627 Iyzi += y * z * ds;
628 Ixzi += x * z * ds;
629 x *= x;
630 y *= y;
631 z *= z;
632 Ixxi += (y + z) * ds;
633 Iyyi += (x + z) * ds;
634 Izzi += (x + y) * ds;
635 }
636 dim += dsi * GaussWV (j);
637 Ix += Ixi * GaussWV (j);
638 Iy += Iyi * GaussWV (j);
639 Iz += Izi * GaussWV (j);
640 Ixx += Ixxi * GaussWV (j);
641 Iyy += Iyyi * GaussWV (j);
642 Izz += Izzi * GaussWV (j);
643 Ixy += Ixyi * GaussWV (j);
644 Iyz += Iyzi * GaussWV (j);
645 Ixz += Ixzi * GaussWV (j);
646 }
647 vr *= ur;
648 Ixx *= vr;
649 Iyy *= vr;
650 Izz *= vr;
651 Ixy *= vr;
652 Ixz *= vr;
653 Iyz *= vr;
654 if (Abs(dim) >= EPS_DIM) {
655 Ix /= dim;
656 Iy /= dim;
657 Iz /= dim;
658 dim *= vr;
659 g.SetCoord (Ix, Iy, Iz);
660 }
661 else {
662 dim =0.;
663 g.SetCoord (0.,0.,0.);
664 }
665 inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
666 gp_XYZ (-Ixy, Iyy, -Iyz),
667 gp_XYZ (-Ixz, -Iyz, Izz));
668}
669
670GProp_SGProps::GProp_SGProps(){}
671
672GProp_SGProps::GProp_SGProps (const Face& S,
673 const gp_Pnt& SLocation
674 )
675{
676 SetLocation(SLocation);
677 Perform(S);
678}
679
680GProp_SGProps::GProp_SGProps (Face& S,
681 Domain& D,
682 const gp_Pnt& SLocation
683 )
684{
685 SetLocation(SLocation);
686 Perform(S,D);
687}
688
689GProp_SGProps::GProp_SGProps(Face& S, const gp_Pnt& SLocation, const Standard_Real Eps){
690 SetLocation(SLocation);
691 Perform(S, Eps);
692}
693
694GProp_SGProps::GProp_SGProps(Face& S, Domain& D, const gp_Pnt& SLocation, const Standard_Real Eps){
695 SetLocation(SLocation);
696 Perform(S, D, Eps);
697}
698
699void GProp_SGProps::SetLocation(const gp_Pnt& SLocation){
700 loc = SLocation;
701}
702
703void GProp_SGProps::Perform(const Face& S){
704 Compute(S,loc,dim,g,inertia);
705 myEpsilon = 1.0;
706 return;
707}
708
709void GProp_SGProps::Perform(Face& S, Domain& D){
710 Compute(S,D,loc,dim,g,inertia);
711 myEpsilon = 1.0;
712 return;
713}
714
715Standard_Real GProp_SGProps::Perform(Face& S, const Standard_Real Eps){
716 return myEpsilon = Compute(S,loc,dim,g,inertia,Eps);
717}
718
719Standard_Real GProp_SGProps::Perform(Face& S, Domain& D, const Standard_Real Eps){
720 return myEpsilon = Compute(S,D,loc,dim,g,inertia,Eps);
721}
722
723
724Standard_Real GProp_SGProps::GetEpsilon(){
725 return myEpsilon;
726}
Open CASCADE Technology repository