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

revision 4660 by sshaw, Fri Feb 7 03:08:49 2014 UTC revision 4668 by jfenwick, Tue Feb 11 03:44:10 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
2915
2916
2917  /* This routine produces a Data object filled with smoothed random data.  /* This routine produces a Data object filled with smoothed random data.
2918  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.
2919  A parameter radius  gives the size of the stencil used for the smoothing.  A parameter radius  gives the size of the stencil used for the smoothing.
# Line 2986  escript::Data Brick::randomFill(long see Line 2990  escript::Data Brick::randomFill(long see
2990      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
2991      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
2992
2994      {      {
2995          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");
2996      }      }
2998        {
2999            throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3000        }
3002        {
3003            throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");
3004        }
3005
3006      double* src=new double[ext[0]*ext[1]*ext[2]];      double* src=new double[ext[0]*ext[1]*ext[2]];
3007      esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]);        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]);
3008

3009  #ifdef ESYS_MPI  #ifdef ESYS_MPI
3010
3011      dim_t X=m_mpiInfo->rank%m_NX[0];      dim_t X=m_mpiInfo->rank%m_NX[0];
3012      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];
3013      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);      dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3014    #endif
3015
3016    /*
3017        // if we wanted to test a repeating pattern
3018        size_t basex=0;
3019        size_t basey=0;
3020        size_t basez=0;
3021    #ifdef ESYS_MPI
3022        basex=X*m_gNE[0]/m_NX[0];
3023        basey=Y*m_gNE[1]/m_NX[1];
3024        basez=Z*m_gNE[2]/m_NX[2];
3025    #endif
3026        if (seed==0)
3027        {
3028        seed=2; // since we are using the seed parameter as the spacing and 0 spacing causes an exception
3029        }
3030        esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, seed, basex, basey, basez);
3031    */
3032
3033
3034    /*
3035    cout << "Pattern:\n";
3036    for (int i=0;i<ext[0]*ext[1]*ext[2];)
3037    {
3038        cout << src[i++] << " ";
3039        if (i%ext[0]==0)
3040          cout << "\n";
3041        if (i%(ext[0]*ext[1])==0)
3042          cout << "\n";
3043    }*/
3044
3045
3046    #ifdef ESYS_MPI
3047
3048
3049
3050      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);
3051      size_t inset=2*radius+1;          size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3052            // there is an overlap of two points both of which need to have "radius" points on either side.
3053
3054      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
3055      size_t ymidlen=ext[1]-2*inset;        size_t ymidlen=ext[1]-2*inset;
# Line 3022  escript::Data Brick::randomFill(long see Line 3068  escript::Data Brick::randomFill(long see
3068      grid.generateOutNeighbours(X, Y, Z, outcoms);      grid.generateOutNeighbours(X, Y, Z, outcoms);
3069
3070
3071      block.copyUsedFromBuffer(src);      block.copyAllToBuffer(src);
3072
3073
3074      int comserr=0;          int comserr=0;
# Line 3053  escript::Data Brick::randomFill(long see Line 3099  escript::Data Brick::randomFill(long see
3099      }      }
3100
3101      block.copyUsedFromBuffer(src);      block.copyUsedFromBuffer(src);
3102
3103
3104
3105
3106  #endif      #endif
3107
3108    /*
3109    cout << "Pattern (post transfer):\n";
3110    for (int i=0;i<ext[0]*ext[1]*ext[2];)
3111    {
3112        cout << src[i++] << " ";
3113        if (i%ext[0]==0)
3114          cout << "\n";
3115        if (i%(ext[0]*ext[1])==0)
3116          cout << "\n";
3117    }   */
3118
3119
3120
3121      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3122      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3123      // don't need to check for exwrite because we just made it      // don't need to check for exwrite because we just made it
3124      escript::DataVector& dv=resdat.getExpandedVectorReference();      escript::DataVector& dv=resdat.getExpandedVectorReference();
3126
3127    // cout << "Convolution matrix\n";
3129    // double ctot=0;
3130    // for (int i=0;i<di*di*di;++i)
3131    // {
3132    //     cout << convolution[i] << " ";
3133    //     ctot+=convolution[i];
3134    //     if ((i+1)%di==0)
3135    //     {
3136    //  cout << "\n";
3137    //     }
3138    //     if ((i+1)%(di*di)==0)
3139    //     {
3140    //  cout << "\n";
3141    //     }
3142    // }
3143    //
3144    // cout << "Sum of matrix is = " << ctot << endl;
3145
3146      for (size_t z=0;z<(internal[2]);++z)      for (size_t z=0;z<(internal[2]);++z)
3147      {      {
3148      for (size_t y=0;y<(internal[1]);++y)          for (size_t y=0;y<(internal[1]);++y)
# Line 3071  escript::Data Brick::randomFill(long see Line 3154  escript::Data Brick::randomFill(long see
3154          }          }
3155      }      }
3156      }      }
3157
3158    // cout << "\nResulting matrix:\n";
3159    //     for (size_t z=0;z<(internal[2]);++z)
3160    //     {
3161    //  for (size_t y=0;y<(internal[1]);++y)
3162    //  {
3163    //      for (size_t x=0;x<(internal[0]);++x)
3164    //      {
3165    //      cout << dv[x+y*(internal[0])+z*internal[0]*internal[1]] << " ";
3166    //      }
3167    //      cout << endl;
3168    //  }
3169    //  cout << endl;
3170    //     }
3171
3172
3173      delete[] convolution;      delete[] convolution;
3174      delete[] src;      delete[] src;
3175      return resdat;      return resdat;

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