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

Diff of /trunk/ripley/src/Brick.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4010 by caltinay, Tue Oct 2 06:57:11 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 321  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

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

  ViewVC Help
Powered by ViewVC 1.1.26