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

revision 4638 by jfenwick, Thu Jan 30 06:25:10 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 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        {
2866            double common=pow(M_1_PI*0.5*1/(sigma*sigma), 3./2);
2867        double total=0;
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;
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        {
2893        size_t sbase=bx+by*width+bz*width*height;
2894        double result=0;
2896        {
2898            {
2900            {
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].
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.
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        }
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
2982
2983        // now we check to see if the radius is acceptable
2984        // That is, would not cross multiple ranks in MPI
2985
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);
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();
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            {
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.4638 changed lines Added in v.4650