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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.682  
changed lines
  Added in v.1811

  ViewVC Help
Powered by ViewVC 1.1.26