/[escript]/trunk/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /trunk/ripley/src/Rectangle.cpp

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

revision 4495 by caltinay, Fri Jul 5 02:19:47 2013 UTC revision 4521 by jfenwick, Mon Aug 26 11:51:30 2013 UTC
# Line 18  Line 18 
18  #include <esysUtils/esysFileWriter.h>  #include <esysUtils/esysFileWriter.h>
19    
20  #include <boost/scoped_array.hpp>  #include <boost/scoped_array.hpp>
21    #include "esysUtils/EsysRandom.h"
22    
23  #ifdef USE_NETCDF  #ifdef USE_NETCDF
24  #include <netcdfcpp.h>  #include <netcdfcpp.h>
# Line 3913  void Rectangle::assemblePDEBoundarySyste Line 3914  void Rectangle::assemblePDEBoundarySyste
3914      } // end of parallel section      } // end of parallel section
3915  }  }
3916    
3917    namespace
3918    {
3919        // Calculates a guassian blur colvolution matrix for 2D
3920        double* get2DGauss(unsigned radius, double sigma)
3921        {
3922            double* arr=new double[2*(radius*2+1)];
3923            double common=M_1_PI*0.5*1/(sigma*sigma);
3924        double total=0;
3925        for (size_t y=-radius;y<=radius;++y)
3926        {
3927            for (size_t x=-radius;x<=radius;++x)
3928            {
3929                arr[x+y*(radius*2+1)]=common*exp(-(x*x+y*y)/(2*sigma*sigma));
3930                total+=arr[x+y*(radius*2+1)];
3931            }
3932        }
3933        double invtotal=1/total;
3934        for (size_t p=0;p<(2*(radius*2+1));++p)
3935        {
3936            arr[p]*=invtotal;
3937        }
3938        return arr;
3939        }
3940        
3941        // applies conv to source to get a point.
3942        // (xp, yp) are the coords in the source matrix not the destination matrix
3943        double Convolve2D(double* conv, double* source, size_t xp, size_t yp, unsigned radius, size_t width)
3944        {
3945        double result=0;
3946        for (size_t y=-radius;y<=radius;++y)
3947        {
3948            for (size_t x=-radius;x<=radius;++x)
3949            {
3950                result+=conv[x+y*(2*radius+1)] * source[xp + yp*width];
3951            }
3952        }
3953            return result;      
3954        }
3955    }
3956    
3957    
3958    
3959    escript::Data Rectangle::randomFill(long seed, const boost::python::tuple& filter) const
3960    {
3961        if (m_mpiInfo->size!=1)
3962        {
3963            throw RipleyException("This type of random does not support MPI yet.");
3964        }
3965        if (m_numDim!=2)
3966        {
3967            throw RipleyException("Only 2D supported at this time.");
3968        }
3969        if (len(filter)!=3) {
3970            throw RipleyException("Unsupported random filter");
3971        }
3972        boost::python::extract<string> ex(filter[0]);
3973        if (!ex.check() || (ex()!="gaussian"))
3974        {
3975            throw RipleyException("Unsupported random filter");
3976        }
3977        boost::python::extract<unsigned int> ex1(filter[1]);
3978        if (!ex1.check())
3979        {
3980            throw RipleyException("Radius of gaussian filter must be a positive integer.");
3981        }
3982        unsigned int radius=ex1();
3983        double sigma=0.5;
3984        boost::python::extract<double> ex2(filter[2]);
3985        if (!ex2.check() || (sigma=ex2())<=0)
3986        {
3987            throw RipleyException("Sigma must be a postive floating point number.");
3988        }    
3989        size_t padding=max((unsigned)max((m_NE[0]-m_ownNE[0])/2, (m_NE[1]-m_ownNE[1])/2), radius);
3990        size_t width=(m_ownNE[0]+2*padding);    
3991        size_t dsize=width * (m_ownNE[1]+2*padding);
3992        
3993        double* src=new double[dsize];
3994        esysUtils::randomFillArray(seed, src, dsize);  
3995      
3996        // Now we need to copy the regions owned by other ranks over here  
3997      
3998        // Lets call that done for now
3999        escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
4000        escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
4001        // don't need to check for exwrite because we just made it
4002        escript::DataVector& dv=resdat.getExpandedVectorReference();
4003        double* convolution=get2DGauss(radius, sigma);
4004        for (size_t y=0;y<m_ownNE[1];++y)    
4005        {
4006            for (size_t x=0;x<m_ownNE[0];++x)
4007        {
4008            dv[x+y*width]=Convolve2D(convolution, src, x, y, radius, width);
4009        }
4010        }
4011        return resdat;
4012    }
4013    
4014    
4015    
4016    
4017  } // end of namespace ripley  } // end of namespace ripley
4018    

Legend:
Removed from v.4495  
changed lines
  Added in v.4521

  ViewVC Help
Powered by ViewVC 1.1.26