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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 950 - (show annotations)
Tue Feb 6 07:01:11 2007 UTC (12 years, 7 months ago) by gross
File size: 11247 byte(s)
escript data objects can now be saved to netCDF files, see http://www.unidata.ucar.edu/software/netcdf/.
Currently only constant data are implemented with expanded and tagged data to follow.
There are two new functions to dump a data object

   s=Data(...)
   s.dump(<filename>)

and to recover it

   s=load(<filename>, domain)

Notice that the function space of s is recovered but domain is still need. 

dump and load will replace archive and extract.

The installation needs now the netCDF installed. 


1 //$Id$
2 /*
3 ************************************************************
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 */
13
14 #include "DataConstant.h"
15 #include "DataException.h"
16 #include "esysUtils/EsysAssert.h"
17
18 #include <iostream>
19 #include <boost/python/extract.hpp>
20 #include <netcdfcpp.h>
21
22 using namespace std;
23
24 namespace escript {
25
26 DataConstant::DataConstant(const boost::python::numeric::array& value,
27 const FunctionSpace& what)
28 : DataAbstract(what)
29 {
30 DataArray temp(value);
31 //
32 // copy the data in the correct format
33 m_data=temp.getData();
34 //
35 // create the view of the data
36 DataArrayView tempView(m_data,temp.getView().getShape());
37 setPointDataView(tempView);
38 }
39
40 DataConstant::DataConstant(const DataArrayView& value,
41 const FunctionSpace& what)
42 : DataAbstract(what)
43 {
44 //
45 // copy the data in the correct format
46 m_data=value.getData();
47 //
48 // create the view of the data
49 DataArrayView tempView(m_data,value.getShape());
50 setPointDataView(tempView);
51 }
52
53 DataConstant::DataConstant(const DataConstant& other)
54 : DataAbstract(other.getFunctionSpace())
55 {
56 //
57 // copy the data in the correct format
58 m_data=other.m_data;
59 //
60 // create the view of the data
61 DataArrayView tempView(m_data,other.getPointDataView().getShape());
62 setPointDataView(tempView);
63 }
64
65 DataConstant::DataConstant(const DataConstant& other,
66 const DataArrayView::RegionType& region)
67 : DataAbstract(other.getFunctionSpace())
68 {
69 //
70 // get the shape of the slice to copy from
71 DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
72 //
73 // allocate space for this new DataConstant's data
74 int len = DataArrayView::noValues(shape);
75 m_data.resize(len,0.,len);
76 //
77 // create a view of the data with the correct shape
78 DataArrayView tempView(m_data,shape);
79 DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
80 //
81 // load the view with the data from the slice
82 tempView.copySlice(other.getPointDataView(),region_loop_range);
83 setPointDataView(tempView);
84 }
85
86 DataConstant::DataConstant(const FunctionSpace& what,
87 const DataArrayView::ShapeType &shape,
88 const DataArrayView::ValueType &data)
89 : DataAbstract(what)
90 {
91 //
92 // copy the data in the correct format
93 m_data=data;
94 //
95 // create the view of the data
96 DataArrayView tempView(m_data,shape);
97 setPointDataView(tempView);
98 }
99
100 string
101 DataConstant::toString() const
102 {
103 return getPointDataView().toString("");
104 }
105
106 DataArrayView::ValueType::size_type
107 DataConstant::getPointOffset(int sampleNo,
108 int dataPointNo) const
109 {
110 EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
111 "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
112 //
113 // Whatever the coord's always return the same value as this is constant data.
114 return 0;
115 }
116
117 DataArrayView::ValueType::size_type
118 DataConstant::getLength() const
119 {
120 return m_data.size();
121 }
122
123 DataArrayView
124 DataConstant::getDataPoint(int sampleNo,
125 int dataPointNo)
126 {
127 EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
128 "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
129 //
130 // Whatever the coord's always return the same value as this is constant data.
131 return getPointDataView();
132 }
133
134 DataAbstract*
135 DataConstant::getSlice(const DataArrayView::RegionType& region) const
136 {
137 return new DataConstant(*this,region);
138 }
139
140 void
141 DataConstant::setSlice(const DataAbstract* value,
142 const DataArrayView::RegionType& region)
143 {
144 const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value);
145 if (tempDataConst==0) {
146 throw DataException("Programming error - casting to DataConstant.");
147 }
148 //
149 DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
150 DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
151 //
152 // check shape:
153 if (getPointDataView().getRank()!=region.size()) {
154 throw DataException("Error - Invalid slice region.");
155 }
156 if (tempDataConst->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {
157 throw DataException (value->getPointDataView().createShapeErrorMessage(
158 "Error - Couldn't copy slice due to shape mismatch.",shape));
159 }
160 //
161 getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region_loop_range);
162 }
163
164 int
165 DataConstant::archiveData(ofstream& archiveFile,
166 const DataArrayView::ValueType::size_type noValues) const
167 {
168 return(m_data.archiveData(archiveFile, noValues));
169 }
170
171 int
172 DataConstant::extractData(ifstream& archiveFile,
173 const DataArrayView::ValueType::size_type noValues)
174 {
175 return(m_data.extractData(archiveFile, noValues));
176 }
177
178 void
179 DataConstant::symmetric(DataAbstract* ev)
180 {
181 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
182 if (temp_ev==0) {
183 throw DataException("Error - DataConstant::symmetric: casting to DataConstant failed (propably a programming error).");
184 }
185 DataArrayView& thisView=getPointDataView();
186 DataArrayView& evView=ev->getPointDataView();
187 DataArrayView::symmetric(thisView,0,evView,0);
188 }
189
190 void
191 DataConstant::nonsymmetric(DataAbstract* ev)
192 {
193 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
194 if (temp_ev==0) {
195 throw DataException("Error - DataConstant::nonsymmetric: casting to DataConstant failed (propably a programming error).");
196 }
197 DataArrayView& thisView=getPointDataView();
198 DataArrayView& evView=ev->getPointDataView();
199 DataArrayView::nonsymmetric(thisView,0,evView,0);
200 }
201
202 void
203 DataConstant::trace(DataAbstract* ev, int axis_offset)
204 {
205 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
206 if (temp_ev==0) {
207 throw DataException("Error - DataConstant::trace: casting to DataConstant failed (propably a programming error).");
208 }
209 DataArrayView& thisView=getPointDataView();
210 DataArrayView& evView=ev->getPointDataView();
211 DataArrayView::trace(thisView,0,evView,0,axis_offset);
212 }
213
214 void
215 DataConstant::swapaxes(DataAbstract* ev, int axis0, int axis1)
216 {
217 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
218 if (temp_ev==0) {
219 throw DataException("Error - DataConstant::swapaxes: casting to DataConstant failed (propably a programming error).");
220 }
221 DataArrayView& thisView=getPointDataView();
222 DataArrayView& evView=ev->getPointDataView();
223 DataArrayView::swapaxes(thisView,0,evView,0,axis0,axis1);
224 }
225
226 void
227 DataConstant::transpose(DataAbstract* ev, int axis_offset)
228 {
229 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
230 if (temp_ev==0) {
231 throw DataException("Error - DataConstant::transpose: casting to DataConstant failed (propably a programming error).");
232 }
233 DataArrayView& thisView=getPointDataView();
234 DataArrayView& evView=ev->getPointDataView();
235 DataArrayView::transpose(thisView,0,evView,0,axis_offset);
236 }
237
238 void
239 DataConstant::eigenvalues(DataAbstract* ev)
240 {
241 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
242 if (temp_ev==0) {
243 throw DataException("Error - DataConstant::eigenvalues: casting to DataConstant failed (propably a programming error).");
244 }
245 DataArrayView& thisView=getPointDataView();
246 DataArrayView& evView=ev->getPointDataView();
247 DataArrayView::eigenvalues(thisView,0,evView,0);
248 }
249 void
250 DataConstant::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
251 {
252 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
253 if (temp_ev==0) {
254 throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
255 }
256 DataConstant* temp_V=dynamic_cast<DataConstant*>(V);
257 if (temp_V==0) {
258 throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
259 }
260 DataArrayView thisView=getPointDataView();
261 DataArrayView evView=ev->getPointDataView();
262 DataArrayView VView=V->getPointDataView();
263
264 DataArrayView::eigenvalues_and_eigenvectors(thisView,0,evView,0,VView,tol);
265 }
266
267 void
268 DataConstant::dump(const std::string fileName) const
269 {
270 #ifdef PASO_MPI
271 throw DataException("Error - DataConstant:: dump is not implemented for MPI yet.")
272 #endif
273 const NcDim* ncdims[DataArrayView::maxRank];
274 NcVar* var;
275 int rank = getPointDataView().getRank();
276 int type= getFunctionSpace().getTypeCode();
277 int ndims =0;
278 long dims[DataArrayView::maxRank];
279 DataArrayView::ShapeType shape = getPointDataView().getShape();
280
281 // netCDF error handler
282 NcError err(NcError::verbose_nonfatal);
283 // Create the file.
284 NcFile dataFile(fileName.c_str(), NcFile::Replace);
285 // check if writing was successful
286 if (!dataFile.is_valid())
287 throw DataException("Error - DataConstant:: opening of netCDF file for output failed.");
288 if (!dataFile.add_att("type","constant") )
289 throw DataException("Error - DataConstant:: appending data type to netCDF file failed.");
290 if (!dataFile.add_att("rank",rank) )
291 throw DataException("Error - DataConstant:: appending rank attribute to netCDF file failed.");
292 if (!dataFile.add_att("function_space_type",type))
293 throw DataException("Error - DataConstant:: appending function space attribute to netCDF file failed.");
294
295 if (rank == 0) {
296 if( ! (ncdims[0] = dataFile.add_dim("l", 1)) )
297 throw DataException("Error - DataConstant:: appending ncdimsion 0 to netCDF file failed.");
298 dims[0]=1,
299 ndims=1;
300 } else {
301 ndims=rank;
302 dims[0]=shape[0];
303 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
304 throw DataException("Error - DataConstant:: appending ncdimsion 0 to netCDF file failed.");
305 if ( rank >1 ) {
306 dims[1]=shape[1];
307 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
308 throw DataException("Error - DataConstant:: appending ncdimsion 1 to netCDF file failed.");
309 }
310 if ( rank >2 ) {
311 dims[2]=shape[2];
312 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
313 throw DataException("Error - DataConstant:: appending ncdimsion 2 to netCDF file failed.");
314 }
315 if ( rank >3 ) {
316 dims[3]=shape[3];
317 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
318 throw DataException("Error - DataConstant:: appending ncdimsion 3 to netCDF file failed.");
319 }
320 }
321 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
322 throw DataException("Error - DataConstant:: appending variable to netCDF file failed.");
323 if (! (var->put(&m_data[0],dims)) )
324 throw DataException("Error - DataConstant:: copy data to netCDF buffer failed.");
325 }
326
327 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26