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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2105 - (show annotations)
Fri Nov 28 01:52:12 2008 UTC (10 years, 9 months ago) by jfenwick
File size: 13154 byte(s)
Data::copySelf() now returns an object instead of a pointer.
Fixed a bug in copyFromArray relating to expanded data.
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 "Data.h"
16 #include "DataConstant.h"
17 #include "DataException.h"
18 #include "esysUtils/EsysAssert.h"
19
20 #include <iostream>
21 #include <boost/python/extract.hpp>
22 #ifdef USE_NETCDF
23 #include <netcdfcpp.h>
24 #endif
25 #ifdef PASO_MPI
26 #include <mpi.h>
27 #endif
28
29 #include <boost/python/extract.hpp>
30 #include "DataMaths.h"
31
32 using namespace std;
33 using namespace boost::python;
34
35 namespace escript {
36
37 DataConstant::DataConstant(const boost::python::numeric::array& value,
38 const FunctionSpace& what)
39 : parent(what,DataTypes::shapeFromNumArray(value))
40 {
41 // // extract the shape of the numarray
42 // DataTypes::ShapeType tempShape;
43 // for (int i=0; i < value.getrank(); i++) {
44 // tempShape.push_back(extract<int>(value.getshape()[i]));
45 // }
46
47 // get the space for the data vector
48 // int len = getNoValues();
49 // DataVector temp_data(len, 0.0, len);
50 // DataArrayView temp_dataView(temp_data, tempShape);
51 // temp_dataView.copy(value);
52
53 m_data.copyFromNumArray(value,1);
54 //
55
56 // copy the data in the correct format
57 // m_data=temp_data;
58 //
59 // create the view of the data
60 // DataArrayView tempView(m_data,temp_dataView.getShape());
61 // setPointDataView(tempView);
62 }
63
64 // DataConstant::DataConstant(const DataArrayView& value,
65 // const FunctionSpace& what)
66 // : DataAbstract(what)
67 // {
68 // //
69 // // copy the data in the correct format
70 // m_data=value.getData();
71 // //
72 // // create the view of the data
73 // DataArrayView tempView(m_data,value.getShape());
74 // setPointDataView(tempView);
75 // }
76
77 DataConstant::DataConstant(const DataConstant& other)
78 : parent(other.getFunctionSpace(),other.getShape())
79 { //
80 // copy the data in the correct format
81 m_data=other.m_data;
82 //
83 // // create the view of the data
84 // DataArrayView tempView(m_data,other.getPointDataView().getShape());
85 // setPointDataView(tempView);
86 }
87
88 DataConstant::DataConstant(const DataConstant& other,
89 const DataTypes::RegionType& region)
90 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
91 {
92 //
93 // get the shape of the slice to copy from
94 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
95 //
96 // allocate space for this new DataConstant's data
97 int len = getNoValues();
98 m_data.resize(len,0.,len);
99 //
100 // create a view of the data with the correct shape
101 // DataArrayView tempView(m_data,shape);
102 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
103 //
104 // load the view with the data from the slice
105 // tempView.copySlice(other.getPointDataView(),region_loop_range);
106 DataTypes::copySlice(m_data,getShape(),0,other.getVector(),other.getShape(),0,region_loop_range);
107 // setPointDataView(tempView);
108 }
109
110 DataConstant::DataConstant(const FunctionSpace& what,
111 const DataTypes::ShapeType &shape,
112 const DataTypes::ValueType &data)
113 : parent(what,shape)
114 {
115 //
116 // copy the data in the correct format
117 m_data=data;
118 //
119 // create the view of the data
120 // DataArrayView tempView(m_data,shape);
121 // setPointDataView(tempView);
122 }
123
124 string
125 DataConstant::toString() const
126 {
127 return DataTypes::pointToString(m_data,getShape(),0,"");
128 }
129
130
131 DataAbstract*
132 DataConstant::deepCopy()
133 {
134 return new DataConstant(*this);
135 }
136
137
138 DataTypes::ValueType::size_type
139 DataConstant::getPointOffset(int sampleNo,
140 int dataPointNo) const
141 {
142 EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
143 "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
144 //
145 // Whatever the coord's always return the same value as this is constant data.
146 return 0;
147 }
148
149 DataTypes::ValueType::size_type
150 DataConstant::getPointOffset(int sampleNo,
151 int dataPointNo)
152 {
153 EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
154 "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
155 //
156 // Whatever the coord's always return the same value as this is constant data.
157 return 0;
158 }
159
160 DataTypes::ValueType::size_type
161 DataConstant::getLength() const
162 {
163 return m_data.size();
164 }
165
166 // DataArrayView
167 // DataConstant::getDataPoint(int sampleNo,
168 // int dataPointNo)
169 // {
170 // EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
171 // "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
172 // //
173 // // Whatever the coord's always return the same value as this is constant data.
174 // return getPointDataView();
175 // }
176
177 DataAbstract*
178 DataConstant::getSlice(const DataTypes::RegionType& region) const
179 {
180 return new DataConstant(*this,region);
181 }
182
183 void
184 DataConstant::setSlice(const DataAbstract* value,
185 const DataTypes::RegionType& region)
186 {
187 const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value);
188 if (tempDataConst==0) {
189 throw DataException("Programming error - casting to DataConstant.");
190 }
191 //
192 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
193 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
194 //
195 // check shape:
196 if (getRank()!=region.size()) {
197 throw DataException("Error - Invalid slice region.");
198 }
199 if (getRank()>0 && !DataTypes::checkShape(value->getShape(),shape)) {
200 throw DataException (DataTypes::createShapeErrorMessage(
201 "Error - Couldn't copy slice due to shape mismatch.",shape,value->getShape()));
202 }
203 // getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region_loop_range);
204 DataTypes::copySliceFrom(m_data,getShape(),0,tempDataConst->getVector(), tempDataConst->getShape(),0,region_loop_range);
205 }
206
207
208
209 void
210 DataConstant::symmetric(DataAbstract* ev)
211 {
212 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
213 if (temp_ev==0) {
214 throw DataException("Error - DataConstant::symmetric: casting to DataConstant failed (propably a programming error).");
215 }
216 /* DataArrayView& thisView=getPointDataView();
217 DataArrayView& evView=ev->getPointDataView();*/
218 DataMaths::symmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
219 }
220
221 void
222 DataConstant::nonsymmetric(DataAbstract* ev)
223 {
224 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
225 if (temp_ev==0) {
226 throw DataException("Error - DataConstant::nonsymmetric: casting to DataConstant failed (propably a programming error).");
227 }
228 /* DataArrayView& thisView=getPointDataView();
229 DataArrayView& evView=ev->getPointDataView();*/
230 DataMaths::nonsymmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
231 }
232
233 void
234 DataConstant::trace(DataAbstract* ev, int axis_offset)
235 {
236 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
237 if (temp_ev==0) {
238 throw DataException("Error - DataConstant::trace: casting to DataConstant failed (propably a programming error).");
239 }
240 /* DataArrayView& thisView=getPointDataView();
241 DataArrayView& evView=ev->getPointDataView();*/
242 ValueType& evVec=temp_ev->getVector();
243 const ShapeType& evShape=temp_ev->getShape();
244 DataMaths::trace(m_data,getShape(),0,evVec,evShape,0,axis_offset);
245 }
246
247 void
248 DataConstant::swapaxes(DataAbstract* ev, int axis0, int axis1)
249 {
250 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
251 if (temp_ev==0) {
252 throw DataException("Error - DataConstant::swapaxes: casting to DataConstant failed (propably a programming error).");
253 }
254 // DataArrayView& thisView=getPointDataView();
255 // DataArrayView& evView=ev->getPointDataView();
256 DataMaths::swapaxes(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,axis0,axis1);
257 }
258
259 void
260 DataConstant::transpose(DataAbstract* ev, int axis_offset)
261 {
262 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
263 if (temp_ev==0) {
264 throw DataException("Error - DataConstant::transpose: casting to DataConstant failed (propably a programming error).");
265 }
266 /* DataArrayView& thisView=getPointDataView();
267 DataArrayView& evView=ev->getPointDataView();*/
268 DataMaths::transpose(m_data, getShape(),0, temp_ev->getVector(),temp_ev->getShape(),0,axis_offset);
269 }
270
271 void
272 DataConstant::eigenvalues(DataAbstract* ev)
273 {
274 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
275 if (temp_ev==0) {
276 throw DataException("Error - DataConstant::eigenvalues: casting to DataConstant failed (propably a programming error).");
277 }
278 /* DataArrayView& thisView=getPointDataView();
279 DataArrayView& evView=ev->getPointDataView();*/
280 DataMaths::eigenvalues(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
281 }
282 void
283 DataConstant::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
284 {
285 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
286 if (temp_ev==0) {
287 throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
288 }
289 DataConstant* temp_V=dynamic_cast<DataConstant*>(V);
290 if (temp_V==0) {
291 throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
292 }
293 // DataArrayView thisView=getPointDataView();
294 // DataArrayView evView=ev->getPointDataView();
295 // DataArrayView VView=V->getPointDataView();
296
297 DataMaths::eigenvalues_and_eigenvectors(m_data, getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,temp_V->getVector(), temp_V->getShape(),0,tol);
298 }
299
300 void
301 DataConstant::setToZero()
302 {
303 DataTypes::ValueType::size_type n=m_data.size();
304 for (int i=0; i<n ;++i) m_data[i]=0.;
305 }
306
307 void
308 DataConstant::dump(const std::string fileName) const
309 {
310 #ifdef USE_NETCDF
311 const NcDim* ncdims[DataTypes::maxRank];
312 NcVar* var;
313 int rank = getRank();
314 int type= getFunctionSpace().getTypeCode();
315 int ndims =0;
316 long dims[DataTypes::maxRank];
317 const double* d_ptr=&(m_data[0]);
318 DataTypes::ShapeType shape = getShape();
319 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
320 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
321 #ifdef PASO_MPI
322 MPI_Status status;
323 #endif
324
325 #ifdef PASO_MPI
326 /* Serialize NetCDF I/O */
327 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81802, MPI_COMM_WORLD, &status);
328 #endif
329
330 // netCDF error handler
331 NcError err(NcError::verbose_nonfatal);
332 // Create the file.
333 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
334 NcFile dataFile(newFileName, NcFile::Replace);
335 // check if writing was successful
336 if (!dataFile.is_valid())
337 throw DataException("Error - DataConstant:: opening of netCDF file for output failed.");
338 if (!dataFile.add_att("type_id",0) )
339 throw DataException("Error - DataConstant:: appending data type to netCDF file failed.");
340 if (!dataFile.add_att("rank",rank) )
341 throw DataException("Error - DataConstant:: appending rank attribute to netCDF file failed.");
342 if (!dataFile.add_att("function_space_type",type))
343 throw DataException("Error - DataConstant:: appending function space attribute to netCDF file failed.");
344
345 if (rank == 0) {
346 if( ! (ncdims[0] = dataFile.add_dim("l", 1)) )
347 throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
348 dims[0]=1,
349 ndims=1;
350 } else {
351 ndims=rank;
352 dims[0]=shape[0];
353 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
354 throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
355 if ( rank >1 ) {
356 dims[1]=shape[1];
357 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
358 throw DataException("Error - DataConstant:: appending ncdimension 1 to netCDF file failed.");
359 }
360 if ( rank >2 ) {
361 dims[2]=shape[2];
362 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
363 throw DataException("Error - DataConstant:: appending ncdimension 2 to netCDF file failed.");
364 }
365 if ( rank >3 ) {
366 dims[3]=shape[3];
367 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
368 throw DataException("Error - DataConstant:: appending ncdimension 3 to netCDF file failed.");
369 }
370 }
371
372 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
373 throw DataException("Error - DataConstant:: appending variable to netCDF file failed.");
374 if (! (var->put(d_ptr,dims)) )
375 throw DataException("Error - DataConstant:: copy data to netCDF buffer failed.");
376 #ifdef PASO_MPI
377 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81802, MPI_COMM_WORLD);
378 #endif
379 #else
380 throw DataException("Error - DataConstant:: dump is not configured with netCDF. Please contact your installation manager.");
381 #endif
382 }
383
384 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26