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

revision 3943 by caltinay, Thu Aug 9 04:43:24 2012 UTC revision 3971 by caltinay, Wed Sep 19 02:55:35 2012 UTC
# Line 51  Brick::Brick(int n0, int n1, int n2, dou Line 51  Brick::Brick(int n0, int n1, int n2, dou
51      // ratio as the number of elements      // ratio as the number of elements
52      if (d0<=0 && d1<=0 && d2<=0) {      if (d0<=0 && d1<=0 && d2<=0) {
53          warn=true;          warn=true;
54          d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1./3));          d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1.f/3));
55          d1=(int)(d0*n1/(float)n0);          d1=(int)(d0*n1/(float)n0);
56          d2=m_mpiInfo->size/(d0*d1);          d2=m_mpiInfo->size/(d0*d1);
57          if (d0*d1*d2 != m_mpiInfo->size) {          if (d0*d1*d2 != m_mpiInfo->size) {
# Line 83  Brick::Brick(int n0, int n1, int n2, dou Line 83  Brick::Brick(int n0, int n1, int n2, dou
83          if (d0*d1*d2 != m_mpiInfo->size) {          if (d0*d1*d2 != m_mpiInfo->size) {
84              // ratios not the same so subdivide side with more elements only              // ratios not the same so subdivide side with more elements only
85              if (n0>n1) {              if (n0>n1) {
86                    d0=0;
87                  d1=1;                  d1=1;
88              } else {              } else {
89                  d0=1;                  d0=1;
90                    d1=0;
91              }              }
92          }          }
93      } else if (d0<=0 && d2<=0) {      } else if (d0<=0 && d2<=0) {
# Line 95  Brick::Brick(int n0, int n1, int n2, dou Line 97  Brick::Brick(int n0, int n1, int n2, dou
97          if (d0*d1*d2 != m_mpiInfo->size) {          if (d0*d1*d2 != m_mpiInfo->size) {
98              // ratios not the same so subdivide side with more elements only              // ratios not the same so subdivide side with more elements only
99              if (n0>n2) {              if (n0>n2) {
100                    d0=0;
101                  d2=1;                  d2=1;
102              } else {              } else {
103                  d0=1;                  d0=1;
104                    d2=0;
105              }              }
106          }          }
107      } else if (d1<=0 && d2<=0) {      } else if (d1<=0 && d2<=0) {
# Line 107  Brick::Brick(int n0, int n1, int n2, dou Line 111  Brick::Brick(int n0, int n1, int n2, dou
111          if (d0*d1*d2 != m_mpiInfo->size) {          if (d0*d1*d2 != m_mpiInfo->size) {
112              // ratios not the same so subdivide side with more elements only              // ratios not the same so subdivide side with more elements only
113              if (n1>n2) {              if (n1>n2) {
114                    d1=0;
115                  d2=1;                  d2=1;
116              } else {              } else {
117                  d1=1;                  d1=1;
118                    d2=0;
119              }              }
120          }          }
121      }      }
# Line 232  bool Brick::operator==(const AbstractDom Line 238  bool Brick::operator==(const AbstractDom
238      return false;      return false;
239  }  }
240
241    void Brick::readBinaryGrid(escript::Data& out, string filename,
242                               const vector<int>& first,
243                               const vector<int>& numValues) const
244    {
245        // check destination function space
246        int myN0, myN1, myN2;
247        if (out.getFunctionSpace().getTypeCode() == Nodes) {
248            myN0 = m_N0;
249            myN1 = m_N1;
250            myN2 = m_N2;
251        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
252                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
253            myN0 = m_NE0;
254            myN1 = m_NE1;
255            myN2 = m_NE2;
256        } else
257            throw RipleyException("readBinaryGrid(): invalid function space for output data object");
258
259        // check file existence and size
260        ifstream f(filename.c_str(), ifstream::binary);
261        if (f.fail()) {
262            throw RipleyException("readBinaryGrid(): cannot open file");
263        }
264        f.seekg(0, ios::end);
265        const int numComp = out.getDataPointSize();
266        const int filesize = f.tellg();
267        const int reqsize = numValues[0]*numValues[1]*numValues[2]*numComp*sizeof(float);
268        if (filesize < reqsize) {
269            f.close();
270            throw RipleyException("readBinaryGrid(): not enough data in file");
271        }
272
273        // check if this rank contributes anything
274        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
275                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1 ||
276                first[2] >= m_offset2+myN2 || first[2]+numValues[2] <= m_offset2) {
277            f.close();
278            return;
279        }
280
281        // now determine how much this rank has to write
282
283        // first coordinates in data object to write to
284        const int first0 = max(0, first[0]-m_offset0);
285        const int first1 = max(0, first[1]-m_offset1);
286        const int first2 = max(0, first[2]-m_offset2);
287        // indices to first value in file
288        const int idx0 = max(0, m_offset0-first[0]);
289        const int idx1 = max(0, m_offset1-first[1]);
290        const int idx2 = max(0, m_offset2-first[2]);
291        // number of values to write
292        const int num0 = min(numValues[0]-idx0, myN0-first0);
293        const int num1 = min(numValues[1]-idx1, myN1-first1);
294        const int num2 = min(numValues[2]-idx2, myN2-first2);
295
296        out.requireWrite();
297        vector<float> values(num0*numComp);
298        const int dpp = out.getNumDataPointsPerSample();
299
300        for (index_t z=0; z<num2; z++) {
301            for (index_t y=0; y<num1; y++) {
302                const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]+(idx2+z)*numValues[0]*numValues[1]);
303                f.seekg(fileofs*sizeof(float));
305                for (index_t x=0; x<num0; x++) {
306                    double* dest = out.getSampleDataRW(first0+x+(first1+y)*myN0+(first2+z)*myN0*myN1);
307                    for (index_t c=0; c<numComp; c++) {
308                        for (index_t q=0; q<dpp; q++) {
309                            *dest++ = static_cast<double>(values[x*numComp+c]);
310                        }
311                    }
312                }
313            }
314        }
315
316        f.close();
317    }
318
319  void Brick::dump(const string& fileName) const  void Brick::dump(const string& fileName) const
320  {  {
321  #if USE_SILO  #if USE_SILO

Legend:
 Removed from v.3943 changed lines Added in v.3971