1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
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.
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.
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.
19 #include <Standard_NotImplemented.hxx>
20 #include <math_Vector.hxx>
22 #include <gp_Pnt2d.hxx>
23 #include <gp_Vec2d.hxx>
27 #include <TColStd_Array1OfReal.hxx>
28 #include <Precision.hxx>
31 void operator=(const math_Vector&){}
33 HMath_Vector(){ pvec = 0;}
34 HMath_Vector(math_Vector* pv){ pvec = pv;}
35 ~HMath_Vector(){ if(pvec != 0) delete pvec;}
36 void operator=(math_Vector* pv){ if(pvec != pv && pvec != 0) delete pvec; pvec = pv;}
37 Standard_Real& operator()(Standard_Integer i){ return (*pvec).operator()(i);}
38 const Standard_Real& operator()(Standard_Integer i) const{ return (*pvec).operator()(i);}
39 const math_Vector* operator->() const{ return pvec;}
40 math_Vector* operator->(){ return pvec;}
41 math_Vector* Init(Standard_Real v, Standard_Integer i = 0, Standard_Integer iEnd = 0){
42 if(pvec == 0) return pvec;
43 if(iEnd - i == 0) pvec->Init(v);
44 else for(; i <= iEnd; i++) pvec->operator()(i) = v;
49 //Minimal value of interval's range for computation | minimal value of "dim" | ...
50 static Standard_Real EPS_PARAM = Precision::Angular(), EPS_DIM = 1.E-30, ERROR_ALGEBR_RATIO = 2.0/3.0;
51 //Maximum of GaussPoints on a subinterval and maximum of subintervals
52 static Standard_Integer GPM = math::GaussPointsMax(), SUBS_POWER = 32, SM = SUBS_POWER*GPM + 1;
53 static Standard_Boolean IS_MIN_DIM = 1; // if the value equal 0 error of algorithm calculted by static moments
55 static math_Vector LGaussP0(1,GPM), LGaussW0(1,GPM),
56 LGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM))), LGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
57 static HMath_Vector L1 = new math_Vector(1,SM), L2 = new math_Vector(1,SM),
58 DimL = new math_Vector(1,SM), ErrL = new math_Vector(1,SM), ErrUL = new math_Vector(1,SM,0.0),
59 IxL = new math_Vector(1,SM), IyL = new math_Vector(1,SM), IzL = new math_Vector(1,SM),
60 IxxL = new math_Vector(1,SM), IyyL = new math_Vector(1,SM), IzzL = new math_Vector(1,SM),
61 IxyL = new math_Vector(1,SM), IxzL = new math_Vector(1,SM), IyzL = new math_Vector(1,SM);
63 static math_Vector* LGaussP[] = {&LGaussP0,&LGaussP1};
64 static math_Vector* LGaussW[] = {&LGaussW0,&LGaussW1};
66 static math_Vector UGaussP0(1,GPM), UGaussW0(1,GPM),
67 UGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM))), UGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
68 static HMath_Vector U1 = new math_Vector(1,SM), U2 = new math_Vector(1,SM),
69 DimU = new math_Vector(1,SM), ErrU = new math_Vector(1,SM,0.0),
70 IxU = new math_Vector(1,SM), IyU = new math_Vector(1,SM), IzU = new math_Vector(1,SM),
71 IxxU = new math_Vector(1,SM), IyyU = new math_Vector(1,SM), IzzU = new math_Vector(1,SM),
72 IxyU = new math_Vector(1,SM), IxzU = new math_Vector(1,SM), IyzU = new math_Vector(1,SM);
74 static math_Vector* UGaussP[] = {&UGaussP0,&UGaussP1};
75 static math_Vector* UGaussW[] = {&UGaussW0,&UGaussW1};
77 static Standard_Integer FillIntervalBounds(Standard_Real A, Standard_Real B, const TColStd_Array1OfReal& Knots,
78 HMath_Vector& VA, HMath_Vector& VB)
80 Standard_Integer i = 1, iEnd = Knots.Upper(), j = 1, k = 1;
82 for(; i <= iEnd; i++){
83 Standard_Real kn = Knots(i);
88 VA(j++) = VB(k++) = kn;
100 static inline Standard_Integer MaxSubs(Standard_Integer n, Standard_Integer coeff = SUBS_POWER){
101 return n = IntegerLast()/coeff < n? IntegerLast(): n*coeff + 1;
104 static Standard_Integer LFillIntervalBounds(Standard_Real A, Standard_Real B, const TColStd_Array1OfReal& Knots,
105 const Standard_Integer NumSubs)
107 Standard_Integer iEnd = Knots.Upper(), jEnd = L1->Upper();
109 // Modified by Sergey KHROMOV - Wed Mar 26 11:22:50 2003
110 iEnd = Max(iEnd, MaxSubs(iEnd-1,NumSubs));
112 // iEnd = MaxSubs(iEnd-1,NumSubs);
113 // Modified by Sergey KHROMOV - Wed Mar 26 11:22:51 2003
114 L1 = new math_Vector(1,iEnd); L2 = new math_Vector(1,iEnd);
115 DimL = new math_Vector(1,iEnd); ErrL = new math_Vector(1,iEnd,0.0); ErrUL = new math_Vector(1,iEnd,0.0);
116 IxL = new math_Vector(1,iEnd); IyL = new math_Vector(1,iEnd); IzL = new math_Vector(1,iEnd);
117 IxxL = new math_Vector(1,iEnd); IyyL = new math_Vector(1,iEnd); IzzL = new math_Vector(1,iEnd);
118 IxyL = new math_Vector(1,iEnd); IxzL = new math_Vector(1,iEnd); IyzL = new math_Vector(1,iEnd);
120 return FillIntervalBounds(A, B, Knots, L1, L2);
123 static Standard_Integer UFillIntervalBounds(Standard_Real A, Standard_Real B, const TColStd_Array1OfReal& Knots,
124 const Standard_Integer NumSubs)
126 Standard_Integer iEnd = Knots.Upper(), jEnd = U1->Upper();
128 // Modified by Sergey KHROMOV - Wed Mar 26 11:22:50 2003
129 iEnd = Max(iEnd, MaxSubs(iEnd-1,NumSubs));
131 // iEnd = MaxSubs(iEnd-1,NumSubs);
132 // Modified by Sergey KHROMOV - Wed Mar 26 11:22:51 2003
133 U1 = new math_Vector(1,iEnd); U2 = new math_Vector(1,iEnd);
134 DimU = new math_Vector(1,iEnd); ErrU = new math_Vector(1,iEnd,0.0);
135 IxU = new math_Vector(1,iEnd); IyU = new math_Vector(1,iEnd); IzU = new math_Vector(1,iEnd);
136 IxxU = new math_Vector(1,iEnd); IyyU = new math_Vector(1,iEnd); IzzU = new math_Vector(1,iEnd);
137 IxyU = new math_Vector(1,iEnd); IxzU = new math_Vector(1,iEnd); IyzU = new math_Vector(1,iEnd);
139 return FillIntervalBounds(A, B, Knots, U1, U2);
142 static Standard_Real CCompute(Face& S, Domain& D, const Standard_Boolean ByPoint, const Standard_Real Coeff[],
143 const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia,
144 const Standard_Real EpsDim,
145 const Standard_Boolean isErrorCalculation, const Standard_Boolean isVerifyComputation)
147 Standard_Boolean isNaturalRestriction = S.NaturalRestriction();
149 Standard_Integer NumSubs = SUBS_POWER;
150 Standard_Boolean isMinDim = IS_MIN_DIM;
152 Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
153 Dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
154 //boundary curve parametrization
155 Standard_Real l1, l2, lm, lr, l;
156 //Face parametrization in U and V direction
157 Standard_Real BV1, BV2, v;
158 Standard_Real BU1, BU2, u1, u2, um, ur, u;
159 S.Bounds (BU1, BU2, BV1, BV2); u1 = BU1;
160 //location point used to compute the inertia
161 Standard_Real xloc, yloc, zloc;
162 loc.Coord (xloc, yloc, zloc);
163 //location point used to compute the inertiard (xloc, yloc, zloc);
164 //Jacobien (x, y, z) -> (u, v) = ||n||
165 Standard_Real xn, yn, zn, s, ds, dDim;
166 Standard_Real x, y, z, xi, px, py, pz, yi, zi, d1, d2, d3;
170 //On the boundary curve u-v
173 Standard_Real Dul; // Dul = Du / Dl
174 Standard_Real CDim[2], CIx, CIy, CIz, CIxx[2], CIyy[2], CIzz[2], CIxy, CIxz, CIyz;
175 Standard_Real LocDim[2], LocIx[2], LocIy[2], LocIz[2], LocIxx[2], LocIyy[2], LocIzz[2], LocIxy[2], LocIxz[2], LocIyz[2];
177 Standard_Integer iD = 0, NbLSubs, iLS, iLSubEnd, iGL, iGLEnd, NbLGaussP[2], LRange[2], iL, kL, kLEnd, IL, JL;
178 Standard_Integer i, NbUSubs, iUS, iUSubEnd, iGU, iGUEnd, NbUGaussP[2], URange[2], iU, kU, kUEnd, IU, JU;
179 Standard_Integer UMaxSubs, LMaxSubs;
181 Standard_Real ErrorU, ErrorL, ErrorLMax = 0.0, Eps=0.0, EpsL=0.0, EpsU=0.0;
182 iGLEnd = isErrorCalculation? 2: 1;
184 for(i = 0; i < 2; i++) {
197 NbUGaussP[0] = S.SIntOrder(EpsDim);
198 NbUGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbUGaussP[0]));
199 math::GaussPoints(NbUGaussP[0],UGaussP0); math::GaussWeights(NbUGaussP[0],UGaussW0);
200 math::GaussPoints(NbUGaussP[1],UGaussP1); math::GaussWeights(NbUGaussP[1],UGaussW1);
202 NbUSubs = S.SUIntSubs();
203 TColStd_Array1OfReal UKnots(1,NbUSubs+1);
206 while (isNaturalRestriction || D.More()) {
207 if(isNaturalRestriction){
208 NbLGaussP[0] = Min(2*NbUGaussP[0],math::GaussPointsMax());
210 S.Load(D.Value()); ++iD;
211 NbLGaussP[0] = S.LIntOrder(EpsDim);
213 NbLGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbLGaussP[0]));
214 math::GaussPoints(NbLGaussP[0],LGaussP0); math::GaussWeights(NbLGaussP[0],LGaussW0);
215 math::GaussPoints(NbLGaussP[1],LGaussP1); math::GaussWeights(NbLGaussP[1],LGaussW1);
217 NbLSubs = isNaturalRestriction? S.SVIntSubs(): S.LIntSubs();
218 TColStd_Array1OfReal LKnots(1,NbLSubs+1);
219 if(isNaturalRestriction){
224 l1 = S.FirstParameter(); l2 = S.LastParameter();
228 //OCC503(apo): if(Abs(l2-l1) < EPS_PARAM) continue;
229 if(Abs(l2-l1) > EPS_PARAM) {
230 iLSubEnd = LFillIntervalBounds(l1, l2, LKnots, NumSubs);
231 LMaxSubs = MaxSubs(iLSubEnd);
232 //-- exception avoiding
233 if(LMaxSubs > SM) LMaxSubs = SM;
234 DimL.Init(0.0,1,LMaxSubs); ErrL.Init(0.0,1,LMaxSubs); ErrUL.Init(0.0,1,LMaxSubs);
237 LRange[0] = IL = ErrL->Max(); LRange[1] = JL;
238 L1(JL) = (L1(IL) + L2(IL))/2.0; L2(JL) = L2(IL); L2(IL) = L1(JL);
239 }else LRange[0] = IL = JL;
240 if(JL == LMaxSubs || Abs(L2(JL) - L1(JL)) < EPS_PARAM)
242 DimL(JL) = ErrL(JL) = IxL(JL) = IyL(JL) = IzL(JL) =
243 IxxL(JL) = IyyL(JL) = IzzL(JL) = IxyL(JL) = IxzL(JL) = IyzL(JL) = 0.0;
246 EpsL = ErrorL; Eps = EpsL/0.9;
250 for(kL=0; kL < kLEnd; kL++){
252 lm = 0.5*(L2(iLS) + L1(iLS));
253 lr = 0.5*(L2(iLS) - L1(iLS));
254 CIx = CIy = CIz = CIxy = CIxz = CIyz = 0.0;
255 for(iGL=0; iGL < iGLEnd; iGL++){//
256 CDim[iGL] = CIxx[iGL] = CIyy[iGL] = CIzz[iGL] = 0.0;
257 for(iL=1; iL<=NbLGaussP[iGL]; iL++){
258 l = lm + lr*(*LGaussP[iGL])(iL);
259 if(isNaturalRestriction){
260 v = l; u2 = BU2; Dul = (*LGaussW[iGL])(iL);
262 S.D12d (l, Puv, Vuv);
263 Dul = Vuv.Y()*(*LGaussW[iGL])(iL); // Dul = Du / Dl
264 if(Abs(Dul) < EPS_PARAM) continue;
265 v = Puv.Y(); u2 = Puv.X();
266 //Check on cause out off bounds of value current parameter
267 if(v < BV1) v = BV1; else if(v > BV2) v = BV2;
268 if(u2 < BU1) u2 = BU1; else if(u2 > BU2) u2 = BU2;
272 if(Abs(u2-u1) < EPS_PARAM) continue;
273 iUSubEnd = UFillIntervalBounds(u1, u2, UKnots, NumSubs);
274 UMaxSubs = MaxSubs(iUSubEnd);
275 //-- exception avoiding
276 if(UMaxSubs > SM) UMaxSubs = SM;
277 DimU.Init(0.0,1,UMaxSubs); ErrU.Init(0.0,1,UMaxSubs); ErrorU = 0.0;
280 URange[0] = IU = ErrU->Max(); URange[1] = JU;
281 U1(JU) = (U1(IU)+U2(IU))/2.0; U2(JU) = U2(IU); U2(IU) = U1(JU);
282 }else URange[0] = IU = JU;
283 if(JU == UMaxSubs || Abs(U2(JU) - U1(JU)) < EPS_PARAM)
285 DimU(JU) = ErrU(JU) = IxU(JU) = IyU(JU) = IzU(JU) =
286 IxxU(JU) = IyyU(JU) = IzzU(JU) = IxyU(JU) = IxzU(JU) = IyzU(JU) = 0.0;
289 EpsU = ErrorU; Eps = EpsU*Abs((u2-u1)*Dul)/0.1; EpsL = 0.9*Eps;
293 for(kU=0; kU < kUEnd; kU++){
295 um = 0.5*(U2(iUS) + U1(iUS));
296 ur = 0.5*(U2(iUS) - U1(iUS));
297 iGUEnd = iGLEnd - iGL;
298 for(iGU=0; iGU < iGUEnd; iGU++){//
300 LocIxx[iGU] = LocIyy[iGU] = LocIzz[iGU] =
301 LocIx[iGU] = LocIy[iGU] = LocIz[iGU] =
302 LocIxy[iGU] = LocIxz[iGU] = LocIyz[iGU] = 0.0;
303 for(iU=1; iU<=NbUGaussP[iGU]; iU++){
304 u = um + ur*(*UGaussP[iGU])(iU);
305 S.Normal(u, v, Ps, VNor);
306 VNor.Coord(xn, yn, zn);
308 x -= xloc; y -= yloc; z -= zloc;
309 xn *= (*UGaussW[iGU])(iU);
310 yn *= (*UGaussW[iGU])(iU);
311 zn *= (*UGaussW[iGU])(iU);
313 //volume of elementary cone
314 dDim = (x*xn+y*yn+z*zn)/3.0;
315 //coordinates of cone's center mass
316 px = 0.75*x; py = 0.75*y; pz = 0.75*z;
318 //if(iGU > 0) continue;
319 LocIx[iGU] += px*dDim;
320 LocIy[iGU] += py*dDim;
321 LocIz[iGU] += pz*dDim;
322 x -= Coeff[0]; y -= Coeff[1]; z -= Coeff[2];
324 LocIxy[iGU] -= x*y*dDim;
325 LocIyz[iGU] -= y*z*dDim;
326 LocIxz[iGU] -= x*z*dDim;
327 xi = x*x; yi = y*y; zi = z*z;
328 LocIxx[iGU] += (yi+zi)*dDim;
329 LocIyy[iGU] += (xi+zi)*dDim;
330 LocIzz[iGU] += (xi+yi)*dDim;
332 s = xn*Coeff[0] + yn*Coeff[1] + zn*Coeff[2];
333 d1 = Coeff[0]*x + Coeff[1]*y + Coeff[2]*z - Coeff[3];
338 //if(iGU > 0) continue;
339 LocIx[iGU] += (x - Coeff[0]*d1/2.0) * ds;
340 LocIy[iGU] += (y - Coeff[1]*d1/2.0) * ds;
341 LocIz[iGU] += (z - Coeff[2]*d1/2.0) * ds;
342 px = x-Coeff[0]*d1; py = y-Coeff[1]*d1; pz = z-Coeff[2]*d1;
343 xi = px*px*d1 + px*Coeff[0]*d2 + Coeff[0]*Coeff[0]*d3;
344 yi = py*py*d1 + py*Coeff[1]*d2 + Coeff[1]*Coeff[1]*d3;
345 zi = pz*pz*d1 + pz*Coeff[2]*d2 + Coeff[2]*Coeff[2]*d3;
346 LocIxx[iGU] += (yi+zi)*s;
347 LocIyy[iGU] += (xi+zi)*s;
348 LocIzz[iGU] += (xi+yi)*s;
350 xi = py*pz*d1 + py*Coeff[2]*d2 + pz*Coeff[1]*d2 + Coeff[1]*Coeff[2]*d3;
351 yi = px*pz*d1 + pz*Coeff[0]*d2 + px*Coeff[2]*d2 + Coeff[0]*Coeff[2]*d3;
352 zi = px*py*d1 + px*Coeff[1]*d2 + py*Coeff[0]*d2 + Coeff[0]*Coeff[1]*d3;
353 LocIxy[iGU] -= zi*s; LocIyz[iGU] -= xi*s; LocIxz[iGU] -= yi*s;
357 DimU(iUS) = LocDim[0]*ur;
358 IxxU(iUS) = LocIxx[0]*ur; IyyU(iUS) = LocIyy[0]*ur; IzzU(iUS) = LocIzz[0]*ur;
359 if(iGL > 0) continue;
360 LocDim[1] = Abs(LocDim[1]-LocDim[0]);
361 LocIxx[1] = Abs(LocIxx[1]-LocIxx[0]);
362 LocIyy[1] = Abs(LocIyy[1]-LocIyy[0]);
363 LocIzz[1] = Abs(LocIzz[1]-LocIzz[0]);
364 ErrU(iUS) = isMinDim? LocDim[1]*ur: (LocIxx[1] + LocIyy[1] + LocIzz[1])*ur;
365 IxU(iUS) = LocIx[0]*ur; IyU(iUS) = LocIy[0]*ur; IzU(iUS) = LocIz[0]*ur;
366 IxyU(iUS) = LocIxy[0]*ur; IxzU(iUS) = LocIxz[0]*ur; IyzU(iUS) = LocIyz[0]*ur;
368 if(JU == iUSubEnd) kUEnd = 2;
370 Standard_Integer imax = ErrU->Max();
371 if(imax > 0) ErrorU = ErrU(imax);
374 }while((ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1);
375 for(i=1; i<=JU; i++) {
376 CDim[iGL] += DimU(i)*Dul;
377 CIxx[iGL] += IxxU(i)*Dul; CIyy[iGL] += IyyU(i)*Dul; CIzz[iGL] += IzzU(i)*Dul;
379 if(iGL > 0) continue;
380 ErrUL(iLS) = ErrorU*Abs((u2-u1)*Dul);
381 for(i=1; i<=JU; i++){
382 CIx += IxU(i)*Dul; CIy += IyU(i)*Dul; CIz += IzU(i)*Dul;
383 //CIxx += IxxU(i)*Dul; CIyy += IyyU(i)*Dul; CIzz += IzzU(i)*Dul;
384 CIxy += IxyU(i)*Dul; CIxz += IxzU(i)*Dul; CIyz += IyzU(i)*Dul;
388 DimL(iLS) = CDim[0]*lr;
389 IxxL(iLS) = CIxx[0]*lr; IyyL(iLS) = CIyy[0]*lr; IzzL(iLS) = CIzz[0]*lr;
391 //ErrL(iLS) = Abs(CDim[1]-CDim[0])*lr + ErrUL(iLS);
392 CDim[1] = Abs(CDim[1]-CDim[0]);
393 CIxx[1] = Abs(CIxx[1]-CIxx[0]); CIyy[1] = Abs(CIyy[1]-CIyy[0]); CIzz[1] = Abs(CIzz[1]-CIzz[0]);
395 ErrL(iLS) = (isMinDim? CDim[1]: (CIxx[1] + CIyy[1] + CIzz[1]))*lr + ErrorU;
397 IxL(iLS) = CIx*lr; IyL(iLS) = CIy*lr; IzL(iLS) = CIz*lr;
398 //IxxL(iLS) = CIxx*lr; IyyL(iLS) = CIyy*lr; IzzL(iLS) = CIzz*lr;
399 IxyL(iLS) = CIxy*lr; IxzL(iLS) = CIxz*lr; IyzL(iLS) = CIyz*lr;
401 // Calculate/correct epsilon of computation by current value of Dim
402 //That is need for not spend time for
405 Standard_Real DDim = 0.0, DIxx = 0.0, DIyy = 0.0, DIzz = 0.0;
406 for(i=1; i<=JL; i++) {
408 DIxx += IxxL(i); DIyy += IyyL(i); DIzz += IzzL(i);
410 DDim = isMinDim? Abs(DDim): Abs(DIxx) + Abs(DIyy) + Abs(DIzz);
411 DDim = Abs(DDim*EpsDim);
413 Eps = DDim; EpsL = 0.9*Eps;
417 Standard_Integer imax = ErrL->Max();
418 if(imax > 0) ErrorL = ErrL(imax);
421 }while((ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1);
422 for(i=1; i<=JL; i++){
424 Ix += IxL(i); Iy += IyL(i); Iz += IzL(i);
425 Ixx += IxxL(i); Iyy += IyyL(i); Izz += IzzL(i);
426 Ixy += IxyL(i); Ixz += IxzL(i); Iyz += IyzL(i);
428 ErrorLMax = Max(ErrorLMax, ErrorL);
430 if(isNaturalRestriction) break;
433 if(Abs(Dim) >= EPS_DIM){
435 Ix = Coeff[0] + Ix/Dim;
436 Iy = Coeff[1] + Iy/Dim;
437 Iz = Coeff[2] + Iz/Dim;
443 g.SetCoord (Ix, Iy, Iz);
446 g.SetCoord(0.,0.,0.);
448 inertia.SetCols (gp_XYZ (Ixx, Ixy, Ixz),
449 gp_XYZ (Ixy, Iyy, Iyz),
450 gp_XYZ (Ixz, Iyz, Izz));
452 Eps = Dim != 0.0? ErrorLMax/(isMinDim? Abs(Dim): (Abs(Ixx) + Abs(Iyy) + Abs(Izz))): 0.0;
457 static Standard_Real Compute(Face& S, const Standard_Boolean ByPoint, const Standard_Real Coeff[],
458 const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, Standard_Real EpsDim)
460 Standard_Boolean isErrorCalculation = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
461 Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
462 EpsDim = Abs(EpsDim);
464 return CCompute(S,D,ByPoint,Coeff,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
467 static Standard_Real Compute(Face& S, Domain& D, const Standard_Boolean ByPoint, const Standard_Real Coeff[],
468 const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, Standard_Real EpsDim)
470 Standard_Boolean isErrorCalculation = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
471 Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
472 EpsDim = Abs(EpsDim);
473 return CCompute(S,D,ByPoint,Coeff,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
476 static void Compute(const Face& S,
477 const Standard_Boolean ByPoint,
478 const Standard_Real Coeff[],
487 Standard_Real dvi, dv;
488 Standard_Real ur, um, u, vr, vm, v;
489 Standard_Real x, y, z, xn, yn, zn, xi, yi, zi;
490 // Standard_Real x, y, z, xn, yn, zn, xi, yi, zi, xyz;
491 Standard_Real px,py,pz,s,d1,d2,d3;
492 Standard_Real Ixi, Iyi, Izi, Ixxi, Iyyi, Izzi, Ixyi, Ixzi, Iyzi;
493 Standard_Real xloc, yloc, zloc;
494 Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
496 Volu = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
497 Loc.Coord (xloc, yloc, zloc);
499 Standard_Real LowerU, UpperU, LowerV, UpperV;
500 S.Bounds ( LowerU, UpperU, LowerV, UpperV);
501 Standard_Integer UOrder = Min(S.UIntegrationOrder (),
502 math::GaussPointsMax());
503 Standard_Integer VOrder = Min(S.VIntegrationOrder (),
504 math::GaussPointsMax());
506 Standard_Integer i, j;
507 math_Vector GaussPU (1, UOrder); //gauss points and weights
508 math_Vector GaussWU (1, UOrder);
509 math_Vector GaussPV (1, VOrder);
510 math_Vector GaussWV (1, VOrder);
512 math::GaussPoints (UOrder,GaussPU);
513 math::GaussWeights (UOrder,GaussWU);
514 math::GaussPoints (VOrder,GaussPV);
515 math::GaussWeights (VOrder,GaussWV);
517 um = 0.5 * (UpperU + LowerU);
518 vm = 0.5 * (UpperV + LowerV);
519 ur = 0.5 * (UpperU - LowerU);
520 vr = 0.5 * (UpperV - LowerV);
522 for (j = 1; j <= VOrder; j++) {
523 v = vm + vr * GaussPV (j);
524 dvi = Ixi = Iyi = Izi = Ixxi = Iyyi = Izzi = Ixyi = Ixzi = Iyzi = 0.0;
526 for (i = 1; i <= UOrder; i++) {
527 u = um + ur * GaussPU (i);
528 S.Normal (u, v, P, VNor);
529 VNor.Coord (xn, yn, zn);
531 x -= xloc; y -= yloc; z -= zloc;
532 xn *= GaussWU (i); yn *= GaussWU (i); zn *= GaussWU (i);
534 ///////////////////// ///////////////////////
535 // OFV code // // Initial code //
536 ///////////////////// ///////////////////////
538 dv = (x*xn+y*yn+z*zn)/3.0; //xyz = x * y * z;
539 dvi += dv; //Ixyi += zn * xyz;
540 Ixi += 0.75*x*dv; //Iyzi += xn * xyz;
541 Iyi += 0.75*y*dv; //Ixzi += yn * xyz;
542 Izi += 0.75*z*dv; //xi = x * x * x * xn / 3.0;
543 x -= Coeff[0]; //yi = y * y * y * yn / 3.0;
544 y -= Coeff[1]; //zi = z * z * z * zn / 3.0;
545 z -= Coeff[2]; //Ixxi += (yi + zi);
546 dv *= 3.0/5.0; //Iyyi += (xi + zi);
547 Ixyi -= x*y*dv; //Izzi += (xi + yi);
548 Iyzi -= y*z*dv; //x -= Coeff[0];
549 Ixzi -= x*z*dv; //y -= Coeff[1];
550 xi = x*x; //z -= Coeff[2];
551 yi = y*y; //dv = x * xn + y * yn + z * zn;
552 zi = z*z; //dvi += dv;
553 Ixxi += (yi + zi)*dv; //Ixi += x * dv;
554 Iyyi += (xi + zi)*dv; //Iyi += y * dv;
555 Izzi += (xi + yi)*dv; //Izi += z * dv;
558 s = xn * Coeff[0] + yn * Coeff[1] + zn * Coeff[2];
559 d1 = Coeff[0] * x + Coeff[1] * y + Coeff[2] * z - Coeff[3];
564 Ixi += (x - (Coeff[0] * d1 / 2.0)) * dv;
565 Iyi += (y - (Coeff[1] * d1 / 2.0)) * dv;
566 Izi += (z - (Coeff[2] * d1 / 2.0)) * dv;
567 px = x - Coeff[0] * d1;
568 py = y - Coeff[1] * d1;
569 pz = z - Coeff[2] * d1;
570 xi = px * px * d1 + px * Coeff[0]* d2 + Coeff[0] * Coeff[0] * d3;
571 yi = py * py * d1 + py * Coeff[1] * d2 + Coeff[1] * Coeff[1] * d3;
572 zi = pz * pz * d1 + pz * Coeff[2] * d2 + Coeff[2] * Coeff[2] * d3;
573 Ixxi += (yi + zi) * s;
574 Iyyi += (xi + zi) * s;
575 Izzi += (xi + yi) * s;
577 xi = (py * pz * d1) + (py * Coeff[2] * d2) + (pz * Coeff[1] * d2) + (Coeff[1] * Coeff[2] * d3);
578 yi = (px * pz * d1) + (pz * Coeff[0] * d2) + (px * Coeff[2] * d2) + (Coeff[0] * Coeff[2] * d3);
579 zi = (px * py * d1) + (px * Coeff[1] * d2) + (py * Coeff[0] * d2) + (Coeff[0] * Coeff[1] * d3);
585 Volu += dvi * GaussWV (j);
586 Ix += Ixi * GaussWV (j);
587 Iy += Iyi * GaussWV (j);
588 Iz += Izi * GaussWV (j);
589 Ixx += Ixxi * GaussWV (j);
590 Iyy += Iyyi * GaussWV (j);
591 Izz += Izzi * GaussWV (j);
592 Ixy += Ixyi * GaussWV (j);
593 Ixz += Ixzi * GaussWV (j);
594 Iyz += Iyzi * GaussWV (j);
603 if (Abs(Volu) >= EPS_DIM ) {
605 Ix = Coeff[0] + Ix/Volu;
606 Iy = Coeff[1] + Iy/Volu;
607 Iz = Coeff[2] + Iz/Volu;
616 G.SetCoord (Ix, Iy, Iz);
619 G.SetCoord(0.,0.,0.);
622 Inertia.SetCols (gp_XYZ (Ixx, Ixy, Ixz),
623 gp_XYZ (Ixy, Iyy, Iyz),
624 gp_XYZ (Ixz, Iyz, Izz));
628 // Last modified by OFV 5.2001:
629 // 1). surface and edge integration order is equal now
630 // 2). "by point" works now rathre correctly (it looks so...)
631 static void Compute(Face& S, Domain& D, const Standard_Boolean ByPoint, const Standard_Real Coeff[],
632 const gp_Pnt& Loc, Standard_Real& Volu, gp_Pnt& G, gp_Mat& Inertia)
635 Standard_Real x, y, z, xi, yi, zi, l1, l2, lm, lr, l, v1, v2, v, u1, u2, um, ur, u, ds, Dul, xloc, yloc, zloc;
636 Standard_Real LocVolu, LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy, LocIxz, LocIyz;
637 Standard_Real CVolu, CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz, Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
638 Standard_Real xn, yn, zn, px, py, pz, s, d1, d2, d3, dSigma;
639 Standard_Integer i, j, vio, sio, max, NbGaussgp_Pnts;
646 Loc.Coord (xloc, yloc, zloc);
647 Volu = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
648 S.Bounds (u1, u2, v1, v2);
649 Standard_Real _u2 = u2; //OCC104
650 vio = S.VIntegrationOrder ();
655 sio = S.IntegrationOrder ();
657 NbGaussgp_Pnts = Min(max,math::GaussPointsMax());
659 math_Vector GaussP (1, NbGaussgp_Pnts);
660 math_Vector GaussW (1, NbGaussgp_Pnts);
661 math::GaussPoints (NbGaussgp_Pnts,GaussP);
662 math::GaussWeights (NbGaussgp_Pnts,GaussW);
664 CVolu = CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
665 l1 = S.FirstParameter();
666 l2 = S.LastParameter();
667 lm = 0.5 * (l2 + l1);
668 lr = 0.5 * (l2 - l1);
670 for (i=1; i<=NbGaussgp_Pnts; i++)
672 l = lm + lr * GaussP(i);
673 S.D12d (l, Puv, Vuv);
680 u2 = u2 < u1? u1: u2;
681 u2 = u2 > _u2? _u2: u2;
683 Dul = Vuv.Y() * GaussW(i);
684 um = 0.5 * (u2 + u1);
685 ur = 0.5 * (u2 - u1);
686 LocVolu = LocIx = LocIy = LocIz = LocIxx = LocIyy = LocIzz = LocIxy = LocIxz = LocIyz = 0.0;
688 for (j=1; j<=NbGaussgp_Pnts; j++)
690 u = um + ur * GaussP(j);
691 S.Normal (u, v, Ps, VNor);
692 VNor.Coord (xn, yn, zn);
697 xn = xn * Dul * GaussW(j);
698 yn = yn * Dul * GaussW(j);
699 zn = zn * Dul * GaussW(j);
702 dSigma = (x*xn+y*yn+z*zn)/3.0;
704 LocIx += 0.75*x*dSigma;
705 LocIy += 0.75*y*dSigma;
706 LocIz += 0.75*z*dSigma;
711 LocIxy -= x*y*dSigma;
712 LocIyz -= y*z*dSigma;
713 LocIxz -= x*z*dSigma;
717 LocIxx += (yi + zi)*dSigma;
718 LocIyy += (xi + zi)*dSigma;
719 LocIzz += (xi + yi)*dSigma;
723 s = xn * Coeff[0] + yn * Coeff[1] + zn * Coeff[2];
724 d1 = Coeff[0] * x + Coeff[1] * y + Coeff[2] * z;
729 LocIx += (x - Coeff[0] * d1 / 2.0) * ds;
730 LocIy += (y - Coeff[1] * d1 / 2.0) * ds;
731 LocIz += (z - Coeff[2] * d1 / 2.0) * ds;
732 px = x - Coeff[0] * d1;
733 py = y - Coeff[1] * d1;
734 pz = z - Coeff[2] * d1;
735 xi = (px * px * d1) + (px * Coeff[0]* d2) + (Coeff[0] * Coeff[0] * d3);
736 yi = (py * py * d1) + (py * Coeff[1] * d2) + (Coeff[1] * Coeff[1] * d3);
737 zi = pz * pz * d1 + pz * Coeff[2] * d2 + (Coeff[2] * Coeff[2] * d3);
738 LocIxx += (yi + zi) * s;
739 LocIyy += (xi + zi) * s;
740 LocIzz += (xi + yi) * s;
742 xi = (py * pz * d1) + (py * Coeff[2] * d2) + (pz * Coeff[1] * d2) + (Coeff[1] * Coeff[2] * d3);
743 yi = (px * pz * d1) + (pz * Coeff[0] * d2) + (px * Coeff[2] * d2) + (Coeff[0] * Coeff[2] * d3);
744 zi = (px * py * d1) + (px * Coeff[1] * d2) + (py * Coeff[0] * d2) + (Coeff[0] * Coeff[1] * d3);
750 CVolu += LocVolu * ur;
774 if(Abs(Volu) >= EPS_DIM)
778 Ix = Coeff[0] + Ix/Volu;
779 Iy = Coeff[1] + Iy/Volu;
780 Iz = Coeff[2] + Iz/Volu;
788 G.SetCoord (Ix, Iy, Iz);
793 G.SetCoord(0.,0.,0.);
796 Inertia.SetCols (gp_XYZ (Ixx, Ixy, Ixz),
797 gp_XYZ (Ixy, Iyy, Iyz),
798 gp_XYZ (Ixz, Iyz, Izz));
802 GProp_VGProps::GProp_VGProps(){}
804 GProp_VGProps::GProp_VGProps(Face& S, const gp_Pnt& VLocation, const Standard_Real Eps){
805 SetLocation(VLocation);
809 GProp_VGProps::GProp_VGProps(Face& S, Domain& D, const gp_Pnt& VLocation, const Standard_Real Eps){
810 SetLocation(VLocation);
814 GProp_VGProps::GProp_VGProps(Face& S, Domain& D, const gp_Pnt& VLocation){
815 SetLocation(VLocation);
819 GProp_VGProps::GProp_VGProps(const Face& S, const gp_Pnt& VLocation){
820 SetLocation(VLocation);
824 GProp_VGProps::GProp_VGProps(Face& S, const gp_Pnt& O, const gp_Pnt& VLocation, const Standard_Real Eps){
825 SetLocation(VLocation);
829 GProp_VGProps::GProp_VGProps(Face& S, Domain& D, const gp_Pnt& O, const gp_Pnt& VLocation, const Standard_Real Eps){
830 SetLocation(VLocation);
834 GProp_VGProps::GProp_VGProps(const Face& S, const gp_Pnt& O, const gp_Pnt& VLocation){
835 SetLocation(VLocation);
839 GProp_VGProps::GProp_VGProps(Face& S, Domain& D, const gp_Pnt& O, const gp_Pnt& VLocation){
840 SetLocation(VLocation);
844 GProp_VGProps::GProp_VGProps(Face& S, const gp_Pln& Pl, const gp_Pnt& VLocation, const Standard_Real Eps){
845 SetLocation(VLocation);
849 GProp_VGProps::GProp_VGProps(Face& S, Domain& D, const gp_Pln& Pl, const gp_Pnt& VLocation, const Standard_Real Eps){
850 SetLocation(VLocation);
854 GProp_VGProps::GProp_VGProps(const Face& S, const gp_Pln& Pl, const gp_Pnt& VLocation){
855 SetLocation(VLocation);
859 GProp_VGProps::GProp_VGProps(Face& S, Domain& D, const gp_Pln& Pl, const gp_Pnt& VLocation){
860 SetLocation(VLocation);
864 void GProp_VGProps::SetLocation(const gp_Pnt& VLocation){
868 Standard_Real GProp_VGProps::Perform(Face& S, const Standard_Real Eps){
869 Standard_Real Coeff[] = {0., 0., 0.};
870 return myEpsilon = Compute(S,Standard_True,Coeff,loc,dim,g,inertia,Eps);
873 Standard_Real GProp_VGProps::Perform(Face& S, Domain& D, const Standard_Real Eps){
874 Standard_Real Coeff[] = {0., 0., 0.};
875 return myEpsilon = Compute(S,D,Standard_True,Coeff,loc,dim,g,inertia,Eps);
878 void GProp_VGProps::Perform(const Face& S){
879 Standard_Real Coeff[] = {0., 0., 0.};
880 Compute(S,Standard_True,Coeff,loc,dim,g,inertia);
885 void GProp_VGProps::Perform(Face& S, Domain& D){
886 Standard_Real Coeff[] = {0., 0., 0.};
887 Compute(S,D,Standard_True,Coeff,loc,dim,g,inertia);
892 Standard_Real GProp_VGProps::Perform(Face& S, const gp_Pnt& O, const Standard_Real Eps){
893 Standard_Real xloc, yloc, zloc;
894 loc.Coord(xloc, yloc, zloc);
895 Standard_Real Coeff[3];
896 O.Coord (Coeff[0], Coeff[1], Coeff[2]);
897 Coeff[0] -= xloc; Coeff[1] -= yloc; Coeff[2] -= zloc;
898 return myEpsilon = Compute(S,Standard_True,Coeff,loc,dim,g,inertia,Eps);
901 Standard_Real GProp_VGProps::Perform(Face& S, Domain& D, const gp_Pnt& O, const Standard_Real Eps){
902 Standard_Real xloc, yloc, zloc;
903 loc.Coord(xloc, yloc, zloc);
904 Standard_Real Coeff[3];
905 O.Coord (Coeff[0], Coeff[1], Coeff[2]);
906 Coeff[0] -= xloc; Coeff[1] -= yloc; Coeff[2] -= zloc;
907 return myEpsilon = Compute(S,D,Standard_True,Coeff,loc,dim,g,inertia,Eps);
910 void GProp_VGProps::Perform(const Face& S, const gp_Pnt& O){
911 Standard_Real xloc, yloc, zloc;
912 loc.Coord(xloc, yloc, zloc);
913 Standard_Real Coeff[3];
914 O.Coord (Coeff[0], Coeff[1], Coeff[2]);
915 Coeff[0] -= xloc; Coeff[1] -= yloc; Coeff[2] -= zloc;
916 Compute(S,Standard_True,Coeff,loc,dim,g,inertia);
921 void GProp_VGProps::Perform(Face& S, Domain& D, const gp_Pnt& O){
922 Standard_Real xloc, yloc, zloc;
923 loc.Coord(xloc, yloc, zloc);
924 Standard_Real Coeff[3];
925 O.Coord (Coeff[0], Coeff[1], Coeff[2]);
926 Coeff[0] -= xloc; Coeff[1] -= yloc; Coeff[2] -= zloc;
927 Compute(S,D,Standard_True,Coeff,loc,dim,g,inertia);
932 Standard_Real GProp_VGProps::Perform(Face& S, const gp_Pln& Pl, const Standard_Real Eps){
933 Standard_Real xloc, yloc, zloc;
934 loc.Coord (xloc, yloc, zloc);
935 Standard_Real Coeff[4];
936 Pl.Coefficients (Coeff[0], Coeff[1],Coeff[2],Coeff[3]);
937 Coeff[3] = Coeff[3] - Coeff[0]*xloc - Coeff[1]*yloc - Coeff[2]*zloc;
938 return myEpsilon = Compute(S,Standard_False,Coeff,loc,dim,g,inertia,Eps);
941 Standard_Real GProp_VGProps::Perform(Face& S, Domain& D, const gp_Pln& Pl, const Standard_Real Eps){
942 Standard_Real xloc, yloc, zloc;
943 loc.Coord (xloc, yloc, zloc);
944 Standard_Real Coeff[4];
945 Pl.Coefficients (Coeff[0], Coeff[1],Coeff[2],Coeff[3]);
946 Coeff[3] = Coeff[3] - Coeff[0]*xloc - Coeff[1]*yloc - Coeff[2]*zloc;
947 return myEpsilon = Compute(S,D,Standard_False,Coeff,loc,dim,g,inertia,Eps);
950 void GProp_VGProps::Perform(const Face& S, const gp_Pln& Pl){
951 Standard_Real xloc, yloc, zloc;
952 loc.Coord (xloc, yloc, zloc);
953 Standard_Real Coeff[4];
954 Pl.Coefficients (Coeff[0], Coeff[1],Coeff[2],Coeff[3]);
955 Coeff[3] = Coeff[3] - Coeff[0]*xloc - Coeff[1]*yloc - Coeff[2]*zloc;
956 Compute(S,Standard_False,Coeff,loc,dim,g,inertia);
961 void GProp_VGProps::Perform(Face& S, Domain& D, const gp_Pln& Pl){
962 Standard_Real xloc, yloc, zloc;
963 loc.Coord (xloc, yloc, zloc);
964 Standard_Real Coeff[4];
965 Pl.Coefficients (Coeff[0], Coeff[1],Coeff[2],Coeff[3]);
966 Coeff[3] = Coeff[3] - Coeff[0]*xloc - Coeff[1]*yloc - Coeff[2]*zloc;
967 Compute(S,D,Standard_False,Coeff,loc,dim,g,inertia);
972 Standard_Real GProp_VGProps::GetEpsilon(){