/[escript]/trunk-mpi-branch/escript/src/DataFactory.cpp
ViewVC logotype

Annotation of /trunk-mpi-branch/escript/src/DataFactory.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1306 - (hide annotations)
Tue Sep 18 05:51:09 2007 UTC (12 years ago) by ksteube
File size: 11152 byte(s)
New Copyright in each .c .h .cpp and .py file

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26