/[escript]/trunk/ripley/src/Rectangle.cpp
ViewVC logotype

Diff of /trunk/ripley/src/Rectangle.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 4765 by sshaw, Wed Mar 19 00:17:16 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    
17    #include <algorithm>
18    
19  #include <ripley/Rectangle.h>  #include <ripley/Rectangle.h>
20  #include <paso/SystemMatrix.h>  #include <paso/SystemMatrix.h>
21  #include <esysUtils/esysFileWriter.h>  #include <esysUtils/esysFileWriter.h>
22  #include <ripley/DefaultAssembler2D.h>  #include <ripley/DefaultAssembler2D.h>
23  #include <ripley/WaveAssembler2D.h>  #include <ripley/WaveAssembler2D.h>
24    #include <ripley/LameAssembler2D.h>
25    #include <ripley/domainhelpers.h>
26  #include <boost/scoped_array.hpp>  #include <boost/scoped_array.hpp>
27  #include "esysUtils/EsysRandom.h"  #include "esysUtils/EsysRandom.h"
28    #include "blocktools.h"
29    
30  #ifdef USE_NETCDF  #ifdef USE_NETCDF
31  #include <netcdfcpp.h>  #include <netcdfcpp.h>
# Line 43  Rectangle::Rectangle(int n0, int n1, dou Line 49  Rectangle::Rectangle(int n0, int n1, dou
49                       double y1, int d0, int d1,                       double y1, int d0, int d1,
50                       const std::vector<double>& points,                       const std::vector<double>& points,
51                       const std::vector<int>& tags,                       const std::vector<int>& tags,
52                       const std::map<std::string, int>& tagnamestonums) :                       const simap_t& tagnamestonums) :
53      RipleyDomain(2)      RipleyDomain(2)
54  {  {
55      // ignore subdivision parameters for serial run      // ignore subdivision parameters for serial run
# Line 53  Rectangle::Rectangle(int n0, int n1, dou Line 59  Rectangle::Rectangle(int n0, int n1, dou
59      }      }
60    
61      bool warn=false;      bool warn=false;
62      // if number of subdivisions is non-positive, try to subdivide by the same      std::vector<int> factors;
63      // ratio as the number of elements      int ranks = m_mpiInfo->size;
64      if (d0<=0 && d1<=0) {      int epr[2] = {n0,n1};
65          warn=true;      int d[2] = {d0,d1};
66          d0=max(1, (int)(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));      if (d0<=0 || d1<=0) {
67          d1=m_mpiInfo->size/d0;          for (int i = 0; i < 2; i++) {
68          if (d0*d1 != m_mpiInfo->size) {              if (d[i] < 1) {
69              // ratios not the same so subdivide side with more elements only                  d[i] = 1;
70              if (n0>n1) {                  continue;
                 d0=0;  
                 d1=1;  
             } else {  
                 d0=1;  
                 d1=0;  
71              }              }
72                epr[i] = -1; // can no longer be max
73                //remove
74                if (ranks % d[i] != 0) {
75                    throw RipleyException("Invalid number of spatial subdivisions");
76                }
77                ranks /= d[i];
78            }
79            factorise(factors, ranks);
80            if (factors.size() != 0) {
81                warn = true;
82          }          }
83      }      }
84      if (d0<=0) {      
85          // d1 is preset, determine d0 - throw further down if result is no good      while (factors.size() > 0) {
86          d0=m_mpiInfo->size/d1;          int i = epr[0] > epr[1] ? 0 : 1;
87      } else if (d1<=0) {          int f = factors.back();
88          // d0 is preset, determine d1 - throw further down if result is no good          factors.pop_back();
89          d1=m_mpiInfo->size/d0;          d[i] *= f;
90            epr[i] /= f;
91      }      }
92        d0 = d[0]; d1 = d[1];
93        
94      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
95      // among number of ranks      // among number of ranks
96      if (d0*d1 != m_mpiInfo->size)      if (d0*d1 != m_mpiInfo->size)
# Line 320  void Rectangle::readBinaryGrid(escript:: Line 333  void Rectangle::readBinaryGrid(escript::
333      }      }
334  }  }
335    
336    void Rectangle::readBinaryGridFromZipped(escript::Data& out, string filename,
337                                   const ReaderParameters& params) const
338    {
339        // the mapping is not universally correct but should work on our
340        // supported platforms
341        switch (params.dataType) {
342            case DATATYPE_INT32:
343                readBinaryGridZippedImpl<int>(out, filename, params);
344                break;
345            case DATATYPE_FLOAT32:
346                readBinaryGridZippedImpl<float>(out, filename, params);
347                break;
348            case DATATYPE_FLOAT64:
349                readBinaryGridZippedImpl<double>(out, filename, params);
350                break;
351            default:
352                throw RipleyException("readBinaryGridFromZipped(): invalid or unsupported datatype");
353        }
354    }
355    
356  template<typename ValueType>  template<typename ValueType>
357  void Rectangle::readBinaryGridImpl(escript::Data& out, const string& filename,  void Rectangle::readBinaryGridImpl(escript::Data& out, const string& filename,
358                                     const ReaderParameters& params) const                                     const ReaderParameters& params) const
# Line 408  void Rectangle::readBinaryGridImpl(escri Line 441  void Rectangle::readBinaryGridImpl(escri
441      f.close();      f.close();
442  }  }
443    
444    template<typename ValueType>
445    void Rectangle::readBinaryGridZippedImpl(escript::Data& out, const string& filename,
446                                       const ReaderParameters& params) const
447    {
448        // check destination function space
449        int myN0, myN1;
450        if (out.getFunctionSpace().getTypeCode() == Nodes) {
451            myN0 = m_NN[0];
452            myN1 = m_NN[1];
453        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
454                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
455            myN0 = m_NE[0];
456            myN1 = m_NE[1];
457        } else
458            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
459    
460        // check file existence and size
461        ifstream f(filename.c_str(), ifstream::binary);
462        if (f.fail()) {
463            throw RipleyException("readBinaryGridFromZipped(): cannot open file");
464        }
465        f.seekg(0, ios::end);
466        const int numComp = out.getDataPointSize();
467        int filesize = f.tellg();
468        f.seekg(0, ios::beg);
469        std::vector<char> compressed(filesize);
470        f.read((char*)&compressed[0], filesize);
471        f.close();
472        std::vector<char> decompressed = unzip(compressed);
473        filesize = decompressed.size();
474        const int reqsize = params.numValues[0]*params.numValues[1]*numComp*sizeof(ValueType);
475        if (filesize < reqsize) {
476            throw RipleyException("readBinaryGridFromZipped(): not enough data in file");
477        }
478    
479        // check if this rank contributes anything
480        if (params.first[0] >= m_offset[0]+myN0 ||
481                params.first[0]+params.numValues[0] <= m_offset[0] ||
482                params.first[1] >= m_offset[1]+myN1 ||
483                params.first[1]+params.numValues[1] <= m_offset[1]) {
484            f.close();
485            return;
486        }
487    
488        // now determine how much this rank has to write
489    
490        // first coordinates in data object to write to
491        const int first0 = max(0, params.first[0]-m_offset[0]);
492        const int first1 = max(0, params.first[1]-m_offset[1]);
493        // indices to first value in file
494        const int idx0 = max(0, m_offset[0]-params.first[0]);
495        const int idx1 = max(0, m_offset[1]-params.first[1]);
496        // number of values to read
497        const int num0 = min(params.numValues[0]-idx0, myN0-first0);
498        const int num1 = min(params.numValues[1]-idx1, myN1-first1);
499    
500        out.requireWrite();
501        vector<ValueType> values(num0*numComp);
502        const int dpp = out.getNumDataPointsPerSample();
503    
504        for (int y=0; y<num1; y++) {
505            const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]);
506                memcpy((char*)&values[0], (char*)&decompressed[fileofs*sizeof(ValueType)], num0*numComp*sizeof(ValueType));
507            for (int x=0; x<num0; x++) {
508                const int baseIndex = first0+x*params.multiplier[0]
509                                        +(first1+y*params.multiplier[1])*myN0;
510                for (int m1=0; m1<params.multiplier[1]; m1++) {
511                    for (int m0=0; m0<params.multiplier[0]; m0++) {
512                        const int dataIndex = baseIndex+m0+m1*myN0;
513                        double* dest = out.getSampleDataRW(dataIndex);
514                        for (int c=0; c<numComp; c++) {
515                            ValueType val = values[x*numComp+c];
516    
517                            if (params.byteOrder != BYTEORDER_NATIVE) {
518                                char* cval = reinterpret_cast<char*>(&val);
519                                // this will alter val!!
520                                byte_swap32(cval);
521                            }
522                            if (!std::isnan(val)) {
523                                for (int q=0; q<dpp; q++) {
524                                    *dest++ = static_cast<double>(val);
525                                }
526                            }
527                        }
528                    }
529                }
530            }
531        }
532    
533        f.close();
534    }
535    
536  void Rectangle::writeBinaryGrid(const escript::Data& in, string filename,  void Rectangle::writeBinaryGrid(const escript::Data& in, string filename,
537                                  int byteOrder, int dataType) const                                  int byteOrder, int dataType) const
538  {  {
# Line 455  void Rectangle::writeBinaryGridImpl(cons Line 580  void Rectangle::writeBinaryGridImpl(cons
580      if (numComp > 1 || dpp > 1)      if (numComp > 1 || dpp > 1)
581          throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");          throw RipleyException("writeBinaryGrid(): only scalar, single-value data supported");
582    
     escript::Data* _in = const_cast<escript::Data*>(&in);  
583      const int fileSize = sizeof(ValueType)*numComp*dpp*totalN0*totalN1;      const int fileSize = sizeof(ValueType)*numComp*dpp*totalN0*totalN1;
584    
585      // from here on we know that each sample consists of one value      // from here on we know that each sample consists of one value
# Line 468  void Rectangle::writeBinaryGridImpl(cons Line 592  void Rectangle::writeBinaryGridImpl(cons
592          ostringstream oss;          ostringstream oss;
593    
594          for (index_t x=0; x<myN0; x++) {          for (index_t x=0; x<myN0; x++) {
595              const double* sample = _in->getSampleDataRO(y*myN0+x);              const double* sample = in.getSampleDataRO(y*myN0+x);
596              ValueType fvalue = static_cast<ValueType>(*sample);              ValueType fvalue = static_cast<ValueType>(*sample);
597              if (byteOrder == BYTEORDER_NATIVE) {              if (byteOrder == BYTEORDER_NATIVE) {
598                  oss.write((char*)&fvalue, sizeof(fvalue));                  oss.write((char*)&fvalue, sizeof(fvalue));
# Line 639  const int* Rectangle::borrowSampleRefere Line 763  const int* Rectangle::borrowSampleRefere
763          case FaceElements:          case FaceElements:
764          case ReducedFaceElements:          case ReducedFaceElements:
765              return &m_faceId[0];              return &m_faceId[0];
766            case Points:
767                return &m_diracPointNodeIDs[0];
768          default:          default:
769              break;              break;
770      }      }
# Line 896  void Rectangle::assembleCoordinates(escr Line 1022  void Rectangle::assembleCoordinates(escr
1022  }  }
1023    
1024  //protected  //protected
1025  void Rectangle::assembleGradient(escript::Data& out, escript::Data& in) const  void Rectangle::assembleGradient(escript::Data& out, const escript::Data& in) const
1026  {  {
1027      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1028      const double cx0 = .21132486540518711775/m_dx[0];      const double cx0 = .21132486540518711775/m_dx[0];
# Line 1101  void Rectangle::assembleGradient(escript Line 1227  void Rectangle::assembleGradient(escript
1227  }  }
1228    
1229  //protected  //protected
1230  void Rectangle::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const  void Rectangle::assembleIntegrate(vector<double>& integrals,
1231                                      const escript::Data& arg) const
1232  {  {
1233      const dim_t numComp = arg.getDataPointSize();      const dim_t numComp = arg.getDataPointSize();
1234      const index_t left = (m_offset[0]==0 ? 0 : 1);      const index_t left = (m_offset[0]==0 ? 0 : 1);
# Line 1287  dim_t Rectangle::insertNeighbourNodes(In Line 1414  dim_t Rectangle::insertNeighbourNodes(In
1414  }  }
1415    
1416  //protected  //protected
1417  void Rectangle::nodesToDOF(escript::Data& out, escript::Data& in) const  void Rectangle::nodesToDOF(escript::Data& out, const escript::Data& in) const
1418  {  {
1419      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1420      out.requireWrite();      out.requireWrite();
# Line 1307  void Rectangle::nodesToDOF(escript::Data Line 1434  void Rectangle::nodesToDOF(escript::Data
1434  }  }
1435    
1436  //protected  //protected
1437  void Rectangle::dofToNodes(escript::Data& out, escript::Data& in) const  void Rectangle::dofToNodes(escript::Data& out, const escript::Data& in) const
1438  {  {
1439      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1440      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);      Paso_Coupler* coupler = Paso_Coupler_alloc(m_connector, numComp);
1441      in.requireWrite();      // expand data object if necessary to be able to grab the whole data
1442      Paso_Coupler_startCollect(coupler, in.getSampleDataRW(0));      const_cast<escript::Data*>(&in)->expand();
1443        Paso_Coupler_startCollect(coupler, in.getSampleDataRO(0));
1444    
1445      const dim_t numDOF = getNumDOF();      const dim_t numDOF = getNumDOF();
1446      out.requireWrite();      out.requireWrite();
# Line 1507  void Rectangle::createPattern() Line 1635  void Rectangle::createPattern()
1635      RankVector neighbour;      RankVector neighbour;
1636      IndexVector offsetInShared(1,0);      IndexVector offsetInShared(1,0);
1637      IndexVector sendShared, recvShared;      IndexVector sendShared, recvShared;
1638      int numShared=0;      int numShared=0, expectedShared=0;
1639      const int x=m_mpiInfo->rank%m_NX[0];      const int x=m_mpiInfo->rank%m_NX[0];
1640      const int y=m_mpiInfo->rank/m_NX[0];      const int y=m_mpiInfo->rank/m_NX[0];
1641        if (x > 0)
1642            expectedShared += nDOF1;
1643        if (x < m_NX[0] - 1)
1644            expectedShared += nDOF1;
1645        if (y > 0)
1646            expectedShared += nDOF0;
1647        if (y < m_NX[1] - 1)
1648            expectedShared += nDOF0;
1649        if (x > 0 && y > 0) expectedShared++;
1650        if (x > 0 && y < m_NX[1] - 1) expectedShared++;
1651        if (x < m_NX[0] - 1 && y > 0) expectedShared++;
1652        if (x < m_NX[0] - 1 && y < m_NX[1] - 1) expectedShared++;
1653        
1654        vector<IndexVector> rowIndices(expectedShared);
1655        
1656      for (int i1=-1; i1<2; i1++) {      for (int i1=-1; i1<2; i1++) {
1657          for (int i0=-1; i0<2; i0++) {          for (int i0=-1; i0<2; i0++) {
1658              // skip this rank              // skip this rank
# Line 1529  void Rectangle::createPattern() Line 1672  void Rectangle::createPattern()
1672                          sendShared.push_back(firstDOF+i);                          sendShared.push_back(firstDOF+i);
1673                          recvShared.push_back(numDOF+numShared);                          recvShared.push_back(numDOF+numShared);
1674                          if (i>0)                          if (i>0)
1675                              colIndices[firstDOF+i-1].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i-1, numShared);
1676                          colIndices[firstDOF+i].push_back(numShared);                          doublyLink(colIndices, rowIndices, firstDOF+i, numShared);
1677                          if (i<nDOF0-1)                          if (i<nDOF0-1)
1678                              colIndices[firstDOF+i+1].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+i+1, numShared);
1679                          m_dofMap[firstNode+i]=numDOF+numShared;                          m_dofMap[firstNode+i]=numDOF+numShared;
1680                      }                      }
1681                  } else if (i1==0) {                  } else if (i1==0) {
# Line 1544  void Rectangle::createPattern() Line 1687  void Rectangle::createPattern()
1687                          sendShared.push_back(firstDOF+i*nDOF0);                          sendShared.push_back(firstDOF+i*nDOF0);
1688                          recvShared.push_back(numDOF+numShared);                          recvShared.push_back(numDOF+numShared);
1689                          if (i>0)                          if (i>0)
1690                              colIndices[firstDOF+(i-1)*nDOF0].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+(i-1)*nDOF0, numShared);
1691                          colIndices[firstDOF+i*nDOF0].push_back(numShared);                          doublyLink(colIndices, rowIndices, firstDOF+i*nDOF0, numShared);
1692                          if (i<nDOF1-1)                          if (i<nDOF1-1)
1693                              colIndices[firstDOF+(i+1)*nDOF0].push_back(numShared);                              doublyLink(colIndices, rowIndices, firstDOF+(i+1)*nDOF0, numShared);
1694                          m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;                          m_dofMap[firstNode+i*m_NN[0]]=numDOF+numShared;
1695                      }                      }
1696                  } else {                  } else {
# Line 1557  void Rectangle::createPattern() Line 1700  void Rectangle::createPattern()
1700                      offsetInShared.push_back(offsetInShared.back()+1);                      offsetInShared.push_back(offsetInShared.back()+1);
1701                      sendShared.push_back(dof);                      sendShared.push_back(dof);
1702                      recvShared.push_back(numDOF+numShared);                      recvShared.push_back(numDOF+numShared);
1703                      colIndices[dof].push_back(numShared);                      doublyLink(colIndices, rowIndices, dof, numShared);
1704                      m_dofMap[node]=numDOF+numShared;                      m_dofMap[node]=numDOF+numShared;
1705                      ++numShared;                      ++numShared;
1706                  }                  }
1707              }              }
1708          }          }
1709      }      }
1710            
1711    #pragma omp parallel for
1712        for (int i = 0; i < numShared; i++) {
1713            std::sort(rowIndices[i].begin(), rowIndices[i].end());
1714        }
1715        
1716      // create connector      // create connector
1717      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
1718              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],              numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
# Line 1579  void Rectangle::createPattern() Line 1727  void Rectangle::createPattern()
1727      // create main and couple blocks      // create main and couple blocks
1728      Paso_Pattern *mainPattern = createMainPattern();      Paso_Pattern *mainPattern = createMainPattern();
1729      Paso_Pattern *colPattern, *rowPattern;      Paso_Pattern *colPattern, *rowPattern;
1730      createCouplePatterns(colIndices, numShared, &colPattern, &rowPattern);      createCouplePatterns(colIndices, rowIndices, numShared, &colPattern, &rowPattern);
1731    
1732      // allocate paso distribution      // allocate paso distribution
1733      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,      Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
# Line 1682  void Rectangle::addToMatrixAndRHS(Paso_S Line 1830  void Rectangle::addToMatrixAndRHS(Paso_S
1830    
1831  //protected  //protected
1832  void Rectangle::interpolateNodesOnElements(escript::Data& out,  void Rectangle::interpolateNodesOnElements(escript::Data& out,
1833                                          escript::Data& in, bool reduced) const                                             const escript::Data& in,
1834                                               bool reduced) const
1835  {  {
1836      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
1837      if (reduced) {      if (reduced) {
# Line 1740  void Rectangle::interpolateNodesOnElemen Line 1889  void Rectangle::interpolateNodesOnElemen
1889  }  }
1890    
1891  //protected  //protected
1892  void Rectangle::interpolateNodesOnFaces(escript::Data& out, escript::Data& in,  void Rectangle::interpolateNodesOnFaces(escript::Data& out,
1893                                            const escript::Data& in,
1894                                          bool reduced) const                                          bool reduced) const
1895  {  {
1896      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
# Line 1864  void Rectangle::interpolateNodesOnFaces( Line 2014  void Rectangle::interpolateNodesOnFaces(
2014  namespace  namespace
2015  {  {
2016      // Calculates a guassian blur colvolution matrix for 2D      // Calculates a guassian blur colvolution matrix for 2D
2017        // See wiki article on the subject    
2018      double* get2DGauss(unsigned radius, double sigma)      double* get2DGauss(unsigned radius, double sigma)
2019      {      {
2020          double* arr=new double[(radius*2+1)*(radius*2+1)];          double* arr=new double[(radius*2+1)*(radius*2+1)];
2021          double common=M_1_PI*0.5*1/(sigma*sigma);          double common=M_1_PI*0.5*1/(sigma*sigma);
2022      double total=0;          double total=0;
2023      int r=static_cast<int>(radius);          int r=static_cast<int>(radius);
2024      for (int y=-r;y<=r;++y)          for (int y=-r;y<=r;++y)
2025      {          {
2026          for (int x=-r;x<=r;++x)              for (int x=-r;x<=r;++x)
2027          {                      {        
2028              arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y)/(2*sigma*sigma));                  arr[(x+r)+(y+r)*(r*2+1)]=common*exp(-(x*x+y*y)/(2*sigma*sigma));
2029  // cout << (x+y*(r*2+1)) << " " << arr[(x+r)+(y+r)*(r*2+1)] << endl;                  total+=arr[(x+r)+(y+r)*(r*2+1)];
2030              total+=arr[(x+r)+(y+r)*(r*2+1)];              }
2031          }          }
2032      }          double invtotal=1/total;
2033      double invtotal=1/total;          for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)
2034  //cout << "Inv total is "    << invtotal << endl;          {
2035      for (size_t p=0;p<(radius*2+1)*(radius*2+1);++p)              arr[p]*=invtotal;
2036      {          }
2037          arr[p]*=invtotal;          return arr;
 //cout << p << "->" << arr[p] << endl;        
     }  
     return arr;  
2038      }      }
2039            
2040      // applies conv to source to get a point.      // applies conv to source to get a point.
# Line 1894  namespace Line 2042  namespace
2042      double Convolve2D(double* conv, double* source, size_t xp, size_t yp, unsigned radius, size_t width)      double Convolve2D(double* conv, double* source, size_t xp, size_t yp, unsigned radius, size_t width)
2043      {      {
2044          size_t bx=xp-radius, by=yp-radius;          size_t bx=xp-radius, by=yp-radius;
2045      size_t sbase=bx+by*width;          size_t sbase=bx+by*width;
2046      double result=0;          double result=0;
2047      for (int y=0;y<2*radius+1;++y)          for (int y=0;y<2*radius+1;++y)
2048      {              {        
2049          for (int x=0;x<2*radius+1;++x)              for (int x=0;x<2*radius+1;++x)
2050          {              {
2051              result+=conv[x+y*(2*radius+1)] * source[sbase + x+y*width];                  result+=conv[x+y*(2*radius+1)] * source[sbase + x+y*width];
2052          }              }
2053      }          }
2054          return result;                return result;      
2055      }      }
2056  }  }
2057    
2058    
2059    /* This is a wrapper for filtered (and non-filtered) randoms
2060     * For detailed doco see randomFillWorker
2061    */
2062    escript::Data Rectangle::randomFill(const escript::DataTypes::ShapeType& shape,
2063           const escript::FunctionSpace& what,
2064           long seed, const boost::python::tuple& filter) const
2065    {
2066        int numvals=escript::DataTypes::noValues(shape);
2067        if (len(filter)>0 && (numvals!=1))
2068        {
2069            throw RipleyException("Ripley only supports filters for scalar data.");
2070        }
2071        escript::Data res=randomFillWorker(shape, seed, filter);
2072        if (res.getFunctionSpace()!=what)
2073        {
2074            escript::Data r=escript::Data(res, what);
2075            return r;
2076        }
2077        return res;
2078    }
2079    
2080    
2081  /* This routine produces a Data object filled with smoothed random data.  /* This routine produces a Data object filled with smoothed random data.
2082  The dimensions of the rectangle being filled are internal[0] x internal[1] points.  The dimensions of the rectangle being filled are internal[0] x internal[1] points.
2083  A parameter radius  dives the size of the stencil used for the smoothing.  A parameter radius  gives the size of the stencil used for the smoothing.
2084  A point on the left hand edge for example, will still require `radius` extra points to the left  A point on the left hand edge for example, will still require `radius` extra points to the left
2085  in order to complete the stencil.  in order to complete the stencil.
2086    
# Line 1948  that ripley has. Line 2117  that ripley has.
2117    
2118    
2119  */  */
2120  escript::Data Rectangle::randomFill(long seed, const boost::python::tuple& filter) const  escript::Data Rectangle::randomFillWorker(const escript::DataTypes::ShapeType& shape,
2121           long seed, const boost::python::tuple& filter) const
2122  {  {
2123      if (m_numDim!=2)      if (m_numDim!=2)
2124      {      {
2125          throw RipleyException("Only 2D supported at this time.");          throw RipleyException("Only 2D supported at this time.");
2126      }      }
2127      if (len(filter)!=3) {  
2128          throw RipleyException("Unsupported random filter");      unsigned int radius=0;      // these are only used by gaussian
2129      }      double sigma=0.5;
2130      boost::python::extract<string> ex(filter[0]);      
2131      if (!ex.check() || (ex()!="gaussian"))      unsigned int numvals=escript::DataTypes::noValues(shape);
2132            
2133        
2134        if (len(filter)==0)
2135      {      {
2136          throw RipleyException("Unsupported random filter");          // nothing special required here yet
2137      }      }    
2138      boost::python::extract<unsigned int> ex1(filter[1]);      else if (len(filter)==3)
     if (!ex1.check())  
2139      {      {
2140          throw RipleyException("Radius of gaussian filter must be a positive integer.");          boost::python::extract<string> ex(filter[0]);
2141            if (!ex.check() || (ex()!="gaussian"))
2142            {
2143                throw RipleyException("Unsupported random filter");
2144            }
2145            boost::python::extract<unsigned int> ex1(filter[1]);
2146            if (!ex1.check())
2147            {
2148                throw RipleyException("Radius of gaussian filter must be a positive integer.");
2149            }
2150            radius=ex1();
2151            sigma=0.5;
2152            boost::python::extract<double> ex2(filter[2]);
2153            if (!ex2.check() || (sigma=ex2())<=0)
2154            {
2155                throw RipleyException("Sigma must be a postive floating point number.");
2156            }        
2157      }      }
2158      unsigned int radius=ex1();      else
     double sigma=0.5;  
     boost::python::extract<double> ex2(filter[2]);  
     if (!ex2.check() || (sigma=ex2())<=0)  
2159      {      {
2160          throw RipleyException("Sigma must be a postive floating point number.");          throw RipleyException("Unsupported random filter for Rectangle.");
2161      }          }
2162          
2163      
2164            
2165      size_t internal[2];      size_t internal[2];
2166      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
2167      internal[1]=m_NE[1]+1;  // that is, the ones we need smoothed versions of      internal[1]=m_NE[1]+1;      // that is, the ones we need smoothed versions of
2168      size_t ext[2];      size_t ext[2];
2169      ext[0]=internal[0]+2*radius;    // includes points we need as input      ext[0]=internal[0]+2*radius;        // includes points we need as input
2170      ext[1]=internal[1]+2*radius;    // for smoothing      ext[1]=internal[1]+2*radius;        // for smoothing
2171            
2172      // now we check to see if the radius is acceptable      // now we check to see if the radius is acceptable
2173      // That is, would not cross multiple ranks in MPI      // That is, would not cross multiple ranks in MPI
2174    
2175      if ((2*radius>=internal[0]) || (2*radius>=internal[1]))      if (2*radius>=internal[0]-4)
     {  
         throw RipleyException("Radius of gaussian filter must be less than half the width/height of a rank");  
     }  
       
     
     double* src=new double[ext[0]*ext[1]];  
     esysUtils::randomFillArray(seed, src, ext[0]*ext[1]);    
       
     
 #ifdef ESYS_MPI      
     size_t inset=2*radius+1;  
     size_t Eheight=ext[1]-2*inset;  // how high the E (shared) region is  
     size_t Swidth=ext[0]-2*inset;  
   
     double* SWin=new double[inset*inset];  memset(SWin, 0, inset*inset*sizeof(double));  
     double* SEin=new double[inset*inset];  memset(SEin, 0, inset*inset*sizeof(double));  
     double* NWin=new double[inset*inset];  memset(NWin, 0, inset*inset*sizeof(double));  
     double* Sin=new double[inset*Swidth];  memset(Sin, 0, inset*Swidth*sizeof(double));  
     double* Win=new double[inset*Eheight];  memset(Win, 0, inset*Eheight*sizeof(double));  
   
     double* NEout=new double[inset*inset];  memset(NEout, 0, inset*inset*sizeof(double));  
     unsigned int base=ext[0]-inset+(ext[1]-inset)*ext[0];  
     for (unsigned int i=0;i<inset;++i)  
2176      {      {
2177      memcpy(NEout+inset*i, src+base, inset*sizeof(double));          throw RipleyException("Radius of gaussian filter is too large for X dimension of a rank");
     base+=ext[0];  
2178      }      }
2179      double* NWout=new double[inset*inset];  memset(NWout, 0, inset*inset*sizeof(double));      if (2*radius>=internal[1]-4)
     base=(ext[1]-inset)*ext[0];  
     for (unsigned int i=0;i<inset;++i)  
2180      {      {
2181      memcpy(NWout+inset*i, src+base, inset*sizeof(double));          throw RipleyException("Radius of gaussian filter is too large for Y dimension of a rank");
     base+=ext[0];  
2182      }      }
2183    
2184        double* src=new double[ext[0]*ext[1]*numvals];
2185        esysUtils::randomFillArray(seed, src, ext[0]*ext[1]*numvals);  
2186            
2187      double* SEout=new double[inset*inset];  memset(SEout, 0, inset*inset*sizeof(double));  
2188      base=ext[0]-inset;  
2189      for (int i=0;i<inset;++i)  #ifdef ESYS_MPI
2190      {      if ((internal[0]<5) || (internal[1]<5))
     memcpy(SEout+inset*i, src+base, inset*sizeof(double));  
     base+=ext[0];  
     }  
     double* Nout=new double[inset*Swidth];  memset(Nout, 0, inset*Swidth*sizeof(double));  
     base=inset+(ext[1]-inset)*ext[0];  
     for (unsigned int i=0;i<inset;++i)  
2191      {      {
2192      memcpy(Nout+Swidth*i, src+base, Swidth*sizeof(double));          // since the dimensions are equal for all ranks, this exception
2193      base+=ext[0];          // will be thrown on all ranks
2194            throw RipleyException("Random Data in Ripley requries at least five elements per side per rank.");
2195      }      }
2196        dim_t X=m_mpiInfo->rank%m_NX[0];
2197        dim_t Y=m_mpiInfo->rank%(m_NX[0]*m_NX[1])/m_NX[0];
2198    #endif      
2199    
2200    /*    
2201        // if we wanted to test a repeating pattern
2202        size_t basex=0;
2203        size_t basey=0;
2204    #ifdef ESYS_MPI    
2205        basex=X*m_gNE[0]/m_NX[0];
2206        basey=Y*m_gNE[1]/m_NX[1];
2207    #endif
2208            
2209        esysUtils::patternFillArray2D(ext[0], ext[1], src, 4, basex, basey, numvals);
2210    */
2211    
2212            
2213      double* Eout=new double[inset*Eheight];  memset(Eout, 0, inset*Eheight*sizeof(double));  #ifdef ESYS_MPI  
2214      base=ext[0]-inset+inset*ext[0];      
2215      for (unsigned int i=0;i<Eheight;++i)      BlockGrid2 grid(m_NX[0]-1, m_NX[1]-1);
2216      {      size_t inset=2*radius+2;    // Its +2 not +1 because a whole element is shared (and hence
2217      memcpy(Eout+i*inset, src+base, inset*sizeof(double));                  // there is an overlap of two points both of which need to have "radius" points on either side.
2218      base+=ext[0];      
2219      }        size_t xmidlen=ext[0]-2*inset;      // how wide is the x-dimension between the two insets
2220        size_t ymidlen=ext[1]-2*inset;      
2221        
2222        Block2 block(ext[0], ext[1], inset, xmidlen, ymidlen, numvals);
2223    
2224      MPI_Request reqs[10];      MPI_Request reqs[40];               // a non-tight upper bound on how many we need
2225      MPI_Status stats[10];      MPI_Status stats[40];
2226      short rused=0;      short rused=0;
2227            
2228      dim_t X=m_mpiInfo->rank%m_NX[0];      messvec incoms;
2229      dim_t Y=m_mpiInfo->rank/m_NX[0];      messvec outcoms;  
     dim_t row=m_NX[0];  
     bool swused=false;      // These vars will be true if data needs to be copied out of them  
     bool seused=false;  
     bool nwused=false;  
     bool sused=false;  
     bool wused=false;      
       
     // Tags:  
     // 10 : EW transfer (middle)  
     // 8 : NS transfer (middle)  
     // 7 : NE corner -> to N, E and NE  
     // 11 : NW corner to SW corner (only used on the left hand edge  
     // 12 : SE corner to SW corner (only used on the bottom edge  
2230            
2231        grid.generateInNeighbours(X, Y, incoms);
2232        grid.generateOutNeighbours(X, Y, outcoms);
     int comserr=0;  
     if (Y!=0)   // not on bottom row,  
     {  
     if (X!=0)   // not on the left hand edge  
     {  
         // recv bottomleft from SW  
         comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, (X-1)+(Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));  
         swused=true;  
         comserr|=MPI_Irecv(Win, Eheight*inset, MPI_DOUBLE, X-1+Y*row, 10, m_mpiInfo->comm, reqs+(rused++));  
         wused=true;  
     }  
     else    // on the left hand edge  
     {  
         comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, (Y-1)*row, 11, m_mpiInfo->comm, reqs+(rused++));  
         swused=true;  
     }  
     comserr|=MPI_Irecv(Sin, Swidth*inset, MPI_DOUBLE, X+(Y-1)*row, 8, m_mpiInfo->comm, reqs+(rused++));  
     sused=true;  
     comserr|=MPI_Irecv(SEin, inset*inset, MPI_DOUBLE, X+(Y-1)*row, 7, m_mpiInfo->comm, reqs+(rused++));  
     seused=true;  
   
         
     }  
     else        // on the bottom row  
     {  
     if (X!=0)  
     {  
         comserr|=MPI_Irecv(Win, Eheight*inset, MPI_DOUBLE, X-1+Y*row, 10, m_mpiInfo->comm, reqs+(rused++));  
         wused=true;  
         // Need to use tag 12 here because SW is coming from the East not South East  
         comserr|=MPI_Irecv(SWin, inset*inset, MPI_DOUBLE, X-1+Y*row, 12, m_mpiInfo->comm, reqs+(rused++));  
         swused=true;  
     }  
     if (X!=(row-1))  
     {  
         comserr|=MPI_Isend(SEout, inset*inset, MPI_DOUBLE, X+1+(Y)*row, 12, m_mpiInfo->comm, reqs+(rused++));    
     }  
     }  
2233            
2234      if (Y!=(m_NX[1]-1)) // not on the top row      block.copyAllToBuffer(src);  
2235      {      
2236      comserr|=MPI_Isend(Nout, inset*Swidth, MPI_DOUBLE, X+(Y+1)*row, 8, m_mpiInfo->comm, reqs+(rused++));      int comserr=0;    
2237      comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+(Y+1)*row, 7, m_mpiInfo->comm, reqs+(rused++));      for (size_t i=0;i<incoms.size();++i)
     if (X!=(row-1)) // not on right hand edge  
     {  
         comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+1+(Y+1)*row, 7, m_mpiInfo->comm, reqs+(rused++));  
     }  
     if (X==0)   // left hand edge  
     {  
         comserr|=MPI_Isend(NWout, inset*inset, MPI_DOUBLE, (Y+1)*row,11, m_mpiInfo->comm, reqs+(rused++));        
     }    
     }  
     if (X!=(row-1)) // not on right hand edge  
2238      {      {
2239      comserr|=MPI_Isend(NEout, inset*inset, MPI_DOUBLE, X+1+(Y)*row, 7, m_mpiInfo->comm, reqs+(rused++));          message& m=incoms[i];
2240      comserr|=MPI_Isend(Eout, Eheight*inset, MPI_DOUBLE, X+1+(Y)*row, 10, m_mpiInfo->comm, reqs+(rused++));          comserr|=MPI_Irecv(block.getInBuffer(m.destbuffid), block.getBuffSize(m.destbuffid) , MPI_DOUBLE, m.sourceID, m.tag, m_mpiInfo->comm, reqs+(rused++));
2241            block.setUsed(m.destbuffid);
2242      }      }
2243      if (X!=0)  
2244        for (size_t i=0;i<outcoms.size();++i)
2245      {      {
2246      comserr|=MPI_Irecv(NWin, inset*inset, MPI_DOUBLE, (X-1)+Y*row, 7, m_mpiInfo->comm, reqs+(rused++));          message& m=outcoms[i];
2247      nwused=true;          comserr|=MPI_Isend(block.getOutBuffer(m.srcbuffid), block.getBuffSize(m.srcbuffid) , MPI_DOUBLE, m.destID, m.tag, m_mpiInfo->comm, reqs+(rused++));
2248              }    
         
     }  
2249            
2250      if (!comserr)      if (!comserr)
2251      {          {    
# Line 2138  escript::Data Rectangle::randomFill(long Line 2254  escript::Data Rectangle::randomFill(long
2254    
2255      if (comserr)      if (comserr)
2256      {      {
2257      // Yes this is throwing an exception as a result of an MPI error.          // Yes this is throwing an exception as a result of an MPI error.
2258      // and no we don't inform the other ranks that we are doing this.          // and no we don't inform the other ranks that we are doing this.
2259      // however, we have no reason to believe coms work at this point anyway          // however, we have no reason to believe coms work at this point anyway
2260          throw RipleyException("Error in coms for randomFill");                throw RipleyException("Error in coms for randomFill");      
2261      }      }
2262            
2263        block.copyUsedFromBuffer(src);    
     // Need to copy the values back from the buffers  
     // Copy from SW  
2264            
2265      if (swused)  #endif    
2266      {      
2267      base=0;      if (radius==0 || numvals>1) // the truth of either should imply the truth of the other but let's be safe
     for (unsigned int i=0;i<inset;++i)  
     {  
         memcpy(src+base, SWin+i*inset, inset*sizeof(double));  
         base+=ext[0];  
     }  
     }  
     if (seused)  
     {  
         base=ext[0]-inset;  
     for (unsigned int i=0;i<inset;++i)  
     {  
         memcpy(src+base, SEin+i*inset, inset*sizeof(double));  
         base+=ext[0];  
     }  
     }  
     if (nwused)  
     {  
         base=(ext[1]-inset)*ext[0];  
     for (unsigned int i=0;i<inset;++i)  
     {  
         memcpy(src+base, NWin+i*inset, inset*sizeof(double));  
         base+=ext[0];  
     }  
     }  
     if (sused)  
     {  
        base=inset;  
        for (unsigned int i=0;i<inset;++i)  
        {  
        memcpy(src+base, Sin+i*Swidth, Swidth*sizeof(double));  
        base+=ext[0];  
        }  
     }  
     if (wused)  
2268      {      {
     base=inset*ext[0];  
     for (unsigned int i=0;i<Eheight;++i)  
     {  
         memcpy(src+base, Win+i*inset, inset*sizeof(double));  
         base+=ext[0];  
     }  
2269                
2270            escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
2271            escript::Data resdat(0, shape, fs , true);
2272            // don't need to check for exwrite because we just made it
2273            escript::DataVector& dv=resdat.getExpandedVectorReference();
2274            
2275            
2276            // now we need to copy values over
2277            for (size_t y=0;y<(internal[1]);++y)    
2278            {
2279                for (size_t x=0;x<(internal[0]);++x)
2280                {
2281                    for (unsigned int i=0;i<numvals;++i)
2282                    {
2283                        dv[i+(x+y*(internal[0]))*numvals]=src[i+(x+y*ext[0])*numvals];
2284                    }
2285                }
2286            }
2287            delete[] src;
2288            return resdat;      
2289        }
2290        else                // filter enabled      
2291        {    
2292            escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());
2293            escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);
2294            // don't need to check for exwrite because we just made it
2295            escript::DataVector& dv=resdat.getExpandedVectorReference();
2296            double* convolution=get2DGauss(radius, sigma);
2297            for (size_t y=0;y<(internal[1]);++y)    
2298            {
2299                for (size_t x=0;x<(internal[0]);++x)
2300                {    
2301                    dv[x+y*(internal[0])]=Convolve2D(convolution, src, x+radius, y+radius, radius, ext[0]);
2302                    
2303                }
2304            }
2305            delete[] convolution;
2306            delete[] src;
2307            return resdat;
2308      }      }
       
     delete[] SWin;  
     delete[] SEin;  
     delete[] NWin;  
     delete[] Sin;  
     delete[] Win;  
   
     delete[] NEout;  
     delete[] NWout;  
     delete[] SEout;  
     delete[] Nout;  
     delete[] Eout;  
 #endif      
     escript::FunctionSpace fs(getPtr(), getContinuousFunctionCode());  
     escript::Data resdat(0, escript::DataTypes::scalarShape, fs , true);  
     // don't need to check for exwrite because we just made it  
     escript::DataVector& dv=resdat.getExpandedVectorReference();  
     double* convolution=get2DGauss(radius, sigma);  
     for (size_t y=0;y<(internal[1]);++y)      
     {  
         for (size_t x=0;x<(internal[0]);++x)  
     {      
         dv[x+y*(internal[0])]=Convolve2D(convolution, src, x+radius, y+radius, radius, ext[0]);  
           
     }  
     }  
     delete[] convolution;  
     delete[] src;  
     return resdat;  
2309  }  }
2310    
2311  int Rectangle::findNode(double *coords) const {  int Rectangle::findNode(const double *coords) const {
2312      const int NOT_MINE = -1;      const int NOT_MINE = -1;
2313      //is the found element even owned by this rank      //is the found element even owned by this rank
2314        // (inside owned or shared elements but will map to an owned element)
2315      for (int dim = 0; dim < m_numDim; dim++) {      for (int dim = 0; dim < m_numDim; dim++) {
2316          if (m_origin[dim] + m_offset[dim] > coords[dim]  || m_origin[dim]          double min = m_origin[dim] + m_offset[dim]* m_dx[dim]
2317                  + m_offset[dim] + m_dx[dim]*m_ownNE[dim] < coords[dim]) {                  - m_dx[dim]/2.; //allows for point outside mapping onto node
2318            double max = m_origin[dim] + (m_offset[dim] + m_NE[dim])*m_dx[dim]
2319                    + m_dx[dim]/2.;
2320            if (min > coords[dim] || max < coords[dim]) {
2321              return NOT_MINE;              return NOT_MINE;
2322          }          }
2323      }      }
# Line 2247  int Rectangle::findNode(double *coords) Line 2334  int Rectangle::findNode(double *coords)
2334          minDist += m_dx[dim]*m_dx[dim];          minDist += m_dx[dim]*m_dx[dim];
2335      }      }
2336      //find the closest node      //find the closest node
2337      for (int dx = 0; dx < 2; dx++) {      for (int dx = 0; dx < 1; dx++) {
2338          double xdist = (x - (ex + dx)*m_dx[0]);          double xdist = (x - (ex + dx)*m_dx[0]);
2339          for (int dy = 0; dy < 2; dy++) {          for (int dy = 0; dy < 1; dy++) {
2340              double ydist = (y - (ey + dy)*m_dx[1]);              double ydist = (y - (ey + dy)*m_dx[1]);
2341              double total = xdist*xdist + ydist*ydist;              double total = xdist*xdist + ydist*ydist;
2342              if (total < minDist) {              if (total < minDist) {
2343                  closest = INDEX2(ex+dx, ey+dy, m_NE[0] + 1);                  closest = INDEX2(ex+dx-m_offset[0], ey+dy-m_offset[1], m_NE[0] + 1);
2344                  minDist = total;                  minDist = total;
2345              }              }
2346          }          }
2347      }      }
2348        //if this happens, we've let a dirac point slip through, which is awful
2349        if (closest == NOT_MINE) {
2350            throw RipleyException("Unable to map appropriate dirac point to a node,"
2351                    " implementation problem in Rectangle::findNode()");
2352        }
2353      return closest;      return closest;
2354  }  }
2355    
2356  void Rectangle::setAssembler(std::string type, std::map<std::string,  void Rectangle::setAssembler(std::string type, std::map<std::string,
2357          escript::Data> constants) {          escript::Data> constants) {
2358      if (type.compare("WaveAssembler") == 0) {      if (type.compare("WaveAssembler") == 0) {
2359            if (assembler_type != WAVE_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
2360                throw RipleyException("Domain already using a different custom assembler");
2361            assembler_type = WAVE_ASSEMBLER;
2362          delete assembler;          delete assembler;
2363          assembler = new WaveAssembler2D(this, m_dx, m_NX, m_NE, m_NN, constants);          assembler = new WaveAssembler2D(this, m_dx, m_NX, m_NE, m_NN, constants);
2364        } else if (type.compare("LameAssembler") == 0) {
2365            if (assembler_type != LAME_ASSEMBLER && assembler_type != DEFAULT_ASSEMBLER)
2366                throw RipleyException("Domain already using a different custom assembler");
2367            assembler_type = LAME_ASSEMBLER;
2368            delete assembler;
2369            assembler = new LameAssembler2D(this, m_dx, m_NX, m_NE, m_NN);
2370      } else { //else ifs would go before this for other types      } else { //else ifs would go before this for other types
2371          throw RipleyException("Ripley::Rectangle does not support the"          throw RipleyException("Ripley::Rectangle does not support the"
2372                                  " requested assembler");                                  " requested assembler");

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

  ViewVC Help
Powered by ViewVC 1.1.26