OCCT

// 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_TFunction.hxx>
#include <gp_Pnt.hxx>
#include <math_KronrodSingleIntegration.hxx>
#include <Precision.hxx>
#include <NCollection_Array1.hxx>
#include <NCollection_HArray1.hxx>

//=================================================================================================

BRepGProp_TFunction::BRepGProp_TFunction(const BRepGProp_Face& theSurface,
                                         const gp_Pnt&         theVertex,
                                         const bool            IsByPoint,
                                         const double*         theCoeffs,
                                         const double          theUMin,
                                         const double          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)
{
}

//=================================================================================================

void BRepGProp_TFunction::Init()
{
  myTolReached = 0.;
  myErrReached = 0.;
  myAbsError   = 0.;
}

//=================================================================================================

bool BRepGProp_TFunction::Value(const double X, double& F)
{
  const double tolU = 1.e-9;

  gp_Pnt2d aP2d;
  gp_Vec2d aV2d;
  double   aUMax;

  mySurface.D12d(X, aP2d, aV2d);
  aUMax = aP2d.X();

  if (aUMax - myUMin < tolU)
  {
    F = 0.;
    return true;
  }

  occ::handle<NCollection_HArray1<double>> aUKnots = mySurface.GetUKnots(myUMin, aUMax);
  myUFunction.SetVParam(aP2d.Y());

  // Compute the integral from myUMin to aUMax of myUFunction.
  int    i;
  double aCoeff = aV2d.Y();
  // int aNbUIntervals = aUKnots->Length() - 1;
  // double    aTol          = myTolerance/aNbUIntervals;
  double 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 false;
  }

  double aAbsCoeff = std::abs(aCoeff);

  if (aAbsCoeff <= Precision::Angular())
  {
    // No need to compute the integral. The value will be equal to 0.
    F = 0.;
    return true;
  }
  // else if (aAbsCoeff > 10.*aTol)
  //   aTol /= aAbsCoeff;
  // else
  //   aTol = 0.1;

  int                           iU = aUKnots->Upper();
  int                           aNbPntsStart;
  int                           aNbMaxIter = 1000;
  math_KronrodSingleIntegration anIntegral;
  double                        aLocalErr = 0.;

  i = aUKnots->Lower();
  F = 0.;

  // Epmirical criterion
  aNbPntsStart = std::min(15, mySurface.UIntegrationOrder() / (aUKnots->Length() - 1) + 1);
  aNbPntsStart = std::max(5, aNbPntsStart);

  while (i < iU)
  {
    double aU1 = aUKnots->Value(i++);
    double aU2 = aUKnots->Value(i);

    if (aU2 - aU1 < tolU)
    {
      continue;
    }

    anIntegral.Perform(myUFunction, aU1, aU2, aNbPntsStart, aTol, aNbMaxIter);

    if (!anIntegral.IsDone())
    {
      return false;
    }

    F += anIntegral.Value();
    aLocalErr += anIntegral.AbsolutError();
    // std::cout << " TFunction : " << anIntegral.NbIterReached() << " " <<
    // anIntegral.AbsolutError() << std::endl;
  }

  F *= aCoeff;
  aLocalErr *= aAbsCoeff;
  myAbsError = std::max(myAbsError, aLocalErr);

  myTolReached += aLocalErr;

  if (std::abs(F) > Epsilon(1.))
  {
    aLocalErr /= std::abs(F);
  }

  myErrReached = std::max(myErrReached, aLocalErr);

  return true;
}

//=================================================================================================

int BRepGProp_TFunction::GetStateNumber()
{
  // myErrReached  = myTolReached;
  // myTolReached  = 0.;
  // myNbPntOuter += RealToInt(0.5*myNbPntOuter);

  // if (myNbPntOuter%2 == 0)
  // myNbPntOuter++;

  return 0;
}