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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1872 - (hide annotations)
Mon Oct 13 00:18:55 2008 UTC (11 years, 1 month ago) by jfenwick
File size: 12786 byte(s)
Closing the moreshared branch

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 1796 : DataAbstract(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 1796 : DataAbstract(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     : DataAbstract(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     : DataAbstract(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 jgs 102 DataConstant::getLength() const
151     {
152     return m_data.size();
153     }
154 jgs 82
155 jfenwick 1796 // DataArrayView
156     // DataConstant::getDataPoint(int sampleNo,
157     // int dataPointNo)
158     // {
159     // EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
160     // "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
161     // //
162     // // Whatever the coord's always return the same value as this is constant data.
163     // return getPointDataView();
164     // }
165 matt 1319
166 jgs 102 DataAbstract*
167 jfenwick 1796 DataConstant::getSlice(const DataTypes::RegionType& region) const
168 jgs 102 {
169     return new DataConstant(*this,region);
170     }
171 jgs 82
172 jgs 102 void
173     DataConstant::setSlice(const DataAbstract* value,
174 jfenwick 1796 const DataTypes::RegionType& region)
175 jgs 102 {
176     const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value);
177     if (tempDataConst==0) {
178     throw DataException("Programming error - casting to DataConstant.");
179 jgs 82 }
180 matt 1319 //
181 jfenwick 1796 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
182     DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
183 jgs 108 //
184     // check shape:
185 jfenwick 1796 if (getRank()!=region.size()) {
186 jgs 108 throw DataException("Error - Invalid slice region.");
187     }
188 jfenwick 1796 if (getRank()>0 && !DataTypes::checkShape(value->getShape(),shape)) {
189     throw DataException (DataTypes::createShapeErrorMessage(
190     "Error - Couldn't copy slice due to shape mismatch.",shape,value->getShape()));
191 jgs 108 }
192 jfenwick 1796 // getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region_loop_range);
193     DataTypes::copySliceFrom(m_data,getShape(),0,tempDataConst->getVector(), tempDataConst->getShape(),0,region_loop_range);
194 jgs 102 }
195 jgs 82
196 jgs 123
197    
198 gross 580 void
199 ksteube 775 DataConstant::symmetric(DataAbstract* ev)
200     {
201     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
202     if (temp_ev==0) {
203     throw DataException("Error - DataConstant::symmetric: casting to DataConstant failed (propably a programming error).");
204     }
205 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
206     DataArrayView& evView=ev->getPointDataView();*/
207     DataMaths::symmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
208 ksteube 775 }
209    
210     void
211     DataConstant::nonsymmetric(DataAbstract* ev)
212     {
213     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
214     if (temp_ev==0) {
215     throw DataException("Error - DataConstant::nonsymmetric: casting to DataConstant failed (propably a programming error).");
216     }
217 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
218     DataArrayView& evView=ev->getPointDataView();*/
219     DataMaths::nonsymmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
220 ksteube 775 }
221    
222     void
223 gross 800 DataConstant::trace(DataAbstract* ev, int axis_offset)
224 ksteube 775 {
225     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
226     if (temp_ev==0) {
227 gross 800 throw DataException("Error - DataConstant::trace: casting to DataConstant failed (propably a programming error).");
228 ksteube 775 }
229 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
230     DataArrayView& evView=ev->getPointDataView();*/
231     ValueType& evVec=temp_ev->getVector();
232     const ShapeType& evShape=temp_ev->getShape();
233     DataMaths::trace(m_data,getShape(),0,evVec,evShape,0,axis_offset);
234 ksteube 775 }
235    
236     void
237 gross 804 DataConstant::swapaxes(DataAbstract* ev, int axis0, int axis1)
238 gross 800 {
239     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
240     if (temp_ev==0) {
241 gross 804 throw DataException("Error - DataConstant::swapaxes: casting to DataConstant failed (propably a programming error).");
242 gross 800 }
243 jfenwick 1796 // DataArrayView& thisView=getPointDataView();
244     // DataArrayView& evView=ev->getPointDataView();
245     DataMaths::swapaxes(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,axis0,axis1);
246 gross 800 }
247    
248     void
249 ksteube 775 DataConstant::transpose(DataAbstract* ev, int axis_offset)
250     {
251     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
252     if (temp_ev==0) {
253     throw DataException("Error - DataConstant::transpose: casting to DataConstant failed (propably a programming error).");
254     }
255 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
256     DataArrayView& evView=ev->getPointDataView();*/
257     DataMaths::transpose(m_data, getShape(),0, temp_ev->getVector(),temp_ev->getShape(),0,axis_offset);
258 ksteube 775 }
259    
260     void
261 gross 580 DataConstant::eigenvalues(DataAbstract* ev)
262     {
263     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
264     if (temp_ev==0) {
265     throw DataException("Error - DataConstant::eigenvalues: casting to DataConstant failed (propably a programming error).");
266     }
267 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
268     DataArrayView& evView=ev->getPointDataView();*/
269     DataMaths::eigenvalues(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
270 gross 580 }
271     void
272     DataConstant::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
273     {
274     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
275     if (temp_ev==0) {
276     throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
277     }
278     DataConstant* temp_V=dynamic_cast<DataConstant*>(V);
279     if (temp_V==0) {
280     throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
281     }
282 jfenwick 1796 // DataArrayView thisView=getPointDataView();
283     // DataArrayView evView=ev->getPointDataView();
284     // DataArrayView VView=V->getPointDataView();
285 gross 580
286 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);
287 gross 580 }
288    
289 gross 950 void
290 gross 1118 DataConstant::setToZero()
291     {
292 jfenwick 1796 DataTypes::ValueType::size_type n=m_data.size();
293 gross 1118 for (int i=0; i<n ;++i) m_data[i]=0.;
294     }
295    
296     void
297 gross 950 DataConstant::dump(const std::string fileName) const
298     {
299 gross 1149 #ifdef USE_NETCDF
300 jfenwick 1796 const NcDim* ncdims[DataTypes::maxRank];
301 gross 950 NcVar* var;
302 jfenwick 1796 int rank = getRank();
303 gross 950 int type= getFunctionSpace().getTypeCode();
304     int ndims =0;
305 jfenwick 1796 long dims[DataTypes::maxRank];
306 gross 1141 const double* d_ptr=&(m_data[0]);
307 jfenwick 1796 DataTypes::ShapeType shape = getShape();
308 jfenwick 1872 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
309     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
310 ksteube 1827 #ifdef PASO_MPI
311     MPI_Status status;
312     #endif
313 matt 1319
314 ksteube 1827 #ifdef PASO_MPI
315     /* Serialize NetCDF I/O */
316     if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81802, MPI_COMM_WORLD, &status);
317     #endif
318    
319 gross 950 // netCDF error handler
320     NcError err(NcError::verbose_nonfatal);
321     // Create the file.
322 ksteube 1827 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
323     NcFile dataFile(newFileName, NcFile::Replace);
324 gross 950 // check if writing was successful
325 matt 1319 if (!dataFile.is_valid())
326 gross 950 throw DataException("Error - DataConstant:: opening of netCDF file for output failed.");
327 matt 1319 if (!dataFile.add_att("type_id",0) )
328 gross 950 throw DataException("Error - DataConstant:: appending data type to netCDF file failed.");
329 matt 1319 if (!dataFile.add_att("rank",rank) )
330 gross 950 throw DataException("Error - DataConstant:: appending rank attribute to netCDF file failed.");
331 matt 1319 if (!dataFile.add_att("function_space_type",type))
332 gross 950 throw DataException("Error - DataConstant:: appending function space attribute to netCDF file failed.");
333 gross 580
334 gross 950 if (rank == 0) {
335 matt 1319 if( ! (ncdims[0] = dataFile.add_dim("l", 1)) )
336 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
337 gross 950 dims[0]=1,
338     ndims=1;
339     } else {
340     ndims=rank;
341     dims[0]=shape[0];
342 matt 1319 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
343 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
344 gross 950 if ( rank >1 ) {
345     dims[1]=shape[1];
346 matt 1319 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
347 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 1 to netCDF file failed.");
348 gross 950 }
349     if ( rank >2 ) {
350     dims[2]=shape[2];
351 matt 1319 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
352 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 2 to netCDF file failed.");
353 gross 950 }
354     if ( rank >3 ) {
355     dims[3]=shape[3];
356 matt 1319 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
357 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 3 to netCDF file failed.");
358 gross 950 }
359     }
360 matt 1319
361 gross 950 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
362     throw DataException("Error - DataConstant:: appending variable to netCDF file failed.");
363 gross 1141 if (! (var->put(d_ptr,dims)) )
364 matt 1319 throw DataException("Error - DataConstant:: copy data to netCDF buffer failed.");
365 ksteube 1827 #ifdef PASO_MPI
366     if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81802, MPI_COMM_WORLD);
367     #endif
368 gross 1023 #else
369     throw DataException("Error - DataConstant:: dump is not configured with netCDF. Please contact your installation manager.");
370     #endif
371 gross 950 }
372    
373 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