/[escript]/branches/lapack2681/escript/src/DataConstant.cpp
ViewVC logotype

Diff of /branches/lapack2681/escript/src/DataConstant.cpp

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

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

Legend:
Removed from v.82  
changed lines
  Added in v.1141

  ViewVC Help
Powered by ViewVC 1.1.26