[occt.git] / src / BRepGProp / BRepGProp_VinertGK.cxx

// Created on: 2005-12-09
// Created by: Sergey KHROMOV
// Copyright (c) 2005-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#include <BRepGProp_VinertGK.ixx>

#include <TColStd_HArray1OfReal.hxx>
#include <TColStd_Array1OfBoolean.hxx>
#include <math_KronrodSingleIntegration.hxx>

#include <BRepGProp_TFunction.hxx>

//==========================================================================
//function : Constructor
//==========================================================================

BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_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
//           
//==========================================================================

BRepGProp_VinertGK::BRepGProp_VinertGK(      BRepGProp_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
//           
//==========================================================================

BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face        &theSurface,
                                       BRepGProp_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
//           
//==========================================================================

BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face        &theSurface,
                                       BRepGProp_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
//           
//==========================================================================

BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_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
//           
//==========================================================================

BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face        &theSurface,
                                       BRepGProp_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 BRepGProp_VinertGK::Perform(BRepGProp_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 BRepGProp_VinertGK::Perform(BRepGProp_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 BRepGProp_VinertGK::Perform(BRepGProp_Face        &theSurface,
                                          BRepGProp_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 BRepGProp_VinertGK::Perform(BRepGProp_Face        &theSurface,
                                          BRepGProp_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 BRepGProp_VinertGK::Perform(BRepGProp_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 BRepGProp_VinertGK::Perform(BRepGProp_Face        &theSurface,
                                          BRepGProp_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 BRepGProp_VinertGK::PrivatePerform
(BRepGProp_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.
  BRepGProp_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 = (BRepGProp_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;
    BRepGProp_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;
}
Commit	Line	Data
b311480e	1	// Created on: 2005-12-09
b311480e	2	// Created by: Sergey KHROMOV
973c2be1	3	// Copyright (c) 2005-2014 OPEN CASCADE SAS
b311480e	4	//
973c2be1	5	// This file is part of Open CASCADE Technology software library.
b311480e	6	//
d5f74e42	7	// This library is free software; you can redistribute it and/or modify it under
d5f74e42	8	// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1	9	// by the Free Software Foundation, with special exception defined in the file
	10	// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
	11	// distribution for complete text of the license and disclaimer of any warranty.
b311480e	12	//
973c2be1	13	// Alternatively, this file may be used under the terms of Open CASCADE
973c2be1	14	// commercial license or contractual agreement.
7fd59977	15
424cd6bb	16	#include <BRepGProp_VinertGK.ixx>
424cd6bb	17
7fd59977	18	#include <TColStd_HArray1OfReal.hxx>
	19	#include <TColStd_Array1OfBoolean.hxx>
	20	#include <math_KronrodSingleIntegration.hxx>
424cd6bb	21
424cd6bb	22	#include <BRepGProp_TFunction.hxx>
7fd59977	23
	24	//==========================================================================
	25	//function : Constructor
7fd59977	26	//==========================================================================
7fd59977	27
424cd6bb	28	BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face &theSurface,
	29	const gp_Pnt &theLocation,
	30	const Standard_Real theTolerance,
	31	const Standard_Boolean theCGFlag,
	32	const Standard_Boolean theIFlag):
	33	myErrorReached(0.)
7fd59977	34	{
	35	SetLocation(theLocation);
	36	Perform(theSurface, theTolerance, theCGFlag, theIFlag);
	37	}
	38
	39	//==========================================================================
	40	//function : Constructor
	41	//
	42	//==========================================================================
	43
424cd6bb	44	BRepGProp_VinertGK::BRepGProp_VinertGK( BRepGProp_Face &theSurface,
	45	const gp_Pnt &thePoint,
	46	const gp_Pnt &theLocation,
	47	const Standard_Real theTolerance,
	48	const Standard_Boolean theCGFlag,
	49	const Standard_Boolean theIFlag):
	50	myErrorReached(0.)
7fd59977	51	{
	52	SetLocation(theLocation);
	53	Perform(theSurface, thePoint, theTolerance, theCGFlag, theIFlag);
	54	}
	55
	56	//==========================================================================
	57	//function : Constructor
	58	//
	59	//==========================================================================
	60
424cd6bb	61	BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face &theSurface,
	62	BRepGProp_Domain &theDomain,
	63	const gp_Pnt &theLocation,
	64	const Standard_Real theTolerance,
	65	const Standard_Boolean theCGFlag,
	66	const Standard_Boolean theIFlag):
	67	myErrorReached(0.)
7fd59977	68	{
	69	SetLocation(theLocation);
	70	Perform(theSurface, theDomain, theTolerance, theCGFlag, theIFlag);
	71	}
	72
	73	//==========================================================================
	74	//function : Constructor
	75	//
	76	//==========================================================================
	77
424cd6bb	78	BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face &theSurface,
	79	BRepGProp_Domain &theDomain,
	80	const gp_Pnt &thePoint,
	81	const gp_Pnt &theLocation,
	82	const Standard_Real theTolerance,
	83	const Standard_Boolean theCGFlag,
	84	const Standard_Boolean theIFlag):
	85	myErrorReached(0.)
7fd59977	86	{
	87	SetLocation(theLocation);
	88	Perform(theSurface, theDomain, thePoint, theTolerance, theCGFlag, theIFlag);
	89	}
	90
	91	//==========================================================================
	92	//function : Constructor
	93	//
	94	//==========================================================================
	95
424cd6bb	96	BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face &theSurface,
	97	const gp_Pln &thePlane,
	98	const gp_Pnt &theLocation,
	99	const Standard_Real theTolerance,
	100	const Standard_Boolean theCGFlag,
	101	const Standard_Boolean theIFlag):
	102	myErrorReached(0.)
7fd59977	103	{
	104	SetLocation(theLocation);
	105	Perform(theSurface, thePlane, theTolerance, theCGFlag, theIFlag);
	106	}
	107
	108	//==========================================================================
	109	//function : Constructor
	110	//
	111	//==========================================================================
	112
424cd6bb	113	BRepGProp_VinertGK::BRepGProp_VinertGK(BRepGProp_Face &theSurface,
	114	BRepGProp_Domain &theDomain,
	115	const gp_Pln &thePlane,
	116	const gp_Pnt &theLocation,
	117	const Standard_Real theTolerance,
	118	const Standard_Boolean theCGFlag,
	119	const Standard_Boolean theIFlag):
	120	myErrorReached(0.)
7fd59977	121	{
	122	SetLocation(theLocation);
	123	Perform(theSurface, theDomain, thePlane, theTolerance, theCGFlag, theIFlag);
	124	}
	125
	126	//==========================================================================
	127	//function : Perform
	128	// Compute the properties.
	129	//==========================================================================
	130
424cd6bb	131	Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face &theSurface,
	132	const Standard_Real theTolerance,
	133	const Standard_Boolean theCGFlag,
	134	const Standard_Boolean theIFlag)
7fd59977	135
	136	{
	137	Standard_Real aShift[] = { 0., 0., 0. };
	138
	139	return PrivatePerform(theSurface, NULL, Standard_True, &aShift, theTolerance,
424cd6bb	140	theCGFlag, theIFlag);
7fd59977	141	}
	142
	143	//==========================================================================
	144	//function : Perform
	145	// Compute the properties.
	146	//==========================================================================
	147
424cd6bb	148	Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face &theSurface,
	149	const gp_Pnt &thePoint,
	150	const Standard_Real theTolerance,
	151	const Standard_Boolean theCGFlag,
	152	const Standard_Boolean theIFlag)
7fd59977	153
	154	{
	155	gp_XYZ aXYZ(thePoint.XYZ().Subtracted(loc.XYZ()));
	156	Standard_Real aShift[3];
	157
	158	aXYZ.Coord(aShift[0], aShift[1], aShift[2]);
	159
	160	return PrivatePerform(theSurface, NULL, Standard_True, &aShift, theTolerance,
424cd6bb	161	theCGFlag, theIFlag);
7fd59977	162	}
	163
	164	//==========================================================================
	165	//function : Perform
	166	// Compute the properties.
	167	//==========================================================================
	168
424cd6bb	169	Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face &theSurface,
	170	BRepGProp_Domain &theDomain,
	171	const Standard_Real theTolerance,
	172	const Standard_Boolean theCGFlag,
	173	const Standard_Boolean theIFlag)
7fd59977	174
	175	{
	176	Standard_Real aShift[] = { 0., 0., 0. };
	177
	178	return PrivatePerform(theSurface, &theDomain,
424cd6bb	179	Standard_True, &aShift, theTolerance,
424cd6bb	180	theCGFlag, theIFlag);
7fd59977	181	}
	182
	183	//==========================================================================
	184	//function : Perform
	185	// Compute the properties.
	186	//==========================================================================
	187
424cd6bb	188	Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face &theSurface,
	189	BRepGProp_Domain &theDomain,
	190	const gp_Pnt &thePoint,
	191	const Standard_Real theTolerance,
	192	const Standard_Boolean theCGFlag,
	193	const Standard_Boolean theIFlag)
7fd59977	194
	195	{
	196	gp_XYZ aXYZ(thePoint.XYZ().Subtracted(loc.XYZ()));
	197	Standard_Real aShift[3];
	198
	199	aXYZ.Coord(aShift[0], aShift[1], aShift[2]);
	200
	201	return PrivatePerform(theSurface, &theDomain,
424cd6bb	202	Standard_True, &aShift, theTolerance,
424cd6bb	203	theCGFlag, theIFlag);
7fd59977	204	}
	205
	206	//==========================================================================
	207	//function : Perform
	208	// Compute the properties.
	209	//==========================================================================
	210
424cd6bb	211	Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face &theSurface,
	212	const gp_Pln &thePlane,
	213	const Standard_Real theTolerance,
	214	const Standard_Boolean theCGFlag,
	215	const Standard_Boolean theIFlag)
7fd59977	216
	217	{
	218	Standard_Real aCoeff[4];
	219	Standard_Real aXLoc;
	220	Standard_Real aYLoc;
	221	Standard_Real aZLoc;
	222
	223	loc.Coord(aXLoc, aYLoc, aZLoc);
	224	thePlane.Coefficients (aCoeff[0], aCoeff[1], aCoeff[2], aCoeff[3]);
	225	aCoeff[3] = aCoeff[3] - aCoeff[0]aXLoc - aCoeff[1]aYLoc - aCoeff[2]*aZLoc;
	226
	227	return PrivatePerform(theSurface, NULL,
424cd6bb	228	Standard_False, &aCoeff, theTolerance,
424cd6bb	229	theCGFlag, theIFlag);
7fd59977	230	}
	231
	232	//==========================================================================
	233	//function : Perform
	234	// Compute the properties.
	235	//==========================================================================
	236
424cd6bb	237	Standard_Real BRepGProp_VinertGK::Perform(BRepGProp_Face &theSurface,
	238	BRepGProp_Domain &theDomain,
	239	const gp_Pln &thePlane,
	240	const Standard_Real theTolerance,
	241	const Standard_Boolean theCGFlag,
	242	const Standard_Boolean theIFlag)
7fd59977	243
	244	{
	245	Standard_Real aCoeff[4];
	246	Standard_Real aXLoc;
	247	Standard_Real aYLoc;
	248	Standard_Real aZLoc;
	249
	250	loc.Coord(aXLoc, aYLoc, aZLoc);
	251	thePlane.Coefficients (aCoeff[0], aCoeff[1], aCoeff[2], aCoeff[3]);
	252	aCoeff[3] = aCoeff[3] - aCoeff[0]aXLoc - aCoeff[1]aYLoc - aCoeff[2]*aZLoc;
	253
	254	return PrivatePerform(theSurface, &theDomain,
424cd6bb	255	Standard_False, &aCoeff, theTolerance,
424cd6bb	256	theCGFlag, theIFlag);
7fd59977	257	}
	258
	259	//==========================================================================
	260	//function : PrivatePerform
	261	// Compute the properties.
	262	//==========================================================================
	263
424cd6bb	264	Standard_Real BRepGProp_VinertGK::PrivatePerform
	265	(BRepGProp_Face &theSurface,
	266	const Standard_Address thePtrDomain,
	267	const Standard_Boolean IsByPoint,
	268	const Standard_Address theCoeffs,
	269	const Standard_Real theTolerance,
	270	const Standard_Boolean theCGFlag,
	271	const Standard_Boolean theIFlag)
7fd59977	272
	273	{
	274
	275	const Standard_Real aTTol = 1.e-9;
	276	Standard_Real aCoeffs = (Standard_Real )theCoeffs;
	277
	278	// Compute the number of 2d bounding curves of the face.
424cd6bb	279	BRepGProp_Domain *aPDomain = NULL;
7fd59977	280	Standard_Integer aNbCurves = 0;
	281
	282	// If the pointer to the domain is NULL, there is only one curve to treat:
	283	// U isoline with the UMax parameter.
	284	if (thePtrDomain == NULL)
	285	aNbCurves = 1;
	286	else {
424cd6bb	287	aPDomain = (BRepGProp_Domain *)thePtrDomain;
7fd59977	288
	289	for (aPDomain->Init(); aPDomain->More(); aPDomain->Next())
	290	aNbCurves++;
	291	}
	292
	293	if (aNbCurves == 0) {
	294	myErrorReached = -1.;
	295
	296	return myErrorReached;
	297	}
	298
	299	//Standard_Real aCrvTol = 0.5*theTolerance/aNbCurves;
	300	Standard_Real aCrvTol = 0.1*theTolerance;
	301	Standard_Real aUMin;
	302	Standard_Real aUMax;
	303	Standard_Real aTMin;
	304	Standard_Real aTMax;
	305	Standard_Integer aNbPnts;
	306	Standard_Integer aNbMaxIter = 1000;
	307	Standard_Integer aNbVal = 10;
	308	Standard_Integer k;
	309	math_Vector aLocalValue(1, aNbVal);
	310	math_Vector aLocalTolReached(1, aNbVal);
	311	math_Vector aValue(1, aNbVal);
	312	math_Vector aTolReached(1, aNbVal);
	313	TColStd_Array1OfBoolean CFlags(1, aNbVal);
	314	CFlags.Init(Standard_False);
	315	Standard_Boolean isMore;
	316
	317	//aNbVal = 1;
	318	aValue.Init(0.);
	319	aTolReached.Init(0.);
	320
	321	CFlags.Init(Standard_False);
	322	CFlags(1) = Standard_True;
	323
	324	if(theCGFlag \|\| theIFlag) {
	325	Standard_Integer i;
	326	for(i = 2; i <= 4; ++i) {CFlags(i) = Standard_True;}
	327	}
	328
	329	if(theIFlag) {
	330	Standard_Integer i;
	331	for(i = 5; i <= 10; ++i) {CFlags(i) = Standard_True;}
	332	}
	333
	334	theSurface.Bounds(aUMin, aUMax, aTMin, aTMax);
	335
	336	if (thePtrDomain == NULL)
	337	isMore = Standard_True;
	338	else {
	339	aPDomain->Init();
	340	isMore = aPDomain->More();
	341	}
	342
	343	while(isMore) {
	344	// If the pointer to the domain is NULL, there is only one curve to treat:
	345	// U isoline with the UMax parameter.
	346
	347	if (thePtrDomain == NULL)
	348	theSurface.Load(Standard_False, GeomAbs_IsoU);
	349	else
	350	theSurface.Load(aPDomain->Value());
	351
352	aTMin = theSurface.FirstParameter();
353	aTMax = theSurface.LastParameter();
354
355
356	// Get the spans on the curve.
357	Handle(TColStd_HArray1OfReal) aTKnots;
424cd6bb	358	BRepGProp_TFunction aTFunc(theSurface, loc, IsByPoint, theCoeffs,
424cd6bb	359	aUMin, aCrvTol);
7fd59977	360
	361	theSurface.GetTKnots(aTMin, aTMax, aTKnots);
	362
	363	Standard_Integer iU = aTKnots->Upper();
	364	Standard_Integer aNbTIntervals = aTKnots->Length() - 1;
	365	//Standard_Real aTolSpan = aCrvTol/aNbTIntervals;
	366	Standard_Real aTolSpan = 0.9*theTolerance; //Relative error
	367	math_KronrodSingleIntegration anIntegral;
	368	GProp_ValueType aValueType;
	369
	370
	371	// Empirical criterion.
	372	aNbPnts = Min(15, theSurface.IntegrationOrder()/aNbTIntervals + 1);
	373	aNbPnts = Max(5, aNbPnts);
424cd6bb	374	// aNbPnts = theSurface.IntegrationOrder();
7fd59977	375
	376	aLocalValue.Init(0.);
	377	aLocalTolReached.Init(0.);
	378
	379	for (k = 1; k <= aNbVal; k++) {
	380
	381	if(!CFlags(k)) continue;
	382
	383	Standard_Integer i = aTKnots->Lower();
	384
	385	switch (k) {
	386	case 1: aValueType = GProp_Mass; break;
	387	case 2: aValueType = GProp_CenterMassX; break;
	388	case 3: aValueType = GProp_CenterMassY; break;
	389	case 4: aValueType = GProp_CenterMassZ; break;
	390	case 5: aValueType = GProp_InertiaXX; break;
	391	case 6: aValueType = GProp_InertiaYY; break;
	392	case 7: aValueType = GProp_InertiaZZ; break;
	393	case 8: aValueType = GProp_InertiaXY; break;
	394	case 9: aValueType = GProp_InertiaXZ; break;
	395	case 10: aValueType = GProp_InertiaYZ; break;
	396
	397	default: myErrorReached = -1.; return myErrorReached;
	398	}
	399	aTFunc.SetValueType(aValueType);
	400
	401	Standard_Real err1 = 0.;
	402	while (i < iU) {
	403
424cd6bb	404	//cout << "-------------- Span " << i << " nbp: " << aNbPnts << endl;
	405	Standard_Real aT1 = aTKnots->Value(i++);
	406	Standard_Real aT2 = aTKnots->Value(i);
7fd59977	407
424cd6bb	408	if(aT2 - aT1 < aTTol) continue;
7fd59977	409
424cd6bb	410	aTFunc.SetNbKronrodPoints(aNbPnts);
	411	aTFunc.Init();
	412	aTFunc.SetTolerance(aCrvTol/(aT2-aT1));
	413	anIntegral.Perform(aTFunc, aT1, aT2, aNbPnts, aTolSpan, aNbMaxIter);
7fd59977	414
424cd6bb	415	if (!anIntegral.IsDone()) {
424cd6bb	416	myErrorReached = -1.;
7fd59977	417
424cd6bb	418	return myErrorReached;
424cd6bb	419	}
7fd59977	420
424cd6bb	421	aLocalValue(k) += anIntegral.Value();
	422	err1 = aTFunc.AbsolutError()*(aT2 - aT1);
	423	//cout << "Errors: " << anIntegral.NbIterReached() << " " << anIntegral.AbsolutError() << " " << err1 << endl;
	424	aLocalTolReached(k) += anIntegral.AbsolutError() + err1;
	425	//cout << "--- Errors: " << anIntegral.NbIterReached() << " " << anIntegral.AbsolutError() << " " << err1 << endl;
7fd59977	426	}
	427
	428	aValue(k) += aLocalValue(k);
	429	aTolReached(k) += aLocalTolReached(k);
	430	}
	431
	432	// If the pointer to the domain is NULL, there is only one curve to treat:
	433	// U isoline with the UMax parameter.
	434	if (thePtrDomain == NULL)
	435	isMore = Standard_False;
	436	else {
	437	aPDomain->Next();
	438	isMore = aPDomain->More();
	439	}
	440	}
	441
	442	// Get volume value.
	443	dim = aValue(1);
	444	myErrorReached = aTolReached(1);
	445	myAbsolutError = myErrorReached;
	446	Standard_Real anAbsDim = Abs(dim);
	447	Standard_Real aVolTol = Epsilon(myAbsolutError);
	448	if(anAbsDim >= aVolTol) myErrorReached /= anAbsDim;
	449
	450	if(theCGFlag \|\| theIFlag) {
	451	// Compute values of center of mass.
	452	if(anAbsDim >= aVolTol) {
	453	if (IsByPoint) {
424cd6bb	454	aValue(2) = aCoeffs[0] + aValue(2)/dim;
	455	aValue(3) = aCoeffs[1] + aValue(3)/dim;
	456	aValue(4) = aCoeffs[2] + aValue(4)/dim;
7fd59977	457	} else {
424cd6bb	458	aValue(2) /= dim;
	459	aValue(3) /= dim;
	460	aValue(4) /= dim;
7fd59977	461	}
	462	} else {
	463	aValue(2) = 0.;
	464	aValue(3) = 0.;
	465	aValue(4) = 0.;
	466	dim = 0.;
	467	}
	468	g.SetCoord(aValue(2), aValue(3), aValue(4));
	469	}
	470
	471	if(theIFlag) {
	472	// Fill the matrix of inertia.
	473	inertia.SetCols (gp_XYZ (aValue(5), aValue(8), aValue(9)),
424cd6bb	474	gp_XYZ (aValue(8), aValue(6), aValue(10)),
424cd6bb	475	gp_XYZ (aValue(9), aValue(10), aValue(7)));
7fd59977	476	}
	477	//return myErrorReached;
	478	return myAbsolutError;
	479	}
424cd6bb	480