/[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 4668 by jfenwick, Tue Feb 11 03:44:10 2014 UTC revision 4687 by jfenwick, Wed Feb 19 00:03:29 2014 UTC
# Line 2911  namespace Line 2911  namespace
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.
# Line 2951  inset=2*radius+1 Line 2969  inset=2*radius+1
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        }
3000        radius=ex1();
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      }      }
3008      unsigned int radius=ex1();      else
     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 3003  escript::Data Brick::randomFill(long see Line 3038  escript::Data Brick::randomFill(long see
3038          throw RipleyException("Radius of gaussian filter is too large for Z dimension of a rank");          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]);
# Line 3055  for (int i=0;i<ext[0]*ext[1]*ext[2];) Line 3097  for (int i=0;i<ext[0]*ext[1]*ext[2];)
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 3116  for (int i=0;i<ext[0]*ext[1]*ext[2];) Line 3158  for (int i=0;i<ext[0]*ext[1]*ext[2];)
3158        cout << "\n";        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          
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 y=0;y<(internal[1]);++y)    
3174            {
3175            for (size_t x=0;x<(internal[0]);++x)
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::FunctionSpace fs(getPtr(), getContinuousFunctionCode());      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3192      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);      // don't need to check for exwrite because we just made it
3193      // don't need to check for exwrite because we just made it      escript::DataVector& dv=resdat.getExpandedVectorReference();
3194      escript::DataVector& dv=resdat.getExpandedVectorReference();      double* convolution=get3DGauss(radius, sigma);
     double* convolution=get3DGauss(radius, sigma);  
3195    
3196  // cout << "Convolution matrix\n";  // cout << "Convolution matrix\n";
3197  // size_t di=(radius*2+1);  // size_t di=(radius*2+1);
# Line 3143  for (int i=0;i<ext[0]*ext[1]*ext[2];) Line 3212  for (int i=0;i<ext[0]*ext[1]*ext[2];)
3212  //  //
3213  // cout << "Sum of matrix is = " << ctot << endl;  // cout << "Sum of matrix is = " << ctot << endl;
3214    
3215      for (size_t z=0;z<(internal[2]);++z)      for (size_t z=0;z<(internal[2]);++z)
     {  
     for (size_t y=0;y<(internal[1]);++y)      
3216      {      {
3217          for (size_t x=0;x<(internal[0]);++x)          for (size_t y=0;y<(internal[1]);++y)    
3218          {              {
3219          dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);          for (size_t x=0;x<(internal[0]);++x)
3220                    {    
3221                dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3222                
3223            }
3224          }          }
3225      }      }
     }  
3226            
3227  // cout << "\nResulting matrix:\n";  // cout << "\nResulting matrix:\n";
3228  //     for (size_t z=0;z<(internal[2]);++z)  //     for (size_t z=0;z<(internal[2]);++z)
# Line 3170  for (int i=0;i<ext[0]*ext[1]*ext[2];) Line 3239  for (int i=0;i<ext[0]*ext[1]*ext[2];)
3239  //     }  //     }
3240    
3241            
3242      delete[] convolution;      delete[] convolution;
3243      delete[] src;      delete[] src;
3244      return resdat;      return resdat;
3245        
3246        }
3247  }  }
3248    
3249    

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

  ViewVC Help
Powered by ViewVC 1.1.26