/[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

trunk/escript/src/Data/DataConstant.cpp revision 155 by jgs, Wed Nov 9 02:02:19 2005 UTC trunk/escript/src/DataConstant.cpp revision 1796 by jfenwick, Wed Sep 17 01:45:46 2008 UTC
# Line 1  Line 1 
 //$Id$  
 /*  
  ******************************************************************************  
  *                                                                            *  
  *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  
  *                                                                            *  
  * This software is the property of ACcESS. No part of this code              *  
  * may be copied in any form or by any means without the expressed written    *  
  * consent of ACcESS.  Copying, use or modification of this software          *  
  * by any unauthorised person is illegal unless that person has a software    *  
  * license agreement with ACcESS.                                             *  
  *                                                                            *  
  ******************************************************************************  
 */  
1    
2  #include "escript/Data/DataConstant.h"  /* $Id$ */
3  #include "escript/Data/DataException.h"  
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    #include "DataConstant.h"
17    #include "DataException.h"
18  #include "esysUtils/EsysAssert.h"  #include "esysUtils/EsysAssert.h"
19    
20  #include <iostream>  #include <iostream>
21  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
22    #ifdef USE_NETCDF
23    #include <netcdfcpp.h>
24    #endif
25    
26    #include <boost/python/extract.hpp>
27    #include "DataMaths.h"
28    
29  using namespace std;  using namespace std;
30    using namespace boost::python;
31    
32  namespace escript {  namespace escript {
33    
34  DataConstant::DataConstant(const boost::python::numeric::array& value,  DataConstant::DataConstant(const boost::python::numeric::array& value,
35                             const FunctionSpace& what)                             const FunctionSpace& what)
36    : DataAbstract(what)    : DataAbstract(what,DataTypes::shapeFromNumArray(value))
37  {  {
38    DataArray temp(value);  //   // extract the shape of the numarray
39    //  //   DataTypes::ShapeType tempShape;
40    // copy the data in the correct format  //   for (int i=0; i < value.getrank(); i++) {
41    m_data=temp.getData();  //     tempShape.push_back(extract<int>(value.getshape()[i]));
42    //  //   }
43    // create the view of the data  
44    DataArrayView tempView(m_data,temp.getView().getShape());    // get the space for the data vector
45    setPointDataView(tempView);  //   int len = getNoValues();
46  }  //   DataVector temp_data(len, 0.0, len);
47    //   DataArrayView temp_dataView(temp_data, tempShape);
48    //   temp_dataView.copy(value);
49    
50  DataConstant::DataConstant(const DataArrayView& value,    m_data.copyFromNumArray(value);
                            const FunctionSpace& what)  
   : DataAbstract(what)  
 {  
51    //    //
52    
53    // copy the data in the correct format    // copy the data in the correct format
54    m_data=value.getData();  //   m_data=temp_data;
55    //    //
56    // create the view of the data    // create the view of the data
57    DataArrayView tempView(m_data,value.getShape());  //   DataArrayView tempView(m_data,temp_dataView.getShape());
58    setPointDataView(tempView);  //   setPointDataView(tempView);
59  }  }
60    
61    // 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  DataConstant::DataConstant(const DataConstant& other)  DataConstant::DataConstant(const DataConstant& other)
75    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace(),other.getShape())
76  {  {  //
   //  
77    // copy the data in the correct format    // copy the data in the correct format
78    m_data=other.m_data;    m_data=other.m_data;
79    //    //
80    // create the view of the data  //   // create the view of the data
81    DataArrayView tempView(m_data,other.getPointDataView().getShape());  //   DataArrayView tempView(m_data,other.getPointDataView().getShape());
82    setPointDataView(tempView);  //   setPointDataView(tempView);
83  }  }
84    
85  DataConstant::DataConstant(const DataConstant& other,  DataConstant::DataConstant(const DataConstant& other,
86                             const DataArrayView::RegionType& region)                             const DataTypes::RegionType& region)
87    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
88  {  {
89    //    //
90    // get the shape of the slice to copy from    // get the shape of the slice to copy from
91    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));  //   DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
92    //    //
93    // allocate space for this new DataConstant's data    // allocate space for this new DataConstant's data
94    int len = DataArrayView::noValues(shape);    int len = getNoValues();
95    m_data.resize(len,0.,len);    m_data.resize(len,0.,len);
96    //    //
97    // create a view of the data with the correct shape    // create a view of the data with the correct shape
98    DataArrayView tempView(m_data,shape);  //   DataArrayView tempView(m_data,shape);
99    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
100    //    //
101    // load the view with the data from the slice    // load the view with the data from the slice
102    tempView.copySlice(other.getPointDataView(),region_loop_range);  //   tempView.copySlice(other.getPointDataView(),region_loop_range);
103    setPointDataView(tempView);    DataTypes::copySlice(m_data,getShape(),0,other.getVector(),other.getShape(),0,region_loop_range);
104    //   setPointDataView(tempView);
105  }  }
106    
107  DataConstant::DataConstant(const FunctionSpace& what,  DataConstant::DataConstant(const FunctionSpace& what,
108                             const DataArrayView::ShapeType &shape,                             const DataTypes::ShapeType &shape,
109                             const DataArrayView::ValueType &data)                             const DataTypes::ValueType &data)
110    : DataAbstract(what)    : DataAbstract(what,shape)
111  {  {
112    //    //
113    // copy the data in the correct format    // copy the data in the correct format
114    m_data=data;    m_data=data;
115    //    //
116    // create the view of the data    // create the view of the data
117    DataArrayView tempView(m_data,shape);  //   DataArrayView tempView(m_data,shape);
118    setPointDataView(tempView);  //   setPointDataView(tempView);
119  }  }
120    
121  string  string
122  DataConstant::toString() const  DataConstant::toString() const
123  {  {
124    return getPointDataView().toString("");    return DataTypes::pointToString(m_data,getShape(),0,"");
125  }  }
126    
127  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
128  DataConstant::getPointOffset(int sampleNo,  DataConstant::getPointOffset(int sampleNo,
129                               int dataPointNo) const                               int dataPointNo) const
130  {  {
# Line 115  DataConstant::getPointOffset(int sampleN Line 135  DataConstant::getPointOffset(int sampleN
135    return 0;    return 0;
136  }  }
137    
138  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
139  DataConstant::getLength() const  DataConstant::getLength() const
140  {  {
141    return m_data.size();    return m_data.size();
142  }  }
143    
144  DataArrayView  // DataArrayView
145  DataConstant::getDataPoint(int sampleNo,  // DataConstant::getDataPoint(int sampleNo,
146                             int dataPointNo)  //                            int dataPointNo)
147  {  // {
148    EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),  //   EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
149               "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);  //              "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
150    //  //   //
151    // 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.
152    return getPointDataView();  //   return getPointDataView();
153  }  // }
154      
155  DataAbstract*  DataAbstract*
156  DataConstant::getSlice(const DataArrayView::RegionType& region) const  DataConstant::getSlice(const DataTypes::RegionType& region) const
157  {  {
158    return new DataConstant(*this,region);    return new DataConstant(*this,region);
159  }  }
160    
161  void  void
162  DataConstant::setSlice(const DataAbstract* value,  DataConstant::setSlice(const DataAbstract* value,
163                         const DataArrayView::RegionType& region)                         const DataTypes::RegionType& region)
164  {  {
165    const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value);    const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value);
166    if (tempDataConst==0) {    if (tempDataConst==0) {
167      throw DataException("Programming error - casting to DataConstant.");      throw DataException("Programming error - casting to DataConstant.");
168    }    }
169    //    //
170    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
171    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
172    //    //
173    // check shape:    // check shape:
174    if (getPointDataView().getRank()!=region.size()) {    if (getRank()!=region.size()) {
175      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
176    }    }
177    if (tempDataConst->getPointDataView().getRank()>0 and !value->getPointDataView().checkShape(shape)) {    if (getRank()>0 && !DataTypes::checkShape(value->getShape(),shape)) {
178      throw DataException (value->getPointDataView().createShapeErrorMessage(      throw DataException (DataTypes::createShapeErrorMessage(
179                  "Error - Couldn't copy slice due to shape mismatch.",shape));                  "Error - Couldn't copy slice due to shape mismatch.",shape,value->getShape()));
180    }    }
181    //    //   getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region_loop_range);
182    getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region_loop_range);    DataTypes::copySliceFrom(m_data,getShape(),0,tempDataConst->getVector(), tempDataConst->getShape(),0,region_loop_range);
183    }
184    
185    
186    
187    void
188    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    /*  DataArrayView& thisView=getPointDataView();
195      DataArrayView& evView=ev->getPointDataView();*/
196      DataMaths::symmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
197    }
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    /*  DataArrayView& thisView=getPointDataView();
207      DataArrayView& evView=ev->getPointDataView();*/
208      DataMaths::nonsymmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
209  }  }
210    
211  void  void
212  DataConstant::reshapeDataPoint(const DataArrayView::ShapeType& shape)  DataConstant::trace(DataAbstract* ev, int axis_offset)
213  {  {
214    if (getPointDataView().getRank()!=0) {    DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
215      stringstream temp;    if (temp_ev==0) {
216      temp << "Error - Can only reshape Data with data points of rank 0. "      throw DataException("Error - DataConstant::trace: casting to DataConstant failed (propably a programming error).");
          << "This Data has data points with rank: " << getPointDataView().getRank();  
     throw DataException(temp.str());  
217    }    }
218    int len = DataArrayView::noValues(shape);  /*  DataArrayView& thisView=getPointDataView();
219    m_data.resize(len,getPointDataView()(),len);    DataArrayView& evView=ev->getPointDataView();*/
220    DataArrayView newView(m_data,shape);    ValueType& evVec=temp_ev->getVector();
221    setPointDataView(newView);    const ShapeType& evShape=temp_ev->getShape();
222      DataMaths::trace(m_data,getShape(),0,evVec,evShape,0,axis_offset);
223  }  }
224    
225  int  void
226  DataConstant::archiveData(ofstream& archiveFile,  DataConstant::swapaxes(DataAbstract* ev, int axis0, int axis1)
                           const DataArrayView::ValueType::size_type noValues) const  
227  {  {
228    return(m_data.archiveData(archiveFile, noValues));    DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
229      if (temp_ev==0) {
230        throw DataException("Error - DataConstant::swapaxes: casting to DataConstant failed (propably a programming error).");
231      }
232    //   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  }  }
236    
237  int  void
238  DataConstant::extractData(ifstream& archiveFile,  DataConstant::transpose(DataAbstract* ev, int axis_offset)
239                            const DataArrayView::ValueType::size_type noValues)  {
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    /*  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    }
248    
249    void
250    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    /*  DataArrayView& thisView=getPointDataView();
257      DataArrayView& evView=ev->getPointDataView();*/
258      DataMaths::eigenvalues(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
259    }
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    //   DataArrayView thisView=getPointDataView();
272    //   DataArrayView evView=ev->getPointDataView();
273    //   DataArrayView VView=V->getPointDataView();
274    
275      DataMaths::eigenvalues_and_eigenvectors(m_data, getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,temp_V->getVector(), temp_V->getShape(),0,tol);
276    }
277    
278    void
279    DataConstant::setToZero()
280    {
281        DataTypes::ValueType::size_type n=m_data.size();
282        for (int i=0; i<n ;++i) m_data[i]=0.;
283    }
284    
285    void
286    DataConstant::dump(const std::string fileName) const
287  {  {
288    return(m_data.extractData(archiveFile, noValues));     #ifdef PASO_MPI
289       throw DataException("Error - DataConstant:: dump is not implemented for MPI yet.");
290       #endif
291       #ifdef USE_NETCDF
292       const NcDim* ncdims[DataTypes::maxRank];
293       NcVar* var;
294       int rank = getRank();
295       int type=  getFunctionSpace().getTypeCode();
296       int ndims =0;
297       long dims[DataTypes::maxRank];
298       const double* d_ptr=&(m_data[0]);
299       DataTypes::ShapeType shape = getShape();
300    
301       // 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       if (!dataFile.is_valid())
307        throw DataException("Error - DataConstant:: opening of netCDF file for output failed.");
308       if (!dataFile.add_att("type_id",0) )
309        throw DataException("Error - DataConstant:: appending data type to netCDF file failed.");
310       if (!dataFile.add_att("rank",rank) )
311        throw DataException("Error - DataConstant:: appending rank attribute to netCDF file failed.");
312       if (!dataFile.add_att("function_space_type",type))
313        throw DataException("Error - DataConstant:: appending function space attribute to netCDF file failed.");
314    
315       if (rank == 0) {
316          if( ! (ncdims[0] = dataFile.add_dim("l", 1)) )
317            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           if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
324            throw DataException("Error - DataConstant:: appending ncdimsion 0 to netCDF file failed.");
325           if ( rank >1 ) {
326               dims[1]=shape[1];
327               if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
328            throw DataException("Error - DataConstant:: appending ncdimsion 1 to netCDF file failed.");
329           }
330           if ( rank >2 ) {
331               dims[2]=shape[2];
332               if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
333            throw DataException("Error - DataConstant:: appending ncdimsion 2 to netCDF file failed.");
334           }
335           if ( rank >3 ) {
336               dims[3]=shape[3];
337               if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
338            throw DataException("Error - DataConstant:: appending ncdimsion 3 to netCDF file failed.");
339           }
340       }
341    
342       if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
343        throw DataException("Error - DataConstant:: appending variable to netCDF file failed.");
344       if (! (var->put(d_ptr,dims)) )
345             throw DataException("Error - DataConstant:: copy data to netCDF buffer failed.");
346       #else
347       throw DataException("Error - DataConstant:: dump is not configured with netCDF. Please contact your installation manager.");
348       #endif
349  }  }
350    
351  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.155  
changed lines
  Added in v.1796

  ViewVC Help
Powered by ViewVC 1.1.26