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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1800 - (hide annotations)
Thu Sep 18 05:28:20 2008 UTC (11 years ago) by ksteube
File size: 12366 byte(s)
Serialized parallel I/O when writing mesh or data to NetCDF file on multiple MPI processors.
Added domain method getMPIComm() to complement getMPISize() and getMPIRank().

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