/[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 3981 by jfenwick, Fri Sep 21 02:47:54 2012 UTC revision 4013 by caltinay, Thu Oct 4 03:13:27 2012 UTC
# Line 18  extern "C" { Line 18  extern "C" {
18  #include <paso/SystemMatrix.h>  #include <paso/SystemMatrix.h>
19  }  }
20    
21    #ifdef USE_NETCDF
22    #include <netcdfcpp.h>
23    #endif
24    
25  #if USE_SILO  #if USE_SILO
26  #include <silo.h>  #include <silo.h>
27  #ifdef ESYS_MPI  #ifdef ESYS_MPI
# Line 54  Brick::Brick(int n0, int n1, int n2, dou Line 58  Brick::Brick(int n0, int n1, int n2, dou
58      if (d0<=0 && d1<=0 && d2<=0) {      if (d0<=0 && d1<=0 && d2<=0) {
59          warn=true;          warn=true;
60          d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1.f/3));          d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1.f/3));
61          d1=(int)(d0*n1/(float)n0);          d0=max(1, d0);
62            d1=max(1, (int)(d0*n1/(float)n0));
63          d2=m_mpiInfo->size/(d0*d1);          d2=m_mpiInfo->size/(d0*d1);
64          if (d0*d1*d2 != m_mpiInfo->size) {          if (d0*d1*d2 != m_mpiInfo->size) {
65              // ratios not the same so leave "smallest" side undivided and try              // ratios not the same so leave "smallest" side undivided and try
# Line 80  Brick::Brick(int n0, int n1, int n2, dou Line 85  Brick::Brick(int n0, int n1, int n2, dou
85      }      }
86      if (d0<=0 && d1<=0) {      if (d0<=0 && d1<=0) {
87          warn=true;          warn=true;
88          d0=int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1)));          d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1))));
89          d1=m_mpiInfo->size/d0;          d1=m_mpiInfo->size/d0;
90          if (d0*d1*d2 != m_mpiInfo->size) {          if (d0*d1*d2 != m_mpiInfo->size) {
91              // ratios not the same so subdivide side with more elements only              // ratios not the same so subdivide side with more elements only
# Line 94  Brick::Brick(int n0, int n1, int n2, dou Line 99  Brick::Brick(int n0, int n1, int n2, dou
99          }          }
100      } else if (d0<=0 && d2<=0) {      } else if (d0<=0 && d2<=0) {
101          warn=true;          warn=true;
102          d0=int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n2+1)));          d0=max(1, int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n2+1))));
103          d2=m_mpiInfo->size/d0;          d2=m_mpiInfo->size/d0;
104          if (d0*d1*d2 != m_mpiInfo->size) {          if (d0*d1*d2 != m_mpiInfo->size) {
105              // ratios not the same so subdivide side with more elements only              // ratios not the same so subdivide side with more elements only
# Line 108  Brick::Brick(int n0, int n1, int n2, dou Line 113  Brick::Brick(int n0, int n1, int n2, dou
113          }          }
114      } else if (d1<=0 && d2<=0) {      } else if (d1<=0 && d2<=0) {
115          warn=true;          warn=true;
116          d1=int(sqrt(m_mpiInfo->size*(n1+1)/(float)(n2+1)));          d1=max(1, int(sqrt(m_mpiInfo->size*(n1+1)/(float)(n2+1))));
117          d2=m_mpiInfo->size/d1;          d2=m_mpiInfo->size/d1;
118          if (d0*d1*d2 != m_mpiInfo->size) {          if (d0*d1*d2 != m_mpiInfo->size) {
119              // ratios not the same so subdivide side with more elements only              // ratios not the same so subdivide side with more elements only
# Line 304  void Brick::readBinaryGrid(escript::Data Line 309  void Brick::readBinaryGrid(escript::Data
309              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]);
310              f.seekg(fileofs*sizeof(float));              f.seekg(fileofs*sizeof(float));
311              f.read((char*)&values[0], num0*numComp*sizeof(float));              f.read((char*)&values[0], num0*numComp*sizeof(float));
312                const int dataIndex = first0+(first1+y)*myN0+(first2+z)*myN0*myN1;
313    
314              for (index_t x=0; x<num0; x++) {              for (index_t x=0; x<num0; x++) {
315                  double* dest = out.getSampleDataRW(first0+x+(first1+y)*myN0+(first2+z)*myN0*myN1);                  double* dest = out.getSampleDataRW(dataIndex+x);
316                  for (index_t c=0; c<numComp; c++) {                  for (index_t c=0; c<numComp; c++) {
317                      for (index_t q=0; q<dpp; q++) {                      for (index_t q=0; q<dpp; q++) {
318                          *dest++ = static_cast<double>(values[x*numComp+c]);                          *dest++ = static_cast<double>(values[x*numComp+c]);
# Line 318  void Brick::readBinaryGrid(escript::Data Line 325  void Brick::readBinaryGrid(escript::Data
325      f.close();      f.close();
326  }  }
327    
328    void Brick::readNcGrid(escript::Data& out, string filename, string varname,
329                const vector<int>& first, const vector<int>& numValues) const
330    {
331    #ifdef USE_NETCDF
332        // check destination function space
333        int myN0, myN1, myN2;
334        if (out.getFunctionSpace().getTypeCode() == Nodes) {
335            myN0 = m_N0;
336            myN1 = m_N1;
337            myN2 = m_N2;
338        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
339                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
340            myN0 = m_NE0;
341            myN1 = m_NE1;
342            myN2 = m_NE2;
343        } else
344            throw RipleyException("readNcGrid(): invalid function space for output data object");
345    
346        if (first.size() != 3)
347            throw RipleyException("readNcGrid(): argument 'first' must have 3 entries");
348    
349        if (numValues.size() != 3)
350            throw RipleyException("readNcGrid(): argument 'numValues' must have 3 entries");
351    
352        // check file existence and size
353        NcFile f(filename.c_str(), NcFile::ReadOnly);
354        if (!f.is_valid())
355            throw RipleyException("readNcGrid(): cannot open file");
356    
357        NcVar* var = f.get_var(varname.c_str());
358        if (!var)
359            throw RipleyException("readNcGrid(): invalid variable name");
360    
361        // TODO: rank>0 data support
362        const int numComp = out.getDataPointSize();
363        if (numComp > 1)
364            throw RipleyException("readNcGrid(): only scalar data supported");
365    
366        const int dims = var->num_dims();
367        const long *edges = var->edges();
368    
369        // is this a slice of the data object (dims!=3)?
370        // note the expected ordering of edges (as in numpy: z,y,x)
371        if ( (dims==3 && (numValues[2] > edges[0] || numValues[1] > edges[1]
372                          || numValues[0] > edges[2]))
373                || (dims==2 && numValues[2]>1)
374                || (dims==1 && (numValues[2]>1 || numValues[1]>1)) ) {
375            throw RipleyException("readNcGrid(): not enough data in file");
376        }
377    
378        // check if this rank contributes anything
379        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
380                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1 ||
381                first[2] >= m_offset2+myN2 || first[2]+numValues[2] <= m_offset2) {
382            return;
383        }
384    
385        // now determine how much this rank has to write
386    
387        // first coordinates in data object to write to
388        const int first0 = max(0, first[0]-m_offset0);
389        const int first1 = max(0, first[1]-m_offset1);
390        const int first2 = max(0, first[2]-m_offset2);
391        // indices to first value in file
392        const int idx0 = max(0, m_offset0-first[0]);
393        const int idx1 = max(0, m_offset1-first[1]);
394        const int idx2 = max(0, m_offset2-first[2]);
395        // number of values to write
396        const int num0 = min(numValues[0]-idx0, myN0-first0);
397        const int num1 = min(numValues[1]-idx1, myN1-first1);
398        const int num2 = min(numValues[2]-idx2, myN2-first2);
399    
400        vector<double> values(num0*num1*num2);
401        if (dims==3) {
402            var->set_cur(idx2, idx1, idx0);
403            var->get(&values[0], num2, num1, num0);
404        } else if (dims==2) {
405            var->set_cur(idx1, idx0);
406            var->get(&values[0], num1, num0);
407        } else {
408            var->set_cur(idx0);
409            var->get(&values[0], num0);
410        }
411    
412        const int dpp = out.getNumDataPointsPerSample();
413        out.requireWrite();
414    
415        for (index_t z=0; z<num2; z++) {
416            for (index_t y=0; y<num1; y++) {
417    #pragma omp parallel for
418                for (index_t x=0; x<num0; x++) {
419                    const int dataIndex = first0+(first1+y)*myN0+(first2+z)*myN0*myN1+x;
420                    const int srcIndex=z*num1*num0+y*num0+x;
421                    double* dest = out.getSampleDataRW(dataIndex);
422                    for (index_t q=0; q<dpp; q++) {
423                        *dest++ = values[srcIndex];
424                    }
425                }
426            }
427        }
428    #else
429        throw RipleyException("readNcGrid(): not compiled with netCDF support");
430    #endif
431    }
432    
433  void Brick::dump(const string& fileName) const  void Brick::dump(const string& fileName) const
434  {  {
435  #if USE_SILO  #if USE_SILO
# Line 1758  void Brick::dofToNodes(escript::Data& ou Line 1870  void Brick::dofToNodes(escript::Data& ou
1870                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);                  : &buffer[(m_dofMap[i]-numDOF)*numComp]);
1871          copy(src, src+numComp, out.getSampleDataRW(i));          copy(src, src+numComp, out.getSampleDataRW(i));
1872      }      }
1873        Paso_Coupler_free(coupler);
1874  }  }
1875    
1876  //private  //private

Legend:
Removed from v.3981  
changed lines
  Added in v.4013

  ViewVC Help
Powered by ViewVC 1.1.26