/[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 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC revision 2628 by jfenwick, Tue Aug 25 03:50:00 2009 UTC
# Line 3028  Data::freeSampleBuffer(BufferGroup* buff Line 3028  Data::freeSampleBuffer(BufferGroup* buff
3028  }  }
3029    
3030    
3031    Data
3032    Data::interpolateFromTable(boost::python::object table, double Amin, double Astep,
3033             double undef, Data& B, double Bmin, double Bstep)
3034    {
3035        WrappedArray t(table);
3036        return interpolateFromTable2D(t, Amin, Astep, undef, B, Bmin, Bstep);
3037    }
3038    
3039            
3040    Data
3041    Data::interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
3042                           double undef, Data& B, double Bmin, double Bstep)
3043    {
3044        int error=0;
3045        if ((getDataPointRank()!=0) || (B.getDataPointRank()!=0))
3046        {
3047            throw DataException("Inputs to 2D interpolation must be scalar");
3048        }
3049        if (table.getRank()!=2)
3050        {
3051        throw DataException("Table for 2D interpolation must be 2D");
3052        }
3053        if (getFunctionSpace()!=B.getFunctionSpace())
3054        {
3055        Data n=B.interpolate(getFunctionSpace());
3056        return interpolateFromTable2D(table, Amin, Astep, undef,
3057            n , Bmin, Bstep);
3058        }
3059        if (!isExpanded())
3060        {
3061        expand();
3062        }
3063        if (!B.isExpanded())
3064        {
3065        B.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            const DataVector& bdat=B.getReady()->getVectorRO();
3075            DataVector& rdat=res.getReady()->getVectorRW();
3076            const DataTypes::ShapeType& ts=table.getShape();
3077            for (int l=0; l<numpts; ++l)
3078            {
3079            double a=adat[l];
3080            double b=bdat[l];
3081            int x=static_cast<int>((a-Amin)/Astep);
3082            int y=static_cast<int>((b-Bmin)/Bstep);
3083            if ((a<Amin) || (b<Bmin))
3084            {
3085                error=1;
3086                break;  
3087            }
3088            if ((x>=(ts[0]-1)) || (y>=(ts[1]-1)))
3089            {
3090                error=1;
3091                break;
3092            }
3093            else        // x and y are in bounds
3094            {
3095                double sw=table.getElt(x,y);
3096                double nw=table.getElt(x,y+1);
3097                double se=table.getElt(x+1,y);
3098                double ne=table.getElt(x+1,y+1);
3099                if ((sw>undef) || (nw>undef) || (se>undef) || (ne>undef))
3100                {
3101                error=2;
3102                break;
3103                }
3104                // map x*Astep <= a << (x+1)*Astep to [-1,1]
3105                // same with b
3106                double la = 2.0*(a-(x*Astep))/Astep-1;
3107                double lb = 2.0*(b-(y*Bstep))/Bstep-1;
3108                rdat[l]=((1-la)*(1-lb)*sw + (1-la)*(1+lb)*nw +
3109                     (1+la)*(1-lb)*se + (1+la)*(1+lb)*ne)/4;
3110            }
3111            }
3112        } catch (DataException d)
3113        {
3114            error=3;
3115            break;
3116        }
3117        } while (false);
3118    #ifdef PASO_MPI
3119        int rerror=0;
3120        MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );
3121        error=rerror;
3122    #endif
3123        if (error)
3124        {
3125        switch (error)
3126        {
3127        case 1: throw DataException("Point out of bounds");
3128        case 2: throw DataException("Interpolated value too large");
3129        default:
3130            throw DataException("Unknown error in interpolation");      
3131        }
3132        }
3133        return res;
3134    }
3135    
3136    
3137  /* Member functions specific to the MPI implementation */  /* Member functions specific to the MPI implementation */
3138    
3139  void  void

Legend:
Removed from v.2548  
changed lines
  Added in v.2628

  ViewVC Help
Powered by ViewVC 1.1.26