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

Legend:
Removed from v.1312  
changed lines
  Added in v.1799

  ViewVC Help
Powered by ViewVC 1.1.26