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

Diff of /trunk/ripley/src/Rectangle.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 159  bool Rectangle::operator==(const Abstrac Line 163  bool Rectangle::operator==(const Abstrac
163      return false;      return false;
164  }  }
165    
166    void Rectangle::readNcGrid(escript::Data& out, string filename, string varname,
167                const vector<int>& first, const vector<int>& numValues) const
168    {
169    #ifdef USE_NETCDF
170        // check destination function space
171        int myN0, myN1;
172        if (out.getFunctionSpace().getTypeCode() == Nodes) {
173            myN0 = m_N0;
174            myN1 = m_N1;
175        } else if (out.getFunctionSpace().getTypeCode() == Elements ||
176                    out.getFunctionSpace().getTypeCode() == ReducedElements) {
177            myN0 = m_NE0;
178            myN1 = m_NE1;
179        } else
180            throw RipleyException("readNcGrid(): invalid function space for output data object");
181    
182        if (first.size() != 2)
183            throw RipleyException("readNcGrid(): argument 'first' must have 2 entries");
184    
185        if (numValues.size() != 2)
186            throw RipleyException("readNcGrid(): argument 'numValues' must have 2 entries");
187    
188        // check file existence and size
189        NcFile f(filename.c_str(), NcFile::ReadOnly);
190        if (!f.is_valid())
191            throw RipleyException("readNcGrid(): cannot open file");
192    
193        NcVar* var = f.get_var(varname.c_str());
194        if (!var)
195            throw RipleyException("readNcGrid(): invalid variable");
196    
197        // TODO: rank>0 data support
198        const int numComp = out.getDataPointSize();
199        if (numComp > 1)
200            throw RipleyException("readNcGrid(): only scalar data supported");
201    
202        const int dims = var->num_dims();
203        const long *edges = var->edges();
204    
205        // is this a slice of the data object (dims!=2)?
206        // note the expected ordering of edges (as in numpy: y,x)
207        if ( (dims==2 && (numValues[1] > edges[0] || numValues[0] > edges[1]))
208                || (dims==1 && numValues[1]>1) ) {
209            throw RipleyException("readNcGrid(): not enough data in file");
210        }
211    
212        // check if this rank contributes anything
213        if (first[0] >= m_offset0+myN0 || first[0]+numValues[0] <= m_offset0 ||
214                first[1] >= m_offset1+myN1 || first[1]+numValues[1] <= m_offset1)
215            return;
216    
217        // now determine how much this rank has to write
218    
219        // first coordinates in data object to write to
220        const int first0 = max(0, first[0]-m_offset0);
221        const int first1 = max(0, first[1]-m_offset1);
222        // indices to first value in file
223        const int idx0 = max(0, m_offset0-first[0]);
224        const int idx1 = max(0, m_offset1-first[1]);
225        // number of values to write
226        const int num0 = min(numValues[0]-idx0, myN0-first0);
227        const int num1 = min(numValues[1]-idx1, myN1-first1);
228    
229        vector<double> values(num0*num1);
230        if (dims==2) {
231            var->set_cur(idx1, idx0);
232            var->get(&values[0], num1, num0);
233        } else {
234            var->set_cur(idx0);
235            var->get(&values[0], num0);
236        }
237    
238        const int dpp = out.getNumDataPointsPerSample();
239        out.requireWrite();
240    
241        for (index_t y=0; y<num1; y++) {
242    #pragma omp parallel for
243            for (index_t x=0; x<num0; x++) {
244                const int dataIndex = (first1+y)*myN0+first0+x;
245                const int srcIndex=y*num0+x;
246                double* dest = out.getSampleDataRW(dataIndex);
247                for (index_t q=0; q<dpp; q++) {
248                    *dest++ = values[srcIndex];
249                }
250            }
251        }
252    #else
253        throw RipleyException("readNcGrid(): not compiled with netCDF support");
254    #endif
255    }
256    
257  void Rectangle::readBinaryGrid(escript::Data& out, string filename,  void Rectangle::readBinaryGrid(escript::Data& out, string filename,
258              const vector<int>& first, const vector<int>& numValues) const              const vector<int>& first, const vector<int>& numValues) const
259  {  {

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

  ViewVC Help
Powered by ViewVC 1.1.26