/[escript]/trunk/escript/src/Data.cpp
ViewVC logotype

Diff of /trunk/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2671 by gross, Thu Sep 17 06:37:22 2009 UTC revision 2673 by jfenwick, Fri Sep 18 05:33:10 2009 UTC
# Line 3055  Data::interpolateFromTable1D(const Wrapp Line 3055  Data::interpolateFromTable1D(const Wrapp
3055      {      {
3056          throw DataException("Table for 1D interpolation must be 1D");          throw DataException("Table for 1D interpolation must be 1D");
3057      }      }
3058        if (Astep<=0)
3059        {
3060            throw DataException("Astep must be positive");
3061        }
3062      if (!isExpanded())      if (!isExpanded())
3063      {      {
3064          expand();          expand();
# Line 3067  Data::interpolateFromTable1D(const Wrapp Line 3071  Data::interpolateFromTable1D(const Wrapp
3071              int numpts=getNumDataPoints();              int numpts=getNumDataPoints();
3072              const DataVector& adat=getReady()->getVectorRO();              const DataVector& adat=getReady()->getVectorRO();
3073              DataVector& rdat=res.getReady()->getVectorRW();              DataVector& rdat=res.getReady()->getVectorRW();
3074              int twidth=table.getShape()[0];              int twidth=table.getShape()[0]-1;
3075              for (int l=0; l<numpts; ++l)              for (int l=0; l<numpts; ++l)
3076              {              {
3077                  double a=adat[l];                  double a=adat[l];
3078                  int x=static_cast<int>(((a-Amin)/Astep));                  int x=static_cast<int>(((a-Amin)/Astep));
3079                                  if (check_boundaries) {                                  if (check_boundaries) {
3080                      if ( (x<0) || (a>Amin) )                      if ((a<Amin) || (x<0))
3081                      {                      {
3082                          error=1;                          error=1;
3083                          break;                            break;  
3084                      }                      }
3085                      if ( (x>=twidth) ||  (a>Amin+Astep*twidth))                      if (a>Amin+Astep*twidth)
3086                      {                      {
3087                          error=4;                          error=4;
3088                          break;                          break;
3089                      }                      }
3090                                  }                                  }
3091                                  if (x<0) x=0;                                  if (x<0) x=0;
3092                                  if (x>=twidth) x=twidth-1;                                  if (x>twidth) x=twidth;
3093    
3094                  if (x==(twidth-1))                  if (x==twidth)          // value is on the far end of the table
3095                  {                  {
3096                      double e=table.getElt(x);                      double e=table.getElt(x);
3097                      if (e>undef)                      if (e>undef)
# Line 3150  Data::interpolateFromTable2D(const Wrapp Line 3154  Data::interpolateFromTable2D(const Wrapp
3154      {      {
3155      throw DataException("Table for 2D interpolation must be 2D");      throw DataException("Table for 2D interpolation must be 2D");
3156      }      }
3157        if ((Astep<=0) || (Bstep<=0))
3158        {
3159        throw DataException("Astep and Bstep must be postive");
3160        }
3161      if (getFunctionSpace()!=B.getFunctionSpace())      if (getFunctionSpace()!=B.getFunctionSpace())
3162      {      {
3163      Data n=B.interpolate(getFunctionSpace());      Data n=B.interpolate(getFunctionSpace());
# Line 3174  Data::interpolateFromTable2D(const Wrapp Line 3182  Data::interpolateFromTable2D(const Wrapp
3182          const DataVector& bdat=B.getReady()->getVectorRO();          const DataVector& bdat=B.getReady()->getVectorRO();
3183          DataVector& rdat=res.getReady()->getVectorRW();          DataVector& rdat=res.getReady()->getVectorRW();
3184          const DataTypes::ShapeType& ts=table.getShape();          const DataTypes::ShapeType& ts=table.getShape();
3185            int twx=ts[0]-1;    // table width x
3186            int twy=ts[1]-1;    // table width y
3187          for (int l=0; l<numpts; ++l)          for (int l=0; l<numpts; ++l)
3188          {          {
3189          double a=adat[l];          double a=adat[l];
# Line 3181  Data::interpolateFromTable2D(const Wrapp Line 3191  Data::interpolateFromTable2D(const Wrapp
3191          int x=static_cast<int>(((a-Amin)/Astep));          int x=static_cast<int>(((a-Amin)/Astep));
3192          int y=static_cast<int>(((b-Bmin)/Bstep));          int y=static_cast<int>(((b-Bmin)/Bstep));
3193                  if (check_boundaries) {                  if (check_boundaries) {
3194                  if ( (x<0) || (a>Amin) || (y<0) || (y>Amin) )                  if ( (a<Amin) || (b<Bmin) || (x<0) || (y<0) )
3195                  {                  {
3196                      error=1;                      error=1;
3197                      break;                        break;  
3198                  }                  }
3199                  if ( (x>=ts[0]) || (a>Amin+Astep*ts[0]) || (y>=ts[1]) || (b>Bmin+Bstep*ts[1]) )                  if ( (a>Amin+Astep*twx) || (b>Bmin+Bstep*twy) )
3200                  {                  {
3201                      error=4;                      error=4;
3202                      break;                      break;
# Line 3194  Data::interpolateFromTable2D(const Wrapp Line 3204  Data::interpolateFromTable2D(const Wrapp
3204                  }                  }
3205                  if (x<0) x=0;                  if (x<0) x=0;
3206                  if (y<0) y=0;                  if (y<0) y=0;
3207                  if (x>=ts[0]) x=ts[0]-1;                  if (x>twx) x=twx;
3208                  if (y>=ts[0]) y=ts[1]-1;                  if (y>twx) y=twy;
3209    
3210                  if (x == ts[0] - 1 ) {                  if (x == twx ) {
3211                       if (y == ts[1] - 1 ) {                       if (y == twy ) {
3212                   double sw=table.getElt(x,y);                   double sw=table.getElt(x,y);
3213                   if ((sw>undef))                   if ((sw>undef))
3214                   {                   {
# Line 3220  Data::interpolateFromTable2D(const Wrapp Line 3230  Data::interpolateFromTable2D(const Wrapp
3230    
3231                       }                       }
3232                  } else {                  } else {
3233                       if (y == ts[1] - 1 ) {                       if (y == twy ) {
3234                   double sw=table.getElt(x,y);                   double sw=table.getElt(x,y);
3235                   double se=table.getElt(x+1,y);                   double se=table.getElt(x+1,y);
3236                   if ((sw>undef) || (se>undef) )                   if ((sw>undef) || (se>undef) )
# Line 3248  Data::interpolateFromTable2D(const Wrapp Line 3258  Data::interpolateFromTable2D(const Wrapp
3258                   rdat[l]=((1-la)*(1-lb)*sw + (1-la)*(1+lb)*nw +                   rdat[l]=((1-la)*(1-lb)*sw + (1-la)*(1+lb)*nw +
3259                        (1+la)*(1-lb)*se + (1+la)*(1+lb)*ne)/4;                        (1+la)*(1-lb)*se + (1+la)*(1+lb)*ne)/4;
3260                       }                       }
   
3261                  }                  }
3262          }          }
3263      } catch (DataException d)      } catch (DataException d)

Legend:
Removed from v.2671  
changed lines
  Added in v.2673

  ViewVC Help
Powered by ViewVC 1.1.26