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

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

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

revision 1358 by gross, Wed Dec 5 03:41:06 2007 UTC revision 1800 by ksteube, Thu Sep 18 05:28:20 2008 UTC
# Line 13  Line 13 
13   *   *
14   *******************************************************/   *******************************************************/
15    
16    #include "Data.h"
17  #include "DataExpanded.h"  #include "DataExpanded.h"
18  #include "DataException.h"  #include "DataException.h"
19  #include "DataConstant.h"  #include "DataConstant.h"
# Line 20  Line 21 
21  #ifdef USE_NETCDF  #ifdef USE_NETCDF
22  #include <netcdfcpp.h>  #include <netcdfcpp.h>
23  #endif  #endif
24    #ifdef PASO_MPI
25    #include <mpi.h>
26    #endif
27    
28  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
29    #include "DataMaths.h"
30    
31  using namespace std;  using namespace std;
32  using namespace boost::python;  using namespace boost::python;
33  using namespace boost;  using namespace boost;
34    using namespace escript::DataTypes;
35    
36  namespace escript {  namespace escript {
37    
38  DataExpanded::DataExpanded(const boost::python::numeric::array& value,  DataExpanded::DataExpanded(const boost::python::numeric::array& value,
39                             const FunctionSpace& what)                             const FunctionSpace& what)
40    : DataAbstract(what)    : DataAbstract(what,DataTypes::shapeFromNumArray(value))
41  {  {
42    DataArrayView::ShapeType tempShape;  /*  DataTypes::ShapeType tempShape;
43    //    //
44    // extract the shape of the python numarray    // extract the shape of the python numarray
45    for (int i=0; i<value.getrank(); i++) {    for (int i=0; i<value.getrank(); i++) {
46      tempShape.push_back(extract<int>(value.getshape()[i]));      tempShape.push_back(extract<int>(value.getshape()[i]));
47    }    }*/
48    //    //
49    // initialise the data array for this object    // initialise the data array for this object
50    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    initialise(what.getNumSamples(),what.getNumDPPSample());
51    //    //
52    // copy the given value to every data point    // copy the given value to every data point
53    copy(value);    copy(value);
54  }  }
55    
56  DataExpanded::DataExpanded(const DataExpanded& other)  DataExpanded::DataExpanded(const DataExpanded& other)
57    : DataAbstract(other.getFunctionSpace()),    : DataAbstract(other.getFunctionSpace(), other.getShape()),
58    m_data(other.m_data)    m_data(other.m_data)
59  {  {
60    /*
61    //    //
62    // create the view for the data    // create the view for the data
63    DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());    DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());
64    setPointDataView(temp);    setPointDataView(temp);
65    
66      // keep local shape in sync
67      setShape(other.getPointDataView().getShape());*/
68  }  }
69    
70  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
71    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace(), other.getShape())
72  {  {
73    //    //
74    // initialise the data array for this object    // initialise the data array for this object
75    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
76    //    //
77    // DataConstant only has one value, copy this to every data point    // DataConstant only has one value, copy this to every data point
78    copy(other.getPointDataView());    copy(other);
79  }  }
80    
81  DataExpanded::DataExpanded(const DataTagged& other)  DataExpanded::DataExpanded(const DataTagged& other)
82    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace(), other.getShape())
83  {  {
84    //    //
85    // initialise the data array for this object    // initialise the data array for this object
86    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
87    //    //
88    // for each data point in this object, extract and copy the corresponding data    // for each data point in this object, extract and copy the corresponding data
89    // value from the given DataTagged object    // value from the given DataTagged object
90    int i,j;    int i,j;
91    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
92    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
93    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
94    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
95      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
96        try {        try {
97          getPointDataView().copy(getPointOffset(i,j),  //         getPointDataView().copy(getPointOffset(i,j),
98                                  other.getPointDataView(),  //                                 other.getPointDataView(),
99    //                                 other.getPointOffset(i,j));
100               DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(),
101                                    other.getVector(),
102                                  other.getPointOffset(i,j));                                  other.getPointOffset(i,j));
103        }        }
104        catch (std::exception& e) {        catch (std::exception& e) {
# Line 96  DataExpanded::DataExpanded(const DataTag Line 109  DataExpanded::DataExpanded(const DataTag
109  }  }
110    
111  DataExpanded::DataExpanded(const DataExpanded& other,  DataExpanded::DataExpanded(const DataExpanded& other,
112                             const DataArrayView::RegionType& region)                             const DataTypes::RegionType& region)
113    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
114  {  {
115    //    //
116    // get the shape of the slice    // get the shape of the slice
117    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));  //   DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
118    //    //
119    // initialise this Data object to the shape of the slice    // initialise this Data object to the shape of the slice
120    initialise(shape,other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
121    //    //
122    // copy the data    // copy the data
123    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
124    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
125    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
126    int i,j;    int i,j;
127    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
128    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
129      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
130        try {        try {
131          getPointDataView().copySlice(getPointOffset(i,j),  //         getPointDataView().copySlice(getPointOffset(i,j),
132                                       other.getPointDataView(),  //                                      other.getPointDataView(),
133    //                                      other.getPointOffset(i,j),
134    //                                      region_loop_range);
135            DataTypes::copySlice(getVector(),getShape(),getPointOffset(i,j),
136                                         other.getVector(),
137                         other.getShape(),
138                                       other.getPointOffset(i,j),                                       other.getPointOffset(i,j),
139                                       region_loop_range);                                       region_loop_range);
140        }        }
# Line 127  DataExpanded::DataExpanded(const DataExp Line 145  DataExpanded::DataExpanded(const DataExp
145    }    }
146  }  }
147    
148  DataExpanded::DataExpanded(const DataArrayView& value,  // DataExpanded::DataExpanded(const DataArrayView& value,
149                             const FunctionSpace& what)  //                            const FunctionSpace& what)
150    : DataAbstract(what)  //   : DataAbstract(what)
151  {  // {
152    //  //   //
153    // get the shape of the given data value  //   // get the shape of the given data value
154    DataArrayView::ShapeType tempShape=value.getShape();  //   DataTypes::ShapeType tempShape=value.getShape();
155    //  //   //
156    // initialise this Data object to the shape of the given data value  //   // initialise this Data object to the shape of the given data value
157    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());  //   initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
158    //  //   //
159    // copy the given value to every data point  //   // copy the given value to every data point
160    copy(value);  //   copy(value);
161  }  // }
162    
163  DataExpanded::DataExpanded(const FunctionSpace& what,  DataExpanded::DataExpanded(const FunctionSpace& what,
164                             const DataArrayView::ShapeType &shape,                             const DataTypes::ShapeType &shape,
165                             const DataArrayView::ValueType &data)                             const DataTypes::ValueType &data)
166    : DataAbstract(what)    : DataAbstract(what,shape)
167  {  {
168    //    EsysAssert(data.size()%getNoValues()==0,
169    // create the view of the data                   "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
170    initialise(shape,what.getNumSamples(),what.getNumDPPSample());  
171    //    if (data.size()==getNoValues())
172    // copy the data in the correct format    {
173    m_data.getData()=data;       ValueType& vec=m_data.getData();
174         //
175         // create the view of the data
176         initialise(what.getNumSamples(),what.getNumDPPSample());
177         // now we copy this value to all elements
178         for (int i=0;i<getLength();)
179         {
180        for (int j=0;j<getNoValues();++j,++i)
181        {
182            vec[i]=data[j];
183        }
184         }
185      }
186      else
187      {
188         //
189         // copy the data in the correct format
190         m_data.getData()=data;
191      }
192    
193    
194  }  }
195    
196  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
# Line 160  DataExpanded::~DataExpanded() Line 198  DataExpanded::~DataExpanded()
198  }  }
199    
200  DataAbstract*  DataAbstract*
201  DataExpanded::getSlice(const DataArrayView::RegionType& region) const  DataExpanded::deepCopy()
202    {
203      return new DataExpanded(*this);
204    }
205    
206    
207    DataAbstract*
208    DataExpanded::getSlice(const DataTypes::RegionType& region) const
209  {  {
210    return new DataExpanded(*this,region);    return new DataExpanded(*this,region);
211  }  }
212    
213  void  void
214  DataExpanded::setSlice(const DataAbstract* value,  DataExpanded::setSlice(const DataAbstract* value,
215                         const DataArrayView::RegionType& region)                         const DataTypes::RegionType& region)
216  {  {
217    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
218    if (tempDataExp==0) {    if (tempDataExp==0) {
# Line 175  DataExpanded::setSlice(const DataAbstrac Line 220  DataExpanded::setSlice(const DataAbstrac
220    }    }
221    //    //
222    // get shape of slice    // get shape of slice
223    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
224    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
225    //    //
226    // check shape    // check shape
227    if (getPointDataView().getRank()!=region.size()) {    if (getRank()!=region.size()) {
228      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
229    }    }
230    if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {    if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
231      throw DataException (value->getPointDataView().createShapeErrorMessage(      throw DataException (DataTypes::createShapeErrorMessage(
232          "Error - Couldn't copy slice due to shape mismatch.",shape));          "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
233    }    }
234    //    //
235    // copy the data from the slice into this object    // copy the data from the slice into this object
236    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
237    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
238    int i, j;    int i, j;
239      ValueType& vec=getVector();
240      const ShapeType& mshape=getShape();
241      const ValueType& tVec=tempDataExp->getVector();
242      const ShapeType& tShape=tempDataExp->getShape();
243    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
244    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
245      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
246        getPointDataView().copySliceFrom(getPointOffset(i,j),  /*      getPointDataView().copySliceFrom(getPointOffset(i,j),
247                                         tempDataExp->getPointDataView(),                                         tempDataExp->getPointDataView(),
248                                         tempDataExp->getPointOffset(i,j),                                         tempDataExp->getPointOffset(i,j),
249                                           region_loop_range);*/
250            DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
251                                           tVec,
252                           tShape,
253                                           tempDataExp->getPointOffset(i,j),
254                                         region_loop_range);                                         region_loop_range);
255    
256      }      }
257    }    }
258  }  }
259    
260  void  void
261  DataExpanded::copy(const DataArrayView& value)  DataExpanded::copy(const DataConstant& value)
262  {  {
263      EsysAssert((checkShape(getShape(), value.getShape())),
264                     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
265    
266    //    //
267    // copy a single value to every data point in this object    // copy a single value to every data point in this object
268    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
# Line 217  DataExpanded::copy(const DataArrayView& Line 275  DataExpanded::copy(const DataArrayView&
275        // DOASSERT is on which of course will play        // DOASSERT is on which of course will play
276        // havoc with the omp threads. Run single threaded        // havoc with the omp threads. Run single threaded
277        // if using DOASSERT.        // if using DOASSERT.
278        getPointDataView().copy(getPointOffset(i,j),value);        //getPointDataView().copy(getPointOffset(i,j),value);
279          DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(), value.getVector(), 0);
280      }      }
281    }    }
282  }  }
283    
284    
285    // void
286    // DataExpanded::copy(const DataArrayView& value)
287    // {
288    //   //
289    //   // copy a single value to every data point in this object
290    //   int nRows=m_data.getNumRows();
291    //   int nCols=m_data.getNumCols();
292    //   int i,j;
293    //   #pragma omp parallel for private(i,j) schedule(static)
294    //   for (i=0;i<nRows;i++) {
295    //     for (j=0;j<nCols;j++) {
296    //       // NOTE: An exception may be thown from this call if
297    //       // DOASSERT is on which of course will play
298    //       // havoc with the omp threads. Run single threaded
299    //       // if using DOASSERT.
300    //       getPointDataView().copy(getPointOffset(i,j),value);
301    //     }
302    //   }
303    // }
304    
305  void  void
306  DataExpanded::copy(const boost::python::numeric::array& value)  DataExpanded::copy(const boost::python::numeric::array& value)
307  {  {
308    
309    // extract the shape of the numarray    // extract the shape of the numarray
310    DataArrayView::ShapeType tempShape;    DataTypes::ShapeType tempShape;
311    for (int i=0; i < value.getrank(); i++) {    for (int i=0; i < value.getrank(); i++) {
312      tempShape.push_back(extract<int>(value.getshape()[i]));      tempShape.push_back(extract<int>(value.getshape()[i]));
313    }    }
314    
315    // get the space for the data vector    // get the space for the data vector
316    int len = DataArrayView::noValues(tempShape);  //   int len = DataTypes::noValues(tempShape);
317    DataVector temp_data(len, 0.0, len);  //   DataVector temp_data(len, 0.0, len);
318    DataArrayView temp_dataView(temp_data, tempShape);  //   DataArrayView temp_dataView(temp_data, tempShape);
319    temp_dataView.copy(value);  //   temp_dataView.copy(value);
320    
321    //    //
322    // check the input shape matches this shape    // check the input shape matches this shape
323    if (!getPointDataView().checkShape(temp_dataView.getShape())) {    if (!DataTypes::checkShape(getShape(),tempShape)) {
324      throw DataException(getPointDataView().createShapeErrorMessage(      throw DataException(DataTypes::createShapeErrorMessage(
325                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
326                          temp_dataView.getShape()));                          tempShape,getShape()));
327    }    }
328    //    //
329    // now copy over the data    // now copy over the data
330    copy(temp_dataView);    //copy(temp_dataView);
331      getVector().copyFromNumArray(value);
332  }  }
333    
334    
335  void  void
336  DataExpanded::initialise(const DataArrayView::ShapeType& shape,  DataExpanded::initialise(int noSamples,
                          int noSamples,  
337                           int noDataPointsPerSample)                           int noDataPointsPerSample)
338  {  {
339    //    //
340    // resize data array to the required size    // resize data array to the required size
341    m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));    m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
342    
343    //    //
344    // create the data view of the data array  //   // create the data view of the data array
345    DataArrayView temp(m_data.getData(),shape);  //   DataArrayView temp(m_data.getData(),shape);
346    setPointDataView(temp);  //   setPointDataView(temp);
347    
348    //   // keep shape in sync
349    //   setShape(shape);
350  }  }
351    
352  string  string
# Line 271  DataExpanded::toString() const Line 356  DataExpanded::toString() const
356    FunctionSpace fs=getFunctionSpace();    FunctionSpace fs=getFunctionSpace();
357    //    //
358    // create a temporary view as the offset will be changed    // create a temporary view as the offset will be changed
359    DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());  //  DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
360    
361      int offset=0;
362    for (int i=0;i<m_data.getNumRows();i++) {    for (int i=0;i<m_data.getNumRows();i++) {
363      for (int j=0;j<m_data.getNumCols();j++) {      for (int j=0;j<m_data.getNumCols();j++) {
364        tempView.setOffset(getPointOffset(i,j));        offset=getPointOffset(i,j);
365        stringstream suffix;        stringstream suffix;
366        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
367        temp << tempView.toString(suffix.str());        temp << DataTypes::pointToString(getVector(),getShape(),offset,suffix.str());
368        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
369          temp << endl;          temp << endl;
370        }        }
# Line 286  DataExpanded::toString() const Line 373  DataExpanded::toString() const
373    return temp.str();    return temp.str();
374  }  }
375    
376  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
377  DataExpanded::getPointOffset(int sampleNo,  DataExpanded::getPointOffset(int sampleNo,
378                               int dataPointNo) const                               int dataPointNo) const
379  {  {
380    return m_data.index(sampleNo,dataPointNo);    return m_data.index(sampleNo,dataPointNo);
381  }  }
382    
383  DataArrayView  // DataArrayView
384  DataExpanded::getDataPoint(int sampleNo,  // DataExpanded::getDataPoint(int sampleNo,
385                             int dataPointNo)  //                            int dataPointNo)
386  {  // {
387    DataArrayView temp(m_data.getData(),getPointDataView().getShape(),getPointOffset(sampleNo,dataPointNo));  //   DataArrayView temp(m_data.getData(),getPointDataView().getShape(),getPointOffset(sampleNo,dataPointNo));
388    return temp;  //   return temp;
389  }  // }
390    
391  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
392  DataExpanded::getLength() const  DataExpanded::getLength() const
393  {  {
394    return m_data.size();    return m_data.size();
395  }  }
396    
 int  
 DataExpanded::archiveData(ofstream& archiveFile,  
                           const DataArrayView::ValueType::size_type noValues) const  
 {  
   return(m_data.archiveData(archiveFile, noValues));  
 }  
397    
 int  
 DataExpanded::extractData(ifstream& archiveFile,  
                           const DataArrayView::ValueType::size_type noValues)  
 {  
   return(m_data.extractData(archiveFile, noValues));  
 }  
398    
399  void  void
400  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
# Line 327  DataExpanded::copyToDataPoint(const int Line 402  DataExpanded::copyToDataPoint(const int
402    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
403    int numSamples = getNumSamples();    int numSamples = getNumSamples();
404    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
405    int dataPointRank = getPointDataView().getRank();    int dataPointRank = getRank();
406    ShapeType dataPointShape = getPointDataView().getShape();    ShapeType dataPointShape = getShape();
407    if (numSamples*numDataPointsPerSample > 0) {    if (numSamples*numDataPointsPerSample > 0) {
408       //TODO: global error handling       //TODO: global error handling
409       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
# Line 337  DataExpanded::copyToDataPoint(const int Line 412  DataExpanded::copyToDataPoint(const int
412       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
413             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
414       }       }
415       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
416         ValueType& vec=getVector();
417       if (dataPointRank==0) {       if (dataPointRank==0) {
418           dataPointView()=value;           vec[0]=value;
419       } else if (dataPointRank==1) {       } else if (dataPointRank==1) {
420          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
421              dataPointView(i)=value;              vec[offset+i]=value;
422          }          }
423       } else if (dataPointRank==2) {       } else if (dataPointRank==2) {
424          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
425             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
426                dataPointView(i,j)=value;                vec[offset+getRelIndex(dataPointShape,i,j)]=value;
427             }             }
428          }          }
429       } else if (dataPointRank==3) {       } else if (dataPointRank==3) {
430          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
431             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
432                for (int k=0; k<dataPointShape[2]; k++) {                for (int k=0; k<dataPointShape[2]; k++) {
433                   dataPointView(i,j,k)=value;                   vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
434                }                }
435             }             }
436          }          }
# Line 363  DataExpanded::copyToDataPoint(const int Line 439  DataExpanded::copyToDataPoint(const int
439             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
440               for (int k=0; k<dataPointShape[2]; k++) {               for (int k=0; k<dataPointShape[2]; k++) {
441                 for (int l=0; l<dataPointShape[3]; l++) {                 for (int l=0; l<dataPointShape[3]; l++) {
442                    dataPointView(i,j,k,l)=value;                    vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
443                 }                 }
444               }               }
445             }             }
# Line 377  DataExpanded::copyToDataPoint(const int Line 453  DataExpanded::copyToDataPoint(const int
453    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
454    int numSamples = getNumSamples();    int numSamples = getNumSamples();
455    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
456    int dataPointRank = getPointDataView().getRank();    int dataPointRank = getRank();
457    ShapeType dataPointShape = getPointDataView().getShape();    const ShapeType& shape = getShape();
458    //    //
459    // check rank:    // check rank:
460    if (value.getrank()!=dataPointRank)    if (value.getrank()!=dataPointRank)
# Line 391  DataExpanded::copyToDataPoint(const int Line 467  DataExpanded::copyToDataPoint(const int
467       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
468             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
469       }       }
470       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);   //    DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
471         ValueType& vec=getVector();
472       if (dataPointRank==0) {       if (dataPointRank==0) {
473           dataPointView()=extract<double>(value[0]);           vec[0]=extract<double>(value[0]);
474       } else if (dataPointRank==1) {       } else if (dataPointRank==1) {
475          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<shape[0]; i++) {
476              dataPointView(i)=extract<double>(value[i]);              vec[i]=extract<double>(value[i]);
477          }          }
478       } else if (dataPointRank==2) {       } else if (dataPointRank==2) {
479          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<shape[0]; i++) {
480             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<shape[1]; j++) {
481                dataPointView(i,j)=extract<double>(value[i][j]);                vec[getRelIndex(shape,i,j)]=extract<double>(value[i][j]);
482             }             }
483          }          }
484       } else if (dataPointRank==3) {       } else if (dataPointRank==3) {
485          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<shape[0]; i++) {
486             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<shape[1]; j++) {
487                for (int k=0; k<dataPointShape[2]; k++) {                for (int k=0; k<shape[2]; k++) {
488                   dataPointView(i,j,k)=extract<double>(value[i][j][k]);                   vec[getRelIndex(shape,i,j,k)]=extract<double>(value[i][j][k]);
489                }                }
490             }             }
491          }          }
492       } else if (dataPointRank==4) {       } else if (dataPointRank==4) {
493           for (int i=0; i<dataPointShape[0]; i++) {           for (int i=0; i<shape[0]; i++) {
494             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<shape[1]; j++) {
495               for (int k=0; k<dataPointShape[2]; k++) {               for (int k=0; k<shape[2]; k++) {
496                 for (int l=0; l<dataPointShape[3]; l++) {                 for (int l=0; l<shape[3]; l++) {
497                    dataPointView(i,j,k,l)=extract<double>(value[i][j][k][l]);                    vec[getRelIndex(shape,i,j,k,l)]=extract<double>(value[i][j][k][l]);
498                 }                 }
499               }               }
500             }             }
# Line 431  DataExpanded::copyAll(const boost::pytho Line 508  DataExpanded::copyAll(const boost::pytho
508    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
509    int numSamples = getNumSamples();    int numSamples = getNumSamples();
510    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
511    int dataPointRank = getPointDataView().getRank();    int dataPointRank = getRank();
512    ShapeType dataPointShape = getPointDataView().getShape();    const ShapeType& dataPointShape = getShape();
513    //    //
514    // check rank:    // check rank:
515    if (value.getrank()!=dataPointRank+1)    if (value.getrank()!=dataPointRank+1)
# Line 440  DataExpanded::copyAll(const boost::pytho Line 517  DataExpanded::copyAll(const boost::pytho
517    if (value.getshape()[0]!=numSamples*numDataPointsPerSample)    if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
518         throw DataException("leading dimension of numarray is too small");         throw DataException("leading dimension of numarray is too small");
519    //    //
520      ValueType& vec=getVector();
521    int dataPoint = 0;    int dataPoint = 0;
522    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
523      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
524        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);  //       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
525          ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
526        if (dataPointRank==0) {        if (dataPointRank==0) {
527           dataPointView()=extract<double>(value[dataPoint]);           vec[offset]=extract<double>(value[dataPoint]);
528        } else if (dataPointRank==1) {        } else if (dataPointRank==1) {
529           for (int i=0; i<dataPointShape[0]; i++) {           for (int i=0; i<dataPointShape[0]; i++) {
530              dataPointView(i)=extract<double>(value[dataPoint][i]);              vec[offset+i]=extract<double>(value[dataPoint][i]);
531           }           }
532        } else if (dataPointRank==2) {        } else if (dataPointRank==2) {
533           for (int i=0; i<dataPointShape[0]; i++) {           for (int i=0; i<dataPointShape[0]; i++) {
534             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
535               dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);           vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
536    //             dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
537             }             }
538           }           }
539         } else if (dataPointRank==3) {         } else if (dataPointRank==3) {
540           for (int i=0; i<dataPointShape[0]; i++) {           for (int i=0; i<dataPointShape[0]; i++) {
541             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
542               for (int k=0; k<dataPointShape[2]; k++) {               for (int k=0; k<dataPointShape[2]; k++) {
543                   dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);           vec[offset+getRelIndex(dataPointShape,i,j,k)]=extract<double>(value[dataPoint][i][j][k]);
544               }               }
545             }             }
546           }           }
# Line 469  DataExpanded::copyAll(const boost::pytho Line 549  DataExpanded::copyAll(const boost::pytho
549             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
550               for (int k=0; k<dataPointShape[2]; k++) {               for (int k=0; k<dataPointShape[2]; k++) {
551                 for (int l=0; l<dataPointShape[3]; l++) {                 for (int l=0; l<dataPointShape[3]; l++) {
552                   dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);                   vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=extract<double>(value[dataPoint][i][j][k][l]);
553                 }                 }
554               }               }
555             }             }
# Line 489  DataExpanded::symmetric(DataAbstract* ev Line 569  DataExpanded::symmetric(DataAbstract* ev
569    if (temp_ev==0) {    if (temp_ev==0) {
570      throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
571    }    }
572    DataArrayView& thisView=getPointDataView();  //  DataArrayView& thisView=getPointDataView();
573    DataArrayView& evView=ev->getPointDataView();  //  DataArrayView& evView=ev->getPointDataView();
574      ValueType& vec=getVector();
575      const ShapeType& shape=getShape();
576      ValueType& evVec=temp_ev->getVector();
577      const ShapeType& evShape=temp_ev->getShape();
578    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
579    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
580      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
581           DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
582                                      evView,ev->getPointOffset(sampleNo,dataPointNo));                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
583      }      }
584    }    }
585  }  }
# Line 509  DataExpanded::nonsymmetric(DataAbstract* Line 593  DataExpanded::nonsymmetric(DataAbstract*
593    if (temp_ev==0) {    if (temp_ev==0) {
594      throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
595    }    }
596    DataArrayView& thisView=getPointDataView();  //   DataArrayView& thisView=getPointDataView();
597    DataArrayView& evView=ev->getPointDataView();  //   DataArrayView& evView=ev->getPointDataView();
598      ValueType& vec=getVector();
599      const ShapeType& shape=getShape();
600      ValueType& evVec=temp_ev->getVector();
601      const ShapeType& evShape=temp_ev->getShape();
602    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
603    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
604      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
605           DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),  //          DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
606                                      evView,ev->getPointOffset(sampleNo,dataPointNo));  //                                     evView,ev->getPointOffset(sampleNo,dataPointNo));
607             DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
608                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
609      }      }
610    }    }
611  }  }
# Line 529  DataExpanded::trace(DataAbstract* ev, in Line 619  DataExpanded::trace(DataAbstract* ev, in
619    if (temp_ev==0) {    if (temp_ev==0) {
620      throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
621    }    }
622    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
623    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
624      ValueType& evVec=temp_ev->getVector();
625      const ShapeType& evShape=temp_ev->getShape();
626    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
627    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
628      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
629           DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),  /*         DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),
630                                      evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);                                      evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);*/
631             DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
632                                        evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
633      }      }
634    }    }
635  }  }
# Line 550  DataExpanded::transpose(DataAbstract* ev Line 644  DataExpanded::transpose(DataAbstract* ev
644    if (temp_ev==0) {    if (temp_ev==0) {
645      throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
646    }    }
647    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
648    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
649      ValueType& evVec=temp_ev->getVector();
650      const ShapeType& evShape=temp_ev->getShape();
651    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
652    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
653      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
654           DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
655                                      evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
656      }      }
657    }    }
658  }  }
# Line 571  DataExpanded::swapaxes(DataAbstract* ev, Line 667  DataExpanded::swapaxes(DataAbstract* ev,
667    if (temp_ev==0) {    if (temp_ev==0) {
668      throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
669    }    }
670    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
671    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
672      ValueType& evVec=temp_ev->getVector();
673      const ShapeType& evShape=temp_ev->getShape();
674    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
675    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
676      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
677           DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
678                                      evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
679      }      }
680    }    }
681  }  }
# Line 591  DataExpanded::eigenvalues(DataAbstract* Line 689  DataExpanded::eigenvalues(DataAbstract*
689    if (temp_ev==0) {    if (temp_ev==0) {
690      throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
691    }    }
692    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
693    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
694      ValueType& evVec=temp_ev->getVector();
695      const ShapeType& evShape=temp_ev->getShape();
696    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
697    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
698      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
699           DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
700                                      evView,ev->getPointOffset(sampleNo,dataPointNo));                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
701      }      }
702    }    }
703  }  }
# Line 615  DataExpanded::eigenvalues_and_eigenvecto Line 715  DataExpanded::eigenvalues_and_eigenvecto
715    if (temp_V==0) {    if (temp_V==0) {
716      throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");      throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
717    }    }
718    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
719    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
720    DataArrayView& VView=V->getPointDataView();    ValueType& evVec=temp_ev->getVector();
721      const ShapeType& evShape=temp_ev->getShape();
722      ValueType& VVec=temp_V->getVector();
723      const ShapeType& VShape=temp_V->getShape();
724    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
725    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
726      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
727           DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
728                                                       evView,ev->getPointOffset(sampleNo,dataPointNo),                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
729                                                       VView,V->getPointOffset(sampleNo,dataPointNo),                                      VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
                                                      tol);  
730      }      }
731    }    }
732  }  }
733    
734  void  void
735  DataExpanded::setToZero(){  DataExpanded::setToZero(){
736    // TODO: Surely there is a more efficient way to do this????
737    // Why is there no memset here? Parallel issues?
738    int numSamples = getNumSamples();    int numSamples = getNumSamples();
739    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
740    DataArrayView& thisView=getPointDataView();    DataTypes::ValueType::size_type n = getNoValues();
   DataArrayView::ValueType::size_type n = thisView.noValues();  
741    double* p;    double* p;
742    int  sampleNo,dataPointNo, i;    int  sampleNo,dataPointNo, i;
743    #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
744    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
745      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
746          p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);          p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
747          for (int i=0; i<n ;++i) p[i]=0.;          for (i=0; i<n ;++i) p[i]=0.;
748      }      }
749    }    }
750  }  }
751    
752    /* Append MPI rank to file name if multiple MPI processes */
753    char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
754      /* Make plenty of room for the mpi_rank number and terminating '\0' */
755      char *newFileName = (char *)malloc(strlen(fileName)+20);
756      strncpy(newFileName, fileName, strlen(fileName)+1);
757      if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
758      return(newFileName);
759    }
760    
761  void  void
762  DataExpanded::dump(const std::string fileName) const  DataExpanded::dump(const std::string fileName) const
763  {  {
    #ifdef PASO_MPI  
    throw DataException("Error - DataExpanded:: dump is not implemented for MPI yet.");  
    #endif  
764     #ifdef USE_NETCDF     #ifdef USE_NETCDF
765     const int ldims=2+DataArrayView::maxRank;     const int ldims=2+DataTypes::maxRank;
766     const NcDim* ncdims[ldims];     const NcDim* ncdims[ldims];
767     NcVar *var, *ids;     NcVar *var, *ids;
768     int rank = getPointDataView().getRank();     int rank = getRank();
769     int type=  getFunctionSpace().getTypeCode();     int type=  getFunctionSpace().getTypeCode();
770     int ndims =0;     int ndims =0;
771     long dims[ldims];     long dims[ldims];
772     const double* d_ptr=&(m_data[0]);     const double* d_ptr=&(m_data[0]);
773     DataArrayView::ShapeType shape = getPointDataView().getShape();     const DataTypes::ShapeType& shape = getShape();
774       int mpi_iam=getFunctionSpace().getDomain().getMPIRank();
775       int mpi_num=getFunctionSpace().getDomain().getMPISize();
776    #ifdef PASO_MPI
777       MPI_Status status;
778    #endif
779    
780     // netCDF error handler     // netCDF error handler
781     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
782     // Create the file.     // Create the file.
783     NcFile dataFile(fileName.c_str(), NcFile::Replace);  #ifdef PASO_MPI
784       if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, getFunctionSpace().getDomain().getMPIComm(), &status);
785    #endif
786       char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
787       NcFile dataFile(newFileName, NcFile::Replace);
788     // check if writing was successful     // check if writing was successful
789     if (!dataFile.is_valid())     if (!dataFile.is_valid())
790          throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");          throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
# Line 681  DataExpanded::dump(const std::string fil Line 798  DataExpanded::dump(const std::string fil
798     if ( rank >0 ) {     if ( rank >0 ) {
799         dims[0]=shape[0];         dims[0]=shape[0];
800         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
801              throw DataException("Error - DataExpanded:: appending ncdimsion 0 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
802     }     }
803     if ( rank >1 ) {     if ( rank >1 ) {
804         dims[1]=shape[1];         dims[1]=shape[1];
805         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
806              throw DataException("Error - DataExpanded:: appending ncdimsion 1 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
807     }     }
808     if ( rank >2 ) {     if ( rank >2 ) {
809         dims[2]=shape[2];         dims[2]=shape[2];
810         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
811              throw DataException("Error - DataExpanded:: appending ncdimsion 2 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
812     }     }
813     if ( rank >3 ) {     if ( rank >3 ) {
814         dims[3]=shape[3];         dims[3]=shape[3];
815         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
816              throw DataException("Error - DataExpanded:: appending ncdimsion 3 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
817     }     }
818     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
819     if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )     if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
# Line 715  DataExpanded::dump(const std::string fil Line 832  DataExpanded::dump(const std::string fil
832          throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");          throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
833     if (! (var->put(d_ptr,dims)) )     if (! (var->put(d_ptr,dims)) )
834          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
835    #ifdef PASO_MPI
836       if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, getFunctionSpace().getDomain().getMPIComm());
837    #endif
838     #else     #else
839     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
840     #endif     #endif
841  }  }
842    
843  void  // void
844    // DataExpanded::setTaggedValue(int tagKey,
845    //                              const DataArrayView& value)
846    // {
847    //   int numSamples = getNumSamples();
848    //   int numDataPointsPerSample = getNumDPPSample();
849    //   int sampleNo,dataPointNo, i;
850    //   DataTypes::ValueType::size_type n = getNoValues();
851    //   double* p,*in=&(value.getData()[0]);
852    //  
853    //   if (value.noValues() != n) {
854    //     throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
855    //   }
856    //
857    //   #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
858    //   for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
859    //     if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
860    //         for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
861    //             p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
862    //             for (i=0; i<n ;++i) p[i]=in[i];
863    //         }
864    //     }
865    //   }
866    // }
867    
868    void  
869  DataExpanded::setTaggedValue(int tagKey,  DataExpanded::setTaggedValue(int tagKey,
870                               const DataArrayView& value)             const DataTypes::ShapeType& pointshape,
871                   const DataTypes::ValueType& value,
872               int dataOffset)
873  {  {
874    int numSamples = getNumSamples();    int numSamples = getNumSamples();
875    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
876    int sampleNo,dataPointNo, i;    int sampleNo,dataPointNo, i;
877    DataArrayView& thisView=getPointDataView();    DataTypes::ValueType::size_type n = getNoValues();
878    DataArrayView::ValueType::size_type n = thisView.noValues();    double* p;
879    double* p,*in=&(value.getData()[0]);    const double* in=&value[0+dataOffset];
880        
881    if (value.noValues() != n) {    if (value.size() != n) {
882      throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");      throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
883    }    }
884    
# Line 740  DataExpanded::setTaggedValue(int tagKey, Line 887  DataExpanded::setTaggedValue(int tagKey,
887      if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {      if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
888          for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {          for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
889              p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);              p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
890              for (int i=0; i<n ;++i) p[i]=in[i];              for (i=0; i<n ;++i) p[i]=in[i];
891          }          }
892      }      }
893    }    }
894  }  }
895    
896    
897    void
898    DataExpanded::reorderByReferenceIDs(int *reference_ids)
899    {
900      int numSamples = getNumSamples();
901      DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
902      int sampleNo, sampleNo2,i;
903      double* p,*p2;
904      register double rtmp;
905      FunctionSpace fs=getFunctionSpace();
906    
907      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
908         const int id_in=reference_ids[sampleNo];
909         const int id=fs.getReferenceIDOfSample(sampleNo);
910         if (id!=id_in) {
911             bool matched=false;
912             for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
913                  if (id == reference_ids[sampleNo2]) {
914                     p=&(m_data[getPointOffset(sampleNo,0)]);
915                     p2=&(m_data[getPointOffset(sampleNo2,0)]);
916                     for (i=0; i<n ;i++) {
917                             rtmp=p[i];
918                             p[i]=p2[i];
919                             p2[i]=rtmp;
920                     }
921                     reference_ids[sampleNo]=id;
922                     reference_ids[sampleNo2]=id_in;
923                     matched=true;
924                     break;
925                  }
926             }
927             if (! matched) {
928                throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
929             }
930         }
931       }
932    }
933    
934    DataTypes::ValueType&
935    DataExpanded::getVector()
936    {
937        return m_data.getData();
938    }
939    
940    const DataTypes::ValueType&
941    DataExpanded::getVector() const
942    {
943        return m_data.getData();
944    }
945    
946  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1358  
changed lines
  Added in v.1800

  ViewVC Help
Powered by ViewVC 1.1.26