# Diff of /branches/diaplayground/ripley/src/Brick.cpp

revision 4660 by sshaw, Fri Feb 7 03:08:49 2014 UTC revision 4687 by jfenwick, Wed Feb 19 00:03:29 2014 UTC
# Line 2875  namespace Line 2875  namespace
2875          {          {
2876          for (int x=-r;x<=r;++x)          for (int x=-r;x<=r;++x)
2877          {                  {
2878              arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y+z*z)/(2*sigma*sigma));              arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)]=common*exp(-(x*x+y*y+z*z)/(2*sigma*sigma));
2879              total+=arr[(x+r)+(y+r)*(r*2+1)];              total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];
2880          }          }
2881          }          }
2882      }      }
2883      double invtotal=1/total;      double invtotal=1/total;
2885      {      {
2886          arr[p]*=invtotal;          arr[p]*=invtotal;
2887      }      }
# Line 2905  namespace Line 2905  namespace
2905          }          }
2906          }          }
2907      }      }
2908          return result;            // use this line for "pass-though" (return the centre point value)
2910        return result;
2911      }      }
2912  }  }
2913
2914    /* This is a wrapper for filtered (and non-filtered) randoms
2915     * For detailed doco see randomFillWorker
2916    */
2917    escript::Data Brick::randomFill(const escript::DataTypes::ShapeType& shape,
2918           const escript::FunctionSpace& what,
2919           long seed, const boost::python::tuple& filter) const
2920    {
2921        int numvals=escript::DataTypes::noValues(shape);
2922        if (len(filter)>0 && (numvals!=1))
2923        {
2924        throw RipleyException("Ripley only supports filters for scalar data.");
2925        }
2926        escript::Data res=randomFillWorker(shape, seed, filter);
2927        if (res.getFunctionSpace()!=what)
2928        {
2929        escript::Data r=escript::Data(res, what);
2930        return r;
2931        }
2932        return res;
2933    }
2934
2935  /* This routine produces a Data object filled with smoothed random data.  /* This routine produces a Data object filled with smoothed random data.
2936  The dimensions of the rectangle being filled are internal[0] x internal[1] x internal[2] points.  The dimensions of the rectangle being filled are internal[0] x internal[1] x internal[2] points.
2969  This is to ensure that values at distance `radius` from the shared/overlapped element  This is to ensure that values at distance `radius` from the shared/overlapped element
2970  that ripley has.  that ripley has.
2971  */  */
2972  escript::Data Brick::randomFill(long seed, const boost::python::tuple& filter) const  escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const
2973  {  {
2974      if (m_numDim!=3)      if (m_numDim!=3)
2975      {      {
2976          throw RipleyException("Brick must be 3D.");          throw RipleyException("Brick must be 3D.");
2977      }      }
2978      if (len(filter)!=3) {
2979          throw RipleyException("Unsupported random filter");      unsigned int radius=0;  // these are only used by gaussian
2980      }      double sigma=0.5;
2981      boost::python::extract<string> ex(filter[0]);
2982      if (!ex.check() || (ex()!="gaussian"))      unsigned int numvals=escript::DataTypes::noValues(shape);
2983
2984        if (len(filter)==0)
2985      {      {
2986          throw RipleyException("Unsupported random filter");      // nothing special required here yet
2987      }      }
2988      boost::python::extract<unsigned int> ex1(filter[1]);      else if (len(filter)==3)
if (!ex1.check())
2989      {      {
2990          throw RipleyException("Radius of gaussian filter must be a positive integer.");      boost::python::extract<string> ex(filter[0]);
2991        if (!ex.check() || (ex()!="gaussian"))
2992        {
2993            throw RipleyException("Unsupported random filter for Brick.");
2994        }
2995        boost::python::extract<unsigned int> ex1(filter[1]);
2996        if (!ex1.check())
2997        {
2998            throw RipleyException("Radius of gaussian filter must be a positive integer.");
2999        }
3001        sigma=0.5;
3002        boost::python::extract<double> ex2(filter[2]);
3003        if (!ex2.check() || (sigma=ex2())<=0)
3004        {
3005            throw RipleyException("Sigma must be a postive floating point number.");
3006        }
3007      }      }
double sigma=0.5;
boost::python::extract<double> ex2(filter[2]);
if (!ex2.check() || (sigma=ex2())<=0)
3009      {      {
3010          throw RipleyException("Sigma must be a postive floating point number.");          throw RipleyException("Unsupported random filter");
3011      }          }
3012
3013
3014
3015
3016      size_t internal[3];      size_t internal[3];
3017      internal[0]=m_NE[0]+1;  // number of points in the internal region      internal[0]=m_NE[0]+1;  // number of points in the internal region
# Line 2986  escript::Data Brick::randomFill(long see Line 3025  escript::Data Brick::randomFill(long see
3025      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
3026      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
3027
3029      {      {
3030          throw RipleyException("Radius of gaussian filter must be less than half the width/height of a rank");          throw RipleyException("Radius of gaussian filter is too large for X dimension of a rank");
3031      }      }
3033        {
3034            throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3035        }
3037        {
3038            throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");
3039        }
3040
3041      double* src=new double[ext[0]*ext[1]*ext[2]];      double* src=new double[ext[0]*ext[1]*ext[2]*numvals];
3042      esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]);        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);
3043

3044  #ifdef ESYS_MPI  #ifdef ESYS_MPI
3045
3046        if ((internal[0]<5) || (internal[1]<5) || (internal[2]<5))
3047        {
3048        // since the dimensions are equal for all ranks, this exception
3049        // will be thrown on all ranks
3050        throw RipleyException("Random Data in Ripley requries at least five elements per side per rank.");
3051
3052        }
3053      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
3054      dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];      dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
3055      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3056    #endif
3057
3058    /*
3059        // if we wanted to test a repeating pattern
3060        size_t basex=0;
3061        size_t basey=0;
3062        size_t basez=0;
3063    #ifdef ESYS_MPI
3064        basex=X*m_gNE[0]/m_NX[0];
3065        basey=Y*m_gNE[1]/m_NX[1];
3066        basez=Z*m_gNE[2]/m_NX[2];
3067    #endif
3068        if (seed==0)
3069        {
3070        seed=2; // since we are using the seed parameter as the spacing and 0 spacing causes an exception
3071        }
3072        esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, seed, basex, basey, basez);
3073    */
3074
3075
3076    /*
3077    cout << "Pattern:\n";
3078    for (int i=0;i<ext[0]*ext[1]*ext[2];)
3079    {
3080        cout << src[i++] << " ";
3081        if (i%ext[0]==0)
3082          cout << "\n";
3083        if (i%(ext[0]*ext[1])==0)
3084          cout << "\n";
3085    }*/
3086
3087
3088    #ifdef ESYS_MPI
3089
3090
3091
3092      BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);      BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);
3093      size_t inset=2*radius+1;          size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3094            // there is an overlap of two points both of which need to have "radius" points on either side.
3095
3096      size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets      size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets
3097      size_t ymidlen=ext[1]-2*inset;        size_t ymidlen=ext[1]-2*inset;
3098      size_t zmidlen=ext[2]-2*inset;      size_t zmidlen=ext[2]-2*inset;
3099
3100      Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen);          Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen, numvals);
3101
3102      MPI_Request reqs[50];       // a non-tight upper bound on how many we need      MPI_Request reqs[50];       // a non-tight upper bound on how many we need
3103      MPI_Status stats[50];      MPI_Status stats[50];
# Line 3022  escript::Data Brick::randomFill(long see Line 3110  escript::Data Brick::randomFill(long see
3110      grid.generateOutNeighbours(X, Y, Z, outcoms);      grid.generateOutNeighbours(X, Y, Z, outcoms);
3111
3112
3113      block.copyUsedFromBuffer(src);      block.copyAllToBuffer(src);
3114
3115
3116      int comserr=0;          int comserr=0;
# Line 3053  escript::Data Brick::randomFill(long see Line 3141  escript::Data Brick::randomFill(long see
3141      }      }
3142
3143      block.copyUsedFromBuffer(src);      block.copyUsedFromBuffer(src);
3144
3145
3146
3147
3148  #endif      #endif
3149      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3150      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);  /*
3151      // don't need to check for exwrite because we just made it  cout << "Pattern (post transfer):\n";
3152      escript::DataVector& dv=resdat.getExpandedVectorReference();  for (int i=0;i<ext[0]*ext[1]*ext[2];)
3154      for (size_t z=0;z<(internal[2]);++z)      cout << src[i++] << " ";
3155        if (i%ext[0]==0)
3156          cout << "\n";
3157        if (i%(ext[0]*ext[1])==0)
3158          cout << "\n";
3159    }   */
3160
3161        if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
3162      {      {
3163      for (size_t y=0;y<(internal[1]);++y)
3164        escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3165        escript::Data resdat(0, shape, fs , true);
3166        // don't need to check for exwrite because we just made it
3167        escript::DataVector& dv=resdat.getExpandedVectorReference();
3168
3169
3170        // now we need to copy values over
3171        for (size_t z=0;z<(internal[2]);++z)
3172      {      {
3173          for (size_t x=0;x<(internal[0]);++x)          for (size_t y=0;y<(internal[1]);++y)
3174          {              {
3176                    {
3177                for (unsigned int i=0;i<numvals;++i)
3178                {
3179                    dv[i+(x+y*(internal[0])+z*internal[0]*internal[1])*numvals]=src[i+(x+y*ext[0]+z*ext[0]*ext[1])*numvals];
3180                }
3181            }
3182            }
3183        }
3184        delete[] src;
3185        return resdat;
3186        }
3187        else        // filter enabled
3188        {
3189
3190        escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3191        escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3192        // don't need to check for exwrite because we just made it
3193        escript::DataVector& dv=resdat.getExpandedVectorReference();
3195
3196    // cout << "Convolution matrix\n";
3198    // double ctot=0;
3199    // for (int i=0;i<di*di*di;++i)
3200    // {
3201    //     cout << convolution[i] << " ";
3202    //     ctot+=convolution[i];
3203    //     if ((i+1)%di==0)
3204    //     {
3205    //  cout << "\n";
3206    //     }
3207    //     if ((i+1)%(di*di)==0)
3208    //     {
3209    //  cout << "\n";
3210    //     }
3211    // }
3212    //
3213    // cout << "Sum of matrix is = " << ctot << endl;
3214
3215        for (size_t z=0;z<(internal[2]);++z)
3216        {
3217            for (size_t y=0;y<(internal[1]);++y)
3218            {
3219            for (size_t x=0;x<(internal[0]);++x)
3220            {
3222
3223            }
3224          }          }
3225      }      }
3226
3227    // cout << "\nResulting matrix:\n";
3228    //     for (size_t z=0;z<(internal[2]);++z)
3229    //     {
3230    //  for (size_t y=0;y<(internal[1]);++y)
3231    //  {
3232    //      for (size_t x=0;x<(internal[0]);++x)
3233    //      {
3234    //      cout << dv[x+y*(internal[0])+z*internal[0]*internal[1]] << " ";
3235    //      }
3236    //      cout << endl;
3237    //  }
3238    //  cout << endl;
3239    //     }
3240
3241
3242        delete[] convolution;
3243        delete[] src;
3244        return resdat;
3245
3246      }      }
delete[] convolution;
delete[] src;
return resdat;
3247  }  }
3248
3249

Legend:
 Removed from v.4660 changed lines Added in v.4687