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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1020 - (hide annotations)
Mon Mar 12 10:12:36 2007 UTC (12 years, 3 months ago) by phornby
File size: 10882 byte(s)
Added explicit destructors to all Exception classes.

Fixed an ifdef in TestCase.cpp

Made the conditional definition of M_PI in LocalOps.h
depend only on M_PI being undefined.

Replace dynamically dimensioned arrays in DataFactory & DataTagged with malloc.

sort() method of list does not take a named argument
(despite the manual's claims to the contary).


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26