dfchaventmobility.cpp

Go to the documentation of this file.

#include "dfchaventmobility.h"


 DFChaventMobility::DFChaventMobility(double vw,double vo,double Srw,double Sro) 
:Function1D(1)
{
  m_M=vw/vo;

  MaxSo=1-Srw;
  MaxSw=1-Sro;
  MinSw=Srw;
  MinSo=Sro;


  //1_Srw2=1_Srw*1_Srw;
  //1_Sro=1-Sro;

  _c = m_M*MaxSo*MaxSo/(MaxSw*MaxSw);
  _d = 2*_c*(MaxSw-MinSw);
  _a1=_c+1;
  _a2=-2*(MinSw + _c*MaxSw);
  _a3=_c*MaxSw*MaxSw+MinSw*MinSw;

  m_inflectPoint = getInflectionDerivativePoint(m_M);
  m_inflectPointDeriv=(*this)(m_inflectPoint);
}




double DFChaventMobility::operator()(double x,unsigned cmp) const
{
   if (x >= MaxSw || x <= MinSw ) 
     return 0.0;
   register double quot=(_a1*x + _a2)*x + _a3;
   quot*=quot;
   return _d*(MaxSw-x)*(x-MinSw)/quot;
}


void DFChaventMobility::getMinMaxValues(double a, double b,double &min,double &max) const 
{
  
  assert(a<=b);
  if (a<0.0) a=0.0;
  if (b>1.0) b=1.0;
  if (b< m_inflectPoint)
  {
    min=(*this)(a,0);
    max=(*this)(b,0);
  }
  else if (a>m_inflectPoint)
  {
    max=(*this)(a,0);
    min=(*this)(b,0);
  }
  else
  {
    max=m_inflectPointDeriv;
    min=std::min( (*this)(a,0), (*this)(b,0));
  }
}


double DFChaventMobility::getMaxNorm(double a, double b) const    
{
  assert(a<=b);
  if (a<0.0) a=0.0;
  if (b>1.0) b=1.0;
  if (b< m_inflectPoint)
  {
    return (*this)(b,0);
  }
  else if (a>m_inflectPoint)
  {
    return (*this)(a,0);
  }
  else
  {
    return m_inflectPointDeriv;
  }

}

 void DFChaventMobility::updateInflectionPoint() 
{
  m_inflectPoint = getInflectionDerivativePoint(m_M);
  m_inflectPointDeriv=(*this)(m_inflectPoint);
}




 double DFChaventMobility::getInflectionDerivativePoint(double M) 
{
  static double sqrt3 = sqrt(3.0);
  

  double dd=MaxSw-MinSw;
  double dd2=dd*dd;
  double dd3=dd*dd2;

  double b = -sqrt(_c)*dd3/(4*(_c+1));
  double a = (1-_c)*dd3/(8*(_c+1));
  double r=pow(hypot(a,b),1.0/3.0);
  double theta=atan2(b,a)/3.0;
  a=r*cos(theta);
  b=r*sin(theta);
  double aux1=(-sqrt3*b-a)/2.0;
    
  return aux1*dd2/(4*(b*b + a*a)) + aux1 + (MaxSw+MinSw)/2.0;
  
  
  
}


 void DFChaventMobility::setParameters(const VecDouble &v) 
{
  assert(v.size() ==2);
  m_M=v(0)/v(1);
}


 DFChaventMobility::~DFChaventMobility() 
{

}

---