/[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/esys2/escript/src/Data/DataConstant.cpp revision 82 by jgs, Tue Oct 26 06:53:54 2004 UTC trunk/escript/src/DataConstant.cpp revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1 
 /*  
  ******************************************************************************  
  *                                                                            *  
  *       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"  /*******************************************************
3  #include "escript/Data/DataException.h"  *
4    * Copyright (c) 2003-2009 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    
14    
15    #include "Data.h"
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    #ifdef PASO_MPI
26    #include <mpi.h>
27    #endif
28    
29    #include <boost/python/extract.hpp>
30    #include "DataMaths.h"
31    
32    // #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {throw DataException("Attempt to modify shared object");}
33    
34    #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {std::ostringstream ss; ss << " Attempt to modify shared object. line " << __LINE__ << " of " << __FILE__; ss << m_owners.size(); cerr << ss << endl; /* *((int*)0)=17; */throw DataException(ss.str());}
35    
36  using namespace std;  using namespace std;
37    using namespace boost::python;
38    
39  namespace escript {  namespace escript {
40    
41    DataConstant::DataConstant(const boost::python::numeric::array& value, const FunctionSpace& what): DataAbstract(what)  DataConstant::DataConstant(const WrappedArray& value,
42    {                             const FunctionSpace& what)
43      //cout << "Calling DataConstant constructor 1." << endl;    : parent(what,value.getShape())
44      DataArray temp(value);  {
45      //    m_data.copyFromArray(value,1);
46      // copy the data in the correct format  }
47      m_data=temp.getData();  
48      DataArrayView tempView(m_data,temp.getView().getShape());  DataConstant::DataConstant(const DataConstant& other)
49      //    : parent(other.getFunctionSpace(),other.getShape())
50      // copy the view of the data  {
51      setPointDataView(tempView);    m_data=other.m_data;
52    }  }
53    
54    DataConstant::DataConstant(const DataArrayView& value, const FunctionSpace& what): DataAbstract(what)  DataConstant::DataConstant(const DataConstant& other,
55    {                             const DataTypes::RegionType& region)
56      //cout << "Calling DataConstant constructor 2." << endl;    : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
57      //  {
58      // copy the data in the correct format    //
59      m_data=value.getData();    // allocate space for this new DataConstant's data
60      DataArrayView tempView(m_data,value.getShape());    int len = getNoValues();
61      //    m_data.resize(len,0.,len);
62      // copy the view of the data    //
63      setPointDataView(tempView);    // create a view of the data with the correct shape
64    }    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
65      //
66    DataConstant::DataConstant(const DataConstant& other): DataAbstract(other.getFunctionSpace()), m_data(other.m_data)    // load the view with the data from the slice
67    {    DataTypes::copySlice(m_data,getShape(),0,other.getVectorRO(),other.getShape(),0,region_loop_range);
68      //cout << "Calling DataConstant copy constructor." << endl;  }
69      //  
70      DataArrayView tempView(m_data,other.getPointDataView().getShape());  DataConstant::DataConstant(const FunctionSpace& what,
71      //                             const DataTypes::ShapeType &shape,
72      // copy the view of the data                             const DataTypes::ValueType &data)
73      setPointDataView(tempView);    : parent(what,shape)
74    }  {
75      //
76    DataConstant::DataConstant(const DataConstant& other, const DataArrayView::RegionType& region): DataAbstract(other.getFunctionSpace())    // copy the data in the correct format
77    {    m_data=data;
78      //cout << "Calling DataConstant slice constructor." << endl;    //
79      //    // create the view of the data
80      //  //   DataArrayView tempView(m_data,shape);
81      // get the shape of the slice  //   setPointDataView(tempView);
82      DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));  }
83      //  
84      // allocate space for this DataConstant  string
85      m_data.resize(DataArrayView::noValues(shape));  DataConstant::toString() const
86      DataArrayView tempView(m_data,shape);  {
87      tempView.copySlice(other.getPointDataView(),region);    return DataTypes::pointToString(m_data,getShape(),0,"");
88      //  }
89      // copy the view of the data  
90      setPointDataView(tempView);  
91    }  DataAbstract*
92    DataConstant::deepCopy()
93    string DataConstant::toString() const  {
94    {    return new DataConstant(*this);
95      return getPointDataView().toString("");  }
96    }  
97    
98    DataArrayView::ValueType::size_type DataConstant::getPointOffset(int sampleNo, int dataPointNo) const  DataTypes::ValueType::size_type
99    {  DataConstant::getPointOffset(int sampleNo,
100      EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),                               int dataPointNo) const
101             "Invalid index, sampleNo: " << sampleNo  {
102             << " dataPointNo: " << dataPointNo);    EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
103      return 0;                "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
104    }    //
105      // Whatever the coord's always return the same value as this is constant data.
106    DataArrayView::ValueType::size_type DataConstant::getLength() const    return 0;
107    {  }
108      return m_data.size();  
109    }  DataTypes::ValueType::size_type
110    DataConstant::getPointOffset(int sampleNo,
111    DataArrayView DataConstant::getDataPoint(int sampleNo, int dataPointNo)                               int dataPointNo)
112    {  {
113      EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),    EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
114             "Invalid index, sampleNo: " << sampleNo                "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
115             << " dataPointNo: " << dataPointNo);    //
116      //    // Whatever the coord's always return the same value as this is constant data.
117      // Whatever the coord's always return the same value    return 0;
118      return getPointDataView();  }
119    }  
120      DataTypes::ValueType::size_type
121    DataAbstract* DataConstant::getSlice(const DataArrayView::RegionType& region) const  DataConstant::getLength() const
122    {  {
123      return new DataConstant(*this,region);    return m_data.size();
124    }  }
125    
126    void DataConstant::setSlice(const DataAbstract* value, const DataArrayView::RegionType& region)  DataAbstract*
127    {  DataConstant::getSlice(const DataTypes::RegionType& region) const
128      const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value);  {
129      if (tempDataConst==0)    return new DataConstant(*this,region);
130      {  }
131        throw DataException("Programming error - casting to DataConstant.");  
132      }  void
133      getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region);  DataConstant::setSlice(const DataAbstract* value,
134    }                         const DataTypes::RegionType& region)
135    {
136    void DataConstant::reshapeDataPoint(const DataArrayView::ShapeType& shape)    const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value);
137    {    if (tempDataConst==0) {
138      //      throw DataException("Programming error - casting to DataConstant.");
139      // reshape a rank zero data point    }
140      if (getPointDataView().getRank()!=0) {    CHECK_FOR_EX_WRITE
141        stringstream temp;    //
142        temp << "Error - Can only reshape Data with data points of rank 0. "    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
143         << "This Data has data points with rank: "    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
144         << getPointDataView().getRank();    //
145        throw DataException(temp.str());    // check shape:
146      }    if (getRank()!=region.size()) {
147      m_data.resize(DataArrayView::noValues(shape),getPointDataView()());      throw DataException("Error - Invalid slice region.");
148      DataArrayView newView(m_data,shape);    }
149      setPointDataView(newView);    if (getRank()>0 && !DataTypes::checkShape(value->getShape(),shape)) {
150    }      throw DataException (DataTypes::createShapeErrorMessage(
151                    "Error - Couldn't copy slice due to shape mismatch.",shape,value->getShape()));
152      }
153      //   getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region_loop_range);
154      DataTypes::copySliceFrom(m_data,getShape(),0,tempDataConst->getVectorRO(), tempDataConst->getShape(),0,region_loop_range);
155    }
156    
157    
158    
159    void
160    DataConstant::symmetric(DataAbstract* ev)
161    {
162      DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
163      if (temp_ev==0) {
164        throw DataException("Error - DataConstant::symmetric: casting to DataConstant failed (propably a programming error).");
165      }
166      DataMaths::symmetric(m_data,getShape(),0,temp_ev->getVectorRW(), temp_ev->getShape(),0);
167    }
168    
169    void
170    DataConstant::nonsymmetric(DataAbstract* ev)
171    {
172      DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
173      if (temp_ev==0) {
174        throw DataException("Error - DataConstant::nonsymmetric: casting to DataConstant failed (propably a programming error).");
175      }
176      DataMaths::nonsymmetric(m_data,getShape(),0,temp_ev->getVectorRW(), temp_ev->getShape(),0);
177    }
178    
179    void
180    DataConstant::trace(DataAbstract* ev, int axis_offset)
181    {
182      DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
183      if (temp_ev==0) {
184        throw DataException("Error - DataConstant::trace: casting to DataConstant failed (propably a programming error).");
185      }
186      ValueType& evVec=temp_ev->getVectorRW();
187      const ShapeType& evShape=temp_ev->getShape();
188      DataMaths::trace(m_data,getShape(),0,evVec,evShape,0,axis_offset);
189    }
190    
191    void
192    DataConstant::swapaxes(DataAbstract* ev, int axis0, int axis1)
193    {
194      DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
195      if (temp_ev==0) {
196        throw DataException("Error - DataConstant::swapaxes: casting to DataConstant failed (propably a programming error).");
197      }
198      DataMaths::swapaxes(m_data,getShape(),0,temp_ev->getVectorRW(), temp_ev->getShape(),0,axis0,axis1);
199    }
200    
201    void
202    DataConstant::transpose(DataAbstract* ev, int axis_offset)
203    {
204      DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
205      if (temp_ev==0) {
206        throw DataException("Error - DataConstant::transpose: casting to DataConstant failed (propably a programming error).");
207      }
208      DataMaths::transpose(m_data, getShape(),0, temp_ev->getVectorRW(),temp_ev->getShape(),0,axis_offset);
209    }
210    
211    void
212    DataConstant::eigenvalues(DataAbstract* ev)
213    {
214      DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
215      if (temp_ev==0) {
216        throw DataException("Error - DataConstant::eigenvalues: casting to DataConstant failed (propably a programming error).");
217      }
218      DataMaths::eigenvalues(m_data,getShape(),0,temp_ev->getVectorRW(), temp_ev->getShape(),0);
219    }
220    void
221    DataConstant::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
222    {
223      DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
224      if (temp_ev==0) {
225        throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
226      }
227      DataConstant* temp_V=dynamic_cast<DataConstant*>(V);
228      if (temp_V==0) {
229        throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
230      }
231      DataMaths::eigenvalues_and_eigenvectors(m_data, getShape(),0,temp_ev->getVectorRW(), temp_ev->getShape(),0,temp_V->getVectorRW(), temp_V->getShape(),0,tol);
232    }
233    
234    void
235    DataConstant::setToZero()
236    {
237        CHECK_FOR_EX_WRITE
238        DataTypes::ValueType::size_type n=m_data.size();
239        for (int i=0; i<n ;++i) m_data[i]=0.;
240    }
241    
242    void
243    DataConstant::dump(const std::string fileName) const
244    {
245       #ifdef USE_NETCDF
246       const NcDim* ncdims[DataTypes::maxRank];
247       NcVar* var;
248       int rank = getRank();
249       int type=  getFunctionSpace().getTypeCode();
250       int ndims =0;
251       long dims[DataTypes::maxRank];
252       const double* d_ptr=&(m_data[0]);
253       DataTypes::ShapeType shape = getShape();
254       int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
255       int mpi_num=getFunctionSpace().getDomain()->getMPISize();
256    #ifdef PASO_MPI
257       MPI_Status status;
258    #endif
259    
260    #ifdef PASO_MPI
261       /* Serialize NetCDF I/O */
262       if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81802, MPI_COMM_WORLD, &status);
263    #endif
264    
265       // netCDF error handler
266       NcError err(NcError::verbose_nonfatal);
267       // Create the file.
268       char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
269       NcFile dataFile(newFileName, NcFile::Replace);
270       // check if writing was successful
271       if (!dataFile.is_valid())
272        throw DataException("Error - DataConstant:: opening of netCDF file for output failed.");
273       if (!dataFile.add_att("type_id",0) )
274        throw DataException("Error - DataConstant:: appending data type to netCDF file failed.");
275       if (!dataFile.add_att("rank",rank) )
276        throw DataException("Error - DataConstant:: appending rank attribute to netCDF file failed.");
277       if (!dataFile.add_att("function_space_type",type))
278        throw DataException("Error - DataConstant:: appending function space attribute to netCDF file failed.");
279    
280       if (rank == 0) {
281          if( ! (ncdims[0] = dataFile.add_dim("l", 1)) )
282            throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
283          dims[0]=1,
284          ndims=1;
285       } else {
286           ndims=rank;
287           dims[0]=shape[0];
288           if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
289            throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
290           if ( rank >1 ) {
291               dims[1]=shape[1];
292               if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
293            throw DataException("Error - DataConstant:: appending ncdimension 1 to netCDF file failed.");
294           }
295           if ( rank >2 ) {
296               dims[2]=shape[2];
297               if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
298            throw DataException("Error - DataConstant:: appending ncdimension 2 to netCDF file failed.");
299           }
300           if ( rank >3 ) {
301               dims[3]=shape[3];
302               if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
303            throw DataException("Error - DataConstant:: appending ncdimension 3 to netCDF file failed.");
304           }
305       }
306    
307       if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
308        throw DataException("Error - DataConstant:: appending variable to netCDF file failed.");
309       if (! (var->put(d_ptr,dims)) )
310             throw DataException("Error - DataConstant:: copy data to netCDF buffer failed.");
311    #ifdef PASO_MPI
312       if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81802, MPI_COMM_WORLD);
313    #endif
314       #else
315       throw DataException("Error - DataConstant:: dump is not configured with netCDF. Please contact your installation manager.");
316       #endif
317    }
318    
319    // These used to be marked as inline in DataConstant.
320    // But they are marked virtual in DataReady
321    DataTypes::ValueType&
322    DataConstant::getVectorRW()
323    {
324      CHECK_FOR_EX_WRITE
325      return m_data;
326    }
327    
328    const DataTypes::ValueType&
329    DataConstant::getVectorRO() const
330    {
331      return m_data;
332    }
333    
334  }  // end of namespace  }  // end of namespace
   
   
   
   

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

  ViewVC Help
Powered by ViewVC 1.1.26