/[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 4529 by caltinay, Fri Oct 25 01:23:27 2013 UTC revision 4615 by caltinay, Mon Jan 13 05:05:33 2014 UTC
# Line 251  bool Brick::operator==(const AbstractDom Line 251  bool Brick::operator==(const AbstractDom
251  }  }
252    
253  void Brick::readNcGrid(escript::Data& out, string filename, string varname,  void Brick::readNcGrid(escript::Data& out, string filename, string varname,
254              const vector<int>& first, const vector<int>& numValues,              const GridParameters& params) const
             const vector<int>& multiplier) const  
255  {  {
256  #ifdef USE_NETCDF  #ifdef USE_NETCDF
257      // check destination function space      // check destination function space
# Line 269  void Brick::readNcGrid(escript::Data& ou Line 268  void Brick::readNcGrid(escript::Data& ou
268      } else      } else
269          throw RipleyException("readNcGrid(): invalid function space for output data object");          throw RipleyException("readNcGrid(): invalid function space for output data object");
270    
271      if (first.size() != 3)      if (params.first.size() != 3)
272          throw RipleyException("readNcGrid(): argument 'first' must have 3 entries");          throw RipleyException("readNcGrid(): argument 'first' must have 3 entries");
273    
274      if (numValues.size() != 3)      if (params.numValues.size() != 3)
275          throw RipleyException("readNcGrid(): argument 'numValues' must have 3 entries");          throw RipleyException("readNcGrid(): argument 'numValues' must have 3 entries");
276    
277      if (multiplier.size() != 3)      if (params.multiplier.size() != 3)
278          throw RipleyException("readNcGrid(): argument 'multiplier' must have 3 entries");          throw RipleyException("readNcGrid(): argument 'multiplier' must have 3 entries");
279      for (size_t i=0; i<multiplier.size(); i++)      for (size_t i=0; i<params.multiplier.size(); i++)
280          if (multiplier[i]<1)          if (params.multiplier[i]<1)
281              throw RipleyException("readNcGrid(): all multipliers must be positive");              throw RipleyException("readNcGrid(): all multipliers must be positive");
282    
283      // check file existence and size      // check file existence and size
# Line 300  void Brick::readNcGrid(escript::Data& ou Line 299  void Brick::readNcGrid(escript::Data& ou
299    
300      // is this a slice of the data object (dims!=3)?      // is this a slice of the data object (dims!=3)?
301      // note the expected ordering of edges (as in numpy: z,y,x)      // note the expected ordering of edges (as in numpy: z,y,x)
302      if ( (dims==3 && (numValues[2] > edges[0] || numValues[1] > edges[1]      if ( (dims==3 && (params.numValues[2] > edges[0] ||
303                        || numValues[0] > edges[2]))                        params.numValues[1] > edges[1] ||
304              || (dims==2 && numValues[2]>1)                        params.numValues[0] > edges[2]))
305              || (dims==1 && (numValues[2]>1 || numValues[1]>1)) ) {              || (dims==2 && params.numValues[2]>1)
306                || (dims==1 && (params.numValues[2]>1 || params.numValues[1]>1)) ) {
307          throw RipleyException("readNcGrid(): not enough data in file");          throw RipleyException("readNcGrid(): not enough data in file");
308      }      }
309    
310      // check if this rank contributes anything      // check if this rank contributes anything
311      if (first[0] >= m_offset[0]+myN0 || first[0]+numValues[0]*multiplier[0] <= m_offset[0] ||      if (params.first[0] >= m_offset[0]+myN0 ||
312              first[1] >= m_offset[1]+myN1 || first[1]+numValues[1]*multiplier[1] <= m_offset[1] ||              params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
313              first[2] >= m_offset[2]+myN2 || first[2]+numValues[2]*multiplier[2] <= m_offset[2]) {              params.first[1] >= m_offset[1]+myN1 ||
314                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1] ||
315                params.first[2] >= m_offset[2]+myN2 ||
316                params.first[2]+params.numValues[2]*params.multiplier[2] <= m_offset[2]) {
317          return;          return;
318      }      }
319    
320      // now determine how much this rank has to write      // now determine how much this rank has to write
321    
322      // first coordinates in data object to write to      // first coordinates in data object to write to
323      const int first0 = max(0, first[0]-m_offset[0]);      const int first0 = max(0, params.first[0]-m_offset[0]);
324      const int first1 = max(0, first[1]-m_offset[1]);      const int first1 = max(0, params.first[1]-m_offset[1]);
325      const int first2 = max(0, first[2]-m_offset[2]);      const int first2 = max(0, params.first[2]-m_offset[2]);
326      // indices to first value in file      // indices to first value in file
327      const int idx0 = max(0, m_offset[0]-first[0]);      const int idx0 = max(0, m_offset[0]-params.first[0]);
328      const int idx1 = max(0, m_offset[1]-first[1]);      const int idx1 = max(0, m_offset[1]-params.first[1]);
329      const int idx2 = max(0, m_offset[2]-first[2]);      const int idx2 = max(0, m_offset[2]-params.first[2]);
330      // number of values to read      // number of values to read
331      const int num0 = min(numValues[0]-idx0, myN0-first0);      const int num0 = min(params.numValues[0]-idx0, myN0-first0);
332      const int num1 = min(numValues[1]-idx1, myN1-first1);      const int num1 = min(params.numValues[1]-idx1, myN1-first1);
333      const int num2 = min(numValues[2]-idx2, myN2-first2);      const int num2 = min(params.numValues[2]-idx2, myN2-first2);
334    
335      vector<double> values(num0*num1*num2);      vector<double> values(num0*num1*num2);
336      if (dims==3) {      if (dims==3) {
# Line 348  void Brick::readNcGrid(escript::Data& ou Line 351  void Brick::readNcGrid(escript::Data& ou
351          for (index_t y=0; y<num1; y++) {          for (index_t y=0; y<num1; y++) {
352  #pragma omp parallel for  #pragma omp parallel for
353              for (index_t x=0; x<num0; x++) {              for (index_t x=0; x<num0; x++) {
354                  const int baseIndex = first0+x*multiplier[0]                  const int baseIndex = first0+x*params.multiplier[0]
355                                          +(first1+y*multiplier[1])*myN0                                       +(first1+y*params.multiplier[1])*myN0
356                                          +(first2+z*multiplier[2])*myN0*myN1;                                       +(first2+z*params.multiplier[2])*myN0*myN1;
357                  const int srcIndex=z*num1*num0+y*num0+x;                  const int srcIndex=z*num1*num0+y*num0+x;
358                  if (!isnan(values[srcIndex])) {                  if (!isnan(values[srcIndex])) {
359                      for (index_t m2=0; m2<multiplier[2]; m2++) {                      for (index_t m2=0; m2<params.multiplier[2]; m2++) {
360                          for (index_t m1=0; m1<multiplier[1]; m1++) {                          for (index_t m1=0; m1<params.multiplier[1]; m1++) {
361                              for (index_t m0=0; m0<multiplier[0]; m0++) {                              for (index_t m0=0; m0<params.multiplier[0]; m0++) {
362                                  const int dataIndex = baseIndex+m0                                  const int dataIndex = baseIndex+m0
363                                                 +m1*myN0                                                 +m1*myN0
364                                                 +m2*myN0*myN1;                                                 +m2*myN0*myN1;
# Line 376  void Brick::readNcGrid(escript::Data& ou Line 379  void Brick::readNcGrid(escript::Data& ou
379  }  }
380    
381  void Brick::readBinaryGrid(escript::Data& out, string filename,  void Brick::readBinaryGrid(escript::Data& out, string filename,
382                             const vector<int>& first,                             const GridParameters& params) const
                            const vector<int>& numValues,  
                            const std::vector<int>& multiplier,  
                            int byteOrder, int dataType) const  
383  {  {
384      // the mapping is not universally correct but should work on our      // the mapping is not universally correct but should work on our
385      // supported platforms      // supported platforms
386      switch (dataType) {      switch (params.dataType) {
387          case DATATYPE_INT32:          case DATATYPE_INT32:
388              readBinaryGridImpl<int>(out, filename, first, numValues,              readBinaryGridImpl<int>(out, filename, params);
                                     multiplier, byteOrder);  
389              break;              break;
390          case DATATYPE_FLOAT32:          case DATATYPE_FLOAT32:
391              readBinaryGridImpl<float>(out, filename, first, numValues,              readBinaryGridImpl<float>(out, filename, params);
                                       multiplier, byteOrder);  
392              break;              break;
393          case DATATYPE_FLOAT64:          case DATATYPE_FLOAT64:
394              readBinaryGridImpl<double>(out, filename, first, numValues,              readBinaryGridImpl<double>(out, filename, params);
                                        multiplier, byteOrder);  
395              break;              break;
396          default:          default:
397              throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");              throw RipleyException("readBinaryGrid(): invalid or unsupported datatype");
# Line 403  void Brick::readBinaryGrid(escript::Data Line 400  void Brick::readBinaryGrid(escript::Data
400    
401  template<typename ValueType>  template<typename ValueType>
402  void Brick::readBinaryGridImpl(escript::Data& out, const string& filename,  void Brick::readBinaryGridImpl(escript::Data& out, const string& filename,
403                                 const vector<int>& first,                                 const GridParameters& params) const
                                const vector<int>& numValues,  
                                const vector<int>& multiplier,  
                                int byteOrder) const  
404  {  {
405      // check destination function space      // check destination function space
406      int myN0, myN1, myN2;      int myN0, myN1, myN2;
# Line 422  void Brick::readBinaryGridImpl(escript:: Line 416  void Brick::readBinaryGridImpl(escript::
416      } else      } else
417          throw RipleyException("readBinaryGrid(): invalid function space for output data object");          throw RipleyException("readBinaryGrid(): invalid function space for output data object");
418    
419      if (first.size() != 3)      if (params.first.size() != 3)
420          throw RipleyException("readBinaryGrid(): argument 'first' must have 3 entries");          throw RipleyException("readBinaryGrid(): argument 'first' must have 3 entries");
421    
422      if (numValues.size() != 3)      if (params.numValues.size() != 3)
423          throw RipleyException("readBinaryGrid(): argument 'numValues' must have 3 entries");          throw RipleyException("readBinaryGrid(): argument 'numValues' must have 3 entries");
424    
425      if (multiplier.size() != 3)      if (params.multiplier.size() != 3)
426          throw RipleyException("readBinaryGrid(): argument 'multiplier' must have 3 entries");          throw RipleyException("readBinaryGrid(): argument 'multiplier' must have 3 entries");
427      for (size_t i=0; i<multiplier.size(); i++)      for (size_t i=0; i<params.multiplier.size(); i++)
428          if (multiplier[i]<1)          if (params.multiplier[i]<1)
429              throw RipleyException("readBinaryGrid(): all multipliers must be positive");              throw RipleyException("readBinaryGrid(): all multipliers must be positive");
430    
431      // check file existence and size      // check file existence and size
# Line 442  void Brick::readBinaryGridImpl(escript:: Line 436  void Brick::readBinaryGridImpl(escript::
436      f.seekg(0, ios::end);      f.seekg(0, ios::end);
437      const int numComp = out.getDataPointSize();      const int numComp = out.getDataPointSize();
438      const int filesize = f.tellg();      const int filesize = f.tellg();
439      const int reqsize = numValues[0]*numValues[1]*numValues[2]*numComp*sizeof(ValueType);      const int reqsize = params.numValues[0]*params.numValues[1]*params.numValues[2]*numComp*sizeof(ValueType);
440      if (filesize < reqsize) {      if (filesize < reqsize) {
441          f.close();          f.close();
442          throw RipleyException("readBinaryGrid(): not enough data in file");          throw RipleyException("readBinaryGrid(): not enough data in file");
443      }      }
444    
445      // check if this rank contributes anything      // check if this rank contributes anything
446      if (first[0] >= m_offset[0]+myN0 || first[0]+numValues[0]*multiplier[0] <= m_offset[0] ||      if (params.first[0] >= m_offset[0]+myN0 ||
447              first[1] >= m_offset[1]+myN1 || first[1]+numValues[1]*multiplier[1] <= m_offset[1] ||              params.first[0]+params.numValues[0]*params.multiplier[0] <= m_offset[0] ||
448              first[2] >= m_offset[2]+myN2 || first[2]+numValues[2]*multiplier[2] <= m_offset[2]) {              params.first[1] >= m_offset[1]+myN1 ||
449                params.first[1]+params.numValues[1]*params.multiplier[1] <= m_offset[1] ||
450                params.first[2] >= m_offset[2]+myN2 ||
451                params.first[2]+params.numValues[2]*params.multiplier[2] <= m_offset[2]) {
452          f.close();          f.close();
453          return;          return;
454      }      }
# Line 459  void Brick::readBinaryGridImpl(escript:: Line 456  void Brick::readBinaryGridImpl(escript::
456      // now determine how much this rank has to write      // now determine how much this rank has to write
457    
458      // first coordinates in data object to write to      // first coordinates in data object to write to
459      const int first0 = max(0, first[0]-m_offset[0]);      const int first0 = max(0, params.first[0]-m_offset[0]);
460      const int first1 = max(0, first[1]-m_offset[1]);      const int first1 = max(0, params.first[1]-m_offset[1]);
461      const int first2 = max(0, first[2]-m_offset[2]);      const int first2 = max(0, params.first[2]-m_offset[2]);
462      // indices to first value in file      // indices to first value in file
463      const int idx0 = max(0, m_offset[0]-first[0]);      const int idx0 = max(0, m_offset[0]-params.first[0]);
464      const int idx1 = max(0, m_offset[1]-first[1]);      const int idx1 = max(0, m_offset[1]-params.first[1]);
465      const int idx2 = max(0, m_offset[2]-first[2]);      const int idx2 = max(0, m_offset[2]-params.first[2]);
466      // number of values to read      // number of values to read
467      const int num0 = min(numValues[0]-idx0, myN0-first0);      const int num0 = min(params.numValues[0]-idx0, myN0-first0);
468      const int num1 = min(numValues[1]-idx1, myN1-first1);      const int num1 = min(params.numValues[1]-idx1, myN1-first1);
469      const int num2 = min(numValues[2]-idx2, myN2-first2);      const int num2 = min(params.numValues[2]-idx2, myN2-first2);
470    
471      out.requireWrite();      out.requireWrite();
472      vector<ValueType> values(num0*numComp);      vector<ValueType> values(num0*numComp);
# Line 477  void Brick::readBinaryGridImpl(escript:: Line 474  void Brick::readBinaryGridImpl(escript::
474    
475      for (int z=0; z<num2; z++) {      for (int z=0; z<num2; z++) {
476          for (int y=0; y<num1; y++) {          for (int y=0; y<num1; y++) {
477              const int fileofs = numComp*(idx0+(idx1+y)*numValues[0]+(idx2+z)*numValues[0]*numValues[1]);              const int fileofs = numComp*(idx0+(idx1+y)*params.numValues[0]
478                                 +(idx2+z)*params.numValues[0]*params.numValues[1]);
479              f.seekg(fileofs*sizeof(ValueType));              f.seekg(fileofs*sizeof(ValueType));
480              f.read((char*)&values[0], num0*numComp*sizeof(ValueType));              f.read((char*)&values[0], num0*numComp*sizeof(ValueType));
481    
482              for (int x=0; x<num0; x++) {              for (int x=0; x<num0; x++) {
483                  const int baseIndex = first0+x*multiplier[0]                  const int baseIndex = first0+x*params.multiplier[0]
484                                          +(first1+y*multiplier[1])*myN0                                       +(first1+y*params.multiplier[1])*myN0
485                                          +(first2+z*multiplier[2])*myN0*myN1;                                       +(first2+z*params.multiplier[2])*myN0*myN1;
486                  for (int m2=0; m2<multiplier[2]; m2++) {                  for (int m2=0; m2<params.multiplier[2]; m2++) {
487                      for (int m1=0; m1<multiplier[1]; m1++) {                      for (int m1=0; m1<params.multiplier[1]; m1++) {
488                          for (int m0=0; m0<multiplier[0]; m0++) {                          for (int m0=0; m0<params.multiplier[0]; m0++) {
489                              const int dataIndex = baseIndex+m0                              const int dataIndex = baseIndex+m0
490                                             +m1*myN0                                             +m1*myN0
491                                             +m2*myN0*myN1;                                             +m2*myN0*myN1;
# Line 495  void Brick::readBinaryGridImpl(escript:: Line 493  void Brick::readBinaryGridImpl(escript::
493                              for (int c=0; c<numComp; c++) {                              for (int c=0; c<numComp; c++) {
494                                  ValueType val = values[x*numComp+c];                                  ValueType val = values[x*numComp+c];
495    
496                                  if (byteOrder != BYTEORDER_NATIVE) {                                  if (params.byteOrder != BYTEORDER_NATIVE) {
497                                      char* cval = reinterpret_cast<char*>(&val);                                      char* cval = reinterpret_cast<char*>(&val);
498                                      // this will alter val!!                                      // this will alter val!!
499                                      byte_swap32(cval);                                      byte_swap32(cval);

Legend:
Removed from v.4529  
changed lines
  Added in v.4615

  ViewVC Help
Powered by ViewVC 1.1.26