/[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 2645 by jfenwick, Wed Sep 2 04:14:03 2009 UTC revision 2646 by jfenwick, Fri Sep 4 00:13:00 2009 UTC
# Line 3031  Data::freeSampleBuffer(BufferGroup* buff Line 3031  Data::freeSampleBuffer(BufferGroup* buff
3031    
3032    
3033  Data  Data
3034  Data::interpolateFromTable(boost::python::object table, double Amin, double Astep,  Data::interpolateFromTable2DP(boost::python::object table, double Amin, double Astep,
3035           double undef, Data& B, double Bmin, double Bstep)          Data& B, double Bmin, double Bstep, double undef)
3036  {  {
3037      WrappedArray t(table);      WrappedArray t(table);
3038      return interpolateFromTable2D(t, Amin, Astep, undef, B, Bmin, Bstep);      return interpolateFromTable2D(t, Amin, Astep, undef, B, Bmin, Bstep);
3039  }  }
3040    
3041    Data
3042    Data::interpolateFromTable1DP(boost::python::object table, double Amin, double Astep,
3043                              double undef)
3044    {
3045        WrappedArray t(table);
3046        return interpolateFromTable1D(t, Amin, Astep, undef);
3047    }
3048    
3049    
3050    Data
3051    Data::interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,
3052                             double undef)
3053    {
3054        int error=0;
3055        if ((getDataPointRank()!=0))
3056        {
3057            throw DataException("Input to 1D interpolation must be scalar");
3058        }
3059        if (table.getRank()!=1)
3060        {
3061            throw DataException("Table for 1D interpolation must be 1D");
3062        }
3063        if (!isExpanded())
3064        {
3065            expand();
3066        }
3067        Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);
3068        do                                   // to make breaks useful
3069        {
3070            try
3071            {
3072                int numpts=getNumDataPoints();
3073                const DataVector& adat=getReady()->getVectorRO();
3074                DataVector& rdat=res.getReady()->getVectorRW();
3075                int twidth=table.getShape()[0];
3076                for (int l=0; l<numpts; ++l)
3077                {
3078                    double a=adat[l];
3079                    int x=static_cast<int>((a-Amin)/Astep);
3080                    if (a<Amin)
3081                    {
3082                        error=1;
3083                        break;  
3084                    }
3085                    if (x>=(twidth-1))
3086                    {
3087                        error=1;
3088                        break;
3089                    }
3090                    else        // x and y are in bounds
3091                    {
3092                        double e=table.getElt(x);
3093                        double w=table.getElt(x+1);
3094                        if ((e>undef) || (w>undef))
3095                        {
3096                            error=2;
3097                            break;
3098                        }
3099                // map x*Astep <= a << (x+1)*Astep to [-1,1]
3100                        double la = 2.0*(a-(x*Astep))/Astep-1;
3101                        rdat[l]=((1-la)*e + (1+la)*w)/2;
3102                    }
3103                }
3104            } catch (DataException d)
3105            {
3106                error=3;
3107                break;
3108            }
3109        } while (false);
3110    #ifdef PASO_MPI
3111        int rerror=0;
3112        MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );
3113        error=rerror;
3114    #endif
3115        if (error)
3116        {
3117            switch (error)
3118            {
3119                case 1: throw DataException("Point out of bounds");
3120                case 2: throw DataException("Interpolated value too large");
3121                default:
3122                    throw DataException("Unknown error in interpolation");      
3123            }
3124        }
3125        return res;
3126    }
3127    
3128                    
3129  Data  Data
3130  Data::interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,  Data::interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,

Legend:
Removed from v.2645  
changed lines
  Added in v.2646

  ViewVC Help
Powered by ViewVC 1.1.26