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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4286 - (show annotations)
Thu Mar 7 04:28:11 2013 UTC (6 years, 7 months ago) by caltinay
File size: 14951 byte(s)
Assorted spelling fixes.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 #include "DataFactory.h"
18 #include "esysUtils/esys_malloc.h"
19 #include "esysUtils/Esys_MPI.h"
20
21 #include <boost/python/extract.hpp>
22 #include <iostream>
23 #include <exception>
24 #ifdef USE_NETCDF
25 #include <netcdfcpp.h>
26 #endif
27
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 DataTypes::ShapeType shape;
40 return Data(value,shape,what,expanded);
41 }
42
43 Data
44 Vector(double value,
45 const FunctionSpace& what,
46 bool expanded)
47 {
48 DataTypes::ShapeType shape(1,what.getDomain()->getDim());
49 return Data(value,shape,what,expanded);
50 }
51
52 Data
53 VectorFromObj(boost::python::object o,
54 const FunctionSpace& what,
55 bool expanded)
56 {
57 double v;
58 try // first try to get a double and route it to the other method
59 {
60 v=boost::python::extract<double>(o);
61 return Vector(v,what,expanded);
62 }
63 catch(...)
64 {
65 PyErr_Clear();
66 }
67 DataTypes::ShapeType shape(1,what.getDomain()->getDim());
68 Data d(o,what,expanded);
69 if (d.getDataPointShape()!=shape)
70 {
71 throw DataException("VectorFromObj: Shape of vector passed to function does not match the dimension of the domain. ");
72 }
73 return d;
74 }
75
76 Data
77 Tensor(double value,
78 const FunctionSpace& what,
79 bool expanded)
80 {
81 DataTypes::ShapeType shape(2,what.getDomain()->getDim());
82 return Data(value,shape,what,expanded);
83 }
84
85
86 // We need to take some care here because this signature trumps the other one from boost's point of view
87 Data
88 TensorFromObj(boost::python::object o,
89 const FunctionSpace& what,
90 bool expanded)
91 {
92 double v;
93 try // first try to get a double and route it to the other method
94 {
95 v=boost::python::extract<double>(o);
96 return Tensor(v,what,expanded);
97 }
98 catch(...)
99 {
100 PyErr_Clear();
101 }
102 DataTypes::ShapeType shape(2,what.getDomain()->getDim());
103 Data d(o,what,expanded);
104 if (d.getDataPointShape()!=shape)
105 {
106 throw DataException("TensorFromObj: Shape of tensor passed to function does not match the dimension of the domain. ");
107 }
108 return d;
109 }
110
111 Data
112 Tensor3(double value,
113 const FunctionSpace& what,
114 bool expanded)
115 {
116 DataTypes::ShapeType shape(3,what.getDomain()->getDim());
117 return Data(value,shape,what,expanded);
118 }
119
120 Data
121 Tensor3FromObj(boost::python::object o,
122 const FunctionSpace& what,
123 bool expanded)
124 {
125 double v;
126 try // first try to get a double and route it to the other method
127 {
128 v=boost::python::extract<double>(o);
129 return Tensor3(v,what,expanded);
130 }
131 catch(...)
132 {
133 PyErr_Clear();
134 }
135 DataTypes::ShapeType shape(3,what.getDomain()->getDim());
136 Data d(o,what,expanded);
137 if (d.getDataPointShape()!=shape)
138 {
139 throw DataException("Tensor3FromObj: Shape of tensor passed to function does not match the dimension of the domain. ");
140 }
141 return d;
142 }
143
144 Data
145 Tensor4(double value,
146 const FunctionSpace& what,
147 bool expanded)
148 {
149 DataTypes::ShapeType shape(4,what.getDomain()->getDim());
150 return Data(value,shape,what,expanded);
151 }
152
153 Data
154 Tensor4FromObj(boost::python::object o,
155 const FunctionSpace& what,
156 bool expanded)
157 {
158 double v;
159 try // first try to get a double and route it to the other method
160 {
161 v=boost::python::extract<double>(o);
162 return Tensor4(v,what,expanded);
163 }
164 catch(...)
165 {
166 PyErr_Clear();
167 }
168 DataTypes::ShapeType shape(4,what.getDomain()->getDim());
169 Data d(o,what,expanded);
170 if (d.getDataPointShape()!=shape)
171 {
172 throw DataException("VectorFromObj: Shape of tensor passed to function does not match the dimension of the domain. ");
173 }
174 return d;
175 }
176
177
178 Data
179 load(const std::string fileName,
180 const AbstractDomain& domain)
181 {
182 #ifdef USE_NETCDF
183 NcAtt *type_att, *rank_att, *function_space_type_att;
184 // netCDF error handler
185 NcError err(NcError::silent_nonfatal);
186 int mpi_iam=0, mpi_num=1;
187 // Create the file.
188 #ifdef ESYS_MPI
189 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_iam);
190 MPI_Comm_size(MPI_COMM_WORLD, &mpi_num);
191 #endif
192 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
193 NcFile dataFile(newFileName, NcFile::ReadOnly);
194 if (!dataFile.is_valid())
195 throw DataException("Error - load:: opening of netCDF file for input failed.");
196 /* recover function space */
197 if (! (function_space_type_att=dataFile.get_att("function_space_type")) )
198 throw DataException("Error - load:: cannot recover function_space_type attribute from escript netCDF file.");
199 int function_space_type = function_space_type_att->as_int(0);
200 delete function_space_type_att;
201 /* test if function space id is valid and create function space instance */
202 if (! domain.isValidFunctionSpaceType(function_space_type) )
203 throw DataException("Error - load:: function space type code in netCDF file is invalid for given domain.");
204 FunctionSpace function_space=FunctionSpace(domain.getPtr(), function_space_type);
205 /* recover rank */
206 if (! (rank_att=dataFile.get_att("rank")) )
207 throw DataException("Error - load:: cannot recover rank attribute from escript netCDF file.");
208 int rank = rank_att->as_int(0);
209 delete rank_att;
210 if (rank<0 || rank>DataTypes::maxRank)
211 throw DataException("Error - load:: rank in escript netCDF file is greater than maximum rank.");
212 /* recover type attribute */
213 int type=-1;
214 if ((type_att=dataFile.get_att("type")) ) {
215 char* type_str = type_att->as_string(0);
216 if (strncmp(type_str, "constant", strlen("constant")) == 0 ) {
217 type =0;
218 } else if (strncmp(type_str, "tagged", strlen("tagged")) == 0 ) {
219 type =1;
220 } else if (strncmp(type_str, "expanded", strlen("expanded")) == 0 ) {
221 type =2;
222 }
223 esysUtils::free(type_str);
224 } else {
225 if (! (type_att=dataFile.get_att("type_id")) )
226 throw DataException("Error - load:: cannot recover type attribute from escript netCDF file.");
227 type=type_att->as_int(0);
228 }
229 delete type_att;
230
231 /* recover dimension */
232 int ndims=dataFile.num_dims();
233 int ntags =0 , num_samples =0 , num_data_points_per_sample =0, d=0, len_data_point=1;
234 NcDim *d_dim, *tags_dim, *num_samples_dim, *num_data_points_per_sample_dim;
235 /* recover shape */
236 DataTypes::ShapeType shape;
237 long dims[DataTypes::maxRank+2];
238 if (rank>0) {
239 if (! (d_dim=dataFile.get_dim("d0")) )
240 throw DataException("Error - load:: unable to recover d0 from netCDF file.");
241 d=d_dim->size();
242 shape.push_back(d);
243 dims[0]=d;
244 len_data_point*=d;
245 }
246 if (rank>1) {
247 if (! (d_dim=dataFile.get_dim("d1")) )
248 throw DataException("Error - load:: unable to recover d1 from netCDF file.");
249 d=d_dim->size();
250 shape.push_back(d);
251 dims[1]=d;
252 len_data_point*=d;
253 }
254 if (rank>2) {
255 if (! (d_dim=dataFile.get_dim("d2")) )
256 throw DataException("Error - load:: unable to recover d2 from netCDF file.");
257 d=d_dim->size();
258 shape.push_back(d);
259 dims[2]=d;
260 len_data_point*=d;
261 }
262 if (rank>3) {
263 if (! (d_dim=dataFile.get_dim("d3")) )
264 throw DataException("Error - load:: unable to recover d3 from netCDF file.");
265 d=d_dim->size();
266 shape.push_back(d);
267 dims[3]=d;
268 len_data_point*=d;
269 }
270 /* recover stuff */
271 Data out;
272 NcVar *var, *ids_var, *tags_var;
273 if (type == 0) {
274 /* constant data */
275 if ( ! ( (ndims == rank && rank >0) || ( ndims ==1 && rank == 0 ) ) )
276 throw DataException("Error - load:: illegal number of dimensions for constant data in netCDF file.");
277 if (rank == 0) {
278 if (! (d_dim=dataFile.get_dim("l")) )
279 throw DataException("Error - load:: unable to recover d0 for scalar constant data in netCDF file.");
280 int d0 = d_dim->size();
281 if (! d0 == 1)
282 throw DataException("Error - load:: d0 is expected to be one for scalar constant data in netCDF file.");
283 dims[0]=1;
284 }
285 out=Data(0,shape,function_space);
286 if (!(var = dataFile.get_var("data")))
287 throw DataException("Error - load:: unable to find data in netCDF file.");
288 if (! var->get(&(out.getDataAtOffsetRW(out.getDataOffset(0,0))), dims) )
289 throw DataException("Error - load:: unable to recover data from netCDF file.");
290 } else if (type == 1) {
291 /* tagged data */
292 if ( ! (ndims == rank + 1) )
293 throw DataException("Error - load:: illegal number of dimensions for tagged data in netCDF file.");
294 if (! (tags_dim=dataFile.get_dim("num_tags")) )
295 throw DataException("Error - load:: unable to recover number of tags from netCDF file.");
296 ntags=tags_dim->size();
297 dims[rank]=ntags;
298 int *tags = (int *) esysUtils::malloc(ntags*sizeof(int));
299 if (! ( tags_var = dataFile.get_var("tags")) )
300 {
301 esysUtils::free(tags);
302 throw DataException("Error - load:: unable to find tags in netCDF file.");
303 }
304 if (! tags_var->get(tags, ntags) )
305 {
306 esysUtils::free(tags);
307 throw DataException("Error - load:: unable to recover tags from netCDF file.");
308 }
309
310 // Current Version
311 /* DataVector data(len_data_point * ntags, 0., len_data_point * ntags);
312 if (!(var = dataFile.get_var("data")))
313 {
314 esysUtils::free(tags);
315 throw DataException("Error - load:: unable to find data in netCDF file.");
316 }
317 if (! var->get(&(data[0]), dims) )
318 {
319 esysUtils::free(tags);
320 throw DataException("Error - load:: unable to recover data from netCDF file.");
321 }
322 out=Data(DataArrayView(data,shape,0),function_space);
323 for (int t=1; t<ntags; ++t) {
324 out.setTaggedValueFromCPP(tags[t],shape, data, t*len_data_point);
325 // out.setTaggedValueFromCPP(tags[t],DataArrayView(data,shape,t*len_data_point));
326 }*/
327 // End current version
328
329 // New version
330
331 // A) create a DataTagged dt
332 // B) Read data from file
333 // C) copy default value into dt
334 // D) copy tagged values into dt
335 // E) create a new Data based on dt
336
337 NcVar* var1;
338 DataVector data1(len_data_point * ntags, 0., len_data_point * ntags);
339 if (!(var1 = dataFile.get_var("data")))
340 {
341 esysUtils::free(tags);
342 throw DataException("Error - load:: unable to find data in netCDF file.");
343 }
344 if (! var1->get(&(data1[0]), dims) )
345 {
346 esysUtils::free(tags);
347 throw DataException("Error - load:: unable to recover data from netCDF file.");
348 }
349 DataTagged* dt=new DataTagged(function_space, shape, tags,data1);
350 out=Data(dt);
351 esysUtils::free(tags);
352 } else if (type == 2) {
353 /* expanded data */
354 if ( ! (ndims == rank + 2) )
355 throw DataException("Error - load:: illegal number of dimensions for expanded data in netCDF file.");
356 if ( ! (num_samples_dim = dataFile.get_dim("num_samples") ) )
357 throw DataException("Error - load:: unable to recover number of samples from netCDF file.");
358 num_samples = num_samples_dim->size();
359 if ( ! (num_data_points_per_sample_dim = dataFile.get_dim("num_data_points_per_sample") ) )
360 throw DataException("Error - load:: unable to recover number of data points per sample from netCDF file.");
361 num_data_points_per_sample=num_data_points_per_sample_dim->size();
362 // check shape:
363 if ( ! (num_samples == function_space.getNumSamples() && num_data_points_per_sample == function_space.getNumDataPointsPerSample()) )
364 throw DataException("Error - load:: data sample layout of file does not match data layout of function space.");
365 if (num_samples==0) {
366 out = Data(0,shape,function_space,true);
367 }
368 else {
369 // get ids
370 if (! ( ids_var = dataFile.get_var("id")) )
371 throw DataException("Error - load:: unable to find reference ids in netCDF file.");
372 const int* ids_p=function_space.borrowSampleReferenceIDs();
373 int *ids_of_nc = (int *)esysUtils::malloc(num_samples*sizeof(int));
374 if (! ids_var->get(ids_of_nc, (long) num_samples) )
375 {
376 esysUtils::free(ids_of_nc);
377 throw DataException("Error - load:: unable to recover ids from netCDF file.");
378 }
379 // check order:
380 int failed=-1, local_failed=-1, i;
381 #pragma omp parallel private(local_failed)
382 {
383 local_failed=-1;
384 #pragma omp for private(i) schedule(static)
385 for (i=0;i < num_samples; ++i) {
386 if (ids_of_nc[i]!=ids_p[i]) local_failed=i;
387 }
388 #pragma omp critical
389 if (local_failed>=0) failed = local_failed;
390 }
391 /* if (failed>=0)
392 {
393 esysUtils::free(ids_of_nc);
394 throw DataException("Error - load:: data ordering in netCDF file does not match ordering of FunctionSpace.");
395 } */
396 // get the data:
397 dims[rank]=num_data_points_per_sample;
398 dims[rank+1]=num_samples;
399 out=Data(0,shape,function_space,true);
400 if (!(var = dataFile.get_var("data")))
401 {
402 esysUtils::free(ids_of_nc);
403 throw DataException("Error - load:: unable to find data in netCDF file.");
404 }
405 if (! var->get(&(out.getDataAtOffsetRW(out.getDataOffset(0,0))), dims) )
406 {
407 esysUtils::free(ids_of_nc);
408 throw DataException("Error - load:: unable to recover data from netCDF file.");
409 }
410 if (failed>=0) {
411 try {
412 std::cout << "Information - load: start reordering data from netCDF file " << fileName << std::endl;
413 out.borrowData()->reorderByReferenceIDs(ids_of_nc);
414 }
415 catch (std::exception&) {
416 esysUtils::free(ids_of_nc);
417 throw DataException("Error - load:: unable to reorder data in netCDF file.");
418 }
419 }
420 }
421 } else {
422 throw DataException("Error - load:: unknown escript data type in netCDF file.");
423 }
424 return out;
425 #else
426 throw DataException("Error - load:: is not compiled with netCDF. Please contact your installation manager.");
427 #endif
428 }
429
430 bool
431 loadConfigured()
432 {
433 #ifdef USE_NETCDF
434 return true;
435 #else
436 return false;
437 #endif
438 }
439
440 Data
441 convertToData(const boost::python::object& value,
442 const FunctionSpace& what)
443 {
444 // first we try to extract a Data object from value
445 extract<Data> value_data(value);
446 if (value_data.check()) {
447 Data extracted_data=value_data();
448 if (extracted_data.isEmpty()) {
449 return extracted_data;
450 } else {
451 return Data(extracted_data,what);
452 }
453 } else {
454 return Data(value,what);
455 }
456 }
457
458 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26