CO2INJECTION: probefunction1d.cpp Source File

00001 #include "probefunction1d.h"
00002 #include "numericmethods.h"
00003 
00004 
00005 
00006 
00007 ProbeFunction1D::ProbeFunction1D(Function1D &f,unsigned cmp,double a,double b,unsigned nPoints) 
00008   :m_f(&f),m_cmp(cmp)
00009 {
00010   assert(a<b);
00011 
00012   m_a = a;
00013   m_b = b;
00014   
00015   double dx = (b-a)/((double) nPoints);
00016   double xPrev = a;
00017   double fxPrev=f(a,m_cmp);
00018 
00019   double x=a+dx;
00020   double fx=f(x,m_cmp);
00021 
00022   FunctionBehavour currType;
00023 
00024   
00025   this->add(xPrev,fxPrev);
00026   if (fxPrev < fx)
00027     currType = INCREASING;
00028   else
00029     currType = DECREASING;
00030   this->setType(currType);
00031   
00032   for (unsigned i=1;i<nPoints;i++)
00033   {
00034     fxPrev = fx;
00035     xPrev = x;
00036     x+=dx;
00037     fx = f(x,m_cmp);
00038     
00039     if (currType == INCREASING && (fx < fxPrev))
00040     {
00041       this->add(xPrev,fxPrev);
00042       currType=DECREASING;
00043       this->setType(currType);
00044     }
00045     if (currType == DECREASING && (fx > fxPrev))
00046     {
00047       this->add(xPrev,fxPrev);
00048       currType=INCREASING;
00049       this->setType(currType);
00050     }
00051   }
00052 }
00053 
00054 void ProbeFunction1D::add(double x,double fx) 
00055 {
00056   m_vF.push_back(fx);
00057   m_vX.push_back(x);
00058 }
00059 
00060 void ProbeFunction1D::setType(FunctionBehavour beh) 
00061 {
00062   m_vType.push_back(beh);
00063 }
00064 
00065 
00072 void  ProbeFunction1D::getMinMaxValues(double c1,double c2,double &min,double &max) const
00073 {
00074   //See if the probe function already got the probe information.
00075   //Users should always call the sampleFunction method first,
00076   //Or use the non default constructor 
00077   assert(m_f);
00078   assert(c1<=c2);
00079   //  assert(c1 >= m_a && c2 <= m_b);
00080   NumericMethods::adjustBounds(c1,m_a,m_b);
00081   NumericMethods::adjustBounds(c2,m_a,m_b);
00082   int indexC1 = this->find(c1);
00083   int indexC2 = this->find(c2);
00084   Function1D &f = *m_f;
00085   if (indexC1 == indexC2)
00086   {
00087     if (this->getType(indexC1) == INCREASING)
00088     {
00089       min = f(c1,m_cmp);
00090       max = f(c2,m_cmp);
00091     }
00092     else 
00093     {
00094       min = f(c2,m_cmp);
00095       max = f(c1,m_cmp);
00096     }
00097   }
00098   else
00099   {
00100     assert(indexC1 < indexC2 && (unsigned) indexC2 < m_vF.size());
00101     double fx;
00102     min = max = f(c1,m_cmp);
00103     for (int i= indexC1+1;i<=indexC2;i++)
00104     {
00105       fx = m_vF[i];
00106       if (fx > max)
00107         max = fx;
00108       else if (fx < min)
00109         min = fx;
00110     }
00111     fx = f(c2,m_cmp);
00112     if (fx > max)
00113       max = fx;
00114     else if (fx < min)
00115       min = fx;
00116   }
00117 }
00118 
00119 
00120 void  ProbeFunction1D::getMaxValues(double c1,double c2,double &max) 
00121 {
00122   assert(c1<=c2);
00123   assert(c1 >= m_a && c2 <= m_b);
00124     //  NumericMethods::adjustBounds(c1,m_a,m_b);
00125     //NumericMethods::adjustBounds(c2,m_a,m_b);
00126 
00127   Function1D &f = *m_f;
00128   int indexC1 = this->find(c1);
00129   int indexC2 = this->find(c2);
00130 
00131   if (indexC1 == indexC2)
00132   {
00133     if (this->getType(indexC1) == INCREASING)
00134     {
00135       max = f(c2,m_cmp);
00136     }
00137     else 
00138     {
00139       max = f(c1,m_cmp);
00140     }
00141   }
00142   else
00143   {
00144     assert(indexC1 < indexC2 && (unsigned) indexC2 < m_vF.size());
00145     double fx;
00146     max = f(c1,m_cmp);
00147     for (int i= indexC1+1;i<=indexC2;i++)
00148     {
00149       fx = m_vF[i];
00150       if (fx > max)
00151         max = fx;
00152     }
00153     fx = f(c2,m_cmp);
00154     if (fx > max)
00155       max = fx;
00156   }
00157 }
00158 
00159 
00160 
00167 unsigned ProbeFunction1D::find(double c) const
00168 {
00169   assert(!m_vX.empty());
00170   if (!((c>=m_a) && (c<=m_b)))
00171        throw new Exception("ProbeFunction1D::getMinMaxValues value %g not contained in the probe interval [%g,%g]",
00172                            c,m_a,m_b);
00173   unsigned i;
00174   for(i=0;i<m_vX.size();i++)
00175   {
00176     if (m_vX[i] > c)
00177       return i-1;
00178   }
00179   return i-1;
00180 }
00181 
00182 
00187 void ProbeFunction1D::print() 
00188 {
00189   for (unsigned i=0;i<m_vX.size();i++)
00190   {
00191     printf("Point: %g,%g ",m_vX[i],m_vF[i]);
00192   if (m_vType[i] == INCREASING)
00193     printf("-- INCREASING --\n");
00194   if (m_vType[i] == DECREASING)
00195     printf("-- DECREASING --\n");
00196   }
00197   printf("Point: %g,%g\n",m_b,(*m_f)(m_b,m_cmp));
00198   
00199 }
00200 
00201 
00202 void ProbeFunction1D::sampleFunction(Function1D &f,unsigned cmp,double a,double b,unsigned nPoints) 
00203 {
00204   m_f=&f;
00205   m_cmp=cmp;
00206 
00207   assert(a<b);
00208 
00209   m_a = a;
00210   m_b = b;
00211   
00212   double dx = (b-a)/((double) nPoints);
00213   double xPrev = a;
00214   double fxPrev=f(a,m_cmp);
00215 
00216   double x=a+dx;
00217   double fx=f(x,m_cmp);
00218 
00219   FunctionBehavour currType;
00220 
00221   
00222   this->add(xPrev,fxPrev);
00223   if (fxPrev < fx)
00224     currType = INCREASING;
00225   else
00226     currType = DECREASING;
00227   this->setType(currType);
00228   
00229   for (unsigned i=1;i<nPoints;i++)
00230   {
00231     fxPrev = fx;
00232     xPrev = x;
00233     x+=dx;
00234     fx = f(x,m_cmp);
00235     
00236     if (currType == INCREASING && (fx < fxPrev))
00237     {
00238       this->add(xPrev,fxPrev);
00239       currType=DECREASING;
00240       this->setType(currType);
00241     }
00242     if (currType == DECREASING && (fx > fxPrev))
00243     {
00244       this->add(xPrev,fxPrev);
00245       currType=INCREASING;
00246       this->setType(currType);
00247     }
00248   }
00249 
00250 }