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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (10 years, 3 months ago) by ksteube
File size: 12386 byte(s)
Copyright updated in all files

1
2 /*******************************************************
3 *
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
14
15 #include "DataFactory.h"
16 #include "esysUtils/esys_malloc.h"
17
18 #include <boost/python/extract.hpp>
19 #include <iostream>
20 #include <exception>
21 #ifdef USE_NETCDF
22 #include <netcdfcpp.h>
23 #endif
24 #ifdef PASO_MPI
25 #include <mpi.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 Tensor(double value,
54 const FunctionSpace& what,
55 bool expanded)
56 {
57 DataTypes::ShapeType shape(2,what.getDomain().getDim());
58 return Data(value,shape,what,expanded);
59 }
60
61 Data
62 Tensor3(double value,
63 const FunctionSpace& what,
64 bool expanded)
65 {
66 DataTypes::ShapeType shape(3,what.getDomain().getDim());
67 return Data(value,shape,what,expanded);
68 }
69
70 Data
71 Tensor4(double value,
72 const FunctionSpace& what,
73 bool expanded)
74 {
75 DataTypes::ShapeType shape(4,what.getDomain().getDim());
76 return Data(value,shape,what,expanded);
77 }
78
79 Data
80 load(const std::string fileName,
81 const AbstractDomain& domain)
82 {
83 #ifdef USE_NETCDF
84 NcAtt *type_att, *rank_att, *function_space_type_att;
85 // netCDF error handler
86 NcError err(NcError::silent_nonfatal);
87 int mpi_iam=0, mpi_num=1;
88 // Create the file.
89 #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 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 FunctionSpace function_space=FunctionSpace(domain, function_space_type);
106 /* 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 if (rank<0 || rank>DataTypes::maxRank)
112 throw DataException("Error - load:: rank in escript netCDF file is greater than maximum rank.");
113 /* recover type attribute */
114 int type=-1;
115 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 esysUtils::free(type_str);
125 } 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 }
130 delete type_att;
131
132 /* recover dimension */
133 int ndims=dataFile.num_dims();
134 int ntags =0 , num_samples =0 , num_data_points_per_sample =0, d=0, len_data_point=1;
135 NcDim *d_dim, *tags_dim, *num_samples_dim, *num_data_points_per_sample_dim;
136 /* recover shape */
137 DataTypes::ShapeType shape;
138 long dims[DataTypes::maxRank+2];
139 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 len_data_point*=d;
146 }
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 len_data_point*=d;
154 }
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 len_data_point*=d;
162 }
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 len_data_point*=d;
170 }
171 /* recover stuff */
172 Data out;
173 NcVar *var, *ids_var, *tags_var;
174 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 if (! var->get(&(out.getDataAtOffset(out.getDataOffset(0,0))), dims) )
190 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 throw DataException("Error - load:: illegal number of dimensions for tagged data in netCDF file.");
195 if (! (tags_dim=dataFile.get_dim("num_tags")) )
196 throw DataException("Error - load:: unable to recover number of tags from netCDF file.");
197 ntags=tags_dim->size();
198 dims[rank]=ntags;
199 int *tags = (int *) esysUtils::malloc(ntags*sizeof(int));
200 if (! ( tags_var = dataFile.get_var("tags")) )
201 {
202 esysUtils::free(tags);
203 throw DataException("Error - load:: unable to find tags in netCDF file.");
204 }
205 if (! tags_var->get(tags, ntags) )
206 {
207 esysUtils::free(tags);
208 throw DataException("Error - load:: unable to recover tags from netCDF file.");
209 }
210
211 // Current Version
212 /* DataVector data(len_data_point * ntags, 0., len_data_point * ntags);
213 if (!(var = dataFile.get_var("data")))
214 {
215 esysUtils::free(tags);
216 throw DataException("Error - load:: unable to find data in netCDF file.");
217 }
218 if (! var->get(&(data[0]), dims) )
219 {
220 esysUtils::free(tags);
221 throw DataException("Error - load:: unable to recover data from netCDF file.");
222 }
223 out=Data(DataArrayView(data,shape,0),function_space);
224 for (int t=1; t<ntags; ++t) {
225 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 }
245 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 esysUtils::free(tags);
253 } 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 if ( ! (num_samples_dim = dataFile.get_dim("num_samples") ) )
258 throw DataException("Error - load:: unable to recover number of samples from netCDF file.");
259 num_samples = num_samples_dim->size();
260 if ( ! (num_data_points_per_sample_dim = dataFile.get_dim("num_data_points_per_sample") ) )
261 throw DataException("Error - load:: unable to recover number of data points per sample from netCDF file.");
262 num_data_points_per_sample=num_data_points_per_sample_dim->size();
263 // check shape:
264 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 if (num_samples==0) {
267 out = Data(0,shape,function_space,true);
268 }
269 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 }
322 } else {
323 throw DataException("Error - load:: unknown escript data type in netCDF file.");
324 }
325 return out;
326 #else
327 throw DataException("Error - load:: is not compiled with netCFD. Please contact your insatllation manager.");
328 #endif
329 }
330
331 bool
332 loadConfigured()
333 {
334 #ifdef USE_NETCDF
335 return true;
336 #else
337 return false;
338 #endif
339 }
340
341 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