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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2005 - (show annotations)
Mon Nov 10 01:21:39 2008 UTC (10 years, 10 months ago) by jfenwick
File size: 13152 byte(s)
Bringing all changes across from schroedinger.
(Note this does not mean development is done, just that it will happen
on the trunk for now).
If anyone notices any problems please contact me.


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);
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