/[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 4650 by jfenwick, Wed Feb 5 04:16:01 2014 UTC revision 4696 by jfenwick, Wed Feb 19 07:29:50 2014 UTC
# Line 1  Line 1 
1    
2  /*****************************************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2013 by University of Queensland  * Copyright (c) 2003-2014 by University of Queensland
5  * http://www.uq.edu.au  * http://www.uq.edu.au
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
# Line 9  Line 9 
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12  * Development since 2012 by School of Earth Sciences  * Development 2012-2013 by School of Earth Sciences
13    * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
# Line 786  const int* Brick::borrowSampleReferenceI Line 787  const int* Brick::borrowSampleReferenceI
787          case FaceElements:          case FaceElements:
788          case ReducedFaceElements:          case ReducedFaceElements:
789              return &m_faceId[0];              return &m_faceId[0];
790            case Points:
791                return &m_diracPointNodeIDs[0];
792          default:          default:
793              break;              break;
794      }      }
# Line 2872  namespace Line 2875  namespace
2875          {          {
2876          for (int x=-r;x<=r;++x)          for (int x=-r;x<=r;++x)
2877          {                  {        
2878              arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y+z*z)/(2*sigma*sigma));              arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)]=common*exp(-(x*x+y*y+z*z)/(2*sigma*sigma));
2879              total+=arr[(x+r)+(y+r)*(r*2+1)];              total+=arr[(x+r)+(y+r)*(r*2+1)+(z+r)*(r*2+1)*(r*2+1)];
2880          }          }
2881          }          }
2882      }      }
2883      double invtotal=1/total;      double invtotal=1/total;
2884      for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)      for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)
2885      {      {
2886          arr[p]*=invtotal;          arr[p]*=invtotal;
2887      }      }
# Line 2902  namespace Line 2905  namespace
2905          }          }
2906          }          }
2907      }      }
2908          return result;            // use this line for "pass-though" (return the centre point value)
2909    //  return source[sbase+(radius)+(radius)*width+(radius)*width*height];
2910        return result;      
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 2944  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 2983  escript::Data Brick::randomFill(long see Line 3025  escript::Data Brick::randomFill(long see
3025      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
3026      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
3027    
3028      if ((2*radius>=internal[0]) || (2*radius>=internal[1]) || (2*radius>=internal[2]))      if (2*radius>=internal[0]-4)
3029      {      {
3030          throw RipleyException("Radius of gaussian filter must be less than half the width/height of a rank");          throw RipleyException("Radius of gaussian filter is too large for X dimension of a rank");
3031      }      }
3032            if (2*radius>=internal[1]-4)
3033        {
3034            throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
3035        }
3036        if (2*radius>=internal[2]-4)
3037        {
3038            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]);
3056    #endif    
3057    
3058    /*    
3059        // if we wanted to test a repeating pattern
3060        size_t basex=0;
3061        size_t basey=0;
3062        size_t basez=0;
3063    #ifdef ESYS_MPI    
3064        basex=X*m_gNE[0]/m_NX[0];
3065        basey=Y*m_gNE[1]/m_NX[1];
3066        basez=Z*m_gNE[2]/m_NX[2];    
3067    #endif    
3068        if (seed==0)
3069        {
3070        seed=2; // since we are using the seed parameter as the spacing and 0 spacing causes an exception
3071        }
3072        esysUtils::patternFillArray(1, ext[0],ext[1],ext[2], src, 4, basex, basey, basez, numvals);
3073    */
3074            
3075        
3076    /*
3077    cout << "Pattern:\n";    
3078    for (int i=0;i<ext[0]*ext[1]*ext[2];)
3079    {
3080        cout << src[i++] << " ";
3081        if (i%ext[0]==0)
3082          cout << "\n";
3083        if (i%(ext[0]*ext[1])==0)
3084          cout << "\n";
3085    }*/
3086        
3087      
3088    #ifdef ESYS_MPI
3089    
3090    
3091    
3092      BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);      BlockGrid grid(m_NX[0]-1, m_NX[1]-1, m_NX[2]-1);
3093      size_t inset=2*radius+1;          size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
3094            // there is an overlap of two points both of which need to have "radius" points on either side.
3095            
3096      size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets      size_t xmidlen=ext[0]-2*inset;  // how wide is the x-dimension between the two insets
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 3019  escript::Data Brick::randomFill(long see Line 3110  escript::Data Brick::randomFill(long see
3110      grid.generateOutNeighbours(X, Y, Z, outcoms);      grid.generateOutNeighbours(X, Y, Z, outcoms);
3111            
3112            
3113      block.copyUsedFromBuffer(src);      block.copyAllToBuffer(src);
3114            
3115            
3116      int comserr=0;          int comserr=0;    
# Line 3030  escript::Data Brick::randomFill(long see Line 3121  escript::Data Brick::randomFill(long see
3121      block.setUsed(m.destbuffid);      block.setUsed(m.destbuffid);
3122      }      }
3123    
3124      for (size_t i=0;i<incoms.size();++i)      for (size_t i=0;i<outcoms.size();++i)
3125      {      {
3126      message& m=incoms[i];      message& m=outcoms[i];
3127      comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));      comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3128      }          }    
3129            
# Line 3050  escript::Data Brick::randomFill(long see Line 3141  escript::Data Brick::randomFill(long see
3141      }      }
3142            
3143      block.copyUsedFromBuffer(src);      block.copyUsedFromBuffer(src);
3144        
3145        
3146        
3147    
3148  #endif      #endif    
3149      escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());      
3150      escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);  /*    
3151      // don't need to check for exwrite because we just made it  cout << "Pattern (post transfer):\n";    
3152      escript::DataVector& dv=resdat.getExpandedVectorReference();  for (int i=0;i<ext[0]*ext[1]*ext[2];)
3153      double* convolution=get3DGauss(radius, sigma);  {
3154      for (size_t z=0;z<(internal[2]);++z)      cout << src[i++] << " ";
3155        if (i%ext[0]==0)
3156          cout << "\n";
3157        if (i%(ext[0]*ext[1])==0)
3158          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      for (size_t y=0;y<(internal[1]);++y)            
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 x=0;x<(internal[0]);++x)          for (size_t y=0;y<(internal[1]);++y)    
3174          {              {
3175          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)
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::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
3192        // don't need to check for exwrite because we just made it
3193        escript::DataVector& dv=resdat.getExpandedVectorReference();
3194        double* convolution=get3DGauss(radius, sigma);
3195    
3196    // cout << "Convolution matrix\n";
3197    // size_t di=(radius*2+1);
3198    // double ctot=0;
3199    // for (int i=0;i<di*di*di;++i)
3200    // {
3201    //     cout << convolution[i] << " ";
3202    //     ctot+=convolution[i];
3203    //     if ((i+1)%di==0)
3204    //     {
3205    //  cout << "\n";
3206    //     }
3207    //     if ((i+1)%(di*di)==0)
3208    //     {
3209    //  cout << "\n";
3210    //     }
3211    // }
3212    //
3213    // cout << "Sum of matrix is = " << ctot << endl;
3214    
3215        for (size_t z=0;z<(internal[2]);++z)
3216        {
3217            for (size_t y=0;y<(internal[1]);++y)    
3218            {
3219            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";
3228    //     for (size_t z=0;z<(internal[2]);++z)
3229    //     {
3230    //  for (size_t y=0;y<(internal[1]);++y)    
3231    //  {
3232    //      for (size_t x=0;x<(internal[0]);++x)
3233    //      {    
3234    //      cout << dv[x+y*(internal[0])+z*internal[0]*internal[1]] << " ";
3235    //      }
3236    //      cout << endl;
3237    //  }
3238    //  cout << endl;
3239    //     }
3240    
3241        
3242        delete[] convolution;
3243        delete[] src;
3244        return resdat;
3245        
3246      }      }
     delete[] convolution;  
     delete[] src;  
     return resdat;  
3247  }  }
3248    
3249    
# Line 3081  escript::Data Brick::randomFill(long see Line 3254  escript::Data Brick::randomFill(long see
3254  int Brick::findNode(const double *coords) const {  int Brick::findNode(const double *coords) const {
3255      const int NOT_MINE = -1;      const int NOT_MINE = -1;
3256      //is the found element even owned by this rank      //is the found element even owned by this rank
3257        // (inside owned or shared elements but will map to an owned element)
3258      for (int dim = 0; dim < m_numDim; dim++) {      for (int dim = 0; dim < m_numDim; dim++) {
3259          if (m_origin[dim] + m_offset[dim] > coords[dim]  || m_origin[dim]          double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
3260                  + m_offset[dim] + m_dx[dim]*m_ownNE[dim] < coords[dim]) {                  - m_dx[dim]/2.; //allows for point outside mapping onto node
3261            double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
3262                    + m_dx[dim]/2.;
3263            if (min > coords[dim] || max < coords[dim]) {
3264              return NOT_MINE;              return NOT_MINE;
3265          }          }
3266      }      }
# Line 3102  int Brick::findNode(const double *coords Line 3279  int Brick::findNode(const double *coords
3279          minDist += m_dx[dim]*m_dx[dim];          minDist += m_dx[dim]*m_dx[dim];
3280      }      }
3281      //find the closest node      //find the closest node
3282      for (int dx = 0; dx < 2; dx++) {      for (int dx = 0; dx < 1; dx++) {
3283          double xdist = x - (ex + dx)*m_dx[0];          double xdist = x - (ex + dx)*m_dx[0];
3284          for (int dy = 0; dy < 2; dy++) {          for (int dy = 0; dy < 1; dy++) {
3285              double ydist = y - (ey + dy)*m_dx[1];              double ydist = y - (ey + dy)*m_dx[1];
3286              for (int dz = 0; dz < 2; dz++) {              for (int dz = 0; dz < 1; dz++) {
3287                  double zdist = z - (ez + dz)*m_dx[2];                  double zdist = z - (ez + dz)*m_dx[2];
3288                  double total = xdist*xdist + ydist*ydist + zdist*zdist;                  double total = xdist*xdist + ydist*ydist + zdist*zdist;
3289                  if (total < minDist) {                  if (total < minDist) {
3290                      closest = INDEX3(ex+dy, ey+dy, ez+dz, m_NE[0]+1, m_NE[1]+1);                      closest = INDEX3(ex+dy-m_offset[0], ey+dy-m_offset[1],
3291                                ez+dz-m_offset[2], m_NE[0]+1, m_NE[1]+1);
3292                      minDist = total;                      minDist = total;
3293                  }                  }
3294              }              }
3295          }          }
3296      }      }
3297        if (closest == NOT_MINE) {
3298            throw RipleyException("Unable to map appropriate dirac point to a node, implementation problem in Brick::findNode()");
3299        }
3300      return closest;      return closest;
3301  }  }
3302    

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

  ViewVC Help
Powered by ViewVC 1.1.26