/[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 4174 by caltinay, Wed Jan 30 03:21:27 2013 UTC revision 4277 by caltinay, Wed Mar 6 01:30:41 2013 UTC
# Line 247  bool Brick::operator==(const AbstractDom Line 247  bool Brick::operator==(const AbstractDom
247    
248  void Brick::readBinaryGrid(escript::Data& out, string filename,  void Brick::readBinaryGrid(escript::Data& out, string filename,
249                             const vector<int>& first,                             const vector<int>& first,
250                             const vector<int>& numValues) const                             const vector<int>& numValues,
251                               const vector<int>& multiplier) const
252  {  {
253      // check destination function space      // check destination function space
254      int myN0, myN1, myN2;      int myN0, myN1, myN2;
# Line 263  void Brick::readBinaryGrid(escript::Data Line 264  void Brick::readBinaryGrid(escript::Data
264      } else      } else
265          throw RipleyException("readBinaryGrid(): invalid function space for output data object");          throw RipleyException("readBinaryGrid(): invalid function space for output data object");
266    
267        if (first.size() != 3)
268            throw RipleyException("readBinaryGrid(): argument 'first' must have 3 entries");
269    
270        if (numValues.size() != 3)
271            throw RipleyException("readBinaryGrid(): argument 'numValues' must have 3 entries");
272    
273        if (multiplier.size() != 3)
274            throw RipleyException("readBinaryGrid(): argument 'multiplier' must have 3 entries");
275        for (size_t i=0; i<multiplier.size(); i++)
276            if (multiplier[i]<1)
277                throw RipleyException("readBinaryGrid(): all multipliers must be positive");
278    
279      // check file existence and size      // check file existence and size
280      ifstream f(filename.c_str(), ifstream::binary);      ifstream f(filename.c_str(), ifstream::binary);
281      if (f.fail()) {      if (f.fail()) {
# Line 278  void Brick::readBinaryGrid(escript::Data Line 291  void Brick::readBinaryGrid(escript::Data
291      }      }
292    
293      // check if this rank contributes anything      // check if this rank contributes anything
294      if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||      if (first[0] >= m_offset0+myN0 || first[0]+numValues[0]*multiplier[0] <= m_offset0 ||
295              first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1 ||              first[1] >= m_offset1+myN1 || first[1]+numValues[1]*multiplier[1] <= m_offset1 ||
296              first[2] >= m_offset2+myN2 || first[2]+numValues[2] <= m_offset2) {              first[2] >= m_offset2+myN2 || first[2]+numValues[2]*multiplier[2] <= m_offset2) {
297          f.close();          f.close();
298          return;          return;
299      }      }
# Line 295  void Brick::readBinaryGrid(escript::Data Line 308  void Brick::readBinaryGrid(escript::Data
308      const int idx0 = max(0, m_offset0-first[0]);      const int idx0 = max(0, m_offset0-first[0]);
309      const int idx1 = max(0, m_offset1-first[1]);      const int idx1 = max(0, m_offset1-first[1]);
310      const int idx2 = max(0, m_offset2-first[2]);      const int idx2 = max(0, m_offset2-first[2]);
311      // number of values to write      // number of values to read
312      const int num0 = min(numValues[0]-idx0, myN0-first0);      const int num0 = min(numValues[0]-idx0, myN0-first0);
313      const int num1 = min(numValues[1]-idx1, myN1-first1);      const int num1 = min(numValues[1]-idx1, myN1-first1);
314      const int num2 = min(numValues[2]-idx2, myN2-first2);      const int num2 = min(numValues[2]-idx2, myN2-first2);
# Line 309  void Brick::readBinaryGrid(escript::Data Line 322  void Brick::readBinaryGrid(escript::Data
322              const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]+(idx2+z)*numValues[0]*numValues[1]);              const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]+(idx2+z)*numValues[0]*numValues[1]);
323              f.seekg(fileofs*sizeof(float));              f.seekg(fileofs*sizeof(float));
324              f.read((char*)&values[0], num0*numComp*sizeof(float));              f.read((char*)&values[0], num0*numComp*sizeof(float));
             const int dataIndex = first0+(first1+y)*myN0+(first2+z)*myN0*myN1;  
325    
326              for (index_t x=0; x<num0; x++) {              for (index_t x=0; x<num0; x++) {
327                  double* dest = out.getSampleDataRW(dataIndex+x);                  const int baseIndex = first0+x*multiplier[0]
328                  for (index_t c=0; c<numComp; c++) {                                          +(first1+y*multiplier[1])*myN0
329                      if (!isnan(values[x*numComp+c])) {                                          +(first2+z*multiplier[2])*myN0*myN1;
330                          for (index_t q=0; q<dpp; q++) {                  for (index_t m2=0; m2<multiplier[2]; m2++) {
331                              *dest++ = static_cast<double>(values[x*numComp+c]);                      for (index_t m1=0; m1<multiplier[1]; m1++) {
332                            for (index_t m0=0; m0<multiplier[0]; m0++) {
333                                const int dataIndex = baseIndex+m0
334                                               +m1*myN0
335                                               +m2*myN0*myN1;
336                                double* dest = out.getSampleDataRW(dataIndex);
337                                for (index_t c=0; c<numComp; c++) {
338                                    if (!isnan(values[x*numComp+c])) {
339                                        for (index_t q=0; q<dpp; q++) {
340                                            *dest++ = static_cast<double>(values[x*numComp+c]);
341                                        }
342                                    }
343                                }
344                          }                          }
345                      }                      }
346                  }                  }
# Line 328  void Brick::readBinaryGrid(escript::Data Line 352  void Brick::readBinaryGrid(escript::Data
352  }  }
353    
354  void Brick::readNcGrid(escript::Data& out, string filename, string varname,  void Brick::readNcGrid(escript::Data& out, string filename, string varname,
355              const vector<int>& first, const vector<int>& numValues) const              const vector<int>& first, const vector<int>& numValues,
356                const vector<int>& multiplier) const
357  {  {
358  #ifdef USE_NETCDF  #ifdef USE_NETCDF
359      // check destination function space      // check destination function space
# Line 351  void Brick::readNcGrid(escript::Data& ou Line 376  void Brick::readNcGrid(escript::Data& ou
376      if (numValues.size() != 3)      if (numValues.size() != 3)
377          throw RipleyException("readNcGrid(): argument 'numValues' must have 3 entries");          throw RipleyException("readNcGrid(): argument 'numValues' must have 3 entries");
378    
379        if (multiplier.size() != 3)
380            throw RipleyException("readNcGrid(): argument 'multiplier' must have 3 entries");
381        for (size_t i=0; i<multiplier.size(); i++)
382            if (multiplier[i]<1)
383                throw RipleyException("readNcGrid(): all multipliers must be positive");
384    
385      // check file existence and size      // check file existence and size
386      NcFile f(filename.c_str(), NcFile::ReadOnly);      NcFile f(filename.c_str(), NcFile::ReadOnly);
387      if (!f.is_valid())      if (!f.is_valid())
# Line 378  void Brick::readNcGrid(escript::Data& ou Line 409  void Brick::readNcGrid(escript::Data& ou
409      }      }
410    
411      // check if this rank contributes anything      // check if this rank contributes anything
412      if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||      if (first[0] >= m_offset0+myN0 || first[0]+numValues[0]*multiplier[0] <= m_offset0 ||
413              first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1 ||              first[1] >= m_offset1+myN1 || first[1]+numValues[1]*multiplier[1] <= m_offset1 ||
414              first[2] >= m_offset2+myN2 || first[2]+numValues[2] <= m_offset2) {              first[2] >= m_offset2+myN2 || first[2]+numValues[2]*multiplier[2] <= m_offset2) {
415          return;          return;
416      }      }
417    
# Line 394  void Brick::readNcGrid(escript::Data& ou Line 425  void Brick::readNcGrid(escript::Data& ou
425      const int idx0 = max(0, m_offset0-first[0]);      const int idx0 = max(0, m_offset0-first[0]);
426      const int idx1 = max(0, m_offset1-first[1]);      const int idx1 = max(0, m_offset1-first[1]);
427      const int idx2 = max(0, m_offset2-first[2]);      const int idx2 = max(0, m_offset2-first[2]);
428      // number of values to write      // number of values to read
429      const int num0 = min(numValues[0]-idx0, myN0-first0);      const int num0 = min(numValues[0]-idx0, myN0-first0);
430      const int num1 = min(numValues[1]-idx1, myN1-first1);      const int num1 = min(numValues[1]-idx1, myN1-first1);
431      const int num2 = min(numValues[2]-idx2, myN2-first2);      const int num2 = min(numValues[2]-idx2, myN2-first2);
# Line 418  void Brick::readNcGrid(escript::Data& ou Line 449  void Brick::readNcGrid(escript::Data& ou
449          for (index_t y=0; y<num1; y++) {          for (index_t y=0; y<num1; y++) {
450  #pragma omp parallel for  #pragma omp parallel for
451              for (index_t x=0; x<num0; x++) {              for (index_t x=0; x<num0; x++) {
452                  const int dataIndex = first0+(first1+y)*myN0+(first2+z)*myN0*myN1+x;                  const int baseIndex = first0+x*multiplier[0]
453                                            +(first1+y*multiplier[1])*myN0
454                                            +(first2+z*multiplier[2])*myN0*myN1;
455                  const int srcIndex=z*num1*num0+y*num0+x;                  const int srcIndex=z*num1*num0+y*num0+x;
456                  if (!isnan(values[srcIndex])) {                  if (!isnan(values[srcIndex])) {
457                      double* dest = out.getSampleDataRW(dataIndex);                      for (index_t m2=0; m2<multiplier[2]; m2++) {
458                      for (index_t q=0; q<dpp; q++) {                          for (index_t m1=0; m1<multiplier[1]; m1++) {
459                          *dest++ = values[srcIndex];                              for (index_t m0=0; m0<multiplier[0]; m0++) {
460                                    const int dataIndex = baseIndex+m0
461                                                   +m1*myN0
462                                                   +m2*myN0*myN1;
463                                    double* dest = out.getSampleDataRW(dataIndex);
464                                    for (index_t q=0; q<dpp; q++) {
465                                        *dest++ = values[srcIndex];
466                                    }
467                                }
468                            }
469                      }                      }
470                  }                  }
471              }              }

Legend:
Removed from v.4174  
changed lines
  Added in v.4277

  ViewVC Help
Powered by ViewVC 1.1.26