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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26