/[escript]/branches/diaplayground/ripley/src/Brick.cpp
ViewVC logotype

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

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

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;
2884      for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)      for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)
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)
2909    //  return source[sbase+(radius)+(radius)*width+(radius)*width*height];
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    
2993      if ((2*radius>=internal[0]) || (2*radius>=internal[1]) || (2*radius>=internal[2]))      if (2*radius>=internal[0]-4)
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      }      }
2997            if (2*radius>=internal[1]-4)
2998        {
2999            throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3000        }
3001        if (2*radius>=internal[2]-4)
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();
3125      double* convolution=get3DGauss(radius, sigma);      double* convolution=get3DGauss(radius, sigma);
3126    
3127    // cout << "Convolution matrix\n";
3128    // size_t di=(radius*2+1);
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

  ViewVC Help
Powered by ViewVC 1.1.26