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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1796 - (hide annotations)
Wed Sep 17 01:45:46 2008 UTC (11 years, 1 month ago) by jfenwick
File size: 12276 byte(s)
Merged noarrayview branch onto trunk.


1 jgs 82
2 ksteube 1312 /* $Id$ */
3    
4     /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
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 jgs 82
26 matt 1319 #include <boost/python/extract.hpp>
27 jfenwick 1796 #include "DataMaths.h"
28 matt 1319
29 jgs 82 using namespace std;
30 matt 1319 using namespace boost::python;
31 jgs 82
32     namespace escript {
33    
34 jgs 102 DataConstant::DataConstant(const boost::python::numeric::array& value,
35     const FunctionSpace& what)
36 jfenwick 1796 : DataAbstract(what,DataTypes::shapeFromNumArray(value))
37 jgs 102 {
38 jfenwick 1796 // // extract the shape of the numarray
39     // DataTypes::ShapeType tempShape;
40     // for (int i=0; i < value.getrank(); i++) {
41     // tempShape.push_back(extract<int>(value.getshape()[i]));
42     // }
43 matt 1319
44     // get the space for the data vector
45 jfenwick 1796 // int len = getNoValues();
46     // DataVector temp_data(len, 0.0, len);
47     // DataArrayView temp_dataView(temp_data, tempShape);
48     // temp_dataView.copy(value);
49 matt 1319
50 jfenwick 1796 m_data.copyFromNumArray(value);
51 jgs 102 //
52 jgs 82
53 jgs 102 // copy the data in the correct format
54 jfenwick 1796 // m_data=temp_data;
55 jgs 119 //
56     // create the view of the data
57 jfenwick 1796 // DataArrayView tempView(m_data,temp_dataView.getShape());
58     // setPointDataView(tempView);
59 jgs 102 }
60 jgs 82
61 jfenwick 1796 // DataConstant::DataConstant(const DataArrayView& value,
62     // const FunctionSpace& what)
63     // : DataAbstract(what)
64     // {
65     // //
66     // // copy the data in the correct format
67     // m_data=value.getData();
68     // //
69     // // create the view of the data
70     // DataArrayView tempView(m_data,value.getShape());
71     // setPointDataView(tempView);
72     // }
73    
74 jgs 102 DataConstant::DataConstant(const DataConstant& other)
75 jfenwick 1796 : DataAbstract(other.getFunctionSpace(),other.getShape())
76 matt 1319 { //
77 jgs 102 // copy the data in the correct format
78     m_data=other.m_data;
79 jgs 119 //
80 jfenwick 1796 // // create the view of the data
81     // DataArrayView tempView(m_data,other.getPointDataView().getShape());
82     // setPointDataView(tempView);
83 jgs 102 }
84 jgs 82
85 jgs 102 DataConstant::DataConstant(const DataConstant& other,
86 jfenwick 1796 const DataTypes::RegionType& region)
87     : DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
88 jgs 102 {
89     //
90     // get the shape of the slice to copy from
91 jfenwick 1796 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
92 jgs 102 //
93     // allocate space for this new DataConstant's data
94 jfenwick 1796 int len = getNoValues();
95 jgs 151 m_data.resize(len,0.,len);
96 jgs 102 //
97     // create a view of the data with the correct shape
98 jfenwick 1796 // DataArrayView tempView(m_data,shape);
99     DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
100 jgs 119 //
101     // load the view with the data from the slice
102 jfenwick 1796 // tempView.copySlice(other.getPointDataView(),region_loop_range);
103     DataTypes::copySlice(m_data,getShape(),0,other.getVector(),other.getShape(),0,region_loop_range);
104     // setPointDataView(tempView);
105 jgs 119 }
106    
107     DataConstant::DataConstant(const FunctionSpace& what,
108 jfenwick 1796 const DataTypes::ShapeType &shape,
109     const DataTypes::ValueType &data)
110     : DataAbstract(what,shape)
111 jgs 119 {
112 jgs 102 //
113 jgs 119 // copy the data in the correct format
114     m_data=data;
115     //
116     // create the view of the data
117 jfenwick 1796 // DataArrayView tempView(m_data,shape);
118     // setPointDataView(tempView);
119 jgs 102 }
120 jgs 82
121 jgs 102 string
122     DataConstant::toString() const
123     {
124 jfenwick 1796 return DataTypes::pointToString(m_data,getShape(),0,"");
125 jgs 102 }
126 jgs 82
127 jfenwick 1796 DataTypes::ValueType::size_type
128 jgs 102 DataConstant::getPointOffset(int sampleNo,
129     int dataPointNo) const
130     {
131     EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
132     "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
133 jgs 122 //
134     // Whatever the coord's always return the same value as this is constant data.
135 jgs 102 return 0;
136     }
137 jgs 82
138 jfenwick 1796 DataTypes::ValueType::size_type
139 jgs 102 DataConstant::getLength() const
140     {
141     return m_data.size();
142     }
143 jgs 82
144 jfenwick 1796 // DataArrayView
145     // DataConstant::getDataPoint(int sampleNo,
146     // int dataPointNo)
147     // {
148     // EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
149     // "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
150     // //
151     // // Whatever the coord's always return the same value as this is constant data.
152     // return getPointDataView();
153     // }
154 matt 1319
155 jgs 102 DataAbstract*
156 jfenwick 1796 DataConstant::getSlice(const DataTypes::RegionType& region) const
157 jgs 102 {
158     return new DataConstant(*this,region);
159     }
160 jgs 82
161 jgs 102 void
162     DataConstant::setSlice(const DataAbstract* value,
163 jfenwick 1796 const DataTypes::RegionType& region)
164 jgs 102 {
165     const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value);
166     if (tempDataConst==0) {
167     throw DataException("Programming error - casting to DataConstant.");
168 jgs 82 }
169 matt 1319 //
170 jfenwick 1796 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
171     DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
172 jgs 108 //
173     // check shape:
174 jfenwick 1796 if (getRank()!=region.size()) {
175 jgs 108 throw DataException("Error - Invalid slice region.");
176     }
177 jfenwick 1796 if (getRank()>0 && !DataTypes::checkShape(value->getShape(),shape)) {
178     throw DataException (DataTypes::createShapeErrorMessage(
179     "Error - Couldn't copy slice due to shape mismatch.",shape,value->getShape()));
180 jgs 108 }
181 jfenwick 1796 // getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region_loop_range);
182     DataTypes::copySliceFrom(m_data,getShape(),0,tempDataConst->getVector(), tempDataConst->getShape(),0,region_loop_range);
183 jgs 102 }
184 jgs 82
185 jgs 123
186    
187 gross 580 void
188 ksteube 775 DataConstant::symmetric(DataAbstract* ev)
189     {
190     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
191     if (temp_ev==0) {
192     throw DataException("Error - DataConstant::symmetric: casting to DataConstant failed (propably a programming error).");
193     }
194 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
195     DataArrayView& evView=ev->getPointDataView();*/
196     DataMaths::symmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
197 ksteube 775 }
198    
199     void
200     DataConstant::nonsymmetric(DataAbstract* ev)
201     {
202     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
203     if (temp_ev==0) {
204     throw DataException("Error - DataConstant::nonsymmetric: casting to DataConstant failed (propably a programming error).");
205     }
206 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
207     DataArrayView& evView=ev->getPointDataView();*/
208     DataMaths::nonsymmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
209 ksteube 775 }
210    
211     void
212 gross 800 DataConstant::trace(DataAbstract* ev, int axis_offset)
213 ksteube 775 {
214     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
215     if (temp_ev==0) {
216 gross 800 throw DataException("Error - DataConstant::trace: casting to DataConstant failed (propably a programming error).");
217 ksteube 775 }
218 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
219     DataArrayView& evView=ev->getPointDataView();*/
220     ValueType& evVec=temp_ev->getVector();
221     const ShapeType& evShape=temp_ev->getShape();
222     DataMaths::trace(m_data,getShape(),0,evVec,evShape,0,axis_offset);
223 ksteube 775 }
224    
225     void
226 gross 804 DataConstant::swapaxes(DataAbstract* ev, int axis0, int axis1)
227 gross 800 {
228     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
229     if (temp_ev==0) {
230 gross 804 throw DataException("Error - DataConstant::swapaxes: casting to DataConstant failed (propably a programming error).");
231 gross 800 }
232 jfenwick 1796 // DataArrayView& thisView=getPointDataView();
233     // DataArrayView& evView=ev->getPointDataView();
234     DataMaths::swapaxes(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,axis0,axis1);
235 gross 800 }
236    
237     void
238 ksteube 775 DataConstant::transpose(DataAbstract* ev, int axis_offset)
239     {
240     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
241     if (temp_ev==0) {
242     throw DataException("Error - DataConstant::transpose: casting to DataConstant failed (propably a programming error).");
243     }
244 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
245     DataArrayView& evView=ev->getPointDataView();*/
246     DataMaths::transpose(m_data, getShape(),0, temp_ev->getVector(),temp_ev->getShape(),0,axis_offset);
247 ksteube 775 }
248    
249     void
250 gross 580 DataConstant::eigenvalues(DataAbstract* ev)
251     {
252     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
253     if (temp_ev==0) {
254     throw DataException("Error - DataConstant::eigenvalues: casting to DataConstant failed (propably a programming error).");
255     }
256 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
257     DataArrayView& evView=ev->getPointDataView();*/
258     DataMaths::eigenvalues(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
259 gross 580 }
260     void
261     DataConstant::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
262     {
263     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
264     if (temp_ev==0) {
265     throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
266     }
267     DataConstant* temp_V=dynamic_cast<DataConstant*>(V);
268     if (temp_V==0) {
269     throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
270     }
271 jfenwick 1796 // DataArrayView thisView=getPointDataView();
272     // DataArrayView evView=ev->getPointDataView();
273     // DataArrayView VView=V->getPointDataView();
274 gross 580
275 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);
276 gross 580 }
277    
278 gross 950 void
279 gross 1118 DataConstant::setToZero()
280     {
281 jfenwick 1796 DataTypes::ValueType::size_type n=m_data.size();
282 gross 1118 for (int i=0; i<n ;++i) m_data[i]=0.;
283     }
284    
285     void
286 gross 950 DataConstant::dump(const std::string fileName) const
287     {
288     #ifdef PASO_MPI
289 ksteube 1312 throw DataException("Error - DataConstant:: dump is not implemented for MPI yet.");
290 gross 950 #endif
291 gross 1149 #ifdef USE_NETCDF
292 jfenwick 1796 const NcDim* ncdims[DataTypes::maxRank];
293 gross 950 NcVar* var;
294 jfenwick 1796 int rank = getRank();
295 gross 950 int type= getFunctionSpace().getTypeCode();
296     int ndims =0;
297 jfenwick 1796 long dims[DataTypes::maxRank];
298 gross 1141 const double* d_ptr=&(m_data[0]);
299 jfenwick 1796 DataTypes::ShapeType shape = getShape();
300 matt 1319
301 gross 950 // netCDF error handler
302     NcError err(NcError::verbose_nonfatal);
303     // Create the file.
304     NcFile dataFile(fileName.c_str(), NcFile::Replace);
305     // check if writing was successful
306 matt 1319 if (!dataFile.is_valid())
307 gross 950 throw DataException("Error - DataConstant:: opening of netCDF file for output failed.");
308 matt 1319 if (!dataFile.add_att("type_id",0) )
309 gross 950 throw DataException("Error - DataConstant:: appending data type to netCDF file failed.");
310 matt 1319 if (!dataFile.add_att("rank",rank) )
311 gross 950 throw DataException("Error - DataConstant:: appending rank attribute to netCDF file failed.");
312 matt 1319 if (!dataFile.add_att("function_space_type",type))
313 gross 950 throw DataException("Error - DataConstant:: appending function space attribute to netCDF file failed.");
314 gross 580
315 gross 950 if (rank == 0) {
316 matt 1319 if( ! (ncdims[0] = dataFile.add_dim("l", 1)) )
317 gross 950 throw DataException("Error - DataConstant:: appending ncdimsion 0 to netCDF file failed.");
318     dims[0]=1,
319     ndims=1;
320     } else {
321     ndims=rank;
322     dims[0]=shape[0];
323 matt 1319 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
324 gross 950 throw DataException("Error - DataConstant:: appending ncdimsion 0 to netCDF file failed.");
325     if ( rank >1 ) {
326     dims[1]=shape[1];
327 matt 1319 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
328 gross 950 throw DataException("Error - DataConstant:: appending ncdimsion 1 to netCDF file failed.");
329     }
330     if ( rank >2 ) {
331     dims[2]=shape[2];
332 matt 1319 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
333 gross 950 throw DataException("Error - DataConstant:: appending ncdimsion 2 to netCDF file failed.");
334     }
335     if ( rank >3 ) {
336     dims[3]=shape[3];
337 matt 1319 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
338 gross 950 throw DataException("Error - DataConstant:: appending ncdimsion 3 to netCDF file failed.");
339     }
340     }
341 matt 1319
342 gross 950 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
343     throw DataException("Error - DataConstant:: appending variable to netCDF file failed.");
344 gross 1141 if (! (var->put(d_ptr,dims)) )
345 matt 1319 throw DataException("Error - DataConstant:: copy data to netCDF buffer failed.");
346 gross 1023 #else
347     throw DataException("Error - DataConstant:: dump is not configured with netCDF. Please contact your installation manager.");
348     #endif
349 gross 950 }
350    
351 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