/[escript]/branches/symbolic_from_3470/escript/src/Data.cpp
ViewVC logotype

Diff of /branches/symbolic_from_3470/escript/src/Data.cpp

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

revision 3788 by caltinay, Wed Sep 21 04:48:06 2011 UTC revision 3789 by caltinay, Tue Jan 31 04:55:05 2012 UTC
# Line 3428  Data::interpolateFromTable1D(const Wrapp Line 3428  Data::interpolateFromTable1D(const Wrapp
3428      return res;      return res;
3429  }  }
3430    
           
3431  Data  Data
3432  Data::interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,  Data::interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
3433                         double undef, Data& B, double Bmin, double Bstep, bool check_boundaries)                         double undef, Data& B, double Bmin, double Bstep, bool check_boundaries)
# Line 3445  Data::interpolateFromTable2D(const Wrapp Line 3444  Data::interpolateFromTable2D(const Wrapp
3444      }      }
3445      if ((Astep<=0) || (Bstep<=0))      if ((Astep<=0) || (Bstep<=0))
3446      {      {
3447      throw DataException("Astep and Bstep must be postive");      throw DataException("All step components must be strictly positive.");
3448      }      }
3449      if (getFunctionSpace()!=B.getFunctionSpace())      if (getFunctionSpace()!=B.getFunctionSpace())
3450      {      {
3451      Data n=B.interpolate(getFunctionSpace());      Data n=B.interpolate(getFunctionSpace());
3452      return interpolateFromTable2D(table, Amin, Astep, undef,      return interpolateFromTable2D(table, Amin, Astep, undef,
3453          n , Bmin, Bstep, check_boundaries);          n , Bmin, Bstep, check_boundaries);      
3454      }      }
3455    
3456      if (!isExpanded())      if (!isExpanded())
3457      {      {
3458      expand();      expand();
# Line 3481  Data::interpolateFromTable2D(const Wrapp Line 3481  Data::interpolateFromTable2D(const Wrapp
3481      }      }
3482      if (!error)      if (!error)
3483      {      {
3484      int twx=ts[0]-1;    // table width x      int twx=ts[1]-1;    // table width x
3485      int twy=ts[1]-1;    // table width y      int twy=ts[0]-1;    // table width y
3486    
3487      bool haserror=false;      bool haserror=false;
3488      int l=0;      int l=0;
3489      #pragma omp parallel for private(l) shared(res,rdat, adat, bdat) schedule(static)      #pragma omp parallel for private(l) shared(res,rdat, adat, bdat) schedule(static)
# Line 3498  Data::interpolateFromTable2D(const Wrapp Line 3499  Data::interpolateFromTable2D(const Wrapp
3499          int y=static_cast<int>(((b-Bmin)/Bstep));          int y=static_cast<int>(((b-Bmin)/Bstep));
3500                  if (check_boundaries)                  if (check_boundaries)
3501          {          {
3502                  if ( (a<Amin) || (b<Bmin) || (x<0) || (y<0) )                  if ( (a<Amin) || (b<Bmin) || (x<0) || (y<0))
3503                  {                  {
3504                  lerror=1;                  lerror=1;
3505                  }                  }
3506                  else if ( (a>Amin+Astep*twx) || (b>Bmin+Bstep*twy) )                  else if ( (a>Amin+Astep*twx) || (b>Bmin+Bstep*twy))
3507                  {                  {
3508                  lerror=4;                  lerror=4;
3509                  }                  }
# Line 3511  Data::interpolateFromTable2D(const Wrapp Line 3512  Data::interpolateFromTable2D(const Wrapp
3512          {          {
3513                      if (x<0) x=0;                      if (x<0) x=0;
3514                      if (y<0) y=0;                      if (y<0) y=0;
3515    
3516                      if (x>twx) x=twx;                      if (x>twx) x=twx;
3517              if (y>twy) y=twy;              if (y>twy) y=twy;
3518              try              try
3519              {              {
3520                  if (x == twx )              int nx=x+1;
3521              {              int ny=y+1;
3522                              if (y == twy )  
3523                  {              double la=0;    /* map position of a between x and nx to [-1,1] */
3524                  double sw=table.getElt(x,y);              double lb=0;
3525                  if ((sw>undef))              double weight=4;
3526                  {  
3527                      lerror=2;              // now we work out which terms we should be considering
3528                  }              bool usex=(x!=twx);
3529                  else              bool usey=(y!=twy);
3530                  {  
3531                      (*rdat)[l]=sw;              la = 2.0*(a-Amin-(x*Astep))/Astep-1;
3532                  }              lb = 2.0*(b-Bmin-(y*Bstep))/Bstep-1;
3533                              }  
3534                  else  
3535                  {              
3536                  double sw=table.getElt(x,y);              double sw=table.getElt(y,x);
3537                  double nw=table.getElt(x,y+1);              double nw=usey?table.getElt(ny,x):0;    // 0 because if !usey ny does not actually exist
3538                  if ((sw>undef) || (nw>undef))              double se=usex?table.getElt(y,nx):0;
3539                  {              double ne=(usex&&usey)?table.getElt(ny,nx):0;
3540                      lerror=2;  
3541                  }  // cout << a << "," << b << " -> " << x << "," << y << "   " <<  sw <<  "," <<
3542                  else  // nw <<  "," <<  se <<  "," <<  ne <<  "\n";          
3543                  {  
3544                      double lb = 2.0*(b-Bmin-(y*Bstep))/Bstep-1;              double ans=(1-la)*(1-lb)*sw +
3545                      (*rdat)[l]=((1-lb)*sw + (1+lb)*nw )/2.;                     (1-la)*(1+lb)*nw +
3546                  }                     (1+la)*(1-lb)*se +
3547                              }                     (1+la)*(1+lb)*ne;
3548                          }              ans/=weight;
3549                  else    // x != twx              (*rdat)[l]=ans;
3550                  {              // this code does not check to see if any of the points used in the interpolation are undef
3551                  if (y == twy )              if (ans>undef)
3552                  {              {
3553                      double sw=table.getElt(x,y);                  lerror=2;
3554                      double se=table.getElt(x+1,y);              }
                     if ((sw>undef) || (se>undef) )  
                     {  
                     lerror=2;  
                     }  
                     else  
                     {  
                     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))  
                     {  
                     lerror=2;  
                             }  
                         else  
                         {  
                     // map x*Astep <= a << (x+1)*Astep to [-1,1]  
                     // same with b  
                     double la = 2.0*(a-Amin-(x*Astep))/Astep-1;  
                     double lb = 2.0*(b-Bmin-(y*Bstep))/Bstep-1;  
                     (*rdat)[l]=((1-la)*(1-lb)*sw + (1-la)*(1+lb)*nw +  
                         (1+la)*(1-lb)*se + (1+la)*(1+lb)*ne)/4;  
                         }  
                 }  //else  (y != twy)  
                 }   // else x != twx  
3555              }              }
3556              catch (DataException d)              catch (DataException d)
3557              {              {
# Line 3591  Data::interpolateFromTable2D(const Wrapp Line 3562  Data::interpolateFromTable2D(const Wrapp
3562          {          {
3563              #pragma omp critical    // Doco says there is a flush associated with critical              #pragma omp critical    // Doco says there is a flush associated with critical
3564              {              {
3565                  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
3566                  error=lerror;   // one                  error=lerror;   // one
3567              }                    }      
3568          }          }
# Line 3617  Data::interpolateFromTable2D(const Wrapp Line 3588  Data::interpolateFromTable2D(const Wrapp
3588       return res;       return res;
3589  }  }
3590    
3591    
3592  Data  Data
3593  Data::interpolateFromTable3D(const WrappedArray& table, double Amin, double Astep,  Data::interpolateFromTable3D(const WrappedArray& table, double Amin, double Astep,
3594                         double undef, Data& B, double Bmin, double Bstep, Data& C,                         double undef, Data& B, double Bmin, double Bstep, Data& C,
# Line 4203  escript::condEval(escript::Data& mask, e Line 4175  escript::condEval(escript::Data& mask, e
4175    
4176      // now we copy the tags from the mask - if the mask does not have it then it doesn't appear      // now we copy the tags from the mask - if the mask does not have it then it doesn't appear
4177      const DataTagged::DataMapType& maskLookup=mdat->getTagLookup();      const DataTagged::DataMapType& maskLookup=mdat->getTagLookup();
4178      DataTagged::DataMapType::const_iterator i;      DataTagged::DataMapType::const_iterator it;
4179      DataTagged::DataMapType::const_iterator thisLookupEnd=maskLookup.end();      DataTagged::DataMapType::const_iterator thisLookupEnd=maskLookup.end();
4180      for (i=maskLookup.begin();i!=thisLookupEnd;i++)      for (it=maskLookup.begin();it!=thisLookupEnd;it++)
4181      {      {
4182          if (mdat->getDataByTagRO(i->first,0)>0)          if (mdat->getDataByTagRO(it->first,0)>0)
4183          {          {
4184          rdat->addTaggedValue(i->first,trueval.getDataPointShape(), tdat->getVectorRO(), tdat->getOffsetForTag(i->first));          rdat->addTaggedValue(it->first,trueval.getDataPointShape(), tdat->getVectorRO(), tdat->getOffsetForTag(it->first));
4185          }          }
4186          else          else
4187          {          {
4188          rdat->addTaggedValue(i->first,falseval.getDataPointShape(), fdat->getVectorRO(), fdat->getOffsetForTag(i->first));          rdat->addTaggedValue(it->first,falseval.getDataPointShape(), fdat->getVectorRO(), fdat->getOffsetForTag(it->first));
4189          }          }
4190      }      }
4191    

Legend:
Removed from v.3788  
changed lines
  Added in v.3789

  ViewVC Help
Powered by ViewVC 1.1.26