/[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 2653 by jfenwick, Tue Sep 8 04:26:30 2009 UTC revision 2668 by gross, Thu Sep 17 04:04:09 2009 UTC
# Line 3027  Data::freeSampleBuffer(BufferGroup* buff Line 3027  Data::freeSampleBuffer(BufferGroup* buff
3027    
3028  Data  Data
3029  Data::interpolateFromTable2DP(boost::python::object table, double Amin, double Astep,  Data::interpolateFromTable2DP(boost::python::object table, double Amin, double Astep,
3030          Data& B, double Bmin, double Bstep, double undef)          Data& B, double Bmin, double Bstep, double undef, bool check_boundaries)
3031  {  {
3032      WrappedArray t(table);      WrappedArray t(table);
3033      return interpolateFromTable2D(t, Amin, Astep, undef, B, Bmin, Bstep);      return interpolateFromTable2D(t, Amin, Astep, undef, B, Bmin, Bstep,check_boundaries);
3034  }  }
3035    
3036  Data  Data
3037  Data::interpolateFromTable1DP(boost::python::object table, double Amin, double Astep,  Data::interpolateFromTable1DP(boost::python::object table, double Amin, double Astep,
3038                            double undef)                            double undef,bool check_boundaries)
3039  {  {
3040      WrappedArray t(table);      WrappedArray t(table);
3041      return interpolateFromTable1D(t, Amin, Astep, undef);      return interpolateFromTable1D(t, Amin, Astep, undef, check_boundaries);
3042  }  }
3043    
3044    
3045  Data  Data
3046  Data::interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,  Data::interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,
3047                           double undef)                   double undef, bool check_boundaries)
3048  {  {
3049      int error=0;      int error=0;
3050      if ((getDataPointRank()!=0))      if ((getDataPointRank()!=0))
# Line 3071  Data::interpolateFromTable1D(const Wrapp Line 3071  Data::interpolateFromTable1D(const Wrapp
3071              for (int l=0; l<numpts; ++l)              for (int l=0; l<numpts; ++l)
3072              {              {
3073                  double a=adat[l];                  double a=adat[l];
3074                  int x=static_cast<int>((a-Amin)/Astep);                  int x=static_cast<int>(((a-Amin)/Astep));
3075                  if (a<Amin)                                  if (check_boundaries) {
3076                  {                      if ( (x<0) || (a>Amin) )
3077                      error=1;                      {
3078                      break;                            error=1;
3079                  }                          break;  
3080                  if (x>=(twidth-1))                      }
3081                        if ( (x>=twidth) ||  (a>Amin+Astep*twidth))
3082                        {
3083                            error=4;
3084                            break;
3085                        }
3086                                    }
3087                                    if (x<0) x=0;
3088                                    if (x>=twidth) x=twidth-1;
3089    
3090                    if (x==(twidth-1))
3091                  {                  {
3092                      error=1;                      double e=table.getElt(x);
3093                      break;                      if (e>undef)
3094                        {
3095                            error=2;
3096                            break;
3097                        }
3098                        rdat[l]=e;
3099                  }                  }
3100                  else        // x and y are in bounds                  else        // x and y are in bounds
3101                  {                  {
# Line 3111  Data::interpolateFromTable1D(const Wrapp Line 3126  Data::interpolateFromTable1D(const Wrapp
3126      {      {
3127          switch (error)          switch (error)
3128          {          {
3129              case 1: throw DataException("Point out of bounds");              case 1: throw DataException("Value below lower table range.");
3130              case 2: throw DataException("Interpolated value too large");              case 2: throw DataException("Interpolated value too large");
3131                case 4: throw DataException("Value greater than upper table range.");
3132              default:              default:
3133                  throw DataException("Unknown error in interpolation");                        throw DataException("Unknown error in interpolation");      
3134          }          }
# Line 3123  Data::interpolateFromTable1D(const Wrapp Line 3139  Data::interpolateFromTable1D(const Wrapp
3139                    
3140  Data  Data
3141  Data::interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,  Data::interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
3142                         double undef, Data& B, double Bmin, double Bstep)                         double undef, Data& B, double Bmin, double Bstep, bool check_boundaries)
3143  {  {
3144      int error=0;      int error=0;
3145      if ((getDataPointRank()!=0) || (B.getDataPointRank()!=0))      if ((getDataPointRank()!=0) || (B.getDataPointRank()!=0))
# Line 3138  Data::interpolateFromTable2D(const Wrapp Line 3154  Data::interpolateFromTable2D(const Wrapp
3154      {      {
3155      Data n=B.interpolate(getFunctionSpace());      Data n=B.interpolate(getFunctionSpace());
3156      return interpolateFromTable2D(table, Amin, Astep, undef,      return interpolateFromTable2D(table, Amin, Astep, undef,
3157          n , Bmin, Bstep);          n , Bmin, Bstep, check_boundaries);
3158      }      }
3159      if (!isExpanded())      if (!isExpanded())
3160      {      {
# Line 3162  Data::interpolateFromTable2D(const Wrapp Line 3178  Data::interpolateFromTable2D(const Wrapp
3178          {          {
3179          double a=adat[l];          double a=adat[l];
3180          double b=bdat[l];          double b=bdat[l];
3181          int x=static_cast<int>((a-Amin)/Astep);          int x=static_cast<int>(((a-Amin)/Astep));
3182          int y=static_cast<int>((b-Bmin)/Bstep);          int y=static_cast<int>(((b-Bmin)/Bstep));
3183          if ((a<Amin) || (b<Bmin))                  if (check_boundaries) {
3184          {                  if ( (x<0) || (a>Amin) || (y<0) || (y>Amin) )
3185              error=1;                  {
3186              break;                        error=1;
3187          }                      break;  
3188          if ((x>=(ts[0]-1)) || (y>=(ts[1]-1)))                  }
3189          {                  if ( (x>=ts[0]) || (a>Amin+Astep*ts[0]) || (y>=ts[1]) || (b>Bmin+Bstep*ts[1]) )
3190              error=1;                  {
3191              break;                      error=4;
3192          }                      break;
3193          else        // x and y are in bounds                  }
3194          {                  }
3195              double sw=table.getElt(x,y);                  if (x<0) x=0;
3196              double nw=table.getElt(x,y+1);                  if (y<0) y=0;
3197              double se=table.getElt(x+1,y);                  if (x>=ts[0]) x=ts[0]-1;
3198              double ne=table.getElt(x+1,y+1);                  if (y>=ts[0]) y=ts[1]-1;
3199              if ((sw>undef) || (nw>undef) || (se>undef) || (ne>undef))  
3200              {                  if (x == ts[0] - 1 ) {
3201              error=2;                       if (y == ts[1] - 1 ) {
3202              break;                   double sw=table.getElt(x,y);
3203              }                   if ((sw>undef))
3204              // map x*Astep <= a << (x+1)*Astep to [-1,1]                   {
3205              // same with b                   error=2;
3206              double la = 2.0*(a-(x*Astep))/Astep-1;                   break;
3207              double lb = 2.0*(b-(y*Bstep))/Bstep-1;                   }
3208              rdat[l]=((1-la)*(1-lb)*sw + (1-la)*(1+lb)*nw +                   rdat[l]=sw;
3209                   (1+la)*(1-lb)*se + (1+la)*(1+lb)*ne)/4;  
3210          }                       } else {
3211                     double sw=table.getElt(x,y);
3212                     double nw=table.getElt(x,y+1);
3213                     if ((sw>undef) || (nw>undef))
3214                     {
3215                     error=2;
3216                     break;
3217                     }
3218                     double lb = 2.0*(b-(y*Bstep))/Bstep-1;
3219                     rdat[l]=((1-lb)*sw + (1+lb)*nw )/2.;
3220    
3221                         }
3222                    } else {
3223                         if (y == ts[1] - 1 ) {
3224                     double sw=table.getElt(x,y);
3225                     double se=table.getElt(x+1,y);
3226                     if ((sw>undef) || (se>undef) )
3227                     {
3228                     error=2;
3229                     break;
3230                     }
3231                     double la = 2.0*(a-(x*Astep))/Astep-1;
3232                     rdat[l]=((1-la)*sw + (1+la)*se )/2;
3233    
3234                         } else {
3235                     double sw=table.getElt(x,y);
3236                     double nw=table.getElt(x,y+1);
3237                     double se=table.getElt(x+1,y);
3238                     double ne=table.getElt(x+1,y+1);
3239                     if ((sw>undef) || (nw>undef) || (se>undef) || (ne>undef))
3240                     {
3241                     error=2;
3242                     break;
3243                     }
3244                     // map x*Astep <= a << (x+1)*Astep to [-1,1]
3245                     // same with b
3246                     double la = 2.0*(a-(x*Astep))/Astep-1;
3247                     double lb = 2.0*(b-(y*Bstep))/Bstep-1;
3248                     rdat[l]=((1-la)*(1-lb)*sw + (1-la)*(1+lb)*nw +
3249                          (1+la)*(1-lb)*se + (1+la)*(1+lb)*ne)/4;
3250    
3251                         }
3252    
3253                    }
3254          }          }
3255      } catch (DataException d)      } catch (DataException d)
3256      {      {
# Line 3208  Data::interpolateFromTable2D(const Wrapp Line 3267  Data::interpolateFromTable2D(const Wrapp
3267      {      {
3268      switch (error)      switch (error)
3269      {      {
3270      case 1: throw DataException("Point out of bounds");              case 1: throw DataException("Value below lower table range.");
3271      case 2: throw DataException("Interpolated value too large");              case 2: throw DataException("Interpolated value too large");
3272      default:              case 4: throw DataException("Value greater than upper table range.");
3273          throw DataException("Unknown error in interpolation");                            default:
3274                          throw DataException("Unknown error in interpolation");        
3275      }      }
3276      }      }
3277      return res;      return res;

Legend:
Removed from v.2653  
changed lines
  Added in v.2668

  ViewVC Help
Powered by ViewVC 1.1.26