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

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        {
3923            double common=M_1_PI*0.5*1/(sigma*sigma);
3924        double total=0;
3926        {
3928            {
3931            }
3932        }
3933        double invtotal=1/total;
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;
3947        {
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        }
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        }
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();
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