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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2005 - (hide 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 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
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 ksteube 1312
14 ksteube 1811
15 ksteube 1827 #include "Data.h"
16 jgs 474 #include "DataConstant.h"
17     #include "DataException.h"
18 robwdcock 682 #include "esysUtils/EsysAssert.h"
19 jgs 82
20     #include <iostream>
21     #include <boost/python/extract.hpp>
22 gross 1023 #ifdef USE_NETCDF
23 ksteube 1312 #include <netcdfcpp.h>
24 gross 1023 #endif
25 ksteube 1827 #ifdef PASO_MPI
26     #include <mpi.h>
27     #endif
28 jgs 82
29 matt 1319 #include <boost/python/extract.hpp>
30 jfenwick 1796 #include "DataMaths.h"
31 matt 1319
32 jgs 82 using namespace std;
33 matt 1319 using namespace boost::python;
34 jgs 82
35     namespace escript {
36    
37 jgs 102 DataConstant::DataConstant(const boost::python::numeric::array& value,
38     const FunctionSpace& what)
39 jfenwick 2005 : parent(what,DataTypes::shapeFromNumArray(value))
40 jgs 102 {
41 jfenwick 1796 // // 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 matt 1319
47     // get the space for the data vector
48 jfenwick 1796 // int len = getNoValues();
49     // DataVector temp_data(len, 0.0, len);
50     // DataArrayView temp_dataView(temp_data, tempShape);
51     // temp_dataView.copy(value);
52 matt 1319
53 jfenwick 1796 m_data.copyFromNumArray(value);
54 jgs 102 //
55 jgs 82
56 jgs 102 // copy the data in the correct format
57 jfenwick 1796 // m_data=temp_data;
58 jgs 119 //
59     // create the view of the data
60 jfenwick 1796 // DataArrayView tempView(m_data,temp_dataView.getShape());
61     // setPointDataView(tempView);
62 jgs 102 }
63 jgs 82
64 jfenwick 1796 // 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 jgs 102 DataConstant::DataConstant(const DataConstant& other)
78 jfenwick 2005 : parent(other.getFunctionSpace(),other.getShape())
79 matt 1319 { //
80 jgs 102 // copy the data in the correct format
81     m_data=other.m_data;
82 jgs 119 //
83 jfenwick 1796 // // create the view of the data
84     // DataArrayView tempView(m_data,other.getPointDataView().getShape());
85     // setPointDataView(tempView);
86 jgs 102 }
87 jgs 82
88 jgs 102 DataConstant::DataConstant(const DataConstant& other,
89 jfenwick 1796 const DataTypes::RegionType& region)
90 jfenwick 2005 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
91 jgs 102 {
92     //
93     // get the shape of the slice to copy from
94 jfenwick 1796 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
95 jgs 102 //
96     // allocate space for this new DataConstant's data
97 jfenwick 1796 int len = getNoValues();
98 jgs 151 m_data.resize(len,0.,len);
99 jgs 102 //
100     // create a view of the data with the correct shape
101 jfenwick 1796 // DataArrayView tempView(m_data,shape);
102     DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
103 jgs 119 //
104     // load the view with the data from the slice
105 jfenwick 1796 // 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 jgs 119 }
109    
110     DataConstant::DataConstant(const FunctionSpace& what,
111 jfenwick 1796 const DataTypes::ShapeType &shape,
112     const DataTypes::ValueType &data)
113 jfenwick 2005 : parent(what,shape)
114 jgs 119 {
115 jgs 102 //
116 jgs 119 // copy the data in the correct format
117     m_data=data;
118     //
119     // create the view of the data
120 jfenwick 1796 // DataArrayView tempView(m_data,shape);
121     // setPointDataView(tempView);
122 jgs 102 }
123 jgs 82
124 jgs 102 string
125     DataConstant::toString() const
126     {
127 jfenwick 1796 return DataTypes::pointToString(m_data,getShape(),0,"");
128 jgs 102 }
129 jgs 82
130 jfenwick 1799
131     DataAbstract*
132     DataConstant::deepCopy()
133     {
134     return new DataConstant(*this);
135     }
136    
137    
138 jfenwick 1796 DataTypes::ValueType::size_type
139 jgs 102 DataConstant::getPointOffset(int sampleNo,
140     int dataPointNo) const
141     {
142     EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
143     "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
144 jgs 122 //
145     // Whatever the coord's always return the same value as this is constant data.
146 jgs 102 return 0;
147     }
148 jgs 82
149 jfenwick 1796 DataTypes::ValueType::size_type
150 jfenwick 2005 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 jgs 102 DataConstant::getLength() const
162     {
163     return m_data.size();
164     }
165 jgs 82
166 jfenwick 1796 // 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 matt 1319
177 jgs 102 DataAbstract*
178 jfenwick 1796 DataConstant::getSlice(const DataTypes::RegionType& region) const
179 jgs 102 {
180     return new DataConstant(*this,region);
181     }
182 jgs 82
183 jgs 102 void
184     DataConstant::setSlice(const DataAbstract* value,
185 jfenwick 1796 const DataTypes::RegionType& region)
186 jgs 102 {
187     const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value);
188     if (tempDataConst==0) {
189     throw DataException("Programming error - casting to DataConstant.");
190 jgs 82 }
191 matt 1319 //
192 jfenwick 1796 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
193     DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
194 jgs 108 //
195     // check shape:
196 jfenwick 1796 if (getRank()!=region.size()) {
197 jgs 108 throw DataException("Error - Invalid slice region.");
198     }
199 jfenwick 1796 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 jgs 108 }
203 jfenwick 1796 // getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region_loop_range);
204     DataTypes::copySliceFrom(m_data,getShape(),0,tempDataConst->getVector(), tempDataConst->getShape(),0,region_loop_range);
205 jgs 102 }
206 jgs 82
207 jgs 123
208    
209 gross 580 void
210 ksteube 775 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 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
217     DataArrayView& evView=ev->getPointDataView();*/
218     DataMaths::symmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
219 ksteube 775 }
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 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
229     DataArrayView& evView=ev->getPointDataView();*/
230     DataMaths::nonsymmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
231 ksteube 775 }
232    
233     void
234 gross 800 DataConstant::trace(DataAbstract* ev, int axis_offset)
235 ksteube 775 {
236     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
237     if (temp_ev==0) {
238 gross 800 throw DataException("Error - DataConstant::trace: casting to DataConstant failed (propably a programming error).");
239 ksteube 775 }
240 jfenwick 1796 /* 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 ksteube 775 }
246    
247     void
248 gross 804 DataConstant::swapaxes(DataAbstract* ev, int axis0, int axis1)
249 gross 800 {
250     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
251     if (temp_ev==0) {
252 gross 804 throw DataException("Error - DataConstant::swapaxes: casting to DataConstant failed (propably a programming error).");
253 gross 800 }
254 jfenwick 1796 // 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 gross 800 }
258    
259     void
260 ksteube 775 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 jfenwick 1796 /* 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 ksteube 775 }
270    
271     void
272 gross 580 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 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
279     DataArrayView& evView=ev->getPointDataView();*/
280     DataMaths::eigenvalues(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
281 gross 580 }
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 jfenwick 1796 // DataArrayView thisView=getPointDataView();
294     // DataArrayView evView=ev->getPointDataView();
295     // DataArrayView VView=V->getPointDataView();
296 gross 580
297 jfenwick 1796 DataMaths::eigenvalues_and_eigenvectors(m_data, getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,temp_V->getVector(), temp_V->getShape(),0,tol);
298 gross 580 }
299    
300 gross 950 void
301 gross 1118 DataConstant::setToZero()
302     {
303 jfenwick 1796 DataTypes::ValueType::size_type n=m_data.size();
304 gross 1118 for (int i=0; i<n ;++i) m_data[i]=0.;
305     }
306    
307     void
308 gross 950 DataConstant::dump(const std::string fileName) const
309     {
310 gross 1149 #ifdef USE_NETCDF
311 jfenwick 1796 const NcDim* ncdims[DataTypes::maxRank];
312 gross 950 NcVar* var;
313 jfenwick 1796 int rank = getRank();
314 gross 950 int type= getFunctionSpace().getTypeCode();
315     int ndims =0;
316 jfenwick 1796 long dims[DataTypes::maxRank];
317 gross 1141 const double* d_ptr=&(m_data[0]);
318 jfenwick 1796 DataTypes::ShapeType shape = getShape();
319 jfenwick 1872 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
320     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
321 ksteube 1827 #ifdef PASO_MPI
322     MPI_Status status;
323     #endif
324 matt 1319
325 ksteube 1827 #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 gross 950 // netCDF error handler
331     NcError err(NcError::verbose_nonfatal);
332     // Create the file.
333 ksteube 1827 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
334     NcFile dataFile(newFileName, NcFile::Replace);
335 gross 950 // check if writing was successful
336 matt 1319 if (!dataFile.is_valid())
337 gross 950 throw DataException("Error - DataConstant:: opening of netCDF file for output failed.");
338 matt 1319 if (!dataFile.add_att("type_id",0) )
339 gross 950 throw DataException("Error - DataConstant:: appending data type to netCDF file failed.");
340 matt 1319 if (!dataFile.add_att("rank",rank) )
341 gross 950 throw DataException("Error - DataConstant:: appending rank attribute to netCDF file failed.");
342 matt 1319 if (!dataFile.add_att("function_space_type",type))
343 gross 950 throw DataException("Error - DataConstant:: appending function space attribute to netCDF file failed.");
344 gross 580
345 gross 950 if (rank == 0) {
346 matt 1319 if( ! (ncdims[0] = dataFile.add_dim("l", 1)) )
347 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
348 gross 950 dims[0]=1,
349     ndims=1;
350     } else {
351     ndims=rank;
352     dims[0]=shape[0];
353 matt 1319 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
354 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
355 gross 950 if ( rank >1 ) {
356     dims[1]=shape[1];
357 matt 1319 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
358 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 1 to netCDF file failed.");
359 gross 950 }
360     if ( rank >2 ) {
361     dims[2]=shape[2];
362 matt 1319 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
363 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 2 to netCDF file failed.");
364 gross 950 }
365     if ( rank >3 ) {
366     dims[3]=shape[3];
367 matt 1319 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
368 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 3 to netCDF file failed.");
369 gross 950 }
370     }
371 matt 1319
372 gross 950 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
373     throw DataException("Error - DataConstant:: appending variable to netCDF file failed.");
374 gross 1141 if (! (var->put(d_ptr,dims)) )
375 matt 1319 throw DataException("Error - DataConstant:: copy data to netCDF buffer failed.");
376 ksteube 1827 #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 gross 1023 #else
380     throw DataException("Error - DataConstant:: dump is not configured with netCDF. Please contact your installation manager.");
381     #endif
382 gross 950 }
383    
384 jgs 82 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26