/[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 2786 by caltinay, Thu Nov 26 04:42:29 2009 UTC revision 2787 by jfenwick, Fri Nov 27 05:03:09 2009 UTC
# Line 3271  Data::interpolateFromTable1D(const Wrapp Line 3271  Data::interpolateFromTable1D(const Wrapp
3271          expand();          expand();
3272      }      }
3273      Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);      Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);
3274        int numpts=getNumDataPoints();
3275        int twidth=table.getShape()[0]-1;  
3276        bool haserror=false;
3277        const DataVector* adat=0;
3278        DataVector* rdat=0;
3279      try      try
3280      {      {
3281          int numpts=getNumDataPoints();          adat=&(getReady()->getVectorRO());
3282          const DataVector& adat=getReady()->getVectorRO();          rdat=&(res.getReady()->getVectorRW());
3283          DataVector& rdat=res.getReady()->getVectorRW();      } catch (DataException d)
3284          int twidth=table.getShape()[0]-1;      {
3285          bool haserror=false;          haserror=true;
3286          int l=0;          error=3;
3287          #pragma omp parallel for private(l) schedule(static)      }
3288          for (l=0;l<numpts; ++l)      if (!haserror)
3289          {      {
3290            int l=0;
3291            #pragma omp parallel for private(l) schedule(static)
3292            for (l=0;l<numpts; ++l)
3293            {
3294              int lerror=0;
3295            #pragma omp flush(haserror)       // In case haserror was in register            #pragma omp flush(haserror)       // In case haserror was in register
3296            if (!haserror)                    if (!haserror)        
3297            {            {
3298             int lerror=0;              double a=(*adat)[l];
            try  
            {  
              do         // so we can use break  
              {  
             double a=adat[l];  
3299              int x=static_cast<int>(((a-Amin)/Astep));              int x=static_cast<int>(((a-Amin)/Astep));
3300              if (check_boundaries) {              if (check_boundaries)
3301                {
3302                  if ((a<Amin) || (x<0))                  if ((a<Amin) || (x<0))
3303                  {                  {
3304                      lerror=1;                      lerror=1;
3305                      break;                    }
3306                  }                  else if (a>Amin+Astep*twidth)
                 if (a>Amin+Astep*twidth)  
3307                  {                  {
3308                          lerror=4;                          lerror=4;
                     break;  
3309                  }                  }
             }  
             if (x<0) x=0;  
             if (x>twidth) x=twidth;  
   
             if (x==twidth)          // value is on the far end of the table  
             {  
                 double e=table.getElt(x);  
                 if (e>undef)  
                 {  
                     lerror=2;  
                     break;  
                 }  
                 rdat[l]=e;  
3310              }              }
3311              else        // x and y are in bounds              if (!lerror)
3312              {              {
3313                  double e=table.getElt(x);                  if (x<0) x=0;
3314                  double w=table.getElt(x+1);                  if (x>twidth) x=twidth;
3315                  if ((e>undef) || (w>undef))                  try{
3316                  {                  if (x==twidth)          // value is on the far end of the table
3317                    {
3318                        double e=table.getElt(x);
3319                        if (e>undef)
3320                        {
3321                      lerror=2;                      lerror=2;
3322                      break;                      }
3323                  }                      else
3324          // map x*Astep <= a << (x+1)*Astep to [-1,1]                      {
3325                  double la = 2.0*(a-Amin-(x*Astep))/Astep-1;                      (*rdat)[l]=e;
3326                  rdat[l]=((1-la)*e + (1+la)*w)/2;                      }
3327              }                    }
3328                } while (false);                  else        // x and y are in bounds
3329              } catch (DataException d)                  {
3330              {                      double e=table.getElt(x);
3331                  lerror=3;                      double w=table.getElt(x+1);
3332              }                      if ((e>undef) || (w>undef))
3333                        {
3334                        lerror=2;
3335                        }
3336                        else
3337                        {
3338                            // map x*Astep <= a << (x+1)*Astep to [-1,1]
3339                        double la = 2.0*(a-Amin-(x*Astep))/Astep-1;
3340                          (*rdat)[l]=((1-la)*e + (1+la)*w)/2;
3341                        }
3342                    }
3343                    }
3344                    catch (DataException d)
3345                    {
3346                    lerror=3;
3347                    }  
3348                }   // if !lerror
3349              if (lerror!=0)              if (lerror!=0)
3350              {              {
3351              #pragma omp critical    // Doco says there is a flush associated with critical              #pragma omp critical    // Doco says there is a flush associated with critical
# Line 3344  Data::interpolateFromTable1D(const Wrapp Line 3355  Data::interpolateFromTable1D(const Wrapp
3355              }              }
3356              }              }
3357            } // if (!error)            } // if (!error)
3358          }   // parallelised for          }   // parallelised for
3359      } catch (DataException d)      } // if !haserror
     {  
         error=3;        // this is outside the parallel region so assign directly  
     }  
3360  #ifdef PASO_MPI  #ifdef PASO_MPI
3361      int rerror=0;      int rerror=0;
3362      MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );      MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );
# Line 3401  Data::interpolateFromTable2D(const Wrapp Line 3409  Data::interpolateFromTable2D(const Wrapp
3409      {      {
3410      B.expand();      B.expand();
3411      }      }
3412    
3413      Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);      Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);
3414    
3415        int numpts=getNumDataPoints();
3416        const DataVector* adat=0;
3417        const DataVector* bdat=0;
3418        DataVector* rdat=0;
3419        const DataTypes::ShapeType& ts=table.getShape();
3420      try      try
3421      {      {
3422      int numpts=getNumDataPoints();      adat=&(getReady()->getVectorRO());
3423      const DataVector& adat=getReady()->getVectorRO();      bdat=&(B.getReady()->getVectorRO());
3424      const DataVector& bdat=B.getReady()->getVectorRO();      rdat=&(res.getReady()->getVectorRW());
3425      DataVector& rdat=res.getReady()->getVectorRW();      }
3426      const DataTypes::ShapeType& ts=table.getShape();      catch (DataException e)
3427        {
3428        error=3;
3429        }
3430        if (!error)
3431        {
3432      int twx=ts[0]-1;    // table width x      int twx=ts[0]-1;    // table width x
3433      int twy=ts[1]-1;    // table width y      int twy=ts[1]-1;    // table width y
3434      bool haserror=false;      bool haserror=false;
3435      int l=0;      int l=0;
3436      #pragma omp parallel for private(l) schedule(static)      #pragma omp parallel for private(l) shared(res,rdat, adat, bdat, ts) schedule(static)
3437      for (l=0; l<numpts; ++l)      for (l=0; l<numpts; ++l)
3438      {      {
3439         #pragma omp flush(haserror)      // In case haserror was in register         #pragma omp flush(haserror)      // In case haserror was in register
3440         if (!haserror)               if (!haserror)      
3441         {         {
3442           int lerror=0;          int lerror=0;
3443           try          double a=(*adat)[l];
3444           {          double b=(*bdat)[l];
            do  
            {  
         double a=adat[l];  
         double b=bdat[l];  
3445          int x=static_cast<int>(((a-Amin)/Astep));          int x=static_cast<int>(((a-Amin)/Astep));
3446          int y=static_cast<int>(((b-Bmin)/Bstep));          int y=static_cast<int>(((b-Bmin)/Bstep));
3447                  if (check_boundaries) {                  if (check_boundaries)
3448            {
3449                  if ( (a<Amin) || (b<Bmin) || (x<0) || (y<0) )                  if ( (a<Amin) || (b<Bmin) || (x<0) || (y<0) )
3450                  {                  {
3451                  #pragma omp critical                  lerror=1;
                 {  
                     lerror=1;  
                 }  
                 break;    
3452                  }                  }
3453                  if ( (a>Amin+Astep*twx) || (b>Bmin+Bstep*twy) )                  else if ( (a>Amin+Astep*twx) || (b>Bmin+Bstep*twy) )
3454                  {                  {
3455                  #pragma omp critical                  lerror=4;
                 {  
                     lerror=4;  
                 }  
                 break;  
3456                  }                  }
3457                  }                  }
3458                  if (x<0) x=0;          if (lerror==0)
3459                  if (y<0) y=0;          {
3460                  if (x>twx) x=twx;                      if (x<0) x=0;
3461                  if (y>twx) y=twy;                      if (y<0) y=0;
3462                        if (x>twx) x=twx;
3463                  if (x == twx ) {              if (y>twx) y=twy;
3464                       if (y == twy ) {              try
3465                   double sw=table.getElt(x,y);              {
3466                   if ((sw>undef))                  if (x == twx )
3467                   {              {
3468                   #pragma omp critical                              if (y == twy )
                  {  
                      lerror=2;  
                  }  
                  break;  
                  }  
                  rdat[l]=sw;  
   
                      } else {  
                  double sw=table.getElt(x,y);  
                  double nw=table.getElt(x,y+1);  
                  if ((sw>undef) || (nw>undef))  
                  {  
                  #pragma omp critical  
                  {  
                     lerror=2;  
                  }  
                  break;  
                  }  
                  double lb = 2.0*(b-Bmin-(y*Bstep))/Bstep-1;  
                  rdat[l]=((1-lb)*sw + (1+lb)*nw )/2.;  
   
                      }  
                 } else {  
                      if (y == twy ) {  
                  double sw=table.getElt(x,y);  
                  double se=table.getElt(x+1,y);  
                  if ((sw>undef) || (se>undef) )  
                  {  
                  #pragma omp critical  
                  {  
                     lerror=2;  
                  }  
                  break;  
                  }  
                  double la = 2.0*(a-Amin-(x*Astep))/Astep-1;  
                  rdat[l]=((1-la)*sw + (1+la)*se )/2;  
   
                      } else {  
                  double sw=table.getElt(x,y);  
                  double nw=table.getElt(x,y+1);  
                  double se=table.getElt(x+1,y);  
                  double ne=table.getElt(x+1,y+1);  
                  if ((sw>undef) || (nw>undef) || (se>undef) || (ne>undef))  
                  {  
                 #pragma omp critical  
3469                  {                  {
3470                       lerror=2;                  double sw=table.getElt(x,y);
3471                  }                  if ((sw>undef))
3472                  break;                  {
3473                   }                      lerror=2;
3474                   // map x*Astep <= a << (x+1)*Astep to [-1,1]                  }
3475                   // same with b                  else
3476                   double la = 2.0*(a-Amin-(x*Astep))/Astep-1;                  {
3477                   double lb = 2.0*(b-Bmin-(y*Bstep))/Bstep-1;                      (*rdat)[l]=sw;
3478                   rdat[l]=((1-la)*(1-lb)*sw + (1-la)*(1+lb)*nw +                  }
3479                        (1+la)*(1-lb)*se + (1+la)*(1+lb)*ne)/4;                              }
3480                       }                  else
3481                  }                  {
3482             } while (false);                  double sw=table.getElt(x,y);
3483            } catch (DataException d)                  double nw=table.getElt(x,y+1);
3484            {                  if ((sw>undef) || (nw>undef))
3485          lerror=3;                  {
3486            }                      lerror=2;
3487            if (lerror!=0)                  }
3488            {                  else
3489                    {
3490                        double lb = 2.0*(b-Bmin-(y*Bstep))/Bstep-1;
3491                        (*rdat)[l]=((1-lb)*sw + (1+lb)*nw )/2.;
3492                    }
3493                                }
3494                            }
3495                    else    // x != twx
3496                    {
3497                    if (y == twy )
3498                    {
3499                        double sw=table.getElt(x,y);
3500                        double se=table.getElt(x+1,y);
3501                        if ((sw>undef) || (se>undef) )
3502                        {
3503                        lerror=2;
3504                        }
3505                        else
3506                        {
3507                        double la = 2.0*(a-Amin-(x*Astep))/Astep-1;
3508                        (*rdat)[l]=((1-la)*sw + (1+la)*se )/2;
3509                        }
3510                    }
3511                    else
3512                    {
3513                            double sw=table.getElt(x,y);
3514                                double nw=table.getElt(x,y+1);
3515                        double se=table.getElt(x+1,y);
3516                        double ne=table.getElt(x+1,y+1);
3517                        if ((sw>undef) || (nw>undef) || (se>undef) || (ne>undef))
3518                        {
3519                        lerror=2;
3520                                }
3521                            else
3522                            {
3523                        // map x*Astep <= a << (x+1)*Astep to [-1,1]
3524                        // same with b
3525                        double la = 2.0*(a-Amin-(x*Astep))/Astep-1;
3526                        double lb = 2.0*(b-Bmin-(y*Bstep))/Bstep-1;
3527                        (*rdat)[l]=((1-la)*(1-lb)*sw + (1-la)*(1+lb)*nw +
3528                            (1+la)*(1-lb)*se + (1+la)*(1+lb)*ne)/4;
3529                            }
3530                    }  //else  (y != twy)
3531                    }   // else x != twx
3532                }
3533                catch (DataException d)
3534                {
3535                lerror=3;
3536                }
3537            }
3538            if (lerror!=0)
3539            {
3540              #pragma omp critical    // Doco says there is a flush associated with critical              #pragma omp critical    // Doco says there is a flush associated with critical
3541              {              {
3542                  haserror=true;  // We only care that one error is recorded. We don't care which                  haserror=true;  // We only care that one error is recorded. We don't care which
3543                  error=lerror;   // one                  error=lerror;   // one
3544              }                    }      
3545            }          }
3546          }  // if (!error)          }  // if (!haserror)
3547      }   // parallel for      }   // parallel for
3548      } catch (DataException d)       } // !error
     {  
         error=3;  
     }  
3549  #ifdef PASO_MPI  #ifdef PASO_MPI
3550      int rerror=0;       int rerror=0;
3551      MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );       MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );
3552      error=rerror;       error=rerror;
3553  #endif  #endif
3554      if (error)       if (error)
3555      {       {
3556      switch (error)      switch (error)
3557      {      {
3558              case 1: throw DataException("Value below lower table range.");              case 1: throw DataException("Value below lower table range.");
# Line 3549  Data::interpolateFromTable2D(const Wrapp Line 3561  Data::interpolateFromTable2D(const Wrapp
3561                      default:                      default:
3562                        throw DataException("Unknown error in interpolation");                                throw DataException("Unknown error in interpolation");        
3563      }      }
3564      }       }
3565      return res;       return res;
3566  }  }
3567    
3568    

Legend:
Removed from v.2786  
changed lines
  Added in v.2787

  ViewVC Help
Powered by ViewVC 1.1.26