/[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 4622 by sshaw, Fri Jan 17 04:55:41 2014 UTC revision 4687 by jfenwick, Wed Feb 19 00:03:29 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 17  Line 18 
18  #include <paso/SystemMatrix.h>  #include <paso/SystemMatrix.h>
19  #include <esysUtils/esysFileWriter.h>  #include <esysUtils/esysFileWriter.h>
20  #include <ripley/DefaultAssembler3D.h>  #include <ripley/DefaultAssembler3D.h>
21    #include <ripley/WaveAssembler3D.h>
22  #include <boost/scoped_array.hpp>  #include <boost/scoped_array.hpp>
23    
24  #ifdef USE_NETCDF  #ifdef USE_NETCDF
# Line 32  Line 34 
34    
35  #include <iomanip>  #include <iomanip>
36    
37    #include "esysUtils/EsysRandom.h"
38    #include "blocktools.h"
39    
40    
41  using namespace std;  using namespace std;
42  using esysUtils::FileWriter;  using esysUtils::FileWriter;
43    
# Line 40  namespace ripley { Line 46  namespace ripley {
46  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,
47               double x1, double y1, double z1, int d0, int d1, int d2,               double x1, double y1, double z1, int d0, int d1, int d2,
48               const std::vector<double>& points, const std::vector<int>& tags,               const std::vector<double>& points, const std::vector<int>& tags,
49               const std::map<std::string, int>& tagnamestonums) :               const simap_t& tagnamestonums) :
50      RipleyDomain(3)      RipleyDomain(3)
51  {  {
52      // ignore subdivision parameters for serial run      // ignore subdivision parameters for serial run
# Line 225  Brick::Brick(int n0, int n1, int n2, dou Line 231  Brick::Brick(int n0, int n1, int n2, dou
231      createPattern();      createPattern();
232            
233      assembler = new DefaultAssembler3D(this, m_dx, m_NX, m_NE, m_NN);      assembler = new DefaultAssembler3D(this, m_dx, m_NX, m_NE, m_NN);
234        for (map<string, int>::const_iterator i = tagnamestonums.begin();
235                i != tagnamestonums.end(); i++) {
236            setTagMap(i->first, i->second);
237        }
238      addPoints(tags.size(), &points[0], &tags[0]);      addPoints(tags.size(), &points[0], &tags[0]);
239  }  }
240    
# Line 590  void Brick::writeBinaryGridImpl(const es Line 600  void Brick::writeBinaryGridImpl(const es
600      if (numComp > 1 || dpp > 1)      if (numComp > 1 || dpp > 1)
601          throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");          throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");
602    
     escript::Data* _in = const_cast<escript::Data*>(&in);  
   
603      // from here on we know that each sample consists of one value      // from here on we know that each sample consists of one value
604      FileWriter fw;      FileWriter fw;
605      fw.openFile(filename, fileSize);      fw.openFile(filename, fileSize);
# Line 604  void Brick::writeBinaryGridImpl(const es Line 612  void Brick::writeBinaryGridImpl(const es
612              ostringstream oss;              ostringstream oss;
613    
614              for (index_t x=0; x<myN0; x++) {              for (index_t x=0; x<myN0; x++) {
615                  const double* sample = _in->getSampleDataRO(z*myN0*myN1+y*myN0+x);                  const double* sample = in.getSampleDataRO(z*myN0*myN1+y*myN0+x);
616                  ValueType fvalue = static_cast<ValueType>(*sample);                  ValueType fvalue = static_cast<ValueType>(*sample);
617                  if (byteOrder == BYTEORDER_NATIVE) {                  if (byteOrder == BYTEORDER_NATIVE) {
618                      oss.write((char*)&fvalue, sizeof(fvalue));                      oss.write((char*)&fvalue, sizeof(fvalue));
# Line 779  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 1151  void Brick::assembleCoordinates(escript: Line 1161  void Brick::assembleCoordinates(escript:
1161  }  }
1162    
1163  //protected  //protected
1164  void Brick::assembleGradient(escript::Data& out, escript::Data& in) const  void Brick::assembleGradient(escript::Data& out, const escript::Data& in) const
1165  {  {
1166      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1167      const double C0 = .044658198738520451079;      const double C0 = .044658198738520451079;
# Line 1622  void Brick::assembleGradient(escript::Da Line 1632  void Brick::assembleGradient(escript::Da
1632  }  }
1633    
1634  //protected  //protected
1635  void Brick::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const  void Brick::assembleIntegrate(vector<double>& integrals, const escript::Data& arg) const
1636  {  {
1637      const dim_t numComp = arg.getDataPointSize();      const dim_t numComp = arg.getDataPointSize();
1638      const index_t left = (m_offset[0]==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
# Line 1910  dim_t Brick::insertNeighbourNodes(IndexV Line 1920  dim_t Brick::insertNeighbourNodes(IndexV
1920  }  }
1921    
1922  //protected  //protected
1923  void Brick::nodesToDOF(escript::Data& out, escript::Data& in) const  void Brick::nodesToDOF(escript::Data& out, const escript::Data& in) const
1924  {  {
1925      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1926      out.requireWrite();      out.requireWrite();
# Line 1934  void Brick::nodesToDOF(escript::Data& ou Line 1944  void Brick::nodesToDOF(escript::Data& ou
1944  }  }
1945    
1946  //protected  //protected
1947  void Brick::dofToNodes(escript::Data& out, escript::Data& in) const  void Brick::dofToNodes(escript::Data& out, const escript::Data& in) const
1948  {  {
1949      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1950      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);
1951      in.requireWrite();      // expand data object if necessary to be able to grab the whole data
1952      Paso_Coupler_startCollect(coupler, in.getSampleDataRW(0));      const_cast<escript::Data*>(&in)->expand();
1953        Paso_Coupler_startCollect(coupler, in.getSampleDataRO(0));
1954    
1955      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
1956      out.requireWrite();      out.requireWrite();
# Line 2530  void Brick::addToMatrixAndRHS(Paso_Syste Line 2541  void Brick::addToMatrixAndRHS(Paso_Syste
2541  }  }
2542    
2543  //protected  //protected
2544  void Brick::interpolateNodesOnElements(escript::Data& out, escript::Data& in,  void Brick::interpolateNodesOnElements(escript::Data& out,
2545                                           const escript::Data& in,
2546                                         bool reduced) const                                         bool reduced) const
2547  {  {
2548      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
# Line 2613  void Brick::interpolateNodesOnElements(e Line 2625  void Brick::interpolateNodesOnElements(e
2625  }  }
2626    
2627  //protected  //protected
2628  void Brick::interpolateNodesOnFaces(escript::Data& out, escript::Data& in,  void Brick::interpolateNodesOnFaces(escript::Data& out, const escript::Data& in,
2629                                      bool reduced) const                                      bool reduced) const
2630  {  {
2631      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
# Line 2847  void Brick::interpolateNodesOnFaces(escr Line 2859  void Brick::interpolateNodesOnFaces(escr
2859      }      }
2860  }  }
2861    
2862  int Brick::findNode(double *coords) const {  namespace
2863    {
2864        // Calculates a guassian blur colvolution matrix for 3D
2865        // See wiki article on the subject
2866        double* get3DGauss(unsigned radius, double sigma)
2867        {
2868            double* arr=new double[(radius*2+1)*(radius*2+1)*(radius*2+1)];
2869            double common=pow(M_1_PI*0.5*1/(sigma*sigma), 3./2);
2870        double total=0;
2871        int r=static_cast<int>(radius);
2872        for (int z=-r;z<=r;++z)
2873        {
2874            for (int y=-r;y<=r;++y)
2875            {
2876            for (int x=-r;x<=r;++x)
2877            {        
2878                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)+(z+r)*(r*2+1)*(r*2+1)];
2880            }
2881            }
2882        }
2883        double invtotal=1/total;
2884        for (size_t p=0;p<(radius*2+1)*(radius*2+1)*(radius*2+1);++p)
2885        {
2886            arr[p]*=invtotal;
2887        }
2888        return arr;
2889        }
2890        
2891        // applies conv to source to get a point.
2892        // (xp, yp) are the coords in the source matrix not the destination matrix
2893        double Convolve3D(double* conv, double* source, size_t xp, size_t yp, size_t zp, unsigned radius, size_t width, size_t height)
2894        {
2895            size_t bx=xp-radius, by=yp-radius, bz=zp-radius;
2896        size_t sbase=bx+by*width+bz*width*height;
2897        double result=0;
2898        for (int z=0;z<2*radius+1;++z)
2899        {
2900            for (int y=0;y<2*radius+1;++y)
2901            {    
2902            for (int x=0;x<2*radius+1;++x)
2903            {
2904                result+=conv[x+y*(2*radius+1)+z*(2*radius+1)*(2*radius+1)] * source[sbase + x+y*width+z*width*height];
2905            }
2906            }
2907        }
2908        // 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.
2936    The dimensions of the rectangle being filled are internal[0] x internal[1] x internal[2] points.
2937    A parameter radius  gives the size of the stencil used for the smoothing.
2938    A point on the left hand edge for example, will still require `radius` extra points to the left
2939    in order to complete the stencil.
2940    
2941    All local calculation is done on an array called `src`, with
2942    dimensions = ext[0] * ext[1] *ext[2].
2943    Where ext[i]= internal[i]+2*radius.
2944    
2945    Now for MPI there is overlap to deal with. We need to share both the overlapping
2946    values themselves but also the external region.
2947    
2948    In a hypothetical 1-D case:
2949    
2950    
2951    1234567
2952    would be split into two ranks thus:
2953    123(4)  (4)567     [4 being a shared element]
2954    
2955    If the radius is 2.   There will be padding elements on the outside:
2956    
2957    pp123(4)  (4)567pp
2958    
2959    To ensure that 4 can be correctly computed on both ranks, values from the other rank
2960    need to be known.
2961    
2962    pp123(4)56   23(4)567pp
2963    
2964    Now in our case, we wout set all the values 23456 on the left rank and send them to the
2965    right hand rank.
2966    
2967    So the edges _may_ need to be shared at a distance `inset` from all boundaries.
2968    inset=2*radius+1    
2969    This is to ensure that values at distance `radius` from the shared/overlapped element
2970    that ripley has.
2971    */
2972    escript::Data Brick::randomFillWorker(const escript::DataTypes::ShapeType& shape, long seed, const boost::python::tuple& filter) const
2973    {
2974        if (m_numDim!=3)
2975        {
2976            throw RipleyException("Brick must be 3D.");
2977        }
2978        
2979        unsigned int radius=0;  // these are only used by gaussian
2980        double sigma=0.5;
2981        
2982        unsigned int numvals=escript::DataTypes::noValues(shape);
2983        
2984        if (len(filter)==0)
2985        {
2986        // nothing special required here yet
2987        }
2988        else if (len(filter)==3)
2989        {
2990        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        else
3009        {
3010            throw RipleyException("Unsupported random filter");
3011        }
3012        
3013        
3014    
3015        
3016        size_t internal[3];
3017        internal[0]=m_NE[0]+1;  // number of points in the internal region
3018        internal[1]=m_NE[1]+1;  // that is, the ones we need smoothed versions of
3019        internal[2]=m_NE[2]+1;  // that is, the ones we need smoothed versions of
3020        size_t ext[3];
3021        ext[0]=internal[0]+2*radius;    // includes points we need as input
3022        ext[1]=internal[1]+2*radius;    // for smoothing
3023        ext[2]=internal[2]+2*radius;    // for smoothing
3024        
3025        // now we check to see if the radius is acceptable
3026        // That is, would not cross multiple ranks in MPI
3027    
3028        if (2*radius>=internal[0]-4)
3029        {
3030            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]*numvals];
3042        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*ext[2]*numvals);      
3043        
3044    #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];
3054        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]);
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, seed, basex, basey, basez);
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);
3093        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
3097        size_t ymidlen=ext[1]-2*inset;  
3098        size_t zmidlen=ext[2]-2*inset;
3099        
3100        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
3103        MPI_Status stats[50];
3104        short rused=0;
3105        
3106        messvec incoms;
3107        messvec outcoms;
3108        
3109        grid.generateInNeighbours(X, Y, Z ,incoms);
3110        grid.generateOutNeighbours(X, Y, Z, outcoms);
3111        
3112        
3113        block.copyAllToBuffer(src);
3114        
3115        
3116        int comserr=0;    
3117        for (size_t i=0;i<incoms.size();++i)
3118        {
3119        message& m=incoms[i];
3120        comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));
3121        block.setUsed(m.destbuffid);
3122        }
3123    
3124        for (size_t i=0;i<outcoms.size();++i)
3125        {
3126        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++));
3128        }    
3129        
3130        if (!comserr)
3131        {    
3132            comserr=MPI_Waitall(rused, reqs, stats);
3133        }
3134    
3135        if (comserr)
3136        {
3137        // Yes this is throwing an exception as a result of an MPI error.
3138        // and no we don't inform the other ranks that we are doing this.
3139        // however, we have no reason to believe coms work at this point anyway
3140            throw RipleyException("Error in coms for randomFill");      
3141        }
3142        
3143        block.copyUsedFromBuffer(src);
3144        
3145        
3146        
3147    
3148    #endif    
3149        
3150    /*    
3151    cout << "Pattern (post transfer):\n";    
3152    for (int i=0;i<ext[0]*ext[1]*ext[2];)
3153    {
3154        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          
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 y=0;y<(internal[1]);++y)    
3174            {
3175            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        }
3247    }
3248    
3249    
3250    
3251    
3252    
3253    
3254    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 2871  int Brick::findNode(double *coords) cons Line 3279  int Brick::findNode(double *coords) cons
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    
3303    void Brick::setAssembler(std::string type, std::map<std::string,
3304            escript::Data> constants) {
3305        if (type.compare("WaveAssembler") == 0) {
3306            delete assembler;
3307            assembler = new WaveAssembler3D(this, m_dx, m_NX, m_NE, m_NN, constants);
3308        } else { //else ifs would go before this for other types
3309            throw RipleyException("Ripley::Rectangle does not support the"
3310                                    " requested assembler");
3311        }
3312    }
3313    
3314  } // end of namespace ripley  } // end of namespace ripley
3315    

Legend:
Removed from v.4622  
changed lines
  Added in v.4687

  ViewVC Help
Powered by ViewVC 1.1.26