/[escript]/trunk/escript/src/DataFactory.cpp
ViewVC logotype

Annotation of /trunk/escript/src/DataFactory.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1872 - (hide annotations)
Mon Oct 13 00:18:55 2008 UTC (11 years ago) by jfenwick
File size: 12399 byte(s)
Closing the moreshared branch

1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 ksteube 1312
14 ksteube 1811
15 jgs 474 #include "DataFactory.h"
16 phornby 1628 #include "esysUtils/esys_malloc.h"
17 jgs 82
18     #include <boost/python/extract.hpp>
19 gross 950 #include <iostream>
20 gross 1487 #include <exception>
21 gross 1023 #ifdef USE_NETCDF
22 ksteube 1312 #include <netcdfcpp.h>
23 gross 1023 #endif
24 ksteube 1748 #ifdef PASO_MPI
25     #include <mpi.h>
26     #endif
27 jgs 82
28     using namespace boost::python;
29    
30     namespace escript {
31    
32     Data
33     Scalar(double value,
34     const FunctionSpace& what,
35     bool expanded)
36     {
37     //
38     // an empty shape is a scalar
39 jfenwick 1796 DataTypes::ShapeType shape;
40 jgs 82 return Data(value,shape,what,expanded);
41     }
42    
43     Data
44     Vector(double value,
45     const FunctionSpace& what,
46     bool expanded)
47     {
48 jfenwick 1872 DataTypes::ShapeType shape(1,what.getDomain()->getDim());
49 jgs 82 return Data(value,shape,what,expanded);
50     }
51    
52     Data
53     Tensor(double value,
54     const FunctionSpace& what,
55     bool expanded)
56     {
57 jfenwick 1872 DataTypes::ShapeType shape(2,what.getDomain()->getDim());
58 jgs 82 return Data(value,shape,what,expanded);
59     }
60    
61     Data
62     Tensor3(double value,
63     const FunctionSpace& what,
64     bool expanded)
65     {
66 jfenwick 1872 DataTypes::ShapeType shape(3,what.getDomain()->getDim());
67 jgs 82 return Data(value,shape,what,expanded);
68     }
69    
70     Data
71     Tensor4(double value,
72     const FunctionSpace& what,
73     bool expanded)
74     {
75 jfenwick 1872 DataTypes::ShapeType shape(4,what.getDomain()->getDim());
76 jgs 82 return Data(value,shape,what,expanded);
77     }
78    
79 gross 950 Data
80     load(const std::string fileName,
81     const AbstractDomain& domain)
82     {
83 gross 1023 #ifdef USE_NETCDF
84 gross 950 NcAtt *type_att, *rank_att, *function_space_type_att;
85     // netCDF error handler
86 gross 1141 NcError err(NcError::silent_nonfatal);
87 ksteube 1748 int mpi_iam=0, mpi_num=1;
88 gross 950 // Create the file.
89 ksteube 1748 #ifdef PASO_MPI
90     MPI_Comm_rank(MPI_COMM_WORLD, &mpi_iam);
91     MPI_Comm_size(MPI_COMM_WORLD, &mpi_num);
92     #endif
93     char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
94     NcFile dataFile(newFileName, NcFile::ReadOnly);
95 gross 950 if (!dataFile.is_valid())
96     throw DataException("Error - load:: opening of netCDF file for input failed.");
97     /* recover function space */
98     if (! (function_space_type_att=dataFile.get_att("function_space_type")) )
99     throw DataException("Error - load:: cannot recover function_space_type attribute from escript netCDF file.");
100     int function_space_type = function_space_type_att->as_int(0);
101     delete function_space_type_att;
102     /* test if function space id is valid and create function space instance */
103     if (! domain.isValidFunctionSpaceType(function_space_type) )
104     throw DataException("Error - load:: function space type code in netCDF file is invalid for given domain.");
105 jfenwick 1872 FunctionSpace function_space=FunctionSpace(domain.getPtr(), function_space_type);
106 gross 950 /* recover rank */
107     if (! (rank_att=dataFile.get_att("rank")) )
108     throw DataException("Error - load:: cannot recover rank attribute from escript netCDF file.");
109     int rank = rank_att->as_int(0);
110     delete rank_att;
111 jfenwick 1796 if (rank<0 || rank>DataTypes::maxRank)
112 gross 950 throw DataException("Error - load:: rank in escript netCDF file is greater than maximum rank.");
113     /* recover type attribute */
114     int type=-1;
115 gross 1141 if ((type_att=dataFile.get_att("type")) ) {
116     char* type_str = type_att->as_string(0);
117     if (strncmp(type_str, "constant", strlen("constant")) == 0 ) {
118     type =0;
119     } else if (strncmp(type_str, "tagged", strlen("tagged")) == 0 ) {
120     type =1;
121     } else if (strncmp(type_str, "expanded", strlen("expanded")) == 0 ) {
122     type =2;
123     }
124 phornby 1628 esysUtils::free(type_str);
125 gross 1141 } else {
126     if (! (type_att=dataFile.get_att("type_id")) )
127     throw DataException("Error - load:: cannot recover type attribute from escript netCDF file.");
128     type=type_att->as_int(0);
129 gross 950 }
130     delete type_att;
131 gross 1141
132 gross 950 /* recover dimension */
133     int ndims=dataFile.num_dims();
134 gross 983 int ntags =0 , num_samples =0 , num_data_points_per_sample =0, d=0, len_data_point=1;
135 gross 967 NcDim *d_dim, *tags_dim, *num_samples_dim, *num_data_points_per_sample_dim;
136 gross 950 /* recover shape */
137 jfenwick 1796 DataTypes::ShapeType shape;
138     long dims[DataTypes::maxRank+2];
139 gross 950 if (rank>0) {
140     if (! (d_dim=dataFile.get_dim("d0")) )
141     throw DataException("Error - load:: unable to recover d0 from netCDF file.");
142     d=d_dim->size();
143     shape.push_back(d);
144     dims[0]=d;
145 gross 983 len_data_point*=d;
146 gross 950 }
147     if (rank>1) {
148     if (! (d_dim=dataFile.get_dim("d1")) )
149     throw DataException("Error - load:: unable to recover d1 from netCDF file.");
150     d=d_dim->size();
151     shape.push_back(d);
152     dims[1]=d;
153 gross 983 len_data_point*=d;
154 gross 950 }
155     if (rank>2) {
156     if (! (d_dim=dataFile.get_dim("d2")) )
157     throw DataException("Error - load:: unable to recover d2 from netCDF file.");
158     d=d_dim->size();
159     shape.push_back(d);
160     dims[2]=d;
161 gross 983 len_data_point*=d;
162 gross 950 }
163     if (rank>3) {
164     if (! (d_dim=dataFile.get_dim("d3")) )
165     throw DataException("Error - load:: unable to recover d3 from netCDF file.");
166     d=d_dim->size();
167     shape.push_back(d);
168     dims[3]=d;
169 gross 983 len_data_point*=d;
170 gross 950 }
171     /* recover stuff */
172     Data out;
173 gross 983 NcVar *var, *ids_var, *tags_var;
174 gross 950 if (type == 0) {
175     /* constant data */
176     if ( ! ( (ndims == rank && rank >0) || ( ndims ==1 && rank == 0 ) ) )
177     throw DataException("Error - load:: illegal number of dimensions for constant data in netCDF file.");
178     if (rank == 0) {
179     if (! (d_dim=dataFile.get_dim("l")) )
180     throw DataException("Error - load:: unable to recover d0 for scalar constant data in netCDF file.");
181     int d0 = d_dim->size();
182     if (! d0 == 1)
183     throw DataException("Error - load:: d0 is expected to be one for scalar constant data in netCDF file.");
184     dims[0]=1;
185     }
186     out=Data(0,shape,function_space);
187     if (!(var = dataFile.get_var("data")))
188     throw DataException("Error - load:: unable to find data in netCDF file.");
189 jfenwick 1796 if (! var->get(&(out.getDataAtOffset(out.getDataOffset(0,0))), dims) )
190 gross 950 throw DataException("Error - load:: unable to recover data from netCDF file.");
191     } else if (type == 1) {
192     /* tagged data */
193     if ( ! (ndims == rank + 1) )
194 phornby 1020 throw DataException("Error - load:: illegal number of dimensions for tagged data in netCDF file.");
195 gross 983 if (! (tags_dim=dataFile.get_dim("num_tags")) )
196 phornby 1020 throw DataException("Error - load:: unable to recover number of tags from netCDF file.");
197 gross 950 ntags=tags_dim->size();
198 gross 983 dims[rank]=ntags;
199 phornby 1628 int *tags = (int *) esysUtils::malloc(ntags*sizeof(int));
200 gross 983 if (! ( tags_var = dataFile.get_var("tags")) )
201 phornby 1020 {
202 phornby 1628 esysUtils::free(tags);
203 gross 983 throw DataException("Error - load:: unable to find tags in netCDF file.");
204 phornby 1020 }
205 gross 983 if (! tags_var->get(tags, ntags) )
206 phornby 1020 {
207 phornby 1628 esysUtils::free(tags);
208 phornby 1020 throw DataException("Error - load:: unable to recover tags from netCDF file.");
209     }
210 gross 983
211 jfenwick 1796 // Current Version
212     /* DataVector data(len_data_point * ntags, 0., len_data_point * ntags);
213 gross 983 if (!(var = dataFile.get_var("data")))
214 phornby 1020 {
215 phornby 1628 esysUtils::free(tags);
216 phornby 1020 throw DataException("Error - load:: unable to find data in netCDF file.");
217     }
218 gross 983 if (! var->get(&(data[0]), dims) )
219 phornby 1020 {
220 phornby 1628 esysUtils::free(tags);
221 phornby 1020 throw DataException("Error - load:: unable to recover data from netCDF file.");
222     }
223 gross 983 out=Data(DataArrayView(data,shape,0),function_space);
224     for (int t=1; t<ntags; ++t) {
225 jfenwick 1796 out.setTaggedValueFromCPP(tags[t],shape, data, t*len_data_point);
226     // out.setTaggedValueFromCPP(tags[t],DataArrayView(data,shape,t*len_data_point));
227     }*/
228     // End current version
229    
230     // New version
231    
232     // A) create a DataTagged dt
233     // B) Read data from file
234     // C) copy default value into dt
235     // D) copy tagged values into dt
236     // E) create a new Data based on dt
237    
238     NcVar* var1;
239     DataVector data1(len_data_point * ntags, 0., len_data_point * ntags);
240     if (!(var1 = dataFile.get_var("data")))
241     {
242     esysUtils::free(tags);
243     throw DataException("Error - load:: unable to find data in netCDF file.");
244 gross 983 }
245 jfenwick 1796 if (! var1->get(&(data1[0]), dims) )
246     {
247     esysUtils::free(tags);
248     throw DataException("Error - load:: unable to recover data from netCDF file.");
249     }
250     DataTagged* dt=new DataTagged(function_space, shape, tags,data1);
251     out=Data(dt);
252 phornby 1628 esysUtils::free(tags);
253 gross 950 } else if (type == 2) {
254     /* expanded data */
255     if ( ! (ndims == rank + 2) )
256     throw DataException("Error - load:: illegal number of dimensions for exanded data in netCDF file.");
257 gross 967 if ( ! (num_samples_dim = dataFile.get_dim("num_samples") ) )
258 gross 950 throw DataException("Error - load:: unable to recover number of samples from netCDF file.");
259 gross 967 num_samples = num_samples_dim->size();
260     if ( ! (num_data_points_per_sample_dim = dataFile.get_dim("num_data_points_per_sample") ) )
261 gross 950 throw DataException("Error - load:: unable to recover number of data points per sample from netCDF file.");
262 gross 967 num_data_points_per_sample=num_data_points_per_sample_dim->size();
263 gross 982 // check shape:
264 gross 967 if ( ! (num_samples == function_space.getNumSamples() && num_data_points_per_sample == function_space.getNumDataPointsPerSample()) )
265     throw DataException("Error - load:: data sample layout of file does not match data layout of function space.");
266 jfenwick 1808 if (num_samples==0) {
267     out = Data(0,shape,function_space,true);
268 phornby 1020 }
269 jfenwick 1808 else {
270     // get ids
271     if (! ( ids_var = dataFile.get_var("id")) )
272     throw DataException("Error - load:: unable to find reference ids in netCDF file.");
273     const int* ids_p=function_space.borrowSampleReferenceIDs();
274     int *ids_of_nc = (int *)esysUtils::malloc(num_samples*sizeof(int));
275     if (! ids_var->get(ids_of_nc, (long) num_samples) )
276     {
277     esysUtils::free(ids_of_nc);
278     throw DataException("Error - load:: unable to recover ids from netCDF file.");
279     }
280     // check order:
281     int failed=-1, local_failed=-1, i;
282     #pragma omp parallel private(local_failed)
283     {
284     local_failed=-1;
285     #pragma omp for private(i) schedule(static)
286     for (i=0;i < num_samples; ++i) {
287     if (ids_of_nc[i]!=ids_p[i]) local_failed=i;
288     }
289     #pragma omp critical
290     if (local_failed>=0) failed = local_failed;
291     }
292     /* if (failed>=0)
293     {
294     esysUtils::free(ids_of_nc);
295     throw DataException("Error - load:: data ordering in netCDF file does not match ordering of FunctionSpace.");
296     } */
297     // get the data:
298     dims[rank]=num_data_points_per_sample;
299     dims[rank+1]=num_samples;
300     out=Data(0,shape,function_space,true);
301     if (!(var = dataFile.get_var("data")))
302     {
303     esysUtils::free(ids_of_nc);
304     throw DataException("Error - load:: unable to find data in netCDF file.");
305     }
306     if (! var->get(&(out.getDataAtOffset(out.getDataOffset(0,0))), dims) )
307     {
308     esysUtils::free(ids_of_nc);
309     throw DataException("Error - load:: unable to recover data from netCDF file.");
310     }
311     if (failed>=0) {
312     try {
313     std::cout << "Information - load: start reordering data from netCDF file " << fileName << std::endl;
314     out.borrowData()->reorderByReferenceIDs(ids_of_nc);
315     }
316     catch (std::exception&) {
317     esysUtils::free(ids_of_nc);
318     throw DataException("Error - load:: unable to reorder data in netCDF file.");
319     }
320     }
321 gross 982 }
322 gross 950 } else {
323     throw DataException("Error - load:: unknown escript data type in netCDF file.");
324     }
325     return out;
326 gross 1023 #else
327     throw DataException("Error - load:: is not compiled with netCFD. Please contact your insatllation manager.");
328     #endif
329     }
330 gross 950
331 gross 1023 bool
332     loadConfigured()
333     {
334     #ifdef USE_NETCDF
335     return true;
336     #else
337     return false;
338     #endif
339 gross 950 }
340    
341 jgs 82 Data
342     convertToData(const boost::python::object& value,
343     const FunctionSpace& what)
344     {
345     // first we try to extract a Data object from value
346     extract<Data> value_data(value);
347     if (value_data.check()) {
348     Data extracted_data=value_data();
349     if (extracted_data.isEmpty()) {
350     return extracted_data;
351     } else {
352     return Data(extracted_data,what);
353     }
354     } else {
355     return Data(value,what);
356     }
357     }
358    
359     } // end of namespace

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26