/[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 4633 by caltinay, Tue Jan 28 01:04:17 2014 UTC revision 4650 by jfenwick, Wed Feb 5 04:16:01 2014 UTC
# Line 33  Line 33 
33    
34  #include <iomanip>  #include <iomanip>
35    
36    #include "esysUtils/EsysRandom.h"
37    #include "blocktools.h"
38    
39    
40  using namespace std;  using namespace std;
41  using esysUtils::FileWriter;  using esysUtils::FileWriter;
42    
# Line 41  namespace ripley { Line 45  namespace ripley {
45  Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,  Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,
46               double x1, double y1, double z1, int d0, int d1, int d2,               double x1, double y1, double z1, int d0, int d1, int d2,
47               const std::vector<double>& points, const std::vector<int>& tags,               const std::vector<double>& points, const std::vector<int>& tags,
48               const std::map<std::string, int>& tagnamestonums) :               const simap_t& tagnamestonums) :
49      RipleyDomain(3)      RipleyDomain(3)
50  {  {
51      // ignore subdivision parameters for serial run      // ignore subdivision parameters for serial run
# Line 2852  void Brick::interpolateNodesOnFaces(escr Line 2856  void Brick::interpolateNodesOnFaces(escr
2856      }      }
2857  }  }
2858    
2859    namespace
2860    {
2861        // Calculates a guassian blur colvolution matrix for 3D
2862        // See wiki article on the subject
2863        double* get3DGauss(unsigned radius, double sigma)
2864        {
2865            double* arr=new double[(radius*2+1)*(radius*2+1)*(radius*2+1)];
2866            double common=pow(M_1_PI*0.5*1/(sigma*sigma), 3./2);
2867        double total=0;
2868        int r=static_cast<int>(radius);
2869        for (int z=-r;z<=r;++z)
2870        {
2871            for (int y=-r;y<=r;++y)
2872            {
2873            for (int x=-r;x<=r;++x)
2874            {        
2875                arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y+z*z)/(2*sigma*sigma));
2876                total+=arr[(x+r)+(y+r)*(r*2+1)];
2877            }
2878            }
2879        }
2880        double invtotal=1/total;
2881        for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)
2882        {
2883            arr[p]*=invtotal;
2884        }
2885        return arr;
2886        }
2887        
2888        // applies conv to source to get a point.
2889        // (xp, yp) are the coords in the source matrix not the destination matrix
2890        double Convolve3D(double* conv, double* source, size_t xp, size_t yp, size_t zp, unsigned radius, size_t width, size_t height)
2891        {
2892            size_t bx=xp-radius, by=yp-radius, bz=zp-radius;
2893        size_t sbase=bx+by*width+bz*width*height;
2894        double result=0;
2895        for (int z=0;z<2*radius+1;++z)
2896        {
2897            for (int y=0;y<2*radius+1;++y)
2898            {    
2899            for (int x=0;x<2*radius+1;++x)
2900            {
2901                result+=conv[x+y*(2*radius+1)+z*(2*radius+1)*(2*radius+1)] * source[sbase + x+y*width+z*width*height];
2902            }
2903            }
2904        }
2905            return result;      
2906        }
2907    }
2908    
2909    
2910    /* This routine produces a Data object filled with smoothed random data.
2911    The dimensions of the rectangle being filled are internal[0] x internal[1] x internal[2] points.
2912    A parameter radius  gives the size of the stencil used for the smoothing.
2913    A point on the left hand edge for example, will still require `radius` extra points to the left
2914    in order to complete the stencil.
2915    
2916    All local calculation is done on an array called `src`, with
2917    dimensions = ext[0] * ext[1] *ext[2].
2918    Where ext[i]= internal[i]+2*radius.
2919    
2920    Now for MPI there is overlap to deal with. We need to share both the overlapping
2921    values themselves but also the external region.
2922    
2923    In a hypothetical 1-D case:
2924    
2925    
2926    1234567
2927    would be split into two ranks thus:
2928    123(4)  (4)567     [4 being a shared element]
2929    
2930    If the radius is 2.   There will be padding elements on the outside:
2931    
2932    pp123(4)  (4)567pp
2933    
2934    To ensure that 4 can be correctly computed on both ranks, values from the other rank
2935    need to be known.
2936    
2937    pp123(4)56   23(4)567pp
2938    
2939    Now in our case, we wout set all the values 23456 on the left rank and send them to the
2940    right hand rank.
2941    
2942    So the edges _may_ need to be shared at a distance `inset` from all boundaries.
2943    inset=2*radius+1    
2944    This is to ensure that values at distance `radius` from the shared/overlapped element
2945    that ripley has.
2946    */
2947    escript::Data Brick::randomFill(long seed, const boost::python::tuple& filter) const
2948    {
2949        if (m_numDim!=3)
2950        {
2951            throw RipleyException("Brick must be 3D.");
2952        }
2953        if (len(filter)!=3) {
2954            throw RipleyException("Unsupported random filter");
2955        }
2956        boost::python::extract<string> ex(filter[0]);
2957        if (!ex.check() || (ex()!="gaussian"))
2958        {
2959            throw RipleyException("Unsupported random filter");
2960        }
2961        boost::python::extract<unsigned int> ex1(filter[1]);
2962        if (!ex1.check())
2963        {
2964            throw RipleyException("Radius of gaussian filter must be a positive integer.");
2965        }
2966        unsigned int radius=ex1();
2967        double sigma=0.5;
2968        boost::python::extract<double> ex2(filter[2]);
2969        if (!ex2.check() || (sigma=ex2())<=0)
2970        {
2971            throw RipleyException("Sigma must be a postive floating point number.");
2972        }    
2973        
2974        size_t internal[3];
2975        internal[0]=m_NE[0]+1;  // number of points in the internal region
2976        internal[1]=m_NE[1]+1;  // that is, the ones we need smoothed versions of
2977        internal[2]=m_NE[2]+1;  // that is, the ones we need smoothed versions of
2978        size_t ext[3];
2979        ext[0]=internal[0]+2*radius;    // includes points we need as input
2980        ext[1]=internal[1]+2*radius;    // for smoothing
2981        ext[2]=internal[2]+2*radius;    // for smoothing
2982        
2983        // now we check to see if the radius is acceptable
2984        // That is, would not cross multiple ranks in MPI
2985    
2986        if ((2*radius>=internal[0]) || (2*radius>=internal[1]) || (2*radius>=internal[2]))
2987        {
2988            throw RipleyException("Radius of gaussian filter must be less than half the width/height of a rank");
2989        }
2990        
2991      
2992        double* src=new double[ext[0]*ext[1]*ext[2]];
2993        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]);  
2994        
2995      
2996    #ifdef ESYS_MPI
2997        
2998        dim_t X=m_mpiInfo->rank%m_NX[0];
2999        dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
3000        dim_t Z=m_mpiInfo->rank/(m_NX[0]*m_NX[1]);
3001        
3002        BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);
3003        size_t inset=2*radius+1;    
3004        
3005        size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets
3006        size_t ymidlen=ext[1]-2*inset;  
3007        size_t zmidlen=ext[2]-2*inset;
3008        
3009        Block block(ext[0], ext[1], ext[2], inset, xmidlen, ymidlen, zmidlen);    
3010        
3011        MPI_Request reqs[50];       // a non-tight upper bound on how many we need
3012        MPI_Status stats[50];
3013        short rused=0;
3014        
3015        messvec incoms;
3016        messvec outcoms;
3017        
3018        grid.generateInNeighbours(X, Y, Z ,incoms);
3019        grid.generateOutNeighbours(X, Y, Z, outcoms);
3020        
3021        
3022        block.copyUsedFromBuffer(src);
3023        
3024        
3025        int comserr=0;    
3026        for (size_t i=0;i<incoms.size();++i)
3027        {
3028        message& m=incoms[i];
3029        comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3030        block.setUsed(m.destbuffid);
3031        }
3032    
3033        for (size_t i=0;i<incoms.size();++i)
3034        {
3035        message& m=incoms[i];
3036        comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3037        }    
3038        
3039        if (!comserr)
3040        {    
3041            comserr=MPI_Waitall(rused, reqs, stats);
3042        }
3043    
3044        if (comserr)
3045        {
3046        // Yes this is throwing an exception as a result of an MPI error.
3047        // and no we don't inform the other ranks that we are doing this.
3048        // however, we have no reason to believe coms work at this point anyway
3049            throw RipleyException("Error in coms for randomFill");      
3050        }
3051        
3052        block.copyUsedFromBuffer(src);
3053    
3054    #endif    
3055        escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
3056        escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3057        // don't need to check for exwrite because we just made it
3058        escript::DataVector& dv=resdat.getExpandedVectorReference();
3059        double* convolution=get3DGauss(radius, sigma);
3060        for (size_t z=0;z<(internal[2]);++z)
3061        {
3062        for (size_t y=0;y<(internal[1]);++y)    
3063        {
3064            for (size_t x=0;x<(internal[0]);++x)
3065            {    
3066            dv[x+y*(internal[0])+z*internal[0]*internal[1]]=Convolve3D(convolution, src, x+radius, y+radius, z+radius, radius, ext[0], ext[1]);
3067            
3068            }
3069        }
3070        }
3071        delete[] convolution;
3072        delete[] src;
3073        return resdat;
3074    }
3075    
3076    
3077    
3078    
3079    
3080    
3081  int Brick::findNode(const double *coords) const {  int Brick::findNode(const double *coords) const {
3082      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3083      //is the found element even owned by this rank      //is the found element even owned by this rank

Legend:
Removed from v.4633  
changed lines
  Added in v.4650

  ViewVC Help
Powered by ViewVC 1.1.26