Integration of OCCT 6.5.0 from SVN

src/GProp/GProp.cdl

@@ -0,0 +1,159 @@
+-- File:	GProp.cdl
+-- Created:	Mon Aug 24 17:41:08 1992
+-- Created:	Tue Mar 12 17:51:12 1991
+-- Author:	Michel CHAUVAT
+--              JCV - January 1992
+---Copyright:	Matra Datavision 1992
+
+
+
+package GProp
+
+        --- Purpose: 
+        --  This package defines algorithmes to compute the global properties 
+        --  of a set of points, a curve, a surface, a solid (non infinite 
+        --  region of space delimited with geometric entities), a compound
+        --  geometric system (heterogeneous composition of the previous 
+        --  entities).
+        --    
+        --  Global properties are :
+        --  . length, area, volume,
+        --  . centre of mass,
+        --  . axis of inertia,
+        --  . moments of inertia,
+        --  . radius of gyration.
+        --  
+        --  It provides  also a class to  compile the average point or
+        --  line of a set of points.
+
+uses 
+    Standard, 
+    TColStd, 
+    TColgp, 
+    gp, 
+    math, 
+    GeomAbs
+
+is
+
+exception UndefinedAxis inherits DomainError;
+        --- Purpose :  
+        --  This exception is raised when a method makes reference to
+        --  an undefined inertia axis of symmetry.
+        
+
+enumeration EquaType
+    is      Plane, Line, Point, Space, None end;
+
+enumeration ValueType
+    is      Mass,  
+            CenterMassX, CenterMassY, CenterMassZ, 
+            InertiaXX, InertiaYY, InertiaZZ, 
+            InertiaXY, InertiaXZ, InertiaYZ, 
+            Unknown 
+    end;
+
+
+
+        --- Purpose : Algorithmes :
+
+
+  class GProps;
+        --- Purpose :
+        --  Computes the global properties of a compound geometric 
+        --  system in 3d space.  It gives facilities to compose the
+        --  properties of hetegogeneous elements of the system 
+        --  (PGProps, CGProps, SGProps, VGProps or GProps). A density 
+        --  can be associated with each component of the system.
+
+
+       class PGProps;
+           --- Purpose :
+           --  Computes the global properties of a set of points in 3d.
+           --  This class inherits GProps.
+       
+         
+       generic class CGProps;
+           ---Purpose :
+           --   Computes the global properties  of a bounded
+           --   curve in 3d. This class inherits GProps.
+
+       class CelGProps;
+           ---Purpose :
+           --   Computes the global properties  of a gp curve in 3d 
+           --   This class inherits GProps.
+
+       generic class SGProps;
+           ---Purpose :
+           --   Computes the global properties and the area of a bounded
+           --   surface in 3d. This class inherits GProps.
+
+
+       class SelGProps;
+           ---Purpose :
+           --   Computes the global properties and the area of a bounded
+           --   elementary surface in 3d. This class inherits GProps.
+       
+       generic class VGProps;
+           ---Purpose :
+           --   Computes the global properties and the volume of a region
+           --   of space. This class inherits GProps.
+
+       generic class VGPropsGK, UFunction, TFunction;
+           ---Purpose :
+           --   Computes the global properties and the volume of a region
+           --   of space by adaptive Gauss-Kronrod integration.  
+           --   This class inherits GProps.
+
+
+       class VelGProps;
+           ---Purpose :
+           --   Computes the global properties and the volume of a region
+           --   of space. the region of space is defined by an elementary 
+           --   surface. This class inherits GProps.
+
+
+  class PrincipalProps;
+        ---Purpose :
+        --  Returns the principal inertia properties of a GProps. 
+
+
+
+
+        --- Purpose :  
+        --  The following abstract classes define templates
+        --  with the minimum of methods required to implement 
+        --  the computation of the global properties for a curve
+        --  or a surface.
+
+
+  deferred generic class CurveTool;
+
+  deferred generic class FaceTool;
+  
+  deferred generic class DomainTool;  
+  
+  --
+  --   Class to compute the average plane or line of a set of points.
+  --   
+  
+  class PEquation;
+    
+        --- Purpose : methods of package
+ 
+  HOperator (G, Q : Pnt from gp; Mass : Real; Operator : out Mat from gp);
+        --- Purpose : Computes the matrix Operator, referred to as the
+-- "Huyghens Operator" of a geometric system at the
+-- point Q of the space, using the following data :
+-- - Mass, i.e. the mass of the system,
+-- - G, the center of mass of the system.
+--   The "Huyghens Operator" is used to compute
+-- Inertia/Q, the matrix of inertia of the system at
+-- the point Q using Huyghens' theorem :
+--    Inertia/Q = Inertia/G + HOperator (Q, G, Mass)
+-- where Inertia/G is the matrix of inertia of the
+-- system relative to its center of mass as returned by
+-- the function MatrixOfInertia on any GProp_GProps object.
+
+
+end GProp;

src/GProp/GProp.cxx

@@ -0,0 +1,28 @@
+#include <GProp.ixx>
+
+#include <Standard_DimensionError.hxx>
+
+#include <gp.hxx>
+#include <gp_XYZ.hxx>
+
+
+void GProp::HOperator (
+
+		       const gp_Pnt& G, 
+		       const gp_Pnt& Q, 
+		       const Standard_Real Mass,
+		       gp_Mat&       Operator
+
+) {
+
+  gp_XYZ QG = G.XYZ() - Q.XYZ();
+  Standard_Real Ixx = QG.Y() * QG.Y() + QG.Z() * QG.Z();
+  Standard_Real Iyy = QG.X() * QG.X() + QG.Z() * QG.Z();
+  Standard_Real Izz = QG.Y() * QG.Y() + QG.X() * QG.X();
+  Standard_Real Ixy =  - QG.X() * QG.Y();
+  Standard_Real Iyz =  - QG.Y() * QG.Z();
+  Standard_Real Ixz =  - QG.X() * QG.Z();
+  Operator.SetCols (gp_XYZ (Ixx, Ixy, Ixz), gp_XYZ (Ixy, Iyy, Iyz),
+                    gp_XYZ (Ixz, Iyz, Izz));
+  Operator.Multiply (Mass);
+}

src/GProp/GProp_CGProps.cdl

@@ -0,0 +1,38 @@
+-- File:	CGProps.cdl
+-- Created:	Thu Aug 27 10:09:41 1992
+-- Created:	Thu Apr 11 17:30:05 1991
+-- Author:	Michel CHAUVAT
+--              Jean-Claude Vauthier January 1992, September 1992
+---Copyright:	Matra Datavision 1992
+
+
+
+generic class CGProps  from GProp (Curve     as any;
+    	    	    	    	   Tool      as any -- as CurveTool(Curve)
+    	    	    	    	  )
+
+inherits GProps from GProp
+
+	--- Purpose  : 
+	--  Computes the  global properties of bounded curves 
+	--  in 3D space. The curve must have at least a continuity C1. 
+	--  It can be a curve as defined in the template CurveTool from 
+	--  package GProp. This template gives the minimum of methods 
+	--  required to evaluate the global properties of a curve 3D with  
+	--  the algorithmes of GProp.
+
+uses  Pnt   from gp
+       
+is
+
+    Create returns CGProps;
+  
+    Create (C : Curve; CLocation : Pnt)  returns CGProps;
+
+    SetLocation(me : in out;CLocation : Pnt) ;
+
+    Perform(me : in out; C : Curve);
+
+end CGProps;
+
+

src/GProp/GProp_CGProps.gxx

@@ -0,0 +1,143 @@
+#include <math.hxx>
+#include <math_Vector.hxx>
+#include <gp.hxx>
+#include <gp_Vec.hxx>
+#include <Standard_NotImplemented.hxx>
+
+#include <TColStd_Array1OfReal.hxx>
+
+
+GProp_CGProps::GProp_CGProps(){}
+
+void GProp_CGProps::SetLocation(const gp_Pnt& CLocation)
+{
+  loc = CLocation;
+}
+
+void GProp_CGProps::Perform (const Curve& C)
+{
+
+  Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
+  dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
+
+  Standard_Real Lower    = Tool::FirstParameter  (C);
+  Standard_Real Upper    = Tool::LastParameter   (C);
+  Standard_Integer Order = Min(Tool::IntegrationOrder (C),
+			       math::GaussPointsMax());
+  
+  gp_Pnt P;    //value on the curve
+  gp_Vec V1;   //first derivative on the curve
+  Standard_Real ds;  //curvilign abscissae
+  Standard_Real ur, um, u;
+  Standard_Real x, y, z; 
+  Standard_Real xloc, yloc, zloc;
+
+  math_Vector GaussP (1, Order);
+  math_Vector GaussW (1, Order);
+  
+  //Recuperation des points de Gauss dans le fichier GaussPoints.
+  math::GaussPoints  (Order,GaussP);
+  math::GaussWeights (Order,GaussW);
+
+  // modified by NIZHNY-MKK  Thu Jun  9 12:13:21 2005.BEGIN
+  Standard_Integer nbIntervals = Tool::NbIntervals(C, GeomAbs_CN);
+  Standard_Boolean bHasIntervals = (nbIntervals > 1);
+  TColStd_Array1OfReal TI(1, nbIntervals + 1);
+
+  if(bHasIntervals) {
+    Tool::Intervals(C, TI, GeomAbs_CN);
+  }
+  else {
+    nbIntervals = 1;
+  }
+  Standard_Integer nIndex = 0;
+  Standard_Real UU1 = Min(Lower, Upper);
+  Standard_Real UU2 = Max(Lower, Upper);
+  
+  for(nIndex = 1; nIndex <= nbIntervals; nIndex++) {
+    if(bHasIntervals) {
+      Lower = Max(TI(nIndex), UU1);
+      Upper = Min(TI(nIndex+1), UU2);
+    }
+    else {
+      Lower = UU1;
+      Upper = UU2;
+    }
+
+    Standard_Real dimLocal, IxLocal, IyLocal, IzLocal, IxxLocal, IyyLocal, IzzLocal, IxyLocal, IxzLocal, IyzLocal;
+    dimLocal = IxLocal = IyLocal = IzLocal = IxxLocal = IyyLocal = IzzLocal = IxyLocal = IxzLocal = IyzLocal = 0.0;
+  // modified by NIZHNY-MKK  Thu Jun  9 12:13:32 2005.END
+
+    loc.Coord (xloc, yloc, zloc);
+
+    Standard_Integer i;
+
+    // Calcul des integrales aux points de gauss :
+    um = 0.5 * (Upper + Lower);
+    ur = 0.5 * (Upper - Lower);
+
+    for (i = 1; i <= Order; i++) {
+      u   = um + ur * GaussP (i);
+      Tool::D1 (C,u, P, V1); 
+      ds  = V1.Magnitude();
+      P.Coord (x, y, z);
+      x   -= xloc;
+      y   -= yloc;
+      z   -= zloc;
+      ds  *= GaussW (i);
+      dimLocal += ds; 
+      IxLocal  += x * ds;  
+      IyLocal  += y * ds;
+      IzLocal  += z * ds;
+      IxyLocal += x * y * ds;
+      IyzLocal += y * z * ds;
+      IxzLocal += x * z * ds;
+      x   *= x;      
+      y   *= y;      
+      z   *= z;      
+      IxxLocal += (y + z) * ds;
+      IyyLocal += (x + z) * ds;
+      IzzLocal += (x + y) * ds;
+    }
+    // modified by NIZHNY-MKK  Thu Jun  9 12:13:47 2005.BEGIN
+    dimLocal *= ur;
+    IxLocal  *= ur;
+    IyLocal  *= ur;
+    IzLocal  *= ur;
+    IxxLocal *= ur;
+    IyyLocal *= ur;
+    IzzLocal *= ur;
+    IxyLocal *= ur;
+    IxzLocal *= ur;
+    IyzLocal *= ur;
+
+    dim += dimLocal;
+    Ix += IxLocal;
+    Iy += IyLocal;
+    Iz += IzLocal;
+    Ixx += IxxLocal;
+    Iyy += IyyLocal;
+    Izz += IzzLocal;
+    Ixy += IxyLocal;
+    Ixz += IxzLocal;
+    Iyz += IyzLocal;
+  }
+  // modified by NIZHNY-MKK  Thu Jun  9 12:13:55 2005.END
+
+  inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
+		    gp_XYZ (-Ixy, Iyy, -Iyz),
+		    gp_XYZ (-Ixz, -Iyz, Izz));
+
+  if (Abs(dim) < gp::Resolution())
+    g = P;
+  else
+    g.SetCoord (Ix/dim, Iy/dim, Iz/dim);
+}
+
+
+GProp_CGProps::GProp_CGProps (const Curve& C, 
+			      const gp_Pnt&   CLocation)
+{
+  SetLocation(CLocation);
+  Perform(C);
+}

src/GProp/GProp_CelGProps.cdl

@@ -0,0 +1,39 @@
+-- File:	GProp_CelGProps.cdl
+-- Created:	Wed Dec  2 16:22:06 1992
+-- Author:	Isabelle GRIGNON
+--		<isg@sdsun2>
+---Copyright:	 Matra Datavision 1992
+
+
+class CelGProps  from GProp inherits GProps from GProp
+
+	--- Purpose  : 
+	--  Computes the  global properties of bounded curves 
+	--  in 3D space. 
+	--  It can be an elementary curve from package gp such as 
+	--  Lin, Circ, Elips, Parab .
+
+uses  Circ  from gp,
+      Lin   from gp,
+      Parab from gp,
+      Pnt   from gp
+       
+is
+
+    Create returns CelGProps;
+  
+    Create (C : Circ from gp; CLocation : Pnt)   returns CelGProps;
+
+    Create (C : Circ from gp; U1, U2 : Real; CLocation : Pnt)returns CelGProps;
+
+    Create (C : Lin from gp; U1, U2 : Real; CLocation : Pnt) returns CelGProps;
+
+    SetLocation(me : in out;CLocation : Pnt) ;
+
+    Perform(me : in out; C :Circ; U1,U2 : Real);
+    
+    Perform(me : in out; C : Lin ; U1,U2 : Real);
+
+end CelGProps;
+
+

src/GProp/GProp_CelGProps.cxx

@@ -0,0 +1,143 @@
+#include <GProp_CelGProps.ixx>
+#include <GProp.hxx>
+#include <gp.hxx>
+#include <gp_Vec.hxx>
+#include <gp.hxx>
+#include <ElCLib.hxx>
+#include <Standard_NotImplemented.hxx>
+#include <math_Matrix.hxx>
+#include <math_Vector.hxx>
+#include <math_Jacobi.hxx>
+
+GProp_CelGProps::GProp_CelGProps(){}
+
+void GProp_CelGProps::SetLocation(const gp_Pnt& CLocation)
+{
+  loc = CLocation;
+}
+
+
+void GProp_CelGProps::Perform(const gp_Circ& C, 			      
+			      const Standard_Real U1,
+			      const Standard_Real U2)
+{
+  Standard_Real X0,Y0,Z0,Xa1,Ya1,Za1,Xa2,Ya2,Za2,Xa3,Ya3,Za3;
+  C.Location().Coord(X0,Y0,Z0);
+  C.XAxis().Direction().Coord(Xa1,Ya1,Za1);
+  C.YAxis().Direction().Coord(Xa2,Ya2,Za2);
+  C.Axis().Direction().Coord(Xa3,Ya3,Za3);
+  Standard_Real Ray = C.Radius();
+
+  dim = Ray*Abs(U2-U1);
+  Standard_Real xloc =Ray*(Sin(U2)-Sin(U1))/(U2-U1);
+  Standard_Real yloc =Ray*(Cos(U1)-Cos(U2))/(U2-U1);
+
+  g.SetCoord(xloc*Xa1+yloc*Xa2+X0,xloc*Ya1+yloc*Ya2+Y0,Z0);
+
+  math_Matrix Dm(1,3,1,3);
+  Dm(1,1)= Ray*Ray*Ray*(U2/2-U1/2-Sin(2*U2)/4+Sin(2*U1)/4);
+  Dm(2,2)= Ray*Ray*Ray*(U2/2-U1/2+Sin(2*U2)/4-Sin(2*U1)/4);
+  Dm(3,3)= Ray*Ray*dim; 
+  Dm(2,1)= Dm(1,2) = - Ray*Ray*Ray*(Cos(2*U1)/4-Cos(2*U2)/4);
+  Dm(3,1) = Dm(1,3) = 0.;
+  Dm(3,2) = Dm(2,3) = 0.;
+  math_Matrix Passage (1,3,1,3);
+  Passage(1,1) = Xa1; Passage(1,2) = Xa2 ;Passage(1,3) = Xa3;
+  Passage(2,1) = Ya1; Passage(2,2) = Ya2 ;Passage(2,3) = Ya3;
+  Passage(3,1) = Za1; Passage(3,2) = Za2 ;Passage(3,3) = Za3;
+
+  math_Jacobi J(Dm);
+  math_Vector V1(1,3),V2(1,3),V3(1,3);
+  J.Vector(1,V1);
+  V1.Multiply(Passage,V1);
+  V1.Multiply(J.Value(1));
+  J.Vector(2,V2);
+  V2.Multiply(Passage,V2);
+  V2.Multiply(J.Value(2));
+  J.Vector(3,V3);
+  V3.Multiply(Passage,V3);
+  V3.Multiply(J.Value(3));
+
+  inertia = gp_Mat (gp_XYZ(V1(1),V2(1),V3(1)),
+		    gp_XYZ(V1(2),V2(2),V3(2)),
+		    gp_XYZ(V1(3),V2(3),V3(3)));
+
+  gp_Mat Hop;
+  GProp::HOperator(g,loc,dim,Hop);
+  inertia = inertia+Hop;
+}
+
+
+void GProp_CelGProps::Perform(const gp_Lin& C, 			      
+			      const Standard_Real U1,
+			      const Standard_Real U2)
+{
+  gp_Ax1 Pos = C.Position();
+  gp_Pnt P1 = ElCLib::LineValue(U1,Pos);
+  dim =Abs(U2-U1);
+  gp_Pnt P2 = ElCLib::LineValue(U2,Pos);
+  g.SetCoord((P1.X()+P2.X())/2.,(P1.Y()+P2.Y())/2.,(P1.Z()+P2.Z())/2.);
+  Standard_Real Vx,Vy,Vz,X0,Y0,Z0;
+  Pos.Direction().Coord(Vx,Vy,Vz);
+  Pos.Location().Coord(X0,Y0,Z0);
+  Standard_Real alfa1 = (Vz*Vz+Vy*Vy)/3.;
+  Standard_Real alfa2 = Vy*Y0+Vz*Z0;
+  Standard_Real alfa3 = Y0*Y0+Z0*Z0;
+  Standard_Real Ixx =  (U2*(U2*(U2*alfa1+alfa2)+alfa3)) -
+              (U1*(U1*(U1*alfa1+alfa2)+alfa3));
+  alfa1 = (Vz*Vz+Vx*Vx)/3.;
+  alfa2 = Vx*X0+Vz*Z0;
+  alfa3 = X0*X0+Z0*Z0;
+  Standard_Real Iyy =  (U2*(U2*(U2*alfa1+alfa2)+alfa3)) -
+              (U1*(U1*(U1*alfa1+alfa2)+alfa3));
+  alfa1 = (Vy*Vy+Vx*Vx)/3.;
+  alfa2 = Vy*Y0+Vz*Z0;
+  alfa3 = Y0*Y0+Z0*Z0;
+  Standard_Real Izz =  (U2*(U2*(U2*alfa1+alfa2)+alfa3)) -
+                       (U1*(U1*(U1*alfa1+alfa2)+alfa3));
+  alfa1 = (Vy*Vx)/3.;
+  alfa2 = (Vy*X0+Vx*Y0)/2.;
+  alfa3 = Y0*X0;
+  Standard_Real Ixy =  (U2*(U2*(U2*alfa1+alfa2)+alfa3)) -
+                       (U1*(U1*(U1*alfa1+alfa2)+alfa3));
+  alfa1 = (Vz*Vx)/3.;
+  alfa2 = (Vz*X0+Vx*Z0)/2;
+  alfa3 = Z0*X0;
+  Standard_Real Ixz =  (U2*(U2*(U2*alfa1+alfa2)+alfa3)) -
+                       (U1*(U1*(U1*alfa1+alfa2)+alfa3));
+  alfa1 = (Vy*Vz)/3.;
+  alfa2 = (Vy*Z0+Vz*Y0)/2.;
+  alfa3 = Y0*Z0;
+  Standard_Real Iyz =  (U2*(U2*(U2*alfa1+alfa2)+alfa3)) -
+                       (U1*(U1*(U1*alfa1+alfa2)+alfa3));
+
+  inertia = gp_Mat (gp_XYZ (Ixx, -Ixy,-Ixz),
+                    gp_XYZ (-Ixy, Iyy,-Iyz),
+                    gp_XYZ (-Ixz,-Iyz, Izz));
+}
+
+
+GProp_CelGProps::GProp_CelGProps (const gp_Circ& C, const gp_Pnt& CLocation) 
+{
+  SetLocation(CLocation);
+  Perform(C,0.,2.*PI);
+}
+
+GProp_CelGProps::GProp_CelGProps (const gp_Circ& C, 
+				  const Standard_Real U1, 
+				  const Standard_Real U2,
+				  const gp_Pnt& CLocation)
+{
+  SetLocation(CLocation);
+  Perform(C,U1,U2);
+}
+
+GProp_CelGProps::GProp_CelGProps (const gp_Lin& C, 
+				  const Standard_Real U1, 
+				  const Standard_Real U2,
+				  const gp_Pnt& CLocation)
+{
+  SetLocation(CLocation);
+  Perform(C,U1,U2);
+}
+

src/GProp/GProp_CurveTool.cdl

@@ -0,0 +1,61 @@
+-- File:	 CurveTool.cdl
+-- Created:	 Wed Aug 26 15:32:21 1992
+-- Author:       Jean-Claude Vauthier
+---Copyright:	 Matra Datavision 1992
+
+
+
+
+
+deferred generic class CurveTool from GProp (Curve as any)   
+
+    --- Purpose :
+    --  This template defines the minimum of methods required
+    --  to compute the global properties of a C1 parametric
+    --  curve in 3d space with the algorithmes of package GProp.
+    --  To compute the global properties of your curves, you 
+    --  have to define your own "CurveTool" using this template.  
+    --
+    --  Curve must be a bounded curve of continuity C1 defined in 3d
+    --  space.
+
+uses Pnt from gp,
+     Vec from gp
+
+is
+
+
+  FirstParameter (myclass; C : Curve)   returns Real;
+        --- Purpose :
+        --  Returns the parametric value of the start point of
+        --  the curve.  The curve is oriented from the start point
+        --  to the end point.
+
+
+  LastParameter (myclass; C : Curve)   returns Real;
+        --- Purpose :
+        --  Returns the parametric value of the end point of
+        --  the curve.  The curve is oriented from the start point
+        --  to the end point.
+
+
+  IntegrationOrder (myclass; C : Curve)    returns Integer;
+        --- Purpose :
+        --  Returns the number of Gauss points required to do
+        --  the integration with a good accuracy using the
+        --  Gauss method.  For a polynomial curve of degree n
+        --  the maxima of accuracy is obtained with an order
+        --  of integration equal to 2*n-1.
+
+
+  Value (myclass; C : Curve; U : Real)  returns Pnt;
+    	--- Purpose : Returns the point of parameter U on the loaded curve.
+
+
+  D1 (myclass; C : Curve; U: Real; P: out Pnt; V1: out Vec);
+    	--- Purpose : 
+    	--  Returns the point of parameter U and the first derivative
+    	--  at this point.
+
+
+end CurveTool;

src/GProp/GProp_DomainTool.cdl

@@ -0,0 +1,29 @@
+-- File:	GProp_DomainTool.cdl
+-- Created:	Fri Nov 27 16:44:12 1992
+-- Author:	Isabelle GRIGNON
+--		<isg@sdsun2>
+---Copyright:	 Matra Datavision 1992
+
+
+deferred generic class DomainTool from GProp (Arc as any)
+
+	---Purpose: Arc iterator
+
+
+is
+  
+  Init(me : in out); 
+        ---Purpose: Initializes the exploration with the parameters already set.
+
+  More(me : in out)  returns Boolean from Standard;
+        --- Purpose :
+        --  Returns True if there is another arc of curve in the list.
+
+  Value(me : in out) returns Arc ;
+  
+  Next(me : in out) ; 
+        --- Purpose :
+        --  Sets the index of the arc iterator to the next arc of
+        --  curve.
+
+end DomainTool;

src/GProp/GProp_FaceTool.cdl

@@ -0,0 +1,96 @@
+-- File:	GProp_FaceTool.cdl
+-- Created:	Fri Nov 27 16:29:37 1992
+-- Author:	Isabelle GRIGNON
+--		<isg@sdsun2>
+---Copyright:	 Matra Datavision 1992
+
+deferred generic class FaceTool  from GProp( Arc as any )
+
+        ---Purpose: This template class defines the minimum of methods required 
+        --          to compute the global properties of Faces with the 
+        --          algorithms of the package GProp. 
+        --          To compute the global properties of a Face, in is necessary
+        --          to define own "FaceTool" and to implement all the methods 
+        --          defined in this template. Note that it is not necessary to 
+        --          inherit this template class.
+	
+uses Pnt           from gp,
+     Pnt2d         from gp,
+     Vec           from gp,
+     Vec2d         from gp, 
+     IsoType       from GeomAbs, 
+     HArray1OfReal from TColStd
+
+is 
+
+  UIntegrationOrder (me) returns Integer;
+        ---Purpose: Returns the number of points required to do the 
+        --          integration in the U parametric direction.
+
+  Bounds (me; U1, U2, V1, V2 : out Real);
+	---Purpose: Returns the parametric bounds of the Face <S>.
+
+  Normal (me;  U, V : Real; P : out Pnt; VNor: out Vec);
+	---Purpose: Computes the point of parameter U, V on the Face <S> and 
+	--          the normal to the face at this point.
+
+  Load(me:in out; A : Arc); 
+        ---Purpose: Loading the boundary arc.
+  
+  Load(me : in out; IsFirstParam: Boolean from Standard; 
+                    theIsoType  : IsoType from GeomAbs);
+        ---Purpose: Loading the boundary arc. This arc is either a top, bottom, 
+        --          left or right bound of a UV rectangle in which the 
+        --          parameters of surface are defined. 
+	--          If IsFirstParam is equal to Standard_True, the face is 
+        --          initialized by either left of bottom bound. Otherwise it is 
+        --          initialized by the top or right one. 
+	--          If theIsoType is equal to GeomAbs_IsoU, the face is 
+        --          initialized with either left or right bound. Otherwise - 
+        --          with either top or bottom one.
+  
+  FirstParameter (me) returns Real ;
+        ---Purpose: Returns the parametric value of the start point of
+        --          the current arc of curve. 
+
+  LastParameter (me) returns Real ;
+        ---Purpose: Returns the parametric value of the end point of
+        --          the current arc of curve. 
+
+  IntegrationOrder (me) returns Integer;
+        ---Purpose: Returns the number of points required to do the 
+        --          integration along the parameter of curve.
+
+  D12d (me; U: Real; P: out Pnt2d; V1: out Vec2d);
+    	---Purpose: Returns the point of parameter U and the first derivative 
+    	--          at this point of a boundary curve.
+
+  GetUKnots(me; theUMin  :        Real          from Standard; 
+                theUMax  :        Real          from Standard; 
+		theUKnots: in out HArray1OfReal from TColStd); 
+	---Purpose: Returns an array of U knots of the face. The first and last 
+        --          elements of the array will be theUMin and theUMax. The 
+        --          middle elements will be the U Knots of the face greater 
+        --          then theUMin and lower then theUMax in increasing order. 
+
+  GetTKnots(me; theTMin  :        Real          from Standard; 
+                theTMax  :        Real          from Standard; 
+		theTKnots: in out HArray1OfReal from TColStd);
+	---Purpose: Returns an array of combination of T knots of the arc and 
+        --          V knots of the face. The first and last elements of the 
+        --          array will be theTMin and theTMax. The middle elements will 
+        --          be the Knots of the arc and the values of parameters of 
+        --          arc on which the value points have V coordinates close to V 
+        --          knots of face. All the parameter will be greater then 
+        --          theTMin and lower then theTMax in increasing order.
+
+end FaceTool;
+
+
+
+
+
+
+
+
+

src/GProp/GProp_GProps.cdl

@@ -0,0 +1,271 @@
+-- File:	GProps.cdl<3>
+-- Created:	Mon Aug 24 18:10:07 1992
+-- Author:	Michel CHAUVAT
+--              JCV January 1992
+---Copyright:	 Matra Datavision 1992
+
+
+class GProps from GProp 
+
+	--- Purpose : 
+	--  Implements a general mechanism to compute the global properties of
+	--  a "compound geometric system" in 3d space    by composition of the
+	--  global properties of "elementary geometric entities"       such as 
+	--  (curve, surface, solid, set of points).  It is possible to compose
+	--  the properties of several "compound geometric systems" too.
+	-- 
+	--  To computes the global properties of a compound geometric
+	--  system you should : 
+        --  . declare the GProps using a constructor which initializes the
+        --    GProps and defines the location point used to compute the inertia
+        --  . compose the global properties of your geometric components with
+        --    the properties of your system using the method Add.
+	--  
+	--  To compute the global properties of the geometric components of
+	--  the system you should  use the services of the following classes :
+	--   - class PGProps for a set of points,
+	--   - class CGProps for a curve,
+	--   - class SGProps for a surface,
+	--   - class VGProps for a "solid".
+	--  The classes CGProps, SGProps, VGProps are generic classes and
+	--  must be instantiated for your application.
+	--  
+	--    
+	--  The global properties computed are :
+	--  - the dimension (length, area or volume)
+	--  - the mass,
+	--  - the centre of mass,
+	--  - the moments of inertia (static moments and quadratic moments),
+	--  - the moment about an axis,
+	--  - the radius of gyration about an axis,
+	--  - the principal properties of inertia  :
+	--    (sea also class PrincipalProps)
+	--    . the principal moments,
+	--    . the principal axis of inertia,
+	--    . the principal radius of gyration,
+        --
+        --         
+        --         
+        --  Example of utilisation in a simplified C++ implementation :
+        --             
+        --    //declares the GProps, the point (0.0, 0.0, 0.0) of the
+        --    //absolute cartesian coordinate system is used as 
+        --    //default reference point to compute the centre of mass
+        --    GProp_GProps System ();
+        --    
+        --    //computes the inertia of a 3d curve
+        --    Your_CGProps Component1 (curve, ....);
+        --    
+        --    //computes the inertia of surfaces
+        --    Your_SGprops Component2 (surface1, ....);        
+        --    Your_SGprops Component3 (surface2,....);        
+        --    
+        --    //composes the global properties of components 1, 2, 3
+        --    //a density can be associated with the components, the
+        --    //density can be defaulted to 1.
+        --    Real Density1 = 2.0;        
+        --    Real Density2 = 3.0;        
+        --    System.Add (Component1, Density1);
+        --    System.Add (Component2, Density2);
+        --    System.Add (Component3);
+        --    
+        --    //returns the centre of mass of the system in the
+        --    //absolute cartesian coordinate system
+        --    gp_Pnt G = System.CentreOfMass ();
+        --    
+        --    //computes the principales inertia of the system
+        --    GProp_PrincipalProps Pp  = System.PrincipalProperties();
+        --    
+        --    //returns the principal moments and radius of gyration
+        --    Real Ixx, Iyy, Izz, Rxx, Ryy, Rzz;
+        --    Pp.Moments (Ixx, Iyy, Izz);
+        --    Pp.RadiusOfGyration (Ixx, Iyy, Izz);
+        --    
+        --    
+
+uses    Ax1 from gp,
+        Mat from gp,
+        Pnt from gp,
+	PrincipalProps from GProp
+
+raises DomainError from Standard
+
+is
+
+  Create  returns GProps;
+        --- Purpose :
+        --  The origin (0, 0, 0) of the absolute cartesian coordinate system
+        --  is used to compute the global properties.
+
+  
+  Create (SystemLocation : Pnt)  returns GProps;
+        --- Purpose :
+        --  The point SystemLocation is used to compute the gobal properties
+        --  of the system. For more accuracy it is better to define this 
+        --  point closed to the location of the system. For example it could
+        --  be a point around the centre of mass of the system.
+--  This point is referred to as the reference point for
+-- this framework. For greater accuracy it is better for
+-- the reference point to be close to the location of the
+-- system. It can, for example, be a point near the
+-- center of mass of the system.
+-- At initialization, the framework is empty; i.e. it
+-- retains no dimensional information such as mass, or
+-- inertia. However, it is now able to bring together
+-- global properties of various other systems, whose
+-- global properties have already been computed
+-- using another framework. To do this, use the
+-- function Add to define the components of the
+-- system. Use it once per component of the system,
+-- and then use the interrogation functions available to
+-- access the computed values.
+       
+
+  Add (me : in out; Item : GProps; Density : Real =1.0)
+    	--- Purpose : Either 
+-- - initializes the global properties retained by this
+--   framework from those retained by the framework Item, or
+-- - brings together the global properties still retained by
+--   this framework with those retained by the framework Item.
+--   The value Density, which is 1.0 by default, is used as
+-- the density of the system analysed by Item.
+-- Sometimes the density will have already been given at
+-- the time of construction of the framework Item. This
+-- may be the case for example, if Item is a
+-- GProp_PGProps framework built to compute the
+-- global properties of a set of points ; or another
+-- GProp_GProps object which already retains
+-- composite global properties. In these cases the real
+-- density was perhaps already taken into account at the
+-- time of construction of Item. Note that this is not
+-- checked: if the density of parts of the system is taken
+-- into account two or more times, results of the
+-- computation will be false.
+-- Notes :
+-- - The point relative to which the inertia of Item is
+--   computed (i.e. the reference point of Item) may be
+--   different from the reference point in this
+--   framework. Huygens' theorem is applied
+--   automatically to transfer inertia values to the
+ --  reference point in this framework.
+-- - The function Add is used once per component of
+--   the system. After that, you use the interrogation
+--   functions available to access values computed for the system.
+-- - The system whose global properties are already
+--   brought together by this framework is referred to
+--   as the current system. However, the current system
+--   is not retained by this framework, which maintains
+--   only its global properties.
+--      Exceptions
+-- Standard_DomainError if Density is less than or
+-- equal to gp::Resolution(). 
+     raises DomainError
+     is static;
+
+
+
+  Mass (me)	returns Real   is static;
+    	--- Purpose: Returns the mass of the current system.
+-- If no density is attached to the components of the
+-- current system the returned value corresponds to :
+-- - the total length of the edges of the current
+--   system if this framework retains only linear
+--   properties, as is the case for example, when
+--   using only the LinearProperties function to
+--   combine properties of lines from shapes, or
+-- - the total area of the faces of the current system if
+--   this framework retains only surface properties,
+--   as is the case for example, when using only the
+--   SurfaceProperties function to combine
+--   properties of surfaces from shapes, or
+-- - the total volume of the solids of the current
+--   system if this framework retains only volume
+--   properties, as is the case for example, when
+--   using only the VolumeProperties function to
+--   combine properties of volumes from solids.
+--   Warning
+-- A length, an area, or a volume is computed in the
+-- current data unit system. The mass of a single
+-- object is obtained by multiplying its length, its area
+-- or its volume by the given density. You must be
+-- consistent with respect to the units used. 
+
+
+  CentreOfMass (me)  returns Pnt  is static;
+    	--- Purpose :
+    	--  Returns the center of mass of the current system. If
+-- the gravitational field is uniform, it is the center of gravity.
+-- The coordinates returned for the center of mass are
+-- expressed in the absolute Cartesian coordinate system.
+	
+  MatrixOfInertia (me)   returns Mat   is static;
+    	--- Purpose: 
+    	--  returns the matrix of inertia. It is a symmetrical matrix. 
+    	--  The coefficients of the matrix are the quadratic moments of 
+    	--  inertia.
+        --  
+        --               | Ixx  Ixy  Ixz | 
+        --   matrix =    | Ixy  Iyy  Iyz |
+        --               | Ixz  Iyz  Izz |
+        --              
+        --  The moments of inertia are denoted by Ixx, Iyy, Izz.
+        --  The products of inertia are denoted by Ixy, Ixz, Iyz.
+        --  The matrix of inertia is returned in the central coordinate
+        --  system (G, Gx, Gy, Gz) where G is the centre of mass of the
+        --  system and Gx, Gy, Gz the directions parallel to the X(1,0,0)
+        --  Y(0,1,0) Z(0,0,1) directions of the absolute cartesian
+        --  coordinate system. It is possible to compute the matrix of
+        --  inertia at another location point using the Huyghens theorem
+        --  (you can use the method of package GProp : HOperator).
+
+
+  StaticMoments (me; Ix, Iy, Iz : out Real) is static;
+    	--- Purpose : Returns Ix, Iy, Iz, the static moments of inertia of the
+-- current system; i.e. the moments of inertia about the
+-- three axes of the Cartesian coordinate system.
+    	
+
+
+  MomentOfInertia (me; A : Ax1)  returns Real is static;
+    	--- Purpose : 
+    	--  computes the moment of inertia of the material system about the
+    	--  axis A.
+
+
+  PrincipalProperties (me) returns PrincipalProps is static;
+        --- Purpose : Computes the principal properties of inertia of the current system.
+-- There is always a set of axes for which the products
+-- of inertia of a geometric system are equal to 0; i.e. the
+-- matrix of inertia of the system is diagonal. These axes
+-- are the principal axes of inertia. Their origin is
+-- coincident with the center of mass of the system. The
+-- associated moments are called the principal moments of inertia.
+-- This function computes the eigen values and the
+-- eigen vectors of the matrix of inertia of the system.
+-- Results are stored by using a presentation framework
+-- of principal properties of inertia
+-- (GProp_PrincipalProps object) which may be
+-- queried to access the value sought.
+      
+
+
+  RadiusOfGyration (me; A : Ax1)  returns Real is static;
+    	--- Purpose : Returns the radius of gyration of the current system about the axis A.
+
+
+
+
+fields
+
+    g       : Pnt is protected;
+        --- Purpose : centre of mass
+    loc     : Pnt is protected;
+        --- Purpose : location point used to compute the inertia
+    dim     : Real is protected;
+        --- Purpose : mass or length or area or volume of the GProps
+    inertia : Mat is protected;    
+        --- Purpose : matrix of inertia
+    	
+end GProps;
+
+

src/GProp/GProp_GProps.cxx

@@ -0,0 +1,178 @@
+#include <GProp_GProps.ixx>
+#include <GProp.hxx>
+#include <math_Jacobi.hxx>
+#include <gp.hxx>
+#include <gp.hxx>
+#include <gp_XYZ.hxx>
+#include <gp_Vec.hxx>
+
+
+GProp_GProps::GProp_GProps () : g (gp::Origin()) , loc (gp::Origin()), dim (0.0) 
+{ 
+  inertia = gp_Mat(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
+}
+
+
+GProp_GProps::GProp_GProps (const gp_Pnt& SystemLocation) : 
+       g (gp::Origin()), loc (SystemLocation), dim (0.0)
+{ 
+  inertia = gp_Mat(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
+}
+
+
+
+void GProp_GProps::Add (const GProp_GProps& Item, const Standard_Real Density) {
+  if (Density <= gp::Resolution())  Standard_DomainError::Raise();
+  if (loc.Distance (Item.loc) <= gp::Resolution ()) {
+    gp_XYZ GXYZ = (Item.g.XYZ()).Multiplied (Item.dim * Density);
+    g.SetXYZ (g.XYZ().Multiplied (dim));
+    GXYZ.Add (g.XYZ());
+    dim = dim + Item.dim * Density;
+    if (Abs(dim) >= 1.e-20 ) {
+      GXYZ.Divide (dim);
+      g.SetXYZ (GXYZ);
+    }
+    else {
+      g.SetCoord(0., 0., 0.);
+    }
+    inertia = inertia + Item.inertia * Density;
+  }else{
+    gp_XYZ Itemloc = loc.XYZ() - Item.loc.XYZ();
+    gp_XYZ Itemg   = Item.loc.XYZ() + Item.g.XYZ();
+    gp_XYZ GXYZ    = Item.g.XYZ() - Itemloc;
+    GXYZ    = GXYZ.Multiplied (Item.dim * Density);
+    g.SetXYZ (g.XYZ().Multiplied (dim));
+    GXYZ.Add (g.XYZ());
+    dim = dim + Item.dim * Density;
+    if (Abs(dim) >= 1.e-20 ) {
+      GXYZ.Divide (dim);
+      g.SetXYZ (GXYZ);
+    }
+    else {
+      g.SetCoord(0., 0., 0.);
+    }
+    //We have to compute the inertia of the Item at the location point
+    //of the system using the Huyghens theorem
+    gp_Mat HMat;
+    gp_Mat ItemInertia = Item.inertia;
+    if (Item.g.XYZ().Modulus() > gp::Resolution()) {
+      //Computes the inertia of Item at its dim centre
+      GProp::HOperator (Itemg, Item.loc, Item.dim, HMat);
+      ItemInertia = ItemInertia - HMat;
+    }
+    //Computes the inertia of Item at the location point of the system
+    GProp::HOperator (Itemg, loc, Item.dim, HMat);
+    ItemInertia = ItemInertia + HMat;
+    inertia = inertia + ItemInertia * Density;
+  }
+}
+
+Standard_Real GProp_GProps::Mass () const { return dim;  }
+
+gp_Pnt GProp_GProps::CentreOfMass () const 
+{ 
+  return gp_Pnt(loc.XYZ() + g.XYZ()); 
+}
+
+gp_Mat GProp_GProps::MatrixOfInertia () const { 
+  gp_Mat HMat;
+  GProp::HOperator (g,gp::Origin(),dim, HMat);
+  return inertia - HMat;
+}
+
+
+
+void GProp_GProps::StaticMoments (Standard_Real& Ix, 
+				  Standard_Real& Iy, 
+				  Standard_Real& Iz) const {
+
+  gp_XYZ G = loc.XYZ() + g.XYZ();
+  Ix = G.X() * dim;
+  Iy = G.Y() * dim;
+  Iz = G.Z() * dim;
+}
+
+
+
+
+Standard_Real GProp_GProps::MomentOfInertia (const gp_Ax1& A) const {
+  // Moment of inertia / axis A
+  // 1] computes the math_Matrix of inertia / A.location()
+  // 2] applies this math_Matrix to A.Direction()
+  // 3] then computes the scalar product between this vector and 
+  //    A.Direction()
+
+      
+  if (loc.Distance (A.Location()) <= gp::Resolution()) {
+    return (A.Direction().XYZ()).Dot (
+           (A.Direction().XYZ()).Multiplied (inertia));
+  }
+  else {
+    gp_Mat HMat;
+    gp_Mat AxisInertia = MatrixOfInertia();
+    GProp::HOperator (gp_Pnt (loc.XYZ() + g.XYZ()), A.Location(), dim, HMat);
+    AxisInertia = AxisInertia + HMat;
+    return (A.Direction().XYZ()).Dot (
+           (A.Direction().XYZ()).Multiplied (AxisInertia));
+  }
+}
+
+
+
+Standard_Real GProp_GProps::RadiusOfGyration (const gp_Ax1& A) const {
+
+  return Sqrt (MomentOfInertia (A) / dim);
+}
+
+
+
+GProp_PrincipalProps GProp_GProps::PrincipalProperties () const {
+  
+  math_Matrix DiagMat (1, 3, 1, 3);
+  Standard_Integer i, j;
+  gp_Mat AxisInertia = MatrixOfInertia();
+  for (j = 1; j <= 3; j++) {
+    for (i = 1; i <= 3; i++) {
+       DiagMat (i, j) = AxisInertia.Value (i, j);
+    }
+  }
+  math_Jacobi J (DiagMat);
+  Standard_Real Ixx = J.Value (1);
+  Standard_Real Iyy = J.Value (2);
+  Standard_Real Izz = J.Value (3);
+  DiagMat  = J.Vectors ();      
+  gp_Vec Vxx (DiagMat (1,1), DiagMat (2, 1), DiagMat (3, 1));
+  gp_Vec Vyy (DiagMat (1,2), DiagMat (2, 2), DiagMat (3, 2));
+  gp_Vec Vzz (DiagMat (1,3), DiagMat (2, 3), DiagMat (3, 3));
+  //
+  // protection contre dim == 0.0e0 au cas ou on aurait rentre qu'un point
+  //
+  Standard_Real Rxx = 0.0e0 ;
+  Standard_Real Ryy = 0.0e0 ;
+  Standard_Real Rzz = 0.0e0 ;  
+  if (0.0e0 != dim) {
+    Rxx = Sqrt (Abs(Ixx / dim));
+    Ryy = Sqrt (Abs(Iyy / dim));
+    Rzz = Sqrt (Abs(Izz / dim));
+  }
+  return GProp_PrincipalProps (Ixx, Iyy, Izz, Rxx, Ryy, Rzz, Vxx, Vyy, Vzz,
+                             gp_Pnt(g.XYZ() + loc.XYZ()));
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

src/GProp/GProp_PEquation.cdl

@@ -0,0 +1,136 @@
+-- File:	GProp_PEquation.cdl
+-- Created:	Wed Jun 16 11:38:18 1993
+-- Author:	Isabelle GRIGNON
+--		<isg@sdsun2>
+---Copyright:	 Matra Datavision 1993
+
+
+class PEquation from GProp 
+
+	---Purpose: A framework to analyze a collection - or cloud
+    	-- - of points and to verify if they are coincident,
+    	-- collinear or coplanar within a given precision. If
+    	-- so, it also computes the mean point, the mean
+    	-- line or the mean plane of the points. If not, it
+    	-- computes the minimal box which includes all the points. 
+
+
+uses Pnt from gp,
+     Pln from gp,
+     Lin from gp,
+     Vec from gp,
+     EquaType from GProp,
+     Array1OfPnt from TColgp
+
+
+raises NoSuchObject from Standard
+
+is
+     
+    Create(Pnts: Array1OfPnt;Tol : Real from Standard) 
+    returns PEquation from GProp;
+    	---Purpose: Constructs a framework to analyze the
+    	-- collection of points Pnts and computes:
+    	-- -   the mean point if the points in question are
+    	--   considered to be coincident within the precision Tol, or
+    	-- -   the mean line if they are considered to be
+    	--   collinear within the precision Tol, or
+    	-- -   the mean plane if they are considered to be
+    	--   coplanar within the precision Tol, or
+    	-- -   the minimal box which contains all the points. Use :
+    	-- -   the functions IsPoint, IsLinear, IsPlanar
+    	--   and IsSpace to find the result of the analysis, and
+    	-- -   the function Point, Line, Plane or Box to
+    	--   access the computed result.
+    
+    IsPlanar(me) returns Boolean
+    	---Purpose: Returns true if, according to the given
+    	-- tolerance, the points analyzed by this framework are  coplanar.
+    	-- Use the function  Plane  to access the computed result.
+    is static;
+    
+    IsLinear(me) returns Boolean
+    	---Purpose:  Returns true if, according to the given
+    	-- tolerance, the points analyzed by this framework are  colinear.
+    	-- Use the function  Line  to access the computed result.
+    is static;
+    
+    IsPoint(me) returns Boolean
+    	---Purpose:  Returns true if, according to the given
+    	-- tolerance, the points analyzed by this framework are  coincident.
+    	--       Use the function  Point  to access the computed result.
+    is static;
+    
+    IsSpace(me) returns Boolean
+    	---Purpose: Returns true if, according to the given
+    	-- tolerance value, the points analyzed by this
+    	-- framework are neither coincident, nor collinear, nor coplanar.
+    	-- Use the function Box to query the smallest box
+    	-- that includes the collection of points.
+    is static;
+    
+    Plane(me) returns Pln from gp
+    raises NoSuchObject
+	---Purpose: Returns the mean plane passing near all the
+    	-- points analyzed by this framework if, according
+    	-- to the given precision, the points are considered to be coplanar.
+    	-- Exceptions
+    	-- Standard_NoSuchObject if, according to the
+    	-- given precision value, the points analyzed by
+    	-- this framework are considered to be:
+    	-- -   coincident, or
+    	-- -   collinear, or
+    	-- -   not coplanar.
+    is static;
+    
+    Line(me) returns Lin from gp
+    raises NoSuchObject
+	---Purpose: Returns the mean line passing near all the
+    	-- points analyzed by this framework if, according
+    	-- to the given precision value, the points are considered to be collinear.
+    	-- Exceptions
+    	-- Standard_NoSuchObject if, according to the
+    	-- given precision, the points analyzed by this
+    	-- framework are considered to be:
+    	-- -   coincident, or
+    	-- -   not collinear.
+    is static;
+    
+    Point(me) returns Pnt from gp
+    raises NoSuchObject
+	---Purpose: Returns the mean point of all the points
+    	-- analyzed by this framework if, according to the
+    	-- given precision, the points are considered to be coincident.
+    	-- Exceptions
+    	-- Standard_NoSuchObject if, according to the
+    	-- given precision, the points analyzed by this
+    	-- framework are not considered to be coincident.
+    is static;
+
+    Box(me; P : out Pnt from gp; V1,V2,V3 : out Vec from gp)
+    	---Purpose: Returns the definition of the smallest box which
+    	-- contains all the points analyzed by this
+    	-- framework if, according to the given precision
+    	-- value, the points are considered to be neither
+    	-- coincident, nor collinear and nor coplanar.
+    	-- This box is centered on the barycenter P of the
+    	-- collection of points. Its sides are parallel to the
+    	-- three vectors V1, V2 and V3, the length of
+    	-- which is the length of the box in the corresponding direction.
+    	-- Note: Vectors V1, V2 and V3 are parallel to
+    	-- the three axes of principal inertia of the system
+    	-- composed of the collection of points where each point is of equal mass.
+    	-- Exceptions
+    	-- Standard_NoSuchObject if, according to the given precision,
+    	-- the points analyzed by this framework are considered to be coincident, collinear or coplanar.
+    is static;
+
+fields
+
+type      : EquaType from GProp;
+g         : Pnt from gp;
+v1        : Vec from gp;
+v2        : Vec from gp;
+v3        : Vec from gp;
+end PEquation;
+

src/GProp/GProp_PEquation.cxx

@@ -0,0 +1,133 @@
+#include <GProp_PEquation.ixx>
+#include <GProp_PrincipalProps.hxx>
+#include <GProp_PGProps.hxx>
+
+GProp_PEquation::GProp_PEquation(const TColgp_Array1OfPnt& Pnts, 
+				       const Standard_Real Tol) 
+{
+  GProp_PGProps Pmat(Pnts);
+  g = Pmat.CentreOfMass(); 
+  Standard_Real Xg,Yg,Zg;
+  g.Coord(Xg,Yg,Zg);
+  GProp_PrincipalProps Pp = Pmat.PrincipalProperties();
+  gp_Vec V1 = Pp.FirstAxisOfInertia();
+  Standard_Real Xv1,Yv1,Zv1;
+  V1.Coord(Xv1,Yv1,Zv1); 
+  gp_Vec V2 = Pp.SecondAxisOfInertia(); 
+  Standard_Real Xv2,Yv2,Zv2;
+  V2.Coord(Xv2,Yv2,Zv2);
+  gp_Vec V3 = Pp.ThirdAxisOfInertia(); 
+  Standard_Real Xv3,Yv3,Zv3;
+  V3.Coord(Xv3,Yv3,Zv3);
+  Standard_Real D,X,Y,Z;
+  Standard_Real Dmx1 = RealFirst();
+  Standard_Real Dmn1 = RealLast();
+  Standard_Real Dmx2 = RealFirst();
+  Standard_Real Dmn2 = RealLast();
+  Standard_Real Dmx3 = RealFirst();
+  Standard_Real Dmn3 = RealLast();
+
+  for (Standard_Integer i = Pnts.Lower(); i <= Pnts.Upper();i++){
+    Pnts(i).Coord(X,Y,Z);
+    D = (X-Xg)*Xv1 +(Y-Yg)*Yv1 + (Z-Zg)*Zv1;
+    if (D > Dmx1) Dmx1 = D;
+    if (D < Dmn1) Dmn1 = D;
+    D = (X-Xg)*Xv2 +(Y-Yg)*Yv2 + (Z-Zg)*Zv2;
+    if (D > Dmx2) Dmx2 = D;
+    if (D < Dmn2) Dmn2 = D;
+    D = (X-Xg)*Xv3 +(Y-Yg)*Yv3 + (Z-Zg)*Zv3;
+    if (D > Dmx3) Dmx3 = D;
+    if (D < Dmn3) Dmn3 = D;
+  }
+  Standard_Integer dimension= 3 ;
+  Standard_Integer It = 0;
+  if (Abs(Dmx1-Dmn1) <= Tol)  {
+    dimension =dimension-1;
+    It =1;
+  }
+  if (Abs(Dmx2-Dmn2) <= Tol)  {
+    dimension =dimension-1;
+    It =2*(It+1);
+  }
+  if (Abs(Dmx3-Dmn3) <= Tol)  {
+    dimension =dimension-1;
+    It = 3*(It+1);
+  }
+  switch (dimension)  {
+  case 0:
+    {
+      type = GProp_Point;
+      break;
+    }
+  case 1:
+    {
+      type = GProp_Line;
+      if (It == 4) v1 = V3;
+      else if (It == 6) v1 = V2;
+      else v1 = V1; 
+      break;
+    }
+  case 2:
+    {
+      type = GProp_Plane;
+      if (It == 1) v1 = V1;
+      else if (It == 2) v1 =V2;
+      else v1 = V3;
+      break;
+    }
+  case 3:
+    {
+      type = GProp_Space;
+      g.SetXYZ(g.XYZ() + Dmn1*V1.XYZ() + Dmn2*V2.XYZ() + Dmn3*V3.XYZ());
+      v1 = (Dmx1-Dmn1)*V1;      
+      v2 = (Dmx2-Dmn2)*V2;
+      v3 = (Dmx3-Dmn3)*V3;
+      break;
+    }
+  }
+}
+Standard_Boolean  GProp_PEquation::IsPlanar() const {
+
+  if (type == GProp_Plane) return Standard_True;
+  else  return Standard_False;
+}
+
+Standard_Boolean  GProp_PEquation::IsLinear() const {
+
+  if (type == GProp_Line) return Standard_True;
+  else  return Standard_False;
+}
+
+Standard_Boolean  GProp_PEquation::IsPoint() const {
+
+  if (type == GProp_Point) return Standard_True;
+  else  return Standard_False;
+}
+
+Standard_Boolean GProp_PEquation::IsSpace() const {
+  if (type == GProp_Space) return Standard_True;
+  else  return Standard_False;
+}
+
+gp_Pln  GProp_PEquation::Plane() const {
+  if (!IsPlanar()) Standard_NoSuchObject::Raise();
+  return gp_Pln(g,v1);
+}
+gp_Lin  GProp_PEquation::Line() const {
+  if (!IsLinear()) Standard_NoSuchObject::Raise();
+  return gp_Lin(g,gp_Dir(v1));  
+}
+
+gp_Pnt  GProp_PEquation::Point() const {
+  if (!IsPoint()) Standard_NoSuchObject::Raise();
+  return g;
+}
+
+void GProp_PEquation::Box(gp_Pnt& P  , gp_Vec& V1,
+			     gp_Vec& V2 , gp_Vec& V3) const {
+  if (!IsSpace()) Standard_NoSuchObject::Raise();
+  P = g;
+  V1 = v1;
+  V2 = v2;
+  V3 = v3;
+}

src/GProp/GProp_PGProps.cdl

@@ -0,0 +1,176 @@
+-- File:	PGProps.cdl
+-- Created:	Fri Feb 14 07:30:21 1992
+-- Author:	Jean Claude VAUTHIER
+---Copyright:	 Matra Datavision 1992
+
+
+
+
+
+
+class PGProps
+
+from GProp
+
+inherits GProps
+---Purpose: A framework for computing the global properties of a
+-- set of points.
+-- A point mass is attached to each point. The global
+-- mass of the system is the sum of each individual
+-- mass. By default, the point mass is equal to 1 and the
+-- mass of a system composed of N points is equal to N.
+-- Warning
+-- A framework of this sort provides functions to handle
+-- sets of points easily. But, like any GProp_GProps
+-- object, by using the Add function, it can theoretically
+-- bring together the computed global properties and
+-- those of a system more complex than a set of points .
+-- The mass of each point and the density of each
+-- component of the composed system must be
+-- coherent. Note that this coherence cannot be checked.
+-- Nonetheless, you are advised to restrict your use of a
+-- GProp_PGProps object to a set of points and to
+-- create a GProp_GProps object in order to bring
+-- together global properties of different systems.
+        
+uses Pnt          from gp,
+     Array1OfPnt  from TColgp,
+     Array2OfPnt  from TColgp,
+     Array1OfReal from TColStd,
+     Array2OfReal from TColStd
+
+     
+raises DimensionError from Standard,
+       DomainError    from Standard
+     
+is
+
+  Create   returns PGProps;
+        --- Purpose : Initializes a framework to compute global properties
+-- on a set of points.
+-- The point relative to which the inertia of the system is
+-- computed will be the origin (0, 0, 0) of the
+-- absolute Cartesian coordinate system.
+-- At initialization, the framework is empty, i.e. it retains
+-- no dimensional information such as mass and inertia.
+-- It is, however, now able to keep global properties of a
+-- set of points while new points are added using the
+-- AddPoint function.
+-- The set of points whose global properties are brought
+-- together by this framework will then be referred to as
+-- the current system. The current system is, however,
+-- not kept by this framework, which only keeps that
+-- system's global properties. Note that the current
+-- system may be more complex than a set of points.
+       
+
+  AddPoint (me : in out; P : Pnt)
+        --- Purpose : Brings together the global properties already
+-- retained by this framework with those induced by
+-- the point Pnt. Pnt may be the first point of the current system.
+-- A point mass is attached to the point Pnt, it is either
+-- equal to 1. or to Density.
+            is static;
+
+
+  AddPoint (me : in out; P : Pnt; Density : Real)
+        --- Purpose :
+        --  Adds a new point P with its density in the system of points
+   -- Exceptions
+-- Standard_DomainError if the mass value Density
+-- is less than gp::Resolution().
+      raises DomainError
+     is static;
+     
+     
+  Create (Pnts : Array1OfPnt) returns PGProps;
+        --- Purpose :
+        --  computes the global properties of the system of points Pnts.
+        --  The density of the points are defaulted to all being 1
+
+
+  Create (Pnts : Array2OfPnt) returns PGProps;
+        --- Purpose :
+        --  computes the global properties of the system of points Pnts.
+        --  The density of the points are defaulted to all being 1
+
+
+  Create (Pnts : Array1OfPnt; Density : Array1OfReal) returns PGProps
+        --- Purpose :
+        --  computes the global properties of the system of points Pnts.
+        --  A density is associated with each point.
+     raises DomainError,
+        --- Purpose : 
+        --  raises if a density is lower or equal to Resolution from package
+        --  gp.
+            DimensionError;
+        --- Purpose : 
+        --  raises if the length of Pnts and the length of Density
+        --  is not the same.
+
+
+
+  Create (Pnts : Array2OfPnt; Density : Array2OfReal) returns PGProps
+        --- Purpose :
+        --  computes the global properties of the system of points Pnts.
+        --  A density is associated with each point.
+     raises DomainError,
+        --- Purpose : 
+        --  Raised if a density is lower or equal to Resolution from package
+        --  gp.
+            DimensionError;
+        --- Purpose : 
+        --  Raised if the length of Pnts and the length of Density
+        --  is not the same.
+
+
+
+
+  Barycentre (myclass; Pnts : Array1OfPnt from TColgp)  returns Pnt;
+        --- Purpose : 
+        --  Computes the barycentre of a set of points. The density of the
+        --  points is defaulted to 1.
+
+
+  Barycentre (myclass; Pnts : Array2OfPnt from TColgp)   returns Pnt;
+        --- Purpose : 
+        --  Computes the barycentre of a set of points. The density of the
+        --  points is defaulted to 1.
+
+
+  Barycentre (myclass; Pnts : Array1OfPnt from TColgp; 
+              Density : Array1OfReal from TColStd;
+              Mass : out Real; G : out Pnt)
+        --- Purpose : 
+        --  Computes the barycentre of a set of points. A density is associated
+        --  with each point.
+     raises DomainError,
+        --- Purpose : 
+        --  raises if a density is lower or equal to Resolution from package
+        --  gp.
+            DimensionError;
+        --- Purpose : 
+        --  Raised if the length of Pnts and the length of Density
+        --  is not the same.
+
+
+  Barycentre (myclass; Pnts : Array2OfPnt from TColgp; 
+              Density : Array2OfReal from TColStd;
+              Mass : out Real; G : out Pnt)
+        --- Purpose : 
+        --  Computes the barycentre of a set of points. A density is associated
+        --  with each point.
+     raises DomainError,
+        --- Purpose : 
+        --  Raised if a density is lower or equal to Resolution from package
+        --  gp.
+            DimensionError;
+        --- Purpose : 
+        --  Raised if the length of Pnts and the length of Density
+        --  is not the same.
+
+
+end PGProps;
+
+
+

src/GProp/GProp_PGProps.cxx

@@ -0,0 +1,213 @@
+#include <GProp_PGProps.ixx>
+#include <Standard_DomainError.hxx>
+#include <Standard_DimensionError.hxx>
+
+#include <gp.hxx>
+#include <gp_XYZ.hxx>
+//#include <gp.hxx>
+
+typedef gp_Pnt Pnt;
+typedef gp_Mat Mat;
+typedef gp_XYZ XYZ;
+typedef TColgp_Array1OfPnt          Array1OfPnt;
+typedef TColgp_Array2OfPnt          Array2OfPnt;
+typedef TColStd_Array1OfReal Array1OfReal;
+typedef TColStd_Array2OfReal Array2OfReal;
+
+
+
+GProp_PGProps::GProp_PGProps ()
+{
+  g = gp::Origin();
+  loc = gp::Origin();
+  dim = 0.0;
+}
+
+void GProp_PGProps::AddPoint (const Pnt& P)
+{
+  Standard_Real Xp, Yp, Zp;
+  P.Coord (Xp, Yp, Zp);
+  Standard_Real Ixy = - Xp * Yp;
+  Standard_Real Ixz = - Xp * Zp;
+  Standard_Real Iyz = - Yp * Zp;
+
+  Standard_Real Ixx = Yp * Yp + Zp * Zp;
+  Standard_Real Iyy = Xp * Xp + Zp * Zp;
+  Standard_Real Izz = Xp * Xp + Yp * Yp;
+  Mat Mp (XYZ (Ixx, Ixy, Ixz), XYZ (Ixy, Iyy, Iyz), XYZ (Ixz, Iyz, Izz));
+  if (dim == 0) {
+    dim = 1;
+    g = P;
+    inertia = Mp;
+  }
+  else {
+    Standard_Real X, Y, Z;
+    g.Coord (X, Y, Z);
+    X = X * dim + Xp;
+    Y = Y * dim + Yp;
+    Z = Z * dim + Zp;
+    dim = dim + 1;
+    X /= dim;
+    Y /= dim;
+    Z /= dim;
+    g.SetCoord (X, Y, Z);
+    inertia = inertia + Mp;
+  }      
+}
+
+void GProp_PGProps::AddPoint (const gp_Pnt& P, const Standard_Real Density)
+{
+  if (Density <= gp::Resolution())  Standard_DomainError::Raise();
+  Standard_Real Xp, Yp, Zp;
+  P.Coord (Xp, Yp, Zp);
+  Standard_Real Ixy = - Xp * Yp;
+  Standard_Real Ixz = - Xp * Zp;
+  Standard_Real Iyz = - Yp * Zp;
+  Standard_Real Ixx = Yp * Yp + Zp * Zp;
+  Standard_Real Iyy = Xp * Xp + Zp * Zp;
+  Standard_Real Izz = Xp * Xp + Yp * Yp;
+  Mat Mp (XYZ (Ixx, Ixy, Ixz), XYZ (Ixy, Iyy, Iyz), XYZ (Ixz, Iyz, Izz));
+  if (dim == 0) {
+    dim = Density;
+    g.SetXYZ (P.XYZ().Multiplied (Density));
+    inertia = Mp * Density;
+  }
+  else {
+    Standard_Real X, Y, Z;
+    g.Coord (X, Y, Z);
+    X = X * dim + Xp * Density;
+    Y = Y * dim + Yp * Density;
+    Z = Z * dim + Zp * Density;
+    dim = dim + Density;
+    X /= dim;
+    Y /= dim;
+    Z /= dim;
+    g.SetCoord (X, Y, Z);
+    inertia = inertia + Mp * Density;
+  }      
+}
+
+GProp_PGProps::GProp_PGProps (const Array1OfPnt& Pnts)
+{
+  for (Standard_Integer i = Pnts.Lower(); i <= Pnts.Upper(); i++) AddPoint(Pnts(i));
+}
+
+GProp_PGProps::GProp_PGProps (const Array2OfPnt& Pnts)
+{
+  for (Standard_Integer j = Pnts.LowerCol(); j <= Pnts.UpperCol(); j++)
+    {
+      for (Standard_Integer i = Pnts.LowerRow(); i <= Pnts.UpperRow(); i++) AddPoint(Pnts (i, j));
+    }      
+}
+
+GProp_PGProps::GProp_PGProps (const Array1OfPnt& Pnts,const Array1OfReal& Density)
+{
+  if (Pnts.Length() != Density.Length())  Standard_DomainError::Raise();
+  Standard_Integer ip = Pnts.Lower();
+  Standard_Integer id = Density.Lower();
+  while (id <= Pnts.Upper()) {
+    Standard_Real D = Density (id);
+    if (D <= gp::Resolution()) Standard_DomainError::Raise();
+    AddPoint(Pnts (ip),D); 
+    ip++;  id++;
+  }      
+}
+
+GProp_PGProps::GProp_PGProps (const Array2OfPnt& Pnts,const Array2OfReal& Density)
+{
+  if (Pnts.ColLength() != Density.ColLength() || Pnts.RowLength() != Density.RowLength()) Standard_DomainError::Raise();
+  Standard_Integer ip = Pnts.LowerRow();
+  Standard_Integer id = Density.LowerRow();
+  Standard_Integer jp = Pnts.LowerCol();
+  Standard_Integer jd = Density.LowerCol();
+  while (jp <= Pnts.UpperCol()) {
+    while (ip <= Pnts.UpperRow()) {
+      Standard_Real D = Density (id, jd);
+      if (D <= gp::Resolution())  Standard_DomainError::Raise();
+      AddPoint(Pnts (ip, jp),D); 
+      ip++; id++;
+    }      
+    jp++; jd++;
+  }
+}
+
+
+
+void GProp_PGProps::Barycentre(const Array1OfPnt&  Pnts,
+			       const Array1OfReal& Density,
+			       Standard_Real&      Mass, 
+			       Pnt&                G)
+{
+  if (Pnts.Length() != Density.Length())  Standard_DimensionError::Raise();
+  Standard_Integer ip = Pnts.Lower();
+  Standard_Integer id = Density.Lower();
+  Mass = Density (id);
+  XYZ Gxyz = Pnts (ip).XYZ();
+  Gxyz.Multiply (Mass);
+  while (ip <= Pnts.Upper()) {
+    Mass = Mass + Density (id);
+    Gxyz.Add ( (Pnts(ip).XYZ()).Multiplied (Density (id)) );
+    ip++;
+    id++;
+  }
+  Gxyz.Divide (Mass);
+  G.SetXYZ (Gxyz);
+}
+
+
+
+
+void GProp_PGProps::Barycentre(const Array2OfPnt&  Pnts, 
+			       const Array2OfReal& Density,
+			       Standard_Real&      Mass,
+			       Pnt&                G)
+{
+  if (Pnts.RowLength() != Density.RowLength() || Pnts.ColLength() != Density.ColLength()) Standard_DimensionError::Raise();
+  Standard_Integer ip = Pnts.LowerRow();
+  Standard_Integer id = Density.LowerRow();
+  Standard_Integer jp = Pnts.LowerCol();
+  Standard_Integer jd = Density.LowerCol();
+  Mass = 0.0;
+  XYZ Gxyz (0.0, 0.0, 0.0);
+  while (jp <= Pnts.UpperCol()) {
+    while (ip <= Pnts.UpperRow()) {
+      Mass = Mass + Density (id, jd);
+      Gxyz.Add ((Pnts(ip, jp).XYZ()).Multiplied (Density (id, jd)));
+      ip++;   id++;
+    }
+    jp++;   jd++;
+  }
+  Gxyz.Divide (Mass);
+  G.SetXYZ (Gxyz);
+}
+
+
+
+
+Pnt GProp_PGProps::Barycentre (const Array1OfPnt& Pnts) {
+
+  XYZ Gxyz = Pnts (Pnts.Lower()).XYZ();
+  for (Standard_Integer i = Pnts.Lower() + 1; i <= Pnts.Upper(); i++) {
+    Gxyz.Add (Pnts(i).XYZ());
+  }
+  Gxyz.Divide (Pnts.Length());
+  return Pnt (Gxyz);
+}
+
+
+
+
+Pnt GProp_PGProps::Barycentre (const Array2OfPnt& Pnts) {
+
+  XYZ Gxyz (0.0, 0.0, 0.0);
+  for (Standard_Integer j = Pnts.LowerCol(); j <= Pnts.UpperCol(); j++) {
+    for (Standard_Integer i = Pnts.LowerRow(); i <= Pnts.UpperRow(); i++) {
+      Gxyz.Add (Pnts(i, j).XYZ());
+    }
+  }
+  Gxyz.Divide (Pnts.RowLength() * Pnts.ColLength());
+  return Pnt (Gxyz);
+}
+
+
+

src/GProp/GProp_PrincipalProps.cdl

@@ -0,0 +1,173 @@
+-- File:	PrincipalProps.cdl
+-- Created:	Mon Feb 17 17:50:44 1992
+-- Author:	Jean Claude VAUTHIER
+---Copyright:	 Matra Datavision 1992
+
+
+class PrincipalProps
+
+
+from GProp
+
+---Purpose:
+-- A framework to present the principal properties of
+-- inertia of a system of which global properties are
+-- computed by a GProp_GProps object.
+-- There is always a set of axes for which the
+-- products of inertia of a geometric system are equal
+-- to 0; i.e. the matrix of inertia of the system is
+-- diagonal. These axes are the principal axes of
+-- inertia. Their origin is coincident with the center of
+-- mass of the system. The associated moments are
+-- called the principal moments of inertia.
+-- This sort of presentation object is created, filled and
+-- returned by the function PrincipalProperties for
+-- any GProp_GProps object, and can be queried to access the result.
+-- Note: The system whose principal properties of
+-- inertia are returned by this framework is referred to
+-- as the current system. The current system,
+-- however, is retained neither by this presentation
+-- framework nor by the GProp_GProps object which activates it.
+uses Vec from gp,
+     Pnt from gp
+     
+ raises UndefinedAxis from GProp
+ 
+ is
+
+
+
+
+  Create   returns PrincipalProps;
+        --- Purpose : creates an undefined PrincipalProps.   
+
+
+  HasSymmetryAxis (me)  returns Boolean is static;
+    	--- Purpose :
+    	--  returns true if the geometric system has an axis of symmetry.  
+	--  For  comparing  moments  relative  tolerance  1.e-10  is  used. 
+	--  Usually  it  is  enough  for  objects,  restricted  by  faces  with 
+	--  analitycal  geometry.
+	
+  HasSymmetryAxis (me; aTol :  Real)  returns Boolean is static;
+    	--- Purpose :
+    	--  returns true if the geometric system has an axis of symmetry. 
+	--  aTol  is  relative  tolerance for  cheking  equality  of  moments
+	--  If  aTol  ==  0,  relative  tolerance  is  ~  1.e-16  (Epsilon(I))
+
+
+  HasSymmetryPoint (me)  returns Boolean is static;
+    	--- Purpose :
+    	--  returns true if the geometric system has a point of symmetry. 
+	--  For  comparing  moments  relative  tolerance  1.e-10  is  used. 
+	--  Usually  it  is  enough  for  objects,  restricted  by  faces  with 
+	--  analitycal  geometry.
+	
+  HasSymmetryPoint (me; aTol :  Real)  returns Boolean is static;
+    	--- Purpose :
+    	--  returns true if the geometric system has a point of symmetry.
+	--  aTol  is  relative  tolerance for  cheking  equality  of  moments 
+	--  If  aTol  ==  0,  relative  tolerance  is  ~  1.e-16  (Epsilon(I))
+
+
+  Moments (me; Ixx, Iyy, Izz: out Real) is static;
+    	--- Purpose : Ixx, Iyy and Izz return the principal moments of inertia
+-- in the current system.
+-- Notes :
+-- - If the current system has an axis of symmetry, two
+--   of the three values Ixx, Iyy and Izz are equal. They
+--   indicate which eigen vectors define an infinity of
+--   axes of principal inertia.
+-- - If the current system has a center of symmetry, Ixx,
+--   Iyy and Izz are equal.
+
+
+  FirstAxisOfInertia (me)   returns Vec
+    	--- Purpose :  returns the first axis of inertia.
+     raises UndefinedAxis
+    	--- Purpose :
+	--  if the system has a point of symmetry there is an infinity of 
+	--  solutions. It is not possible to defines the three axis of 
+	--  inertia. 
+        ---C++: return const&
+     is static;
+
+  SecondAxisOfInertia (me)	returns Vec
+    	--- Purpose :  returns the second axis of inertia.
+     raises  UndefinedAxis
+	--- Purpose :
+	--  if the system has a point of symmetry or an axis of symmetry the 
+	--  second and the third axis of symmetry are undefined.
+        ---C++: return const&
+     is static;
+
+
+  ThirdAxisOfInertia (me)   returns Vec
+    --- Purpose :  returns the third axis of inertia.
+    --     This and the above functions return the first, second or third eigen vector of the
+    -- matrix of inertia of the current system.
+    -- The first, second and third principal axis of inertia
+    -- pass through the center of mass of the current
+    -- system. They are respectively parallel to these three eigen vectors.
+    -- Note that:
+    -- - If the current system has an axis of symmetry, any
+    --   axis is an axis of principal inertia if it passes
+    --   through the center of mass of the system, and runs
+    --   parallel to a linear combination of the two eigen
+    --   vectors of the matrix of inertia, corresponding to the
+    --  two eigen values which are equal. If the current
+    --  system has a center of symmetry, any axis passing
+    --  through the center of mass of the system is an axis
+    --  of principal inertia. Use the functions
+    --  HasSymmetryAxis and HasSymmetryPoint to
+    --  check these particular cases, where the returned
+    --  eigen vectors define an infinity of principal axis of inertia.
+    -- - The Moments function can be used to know which
+    --   of the three eigen vectors corresponds to the two
+    --   eigen values which are equal.
+     raises  UndefinedAxis
+	--- Purpose :
+	--  if the system has a point of symmetry or an axis of symmetry the
+	--  second and the third axis of symmetry are undefined.
+        ---C++: return const&
+   is static;
+
+
+  RadiusOfGyration (me; Rxx, Ryy, Rzz : out Real) is static;
+    	--- Purpose :  Returns the principal radii of gyration  Rxx, Ryy
+-- and Rzz are the radii of gyration of the current
+-- system about its three principal axes of inertia.
+-- Note that:
+-- - If the current system has an axis of symmetry,
+--   two of the three values Rxx, Ryy and Rzz are equal.
+-- - If the current system has a center of symmetry,
+--   Rxx, Ryy and Rzz are equal.
+
+
+
+
+
+  Create (Ixx, Iyy, Izz, Rxx, Ryy, Rzz : Real; Vxx, Vyy, Vzz : Vec; G : Pnt)
+     returns PrincipalProps
+     is private;
+
+
+
+fields
+
+   i1 : Real;
+   i2 : Real;
+   i3 : Real;
+   r1 : Real;
+   r2 : Real;
+   r3 : Real;
+   v1 : Vec;
+   v2 : Vec;
+   v3 : Vec;
+   g  : Pnt;
+   
+friends
+
+  PrincipalProperties from GProps (me)
+  
+end PrincipalProps;

src/GProp/GProp_PrincipalProps.cxx

@@ -0,0 +1,96 @@
+#include <GProp_PrincipalProps.ixx>
+
+
+typedef gp_Vec Vec;
+typedef gp_Pnt Pnt;
+
+
+
+
+
+   GProp_PrincipalProps::GProp_PrincipalProps () {
+
+      i1= i2 = i3 = RealLast();
+      r1 = r2 = r3 = RealLast();
+      v1 = Vec (1.0, 0.0, 0.0);
+      v2 = Vec (0.0, 1.0, 0.0);
+      v3 = Vec (0.0, 0.0, 1.0);
+      g = Pnt (RealLast(), RealLast(), RealLast());
+   }
+
+
+   GProp_PrincipalProps::GProp_PrincipalProps (
+   const Standard_Real Ixx, const Standard_Real Iyy, const Standard_Real Izz, 
+   const Standard_Real Rxx, const Standard_Real Ryy, const Standard_Real Rzz,
+   const gp_Vec& Vxx, const gp_Vec& Vyy, const gp_Vec& Vzz, const gp_Pnt& G) :
+   i1 (Ixx), i2 (Iyy), i3 (Izz), r1 (Rxx), r2 (Ryy), r3 (Rzz),
+   v1 (Vxx), v2 (Vyy), v3 (Vzz), g (G) { }
+
+
+   Standard_Boolean GProp_PrincipalProps::HasSymmetryAxis () const {
+
+//     Standard_Real Eps1 = Abs(Epsilon (i1));
+//     Standard_Real Eps2 = Abs(Epsilon (i2));
+     const Standard_Real aRelTol = 1.e-10;
+     Standard_Real Eps1 = Abs(i1)*aRelTol;
+     Standard_Real Eps2 = Abs(i2)*aRelTol;
+     return (Abs (i1 - i2) <= Eps1 || Abs (i1 - i3) <= Eps1 ||
+             Abs (i2 - i3) <= Eps2);
+   }
+
+   Standard_Boolean GProp_PrincipalProps::HasSymmetryAxis (const Standard_Real aTol) const {
+
+
+     Standard_Real Eps1 = Abs(i1*aTol) + Abs(Epsilon(i1));
+     Standard_Real Eps2 = Abs(i2*aTol) + Abs(Epsilon(i2));
+     return (Abs (i1 - i2) <= Eps1 || Abs (i1 - i3) <= Eps1 ||
+             Abs (i2 - i3) <= Eps2);
+   }
+
+
+   Standard_Boolean GProp_PrincipalProps::HasSymmetryPoint () const {
+
+//     Standard_Real Eps1 = Abs(Epsilon (i1));
+     const Standard_Real aRelTol = 1.e-10;
+     Standard_Real Eps1 = Abs(i1)*aRelTol;
+     return (Abs (i1 - i2) <= Eps1 && Abs (i1 - i3) <= Eps1);
+   }
+
+   Standard_Boolean GProp_PrincipalProps::HasSymmetryPoint (const Standard_Real aTol) const {
+
+     Standard_Real Eps1 = Abs(i1*aTol) + Abs(Epsilon(i1));
+     return (Abs (i1 - i2) <= Eps1 && Abs (i1 - i3) <= Eps1);
+   }
+
+
+   void GProp_PrincipalProps::Moments (Standard_Real& Ixx, Standard_Real& Iyy, Standard_Real& Izz) const {
+
+      Ixx = i1;
+      Iyy = i2;
+      Izz = i3;
+   }
+
+
+   const Vec& GProp_PrincipalProps::FirstAxisOfInertia () const {return  v1;}
+
+
+   const Vec& GProp_PrincipalProps::SecondAxisOfInertia () const {return  v2;}
+
+
+   const Vec& GProp_PrincipalProps::ThirdAxisOfInertia () const {return v3;}
+
+
+   void GProp_PrincipalProps::RadiusOfGyration (
+   Standard_Real& Rxx, Standard_Real& Ryy, Standard_Real& Rzz) const {
+   
+     Rxx = r1;
+     Ryy = r2;
+     Rzz = r3;
+   }
+
+
+
+
+
+
+

src/GProp/GProp_SGProps.cdl

@@ -0,0 +1,51 @@
+-- File:         SGProps.cdl
+-- Created:	 Fri Apr 12 09:43:09 1991
+-- Author:	 Michel CHAUVAT
+--		 Jean-Claude VAUTHIER January 1992
+---Copyright:	 Matra Datavision 1992
+
+
+generic class SGProps from GProp ( Arc      as any;
+    	    	    	    	   Face     as any; --as FaceTool (Arc)
+    	    	    	    	   Domain   as any  --as DomainTool(Arc)
+    	    	    	    	  )
+inherits GProps
+
+    	--- Purpose : 
+        --  Computes the global properties of a face in 3D space. 
+        --  The face 's requirements to evaluate the global properties
+        --  are defined in the template FaceTool from package GProp.
+	
+uses  Pnt  from gp
+is
+
+  Create  returns SGProps;
+
+  Create (S: Face; SLocation: Pnt) returns SGProps;
+  Create (S : in out Face;  D : in out Domain; SLocation : Pnt) returns SGProps; 
+        --- Purpose :
+        --  Builds a SGProps to evaluate the global properties of
+        --  the face <S>. If isNaturalRestriction is true the domain of S is defined  
+    	--  with the natural bounds, else it defined with an iterator
+        --  of Arc (see DomainTool from GProp) 
+  Create (S: in out Face; SLocation: Pnt; Eps: Real) returns SGProps;
+  Create (S: in out Face; D : in out Domain; SLocation: Pnt; Eps: Real) returns SGProps;
+    	--  --"--
+	--  Parameter Eps sets maximal relative error of computed area.
+
+  SetLocation(me: in out; SLocation: Pnt);
+
+  Perform(me: in out; S: Face); 
+  Perform(me : in out; S : in out Face ; D : in out Domain);
+  Perform(me: in out; S: in out Face; Eps: Real) returns Real;
+  Perform(me: in out; S: in out Face; D : in out Domain; Eps: Real) returns Real;
+ 
+  GetEpsilon(me: out) returns Real;
+        --- Purpose :  
+        --  If previously used method contained Eps parameter  
+    	--  get actual relative error of the computation, else return  1.0.
+fields
+
+    myEpsilon: Real from Standard;
+
+end SGProps;

src/GProp/GProp_SGProps.gxx

@@ -0,0 +1,708 @@
+#include <Standard_NotImplemented.hxx>
+#include <math_Vector.hxx>
+#include <math.hxx>
+#include <gp_Pnt2d.hxx>
+#include <gp_Vec2d.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Vec.hxx>
+
+#include <TColStd_Array1OfReal.hxx>
+#include <Precision.hxx>
+
+class HMath_Vector{
+  math_Vector *pvec;
+  void operator=(const math_Vector&){}
+ public:
+  HMath_Vector(){ pvec = 0;}
+  HMath_Vector(math_Vector* pv){ pvec = pv;}
+  ~HMath_Vector(){ if(pvec != 0) delete pvec;}
+  void operator=(math_Vector* pv){ if(pvec != pv && pvec != 0) delete pvec;  pvec = pv;}
+  Standard_Real& operator()(Standard_Integer i){ return (*pvec).operator()(i);}
+  const Standard_Real& operator()(Standard_Integer i) const{ return (*pvec).operator()(i);}
+  const math_Vector* operator->() const{ return pvec;}
+  math_Vector* operator->(){ return pvec;}
+  math_Vector* Vector(){ return pvec;}
+  math_Vector* Init(Standard_Real v, Standard_Integer i = 0, Standard_Integer iEnd = 0){ 
+    if(pvec == 0) return pvec;
+    if(iEnd - i == 0) pvec->Init(v);
+    else { Standard_Integer End = (iEnd <= pvec->Upper()) ? iEnd : pvec->Upper();
+           for(; i <= End; i++) pvec->operator()(i) = v; }
+    return pvec;
+  }
+};
+
+static Standard_Real EPS_PARAM          = 1.e-12;
+static Standard_Real EPS_DIM            = 1.e-20;
+static Standard_Real ERROR_ALGEBR_RATIO = 2.0/3.0;
+
+static Standard_Integer  GPM        = 61;
+static Standard_Integer  SUBS_POWER = 32;
+static Standard_Integer  SM         = 1953;
+
+static math_Vector LGaussP0(1,GPM);
+static math_Vector LGaussW0(1,GPM);
+static math_Vector LGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
+static math_Vector LGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
+
+static math_Vector* LGaussP[] = {&LGaussP0,&LGaussP1};
+static math_Vector* LGaussW[] = {&LGaussW0,&LGaussW1};
+
+static HMath_Vector L1    = new math_Vector(1,SM,0.0);
+static HMath_Vector L2    = new math_Vector(1,SM,0.0);
+static HMath_Vector DimL  = new math_Vector(1,SM,0.0);
+static HMath_Vector ErrL  = new math_Vector(1,SM,0.0);
+static HMath_Vector ErrUL = new math_Vector(1,SM,0.0);
+static HMath_Vector IxL   = new math_Vector(1,SM,0.0);
+static HMath_Vector IyL   = new math_Vector(1,SM,0.0);
+static HMath_Vector IzL   = new math_Vector(1,SM,0.0);
+static HMath_Vector IxxL  = new math_Vector(1,SM,0.0);
+static HMath_Vector IyyL  = new math_Vector(1,SM,0.0);
+static HMath_Vector IzzL  = new math_Vector(1,SM,0.0);
+static HMath_Vector IxyL  = new math_Vector(1,SM,0.0);
+static HMath_Vector IxzL  = new math_Vector(1,SM,0.0);
+static HMath_Vector IyzL  = new math_Vector(1,SM,0.0);
+
+static math_Vector UGaussP0(1,GPM);
+static math_Vector UGaussW0(1,GPM);
+static math_Vector UGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
+static math_Vector UGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
+
+static math_Vector* UGaussP[] = {&UGaussP0,&UGaussP1};
+static math_Vector* UGaussW[] = {&UGaussW0,&UGaussW1};
+
+static HMath_Vector U1   = new math_Vector(1,SM,0.0);
+static HMath_Vector U2   = new math_Vector(1,SM,0.0);
+static HMath_Vector DimU = new math_Vector(1,SM,0.0);
+static HMath_Vector ErrU = new math_Vector(1,SM,0.0);
+static HMath_Vector IxU  = new math_Vector(1,SM,0.0);
+static HMath_Vector IyU  = new math_Vector(1,SM,0.0);
+static HMath_Vector IzU  = new math_Vector(1,SM,0.0);
+static HMath_Vector IxxU = new math_Vector(1,SM,0.0);
+static HMath_Vector IyyU = new math_Vector(1,SM,0.0);
+static HMath_Vector IzzU = new math_Vector(1,SM,0.0);
+static HMath_Vector IxyU = new math_Vector(1,SM,0.0);
+static HMath_Vector IxzU = new math_Vector(1,SM,0.0);
+static HMath_Vector IyzU = new math_Vector(1,SM,0.0);
+
+static Standard_Integer FillIntervalBounds(Standard_Real               A,
+                                           Standard_Real               B,
+                                           const TColStd_Array1OfReal& Knots, 
+					   HMath_Vector&               VA,
+                                           HMath_Vector&               VB)
+{
+  Standard_Integer i = 1, iEnd = Knots.Upper(), j = 1, k = 1;
+  VA(j++) = A;
+  for(; i <= iEnd; i++){
+    Standard_Real kn = Knots(i);
+    if(A < kn)
+      if(kn < B) VA(j++) = VB(k++) = kn; else break;
+  }
+  VB(k) = B;
+  return k;
+}
+
+static inline Standard_Integer MaxSubs(Standard_Integer n, Standard_Integer coeff = SUBS_POWER){
+//  return n = IntegerLast()/coeff < n? IntegerLast(): n*coeff + 1;
+  return Min((n * coeff + 1),SM);
+}
+
+static Standard_Integer LFillIntervalBounds(Standard_Real               A,
+                                            Standard_Real               B,
+                                            const TColStd_Array1OfReal& Knots, 
+					    const Standard_Integer      NumSubs)
+{
+  Standard_Integer iEnd = Knots.Upper(), jEnd = L1->Upper();
+  if(iEnd - 1 > jEnd){
+    iEnd = MaxSubs(iEnd-1,NumSubs); 
+    L1    = new math_Vector(1,iEnd);
+    L2    = new math_Vector(1,iEnd);
+    DimL  = new math_Vector(1,iEnd);
+    ErrL  = new math_Vector(1,iEnd,0.0);
+    ErrUL = new math_Vector(1,iEnd,0.0);
+    IxL   = new math_Vector(1,iEnd);
+    IyL   = new math_Vector(1,iEnd);
+    IzL   = new math_Vector(1,iEnd);
+    IxxL  = new math_Vector(1,iEnd);
+    IyyL  = new math_Vector(1,iEnd);
+    IzzL  = new math_Vector(1,iEnd);
+    IxyL  = new math_Vector(1,iEnd);
+    IxzL  = new math_Vector(1,iEnd);
+    IyzL  = new math_Vector(1,iEnd);
+  }
+  return FillIntervalBounds(A, B, Knots, L1, L2);
+}
+
+static Standard_Integer UFillIntervalBounds(Standard_Real               A,
+                                            Standard_Real               B,
+                                            const TColStd_Array1OfReal& Knots, 
+					    const Standard_Integer      NumSubs)
+{
+  Standard_Integer iEnd = Knots.Upper(), jEnd = U1->Upper();
+  if(iEnd - 1 > jEnd){
+    iEnd = MaxSubs(iEnd-1,NumSubs); 
+    U1   = new math_Vector(1,iEnd);
+    U2   = new math_Vector(1,iEnd);
+    DimU = new math_Vector(1,iEnd);
+    ErrU = new math_Vector(1,iEnd,0.0);
+    IxU  = new math_Vector(1,iEnd);
+    IyU  = new math_Vector(1,iEnd);
+    IzU  = new math_Vector(1,iEnd);
+    IxxU = new math_Vector(1,iEnd);
+    IyyU = new math_Vector(1,iEnd);
+    IzzU = new math_Vector(1,iEnd);
+    IxyU = new math_Vector(1,iEnd);
+    IxzU = new math_Vector(1,iEnd);
+    IyzU = new math_Vector(1,iEnd);
+  }
+  return FillIntervalBounds(A, B, Knots, U1, U2);
+}
+
+static Standard_Real CCompute(Face&                  S,
+                              Domain&                D,
+                              const gp_Pnt&          loc,
+                              Standard_Real&         Dim,
+                              gp_Pnt&                g,
+                              gp_Mat&                inertia,
+			      const Standard_Real    EpsDim,
+			      const Standard_Boolean isErrorCalculation,
+                              const Standard_Boolean isVerifyComputation) 
+{
+  Standard_Boolean isNaturalRestriction = S.NaturalRestriction();
+
+  Standard_Integer NumSubs = SUBS_POWER;
+
+  Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
+  Dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
+  Standard_Real x, y, z;
+  //boundary curve parametrization
+  Standard_Real l1, l2, lm, lr, l;   
+  //Face parametrization in U and V direction
+  Standard_Real BV1, BV2, v;         
+  Standard_Real BU1, BU2, u1, u2, um, ur, u;
+  S.Bounds (BU1, BU2, BV1, BV2);  u1 = BU1;
+  //location point used to compute the inertia
+  Standard_Real xloc, yloc, zloc;
+  loc.Coord (xloc, yloc, zloc); // use member of parent class
+  //Jacobien (x, y, z) -> (u, v) = ||n||
+  Standard_Real ds;                  
+  //On the Face
+  gp_Pnt Ps;                    
+  gp_Vec VNor;
+  //On the boundary curve u-v
+  gp_Pnt2d Puv;                
+  gp_Vec2d Vuv;
+  Standard_Real Dul;  // Dul = Du / Dl
+  Standard_Real CDim[2], CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz;
+  Standard_Real LocDim[2], LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy, LocIxz, LocIyz;
+  
+  Standard_Real ErrorU, ErrorL, ErrorLMax = 0.0, Eps=0.0, EpsL=0.0, EpsU=0.0;
+
+  Standard_Integer iD = 0, NbLSubs, iLS, iLSubEnd, iGL, iGLEnd, NbLGaussP[2], LRange[2], iL, kL, kLEnd, IL, JL;
+  Standard_Integer i, NbUSubs, iUS, iUSubEnd, iGU, iGUEnd, NbUGaussP[2], URange[2], iU, kU, kUEnd, IU, JU;
+  Standard_Integer UMaxSubs, LMaxSubs;
+  iGLEnd = isErrorCalculation? 2: 1; 
+  for(i = 0; i < 2; i++) {
+    LocDim[i] = 0.0;
+    CDim[i] = 0.0;
+  }
+
+  NbUGaussP[0] = S.SIntOrder(EpsDim);  
+  NbUGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbUGaussP[0]));
+  math::GaussPoints(NbUGaussP[0],UGaussP0);  math::GaussWeights(NbUGaussP[0],UGaussW0);
+  math::GaussPoints(NbUGaussP[1],UGaussP1);  math::GaussWeights(NbUGaussP[1],UGaussW1);
+  
+  NbUSubs = S.SUIntSubs();
+  TColStd_Array1OfReal UKnots(1,NbUSubs+1);
+  S.UKnots(UKnots);
+  
+
+  while (isNaturalRestriction || D.More()) {
+    if(isNaturalRestriction){ 
+      NbLGaussP[0] = Min(2*NbUGaussP[0],math::GaussPointsMax());
+    }else{
+      S.Load(D.Value());  ++iD;
+      NbLGaussP[0] = S.LIntOrder(EpsDim);  
+    }
+
+
+    NbLGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbLGaussP[0]));
+    math::GaussPoints(NbLGaussP[0],LGaussP0);  math::GaussWeights(NbLGaussP[0],LGaussW0);
+    math::GaussPoints(NbLGaussP[1],LGaussP1);  math::GaussWeights(NbLGaussP[1],LGaussW1);
+    
+    NbLSubs = isNaturalRestriction? S.SVIntSubs(): S.LIntSubs();
+
+    TColStd_Array1OfReal LKnots(1,NbLSubs+1);
+    if(isNaturalRestriction){
+      S.VKnots(LKnots); 
+      l1 = BV1; l2 = BV2;
+    }else{
+      S.LKnots(LKnots);
+      l1 = S.FirstParameter();  l2 = S.LastParameter();
+    }
+    ErrorL = 0.0;
+    kLEnd = 1; JL = 0;
+    //OCC503(apo): if(Abs(l2-l1) < EPS_PARAM) continue;
+    if(Abs(l2-l1) > EPS_PARAM) {
+      iLSubEnd = LFillIntervalBounds(l1, l2, LKnots, NumSubs);
+      LMaxSubs = MaxSubs(iLSubEnd);
+      if(LMaxSubs > DimL.Vector()->Upper()) LMaxSubs = DimL.Vector()->Upper();
+      DimL.Init(0.0,1,LMaxSubs);  ErrL.Init(0.0,1,LMaxSubs);  ErrUL.Init(0.0,1,LMaxSubs);
+      do{// while: L
+	if(++JL > iLSubEnd){
+	  LRange[0] = IL = ErrL->Max();  LRange[1] = JL;
+	  L1(JL) = (L1(IL) + L2(IL))/2.0;  L2(JL) = L2(IL);  L2(IL) = L1(JL);
+	}else  LRange[0] = IL = JL;
+	if(JL == LMaxSubs || Abs(L2(JL) - L1(JL)) < EPS_PARAM)
+	if(kLEnd == 1){
+	  DimL(JL) = ErrL(JL) = IxL(JL) = IyL(JL) = IzL(JL) = 
+	    IxxL(JL) = IyyL(JL) = IzzL(JL) = IxyL(JL) = IxzL(JL) = IyzL(JL) = 0.0;
+	}else{
+	  JL--;
+	  EpsL = ErrorL;  Eps = EpsL/0.9;
+	  break;
+	}
+	else
+	  for(kL=0; kL < kLEnd; kL++){
+	    iLS = LRange[kL];
+	    lm = 0.5*(L2(iLS) + L1(iLS));         
+	    lr = 0.5*(L2(iLS) - L1(iLS));
+	    CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
+	    for(iGL=0; iGL < iGLEnd; iGL++){//
+	      CDim[iGL] = 0.0;
+	      for(iL=1; iL<=NbLGaussP[iGL]; iL++){
+		l = lm + lr*(*LGaussP[iGL])(iL);
+		if(isNaturalRestriction){ 
+		  v = l; u2 = BU2; Dul = (*LGaussW[iGL])(iL);
+		}else{
+		  S.D12d (l, Puv, Vuv);
+		  Dul = Vuv.Y()*(*LGaussW[iGL])(iL);  // Dul = Du / Dl
+		  if(Abs(Dul) < EPS_PARAM) continue;
+		  v  = Puv.Y();  u2 = Puv.X();
+		  //Check on cause out off bounds of value current parameter
+		  if(v < BV1) v = BV1; else if(v > BV2) v = BV2;
+		  if(u2 < BU1) u2 = BU1; else if(u2 > BU2) u2 = BU2; 
+		}
+		ErrUL(iLS) = 0.0;
+		kUEnd = 1; JU = 0;
+		if(Abs(u2-u1) < EPS_PARAM) continue;
+		iUSubEnd = UFillIntervalBounds(u1, u2, UKnots, NumSubs);
+		UMaxSubs = MaxSubs(iUSubEnd);
+                if(UMaxSubs > DimU.Vector()->Upper()) UMaxSubs = DimU.Vector()->Upper();
+		DimU.Init(0.0,1,UMaxSubs);  ErrU.Init(0.0,1,UMaxSubs);  ErrorU = 0.0;
+		do{//while: U
+		  if(++JU > iUSubEnd){
+		    URange[0] = IU = ErrU->Max();  URange[1] = JU;  
+		    U1(JU) = (U1(IU)+U2(IU))/2.0;  U2(JU) = U2(IU);  U2(IU) = U1(JU);
+		  }else  URange[0] = IU = JU;
+		  if(JU == UMaxSubs || Abs(U2(JU) - U1(JU)) < EPS_PARAM)
+		    if(kUEnd == 1){
+		      DimU(JU) = ErrU(JU) = IxU(JU) = IyU(JU) = IzU(JU) = 
+			IxxU(JU) = IyyU(JU) = IzzU(JU) = IxyU(JU) = IxzU(JU) = IyzU(JU) = 0.0;
+		    }else{
+		      JU--;  
+		      EpsU = ErrorU;  Eps = EpsU*Abs((u2-u1)*Dul)/0.1;  EpsL = 0.9*Eps;
+		      break;
+		    }
+		  else
+		    for(kU=0; kU < kUEnd; kU++){
+		      iUS = URange[kU];
+		      um = 0.5*(U2(iUS) + U1(iUS));
+		      ur = 0.5*(U2(iUS) - U1(iUS));
+		      LocIx = LocIy = LocIz = LocIxx = LocIyy = LocIzz = LocIxy = LocIxz = LocIyz = 0.0;
+		      iGUEnd = iGLEnd - iGL;
+		      for(iGU=0; iGU < iGUEnd; iGU++){//
+			LocDim[iGU] = 0.0;
+			for(iU=1; iU<=NbUGaussP[iGU]; iU++){
+			  u = um + ur*(*UGaussP[iGU])(iU);
+			  S.Normal(u, v, Ps, VNor);
+			  ds = VNor.Magnitude();    //Jacobien(x,y,z) -> (u,v)=||n||
+			  ds *= (*UGaussW[iGU])(iU); 
+			  LocDim[iGU] += ds; 
+			  if(iGU > 0) continue;
+			  Ps.Coord(x, y, z);  
+			  x -= xloc;  y -= yloc;  z -= zloc;
+			  LocIx += x*ds;  LocIy += y*ds;  LocIz += z*ds;
+			  LocIxy += x*y*ds;  LocIyz += y*z*ds;  LocIxz += x*z*ds;
+			  x *= x;  y *= y;  z *= z;
+			  LocIxx += (y+z)*ds;  LocIyy += (x+z)*ds;  LocIzz += (x+y)*ds;
+			}//for: iU
+		      }//for: iGU
+		      DimU(iUS) = LocDim[0]*ur;
+		      if(iGL > 0) continue;
+		      ErrU(iUS) = Abs(LocDim[1]-LocDim[0])*ur;
+		      IxU(iUS) = LocIx*ur; IyU(iUS) = LocIy*ur; IzU(iUS) = LocIz*ur;
+		      IxxU(iUS) = LocIxx*ur; IyyU(iUS) = LocIyy*ur; IzzU(iUS) = LocIzz*ur;
+		      IxyU(iUS) = LocIxy*ur; IxzU(iUS) = LocIxz*ur; IyzU(iUS) = LocIyz*ur;
+		    }//for: kU (iUS)
+		  if(JU == iUSubEnd)  kUEnd = 2;
+		  if(kUEnd == 2)  ErrorU = ErrU(ErrU->Max());
+		}while((ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1);
+		for(i=1; i<=JU; i++)  CDim[iGL] += DimU(i)*Dul;
+		if(iGL > 0) continue;
+		ErrUL(iLS) = ErrorU*Abs((u2-u1)*Dul);
+		for(i=1; i<=JU; i++){
+		  CIx += IxU(i)*Dul; CIy += IyU(i)*Dul; CIz += IzU(i)*Dul;
+		  CIxx += IxxU(i)*Dul; CIyy += IyyU(i)*Dul; CIzz += IzzU(i)*Dul;
+		  CIxy += IxyU(i)*Dul; CIxz += IxzU(i)*Dul; CIyz += IyzU(i)*Dul;
+		}
+	      }//for: iL 
+	    }//for: iGL
+	    DimL(iLS) = CDim[0]*lr;  
+	    if(iGLEnd == 2) ErrL(iLS) = Abs(CDim[1]-CDim[0])*lr + ErrUL(iLS);
+	    IxL(iLS) = CIx*lr; IyL(iLS) = CIy*lr; IzL(iLS) = CIz*lr; 
+	    IxxL(iLS) = CIxx*lr; IyyL(iLS) = CIyy*lr; IzzL(iLS) = CIzz*lr; 
+	    IxyL(iLS) = CIxy*lr; IxzL(iLS) = CIxz*lr; IyzL(iLS) = CIyz*lr; 
+	  }//for: (kL)iLS
+	//  Calculate/correct epsilon of computation by current value of Dim
+	//That is need for not spend time for 
+	if(JL == iLSubEnd){  
+	  kLEnd = 2; 
+	  Standard_Real DDim = 0.0;
+	  for(i=1; i<=JL; i++) DDim += DimL(i);
+	  DDim = Abs(DDim*EpsDim);
+	  if(DDim > Eps) { 
+	    Eps = DDim;  EpsL = 0.9*Eps;
+	  }
+	}
+	if(kLEnd == 2) ErrorL = ErrL(ErrL->Max());
+      }while((ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1);
+      for(i=1; i<=JL; i++){
+	Dim += DimL(i); 
+	Ix += IxL(i); Iy += IyL(i); Iz += IzL(i);
+	Ixx += IxxL(i); Iyy += IyyL(i); Izz += IzzL(i);
+	Ixy += IxyL(i); Ixz += IxzL(i); Iyz += IyzL(i);
+      }
+      ErrorLMax = Max(ErrorLMax, ErrorL);
+    }
+    if(isNaturalRestriction) break;
+    D.Next();
+  }
+  if(Abs(Dim) >= EPS_DIM){
+    Ix /= Dim;  Iy /= Dim;  Iz /= Dim;
+    g.SetCoord (Ix, Iy, Iz);   
+  }else{
+    Dim =0.0;
+    g.SetCoord (0., 0.,0.);
+  }
+  inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
+		    gp_XYZ (-Ixy, Iyy, -Iyz),
+		    gp_XYZ (-Ixz, -Iyz, Izz));
+
+  if(iGLEnd == 2) Eps = Dim != 0.0? ErrorLMax/Abs(Dim): 0.0;
+  else Eps = EpsDim;
+  return Eps;
+}
+
+static Standard_Real Compute(Face& S, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, 
+			     Standard_Real EpsDim) 
+{
+  Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
+  Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
+  EpsDim = Abs(EpsDim);
+  Domain D;
+  return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
+}
+
+static Standard_Real Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, 
+			     Standard_Real EpsDim) 
+{
+  Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
+  Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
+  EpsDim = Abs(EpsDim);
+  return CCompute(S,D,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
+}
+
+static void Compute(Face& S, Domain& D, const gp_Pnt& loc, Standard_Real& dim, gp_Pnt& g, gp_Mat& inertia){
+   Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
+   dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
+
+   Standard_Real x, y, z;
+   Standard_Integer NbCGaussgp_Pnts = 0;
+
+   Standard_Real l1, l2, lm, lr, l;   //boundary curve parametrization
+   Standard_Real v1, v2, vm, vr, v;   //Face parametrization in v direction
+   Standard_Real u1, u2, um, ur, u;
+   Standard_Real ds;                  //Jacobien (x, y, z) -> (u, v) = ||n||
+
+   gp_Pnt P;                    //On the Face
+   gp_Vec VNor;
+
+   gp_Pnt2d Puv;                //On the boundary curve u-v
+   gp_Vec2d Vuv;
+   Standard_Real Dul;                 // Dul = Du / Dl
+   Standard_Real CArea, CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz;
+   Standard_Real LocArea, LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy,
+        LocIxz, LocIyz;
+
+
+   S.Bounds (u1, u2, v1, v2);
+
+   Standard_Integer NbUGaussgp_Pnts = Min(S.UIntegrationOrder (),
+					  math::GaussPointsMax());
+   Standard_Integer NbVGaussgp_Pnts = Min(S.VIntegrationOrder (),
+					  math::GaussPointsMax());
+
+   Standard_Integer NbGaussgp_Pnts = Max(NbUGaussgp_Pnts, NbVGaussgp_Pnts);
+
+   //Number of Gauss points for the integration
+   //on the Face
+   math_Vector GaussSPV (1, NbGaussgp_Pnts);
+   math_Vector GaussSWV (1, NbGaussgp_Pnts);
+   math::GaussPoints  (NbGaussgp_Pnts,GaussSPV);
+   math::GaussWeights (NbGaussgp_Pnts,GaussSWV);
+
+
+   //location point used to compute the inertia
+   Standard_Real xloc, yloc, zloc;
+   loc.Coord (xloc, yloc, zloc);
+
+   while (D.More()) {
+
+      S.Load(D.Value());
+      NbCGaussgp_Pnts =  Min(S.IntegrationOrder (), math::GaussPointsMax());        
+      
+      math_Vector GaussCP (1, NbCGaussgp_Pnts);
+      math_Vector GaussCW (1, NbCGaussgp_Pnts);
+      math::GaussPoints  (NbCGaussgp_Pnts,GaussCP);
+      math::GaussWeights (NbCGaussgp_Pnts,GaussCW);
+
+      CArea = 0.0;
+      CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
+      l1 = S.FirstParameter ();
+      l2 = S.LastParameter  ();
+      lm = 0.5 * (l2 + l1);         
+      lr = 0.5 * (l2 - l1);
+
+      Puv = S.Value2d (lm);
+      vm  = Puv.Y();
+      Puv = S.Value2d (lr);
+      vr  = Puv.Y();
+
+      for (Standard_Integer i = 1; i <= NbCGaussgp_Pnts; i++) {
+        l = lm + lr * GaussCP (i);
+        S.D12d(l, Puv, Vuv);
+        v   = Puv.Y();
+        u2  = Puv.X();
+        Dul = Vuv.Y();
+        Dul *= GaussCW (i);
+        um  = 0.5 * (u2 + u1);
+        ur  = 0.5 * (u2 - u1);
+        LocArea = LocIx  = LocIy  = LocIz = LocIxx = LocIyy = LocIzz = 
+        LocIxy  = LocIxz = LocIyz = 0.0;
+        for (Standard_Integer j = 1; j <= NbGaussgp_Pnts; j++) {
+          u = um + ur * GaussSPV (j);
+          S.Normal (u, v, P, VNor);
+          ds = VNor.Magnitude();    //normal.Magnitude
+          ds = ds * Dul * GaussSWV (j); 
+          LocArea +=  ds; 
+          P.Coord (x, y, z);
+          x      -= xloc;
+          y      -= yloc;
+          z      -= zloc;
+          LocIx  += x * ds;  
+          LocIy  += y * ds;
+          LocIz  += z * ds;
+          LocIxy += x * y * ds;
+          LocIyz += y * z * ds;
+          LocIxz += x * z * ds;
+          x *= x;
+          y *= y;
+          z *= z;
+          LocIxx += (y + z) * ds;
+          LocIyy += (x + z) * ds;
+          LocIzz += (x + y) * ds;
+        }
+        CArea += LocArea * ur;
+        CIx   += LocIx * ur;
+        CIy   += LocIy * ur;
+        CIz   += LocIz * ur;
+        CIxx  += LocIxx * ur;
+        CIyy  += LocIyy * ur;
+        CIzz  += LocIzz * ur;
+        CIxy  += LocIxy * ur;
+        CIxz  += LocIxz * ur;
+        CIyz  += LocIyz * ur;
+      }
+      dim += CArea * lr;
+      Ix  += CIx * lr;
+      Iy  += CIy * lr;
+      Iz  += CIz * lr;
+      Ixx += CIxx * lr;
+      Iyy += CIyy * lr;
+      Izz += CIzz * lr;
+      Ixy += CIxy * lr;
+      Ixz += CIxz * lr;
+      Iyz += CIyz * lr;
+      D.Next();
+   }
+   if (Abs(dim) >= EPS_DIM) {
+     Ix /= dim;
+     Iy /= dim;
+     Iz /= dim;
+     g.SetCoord (Ix, Iy, Iz);
+   }
+   else {
+     dim =0.;
+     g.SetCoord (0., 0.,0.);
+   }
+   inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
+		     gp_XYZ (-Ixy, Iyy, -Iyz),
+		     gp_XYZ (-Ixz, -Iyz, Izz));
+}
+
+
+ 
+static void Compute(const Face& S, const gp_Pnt& loc, Standard_Real& dim, gp_Pnt& g, gp_Mat& inertia){
+   Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
+   dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
+
+   Standard_Real LowerU, UpperU, LowerV, UpperV;
+   S.Bounds (LowerU, UpperU, LowerV, UpperV);
+   Standard_Integer UOrder = Min(S.UIntegrationOrder (),
+				 math::GaussPointsMax());
+   Standard_Integer VOrder = Min(S.VIntegrationOrder (),
+				 math::GaussPointsMax());   
+   gp_Pnt P;          
+   gp_Vec VNor;   
+   Standard_Real dsi, ds;        
+   Standard_Real ur, um, u, vr, vm, v;
+   Standard_Real x, y, z; 
+   Standard_Real Ixi, Iyi, Izi, Ixxi, Iyyi, Izzi, Ixyi, Ixzi, Iyzi;
+   Standard_Real xloc, yloc, zloc;
+   loc.Coord (xloc, yloc, zloc);
+
+   Standard_Integer i, j;
+   math_Vector GaussPU (1, UOrder);     //gauss points and weights
+   math_Vector GaussWU (1, UOrder);
+   math_Vector GaussPV (1, VOrder);
+   math_Vector GaussWV (1, VOrder);
+
+   //Recuperation des points de Gauss dans le fichier GaussPoints.
+   math::GaussPoints  (UOrder,GaussPU);
+   math::GaussWeights (UOrder,GaussWU);
+   math::GaussPoints  (VOrder,GaussPV);
+   math::GaussWeights (VOrder,GaussWV);
+
+   // Calcul des integrales aux points de gauss :
+   um = 0.5 * (UpperU + LowerU);
+   vm = 0.5 * (UpperV + LowerV);
+   ur = 0.5 * (UpperU - LowerU);
+   vr = 0.5 * (UpperV - LowerV);
+
+   for (j = 1; j <= VOrder; j++) {
+     v = vm + vr * GaussPV (j);
+     dsi = Ixi = Iyi = Izi = Ixxi = Iyyi = Izzi = Ixyi = Ixzi = Iyzi = 0.0;
+
+     for (i = 1; i <= UOrder; i++) {
+       u = um + ur * GaussPU (i);
+       S.Normal (u, v, P, VNor); 
+       ds = VNor.Magnitude() * GaussWU (i);
+       P.Coord (x, y, z);
+       x    -=  xloc;
+       y    -=  yloc;
+       z    -=  zloc;
+       dsi  +=  ds; 
+       Ixi  += x * ds;  
+       Iyi  += y * ds;
+       Izi  += z * ds;
+       Ixyi += x * y * ds;
+       Iyzi += y * z * ds;
+       Ixzi += x * z * ds;
+       x    *= x;
+       y    *= y;
+       z    *= z;
+       Ixxi += (y + z) * ds;
+       Iyyi += (x + z) * ds;
+       Izzi += (x + y) * ds;
+     }
+     dim  += dsi  * GaussWV (j);
+     Ix    += Ixi  * GaussWV (j);
+     Iy    += Iyi  * GaussWV (j);
+     Iz    += Izi  * GaussWV (j);
+     Ixx   += Ixxi * GaussWV (j);
+     Iyy   += Iyyi * GaussWV (j);
+     Izz   += Izzi * GaussWV (j);
+     Ixy   += Ixyi * GaussWV (j);
+     Iyz   += Iyzi * GaussWV (j);
+     Ixz   += Ixzi * GaussWV (j);
+   }
+   vr    *= ur;
+   Ixx   *= vr;
+   Iyy   *= vr;
+   Izz   *= vr;
+   Ixy   *= vr;
+   Ixz   *= vr;
+   Iyz   *= vr;
+   if (Abs(dim) >= EPS_DIM) {
+     Ix    /= dim;
+     Iy    /= dim;
+     Iz    /= dim;
+     dim   *= vr;
+     g.SetCoord (Ix, Iy, Iz);
+   }
+   else {
+     dim =0.;
+     g.SetCoord (0.,0.,0.);
+   }
+   inertia = gp_Mat (gp_XYZ (Ixx, -Ixy, -Ixz),
+		     gp_XYZ (-Ixy, Iyy, -Iyz),
+		     gp_XYZ (-Ixz, -Iyz, Izz));
+}
+
+GProp_SGProps::GProp_SGProps(){}
+
+GProp_SGProps::GProp_SGProps (const Face&   S,
+			      const gp_Pnt& SLocation
+			      ) 
+{
+   SetLocation(SLocation);
+   Perform(S);
+}
+
+GProp_SGProps::GProp_SGProps (Face&   S,
+                              Domain& D,
+			      const gp_Pnt& SLocation
+			      ) 
+{
+   SetLocation(SLocation);
+   Perform(S,D);
+}
+
+GProp_SGProps::GProp_SGProps(Face& S, const gp_Pnt& SLocation, const Standard_Real Eps){
+  SetLocation(SLocation);
+  Perform(S, Eps);
+}
+
+GProp_SGProps::GProp_SGProps(Face& S, Domain& D, const gp_Pnt& SLocation, const Standard_Real Eps){
+  SetLocation(SLocation);
+  Perform(S, D, Eps);
+}
+
+void GProp_SGProps::SetLocation(const gp_Pnt& SLocation){
+  loc = SLocation;
+}
+
+void GProp_SGProps::Perform(const Face& S){
+  Compute(S,loc,dim,g,inertia);
+  myEpsilon = 1.0;
+  return;
+}
+
+void GProp_SGProps::Perform(Face& S, Domain& D){
+  Compute(S,D,loc,dim,g,inertia);
+  myEpsilon = 1.0;
+  return;
+}
+
+Standard_Real GProp_SGProps::Perform(Face& S, const Standard_Real Eps){
+  return myEpsilon = Compute(S,loc,dim,g,inertia,Eps);
+}
+
+Standard_Real GProp_SGProps::Perform(Face& S, Domain& D, const Standard_Real Eps){
+  return myEpsilon = Compute(S,D,loc,dim,g,inertia,Eps);
+}
+
+
+Standard_Real GProp_SGProps::GetEpsilon(){
+  return myEpsilon;
+}

src/GProp/GProp_SelGProps.cdl

@@ -0,0 +1,51 @@
+-- File:	GProp_SelGProps.cdl
+-- Created:	Wed Dec  2 15:46:16 1992
+-- Author:	Isabelle GRIGNON
+--		<isg@sdsun2>
+---Copyright:	 Matra Datavision 1992
+
+
+class SelGProps from GProp inherits GProps
+
+    	    ---Purpose: 
+    	    --          Computes the global properties of an elementary 
+    	    --          surface (surface of the gp package)
+
+uses  Cylinder from gp,
+      Cone     from gp,
+      Pnt      from gp,
+      Sphere   from gp,
+      Torus    from gp
+
+is
+
+  Create returns SelGProps;
+
+  Create (S : Cylinder; Alpha1, Alpha2, Z1, Z2 : Real; SLocation : Pnt)   
+     returns SelGProps;
+
+
+  Create (S : Cone; Alpha1, Alpha2, Z1, Z2 : Real; SLocation : Pnt)   
+     returns SelGProps;
+
+
+  Create (S : Sphere from gp; Teta1, Teta2, Alpha1, Alpha2 : Real;
+          SLocation : Pnt)   
+     returns SelGProps;
+
+
+  Create (S : Torus from gp; Teta1, Teta2, Alpha1, Alpha2 : Real;
+          SLocation : Pnt)   
+     returns SelGProps;
+
+  SetLocation(me : in out ;SLocation :Pnt);
+
+  Perform(me : in out;S : Cylinder; Alpha1, Alpha2, Z1, Z2 : Real);
+
+  Perform(me : in out;S : Cone; Alpha1, Alpha2, Z1, Z2 : Real);
+
+  Perform(me : in out;S : Sphere; Teta1, Teta2, Alpha1, Alpha2 : Real);
+
+  Perform(me : in out;S : Torus; Teta1, Teta2, Alpha1, Alpha2 : Real);
+
+end SelGProps;

src/GProp/GProp_SelGProps.cxx

@@ -0,0 +1,371 @@
+#include <GProp_SelGProps.ixx>
+#include <Standard_NotImplemented.hxx>
+#include <math_Matrix.hxx>
+#include <math_Vector.hxx>
+#include <math_Jacobi.hxx>
+#include <GProp.hxx>
+
+GProp_SelGProps::GProp_SelGProps(){};
+
+void GProp_SelGProps::SetLocation(const gp_Pnt& SLocation )
+{
+  loc = SLocation;
+}
+
+
+void GProp_SelGProps::Perform(const gp_Cylinder& S, 
+			      const Standard_Real         Alpha1,
+			      const Standard_Real         Alpha2, 
+			      const Standard_Real         Z1,
+			      const Standard_Real         Z2)
+ {
+  Standard_Real X0,Y0,Z0,Xa1,Ya1,Za1,Xa2,Ya2,Za2,Xa3,Ya3,Za3;
+  S.Location().Coord(X0,Y0,Z0);
+  Standard_Real R = S.Radius();
+  S.Position().XDirection().Coord(Xa1,Ya1,Za1);
+  S.Position().YDirection().Coord(Xa2,Ya2,Za2);
+  S.Position().Direction().Coord(Xa3,Ya3,Za3);
+  dim = R*(Z2-Z1)*(Alpha2-Alpha1);
+  Standard_Real SA2 = Sin(Alpha2);
+  Standard_Real SA1 = Sin(Alpha1);
+  Standard_Real CA2 = Cos(Alpha2);
+  Standard_Real CA1 = Cos(Alpha1);
+  Standard_Real Ix = R*(SA2-SA1)/(Alpha2-Alpha1);
+  Standard_Real Iy = R*(CA1-CA2)/(Alpha2-Alpha1);
+  g.SetCoord( X0 + Ix*Xa1+Iy*Xa2 + Xa3*(Z2+Z1)/2.,
+	      Y0 + Ix*Ya1+Iy*Ya2 + Ya3*(Z2+Z1)/2.,
+	      Z0 + Ix*Za1+Iy*Za2 + Za3*(Z2+Z1)/2.);
+  Standard_Real ICn2 =R*R*( Alpha2-Alpha1 + SA2*CA2 - SA1*CA1 )/2.;  
+  Standard_Real ISn2 =R*R*( Alpha2-Alpha1 - SA2*CA2 + SA1*CA1 )/2.;  
+  Standard_Real IZ2 = (Alpha2-Alpha1)*(Z2*Z2+Z2*Z1+Z1*Z1)/3.;
+  Standard_Real ICnSn= R*R*( SA2*SA2-SA1*SA1)/2.;
+  Standard_Real ICnz = (Z2+Z1)*(SA2-SA1)/2.;
+  Standard_Real ISnz = (Z2+Z1)*(CA1-CA2)/2.;
+
+  math_Matrix Dm(1,3,1,3);
+
+  Dm(1,1) = ISn2 + IZ2;
+  Dm(2,2) = ICn2 + IZ2;
+  Dm(3,3) = Alpha2-Alpha1;
+  Dm(1,2) = Dm(2,1) = -ICnSn;
+  Dm(1,3) = Dm(3,1) = -ICnz;
+  Dm(3,2) = Dm(2,3) = -ISnz;
+
+  math_Matrix Passage (1,3,1,3);
+  Passage(1,1) = Xa1; Passage(1,2) = Xa2 ;Passage(1,3) = Xa3;
+  Passage(2,1) = Ya1; Passage(2,2) = Ya2 ;Passage(2,3) = Ya3;
+  Passage(3,1) = Za1; Passage(3,2) = Za2 ;Passage(3,3) = Za3;
+
+  math_Jacobi J(Dm);
+  R = R*(Z2-Z1);
+  math_Vector V1(1,3), V2(1,3), V3(1,3);
+  J.Vector(1,V1);
+  V1.Multiply(Passage,V1);
+  V1.Multiply(R*J.Value(1));
+  J.Vector(2,V2);
+  V2.Multiply(Passage,V2);
+  V2.Multiply(R*J.Value(2));
+  J.Vector(3,V3);
+  V3.Multiply(Passage,V3);
+  V3.Multiply(R*J.Value(3));
+
+  inertia = gp_Mat (gp_XYZ(V1(1),V2(1),V3(1)),
+		    gp_XYZ(V1(2),V2(2),V3(2)),
+		    gp_XYZ(V1(3),V2(3),V3(3)));
+  gp_Mat Hop;
+  GProp::HOperator(g,loc,dim,Hop);
+  inertia = inertia+Hop;
+}
+
+void GProp_SelGProps::Perform (const gp_Cone& S, 
+			       const Standard_Real     Alpha1,
+			       const Standard_Real     Alpha2,
+			       const Standard_Real     Z1,
+			       const Standard_Real     Z2)
+
+{
+  Standard_Real X0,Y0,Z0,Xa1,Ya1,Za1,Xa2,Ya2,Za2,Xa3,Ya3,Za3;
+  S.Location().Coord(X0,Y0,Z0);
+  S.Position().XDirection().Coord(Xa1,Ya1,Za1);
+  S.Position().YDirection().Coord(Xa2,Ya2,Za2);
+  S.Position().Direction().Coord(Xa3,Ya3,Za3);
+  Standard_Real t =S.SemiAngle();
+  Standard_Real Cnt = Cos(t);
+  Standard_Real Snt = Sin(t); 
+  Standard_Real R = S.RefRadius();
+  Standard_Real Sn2 = Sin(Alpha2);
+  Standard_Real Sn1 = Sin(Alpha1);
+  Standard_Real Cn2 = Cos(Alpha2);
+  Standard_Real Cn1 = Cos(Alpha1);
+
+  Standard_Real Auxi1 = R + (Z2+Z1)*Snt/2.;
+  Standard_Real Auxi2 = (Z2*Z2+Z1*Z2+Z1*Z1)/3.;
+  dim = (Alpha2-Alpha1)*Cnt*(Z2-Z1)*Auxi1; 
+
+  Standard_Real Ix  =
+    (R*R+R*(Z2+Z1)*Snt + Snt*Auxi2)/Auxi1;
+  Standard_Real Iy = Ix*(Cn1-Cn2)/(Alpha2-Alpha1);
+  Ix = Ix*(Sn2-Sn1)/(Alpha2-Alpha1);
+  Standard_Real Iz = Cnt*(R*(Z2+Z1)/2.+Snt*Auxi2)/Auxi1;
+
+  g.SetCoord(X0 + Xa1*Ix+ Xa2*Iy+ Xa3*Iz,
+	     Y0 + Ya1*Ix+ Ya2*Iy+ Ya3*Iz,
+	     Z0 + Za1*Ix+ Za2*Iy+ Za3*Iz);
+ 
+ Standard_Real R1 = R+Z1*Snt;
+ Standard_Real R2 = R+Z2*Snt;
+ Standard_Real ZZ = (Z2-Z1)*Cnt;
+ Standard_Real IR2 = ZZ*Snt*(R1*R1*R1+R1*R1*R2+R1*R2*R2+R2*R2*R2)/4.;
+ Standard_Real ICn2  = IR2*(Alpha2-Alpha1+Cn2*Sn2-Cn1*Sn1)/2.;
+ Standard_Real ISn2 = IR2*(Alpha2-Alpha1+Cn2*Sn2-Cn1*Sn1)/2.;
+ Standard_Real IZ2 = ZZ*Cnt*Cnt*(Z2-Z1)*(Alpha2-Alpha1)*
+   (R*Auxi2 +Snt*(Z2*Z2*Z2+Z2*Z2*Z1+Z2*Z1*Z1+Z1*Z1*Z1))/4.;
+ Standard_Real ICnSn = IR2*(Cn2*Cn2-Cn1*Cn1);
+ Standard_Real ICnz = Cnt*Snt*ZZ*(R*(Z1+Z2)/2.+Auxi2)*(Sn2-Sn1);
+ Standard_Real ISnz = Cnt*Snt*ZZ*(R*(Z1+Z2)/2.+Auxi2)*(Cn1-Cn2);    
+
+  math_Matrix Dm(1,3,1,3);
+  Dm(1,1) = ISn2 + IZ2;
+  Dm(2,2) = ICn2 + IZ2;
+  Dm(3,3) = IR2*(Alpha2-Alpha1);
+  Dm(1,2) = Dm(2,1) = -ICnSn;
+  Dm(1,3) = Dm(3,1) = -ICnz;
+  Dm(3,2) = Dm(2,3) = -ISnz;
+
+  math_Matrix Passage (1,3,1,3);
+  Passage(1,1) = Xa1; Passage(1,2) = Xa2 ;Passage(1,3) = Xa3;
+  Passage(2,1) = Ya1; Passage(2,2) = Ya2 ;Passage(2,3) = Ya3;
+  Passage(3,1) = Za1; Passage(3,2) = Za2 ;Passage(3,3) = Za3;
+
+  math_Jacobi J(Dm);
+  math_Vector V1(1,3),V2(1,3),V3(1,3);
+  J.Vector(1,V1);
+  V1.Multiply(Passage,V1);
+  V1.Multiply(J.Value(1));
+  J.Vector(2,V2);
+  V2.Multiply(Passage,V2);
+  V2.Multiply(J.Value(2));
+  J.Vector(3,V3);
+  V3.Multiply(Passage,V3);
+  V3.Multiply(J.Value(3));
+
+  inertia = gp_Mat (gp_XYZ(V1(1),V2(1),V3(1)),
+		    gp_XYZ(V1(2),V2(2),V3(2)),
+		    gp_XYZ(V1(3),V2(3),V3(3)));
+  gp_Mat Hop;
+  GProp::HOperator(g,loc,dim,Hop);
+  inertia = inertia+Hop;
+}
+
+
+void GProp_SelGProps::Perform(const gp_Sphere& S,
+			      const Standard_Real       Teta1,
+			      const Standard_Real       Teta2,
+			      const Standard_Real       Alpha1,
+			      const Standard_Real       Alpha2)
+{
+  Standard_Real X0,Y0,Z0,Xa1,Ya1,Za1,Xa2,Ya2,Za2,Xa3,Ya3,Za3;
+  S.Location().Coord(X0,Y0,Z0);
+  S.Position().XDirection().Coord(Xa1,Ya1,Za1);
+  S.Position().YDirection().Coord(Xa2,Ya2,Za2);
+  S.Position().Direction().Coord(Xa3,Ya3,Za3);
+  Standard_Real R = S.Radius();
+  Standard_Real Cnt1 = Cos(Teta1);
+  Standard_Real Snt1 = Sin(Teta1);
+  Standard_Real Cnt2 = Cos(Teta2);
+  Standard_Real Snt2 = Sin(Teta2);
+  Standard_Real Cnf1 = Cos(Alpha1);
+  Standard_Real Snf1 = Sin(Alpha1);
+  Standard_Real Cnf2 = Cos(Alpha2);
+  Standard_Real Snf2 = Sin(Alpha2);
+  dim = R*R*(Teta2-Teta1)*(Snf2-Snf1);
+  Standard_Real Ix = 
+    R*(Snt2-Snt1)/(Teta2-Teta1)*
+    (Alpha2-Alpha1+Snf2*Cnf2-Snf1*Cnf1)/(Snf2-Snf1)/2.;
+  Standard_Real Iy = 
+    R*(Cnt1-Cnt2)/(Teta2-Teta1)*
+    (Alpha2-Alpha1+Snf2*Cnf2-Snf1*Cnf1)/(Snf2-Snf1)/2.;
+  Standard_Real Iz = R*(Snf2+Snf1)/2.;
+  g.SetCoord(
+   X0 + Ix*Xa1 + Iy*Xa2 + Iz*Xa3,
+   Y0 + Ix*Ya1 + Iy*Ya2 + Iz*Ya3,
+   Z0 + Ix*Za1 + Iy*Za2 + Iz*Za3);
+
+  Standard_Real IR2 = ( Cnf2*Snf2*(Cnf2+1.)- Cnf1*Snf1*(Cnf1+1.) +
+                       Alpha2-Alpha1 )/3.;
+  Standard_Real ICn2 = (Teta2-Teta1+ Cnt2*Snt2-Cnt1*Snt1)*IR2/2.;
+  Standard_Real ISn2 = (Teta2-Teta1-Cnt2*Snt2+Cnt1*Snt1)*IR2/2.;
+  Standard_Real ICnSn = ( Snt2*Snt2-Snt1*Snt1)*IR2/2.;
+  Standard_Real IZ2 = (Teta2-Teta1)*(Snf2*Snf2*Snf2-Snf1*Snf1*Snf1)/3.;
+  Standard_Real ICnz =(Snt2-Snt1)*(Cnf1*Cnf1*Cnf1-Cnf2*Cnf2*Cnf2)/3.;
+  Standard_Real ISnz =(Cnt1-Cnt2)*(Cnf1*Cnf1*Cnf1-Cnf2*Cnf2*Cnf2)/3.;
+
+  math_Matrix Dm(1,3,1,3);
+  Dm(1,1) = ISn2 +IZ2; 
+  Dm(2,2) = ICn2 +IZ2;
+  Dm(3,3) = IR2*(Teta2-Teta1);
+  Dm(1,2) = Dm(2,1) = -ICnSn;
+  Dm(1,3) = Dm(3,1) = -ICnz;
+  Dm(3,2) = Dm(2,3) = -ISnz;
+
+  math_Matrix Passage (1,3,1,3);
+  Passage(1,1) = Xa1; Passage(1,2) = Xa2 ;Passage(1,3) = Xa3;
+  Passage(2,1) = Ya1; Passage(2,2) = Ya2 ;Passage(2,3) = Ya3;
+  Passage(3,1) = Za1; Passage(3,2) = Za2 ;Passage(3,3) = Za3;
+
+  math_Jacobi J(Dm);
+  R = R*R*R*R;
+  math_Vector V1(1,3), V2(1,3), V3(1,3);
+  J.Vector(1,V1);
+  V1.Multiply(Passage,V1);
+  V1.Multiply(R*J.Value(1));
+  J.Vector(2,V2);
+  V2.Multiply(Passage,V2);
+  V2.Multiply(R*J.Value(2));
+  J.Vector(3,V3);
+  V3.Multiply(Passage,V3);
+  V3.Multiply(R*J.Value(3));
+
+  inertia = gp_Mat (gp_XYZ(V1(1),V2(1),V3(1)),
+		    gp_XYZ(V1(2),V2(2),V3(2)),
+		    gp_XYZ(V1(3),V2(3),V3(3)));
+  gp_Mat Hop;
+  GProp::HOperator(g,loc,dim,Hop);
+  inertia = inertia+Hop;
+
+ 
+}
+
+
+void GProp_SelGProps::Perform (const gp_Torus& S, 
+			       const Standard_Real      Teta1,
+			       const Standard_Real      Teta2,
+			       const Standard_Real      Alpha1,
+			       const Standard_Real      Alpha2)
+{
+  Standard_Real X0,Y0,Z0,Xa1,Ya1,Za1,Xa2,Ya2,Za2,Xa3,Ya3,Za3;
+  S.Location().Coord(X0,Y0,Z0);
+  S.XAxis().Direction().Coord(Xa1,Ya1,Za1);
+  S.YAxis().Direction().Coord(Xa2,Ya2,Za2);
+  S.Axis().Direction().Coord(Xa3,Ya3,Za3);
+  Standard_Real RMax = S.MajorRadius();
+  Standard_Real Rmin = S.MinorRadius();
+  Standard_Real Cnt1 = Cos(Teta1);
+  Standard_Real Snt1 = Sin(Teta1);
+  Standard_Real Cnt2 = Cos(Teta2);
+  Standard_Real Snt2 = Sin(Teta2);
+  Standard_Real Cnf1 = Cos(Alpha1);
+  Standard_Real Snf1 = Sin(Alpha1);
+  Standard_Real Cnf2 = Cos(Alpha2);
+  Standard_Real Snf2 = Sin(Alpha2);
+
+
+  dim = RMax*Rmin*(Teta2-Teta1)*(Alpha2-Alpha1);
+  Standard_Real Ix = 
+    (Snt2-Snt1)/(Teta2-Teta1)*(Rmin*(Snf2-Snf1)/(Alpha2-Alpha1) + RMax);
+  Standard_Real Iy = 
+    (Cnt1-Cnt2)/(Teta2-Teta1)*(Rmin*(Snf2-Snf1)/(Alpha2-Alpha1) + RMax);
+   Standard_Real Iz = Rmin*(Cnf1-Cnf2)/(Alpha2-Alpha1);
+  
+   g.SetCoord(
+	     X0+Ix*Xa1+Iy*Xa2+Iz*Xa3,
+	     Y0+Ix*Ya1+Iy*Ya2+Iz*Ya3,
+	     Z0+Ix*Za1+Iy*Za2+Iz*Za3);
+
+  Standard_Real IR2 = RMax*RMax + 2.*RMax*Rmin*(Snf2-Snf1) +
+		      Rmin*Rmin/2.*(Snf2*Cnf2-Snf1*Cnf1);
+  Standard_Real ICn2 = IR2*(Teta2-Teta1 +Snt2*Cnt2-Snt1*Cnt1)/2.;
+  Standard_Real ISn2 = IR2*(Teta2-Teta1 -Snt2*Cnt2+Snt1*Cnt1)/2.;
+  Standard_Real ICnSn = IR2*(Snt2*Snt2-Snt1*Snt1)/2.;
+  Standard_Real IZ2 = 
+     (Teta2-Teta1)*Rmin*Rmin*(Alpha2-Alpha1-Snf2*Cnf2+Snf1*Cnf1)/2.;
+  Standard_Real ICnz = Rmin*(Snt2-Snt1)*(Cnf1-Cnf2)*(RMax+Rmin*(Cnf1+Cnf2)/2.);
+  Standard_Real ISnz = Rmin*(Cnt2-Cnt1)*(Cnf1-Cnf2)*(RMax+Rmin*(Cnf1+Cnf2)/2.);
+
+  math_Matrix Dm(1,3,1,3);
+  Dm(1,1) = ISn2 + IZ2; 
+  Dm(2,2) = ICn2 + IZ2;
+  Dm(3,3) = IR2*(Teta2-Teta1);
+  Dm(1,2) = Dm(2,1) = -ICnSn;
+  Dm(1,3) = Dm(3,1) = -ICnz;
+  Dm(3,2) = Dm(2,3) = -ISnz;
+
+  math_Matrix Passage (1,3,1,3);
+  Passage(1,1) = Xa1; Passage(1,2) = Xa2 ;Passage(1,3) = Xa3;
+  Passage(2,1) = Ya1; Passage(2,2) = Ya2 ;Passage(2,3) = Ya3;
+  Passage(3,1) = Za1; Passage(3,2) = Za2 ;Passage(3,3) = Za3;
+
+  math_Jacobi J(Dm);
+  RMax = RMax*Rmin;
+  math_Vector V1(1,3), V2(1,3), V3(1,3);
+  J.Vector(1,V1);
+  V1.Multiply(Passage,V1);
+  V1.Multiply(RMax*J.Value(1));
+  J.Vector(2,V2);
+  V2.Multiply(Passage,V2);
+  V2.Multiply(RMax*J.Value(2));
+  J.Vector(3,V3);
+  V3.Multiply(Passage,V3);
+  V3.Multiply(RMax*J.Value(3));
+
+  inertia = gp_Mat (gp_XYZ(V1(1),V2(1),V3(1)),
+		    gp_XYZ(V1(2),V2(2),V3(2)),
+		    gp_XYZ(V1(3),V2(3),V3(3)));
+  gp_Mat Hop;
+  GProp::HOperator(g,loc,dim,Hop);
+  inertia = inertia+Hop;
+
+
+}
+
+   
+
+GProp_SelGProps::GProp_SelGProps (const gp_Cone& S, 
+				  const Standard_Real     Alpha1,
+				  const Standard_Real     Alpha2,
+				  const Standard_Real     Z1,
+				  const Standard_Real     Z2,
+				  const gp_Pnt& SLocation )
+{
+  SetLocation(SLocation);
+  Perform(S,Alpha1,Alpha2,Z1,Z2);
+}
+
+
+GProp_SelGProps::GProp_SelGProps (const gp_Cylinder& S, 
+				  const Standard_Real         Alpha1,
+				  const Standard_Real         Alpha2, 
+				  const Standard_Real         Z1,
+				  const Standard_Real         Z2,
+				  const gp_Pnt& SLocation)
+ {
+  SetLocation(SLocation);
+  Perform(S,Alpha1,Alpha2,Z1,Z2);
+}
+
+
+GProp_SelGProps::GProp_SelGProps (const gp_Sphere& S,
+				  const Standard_Real       Teta1,
+				  const Standard_Real       Teta2,
+				  const Standard_Real       Alpha1,
+				  const Standard_Real       Alpha2,
+				  const gp_Pnt& SLocation)
+ {
+  SetLocation(SLocation);
+  Perform(S,Teta1,Teta2,Alpha1,Alpha2);
+}
+
+
+GProp_SelGProps::GProp_SelGProps (const gp_Torus& S, 
+				  const Standard_Real      Teta1,
+				  const Standard_Real      Teta2,
+				  const Standard_Real      Alpha1,
+				  const Standard_Real      Alpha2,
+				  const gp_Pnt& SLocation)
+ {
+  SetLocation(SLocation);
+  Perform(S,Teta1,Teta2,Alpha1,Alpha2);
+}
+
+   
+

src/GProp/GProp_TFunction.gxx

@@ -0,0 +1,173 @@
+// File:	GProp_TFunction.gxx
+// Created:	Fri Dec  9 16:23:25 2005
+// Author:	Sergey KHROMOV
+//		<skv@dimox>
+
+
+#include <TColStd_HArray1OfReal.hxx>
+#include <math_KronrodSingleIntegration.hxx>
+#include <Precision.hxx>
+#include <math.hxx>
+
+//=======================================================================
+//function : Constructor.
+//purpose  : 
+//=======================================================================
+
+GProp_TFunction::GProp_TFunction(const Face             &theSurface,
+				 const gp_Pnt           &theVertex,
+				 const Standard_Boolean  IsByPoint,
+				 const Standard_Address  theCoeffs,
+				 const Standard_Real     theUMin,
+				 const Standard_Real     theTolerance)
+     : mySurface(theSurface),
+       myUFunction(mySurface, theVertex, IsByPoint, theCoeffs),
+       myUMin(theUMin),
+       myTolerance(theTolerance),
+       myTolReached(0.),
+       myErrReached(0.),
+       myAbsError(0.),
+       myValueType(GProp_Unknown),
+       myIsByPoint(IsByPoint),
+       myNbPntOuter(3)
+{
+}
+
+//=======================================================================
+//function : Init
+//purpose  : 
+//=======================================================================
+
+void GProp_TFunction::Init() 
+{
+  myTolReached = 0.;
+  myErrReached = 0.;
+  myAbsError = 0.;
+}
+
+//=======================================================================
+//function : Value
+//purpose  : 
+//=======================================================================
+
+Standard_Boolean GProp_TFunction::Value(const Standard_Real  X,
+					      Standard_Real &F)
+{
+
+  const Standard_Real tolU = 1.e-9;
+
+  gp_Pnt2d                      aP2d;
+  gp_Vec2d                      aV2d;
+  Standard_Real                 aUMax;
+  Handle(TColStd_HArray1OfReal) anUKnots;
+
+  mySurface.D12d(X, aP2d, aV2d);
+  aUMax = aP2d.X();
+
+  if(aUMax - myUMin < tolU) {
+    F = 0.;
+    return Standard_True;
+  } 
+
+  mySurface.GetUKnots(myUMin, aUMax, anUKnots);
+  myUFunction.SetVParam(aP2d.Y());
+
+  // Compute the integral from myUMin to aUMax of myUFunction.
+  Standard_Integer i;
+  Standard_Real    aCoeff        = aV2d.Y();
+  Standard_Integer aNbUIntervals = anUKnots->Length() - 1;
+  //Standard_Real    aTol          = myTolerance/aNbUIntervals;
+  Standard_Real    aTol          = myTolerance;
+
+  //aTol /= myNbPntOuter;
+
+  if (myValueType == GProp_Mass) {
+    if (myIsByPoint)
+      aCoeff /= 3.;
+  } else if (myValueType == GProp_CenterMassX ||
+	     myValueType == GProp_CenterMassY ||
+	     myValueType == GProp_CenterMassZ) {
+    if (myIsByPoint)
+      aCoeff *= 0.25;
+  } else if (myValueType == GProp_InertiaXX ||
+	     myValueType == GProp_InertiaYY ||
+	     myValueType == GProp_InertiaZZ ||
+	     myValueType == GProp_InertiaXY ||
+	     myValueType == GProp_InertiaXZ ||
+	     myValueType == GProp_InertiaYZ) {
+    if (myIsByPoint)
+      aCoeff *= 0.2;
+  } else
+    return Standard_False;
+
+  Standard_Real aAbsCoeff = Abs(aCoeff);
+
+  if (aAbsCoeff <= Precision::Angular()) {
+    // No need to compute the integral. The value will be equal to 0.
+    F = 0.;
+    return Standard_True;
+  } 
+  //else if (aAbsCoeff > 10.*aTol)
+  //  aTol /= aAbsCoeff;
+  //else
+  //  aTol = 0.1;
+
+  Standard_Integer              iU           = anUKnots->Upper();
+  Standard_Integer              aNbPntsStart;
+  Standard_Integer              aNbMaxIter   = 1000;
+  math_KronrodSingleIntegration anIntegral;
+  Standard_Real                 aLocalErr    = 0.;
+
+  i            = anUKnots->Lower();
+  F            = 0.;
+ 
+  // Epmirical criterion
+  aNbPntsStart = Min(15, mySurface.UIntegrationOrder()/(anUKnots->Length() - 1)+1);
+  aNbPntsStart = Max(5, aNbPntsStart);
+
+   while (i < iU) {
+    Standard_Real aU1 = anUKnots->Value(i++);
+    Standard_Real aU2 = anUKnots->Value(i);
+
+    if(aU2 - aU1 < tolU) continue;
+
+    anIntegral.Perform(myUFunction, aU1, aU2, aNbPntsStart, aTol, aNbMaxIter);
+
+    if (!anIntegral.IsDone())
+      return Standard_False;
+
+    F         += anIntegral.Value();
+    aLocalErr += anIntegral.AbsolutError();
+    //cout << " TFunction : " << anIntegral.NbIterReached() << " " << anIntegral.AbsolutError() << endl;
+  }
+
+  F            *= aCoeff;
+  aLocalErr *= aAbsCoeff;
+  myAbsError = Max(myAbsError, aLocalErr);
+
+  myTolReached += aLocalErr;
+
+  if(Abs(F) > Epsilon(1.)) aLocalErr /= Abs(F);
+
+  myErrReached = Max(myErrReached, aLocalErr);
+
+
+  return Standard_True;
+}
+
+//=======================================================================
+//function : GetStateNumber
+//purpose  : 
+//=======================================================================
+
+Standard_Integer GProp_TFunction::GetStateNumber()
+{
+  //myErrReached  = myTolReached;
+  //myTolReached  = 0.;
+  //myNbPntOuter += RealToInt(0.5*myNbPntOuter);
+
+  //if (myNbPntOuter%2 == 0)
+    //myNbPntOuter++;
+
+  return 0;
+}

src/GProp/GProp_TFunction.lxx

@@ -0,0 +1,57 @@
+// File:	GProp_TFunction.lxx
+// Created:	Wed Dec 21 10:47:36 2005
+// Author:	Sergey KHROMOV
+//		<skv@dimox>
+
+
+//=======================================================================
+//function : SetNbKronrodPoints
+//purpose  : 
+//=======================================================================
+
+inline void GProp_TFunction::SetNbKronrodPoints
+                   (const Standard_Integer theNbPoints)
+{
+  myNbPntOuter = (theNbPoints%2 == 0) ? theNbPoints + 1 : theNbPoints;
+}
+ 
+//=======================================================================
+//function : SetValueType
+//purpose  : 
+//=======================================================================
+
+inline void GProp_TFunction::SetValueType(const GProp_ValueType theType)
+{
+  myValueType = theType;
+  myUFunction.SetValueType(myValueType);
+}
+
+//=======================================================================
+//function : SetTolerance
+//purpose  : 
+//=======================================================================
+
+inline void GProp_TFunction::SetTolerance(const Standard_Real theTolerance)
+{
+  myTolerance = theTolerance;
+}
+
+//=======================================================================
+//function : TolReached
+//purpose  : 
+//=======================================================================
+
+inline Standard_Real GProp_TFunction::ErrorReached() const
+{
+  return myErrReached;
+}
+
+//=======================================================================
+//function : ErrorReached
+//purpose  : 
+//=======================================================================
+
+inline Standard_Real GProp_TFunction::AbsolutError() const
+{
+  return myAbsError;
+}

src/GProp/GProp_UFunction.gxx

@@ -0,0 +1,257 @@
+// File:	GProp_UFunction.gxx
+// Created:	Fri Dec  9 16:18:05 2005
+// Author:	Sergey KHROMOV
+//		<skv@dimox>
+
+
+//=======================================================================
+//function : Constructor.
+//purpose  : 
+//=======================================================================
+
+GProp_UFunction::GProp_UFunction(const Face             &theSurface,
+				 const gp_Pnt           &theVertex,
+				 const Standard_Boolean  IsByPoint,
+				 const Standard_Address  theCoeffs)
+     : mySurface(theSurface),
+       myVertex(theVertex),
+       myCoeffs(theCoeffs),
+       myVParam(0.),
+       myValueType(GProp_Unknown),
+       myIsByPoint(IsByPoint)
+{
+}
+
+//=======================================================================
+//function : Value
+//purpose  : Returns a value of the function.
+//=======================================================================
+
+Standard_Boolean GProp_UFunction::Value(const Standard_Real  X,
+					      Standard_Real &F)
+{
+  // Volume computation
+  if (myValueType == GProp_Mass) {
+    gp_XYZ        aPMP0;
+    Standard_Real aTmpPar1;
+    Standard_Real aTmpPar2;
+
+    F = VolumeValue(X, aPMP0, aTmpPar1, aTmpPar2);
+
+    return Standard_True;
+  }
+
+  // Center of mass computation
+  if (myValueType == GProp_CenterMassX ||
+      myValueType == GProp_CenterMassY ||
+      myValueType == GProp_CenterMassZ)
+    return CenterMassValue(X, F);
+
+  // Inertia computation
+  if (myValueType == GProp_InertiaXX ||
+      myValueType == GProp_InertiaYY ||
+      myValueType == GProp_InertiaZZ ||
+      myValueType == GProp_InertiaXY ||
+      myValueType == GProp_InertiaXZ ||
+      myValueType == GProp_InertiaYZ)
+    return InertiaValue(X, F);
+
+  return Standard_False;
+}
+
+//=======================================================================
+//function : VolumeValue
+//purpose  : Returns the value for volume computation.
+//=======================================================================
+
+Standard_Real GProp_UFunction::VolumeValue(const Standard_Real  X,
+					         gp_XYZ        &thePMP0,
+					         Standard_Real &theS,
+					         Standard_Real &theD1)
+{
+  gp_Pnt aPnt;
+  gp_Vec aNorm;
+
+  mySurface.Normal(X, myVParam, aPnt, aNorm);
+
+  thePMP0 = aPnt.XYZ().Subtracted(myVertex.XYZ());
+
+  // Volume computation for ByPoint mode.
+  if (myIsByPoint)
+    return thePMP0.Dot(aNorm.XYZ());
+
+  // Volume and additional coefficients computation for ByPlane mode.
+  Standard_Real *aCoeff = (Standard_Real *)myCoeffs;
+
+  theS  =   aNorm.X()*aCoeff[0] + aNorm.Y()*aCoeff[1] + aNorm.Z()*aCoeff[2];
+  theD1 =   thePMP0.X()*aCoeff[0] + thePMP0.Y()*aCoeff[1]
+          + thePMP0.Z()*aCoeff[2] - aCoeff[3];
+
+  return theS*theD1;
+}
+
+//=======================================================================
+//function : CenterMassValue
+//purpose  : Returns a value for the center of mass computation.
+//=======================================================================
+
+Standard_Boolean GProp_UFunction::CenterMassValue(const Standard_Real  X,
+						        Standard_Real &F)
+{
+  gp_XYZ        aPmP0;
+  Standard_Real aS;
+  Standard_Real aD1;
+
+  F = VolumeValue(X, aPmP0, aS, aD1);
+
+  // Center of mass computation for ByPoint mode.
+  if (myIsByPoint) {
+    switch (myValueType) {
+    case GProp_CenterMassX:   F *= aPmP0.X(); break;
+    case GProp_CenterMassY:   F *= aPmP0.Y(); break;
+    case GProp_CenterMassZ:   F *= aPmP0.Z(); break;
+    default:
+      return Standard_False;
+    }
+
+    return Standard_True;
+  }
+
+  // Center of mass computation for ByPlane mode.
+  Standard_Real *aCoeff = (Standard_Real *)myCoeffs;
+
+  switch (myValueType) {
+  case GProp_CenterMassX:   F *= (aPmP0.X() - 0.5*aCoeff[0]*aD1); break;
+  case GProp_CenterMassY:   F *= (aPmP0.Y() - 0.5*aCoeff[1]*aD1); break;
+  case GProp_CenterMassZ:   F *= (aPmP0.Z() - 0.5*aCoeff[2]*aD1); break;
+  default:
+    return Standard_False;
+  }
+
+  return Standard_True;
+}
+
+//=======================================================================
+//function : InertiaValue
+//purpose  : Compute the value of intertia.
+//=======================================================================
+
+Standard_Boolean GProp_UFunction::InertiaValue(const Standard_Real  X,
+					             Standard_Real &F)
+{
+  gp_XYZ        aPmP0;
+  Standard_Real aS;
+  Standard_Real aD1;
+  Standard_Real aParam1;
+  Standard_Real aParam2;
+  Standard_Real *aCoeffs = (Standard_Real *)myCoeffs;
+
+  F = VolumeValue(X, aPmP0, aS, aD1);
+
+  // Inertia computation for ByPoint mode.
+  if (myIsByPoint) {
+    switch(myValueType) {
+    case GProp_InertiaXX:
+    case GProp_InertiaYZ:
+      aParam1 = aPmP0.Y() - aCoeffs[1];
+      aParam2 = aPmP0.Z() - aCoeffs[2];
+      break;
+    case GProp_InertiaYY:
+    case GProp_InertiaXZ:
+      aParam1 = aPmP0.X() - aCoeffs[0];
+      aParam2 = aPmP0.Z() - aCoeffs[2];
+      break;
+    case GProp_InertiaZZ:
+    case GProp_InertiaXY:
+      aParam1 = aPmP0.X() - aCoeffs[0];
+      aParam2 = aPmP0.Y() - aCoeffs[1];
+      break;
+    default:
+      return Standard_False;
+    }
+
+    if (myValueType == GProp_InertiaXX ||
+	myValueType == GProp_InertiaYY ||
+	myValueType == GProp_InertiaZZ)
+      F *=  aParam1*aParam1 + aParam2*aParam2;
+    else
+      F *= -aParam1*aParam2;
+
+    return Standard_True;
+  }
+
+  // Inertia computation for ByPlane mode.
+  Standard_Real aD2 = aD1*aD1;
+  Standard_Real aD3 = aD1*aD2/3.;
+  Standard_Real aPPar1;
+  Standard_Real aPPar2;
+  Standard_Real aCoeff1;
+  Standard_Real aCoeff2;
+
+  // Inertia computation for XX, YY and ZZ.
+  if (myValueType == GProp_InertiaXX ||
+      myValueType == GProp_InertiaYY ||
+      myValueType == GProp_InertiaZZ) {
+
+    if (myValueType == GProp_InertiaXX) {
+      aPPar1  = aPmP0.Y();
+      aPPar2  = aPmP0.Z();
+      aCoeff1 = aCoeffs[1];
+      aCoeff2 = aCoeffs[2];
+    } else if (myValueType == GProp_InertiaYY) {
+      aPPar1  = aPmP0.X();
+      aPPar2  = aPmP0.Z();
+      aCoeff1 = aCoeffs[0];
+      aCoeff2 = aCoeffs[2];
+    } else { // myValueType == GProp_InertiaZZ
+      aPPar1  = aPmP0.X();
+      aPPar2  = aPmP0.Y();
+      aCoeff1 = aCoeffs[0];
+      aCoeff2 = aCoeffs[1];
+    }
+
+    aPPar1  -= aCoeff1*aD1;
+    aPPar2  -= aCoeff2*aD1;
+    aParam1  = aPPar1*aPPar1*aD1 + aPPar1*aCoeff1*aD2 + aCoeff1*aCoeff1*aD3;
+    aParam2  = aPPar2*aPPar2*aD1 + aPPar2*aCoeff2*aD2 + aCoeff2*aCoeff2*aD3;
+
+    F = (aParam1 + aParam2)*aS;
+
+    return Standard_True;
+  }
+
+  // Inertia computation for XY, YZ and XZ.
+  if (myValueType == GProp_InertiaXY ||
+      myValueType == GProp_InertiaYZ ||
+      myValueType == GProp_InertiaXZ) {
+
+    if (myValueType == GProp_InertiaXY) {
+      aPPar1  = aPmP0.X();
+      aPPar2  = aPmP0.Y();
+      aCoeff1 = aCoeffs[0];
+      aCoeff2 = aCoeffs[1];
+    } else if (myValueType == GProp_InertiaYZ) {
+      aPPar1  = aPmP0.Y();
+      aPPar2  = aPmP0.Z();
+      aCoeff1 = aCoeffs[1];
+      aCoeff2 = aCoeffs[2];
+    } else { // myValueType == GProp_InertiaXZ
+      aPPar1  = aPmP0.X();
+      aPPar2  = aPmP0.Z();
+      aCoeff1 = aCoeffs[0];
+      aCoeff2 = aCoeffs[2];
+    }
+
+    aD2     *=   0.5;
+    aPPar1  -=   aCoeff1*aD1;
+    aPPar2  -=   aCoeff2*aD1;
+    aParam1  =   aPPar1*aPPar2*aD1
+               + (aPPar1*aCoeff2 + aPPar2*aCoeff1)*aD2 + aCoeff1*aCoeff2*aD3;
+
+    F = -aParam1*aS;
+
+    return Standard_True;
+  }
+
+  return Standard_False;
+}

src/GProp/GProp_UFunction.lxx

@@ -0,0 +1,27 @@
+// File:	GProp_UFunction.lxx
+// Created:	Wed Dec 21 10:00:26 2005
+// Author:	Sergey KHROMOV
+//		<skv@dimox>
+
+
+
+//=======================================================================
+//function : SetValueType
+//purpose  : Setting the type of the value to be returned.
+//=======================================================================
+
+inline void GProp_UFunction::SetValueType(const GProp_ValueType theType)
+{
+  myValueType = theType;
+}
+
+//=======================================================================
+//function : SetVParam
+//purpose  : Setting the V parameter that is constant during the 
+//           integral computation.
+//=======================================================================
+
+inline void GProp_UFunction::SetVParam(const Standard_Real theVParam)
+{
+  myVParam = theVParam;
+}

src/GProp/GProp_VGProps.cdl

@@ -0,0 +1,184 @@
+-- File:	VGProps.cdl<2>
+-- Created:	Fri Apr 12 10:01:47 1991
+-- Author:	Michel CHAUVAT
+--		Jean-Claude VAUTHIER January 1992
+---Copyright:	 Matra Datavision 1992
+
+
+generic class VGProps from GProp (Arc as any;
+                                  Face as any;  -- as FaceTool(Arc)
+    	    	    	    	  Domain as any -- as DomainTool(Arc)
+    	    	    	    	 )
+inherits GProps
+
+        --- Purpose :
+        --  Computes the global properties of a geometric solid 
+        --  (3D closed region of space) delimited with :
+        --  . a surface
+        --  . a point and a surface
+        --  . a plane and a surface
+        --  
+        --  The surface can be :
+        --  . a surface limited with its parametric values U-V,
+        --  . a surface limited in U-V space with its curves of restriction,
+        --  
+        --  The surface 's requirements to evaluate the global properties
+        --  are defined in the template SurfaceTool from package GProp.
+
+uses   	Pnt      from gp,
+        Pln      from gp
+is
+
+  Create returns VGProps;
+  
+  Create (S: Face; VLocation: Pnt from gp) returns VGProps;
+        --- Purpose : 
+        --  Computes the global properties of a region of 3D space
+        --  delimited with the surface <S> and the point VLocation. S can be closed 
+	--  The method is quick and its precision is enough for many cases of analytical 
+    	--  surfaces. 
+        --  Non-adaptive 2D Gauss integration with predefined numbers of Gauss points 
+	--  is used. Numbers of points depend on types of surfaces and  curves.
+    	--  Errror of the computation is not calculated. 
+	
+  Create (S: in out Face; VLocation: Pnt from gp; Eps: Real) returns VGProps; 
+        --- Purpose : 
+        --  Computes the global properties of a region of 3D space
+        --  delimited with the surface <S> and the point VLocation. S can be closed 
+	--  Adaptive 2D Gauss integration is used. 
+	--  Parameter Eps sets maximal relative error of computed mass (volume) for face. 
+	--  Error is calculated as Abs((M(i+1)-M(i))/M(i+1)), M(i+1) and M(i) are values 
+	--  for two successive steps of adaptive integration. 
+
+  Create (S: Face; O: Pnt from gp; VLocation: Pnt from gp) returns VGProps;
+        --- Purpose :  
+        --  Computes the global properties of the region of 3D space
+        --  delimited with the surface <S> and the point VLocation.
+	--  The method is quick and its precision is enough for many cases of analytical 
+    	--  surfaces. 
+        --  Non-adaptive 2D Gauss integration with predefined numbers of Gauss points 
+	--  is used. Numbers of points depend on types of surfaces and  curves.
+    	--  Error of the computation is not calculated. 
+	
+  Create (S: in out Face; O: Pnt from gp; VLocation: Pnt from gp; Eps: Real) returns VGProps;
+        --- Purpose :  
+        --  Computes the global properties of the region of 3D space
+        --  delimited with the surface <S> and the point VLocation.
+	--  Adaptive 2D Gauss integration is used. 
+	--  Parameter Eps sets maximal relative error of computed mass (volume) for face. 
+	--  Error is calculated as Abs((M(i+1)-M(i))/M(i+1)), M(i+1) and M(i) are values 
+	--  for two successive steps of adaptive integration. 
+	--          WARNING: if Eps > 0.001 algorithm performs non-adaptive integration.	  
+     
+  Create (S: Face; Pl: Pln from gp; VLocation: Pnt from gp) returns VGProps;
+        --- Purpose :  
+        --  Computes the global properties of the region of 3D space
+        --  delimited with the surface <S> and the plane Pln.
+	--  The method is quick and its precision is enough for many cases of analytical 
+    	--  surfaces.  
+        --  Non-adaptive 2D Gauss integration with predefined numbers of Gauss points 
+	--  is used. Numbers of points depend on types of surfaces and  curves.
+    	--  Error of the computation is not calculated. 
+	
+  Create (S: in out Face; Pl: Pln from gp; VLocation: Pnt from gp; Eps: Real) returns VGProps;
+        --- Purpose :  
+        --  Computes the global properties of the region of 3D space
+        --  delimited with the surface <S> and the plane Pln.
+	--  Adaptive 2D Gauss integration is used. 
+	--  Parameter Eps sets maximal relative error of computed mass (volume) for face. 
+	--  Error is calculated as Abs((M(i+1)-M(i))/M(i+1)), M(i+1) and M(i) are values 
+	--  for two successive steps of adaptive integration. 
+	--          WARNING: if Eps > 0.001 algorithm performs non-adaptive integration.	  
+
+    	--  With  Domain  -- 
+
+  Create (S: in out Face; D : in out Domain; VLocation: Pnt from gp) returns VGProps;
+        --- Purpose : 
+        --  Computes the global properties of a region of 3D space
+        --  delimited with the surface <S> and the point VLocation. S can be closed 
+	--  The method is quick and its precision is enough for many cases of analytical 
+    	--  surfaces. 
+        --  Non-adaptive 2D Gauss integration with predefined numbers of Gauss points 
+	--  is used. Numbers of points depend on types of surfaces and  curves.
+    	--  Errror of the computation is not calculated. 
+	
+  Create (S: in out Face; D : in out Domain; VLocation: Pnt from gp; Eps: Real) returns VGProps; 
+        --- Purpose : 
+        --  Computes the global properties of a region of 3D space
+        --  delimited with the surface <S> and the point VLocation. S can be closed 
+	--  Adaptive 2D Gauss integration is used. 
+	--  Parameter Eps sets maximal relative error of computed mass (volume) for face. 
+	--  Error is calculated as Abs((M(i+1)-M(i))/M(i+1)), M(i+1) and M(i) are values 
+	--  for two successive steps of adaptive integration. 
+
+  Create (S: in out Face; D : in out Domain; O: Pnt from gp; VLocation: Pnt from gp) returns VGProps;
+        --- Purpose :  
+        --  Computes the global properties of the region of 3D space
+        --  delimited with the surface <S> and the point VLocation.
+	--  The method is quick and its precision is enough for many cases of analytical 
+    	--  surfaces. 
+        --  Non-adaptive 2D Gauss integration with predefined numbers of Gauss points 
+	--  is used. Numbers of points depend on types of surfaces and  curves.
+    	--  Error of the computation is not calculated. 
+	
+  Create (S: in out Face; D : in out Domain; O: Pnt from gp; VLocation: Pnt from gp; Eps: Real) returns VGProps;
+        --- Purpose :  
+        --  Computes the global properties of the region of 3D space
+        --  delimited with the surface <S> and the point VLocation.
+	--  Adaptive 2D Gauss integration is used. 
+	--  Parameter Eps sets maximal relative error of computed mass (volume) for face. 
+	--  Error is calculated as Abs((M(i+1)-M(i))/M(i+1)), M(i+1) and M(i) are values 
+	--  for two successive steps of adaptive integration. 
+	--          WARNING: if Eps > 0.001 algorithm performs non-adaptive integration.	  
+     
+  Create (S: in out Face; D : in out Domain; Pl: Pln from gp; VLocation: Pnt from gp) returns VGProps;
+        --- Purpose :  
+        --  Computes the global properties of the region of 3D space
+        --  delimited with the surface <S> and the plane Pln.
+	--  The method is quick and its precision is enough for many cases of analytical 
+    	--  surfaces.  
+        --  Non-adaptive 2D Gauss integration with predefined numbers of Gauss points 
+	--  is used. Numbers of points depend on types of surfaces and  curves.
+    	--  Error of the computation is not calculated. 
+	
+  Create (S: in out Face; D : in out Domain; Pl: Pln from gp; VLocation: Pnt from gp; Eps: Real) returns VGProps;
+        --- Purpose :  
+        --  Computes the global properties of the region of 3D space
+        --  delimited with the surface <S> and the plane Pln.
+	--  Adaptive 2D Gauss integration is used. 
+	--  Parameter Eps sets maximal relative error of computed mass (volume) for face. 
+	--  Error is calculated as Abs((M(i+1)-M(i))/M(i+1)), M(i+1) and M(i) are values 
+	--  for two successive steps of adaptive integration. 
+	--          WARNING: if Eps > 0.001 algorithm performs non-adaptive integration.	  
+
+  SetLocation(me: in out; VLocation: Pnt from gp);
+  
+  Perform(me: in out; S: Face);
+  Perform(me: in out; S: in out Face; Eps: Real) returns Real;
+
+  Perform(me: in out; S: Face; O : Pnt from gp);
+  Perform(me: in out; S: in out Face; O : Pnt from gp; Eps: Real) returns Real;
+
+  Perform(me: in out; S: Face; Pl : Pln from gp);
+  Perform(me: in out; S: in out Face; Pl : Pln from gp; Eps: Real) returns Real;
+
+  Perform(me: in out; S: in out  Face; D : in out Domain);
+  Perform(me: in out; S: in out Face; D : in out Domain;  Eps: Real) returns Real;
+
+  Perform(me: in out; S: in out Face; D : in out Domain;  O : Pnt from gp);
+  Perform(me: in out; S: in out Face; D : in out Domain;  O : Pnt from gp; Eps: Real) returns Real;
+
+  Perform(me: in out; S: in out Face; D : in out Domain;  Pl : Pln from gp);
+  Perform(me: in out; S: in out Face; D : in out Domain;  Pl : Pln from gp; Eps: Real) returns Real;
+
+  GetEpsilon(me: out) returns Real;
+        --- Purpose :  
+        --  If previously used methods containe Eps parameter  
+    	--  gets actual relative error of the computation, else returns  1.0.
+fields
+ 
+    myEpsilon: Real from Standard;
+
+end VGProps;
+
+

src/GProp/GProp_VGProps.gxx

@@ -0,0 +1,947 @@
+#include <Standard_NotImplemented.hxx>
+#include <math_Vector.hxx>
+#include <math.hxx>
+#include <gp_Pnt2d.hxx>
+#include <gp_Vec2d.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Vec.hxx>
+
+#include <TColStd_Array1OfReal.hxx>
+#include <Precision.hxx>
+class HMath_Vector{
+  math_Vector *pvec;
+  void operator=(const math_Vector&){}
+ public:
+  HMath_Vector(){ pvec = 0;}
+  HMath_Vector(math_Vector* pv){ pvec = pv;}
+  ~HMath_Vector(){ if(pvec != 0) delete pvec;}
+  void operator=(math_Vector* pv){ if(pvec != pv && pvec != 0) delete pvec;  pvec = pv;}
+  Standard_Real& operator()(Standard_Integer i){ return (*pvec).operator()(i);}
+  const Standard_Real& operator()(Standard_Integer i) const{ return (*pvec).operator()(i);}
+  const math_Vector* operator->() const{ return pvec;}
+  math_Vector* operator->(){ return pvec;}
+  math_Vector* Init(Standard_Real v, Standard_Integer i = 0, Standard_Integer iEnd = 0){ 
+    if(pvec == 0) return pvec;
+    if(iEnd - i == 0) pvec->Init(v);
+    else for(; i <= iEnd; i++) pvec->operator()(i) = v;
+    return pvec;
+  }
+};
+
+//Minimal value of interval's range for computation | minimal value of "dim" | ... 
+static Standard_Real     EPS_PARAM = Precision::Angular(), EPS_DIM = 1.E-30, ERROR_ALGEBR_RATIO = 2.0/3.0;
+//Maximum of GaussPoints on a subinterval and maximum of subintervals
+static Standard_Integer  GPM = math::GaussPointsMax(), SUBS_POWER = 32, SM = SUBS_POWER*GPM + 1; 
+static Standard_Boolean  IS_MIN_DIM = 1; // if the value equal 0 error of algorithm calculted by static moments
+
+static math_Vector  LGaussP0(1,GPM), LGaussW0(1,GPM),
+  LGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM))), LGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
+static HMath_Vector  L1 = new math_Vector(1,SM), L2 = new math_Vector(1,SM),
+  DimL = new math_Vector(1,SM), ErrL = new math_Vector(1,SM), ErrUL = new math_Vector(1,SM,0.0),
+  IxL = new math_Vector(1,SM), IyL = new math_Vector(1,SM), IzL = new math_Vector(1,SM),
+  IxxL = new math_Vector(1,SM), IyyL = new math_Vector(1,SM), IzzL = new math_Vector(1,SM),
+  IxyL = new math_Vector(1,SM), IxzL = new math_Vector(1,SM), IyzL = new math_Vector(1,SM);
+
+static math_Vector* LGaussP[] = {&LGaussP0,&LGaussP1};
+static math_Vector* LGaussW[] = {&LGaussW0,&LGaussW1};
+
+static math_Vector  UGaussP0(1,GPM), UGaussW0(1,GPM),
+  UGaussP1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM))), UGaussW1(1,RealToInt(Ceiling(ERROR_ALGEBR_RATIO*GPM)));
+static HMath_Vector  U1 = new math_Vector(1,SM), U2 = new math_Vector(1,SM), 
+  DimU = new math_Vector(1,SM), ErrU = new math_Vector(1,SM,0.0),
+  IxU = new math_Vector(1,SM), IyU = new math_Vector(1,SM), IzU = new math_Vector(1,SM),
+  IxxU = new math_Vector(1,SM), IyyU = new math_Vector(1,SM), IzzU = new math_Vector(1,SM),
+  IxyU = new math_Vector(1,SM), IxzU = new math_Vector(1,SM), IyzU = new math_Vector(1,SM);
+
+static math_Vector* UGaussP[] = {&UGaussP0,&UGaussP1};
+static math_Vector* UGaussW[] = {&UGaussW0,&UGaussW1};
+
+static Standard_Integer FillIntervalBounds(Standard_Real A, Standard_Real B, const TColStd_Array1OfReal& Knots, 
+					   HMath_Vector& VA, HMath_Vector& VB)
+{
+  Standard_Integer i = 1, iEnd = Knots.Upper(), j = 1, k = 1;
+  VA(j++) = A;
+  for(; i <= iEnd; i++){
+    Standard_Real kn = Knots(i);
+    if(A < kn)
+      if(kn < B) VA(j++) = VB(k++) = kn; else break;
+  }
+  VB(k) = B;
+  return k;
+}
+
+static inline Standard_Integer MaxSubs(Standard_Integer n, Standard_Integer coeff = SUBS_POWER){
+  return n = IntegerLast()/coeff < n? IntegerLast(): n*coeff + 1; 
+}
+
+static Standard_Integer LFillIntervalBounds(Standard_Real A, Standard_Real B, const TColStd_Array1OfReal& Knots, 
+					    const Standard_Integer NumSubs)
+{
+  Standard_Integer iEnd = Knots.Upper(), jEnd = L1->Upper();
+
+//  Modified by Sergey KHROMOV - Wed Mar 26 11:22:50 2003
+  iEnd = Max(iEnd, MaxSubs(iEnd-1,NumSubs));
+  if(iEnd - 1 > jEnd){ 
+//     iEnd = MaxSubs(iEnd-1,NumSubs); 
+//  Modified by Sergey KHROMOV - Wed Mar 26 11:22:51 2003
+    L1 = new math_Vector(1,iEnd);  L2 = new math_Vector(1,iEnd);
+    DimL = new math_Vector(1,iEnd); ErrL = new math_Vector(1,iEnd,0.0); ErrUL = new math_Vector(1,iEnd,0.0);
+    IxL = new math_Vector(1,iEnd); IyL = new math_Vector(1,iEnd); IzL = new math_Vector(1,iEnd);
+    IxxL = new math_Vector(1,iEnd); IyyL = new math_Vector(1,iEnd); IzzL = new math_Vector(1,iEnd);
+    IxyL = new math_Vector(1,iEnd); IxzL = new math_Vector(1,iEnd); IyzL = new math_Vector(1,iEnd);
+  }
+  return FillIntervalBounds(A, B, Knots, L1, L2);
+}
+
+static Standard_Integer UFillIntervalBounds(Standard_Real A, Standard_Real B, const TColStd_Array1OfReal& Knots,
+					    const Standard_Integer NumSubs)
+{
+  Standard_Integer iEnd = Knots.Upper(), jEnd = U1->Upper();
+
+//  Modified by Sergey KHROMOV - Wed Mar 26 11:22:50 2003
+  iEnd = Max(iEnd, MaxSubs(iEnd-1,NumSubs));
+  if(iEnd - 1 > jEnd){
+//     iEnd = MaxSubs(iEnd-1,NumSubs); 
+//  Modified by Sergey KHROMOV - Wed Mar 26 11:22:51 2003
+    U1 = new math_Vector(1,iEnd);  U2 = new math_Vector(1,iEnd);
+    DimU = new math_Vector(1,iEnd); ErrU = new math_Vector(1,iEnd,0.0);
+    IxU = new math_Vector(1,iEnd); IyU = new math_Vector(1,iEnd); IzU = new math_Vector(1,iEnd);
+    IxxU = new math_Vector(1,iEnd); IyyU = new math_Vector(1,iEnd); IzzU = new math_Vector(1,iEnd);
+    IxyU = new math_Vector(1,iEnd); IxzU = new math_Vector(1,iEnd); IyzU = new math_Vector(1,iEnd);
+  }
+  return FillIntervalBounds(A, B, Knots, U1, U2);
+}
+
+static Standard_Real CCompute(Face& S, Domain& D, const Standard_Boolean ByPoint, const Standard_Real Coeff[],
+			      const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, 
+			      const Standard_Real EpsDim, 
+			      const Standard_Boolean isErrorCalculation, const Standard_Boolean isVerifyComputation) 
+{
+  Standard_Boolean isNaturalRestriction = S.NaturalRestriction();
+
+  Standard_Integer NumSubs = SUBS_POWER;
+  Standard_Boolean isMinDim = IS_MIN_DIM;
+
+  Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
+  Dim = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
+  //boundary curve parametrization
+  Standard_Real l1, l2, lm, lr, l;   
+  //Face parametrization in U and V direction
+  Standard_Real BV1, BV2, v;         
+  Standard_Real BU1, BU2, u1, u2, um, ur, u;
+  S.Bounds (BU1, BU2, BV1, BV2);  u1 = BU1;
+  //location point used to compute the inertia
+  Standard_Real xloc, yloc, zloc;
+  loc.Coord (xloc, yloc, zloc);
+  //location point used to compute the inertiard (xloc, yloc, zloc);
+  //Jacobien (x, y, z) -> (u, v) = ||n||
+  Standard_Real xn, yn, zn, s, ds, dDim;
+  Standard_Real x, y, z, xi, px, py, pz, yi, zi, d1, d2, d3;
+  //On the Face
+  gp_Pnt Ps;                    
+  gp_Vec VNor;
+  //On the boundary curve u-v
+  gp_Pnt2d Puv;                
+  gp_Vec2d Vuv;
+  Standard_Real Dul;  // Dul = Du / Dl
+  Standard_Real CDim[2], CIx, CIy, CIz, CIxx[2], CIyy[2], CIzz[2], CIxy, CIxz, CIyz;
+  Standard_Real LocDim[2], LocIx[2], LocIy[2], LocIz[2], LocIxx[2], LocIyy[2], LocIzz[2], LocIxy[2], LocIxz[2], LocIyz[2];
+  
+  Standard_Integer iD = 0, NbLSubs, iLS, iLSubEnd, iGL, iGLEnd, NbLGaussP[2], LRange[2], iL, kL, kLEnd, IL, JL;
+  Standard_Integer i, NbUSubs, iUS, iUSubEnd, iGU, iGUEnd, NbUGaussP[2], URange[2], iU, kU, kUEnd, IU, JU;
+  Standard_Integer UMaxSubs, LMaxSubs;
+
+  Standard_Real ErrorU, ErrorL, ErrorLMax = 0.0, Eps=0.0, EpsL=0.0, EpsU=0.0;
+  iGLEnd = isErrorCalculation? 2: 1; 
+
+  for(i = 0; i < 2; i++) {
+    LocDim[i] = 0.0;
+    LocIx[i] = 0.0;
+    LocIy[i] = 0.0;
+    LocIz[i] = 0.0;
+    LocIxx[i] = 0.0;
+    LocIyy[i] = 0.0;
+    LocIzz[i] = 0.0;
+    LocIxy[i] = 0.0;
+    LocIyz[i] = 0.0;
+    LocIxz[i] = 0.0;
+  }
+
+  NbUGaussP[0] = S.SIntOrder(EpsDim);  
+  NbUGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbUGaussP[0]));
+  math::GaussPoints(NbUGaussP[0],UGaussP0);  math::GaussWeights(NbUGaussP[0],UGaussW0);
+  math::GaussPoints(NbUGaussP[1],UGaussP1);  math::GaussWeights(NbUGaussP[1],UGaussW1);
+  
+  NbUSubs = S.SUIntSubs();
+  TColStd_Array1OfReal UKnots(1,NbUSubs+1);
+  S.UKnots(UKnots);
+  
+  while (isNaturalRestriction || D.More()) {
+    if(isNaturalRestriction){ 
+      NbLGaussP[0] = Min(2*NbUGaussP[0],math::GaussPointsMax());
+    }else{
+      S.Load(D.Value());  ++iD;
+      NbLGaussP[0] = S.LIntOrder(EpsDim);  
+    }
+    NbLGaussP[1] = RealToInt(Ceiling(ERROR_ALGEBR_RATIO*NbLGaussP[0]));
+    math::GaussPoints(NbLGaussP[0],LGaussP0);  math::GaussWeights(NbLGaussP[0],LGaussW0);
+    math::GaussPoints(NbLGaussP[1],LGaussP1);  math::GaussWeights(NbLGaussP[1],LGaussW1);
+    
+    NbLSubs = isNaturalRestriction? S.SVIntSubs(): S.LIntSubs();
+    TColStd_Array1OfReal LKnots(1,NbLSubs+1);
+    if(isNaturalRestriction){
+      S.VKnots(LKnots); 
+      l1 = BV1; l2 = BV2;
+    }else{
+      S.LKnots(LKnots);
+      l1 = S.FirstParameter();  l2 = S.LastParameter();
+    }
+    ErrorL = 0.0;
+    kLEnd = 1; JL = 0;
+    //OCC503(apo): if(Abs(l2-l1) < EPS_PARAM) continue;
+    if(Abs(l2-l1) > EPS_PARAM) {
+      iLSubEnd = LFillIntervalBounds(l1, l2, LKnots, NumSubs);
+      LMaxSubs = MaxSubs(iLSubEnd);
+      //-- exception avoiding
+      if(LMaxSubs > SM) LMaxSubs = SM;
+      DimL.Init(0.0,1,LMaxSubs);  ErrL.Init(0.0,1,LMaxSubs);  ErrUL.Init(0.0,1,LMaxSubs);
+      do{// while: L
+	if(++JL > iLSubEnd){
+	  LRange[0] = IL = ErrL->Max();  LRange[1] = JL;
+	  L1(JL) = (L1(IL) + L2(IL))/2.0;  L2(JL) = L2(IL);  L2(IL) = L1(JL);
+	}else  LRange[0] = IL = JL;
+	if(JL == LMaxSubs || Abs(L2(JL) - L1(JL)) < EPS_PARAM)
+	  if(kLEnd == 1){
+	    DimL(JL) = ErrL(JL) = IxL(JL) = IyL(JL) = IzL(JL) = 
+	      IxxL(JL) = IyyL(JL) = IzzL(JL) = IxyL(JL) = IxzL(JL) = IyzL(JL) = 0.0;
+	  }else{
+	    JL--;
+	    EpsL = ErrorL;  Eps = EpsL/0.9;
+	    break;
+	  }
+	else
+	  for(kL=0; kL < kLEnd; kL++){
+	    iLS = LRange[kL];
+	    lm = 0.5*(L2(iLS) + L1(iLS));         
+	    lr = 0.5*(L2(iLS) - L1(iLS));
+	    CIx = CIy = CIz = CIxy = CIxz = CIyz = 0.0;
+	    for(iGL=0; iGL < iGLEnd; iGL++){//
+	      CDim[iGL] = CIxx[iGL] = CIyy[iGL] = CIzz[iGL] = 0.0;
+	      for(iL=1; iL<=NbLGaussP[iGL]; iL++){
+		l = lm + lr*(*LGaussP[iGL])(iL);
+		if(isNaturalRestriction){ 
+		  v = l; u2 = BU2; Dul = (*LGaussW[iGL])(iL);
+		}else{
+		  S.D12d (l, Puv, Vuv);
+		  Dul = Vuv.Y()*(*LGaussW[iGL])(iL);  // Dul = Du / Dl
+		  if(Abs(Dul) < EPS_PARAM) continue;
+		  v  = Puv.Y();  u2 = Puv.X();
+		  //Check on cause out off bounds of value current parameter
+		  if(v < BV1) v = BV1; else if(v > BV2) v = BV2;
+		  if(u2 < BU1) u2 = BU1; else if(u2 > BU2) u2 = BU2; 
+		}
+		ErrUL(iLS) = 0.0;
+		kUEnd = 1; JU = 0;
+		if(Abs(u2-u1) < EPS_PARAM) continue;
+		iUSubEnd = UFillIntervalBounds(u1, u2, UKnots, NumSubs);
+		UMaxSubs = MaxSubs(iUSubEnd);
+		//-- exception avoiding
+		if(UMaxSubs > SM) UMaxSubs = SM;
+		DimU.Init(0.0,1,UMaxSubs);  ErrU.Init(0.0,1,UMaxSubs);  ErrorU = 0.0;
+		do{//while: U
+		  if(++JU > iUSubEnd){
+		    URange[0] = IU = ErrU->Max();  URange[1] = JU;  
+		    U1(JU) = (U1(IU)+U2(IU))/2.0;  U2(JU) = U2(IU);  U2(IU) = U1(JU);
+		  }else  URange[0] = IU = JU;
+		  if(JU == UMaxSubs || Abs(U2(JU) - U1(JU)) < EPS_PARAM)
+		    if(kUEnd == 1){
+		      DimU(JU) = ErrU(JU) = IxU(JU) = IyU(JU) = IzU(JU) = 
+			IxxU(JU) = IyyU(JU) = IzzU(JU) = IxyU(JU) = IxzU(JU) = IyzU(JU) = 0.0;
+		    }else{
+		      JU--;  
+		      EpsU = ErrorU;  Eps = EpsU*Abs((u2-u1)*Dul)/0.1;  EpsL = 0.9*Eps;
+		      break;
+		    }
+		  else
+		    for(kU=0; kU < kUEnd; kU++){
+		      iUS = URange[kU];
+		      um = 0.5*(U2(iUS) + U1(iUS));
+		      ur = 0.5*(U2(iUS) - U1(iUS));
+		      iGUEnd = iGLEnd - iGL;
+		      for(iGU=0; iGU < iGUEnd; iGU++){//
+			LocDim[iGU] = 
+			  LocIxx[iGU] = LocIyy[iGU] = LocIzz[iGU] = 
+			    LocIx[iGU] = LocIy[iGU] = LocIz[iGU] = 
+			      LocIxy[iGU] = LocIxz[iGU] = LocIyz[iGU] = 0.0;
+			for(iU=1; iU<=NbUGaussP[iGU]; iU++){
+			  u = um + ur*(*UGaussP[iGU])(iU);
+			  S.Normal(u, v, Ps, VNor);
+			  VNor.Coord(xn, yn, zn);
+			  Ps.Coord(x, y, z);  
+			  x -= xloc;  y -= yloc;  z -= zloc;
+			  xn *= (*UGaussW[iGU])(iU);
+			  yn *= (*UGaussW[iGU])(iU);
+			  zn *= (*UGaussW[iGU])(iU);
+			  if(ByPoint){
+			    //volume of elementary cone
+			    dDim = (x*xn+y*yn+z*zn)/3.0;
+			    //coordinates of cone's center mass 
+			    px = 0.75*x;  py = 0.75*y;  pz = 0.75*z;  
+			    LocDim[iGU] +=  dDim;
+			    //if(iGU > 0) continue;
+			    LocIx[iGU] += px*dDim;  
+			    LocIy[iGU] += py*dDim;  
+			    LocIz[iGU] += pz*dDim;
+			    x -= Coeff[0];  y -= Coeff[1];  z -= Coeff[2];
+			    dDim *= 3.0/5.0;
+			    LocIxy[iGU] -= x*y*dDim;  
+			    LocIyz[iGU] -= y*z*dDim;  
+			    LocIxz[iGU] -= x*z*dDim;
+			    xi = x*x;  yi = y*y;  zi = z*z;
+			    LocIxx[iGU] += (yi+zi)*dDim;  
+			    LocIyy[iGU] += (xi+zi)*dDim;  
+			    LocIzz[iGU] += (xi+yi)*dDim;
+			  }else{ // by plane
+			    s = xn*Coeff[0] + yn*Coeff[1] + zn*Coeff[2];
+			    d1 = Coeff[0]*x + Coeff[1]*y + Coeff[2]*z - Coeff[3];
+			    d2 = d1*d1;
+			    d3 = d1*d2/3.0;
+			    ds = s*d1;
+			    LocDim[iGU] += ds;
+			    //if(iGU > 0) continue;
+			    LocIx[iGU] += (x - Coeff[0]*d1/2.0) * ds;  
+			    LocIy[iGU] += (y - Coeff[1]*d1/2.0) * ds;
+			    LocIz[iGU] += (z - Coeff[2]*d1/2.0) * ds;
+			    px = x-Coeff[0]*d1;  py = y-Coeff[1]*d1;  pz = z-Coeff[2]*d1;
+			    xi = px*px*d1 + px*Coeff[0]*d2 + Coeff[0]*Coeff[0]*d3;
+			    yi = py*py*d1 + py*Coeff[1]*d2 + Coeff[1]*Coeff[1]*d3;
+			    zi = pz*pz*d1 + pz*Coeff[2]*d2 + Coeff[2]*Coeff[2]*d3;
+			    LocIxx[iGU] += (yi+zi)*s;  
+			    LocIyy[iGU] += (xi+zi)*s;  
+			    LocIzz[iGU] += (xi+yi)*s;
+			    d2 /= 2.0;
+			    xi = py*pz*d1 + py*Coeff[2]*d2 + pz*Coeff[1]*d2 + Coeff[1]*Coeff[2]*d3;
+			    yi = px*pz*d1 + pz*Coeff[0]*d2 + px*Coeff[2]*d2 + Coeff[0]*Coeff[2]*d3;
+			    zi = px*py*d1 + px*Coeff[1]*d2 + py*Coeff[0]*d2 + Coeff[0]*Coeff[1]*d3;
+			    LocIxy[iGU] -= zi*s;  LocIyz[iGU] -= xi*s;  LocIxz[iGU] -= yi*s;
+			  }
+			}//for: iU
+		      }//for: iGU
+		      DimU(iUS) = LocDim[0]*ur;
+		      IxxU(iUS) = LocIxx[0]*ur; IyyU(iUS) = LocIyy[0]*ur; IzzU(iUS) = LocIzz[0]*ur;
+		      if(iGL > 0) continue;
+		      LocDim[1] = Abs(LocDim[1]-LocDim[0]); 
+		      LocIxx[1] = Abs(LocIxx[1]-LocIxx[0]); 
+		      LocIyy[1] = Abs(LocIyy[1]-LocIyy[0]); 
+		      LocIzz[1] = Abs(LocIzz[1]-LocIzz[0]);
+		      ErrU(iUS) = isMinDim? LocDim[1]*ur: (LocIxx[1] + LocIyy[1] + LocIzz[1])*ur;
+		      IxU(iUS) = LocIx[0]*ur; IyU(iUS) = LocIy[0]*ur; IzU(iUS) = LocIz[0]*ur;
+		      IxyU(iUS) = LocIxy[0]*ur; IxzU(iUS) = LocIxz[0]*ur; IyzU(iUS) = LocIyz[0]*ur;
+		    }//for: kU (iUS)
+		  if(JU == iUSubEnd)  kUEnd = 2; 
+		  if(kUEnd == 2)  {
+		    Standard_Integer imax = ErrU->Max();
+		    if(imax > 0) ErrorU = ErrU(imax);
+		    else ErrorU = 0.0;
+		  }
+		}while((ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1);
+		for(i=1; i<=JU; i++) {
+		  CDim[iGL] += DimU(i)*Dul;
+		  CIxx[iGL] += IxxU(i)*Dul; CIyy[iGL] += IyyU(i)*Dul; CIzz[iGL] += IzzU(i)*Dul;
+		}
+		if(iGL > 0) continue;
+		ErrUL(iLS) = ErrorU*Abs((u2-u1)*Dul);
+		for(i=1; i<=JU; i++){
+		  CIx += IxU(i)*Dul; CIy += IyU(i)*Dul; CIz += IzU(i)*Dul;
+		  //CIxx += IxxU(i)*Dul; CIyy += IyyU(i)*Dul; CIzz += IzzU(i)*Dul;
+		  CIxy += IxyU(i)*Dul; CIxz += IxzU(i)*Dul; CIyz += IyzU(i)*Dul;
+		}
+	      }//for: iL 
+	    }//for: iGL
+	    DimL(iLS) = CDim[0]*lr;  
+	    IxxL(iLS) = CIxx[0]*lr; IyyL(iLS) = CIyy[0]*lr; IzzL(iLS) = CIzz[0]*lr; 
+	    if(iGLEnd == 2) {
+	      //ErrL(iLS) = Abs(CDim[1]-CDim[0])*lr + ErrUL(iLS);
+	      CDim[1] = Abs(CDim[1]-CDim[0]);
+	      CIxx[1] = Abs(CIxx[1]-CIxx[0]); CIyy[1] = Abs(CIyy[1]-CIyy[0]); CIzz[1] = Abs(CIzz[1]-CIzz[0]);
+	      ErrorU = ErrUL(iLS);
+	      ErrL(iLS) = (isMinDim? CDim[1]: (CIxx[1] + CIyy[1] + CIzz[1]))*lr + ErrorU;
+	    }
+	    IxL(iLS) = CIx*lr; IyL(iLS) = CIy*lr; IzL(iLS) = CIz*lr; 
+	    //IxxL(iLS) = CIxx*lr; IyyL(iLS) = CIyy*lr; IzzL(iLS) = CIzz*lr; 
+	    IxyL(iLS) = CIxy*lr; IxzL(iLS) = CIxz*lr; IyzL(iLS) = CIyz*lr; 
+	  }//for: (kL)iLS
+	//  Calculate/correct epsilon of computation by current value of Dim
+	//That is need for not spend time for 
+	if(JL == iLSubEnd){  
+	  kLEnd = 2; 
+	  Standard_Real DDim = 0.0, DIxx = 0.0, DIyy = 0.0, DIzz = 0.0;
+	  for(i=1; i<=JL; i++) {
+	    DDim += DimL(i);
+	    DIxx += IxxL(i);  DIyy += IyyL(i);  DIzz += IzzL(i);
+	  }
+	  DDim = isMinDim? Abs(DDim): Abs(DIxx) + Abs(DIyy) + Abs(DIzz);
+	  DDim = Abs(DDim*EpsDim);
+	  if(DDim > Eps) { 
+	    Eps = DDim;  EpsL = 0.9*Eps;
+	  }
+	}
+	if(kLEnd == 2) {
+	  Standard_Integer imax = ErrL->Max();
+	  if(imax > 0) ErrorL = ErrL(imax);
+	  else ErrorL = 0.0;
+	}
+      }while((ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1);
+      for(i=1; i<=JL; i++){
+	Dim += DimL(i); 
+	Ix += IxL(i); Iy += IyL(i); Iz += IzL(i);
+	Ixx += IxxL(i); Iyy += IyyL(i); Izz += IzzL(i);
+      Ixy += IxyL(i); Ixz += IxzL(i); Iyz += IyzL(i);
+      }
+      ErrorLMax = Max(ErrorLMax, ErrorL);
+    }
+    if(isNaturalRestriction) break;
+    D.Next();
+  }
+  if(Abs(Dim) >= EPS_DIM){
+    if(ByPoint){
+      Ix = Coeff[0] + Ix/Dim;
+      Iy = Coeff[1] + Iy/Dim;
+      Iz = Coeff[2] + Iz/Dim;
+    }else{
+      Ix /= Dim;
+      Iy /= Dim;
+      Iz /= Dim;
+    }
+    g.SetCoord (Ix, Iy, Iz);
+  }else{
+    Dim =0.;
+    g.SetCoord(0.,0.,0.);
+  }
+  inertia.SetCols (gp_XYZ (Ixx, Ixy, Ixz),
+		   gp_XYZ (Ixy, Iyy, Iyz),
+		   gp_XYZ (Ixz, Iyz, Izz));
+  if(iGLEnd == 2) 
+    Eps = Dim != 0.0? ErrorLMax/(isMinDim? Abs(Dim): (Abs(Ixx) + Abs(Iyy) + Abs(Izz))): 0.0;
+  else Eps = EpsDim;
+  return Eps;
+}
+
+static Standard_Real Compute(Face& S, const Standard_Boolean ByPoint, const Standard_Real Coeff[],
+			     const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, Standard_Real EpsDim) 
+{
+  Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
+  Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
+  EpsDim = Abs(EpsDim);
+  Domain D;
+  return CCompute(S,D,ByPoint,Coeff,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
+}
+
+static Standard_Real Compute(Face& S, Domain& D, const Standard_Boolean ByPoint, const Standard_Real Coeff[],
+			     const gp_Pnt& loc, Standard_Real& Dim, gp_Pnt& g, gp_Mat& inertia, Standard_Real EpsDim) 
+{
+  Standard_Boolean isErrorCalculation  = 0.0 > EpsDim || EpsDim < 0.001? 1: 0;
+  Standard_Boolean isVerifyComputation = 0.0 < EpsDim && EpsDim < 0.001? 1: 0;
+  EpsDim = Abs(EpsDim);
+  return CCompute(S,D,ByPoint,Coeff,loc,Dim,g,inertia,EpsDim,isErrorCalculation,isVerifyComputation);
+}
+
+static void Compute(const Face& S,
+		    const Standard_Boolean ByPoint,
+		    const Standard_Real  Coeff[],
+                    const gp_Pnt& Loc,
+		    Standard_Real& Volu,
+                    gp_Pnt& G,
+                    gp_Mat& Inertia) 
+{
+
+  gp_Pnt P;            
+  gp_Vec VNor;     
+  Standard_Real dvi, dv;
+  Standard_Real ur, um, u, vr, vm, v;
+  Standard_Real x, y, z, xn, yn, zn, xi, yi, zi;
+//  Standard_Real x, y, z, xn, yn, zn, xi, yi, zi, xyz;
+  Standard_Real px,py,pz,s,d1,d2,d3;
+  Standard_Real Ixi, Iyi, Izi, Ixxi, Iyyi, Izzi, Ixyi, Ixzi, Iyzi;
+  Standard_Real xloc, yloc, zloc;
+  Standard_Real Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
+
+  Volu = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
+  Loc.Coord (xloc, yloc, zloc);
+
+  Standard_Real LowerU, UpperU, LowerV, UpperV;
+  S.Bounds ( LowerU, UpperU, LowerV, UpperV);
+  Standard_Integer UOrder = Min(S.UIntegrationOrder (),
+				math::GaussPointsMax());
+  Standard_Integer VOrder = Min(S.VIntegrationOrder (),
+				math::GaussPointsMax());
+
+  Standard_Integer i, j;
+  math_Vector GaussPU (1, UOrder);     //gauss points and weights
+  math_Vector GaussWU (1, UOrder);
+  math_Vector GaussPV (1, VOrder);
+  math_Vector GaussWV (1, VOrder);
+
+  math::GaussPoints  (UOrder,GaussPU);
+  math::GaussWeights (UOrder,GaussWU);
+  math::GaussPoints  (VOrder,GaussPV);
+  math::GaussWeights (VOrder,GaussWV);
+
+  um = 0.5 * (UpperU + LowerU);
+  vm = 0.5 * (UpperV + LowerV);
+  ur = 0.5 * (UpperU - LowerU);
+  vr = 0.5 * (UpperV - LowerV);
+
+  for (j = 1; j <= VOrder; j++) {
+    v = vm + vr * GaussPV (j);
+    dvi = Ixi = Iyi = Izi = Ixxi = Iyyi = Izzi = Ixyi = Ixzi = Iyzi = 0.0;
+
+    for (i = 1; i <= UOrder; i++) {
+      u = um + ur * GaussPU (i);
+      S.Normal (u, v, P, VNor); 
+      VNor.Coord (xn, yn, zn);
+      P.Coord (x, y, z);
+      x -= xloc;  y -= yloc;  z -= zloc;
+      xn *= GaussWU (i);  yn *= GaussWU (i);  zn *= GaussWU (i);
+      if (ByPoint) {   
+	/////////////////////           ///////////////////////
+	//    OFV code     //           //    Initial code   //
+	/////////////////////           ///////////////////////
+	// modified by APO 
+	dv = (x*xn+y*yn+z*zn)/3.0;     //xyz  = x * y * z;
+	dvi += dv;                     //Ixyi += zn * xyz; 
+	Ixi += 0.75*x*dv;              //Iyzi += xn * xyz;
+	Iyi += 0.75*y*dv;              //Ixzi += yn * xyz; 
+	Izi += 0.75*z*dv;              //xi = x * x * x * xn / 3.0;
+	x -= Coeff[0];                 //yi = y * y * y * yn / 3.0;
+	y -= Coeff[1];                 //zi = z * z * z * zn / 3.0;
+	z -= Coeff[2];                 //Ixxi += (yi + zi);
+	dv *= 3.0/5.0;                 //Iyyi += (xi + zi);
+	Ixyi -= x*y*dv;                //Izzi += (xi + yi);
+	Iyzi -= y*z*dv;                //x -= Coeff[0];
+	Ixzi -= x*z*dv;                //y -= Coeff[1];
+	xi = x*x;                      //z -= Coeff[2];
+	yi = y*y;                      //dv = x * xn + y * yn + z * zn;
+	zi = z*z;                      //dvi +=  dv; 
+	Ixxi += (yi + zi)*dv;          //Ixi += x * dv;
+	Iyyi += (xi + zi)*dv;          //Iyi += y * dv;
+	Izzi += (xi + yi)*dv;          //Izi += z * dv;
+      }
+      else { // by plane
+	s   = xn * Coeff[0] + yn * Coeff[1] + zn * Coeff[2];
+	d1  = Coeff[0] * x + Coeff[1] * y + Coeff[2] * z - Coeff[3];
+	d2  = d1 * d1;
+	d3  = d1 * d2 / 3.0;
+	dv  = s * d1;
+	dvi += dv;
+	Ixi += (x - (Coeff[0] * d1 / 2.0)) * dv;
+	Iyi += (y - (Coeff[1] * d1 / 2.0)) * dv;
+	Izi += (z - (Coeff[2] * d1 / 2.0)) * dv;
+	px  = x - Coeff[0] * d1;
+	py  = y - Coeff[1] * d1;
+	pz  = z - Coeff[2] * d1;
+	xi  = px * px * d1 + px * Coeff[0]* d2 + Coeff[0] * Coeff[0] * d3;
+	yi  = py * py * d1 + py * Coeff[1] * d2 + Coeff[1] * Coeff[1] * d3;
+	zi  = pz * pz * d1 + pz * Coeff[2] * d2 + Coeff[2] * Coeff[2] * d3;
+	Ixxi += (yi + zi) * s;
+	Iyyi += (xi + zi) * s;
+	Izzi += (xi + yi) * s;
+	d2  /= 2.0;
+	xi  = (py * pz * d1) + (py * Coeff[2] * d2) + (pz * Coeff[1] * d2) + (Coeff[1] * Coeff[2] * d3);
+	yi  = (px * pz * d1) + (pz * Coeff[0] * d2) + (px * Coeff[2] * d2) + (Coeff[0] * Coeff[2] * d3);
+	zi  = (px * py * d1) + (px * Coeff[1] * d2) + (py * Coeff[0] * d2) + (Coeff[0] * Coeff[1] * d3);
+	Ixyi -=  zi * s;
+	Iyzi -=  xi * s;
+	Ixzi -=  yi * s;
+      }
+    }
+    Volu += dvi  * GaussWV (j);
+    Ix   += Ixi  * GaussWV (j);
+    Iy   += Iyi  * GaussWV (j);
+    Iz   += Izi  * GaussWV (j);
+    Ixx  += Ixxi * GaussWV (j);
+    Iyy  += Iyyi * GaussWV (j);
+    Izz  += Izzi * GaussWV (j);
+    Ixy  += Ixyi * GaussWV (j);
+    Ixz  += Ixzi * GaussWV (j);
+    Iyz  += Iyzi * GaussWV (j);
+  }
+  vr  *= ur;
+  Ixx *= vr;
+  Iyy *= vr;
+  Izz *= vr;
+  Ixy *= vr;
+  Ixz *= vr;
+  Iyz *= vr;
+  if (Abs(Volu) >= EPS_DIM ) {
+    if (ByPoint) {
+      Ix  = Coeff[0] + Ix/Volu;
+      Iy  = Coeff[1] + Iy/Volu;
+      Iz  = Coeff[2] + Iz/Volu;
+      Volu *= vr;
+    }
+    else { //by plane
+     Ix /= Volu;
+     Iy /= Volu;
+     Iz /= Volu;
+     Volu *= vr;
+    }
+    G.SetCoord (Ix, Iy, Iz);
+  } 
+  else {
+    G.SetCoord(0.,0.,0.);
+    Volu =0.;
+  }
+  Inertia.SetCols (gp_XYZ (Ixx, Ixy, Ixz),
+		   gp_XYZ (Ixy, Iyy, Iyz),
+		   gp_XYZ (Ixz, Iyz, Izz));
+
+}
+
+// Last modified by OFV 5.2001:
+// 1). surface and edge integration order is equal now
+// 2). "by point" works now rathre correctly (it looks so...)
+static void Compute(Face& S, Domain& D, const Standard_Boolean ByPoint, const Standard_Real  Coeff[],
+                    const gp_Pnt& Loc, Standard_Real& Volu, gp_Pnt& G, gp_Mat& Inertia) 
+
+{
+  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;
+  Standard_Real LocVolu, LocIx, LocIy, LocIz, LocIxx, LocIyy, LocIzz, LocIxy, LocIxz, LocIyz;
+  Standard_Real CVolu, CIx, CIy, CIz, CIxx, CIyy, CIzz, CIxy, CIxz, CIyz, Ix, Iy, Iz, Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
+  Standard_Real xn, yn, zn, px, py, pz, s, d1, d2, d3, dSigma;
+  Standard_Integer i, j, vio, sio, max, NbGaussgp_Pnts;
+  
+  gp_Pnt Ps;
+  gp_Vec VNor; 
+  gp_Pnt2d Puv;
+  gp_Vec2d Vuv;
+      
+  Loc.Coord (xloc, yloc, zloc);
+  Volu = Ix = Iy = Iz = Ixx = Iyy = Izz = Ixy = Ixz = Iyz = 0.0;
+  S.Bounds (u1, u2, v1, v2);
+  Standard_Real _u2 = u2;  //OCC104
+  vio = S.VIntegrationOrder ();
+
+  while (D.More())
+    {
+      S.Load(D.Value());
+      sio = S.IntegrationOrder ();
+      max = Max(vio,sio);
+      NbGaussgp_Pnts = Min(max,math::GaussPointsMax());
+            
+      math_Vector GaussP (1, NbGaussgp_Pnts);
+      math_Vector GaussW (1, NbGaussgp_Pnts);
+      math::GaussPoints  (NbGaussgp_Pnts,GaussP);
+      math::GaussWeights (NbGaussgp_Pnts,GaussW);
+      
+      CVolu = CIx = CIy = CIz = CIxx = CIyy = CIzz = CIxy = CIxz = CIyz = 0.0;
+      l1  = S.FirstParameter();
+      l2  = S.LastParameter();
+      lm  = 0.5 * (l2 + l1);         
+      lr  = 0.5 * (l2 - l1);
+
+      for (i=1; i<=NbGaussgp_Pnts; i++)
+	{
+	  l = lm + lr * GaussP(i);
+	  S.D12d (l, Puv, Vuv);
+	  v   = Puv.Y();
+	  u2  = Puv.X();
+
+	  //OCC104
+	  v = v < v1? v1: v; 
+	  v = v > v2? v2: v; 
+	  u2 = u2 < u1? u1: u2; 
+	  u2 = u2 > _u2? _u2: u2; 
+
+	  Dul = Vuv.Y() * GaussW(i);
+	  um  = 0.5 * (u2 + u1);
+	  ur  = 0.5 * (u2 - u1);
+	  LocVolu = LocIx = LocIy = LocIz = LocIxx = LocIyy = LocIzz = LocIxy = LocIxz = LocIyz = 0.0;
+	  
+	  for (j=1; j<=NbGaussgp_Pnts; j++)
+	    {
+	      u = um + ur * GaussP(j);
+	      S.Normal (u, v, Ps, VNor);
+	      VNor.Coord (xn, yn, zn);
+	      Ps.Coord (x, y, z);
+	      x      -= xloc; 
+	      y      -= yloc;
+	      z      -= zloc;
+	      xn = xn * Dul * GaussW(j);
+	      yn = yn * Dul * GaussW(j);
+	      zn = zn * Dul * GaussW(j);
+	      if(ByPoint)
+		{
+		  dSigma = (x*xn+y*yn+z*zn)/3.0;
+		  LocVolu +=  dSigma; 
+		  LocIx   += 0.75*x*dSigma;  
+		  LocIy   += 0.75*y*dSigma;
+		  LocIz   += 0.75*z*dSigma;
+		  x      -= Coeff[0];
+		  y      -= Coeff[1];
+		  z      -= Coeff[2];
+		  dSigma *= 3.0/5.0; 
+		  LocIxy -= x*y*dSigma;
+		  LocIyz -= y*z*dSigma;
+		  LocIxz -= x*z*dSigma;
+		  xi      = x*x;
+		  yi      = y*y;
+		  zi      = z*z;
+		  LocIxx += (yi + zi)*dSigma;
+		  LocIyy += (xi + zi)*dSigma;
+		  LocIzz += (xi + yi)*dSigma;
+		}
+	      else
+		{
+		  s   = xn * Coeff[0] + yn * Coeff[1] + zn * Coeff[2];
+		  d1  = Coeff[0] * x + Coeff[1] * y + Coeff[2] * z;
+		  d2  = d1 * d1;
+		  d3  = d1 * d2 / 3.0;
+		  ds  = s * d1;
+		  LocVolu += ds;
+		  LocIx += (x - Coeff[0] * d1 / 2.0) * ds;
+		  LocIy += (y - Coeff[1] * d1 / 2.0) * ds;
+		  LocIz += (z - Coeff[2] * d1 / 2.0) * ds;
+		  px  = x - Coeff[0] * d1;
+		  py  = y - Coeff[1] * d1;
+		  pz  = z - Coeff[2] * d1;
+		  xi  = (px * px * d1) + (px * Coeff[0]* d2) + (Coeff[0] * Coeff[0] * d3);
+		  yi  = (py * py * d1) + (py * Coeff[1] * d2) + (Coeff[1] * Coeff[1] * d3);
+		  zi  = pz * pz * d1 + pz * Coeff[2] * d2 + (Coeff[2] * Coeff[2] * d3);
+		  LocIxx += (yi + zi) * s;
+		  LocIyy += (xi + zi) * s;
+		  LocIzz += (xi + yi) * s;
+		  d2  /= 2.0;
+		  xi  = (py * pz * d1) + (py * Coeff[2] * d2) + (pz * Coeff[1] * d2) + (Coeff[1] * Coeff[2] * d3);
+		  yi  = (px * pz * d1) + (pz * Coeff[0] * d2) + (px * Coeff[2] * d2) + (Coeff[0] * Coeff[2] * d3);
+		  zi  = (px * py * d1) + (px * Coeff[1] * d2) + (py * Coeff[0] * d2) + (Coeff[0] * Coeff[1] * d3);
+		  LocIxy -=  zi * s;
+		  LocIyz -=  xi * s;
+		  LocIxz -=  yi * s;
+		}
+	    }
+	  CVolu += LocVolu * ur;
+	  CIx   += LocIx   * ur;
+	  CIy   += LocIy   * ur;
+	  CIz   += LocIz   * ur;
+	  CIxx  += LocIxx  * ur;
+	  CIyy  += LocIyy  * ur;
+	  CIzz  += LocIzz  * ur;
+	  CIxy  += LocIxy  * ur;
+	  CIxz  += LocIxz  * ur;
+	  CIyz  += LocIyz  * ur;
+	}
+      Volu   += CVolu * lr;
+      Ix     += CIx   * lr;
+      Iy     += CIy   * lr;
+      Iz     += CIz   * lr;
+      Ixx    += CIxx  * lr;
+      Iyy    += CIyy  * lr;
+      Izz    += CIzz  * lr;
+      Ixy    += CIxy  * lr;
+      Ixz    += CIxz  * lr;
+      Iyz    += CIyz  * lr;
+      D.Next();
+    }
+
+  if(Abs(Volu) >= EPS_DIM)
+    {
+      if(ByPoint)
+	{
+	  Ix = Coeff[0] + Ix/Volu;
+	  Iy = Coeff[1] + Iy/Volu;
+	  Iz = Coeff[2] + Iz/Volu;
+	}
+      else
+	{
+	  Ix /= Volu;
+	  Iy /= Volu;
+	  Iz /= Volu;
+	}
+      G.SetCoord (Ix, Iy, Iz);
+    }
+  else
+    {
+      Volu =0.;
+      G.SetCoord(0.,0.,0.);
+    }
+
+  Inertia.SetCols (gp_XYZ (Ixx, Ixy, Ixz),
+		   gp_XYZ (Ixy, Iyy, Iyz),
+		   gp_XYZ (Ixz, Iyz, Izz));
+  
+}
+
+GProp_VGProps::GProp_VGProps(){}
+
+GProp_VGProps::GProp_VGProps(Face& S, const gp_Pnt& VLocation, const Standard_Real Eps){
+  SetLocation(VLocation);
+  Perform(S,Eps);
+}
+
+GProp_VGProps::GProp_VGProps(Face& S, Domain& D, const gp_Pnt& VLocation, const Standard_Real Eps){
+  SetLocation(VLocation);
+  Perform(S,D,Eps);
+}
+
+GProp_VGProps::GProp_VGProps(Face& S, Domain& D, const gp_Pnt& VLocation){
+  SetLocation(VLocation);
+  Perform(S,D);
+}
+
+GProp_VGProps::GProp_VGProps(const Face& S, const gp_Pnt& VLocation){
+  SetLocation(VLocation);
+  Perform(S);
+}
+
+GProp_VGProps::GProp_VGProps(Face& S, const gp_Pnt& O, const gp_Pnt& VLocation, const Standard_Real Eps){
+  SetLocation(VLocation);
+  Perform(S,O,Eps);
+}
+
+GProp_VGProps::GProp_VGProps(Face& S, Domain& D, const gp_Pnt& O, const gp_Pnt& VLocation, const Standard_Real Eps){
+  SetLocation(VLocation);
+  Perform(S,D,O,Eps);
+}
+
+GProp_VGProps::GProp_VGProps(const Face& S, const gp_Pnt& O, const gp_Pnt& VLocation){
+  SetLocation(VLocation);
+  Perform(S,O);
+}
+
+GProp_VGProps::GProp_VGProps(Face& S, Domain& D, const gp_Pnt& O, const gp_Pnt& VLocation){
+  SetLocation(VLocation);
+  Perform(S,D,O);
+}
+
+GProp_VGProps::GProp_VGProps(Face& S, const gp_Pln& Pl, const gp_Pnt& VLocation, const Standard_Real Eps){
+  SetLocation(VLocation);
+  Perform(S,Pl,Eps);
+}
+
+GProp_VGProps::GProp_VGProps(Face& S, Domain& D, const gp_Pln& Pl, const gp_Pnt& VLocation, const Standard_Real Eps){
+  SetLocation(VLocation);
+  Perform(S,D,Pl,Eps);
+}
+
+GProp_VGProps::GProp_VGProps(const Face& S, const gp_Pln& Pl, const gp_Pnt& VLocation){
+  SetLocation(VLocation);
+  Perform(S,Pl);
+}
+
+GProp_VGProps::GProp_VGProps(Face& S, Domain& D, const gp_Pln& Pl, const gp_Pnt& VLocation){
+  SetLocation(VLocation);
+  Perform(S,D,Pl);
+}
+
+void GProp_VGProps::SetLocation(const gp_Pnt& VLocation){
+  loc = VLocation;
+}
+
+Standard_Real GProp_VGProps::Perform(Face& S, const Standard_Real Eps){
+  Standard_Real Coeff[] = {0., 0., 0.};
+  return myEpsilon = Compute(S,Standard_True,Coeff,loc,dim,g,inertia,Eps);
+}
+
+Standard_Real GProp_VGProps::Perform(Face& S, Domain& D, const Standard_Real Eps){
+  Standard_Real Coeff[] = {0., 0., 0.};
+  return myEpsilon = Compute(S,D,Standard_True,Coeff,loc,dim,g,inertia,Eps);
+}
+
+void GProp_VGProps::Perform(const Face& S){
+  Standard_Real Coeff[] = {0., 0., 0.};
+  Compute(S,Standard_True,Coeff,loc,dim,g,inertia);
+  myEpsilon = 1.0;
+  return;
+}
+
+void GProp_VGProps::Perform(Face& S, Domain& D){
+  Standard_Real Coeff[] = {0., 0., 0.};
+  Compute(S,D,Standard_True,Coeff,loc,dim,g,inertia);
+  myEpsilon = 1.0;
+  return;
+}
+
+Standard_Real GProp_VGProps::Perform(Face& S, const gp_Pnt&  O, const Standard_Real Eps){
+  Standard_Real xloc, yloc, zloc;
+  loc.Coord(xloc, yloc, zloc);
+  Standard_Real Coeff[3];
+  O.Coord (Coeff[0], Coeff[1], Coeff[2]);
+  Coeff[0] -= xloc;  Coeff[1] -= yloc;   Coeff[2] -= zloc;
+  return myEpsilon = Compute(S,Standard_True,Coeff,loc,dim,g,inertia,Eps);
+}
+
+Standard_Real GProp_VGProps::Perform(Face& S, Domain& D, const gp_Pnt&  O, const Standard_Real Eps){
+  Standard_Real xloc, yloc, zloc;
+  loc.Coord(xloc, yloc, zloc);
+  Standard_Real Coeff[3];
+  O.Coord (Coeff[0], Coeff[1], Coeff[2]);
+  Coeff[0] -= xloc;  Coeff[1] -= yloc;   Coeff[2] -= zloc;
+  return myEpsilon = Compute(S,D,Standard_True,Coeff,loc,dim,g,inertia,Eps);
+}
+
+void GProp_VGProps::Perform(const Face& S, const gp_Pnt&  O){
+  Standard_Real xloc, yloc, zloc;
+  loc.Coord(xloc, yloc, zloc);
+  Standard_Real Coeff[3];
+  O.Coord (Coeff[0], Coeff[1], Coeff[2]);
+  Coeff[0] -= xloc;  Coeff[1] -= yloc;   Coeff[2] -= zloc;
+  Compute(S,Standard_True,Coeff,loc,dim,g,inertia);
+  myEpsilon = 1.0;
+  return;
+}
+
+void GProp_VGProps::Perform(Face& S, Domain& D, const gp_Pnt&  O){
+  Standard_Real xloc, yloc, zloc;
+  loc.Coord(xloc, yloc, zloc);
+  Standard_Real Coeff[3];
+  O.Coord (Coeff[0], Coeff[1], Coeff[2]);
+  Coeff[0] -= xloc;  Coeff[1] -= yloc;   Coeff[2] -= zloc;
+  Compute(S,D,Standard_True,Coeff,loc,dim,g,inertia);
+  myEpsilon = 1.0;
+  return;
+}
+
+Standard_Real GProp_VGProps::Perform(Face& S, const gp_Pln& Pl, const Standard_Real Eps){
+  Standard_Real xloc, yloc, zloc;
+  loc.Coord (xloc, yloc, zloc);
+  Standard_Real Coeff[4];
+  Pl.Coefficients (Coeff[0], Coeff[1],Coeff[2],Coeff[3]); 
+  Coeff[3] = Coeff[3] - Coeff[0]*xloc - Coeff[1]*yloc - Coeff[2]*zloc;
+  return myEpsilon = Compute(S,Standard_False,Coeff,loc,dim,g,inertia,Eps);
+}
+
+Standard_Real GProp_VGProps::Perform(Face& S, Domain& D, const gp_Pln& Pl, const Standard_Real Eps){
+  Standard_Real xloc, yloc, zloc;
+  loc.Coord (xloc, yloc, zloc);
+  Standard_Real Coeff[4];
+  Pl.Coefficients (Coeff[0], Coeff[1],Coeff[2],Coeff[3]); 
+  Coeff[3] = Coeff[3] - Coeff[0]*xloc - Coeff[1]*yloc - Coeff[2]*zloc;
+  return myEpsilon = Compute(S,D,Standard_False,Coeff,loc,dim,g,inertia,Eps);
+}
+
+void GProp_VGProps::Perform(const Face& S, const gp_Pln& Pl){
+  Standard_Real xloc, yloc, zloc;
+  loc.Coord (xloc, yloc, zloc);
+  Standard_Real Coeff[4];
+  Pl.Coefficients (Coeff[0], Coeff[1],Coeff[2],Coeff[3]); 
+  Coeff[3] = Coeff[3] - Coeff[0]*xloc - Coeff[1]*yloc - Coeff[2]*zloc;
+  Compute(S,Standard_False,Coeff,loc,dim,g,inertia);
+  myEpsilon = 1.0;
+  return;
+}
+
+void GProp_VGProps::Perform(Face& S, Domain& D, const gp_Pln& Pl){
+  Standard_Real xloc, yloc, zloc;
+  loc.Coord (xloc, yloc, zloc);
+  Standard_Real Coeff[4];
+  Pl.Coefficients (Coeff[0], Coeff[1],Coeff[2],Coeff[3]); 
+  Coeff[3] = Coeff[3] - Coeff[0]*xloc - Coeff[1]*yloc - Coeff[2]*zloc;
+  Compute(S,D,Standard_False,Coeff,loc,dim,g,inertia);
+  myEpsilon = 1.0;
+  return;
+}
+
+Standard_Real GProp_VGProps::GetEpsilon(){
+  return myEpsilon;
+}

src/GProp/GProp_VGPropsGK.cdl

@@ -0,0 +1,464 @@
+-- File:	GProp_VGPropsGK.cdl
+-- Created:	Wed Dec 21 17:35:57 2005
+-- Author:	Sergey KHROMOV
+--		<skv@dimox>
+
+
+generic class VGPropsGK from GProp (Arc as any;
+                                  Face as any;  -- as FaceTool(Arc)
+    	    	    	    	  Domain as any -- as DomainTool(Arc)
+    	    	    	    	 )
+inherits GProps from GProp
+
+        ---Purpose: Computes the global properties of a geometric solid 
+        --          (3D closed region of space) delimited with :
+        --            -  a point and a surface
+        --            -  a plane and a surface
+        --  
+        --          The surface can be :
+        --            -  a surface limited with its parametric values U-V, 
+        --               (naturally restricted)
+        --            -  a surface limited in U-V space with its boundary 
+        --               curves.
+        --  
+        --          The surface's requirements to evaluate the global 
+        --          properties are defined in the template FaceTool class from 
+        --          the package GProp. 
+	--  
+	--          The adaptive 2D algorithm of Gauss-Kronrod integration of 
+        --          double integral is used. 
+        --  
+        --          The inner integral is computed along U parameter of 
+        --          surface. The integrand function is encapsulated in the 
+        --          support class UFunction that is defined below. 
+	--  
+	--          The outer integral is computed along T parameter of a 
+        --          bounding curve. The integrand function is encapsulated in 
+        --          the support class TFunction that is defined below.
+
+uses 
+ 
+    Pnt     from gp, 
+    XYZ     from gp, 
+    Pln     from gp, 
+    Address from Standard, 
+    Boolean from Standard, 
+    Real    from Standard 
+
+
+--  Template class functions. Used for integration. Begin
+
+    class UFunction from GProp inherits Function from math 
+        ---Purpose: This class represents the integrand function for 
+        --          computation of an inner integral. The returned value 
+        --          depends on the value type and the flag IsByPoint. 
+	--  
+	--          The type of returned value is the one of the following 
+        --          values: 
+	--            -  GProp_Mass - volume computation. 
+	--            -  GProp_CenterMassX, GProp_CenterMassY, 
+        --               GProp_CenterMassZ - X, Y and Z coordinates of center 
+        --                                   of mass computation. 
+	--            -  GProp_InertiaXX, GProp_InertiaYY, GProp_InertiaZZ, 
+        --               GProp_InertiaXY, GProp_InertiaXZ, GProp_InertiaYZ 
+        --                                 - moments of inertia computation. 
+	--  
+	--          If the flag IsByPoint is set to Standard_True, the value is 
+        --          returned for the region of space that is delimited by a 
+        --          surface and a point. Otherwise all computations are 
+        --          performed for the region of space delimited by a surface 
+        --          and a plane.
+     
+    uses 
+     
+        Pnt       from gp, 
+	XYZ       from gp, 
+        Address   from Standard, 
+        Boolean   from Standard, 
+	Real      from Standard, 
+	ValueType from GProp
+    
+    is 
+     
+        Create(theSurface: Face; 
+               theVertex : Pnt     from gp; 
+	       IsByPoint : Boolean from Standard; 
+               theCoeffs : Address from Standard) 
+	---Purpose: Constructor. Initializes the function with the face, the 
+        --          location point, the flag IsByPoint and the coefficients 
+        --          theCoeff that have different meaning depending on the value 
+        --          of IsByPoint. 
+        --          If IsByPoint is equal to Standard_True, the number of the 
+        --          coefficients is equal to 3 and they represent X, Y and Z 
+        --          coordinates (theCoeff[0], theCoeff[1] and theCoeff[2] 
+        --          correspondingly) of the shift, if the inertia is computed 
+        --          with respect to the point different then the location. 
+        --          If IsByPoint is equal to Standard_False, the number of the 
+        --          coefficients is 4 and they represent the combination of 
+        --          plane parameters and shift values.
+        returns UFunction from GProp;  
+     
+        SetValueType(me: in out; theType: ValueType from GProp); 
+	---Purpose: Setting the type of the value to be returned. 
+	---C++: inline 
+     
+	SetVParam(me: in out; theVParam: Real from Standard);
+	---Purpose: Setting the V parameter that is constant during the 
+        --          integral computation.
+	---C++: inline 
+     
+        Value(me: in out; X:     Real from Standard; 
+                          F: out Real from Standard) 
+        ---Purpose: Returns a value of the function.
+	returns Boolean from Standard 
+	is redefined; 
+     
+        -----------------------
+        --  Private methods  --
+        -----------------------
+
+        VolumeValue(me: in out; X      :     Real from Standard; 
+                                thePMP0: out XYZ  from gp; 
+                                theS   : out Real from Standard; 
+                                theD1  : out Real from Standard) 
+        ---Purpose: Private method. Returns the value for volume computation. 
+        --          Other returned values are: 
+	--            -  thePMP0 - PSurf(X,Y) minus Location. 
+	--            -  theS and theD1 coeffitients that are computed and used 
+        --               for computation of center of mass and inertia values 
+        --               by plane.
+	returns Real from Standard 
+        is private; 
+     
+        CenterMassValue(me: in out; X:     Real from Standard; 
+                                    F: out Real from Standard) 
+        ---Purpose: Private method. Returns a value for the center of mass 
+        --          computation. If the value type other then GProp_CenterMassX,
+        --          GProp_CenterMassY or GProp_CenterMassZ this method returns 
+        --          Standard_False. Returns Standard_True in case of successful 
+        --          computation of a value.
+	returns Boolean from Standard 
+        is private; 
+     
+        InertiaValue(me: in out; X:     Real from Standard; 
+                                 F: out Real from Standard) 
+        ---Purpose: Private method. Computes the value of intertia. The type of 
+        --          a value returned is defined by the value type. If it is 
+        --          other then GProp_InertiaXX, GProp_InertiaYY, 
+        --          GProp_InertiaZZ, GProp_InertiaXY, GProp_InertiaXZ or 
+        --          GProp_InertiaYZ, the method returns Standard_False. Returns 
+        --          Standard_True in case of successful computation of a value.
+	returns Boolean from Standard 
+        is private; 
+     
+    fields 
+     
+        mySurface  : Face; 
+	myVertex   : Pnt       from gp; 
+	myCoeffs   : Address   from Standard; 
+	myVParam   : Real      from Standard; 
+	myValueType: ValueType from GProp;
+	myIsByPoint: Boolean   from Standard; 
+     
+    end UFunction;
+
+
+    --  Class  TFunction.
+
+    class TFunction from GProp inherits Function from math
+        ---Purpose: This class represents the integrand function for the outer 
+        --          integral computation. The returned value represents the 
+        --          integral of UFunction. It depends on the value type and the 
+        --          flag IsByPoint. 
+     
+    uses 
+     
+        Pnt       from gp, 
+        Address   from Standard, 
+        Boolean   from Standard, 
+        Integer   from Standard, 
+	Real      from Standard, 
+	ValueType from GProp
+    
+    is 
+     
+        Create(theSurface  : Face; 
+               theVertex   : Pnt     from gp; 
+	       IsByPoint   : Boolean from Standard; 
+               theCoeffs   : Address from Standard; 
+               theUMin     : Real    from Standard; 
+               theTolerance: Real    from Standard) 
+	---Purpose: Constructor. Initializes the function with the face, the 
+        --          location point, the flag IsByPoint, the coefficients 
+        --          theCoeff that have different meaning depending on the value 
+        --          of IsByPoint. The last two parameters are theUMin - the 
+        --          lower bound of the inner integral. This value is fixed for 
+        --          any integral. And the value of tolerance of inner integral 
+        --          computation.
+        --          If IsByPoint is equal to Standard_True, the number of the 
+        --          coefficients is equal to 3 and they represent X, Y and Z 
+        --          coordinates (theCoeff[0], theCoeff[1] and theCoeff[2] 
+        --          correspondingly) of the shift if the inertia is computed 
+        --          with respect to the point different then the location. 
+        --          If IsByPoint is equal to Standard_False, the number of the 
+        --          coefficients is 4 and they represent the compbination of 
+        --          plane parameters and shift values.
+        returns TFunction from GProp;   
+	 
+	Init(me: in out);
+      
+        SetNbKronrodPoints(me: in out; theNbPoints: Integer from Standard);
+        ---Purpose: Setting the expected number of Kronrod points for the outer 
+        --          integral computation. This number is required for 
+        --          computation of a value of tolerance for inner integral 
+        --          computation. After GetStateNumber method call, this number 
+        --          is recomputed by the same law as in 
+        --          math_KronrodSingleIntegration, i.e. next number of points 
+        --          is equal to the current number plus a square root of the 
+        --          current number. If the law in math_KronrodSingleIntegration
+        --          is changed, the modification algo should be modified 
+        --          accordingly.
+	---C++: inline 
+     
+        SetValueType(me: in out; aType: ValueType from GProp);
+	---Purpose: Setting the type of the value to be returned. This 
+        --          parameter is directly passed to the UFunction. 
+	---C++: inline 
+
+        SetTolerance(me: in out; aTol: Real from Standard);
+	---Purpose: Setting the tolerance  for  inner integration
+	---C++: inline 
+     
+        ErrorReached(me) 
+	---Purpose: Returns the relative reached error of all values computation since 
+        --          the last call of GetStateNumber method.
+	---C++: inline 
+	returns Real from Standard;
+
+        AbsolutError(me) 
+	---Purpose: Returns the absolut reached error of all values computation since 
+        --          the last call of GetStateNumber method.
+	---C++: inline 
+	returns Real from Standard;
+     
+        Value(me: in out; X:     Real from Standard; 
+                          F: out Real from Standard) 
+        ---Purpose: Returns a value of the function. The value represents an 
+        --          integral of UFunction. It is computed with the predefined 
+        --          tolerance using the adaptive Gauss-Kronrod method.
+	returns Boolean from Standard 
+	is redefined; 
+     
+        GetStateNumber(me: in out) 
+	---Purpose:  Redefined  method. Remembers the error reached during 
+        --           computation of integral values since the object creation 
+        --           or the last call of GetStateNumber. It is invoked in each 
+        --           algorithm from the package math. Particularly in the 
+        --           algorithm math_KronrodSingleIntegration that is used to 
+        --           compute the integral of TFunction.
+        returns Integer
+        is redefined;
+     
+    fields 
+     
+        mySurface   : Face; 
+	myUFunction : UFunction; 
+	myUMin      : Real      from Standard; 
+	myTolerance : Real      from Standard; 
+	myTolReached: Real      from Standard; 
+	myErrReached: Real      from Standard; 
+	myAbsError  : Real      from Standard; 
+	myValueType : ValueType from GProp; 
+	myIsByPoint : Boolean   from Standard; 
+	myNbPntOuter: Integer   from Standard;
+    
+    end TFunction;
+
+--  Template class functions. Used for integration. End
+
+is
+ 
+    Create
+        ---Purpose: Empty constructor. 
+	---C++: inline 
+    returns VGPropsGK; 
+
+    Create(theSurface  : in out Face; 
+           theLocation :        Pnt  from gp; 
+           theTolerance:        Real from Standard = 0.001; 
+           theCGFlag: Boolean from Standard = Standard_False;  
+           theIFlag: Boolean from Standard = Standard_False)
+        ---Purpose: Constructor. Computes the global properties of a region of 
+        --          3D space delimited with the naturally restricted surface 
+        --          and the point VLocation.
+    returns VGPropsGK; 
+
+    Create(theSurface  : in out Face; 
+           thePoint    :        Pnt  from gp;
+           theLocation :        Pnt  from gp; 
+           theTolerance:        Real from Standard = 0.001; 
+           theCGFlag: Boolean from Standard = Standard_False;  
+           theIFlag: Boolean from Standard = Standard_False)
+
+        ---Purpose: Constructor. Computes the global properties of a region of 
+        --          3D space delimited with the naturally restricted surface 
+        --          and the point VLocation. The inertia is computed with 
+        --          respect to thePoint.
+    returns VGPropsGK; 
+
+    Create(theSurface  : in out Face; 
+           theDomain   : in out Domain; 
+           theLocation :        Pnt  from gp; 
+           theTolerance:        Real from Standard = 0.001; 
+           theCGFlag: Boolean from Standard = Standard_False;  
+           theIFlag: Boolean from Standard = Standard_False)
+	   
+        ---Purpose: Constructor. Computes the global properties of a region of 
+        --          3D space delimited with the surface bounded by the domain 
+        --          and the point VLocation.
+    returns VGPropsGK; 
+
+    Create(theSurface  : in out Face; 
+           theDomain   : in out Domain; 
+           thePoint    :        Pnt  from gp; 
+           theLocation :        Pnt  from gp; 
+           theTolerance:        Real from Standard = 0.001;
+           theCGFlag: Boolean from Standard = Standard_False;  
+           theIFlag: Boolean from Standard = Standard_False)
+        ---Purpose: Constructor. Computes the global properties of a region of 
+        --          3D space delimited with the surface bounded by the domain 
+        --          and the point VLocation. The inertia is computed with 
+        --          respect to thePoint.
+    returns VGPropsGK; 
+
+    Create(theSurface  : in out Face; 
+           thePlane    :        Pln  from gp; 
+           theLocation :        Pnt  from gp; 
+           theTolerance:        Real from Standard = 0.001; 
+           theCGFlag: Boolean from Standard = Standard_False;  
+           theIFlag: Boolean from Standard = Standard_False)
+	   
+        ---Purpose: Constructor. Computes the global properties of a region of 
+        --          3D space delimited with the naturally restricted surface 
+        --          and the plane.
+    returns VGPropsGK; 
+
+    Create(theSurface  : in out Face; 
+           theDomain   : in out Domain; 
+           thePlane    :        Pln  from gp; 
+           theLocation :        Pnt  from gp; 
+           theTolerance:        Real from Standard = 0.001; 
+           theCGFlag: Boolean from Standard = Standard_False;  
+           theIFlag: Boolean from Standard = Standard_False)
+	   
+        ---Purpose: Constructor. Computes the global properties of a region of 
+        --          3D space delimited with the surface bounded by the domain 
+        --          and the plane.
+    returns VGPropsGK; 
+
+    SetLocation(me: in out; theLocation: Pnt from gp);
+        ---Purpose:  Sets the vertex that delimit 3D closed region of space.
+	---C++: inline 
+
+    Perform(me: in out; theSurface  : in out Face; 
+                        theTolerance:        Real from Standard = 0.001; 
+                        theCGFlag: Boolean from Standard = Standard_False;  
+                        theIFlag: Boolean from Standard = Standard_False)
+			
+        ---Purpose: Computes the global properties of a region of 3D space 
+        --          delimited with the naturally restricted surface and the 
+        --          point VLocation.
+    returns Real from Standard; 
+
+    Perform(me: in out; theSurface  : in out Face; 
+			thePoint    :        Pnt  from gp;
+                        theTolerance:        Real from Standard = 0.001; 
+                        theCGFlag: Boolean from Standard = Standard_False;  
+                        theIFlag: Boolean from Standard = Standard_False)
+			
+        ---Purpose: Computes the global properties of a region of 3D space 
+        --          delimited with the naturally restricted surface and the 
+        --          point VLocation. The inertia is computed with respect to 
+        --          thePoint.
+    returns Real from Standard; 
+
+    Perform(me: in out; theSurface  : in out Face; 
+                        theDomain   : in out Domain; 
+                        theTolerance:        Real from Standard = 0.001; 
+                        theCGFlag: Boolean from Standard = Standard_False;  
+                        theIFlag: Boolean from Standard = Standard_False)
+			
+        ---Purpose: Computes the global properties of a region of 3D space 
+        --          delimited with the surface bounded by the domain and the 
+        --          point VLocation.
+    returns Real from Standard; 
+
+    Perform(me: in out; theSurface  : in out Face; 
+                        theDomain   : in out Domain; 
+			thePoint    :        Pnt  from gp;
+                        theTolerance:        Real from Standard = 0.001; 
+                        theCGFlag: Boolean from Standard = Standard_False;  
+                        theIFlag: Boolean from Standard = Standard_False)
+        ---Purpose: Computes the global properties of a region of 3D space 
+        --          delimited with the surface bounded by the domain and the 
+        --          point VLocation. The inertia is computed with respect to 
+        --          thePoint.
+    returns Real from Standard; 
+
+    Perform(me: in out; theSurface  : in out Face; 
+                        thePlane    :        Pln  from gp;
+                        theTolerance:        Real from Standard = 0.001; 
+                        theCGFlag: Boolean from Standard = Standard_False;  
+                        theIFlag: Boolean from Standard = Standard_False)
+			
+        ---Purpose: Computes the global properties of a region of 3D space 
+        --          delimited with the naturally restricted surface and the 
+        --          plane.
+    returns Real from Standard; 
+
+    Perform(me: in out; theSurface  : in out Face; 
+                        theDomain   : in out Domain; 
+                        thePlane    :        Pln  from gp;
+                        theTolerance:        Real from Standard = 0.001; 
+                        theCGFlag: Boolean from Standard = Standard_False;  
+                        theIFlag: Boolean from Standard = Standard_False)
+			
+        ---Purpose: Computes the global properties of a region of 3D space 
+        --          delimited with the surface bounded by the domain and the 
+        --          plane.
+    returns Real from Standard; 
+
+    GetErrorReached(me)
+        ---Purpose: Returns the relative reached computation error. 
+	---C++: inline 
+    returns Real from Standard;
+
+    GetAbsolutError(me)
+        ---Purpose: Returns the absolut reached computation error. 
+	---C++: inline 
+    returns Real from Standard;
+  
+-----------------------
+--  Private methods  --
+-----------------------
+
+    PrivatePerform(me: in out; 
+           theSurface  : in out Face; 
+           thePtrDomain:        Address from Standard; --  pointer to Domain.
+           IsByPoint   :        Boolean from Standard; 
+	   theCoeffs   :        Address from Standard; 
+           theTolerance:        Real    from Standard; 
+           theCGFlag   :        Boolean from Standard;  
+           theIFlag    :        Boolean from Standard)
+	   
+        ---Purpose: Main method for computation of the global properties that 
+        --          is invoked by each Perform method.
+    returns Real from Standard 
+    is private; 
+
+fields
+ 
+    myErrorReached: Real from Standard;
+    myAbsolutError: Real from Standard;
+
+end VGPropsGK;

src/GProp/GProp_VGPropsGK.gxx

@@ -0,0 +1,475 @@
+// File:	GProp_VGPropsGK.gxx
+// Created:	Fri Dec  9 14:00:04 2005
+// Author:	Sergey KHROMOV
+//		<skv@dimox>
+
+
+#include <TColStd_HArray1OfReal.hxx>
+#include <TColStd_Array1OfBoolean.hxx>
+#include <math_KronrodSingleIntegration.hxx>
+#include <math_Vector.hxx>
+#include <math.hxx>
+
+//==========================================================================
+//function : Constructor
+//           
+//==========================================================================
+
+GProp_VGPropsGK::GProp_VGPropsGK(      Face          &theSurface,
+				 const gp_Pnt        &theLocation,
+				 const Standard_Real  theTolerance,
+				 const Standard_Boolean theCGFlag,
+				 const Standard_Boolean theIFlag)
+     : myErrorReached(0.)
+{
+  SetLocation(theLocation);
+  Perform(theSurface, theTolerance, theCGFlag, theIFlag);
+}
+
+//==========================================================================
+//function : Constructor
+//           
+//==========================================================================
+
+GProp_VGPropsGK::GProp_VGPropsGK(      Face          &theSurface,
+				 const gp_Pnt        &thePoint,
+				 const gp_Pnt        &theLocation,
+				 const Standard_Real  theTolerance,
+				 const Standard_Boolean theCGFlag,
+				 const Standard_Boolean theIFlag)
+
+     : myErrorReached(0.)
+{
+  SetLocation(theLocation);
+  Perform(theSurface, thePoint, theTolerance, theCGFlag, theIFlag);
+}
+
+//==========================================================================
+//function : Constructor
+//           
+//==========================================================================
+
+GProp_VGPropsGK::GProp_VGPropsGK(      Face          &theSurface,
+				 Domain        &theDomain,
+				 const gp_Pnt        &theLocation,
+				 const Standard_Real  theTolerance,
+				 const Standard_Boolean theCGFlag,
+				 const Standard_Boolean theIFlag)
+
+     : myErrorReached(0.)
+{
+  SetLocation(theLocation);
+  Perform(theSurface, theDomain, theTolerance, theCGFlag, theIFlag);
+}
+
+//==========================================================================
+//function : Constructor
+//           
+//==========================================================================
+
+GProp_VGPropsGK::GProp_VGPropsGK(      Face          &theSurface,
+			           Domain        &theDomain,
+				 const gp_Pnt        &thePoint,
+				 const gp_Pnt        &theLocation,
+				 const Standard_Real  theTolerance,
+				 const Standard_Boolean theCGFlag,
+				 const Standard_Boolean theIFlag)
+
+     : myErrorReached(0.)
+{
+  SetLocation(theLocation);
+  Perform(theSurface, theDomain, thePoint, theTolerance, theCGFlag, theIFlag);
+}
+
+//==========================================================================
+//function : Constructor
+//           
+//==========================================================================
+
+GProp_VGPropsGK::GProp_VGPropsGK(      Face          &theSurface,
+				 const gp_Pln        &thePlane,
+				 const gp_Pnt        &theLocation,
+				 const Standard_Real  theTolerance,
+				 const Standard_Boolean theCGFlag,
+				 const Standard_Boolean theIFlag)
+
+     : myErrorReached(0.)
+{
+  SetLocation(theLocation);
+  Perform(theSurface, thePlane, theTolerance, theCGFlag, theIFlag);
+}
+
+//==========================================================================
+//function : Constructor
+//           
+//==========================================================================
+
+GProp_VGPropsGK::GProp_VGPropsGK(      Face          &theSurface,
+			           Domain        &theDomain,
+				 const gp_Pln        &thePlane,
+				 const gp_Pnt        &theLocation,
+				 const Standard_Real  theTolerance,
+				 const Standard_Boolean theCGFlag,
+				 const Standard_Boolean theIFlag)
+
+     : myErrorReached(0.)
+{
+  SetLocation(theLocation);
+  Perform(theSurface, theDomain, thePlane, theTolerance, theCGFlag, theIFlag);
+}
+
+//==========================================================================
+//function : Perform
+//           Compute the properties.
+//==========================================================================
+
+Standard_Real GProp_VGPropsGK::Perform(      Face          &theSurface,
+				       const Standard_Real  theTolerance,
+				       const Standard_Boolean theCGFlag,
+				       const Standard_Boolean theIFlag)
+
+{
+  Standard_Real aShift[] = { 0., 0., 0. };
+
+  return PrivatePerform(theSurface, NULL, Standard_True, &aShift, theTolerance, 
+			theCGFlag, theIFlag);
+}
+
+//==========================================================================
+//function : Perform
+//           Compute the properties.
+//==========================================================================
+
+Standard_Real GProp_VGPropsGK::Perform(      Face          &theSurface,
+				       const gp_Pnt        &thePoint,
+				       const Standard_Real  theTolerance,
+				       const Standard_Boolean theCGFlag,
+				       const Standard_Boolean theIFlag)
+
+{
+  gp_XYZ        aXYZ(thePoint.XYZ().Subtracted(loc.XYZ()));
+  Standard_Real aShift[3];
+
+  aXYZ.Coord(aShift[0], aShift[1], aShift[2]);
+
+  return PrivatePerform(theSurface, NULL, Standard_True, &aShift, theTolerance, 
+			theCGFlag, theIFlag);
+}
+
+//==========================================================================
+//function : Perform
+//           Compute the properties.
+//==========================================================================
+
+Standard_Real GProp_VGPropsGK::Perform(      Face          &theSurface,
+				       Domain        &theDomain,
+				       const Standard_Real  theTolerance,
+				       const Standard_Boolean theCGFlag,
+				       const Standard_Boolean theIFlag)
+
+{
+  Standard_Real aShift[] = { 0., 0., 0. };
+
+  return PrivatePerform(theSurface, &theDomain,
+			Standard_True, &aShift, theTolerance, 
+			theCGFlag, theIFlag);
+}
+
+//==========================================================================
+//function : Perform
+//           Compute the properties.
+//==========================================================================
+
+Standard_Real GProp_VGPropsGK::Perform(      Face          &theSurface,
+				       Domain        &theDomain,
+				       const gp_Pnt        &thePoint,
+				       const Standard_Real  theTolerance,
+				       const Standard_Boolean theCGFlag,
+				       const Standard_Boolean theIFlag)
+
+{
+  gp_XYZ        aXYZ(thePoint.XYZ().Subtracted(loc.XYZ()));
+  Standard_Real aShift[3];
+
+  aXYZ.Coord(aShift[0], aShift[1], aShift[2]);
+
+  return PrivatePerform(theSurface, &theDomain,
+			Standard_True, &aShift, theTolerance, 
+			theCGFlag, theIFlag);
+}
+
+//==========================================================================
+//function : Perform
+//           Compute the properties.
+//==========================================================================
+
+Standard_Real GProp_VGPropsGK::Perform(      Face          &theSurface,
+				       const gp_Pln        &thePlane,
+				       const Standard_Real  theTolerance,
+				       const Standard_Boolean theCGFlag,
+				       const Standard_Boolean theIFlag)
+
+{
+  Standard_Real aCoeff[4];
+  Standard_Real aXLoc;
+  Standard_Real aYLoc;
+  Standard_Real aZLoc;
+
+  loc.Coord(aXLoc, aYLoc, aZLoc);
+  thePlane.Coefficients (aCoeff[0], aCoeff[1], aCoeff[2], aCoeff[3]);
+  aCoeff[3] = aCoeff[3] - aCoeff[0]*aXLoc - aCoeff[1]*aYLoc - aCoeff[2]*aZLoc;
+
+  return PrivatePerform(theSurface, NULL,
+			Standard_False, &aCoeff, theTolerance, 
+			theCGFlag, theIFlag);
+}
+
+//==========================================================================
+//function : Perform
+//           Compute the properties.
+//==========================================================================
+
+Standard_Real GProp_VGPropsGK::Perform(      Face          &theSurface,
+				       Domain        &theDomain,
+				       const gp_Pln        &thePlane,
+				       const Standard_Real  theTolerance,
+				       const Standard_Boolean theCGFlag,
+				       const Standard_Boolean theIFlag)
+
+{
+  Standard_Real aCoeff[4];
+  Standard_Real aXLoc;
+  Standard_Real aYLoc;
+  Standard_Real aZLoc;
+
+  loc.Coord(aXLoc, aYLoc, aZLoc);
+  thePlane.Coefficients (aCoeff[0], aCoeff[1], aCoeff[2], aCoeff[3]);
+  aCoeff[3] = aCoeff[3] - aCoeff[0]*aXLoc - aCoeff[1]*aYLoc - aCoeff[2]*aZLoc;
+
+  return PrivatePerform(theSurface, &theDomain,
+			Standard_False, &aCoeff, theTolerance, 
+			theCGFlag, theIFlag);
+}
+
+//==========================================================================
+//function : PrivatePerform
+//           Compute the properties.
+//==========================================================================
+
+Standard_Real GProp_VGPropsGK::PrivatePerform
+              (      Face             &theSurface,
+	       const Standard_Address  thePtrDomain,
+	       const Standard_Boolean  IsByPoint,
+	       const Standard_Address  theCoeffs,
+	       const Standard_Real     theTolerance,
+	       const Standard_Boolean theCGFlag,
+	       const Standard_Boolean theIFlag)
+
+{
+
+  const Standard_Real aTTol = 1.e-9;
+  Standard_Real *aCoeffs = (Standard_Real *)theCoeffs;
+
+  // Compute the number of 2d bounding curves of the face.
+  Domain           *aPDomain = NULL;
+  Standard_Integer  aNbCurves = 0;
+
+  // If the pointer to the domain is NULL, there is only one curve to treat:
+  // U isoline with the UMax parameter.
+  if (thePtrDomain == NULL)
+    aNbCurves = 1;
+  else {
+    aPDomain = (Domain *)thePtrDomain;
+
+    for (aPDomain->Init(); aPDomain->More(); aPDomain->Next())
+      aNbCurves++;
+  }
+
+  if (aNbCurves == 0) {
+    myErrorReached = -1.;
+
+    return myErrorReached;
+  }
+
+  //Standard_Real    aCrvTol = 0.5*theTolerance/aNbCurves;
+  Standard_Real    aCrvTol = 0.1*theTolerance;
+  Standard_Real    aUMin;
+  Standard_Real    aUMax;
+  Standard_Real    aTMin;
+  Standard_Real    aTMax;
+  Standard_Integer aNbPnts;
+  Standard_Integer aNbMaxIter = 1000;
+  Standard_Integer aNbVal = 10;
+  Standard_Integer k;
+  math_Vector      aLocalValue(1, aNbVal);
+  math_Vector      aLocalTolReached(1, aNbVal);
+  math_Vector      aValue(1, aNbVal);
+  math_Vector      aTolReached(1, aNbVal);
+  TColStd_Array1OfBoolean CFlags(1, aNbVal); 
+  CFlags.Init(Standard_False);
+  Standard_Boolean isMore;
+
+  //aNbVal = 1;
+  aValue.Init(0.);
+  aTolReached.Init(0.);
+
+  CFlags.Init(Standard_False);
+  CFlags(1) = Standard_True;
+
+  if(theCGFlag || theIFlag) {
+    Standard_Integer i;
+    for(i = 2; i <= 4; ++i) {CFlags(i) = Standard_True;}
+  }
+
+  if(theIFlag) {
+    Standard_Integer i;
+    for(i = 5; i <= 10; ++i) {CFlags(i) = Standard_True;}
+  }
+
+  theSurface.Bounds(aUMin, aUMax, aTMin, aTMax);
+
+  if (thePtrDomain == NULL)
+    isMore = Standard_True;
+  else {
+    aPDomain->Init();
+    isMore = aPDomain->More();
+  }
+
+  while(isMore) {
+    // If the pointer to the domain is NULL, there is only one curve to treat:
+    // U isoline with the UMax parameter.
+
+    if (thePtrDomain == NULL)
+      theSurface.Load(Standard_False, GeomAbs_IsoU);
+    else
+      theSurface.Load(aPDomain->Value());
+
+    aTMin = theSurface.FirstParameter();
+    aTMax = theSurface.LastParameter();
+
+
+    // Get the spans on the curve.
+    Handle(TColStd_HArray1OfReal) aTKnots;
+    GProp_TFunction               aTFunc(theSurface, loc, IsByPoint, theCoeffs,
+					 aUMin, aCrvTol);
+
+    theSurface.GetTKnots(aTMin, aTMax, aTKnots);
+
+    Standard_Integer              iU             = aTKnots->Upper();
+    Standard_Integer              aNbTIntervals  = aTKnots->Length() - 1;
+    //Standard_Real                 aTolSpan       = aCrvTol/aNbTIntervals;
+    Standard_Real                 aTolSpan       = 0.9*theTolerance; //Relative error
+    math_KronrodSingleIntegration anIntegral;
+    GProp_ValueType               aValueType;
+
+
+    // Empirical criterion.
+    aNbPnts = Min(15, theSurface.IntegrationOrder()/aNbTIntervals + 1);
+    aNbPnts = Max(5, aNbPnts);
+//     aNbPnts = theSurface.IntegrationOrder();
+
+    aLocalValue.Init(0.);
+    aLocalTolReached.Init(0.);
+
+    for (k = 1; k <= aNbVal; k++) {
+
+      if(!CFlags(k)) continue;
+
+      Standard_Integer i = aTKnots->Lower();
+
+      switch (k) {
+      case 1:   aValueType = GProp_Mass;        break;
+      case 2:   aValueType = GProp_CenterMassX; break;
+      case 3:   aValueType = GProp_CenterMassY; break;
+      case 4:   aValueType = GProp_CenterMassZ; break;
+      case 5:   aValueType = GProp_InertiaXX;   break;
+      case 6:   aValueType = GProp_InertiaYY;   break;
+      case 7:   aValueType = GProp_InertiaZZ;   break;
+      case 8:   aValueType = GProp_InertiaXY;   break;
+      case 9:   aValueType = GProp_InertiaXZ;   break;
+      case 10:  aValueType = GProp_InertiaYZ;   break;
+
+      default: myErrorReached = -1.; return myErrorReached;
+      }
+      aTFunc.SetValueType(aValueType);
+
+      Standard_Real err1 = 0.;
+      while (i < iU) {
+
+	//cout << "-------------- Span " << i << " nbp: " << aNbPnts << endl;
+	Standard_Real aT1 = aTKnots->Value(i++);
+	Standard_Real aT2 = aTKnots->Value(i);
+
+	if(aT2 - aT1 < aTTol) continue;
+
+	aTFunc.SetNbKronrodPoints(aNbPnts);
+	aTFunc.Init();
+	aTFunc.SetTolerance(aCrvTol/(aT2-aT1));
+	anIntegral.Perform(aTFunc, aT1, aT2, aNbPnts, aTolSpan, aNbMaxIter);
+
+	if (!anIntegral.IsDone()) {
+	  myErrorReached = -1.;
+
+	  return myErrorReached;
+	}
+
+	aLocalValue(k)      += anIntegral.Value();
+	err1 = aTFunc.AbsolutError()*(aT2 - aT1); 
+	//cout << "Errors: " << anIntegral.NbIterReached() << " " << anIntegral.AbsolutError() << " " << err1 << endl;
+	aLocalTolReached(k) += anIntegral.AbsolutError() + err1; 
+	//cout << "--- Errors: " << anIntegral.NbIterReached() << " " << anIntegral.AbsolutError() << " " << err1 << endl;
+      }
+
+      aValue(k)      += aLocalValue(k);
+      aTolReached(k) += aLocalTolReached(k);
+    }
+
+    // If the pointer to the domain is NULL, there is only one curve to treat:
+    // U isoline with the UMax parameter.
+    if (thePtrDomain == NULL)
+      isMore = Standard_False;
+    else {
+      aPDomain->Next();
+      isMore = aPDomain->More();
+    }
+  }
+
+  // Get volume value.
+  dim            = aValue(1);
+  myErrorReached = aTolReached(1);
+  myAbsolutError = myErrorReached;
+  Standard_Real anAbsDim = Abs(dim);
+  Standard_Real aVolTol  = Epsilon(myAbsolutError);
+  if(anAbsDim >= aVolTol) myErrorReached /= anAbsDim;
+
+  if(theCGFlag || theIFlag) {
+    // Compute values of center of mass.
+    if(anAbsDim >= aVolTol) {
+      if (IsByPoint) {
+	aValue(2) = aCoeffs[0] + aValue(2)/dim;
+	aValue(3) = aCoeffs[1] + aValue(3)/dim;
+	aValue(4) = aCoeffs[2] + aValue(4)/dim;
+      } else {
+	aValue(2) /= dim;
+	aValue(3) /= dim;
+	aValue(4) /= dim;
+      }
+    } else {
+      aValue(2) = 0.;
+      aValue(3) = 0.;
+      aValue(4) = 0.;
+      dim = 0.;
+    }
+    g.SetCoord(aValue(2), aValue(3), aValue(4));
+  }
+
+  if(theIFlag) {
+    // Fill the matrix of inertia.
+    inertia.SetCols (gp_XYZ (aValue(5), aValue(8),  aValue(9)),
+		     gp_XYZ (aValue(8), aValue(6),  aValue(10)),
+		     gp_XYZ (aValue(9), aValue(10), aValue(7)));
+  }
+  //return myErrorReached;
+  return myAbsolutError;
+}
+ 

src/GProp/GProp_VGPropsGK.lxx

@@ -0,0 +1,35 @@
+// File:	GProp_VGPropsGK.lxx
+// Created:	Wed Dec 21 11:31:17 2005
+// Author:	Sergey KHROMOV
+//		<skv@dimox>
+
+
+//==========================================================================
+//function : Constructor
+//           Empty constructor.
+//==========================================================================
+
+inline GProp_VGPropsGK::GProp_VGPropsGK()
+     : myErrorReached(0.)
+{
+}
+
+//==========================================================================
+//function : SetLocation
+//           Sets the vertex that delimit 3D closed region of space.
+//==========================================================================
+
+inline void GProp_VGPropsGK::SetLocation(const gp_Pnt &theVertex)
+{
+  loc = theVertex;
+}
+
+//==========================================================================
+//function : GetErrorReached
+//           Returns the reached Error.
+//==========================================================================
+
+inline Standard_Real GProp_VGPropsGK::GetErrorReached() const
+{
+  return myErrorReached;
+}

src/GProp/GProp_VelGProps.cdl

@@ -0,0 +1,56 @@
+-- File:	GProp_VelGProps.cdl
+-- Created:	Wed Dec  2 16:14:27 1992
+-- Author:	Isabelle GRIGNON
+--		<isg@sdsun2>
+---Copyright:	 Matra Datavision 1992
+
+
+class VelGProps from GProp  inherits GProps
+
+        --- Purpose :
+        --  Computes the global properties of a geometric solid 
+        --  (3D closed region of space) 
+        --  The solid can be elementary(definition in the gp package) 
+
+
+uses   	Cone     from gp,
+	Cylinder from gp,
+        Pnt      from gp,
+	Sphere   from gp,
+	Torus    from gp
+
+
+is
+
+
+  Create returns VelGProps;
+
+  
+  Create (S : Cylinder; Alpha1, Alpha2, Z1, Z2 : Real; VLocation : Pnt)   
+     returns VelGProps;
+
+
+  Create (S : Cone; Alpha1, Alpha2, Z1, Z2 : Real; VLocation : Pnt) 
+     returns VelGProps;
+
+
+  Create (S : Sphere; Teta1, Teta2, Alpha1, Alpha2 : Real; VLocation : Pnt)
+     returns VelGProps;
+
+
+  Create (S : Torus; Teta1, Teta2, Alpha1, Alpha2 : Real; VLocation : Pnt)
+     returns VelGProps;
+
+
+  SetLocation(me : in out ;VLocation :Pnt);
+  
+  Perform(me : in out;S : Cylinder; Alpha1, Alpha2, Z1, Z2 : Real);
+
+  Perform(me : in out;S : Cone; Alpha1, Alpha2, Z1, Z2 : Real);
+
+  Perform(me : in out;S : Sphere; Teta1, Teta2, Alpha1, Alpha2 : Real);
+
+  Perform(me : in out;S : Torus; Teta1, Teta2, Alpha1, Alpha2 : Real);
+
+end VelGProps;
+

src/GProp/GProp_VelGProps.cxx

@@ -0,0 +1,369 @@
+#include <GProp_VelGProps.ixx>
+#include <Standard_NotImplemented.hxx>
+#include <gp.hxx>
+#include <GProp.hxx>
+#include <math_Matrix.hxx>
+#include <math_Vector.hxx>
+#include <math_Jacobi.hxx>
+
+
+GProp_VelGProps::GProp_VelGProps(){}
+
+void GProp_VelGProps::SetLocation(const gp_Pnt& VLocation)
+{
+  loc =VLocation;
+}
+
+
+GProp_VelGProps::GProp_VelGProps(const gp_Cylinder&    S,
+				 const Standard_Real Alpha1,
+				 const Standard_Real Alpha2,
+				 const Standard_Real Z1,
+				 const Standard_Real Z2,
+				 const gp_Pnt&        VLocation)
+{
+  SetLocation(VLocation);
+  Perform(S,Alpha1,Alpha2,Z1,Z2);
+}
+
+GProp_VelGProps::GProp_VelGProps(const gp_Cone&    S,
+				 const Standard_Real  Alpha1,
+				 const Standard_Real Alpha2,
+				 const Standard_Real Z1,
+				 const Standard_Real Z2,
+				 const gp_Pnt&        VLocation)
+{
+  SetLocation(VLocation);
+  Perform(S,Alpha1,Alpha2,Z1,Z2);
+}
+
+GProp_VelGProps::GProp_VelGProps(const gp_Sphere&    S,
+				 const Standard_Real Teta1,
+				 const Standard_Real Teta2,
+				 const Standard_Real Alpha1,
+				 const Standard_Real Alpha2,
+				 const gp_Pnt&        VLocation)
+{
+  SetLocation(VLocation);
+  Perform(S,Teta1,Teta2,Alpha1,Alpha2);
+}
+
+
+GProp_VelGProps::GProp_VelGProps(const gp_Torus&    S,
+				 const Standard_Real Teta1,
+				 const Standard_Real Teta2,
+				 const Standard_Real Alpha1,
+				 const Standard_Real Alpha2,
+				 const gp_Pnt&        VLocation)
+{
+  SetLocation(VLocation);
+  Perform(S,Teta1,Teta2,Alpha1,Alpha2);
+}
+
+void GProp_VelGProps::Perform(const gp_Cylinder&    S,
+			      const Standard_Real  Alpha1,
+			      const Standard_Real Alpha2,
+			      const Standard_Real Z1,
+			      const Standard_Real Z2)
+{
+  Standard_Real X0,Y0,Z0,Xa1,Ya1,Za1,Xa2,Ya2,Za2,Xa3,Ya3,Za3;
+  S.Location().Coord(X0,Y0,Z0);
+  Standard_Real Rayon = S.Radius();
+  S.Position().XDirection().Coord(Xa1,Ya1,Za1);
+  S.Position().YDirection().Coord(Xa2,Ya2,Za2);
+  S.Position().Direction().Coord(Xa3,Ya3,Za3);
+  dim = Rayon*Rayon*(Z2-Z1)/2.;
+  Standard_Real SA2 = Sin(Alpha2);
+  Standard_Real SA1 = Sin(Alpha1);
+  Standard_Real CA2 = Cos(Alpha2);
+  Standard_Real CA1 = Cos(Alpha1);
+  Standard_Real Dsin = SA2-SA1; 
+  Standard_Real Dcos = CA1-CA2; 
+  Standard_Real Coef = Rayon/(Alpha2-Alpha1);
+
+  g.SetCoord(X0+(Coef*(Xa1*Dsin+Xa2*Dcos) ) + (Xa3*(Z2+Z1)/2.),
+	     Y0+(Coef*(Ya1*Dsin+Ya2*Dcos) ) + (Ya3*(Z2+Z1)/2.),
+	     Z0+(Coef*(Za1*Dsin+Za2*Dcos) ) + (Za3*(Z2+Z1)/2.) );
+  
+  Standard_Real ICn2 = dim/2. *( Alpha2-Alpha1 + SA2*CA2 - SA1*CA1 );  
+  Standard_Real ISn2 = dim/2. *( Alpha2-Alpha1 - SA2*CA2 + SA1*CA1 );
+  Standard_Real IZ2 = dim * (Alpha2-Alpha1)*(Z2*Z2+Z1*Z2+Z1*Z1);
+  Standard_Real ICnSn = dim *(CA2*CA2-CA1*CA1)/2.;
+  Standard_Real ICnz = dim *(Z2+Z1)/2.*Dsin;
+  Standard_Real ISnz = dim *(Z2+Z1)/2.*Dcos;
+  dim =(Alpha2-Alpha1)*dim;  
+
+  math_Matrix Dm(1,3,1,3);
+
+  Dm(1,1) = Rayon*Rayon*ISn2 + IZ2;
+  Dm(2,2) = Rayon*Rayon*ICn2 + IZ2;
+  Dm(3,3) = Rayon*Rayon*dim;
+  Dm(1,2) = Dm(2,1) = -Rayon*Rayon*ICnSn;
+  Dm(1,3) = Dm(3,1) = -Rayon*ICnz;
+  Dm(3,2) = Dm(2,3) = -Rayon*ISnz;
+
+  math_Matrix Passage (1,3,1,3);
+  Passage(1,1) = Xa1; Passage(1,2) = Xa2 ;Passage(1,3) = Xa3;
+  Passage(2,1) = Ya1; Passage(2,2) = Ya2 ;Passage(2,3) = Ya3;
+  Passage(3,1) = Za1; Passage(3,2) = Za2 ;Passage(3,3) = Za3;
+
+  math_Jacobi J(Dm);
+  math_Vector V1(1,3),V2(1,3),V3(1,3);
+  J.Vector(1,V1);
+  V1.Multiply(Passage,V1);
+  V1.Multiply(J.Value(1));
+  J.Vector(2,V2);
+  V2.Multiply(Passage,V2);
+  V2.Multiply(J.Value(2));
+  J.Vector(3,V3);
+  V3.Multiply(Passage,V3);
+  V3.Multiply(J.Value(3));
+
+  inertia = gp_Mat (gp_XYZ(V1(1),V2(1),V3(1)),
+		    gp_XYZ(V1(2),V2(2),V3(2)),
+		    gp_XYZ(V1(3),V2(3),V3(3)));
+  gp_Mat Hop;
+  GProp::HOperator(g,loc,dim,Hop);
+  inertia = inertia+Hop;
+}
+
+void GProp_VelGProps::Perform(const gp_Cone&    S,
+			      const Standard_Real Alpha1,
+			      const Standard_Real Alpha2,
+			      const Standard_Real Z1,
+			      const Standard_Real Z2)
+{
+  Standard_Real X0,Y0,Z0,Xa1,Ya1,Za1,Xa2,Ya2,Za2,Xa3,Ya3,Za3;
+  S.Location().Coord(X0,Y0,Z0);
+  S.Position().XDirection().Coord(Xa1,Ya1,Za1);
+  S.Position().YDirection().Coord(Xa2,Ya2,Za2);
+  S.Position().Direction().Coord(Xa3,Ya3,Za3);
+  Standard_Real t =S.SemiAngle();
+  Standard_Real Cnt = Cos(t);
+  Standard_Real Snt = Sin(t); 
+  Standard_Real R = S.RefRadius();
+  Standard_Real Sn2 = Sin(Alpha2);
+  Standard_Real Sn1 = Sin(Alpha1);
+  Standard_Real Cn2 = Cos(Alpha2);
+  Standard_Real Cn1 = Cos(Alpha1);
+  Standard_Real ZZ = (Z2-Z1)*(Z2-Z1)*Cnt*Snt;
+  Standard_Real Auxi1=  2*R +(Z2+Z1)*Snt;
+
+  dim = ZZ*(Alpha2-Alpha1)*Auxi1/2.;
+
+  Standard_Real R1 = R + Z1*Snt;
+  Standard_Real R2 = R + Z2*Snt;
+  Standard_Real Coef0 = (R1*R1+R1*R2+R2*R2);
+  Standard_Real Iz = Cnt*(R*(Z2+Z1) + 2*Snt*(Z1*Z1+Z1*Z2+Z2*Z2)/3.)/Auxi1;
+  Standard_Real Ix = Coef0*(Sn2-Sn1)/(Alpha2-Alpha1)/Auxi1;
+  Standard_Real Iy = Coef0*(Cn1-Cn2)/(Alpha2-Alpha1)/Auxi1;
+
+  g.SetCoord(X0 + Xa1*Ix + Xa2*Iy  + Xa3*Iz,
+	     Y0 + Ya1*Ix + Ya2*Iy  + Ya3*Iz,
+	     Z0 + Za1*Ix + Za2*Iy  + Za3*Iz);
+
+  Standard_Real IR2 = ZZ*(R2*R2*R2+R2*R2*R1+R1*R1*R2+R1*R1*R1)/4.;
+  Standard_Real ICn2  = IR2*(Alpha2-Alpha1+Cn2*Sn2-Cn1*Sn1)/2.;
+  Standard_Real ISn2 = IR2*(Alpha2-Alpha1+Cn2*Sn2-Cn1*Sn1)/2.;
+  Standard_Real IZ2  = ZZ*Cnt*Cnt*(Alpha2-Alpha1)*
+                       (Z1*Z1*(R/3 + Z1*Snt/4) + 
+			Z2*Z2*(R/3 + Z2*Snt/4) +
+			Z1*Z2*(R/3 +Z1*Snt/4 +Z2*Snt/4));
+  Standard_Real ICnSn = IR2*(Cn2*Cn2-Cn1*Cn1);
+  Standard_Real ICnz = (Z1+Z2)*ZZ*Coef0*(Sn2-Sn1)/3;
+  Standard_Real ISnz = (Z1+Z2)*ZZ*Coef0*(Cn1-Cn2)/3;    
+
+   
+  math_Matrix Dm(1,3,1,3);
+  Dm(1,1) = ISn2 + IZ2;
+  Dm(2,2) = ICn2 + IZ2;
+  Dm(3,3) = IR2*(Alpha2-Alpha1);
+  Dm(1,2) = Dm(2,1) = -ICnSn;
+  Dm(1,3) = Dm(3,1) = -ICnz;
+  Dm(3,2) = Dm(2,3) = -ISnz;
+
+  math_Matrix Passage (1,3,1,3);
+  Passage(1,1) = Xa1; Passage(1,2) = Xa2 ;Passage(1,3) = Xa3;
+  Passage(2,1) = Ya1; Passage(2,2) = Ya2 ;Passage(2,3) = Ya3;
+  Passage(3,1) = Za1; Passage(3,2) = Za2 ;Passage(3,3) = Za3;
+
+  math_Jacobi J(Dm);
+  math_Vector V1(1,3),V2(1,3),V3(1,3);
+  J.Vector(1,V1);
+  V1.Multiply(Passage,V1);
+  V1.Multiply(J.Value(1));
+  J.Vector(2,V2);
+  V2.Multiply(Passage,V2);
+  V2.Multiply(J.Value(2));
+  J.Vector(3,V3);
+  V3.Multiply(Passage,V3);
+  V3.Multiply(J.Value(3));
+
+  inertia = gp_Mat (gp_XYZ(V1(1),V2(1),V3(1)),
+		    gp_XYZ(V1(2),V2(2),V3(2)),
+		    gp_XYZ(V1(3),V2(3),V3(3)));
+  gp_Mat Hop;
+  GProp::HOperator(g,loc,dim,Hop);
+  inertia = inertia+Hop;
+}
+
+void GProp_VelGProps::Perform(const gp_Sphere&    S,
+			      const Standard_Real Teta1,
+			      const Standard_Real Teta2,
+			      const Standard_Real Alpha1,
+			      const Standard_Real Alpha2)
+{
+  Standard_Real X0,Y0,Z0,Xa1,Ya1,Za1,Xa2,Ya2,Za2,Xa3,Ya3,Za3;
+  S.Location().Coord(X0,Y0,Z0);
+  S.Position().XDirection().Coord(Xa1,Ya1,Za1);
+  S.Position().YDirection().Coord(Xa2,Ya2,Za2);
+  S.Position().Direction().Coord(Xa3,Ya3,Za3);
+  Standard_Real R = S.Radius();
+  Standard_Real Cnt1 = Cos(Teta1);
+  Standard_Real Snt1 = Sin(Teta1);
+  Standard_Real Cnt2 = Cos(Teta2);
+  Standard_Real Snt2 = Sin(Teta2);
+  Standard_Real Cnf1 = Cos(Alpha1);
+  Standard_Real Snf1 = Sin(Alpha1);
+  Standard_Real Cnf2 = Cos(Alpha2);
+  Standard_Real Snf2 = Sin(Alpha2);
+   
+  dim = (Teta2-Teta1)*R*R*R*(Snf2-Snf1)/3.;   
+
+  Standard_Real Ix = 
+    R*(Snt2-Snt1)/(Teta2-Teta1)*
+    (Alpha2-Alpha1+Snf2*Cnf2-Snf1*Cnf1)/(Snf2-Snf1)/2.;
+  Standard_Real Iy = 
+    R*(Cnt1-Cnt2)/(Teta2-Teta1)*
+    (Alpha2-Alpha1+Snf2*Cnf2-Snf1*Cnf1)/(Snf2-Snf1)/2.;
+  Standard_Real Iz = R*(Snf2+Snf1)/2.;
+  g.SetCoord(
+   X0 + Ix*Xa1 + Iy*Xa2 + Iz*Xa3,
+   Y0 + Ix*Ya1 + Iy*Ya2 + Iz*Ya3,
+   Z0 + Ix*Za1 + Iy*Za2 + Iz*Za3);
+
+  Standard_Real IR2 = ( Cnf2*Snf2*(Cnf2+1.)- Cnf1*Snf1*(Cnf1+1.) +
+                       Alpha2-Alpha1 )/9.;
+  Standard_Real ICn2 = (Teta2-Teta1+ Cnt2*Snt2-Cnt1*Snt1)*IR2/2.;
+  Standard_Real ISn2 = (Teta2-Teta1-Cnt2*Snt2+Cnt1*Snt1)*IR2/2.;
+  Standard_Real ICnSn = ( Snt2*Snt2-Snt1*Snt1)*IR2/2.;
+  Standard_Real IZ2 = (Teta2-Teta1)*(Snf2*Snf2*Snf2-Snf1*Snf1*Snf1)/9.;
+  Standard_Real ICnz =(Snt2-Snt1)*(Cnf1*Cnf1*Cnf1-Cnf2*Cnf2*Cnf2)/9.;
+  Standard_Real ISnz =(Cnt1-Cnt2)*(Cnf1*Cnf1*Cnf1-Cnf2*Cnf2*Cnf2)/9.;
+
+  math_Matrix Dm(1,3,1,3);
+  Dm(1,1) = ISn2 +IZ2; 
+  Dm(2,2) = ICn2 +IZ2;
+  Dm(3,3) = IR2*(Teta2-Teta1);
+  Dm(1,2) = Dm(2,1) = -ICnSn;
+  Dm(1,3) = Dm(3,1) = -ICnz;
+  Dm(3,2) = Dm(2,3) = -ISnz;
+
+  math_Matrix Passage (1,3,1,3);
+  Passage(1,1) = Xa1; Passage(1,2) = Xa2 ;Passage(1,3) = Xa3;
+  Passage(2,1) = Ya1; Passage(2,2) = Ya2 ;Passage(2,3) = Ya3;
+  Passage(3,1) = Za1; Passage(3,2) = Za2 ;Passage(3,3) = Za3;
+
+  math_Jacobi J(Dm);
+  R = R*R*R*R*R;
+  math_Vector V1(1,3), V2(1,3), V3(1,3);
+  J.Vector(1,V1);
+  V1.Multiply(Passage,V1);
+  V1.Multiply(R*J.Value(1));
+  J.Vector(2,V2);
+  V2.Multiply(Passage,V2);
+  V2.Multiply(R*J.Value(2));
+  J.Vector(3,V3);
+  V3.Multiply(Passage,V3);
+  V3.Multiply(R*J.Value(3));
+
+  inertia = gp_Mat (gp_XYZ(V1(1),V2(1),V3(1)),
+		    gp_XYZ(V1(2),V2(2),V3(2)),
+		    gp_XYZ(V1(3),V2(3),V3(3)));
+  gp_Mat Hop;
+  GProp::HOperator(g,loc,dim,Hop);
+  inertia = inertia+Hop;
+
+}
+
+
+void GProp_VelGProps::Perform(const gp_Torus&    S,
+			      const Standard_Real Teta1,
+			      const Standard_Real Teta2,
+			      const Standard_Real Alpha1,
+			      const Standard_Real Alpha2)
+{
+  Standard_Real X0,Y0,Z0,Xa1,Ya1,Za1,Xa2,Ya2,Za2,Xa3,Ya3,Za3;
+  S.Location().Coord(X0,Y0,Z0);
+  S.Position().XDirection().Coord(Xa1,Ya1,Za1);
+  S.Position().YDirection().Coord(Xa2,Ya2,Za2);
+  S.Position().Direction().Coord(Xa3,Ya3,Za3);
+  Standard_Real RMax = S.MajorRadius();
+  Standard_Real Rmin = S.MinorRadius();
+  Standard_Real Cnt1 = Cos(Teta1);
+  Standard_Real Snt1 = Sin(Teta1);
+  Standard_Real Cnt2 = Cos(Alpha2);
+  Standard_Real Snt2 = Sin(Alpha2);
+  Standard_Real Cnf1 = Cos(Alpha1);
+  Standard_Real Snf1 = Sin(Alpha1);
+  Standard_Real Cnf2 = Cos(Alpha2);
+  Standard_Real Snf2 = Sin(Alpha2);
+
+  dim = RMax*Rmin*Rmin*(Teta2-Teta1)*(Alpha2-Alpha1)/2.;
+  Standard_Real Ix = 
+    (Snt2-Snt1)/(Teta2-Teta1)*(Rmin*(Snf2-Snf1)/(Alpha2-Alpha1) + RMax);
+  Standard_Real Iy = 
+    (Cnt1-Cnt2)/(Teta2-Teta1)*(Rmin*(Snf2-Snf1)/(Alpha2-Alpha1) + RMax);
+  Standard_Real Iz = Rmin*(Cnf1-Cnf2)/(Alpha2-Alpha1);
+  
+  g.SetCoord(
+	     X0+Ix*Xa1+Iy*Xa2+Iz*Xa3,
+	     Y0+Ix*Ya1+Iy*Ya2+Iz*Ya3,
+	     Z0+Ix*Za1+Iy*Za2+Iz*Za3);
+
+  Standard_Real IR2 =  RMax*RMax+Rmin*Rmin/2. +2.*RMax*Rmin*(Snf2-Snf1) +
+   		       Rmin*Rmin/2.*(Snf2*Cnf2-Snf1*Cnf1);
+  Standard_Real ICn2 = IR2*(Teta2-Teta1 +Snt2*Cnt2-Snt1*Cnt1)/2.;
+  Standard_Real ISn2 = IR2*(Teta2-Teta1 -Snt2*Cnt2+Snt1*Cnt1)/2.;
+  Standard_Real ICnSn = IR2*(Snt2*Snt2-Snt1*Snt1)/2.;
+  Standard_Real IZ2 = 
+     (Teta2-Teta1)*Rmin*Rmin*(Alpha2-Alpha1-Snf2*Cnf2+Snf1*Cnf1)/2.;
+  Standard_Real ICnz = Rmin*(Snt2-Snt1)*(Cnf1-Cnf2)*(RMax+Rmin*(Cnf1+Cnf2)/2.);
+  Standard_Real ISnz = Rmin*(Cnt2-Cnt1)*(Cnf1-Cnf2)*(RMax+Rmin*(Cnf1+Cnf2)/2.);
+
+  math_Matrix Dm(1,3,1,3);
+  Dm(1,1) = ISn2 + IZ2; 
+  Dm(2,2) = ICn2 + IZ2;
+  Dm(3,3) = IR2*(Teta2-Teta1);
+  Dm(1,2) = Dm(2,1) = -ICnSn;
+  Dm(1,3) = Dm(3,1) = -ICnz;
+  Dm(3,2) = Dm(2,3) = -ISnz;
+
+  math_Matrix Passage (1,3,1,3);
+  Passage(1,1) = Xa1; Passage(1,2) = Xa2 ;Passage(1,3) = Xa3;
+  Passage(2,1) = Ya1; Passage(2,2) = Ya2 ;Passage(2,3) = Ya3;
+  Passage(3,1) = Za1; Passage(3,2) = Za2 ;Passage(3,3) = Za3;
+
+  math_Jacobi J(Dm);
+  RMax = RMax*Rmin*Rmin/2.;
+  math_Vector V1(1,3), V2(1,3), V3(1,3);
+  J.Vector(1,V1);
+  V1.Multiply(Passage,V1);
+  V1.Multiply(RMax*J.Value(1));
+  J.Vector(2,V2);
+  V2.Multiply(Passage,V2);
+  V2.Multiply(RMax*J.Value(2));
+  J.Vector(3,V3);
+  V3.Multiply(Passage,V3);
+  V3.Multiply(RMax*J.Value(3));
+
+  inertia = gp_Mat (gp_XYZ(V1(1),V2(1),V3(1)),
+		    gp_XYZ(V1(2),V2(2),V3(2)),
+		    gp_XYZ(V1(3),V2(3),V3(3)));
+  gp_Mat Hop;
+  GProp::HOperator(g,loc,dim,Hop);
+  inertia = inertia+Hop;
+}
+
+
+