/[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 1808 by jfenwick, Thu Sep 25 03:14:56 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  {  {
   DataArrayView::ShapeType tempShape;  
   //  
   // extract the shape of the python numarray  
   for (int i=0; i<value.getrank(); i++) {  
     tempShape.push_back(extract<int>(value.getshape()[i]));  
   }  
42    //    //
43    // initialise the data array for this object    // initialise the data array for this object
44    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    initialise(what.getNumSamples(),what.getNumDPPSample());
45    //    //
46    // copy the given value to every data point    // copy the given value to every data point
47    copy(value);    copy(value);
48  }  }
49    
50  DataExpanded::DataExpanded(const DataExpanded& other)  DataExpanded::DataExpanded(const DataExpanded& other)
51    : DataAbstract(other.getFunctionSpace()),    : DataAbstract(other.getFunctionSpace(), other.getShape()),
52    m_data(other.m_data)    m_data(other.m_data)
53  {  {
   //  
   // create the view for the data  
   DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());  
   setPointDataView(temp);  
54  }  }
55    
56  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
57    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace(), other.getShape())
58  {  {
59    //    //
60    // initialise the data array for this object    // initialise the data array for this object
61    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
62    //    //
63    // DataConstant only has one value, copy this to every data point    // DataConstant only has one value, copy this to every data point
64    copy(other.getPointDataView());    copy(other);
65  }  }
66    
67  DataExpanded::DataExpanded(const DataTagged& other)  DataExpanded::DataExpanded(const DataTagged& other)
68    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace(), other.getShape())
69  {  {
70    //    //
71    // initialise the data array for this object    // initialise the data array for this object
72    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
73    //    //
74    // 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
75    // value from the given DataTagged object    // value from the given DataTagged object
76    int i,j;    int i,j;
77    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
78    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
79    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
80    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
81      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
82        try {        try {
83          getPointDataView().copy(getPointOffset(i,j),             DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(),
84                                  other.getPointDataView(),                                  other.getVector(),
85                                  other.getPointOffset(i,j));                                  other.getPointOffset(i,j));
86        }        }
87        catch (std::exception& e) {        catch (std::exception& e) {
# Line 96  DataExpanded::DataExpanded(const DataTag Line 92  DataExpanded::DataExpanded(const DataTag
92  }  }
93    
94  DataExpanded::DataExpanded(const DataExpanded& other,  DataExpanded::DataExpanded(const DataExpanded& other,
95                             const DataArrayView::RegionType& region)                             const DataTypes::RegionType& region)
96    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
97  {  {
98    //    //
99    // get the shape of the slice    // get the shape of the slice
100    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));  //   DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
101    //    //
102    // initialise this Data object to the shape of the slice    // initialise this Data object to the shape of the slice
103    initialise(shape,other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
104    //    //
105    // copy the data    // copy the data
106    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
107    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
108    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
109    int i,j;    int i,j;
110    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
111    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
112      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
113        try {        try {
114          getPointDataView().copySlice(getPointOffset(i,j),  //         getPointDataView().copySlice(getPointOffset(i,j),
115                                       other.getPointDataView(),  //                                      other.getPointDataView(),
116    //                                      other.getPointOffset(i,j),
117    //                                      region_loop_range);
118            DataTypes::copySlice(getVector(),getShape(),getPointOffset(i,j),
119                                         other.getVector(),
120                         other.getShape(),
121                                       other.getPointOffset(i,j),                                       other.getPointOffset(i,j),
122                                       region_loop_range);                                       region_loop_range);
123        }        }
# Line 127  DataExpanded::DataExpanded(const DataExp Line 128  DataExpanded::DataExpanded(const DataExp
128    }    }
129  }  }
130    
131  DataExpanded::DataExpanded(const DataArrayView& value,  // DataExpanded::DataExpanded(const DataArrayView& value,
132                             const FunctionSpace& what)  //                            const FunctionSpace& what)
133    : DataAbstract(what)  //   : DataAbstract(what)
134  {  // {
135    //  //   //
136    // get the shape of the given data value  //   // get the shape of the given data value
137    DataArrayView::ShapeType tempShape=value.getShape();  //   DataTypes::ShapeType tempShape=value.getShape();
138    //  //   //
139    // initialise this Data object to the shape of the given data value  //   // initialise this Data object to the shape of the given data value
140    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());  //   initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
141    //  //   //
142    // copy the given value to every data point  //   // copy the given value to every data point
143    copy(value);  //   copy(value);
144  }  // }
145    
146  DataExpanded::DataExpanded(const FunctionSpace& what,  DataExpanded::DataExpanded(const FunctionSpace& what,
147                             const DataArrayView::ShapeType &shape,                             const DataTypes::ShapeType &shape,
148                             const DataArrayView::ValueType &data)                             const DataTypes::ValueType &data)
149    : DataAbstract(what)    : DataAbstract(what,shape)
150  {  {
151    //    EsysAssert(data.size()%getNoValues()==0,
152    // create the view of the data                   "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
153    initialise(shape,what.getNumSamples(),what.getNumDPPSample());  
154    //    if (data.size()==getNoValues())
155    // copy the data in the correct format    {
156    m_data.getData()=data;       ValueType& vec=m_data.getData();
157         //
158         // create the view of the data
159         initialise(what.getNumSamples(),what.getNumDPPSample());
160         // now we copy this value to all elements
161         for (int i=0;i<getLength();)
162         {
163        for (int j=0;j<getNoValues();++j,++i)
164        {
165            vec[i]=data[j];
166        }
167         }
168      }
169      else
170      {
171         //
172         // copy the data in the correct format
173         m_data.getData()=data;
174      }
175    
176    
177  }  }
178    
179  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
# Line 160  DataExpanded::~DataExpanded() Line 181  DataExpanded::~DataExpanded()
181  }  }
182    
183  DataAbstract*  DataAbstract*
184  DataExpanded::getSlice(const DataArrayView::RegionType& region) const  DataExpanded::deepCopy()
185    {
186      return new DataExpanded(*this);
187    }
188    
189    
190    DataAbstract*
191    DataExpanded::getSlice(const DataTypes::RegionType& region) const
192  {  {
193    return new DataExpanded(*this,region);    return new DataExpanded(*this,region);
194  }  }
195    
196  void  void
197  DataExpanded::setSlice(const DataAbstract* value,  DataExpanded::setSlice(const DataAbstract* value,
198                         const DataArrayView::RegionType& region)                         const DataTypes::RegionType& region)
199  {  {
200    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
201    if (tempDataExp==0) {    if (tempDataExp==0) {
# Line 175  DataExpanded::setSlice(const DataAbstrac Line 203  DataExpanded::setSlice(const DataAbstrac
203    }    }
204    //    //
205    // get shape of slice    // get shape of slice
206    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
207    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
208    //    //
209    // check shape    // check shape
210    if (getPointDataView().getRank()!=region.size()) {    if (getRank()!=region.size()) {
211      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
212    }    }
213    if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {    if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
214      throw DataException (value->getPointDataView().createShapeErrorMessage(      throw DataException (DataTypes::createShapeErrorMessage(
215          "Error - Couldn't copy slice due to shape mismatch.",shape));          "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
216    }    }
217    //    //
218    // copy the data from the slice into this object    // copy the data from the slice into this object
219    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
220    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
221    int i, j;    int i, j;
222      ValueType& vec=getVector();
223      const ShapeType& mshape=getShape();
224      const ValueType& tVec=tempDataExp->getVector();
225      const ShapeType& tShape=tempDataExp->getShape();
226    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
227    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
228      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
229        getPointDataView().copySliceFrom(getPointOffset(i,j),  /*      getPointDataView().copySliceFrom(getPointOffset(i,j),
230                                         tempDataExp->getPointDataView(),                                         tempDataExp->getPointDataView(),
231                                         tempDataExp->getPointOffset(i,j),                                         tempDataExp->getPointOffset(i,j),
232                                           region_loop_range);*/
233            DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
234                                           tVec,
235                           tShape,
236                                           tempDataExp->getPointOffset(i,j),
237                                         region_loop_range);                                         region_loop_range);
238    
239      }      }
240    }    }
241  }  }
242    
243  void  void
244  DataExpanded::copy(const DataArrayView& value)  DataExpanded::copy(const DataConstant& value)
245  {  {
246      EsysAssert((checkShape(getShape(), value.getShape())),
247                     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
248    
249    //    //
250    // copy a single value to every data point in this object    // copy a single value to every data point in this object
251    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
# Line 217  DataExpanded::copy(const DataArrayView& Line 258  DataExpanded::copy(const DataArrayView&
258        // DOASSERT is on which of course will play        // DOASSERT is on which of course will play
259        // havoc with the omp threads. Run single threaded        // havoc with the omp threads. Run single threaded
260        // if using DOASSERT.        // if using DOASSERT.
261        getPointDataView().copy(getPointOffset(i,j),value);        //getPointDataView().copy(getPointOffset(i,j),value);
262          DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(), value.getVector(), 0);
263      }      }
264    }    }
265  }  }
266    
267    
268    // void
269    // DataExpanded::copy(const DataArrayView& value)
270    // {
271    //   //
272    //   // copy a single value to every data point in this object
273    //   int nRows=m_data.getNumRows();
274    //   int nCols=m_data.getNumCols();
275    //   int i,j;
276    //   #pragma omp parallel for private(i,j) schedule(static)
277    //   for (i=0;i<nRows;i++) {
278    //     for (j=0;j<nCols;j++) {
279    //       // NOTE: An exception may be thown from this call if
280    //       // DOASSERT is on which of course will play
281    //       // havoc with the omp threads. Run single threaded
282    //       // if using DOASSERT.
283    //       getPointDataView().copy(getPointOffset(i,j),value);
284    //     }
285    //   }
286    // }
287    
288  void  void
289  DataExpanded::copy(const boost::python::numeric::array& value)  DataExpanded::copy(const boost::python::numeric::array& value)
290  {  {
291    
292    // extract the shape of the numarray    // extract the shape of the numarray
293    DataArrayView::ShapeType tempShape;    DataTypes::ShapeType tempShape;
294    for (int i=0; i < value.getrank(); i++) {    for (int i=0; i < value.getrank(); i++) {
295      tempShape.push_back(extract<int>(value.getshape()[i]));      tempShape.push_back(extract<int>(value.getshape()[i]));
296    }    }
297    
298    // get the space for the data vector    // get the space for the data vector
299    int len = DataArrayView::noValues(tempShape);  //   int len = DataTypes::noValues(tempShape);
300    DataVector temp_data(len, 0.0, len);  //   DataVector temp_data(len, 0.0, len);
301    DataArrayView temp_dataView(temp_data, tempShape);  //   DataArrayView temp_dataView(temp_data, tempShape);
302    temp_dataView.copy(value);  //   temp_dataView.copy(value);
303    
304    //    //
305    // check the input shape matches this shape    // check the input shape matches this shape
306    if (!getPointDataView().checkShape(temp_dataView.getShape())) {    if (!DataTypes::checkShape(getShape(),tempShape)) {
307      throw DataException(getPointDataView().createShapeErrorMessage(      throw DataException(DataTypes::createShapeErrorMessage(
308                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
309                          temp_dataView.getShape()));                          tempShape,getShape()));
310    }    }
311    //    //
312    // now copy over the data    // now copy over the data
313    copy(temp_dataView);    //copy(temp_dataView);
314      getVector().copyFromNumArray(value);
315  }  }
316    
317    
318  void  void
319  DataExpanded::initialise(const DataArrayView::ShapeType& shape,  DataExpanded::initialise(int noSamples,
                          int noSamples,  
320                           int noDataPointsPerSample)                           int noDataPointsPerSample)
321  {  {
322      if (noSamples==0)     //retain the default empty object
323      {
324         return;
325      }
326    //    //
327    // resize data array to the required size    // resize data array to the required size
328    m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));    m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
   //  
   // create the data view of the data array  
   DataArrayView temp(m_data.getData(),shape);  
   setPointDataView(temp);  
329  }  }
330    
331  string  string
# Line 269  DataExpanded::toString() const Line 333  DataExpanded::toString() const
333  {  {
334    stringstream temp;    stringstream temp;
335    FunctionSpace fs=getFunctionSpace();    FunctionSpace fs=getFunctionSpace();
336    //  
337    // create a temporary view as the offset will be changed    int offset=0;
   DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());  
338    for (int i=0;i<m_data.getNumRows();i++) {    for (int i=0;i<m_data.getNumRows();i++) {
339      for (int j=0;j<m_data.getNumCols();j++) {      for (int j=0;j<m_data.getNumCols();j++) {
340        tempView.setOffset(getPointOffset(i,j));        offset=getPointOffset(i,j);
341        stringstream suffix;        stringstream suffix;
342        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
343        temp << tempView.toString(suffix.str());        temp << DataTypes::pointToString(getVector(),getShape(),offset,suffix.str());
344        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
345          temp << endl;          temp << endl;
346        }        }
347      }      }
348    }    }
349      string result=temp.str();
350      if (result.empty())
351      {
352        return "(data contains no samples)\n";
353      }
354    return temp.str();    return temp.str();
355  }  }
356    
357  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
358  DataExpanded::getPointOffset(int sampleNo,  DataExpanded::getPointOffset(int sampleNo,
359                               int dataPointNo) const                               int dataPointNo) const
360  {  {
361    return m_data.index(sampleNo,dataPointNo);    return m_data.index(sampleNo,dataPointNo);
362  }  }
363    
364  DataArrayView  DataTypes::ValueType::size_type
 DataExpanded::getDataPoint(int sampleNo,  
                            int dataPointNo)  
 {  
   DataArrayView temp(m_data.getData(),getPointDataView().getShape(),getPointOffset(sampleNo,dataPointNo));  
   return temp;  
 }  
   
 DataArrayView::ValueType::size_type  
365  DataExpanded::getLength() const  DataExpanded::getLength() const
366  {  {
367    return m_data.size();    return m_data.size();
368  }  }
369    
 int  
 DataExpanded::archiveData(ofstream& archiveFile,  
                           const DataArrayView::ValueType::size_type noValues) const  
 {  
   return(m_data.archiveData(archiveFile, noValues));  
 }  
370    
 int  
 DataExpanded::extractData(ifstream& archiveFile,  
                           const DataArrayView::ValueType::size_type noValues)  
 {  
   return(m_data.extractData(archiveFile, noValues));  
 }  
371    
372  void  void
373  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 375  DataExpanded::copyToDataPoint(const int
375    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
376    int numSamples = getNumSamples();    int numSamples = getNumSamples();
377    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
378    int dataPointRank = getPointDataView().getRank();    int dataPointRank = getRank();
379    ShapeType dataPointShape = getPointDataView().getShape();    ShapeType dataPointShape = getShape();
380    if (numSamples*numDataPointsPerSample > 0) {    if (numSamples*numDataPointsPerSample > 0) {
381       //TODO: global error handling       //TODO: global error handling
382       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
# Line 337  DataExpanded::copyToDataPoint(const int Line 385  DataExpanded::copyToDataPoint(const int
385       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
386             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
387       }       }
388       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
389         ValueType& vec=getVector();
390       if (dataPointRank==0) {       if (dataPointRank==0) {
391           dataPointView()=value;           vec[0]=value;
392       } else if (dataPointRank==1) {       } else if (dataPointRank==1) {
393          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
394              dataPointView(i)=value;              vec[offset+i]=value;
395          }          }
396       } else if (dataPointRank==2) {       } else if (dataPointRank==2) {
397          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
398             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
399                dataPointView(i,j)=value;                vec[offset+getRelIndex(dataPointShape,i,j)]=value;
400             }             }
401          }          }
402       } else if (dataPointRank==3) {       } else if (dataPointRank==3) {
403          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
404             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
405                for (int k=0; k<dataPointShape[2]; k++) {                for (int k=0; k<dataPointShape[2]; k++) {
406                   dataPointView(i,j,k)=value;                   vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
407                }                }
408             }             }
409          }          }
# Line 363  DataExpanded::copyToDataPoint(const int Line 412  DataExpanded::copyToDataPoint(const int
412             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
413               for (int k=0; k<dataPointShape[2]; k++) {               for (int k=0; k<dataPointShape[2]; k++) {
414                 for (int l=0; l<dataPointShape[3]; l++) {                 for (int l=0; l<dataPointShape[3]; l++) {
415                    dataPointView(i,j,k,l)=value;                    vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
416                 }                 }
417               }               }
418             }             }
# Line 377  DataExpanded::copyToDataPoint(const int Line 426  DataExpanded::copyToDataPoint(const int
426    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
427    int numSamples = getNumSamples();    int numSamples = getNumSamples();
428    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
429    int dataPointRank = getPointDataView().getRank();    int dataPointRank = getRank();
430    ShapeType dataPointShape = getPointDataView().getShape();    const ShapeType& shape = getShape();
431    //    //
432    // check rank:    // check rank:
433    if (value.getrank()!=dataPointRank)    if (value.getrank()!=dataPointRank)
# Line 391  DataExpanded::copyToDataPoint(const int Line 440  DataExpanded::copyToDataPoint(const int
440       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
441             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
442       }       }
443       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);       ValueType& vec=getVector();
444       if (dataPointRank==0) {       if (dataPointRank==0) {
445           dataPointView()=extract<double>(value[0]);           vec[0]=extract<double>(value[0]);
446       } else if (dataPointRank==1) {       } else if (dataPointRank==1) {
447          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<shape[0]; i++) {
448              dataPointView(i)=extract<double>(value[i]);              vec[i]=extract<double>(value[i]);
449          }          }
450       } else if (dataPointRank==2) {       } else if (dataPointRank==2) {
451          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<shape[0]; i++) {
452             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<shape[1]; j++) {
453                dataPointView(i,j)=extract<double>(value[i][j]);                vec[getRelIndex(shape,i,j)]=extract<double>(value[i][j]);
454             }             }
455          }          }
456       } else if (dataPointRank==3) {       } else if (dataPointRank==3) {
457          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<shape[0]; i++) {
458             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<shape[1]; j++) {
459                for (int k=0; k<dataPointShape[2]; k++) {                for (int k=0; k<shape[2]; k++) {
460                   dataPointView(i,j,k)=extract<double>(value[i][j][k]);                   vec[getRelIndex(shape,i,j,k)]=extract<double>(value[i][j][k]);
461                }                }
462             }             }
463          }          }
464       } else if (dataPointRank==4) {       } else if (dataPointRank==4) {
465           for (int i=0; i<dataPointShape[0]; i++) {           for (int i=0; i<shape[0]; i++) {
466             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<shape[1]; j++) {
467               for (int k=0; k<dataPointShape[2]; k++) {               for (int k=0; k<shape[2]; k++) {
468                 for (int l=0; l<dataPointShape[3]; l++) {                 for (int l=0; l<shape[3]; l++) {
469                    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]);
470                 }                 }
471               }               }
472             }             }
# Line 431  DataExpanded::copyAll(const boost::pytho Line 480  DataExpanded::copyAll(const boost::pytho
480    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
481    int numSamples = getNumSamples();    int numSamples = getNumSamples();
482    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
483    int dataPointRank = getPointDataView().getRank();    int dataPointRank = getRank();
484    ShapeType dataPointShape = getPointDataView().getShape();    const ShapeType& dataPointShape = getShape();
485    //    //
486    // check rank:    // check rank:
487    if (value.getrank()!=dataPointRank+1)    if (value.getrank()!=dataPointRank+1)
# Line 440  DataExpanded::copyAll(const boost::pytho Line 489  DataExpanded::copyAll(const boost::pytho
489    if (value.getshape()[0]!=numSamples*numDataPointsPerSample)    if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
490         throw DataException("leading dimension of numarray is too small");         throw DataException("leading dimension of numarray is too small");
491    //    //
492      ValueType& vec=getVector();
493    int dataPoint = 0;    int dataPoint = 0;
494    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
495      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
496        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);        ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
497        if (dataPointRank==0) {        if (dataPointRank==0) {
498           dataPointView()=extract<double>(value[dataPoint]);           vec[offset]=extract<double>(value[dataPoint]);
499        } else if (dataPointRank==1) {        } else if (dataPointRank==1) {
500           for (int i=0; i<dataPointShape[0]; i++) {           for (int i=0; i<dataPointShape[0]; i++) {
501              dataPointView(i)=extract<double>(value[dataPoint][i]);              vec[offset+i]=extract<double>(value[dataPoint][i]);
502           }           }
503        } else if (dataPointRank==2) {        } else if (dataPointRank==2) {
504           for (int i=0; i<dataPointShape[0]; i++) {           for (int i=0; i<dataPointShape[0]; i++) {
505             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
506               dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);           vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
507             }             }
508           }           }
509         } else if (dataPointRank==3) {         } else if (dataPointRank==3) {
510           for (int i=0; i<dataPointShape[0]; i++) {           for (int i=0; i<dataPointShape[0]; i++) {
511             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
512               for (int k=0; k<dataPointShape[2]; k++) {               for (int k=0; k<dataPointShape[2]; k++) {
513                   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]);
514               }               }
515             }             }
516           }           }
# Line 469  DataExpanded::copyAll(const boost::pytho Line 519  DataExpanded::copyAll(const boost::pytho
519             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
520               for (int k=0; k<dataPointShape[2]; k++) {               for (int k=0; k<dataPointShape[2]; k++) {
521                 for (int l=0; l<dataPointShape[3]; l++) {                 for (int l=0; l<dataPointShape[3]; l++) {
522                   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]);
523                 }                 }
524               }               }
525             }             }
# Line 489  DataExpanded::symmetric(DataAbstract* ev Line 539  DataExpanded::symmetric(DataAbstract* ev
539    if (temp_ev==0) {    if (temp_ev==0) {
540      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).");
541    }    }
542    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
543    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
544      ValueType& evVec=temp_ev->getVector();
545      const ShapeType& evShape=temp_ev->getShape();
546    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
547    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
548      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
549           DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
550                                      evView,ev->getPointOffset(sampleNo,dataPointNo));                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
551      }      }
552    }    }
553  }  }
# Line 509  DataExpanded::nonsymmetric(DataAbstract* Line 561  DataExpanded::nonsymmetric(DataAbstract*
561    if (temp_ev==0) {    if (temp_ev==0) {
562      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).");
563    }    }
564    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
565    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
566      ValueType& evVec=temp_ev->getVector();
567      const ShapeType& evShape=temp_ev->getShape();
568    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
569    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
570      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
571           DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
572                                      evView,ev->getPointOffset(sampleNo,dataPointNo));                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
573      }      }
574    }    }
575  }  }
# Line 529  DataExpanded::trace(DataAbstract* ev, in Line 583  DataExpanded::trace(DataAbstract* ev, in
583    if (temp_ev==0) {    if (temp_ev==0) {
584      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).");
585    }    }
586    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
587    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
588      ValueType& evVec=temp_ev->getVector();
589      const ShapeType& evShape=temp_ev->getShape();
590    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
591    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
592      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
593           DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
594                                      evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
595      }      }
596    }    }
597  }  }
# Line 550  DataExpanded::transpose(DataAbstract* ev Line 606  DataExpanded::transpose(DataAbstract* ev
606    if (temp_ev==0) {    if (temp_ev==0) {
607      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).");
608    }    }
609    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
610    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
611      ValueType& evVec=temp_ev->getVector();
612      const ShapeType& evShape=temp_ev->getShape();
613    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
614    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
615      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
616           DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
617                                      evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
618      }      }
619    }    }
620  }  }
# Line 571  DataExpanded::swapaxes(DataAbstract* ev, Line 629  DataExpanded::swapaxes(DataAbstract* ev,
629    if (temp_ev==0) {    if (temp_ev==0) {
630      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).");
631    }    }
632    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
633    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
634      ValueType& evVec=temp_ev->getVector();
635      const ShapeType& evShape=temp_ev->getShape();
636    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
637    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
638      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
639           DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
640                                      evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
641      }      }
642    }    }
643  }  }
# Line 591  DataExpanded::eigenvalues(DataAbstract* Line 651  DataExpanded::eigenvalues(DataAbstract*
651    if (temp_ev==0) {    if (temp_ev==0) {
652      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).");
653    }    }
654    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
655    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
656      ValueType& evVec=temp_ev->getVector();
657      const ShapeType& evShape=temp_ev->getShape();
658    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
659    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
660      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
661           DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
662                                      evView,ev->getPointOffset(sampleNo,dataPointNo));                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
663      }      }
664    }    }
665  }  }
# Line 615  DataExpanded::eigenvalues_and_eigenvecto Line 677  DataExpanded::eigenvalues_and_eigenvecto
677    if (temp_V==0) {    if (temp_V==0) {
678      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).");
679    }    }
680    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
681    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
682    DataArrayView& VView=V->getPointDataView();    ValueType& evVec=temp_ev->getVector();
683      const ShapeType& evShape=temp_ev->getShape();
684      ValueType& VVec=temp_V->getVector();
685      const ShapeType& VShape=temp_V->getShape();
686    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
687    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
688      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
689           DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
690                                                       evView,ev->getPointOffset(sampleNo,dataPointNo),                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
691                                                       VView,V->getPointOffset(sampleNo,dataPointNo),                                      VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
                                                      tol);  
692      }      }
693    }    }
694  }  }
695    
696  void  void
697  DataExpanded::setToZero(){  DataExpanded::setToZero(){
698    // TODO: Surely there is a more efficient way to do this????
699    // Why is there no memset here? Parallel issues?
700    int numSamples = getNumSamples();    int numSamples = getNumSamples();
701    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
702    DataArrayView& thisView=getPointDataView();    DataTypes::ValueType::size_type n = getNoValues();
   DataArrayView::ValueType::size_type n = thisView.noValues();  
703    double* p;    double* p;
704    int  sampleNo,dataPointNo, i;    int  sampleNo,dataPointNo, i;
705    #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
706    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
707      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
708          p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);          p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
709          for (int i=0; i<n ;++i) p[i]=0.;          for (i=0; i<n ;++i) p[i]=0.;
710      }      }
711    }    }
712  }  }
713    
714    /* Append MPI rank to file name if multiple MPI processes */
715    char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
716      /* Make plenty of room for the mpi_rank number and terminating '\0' */
717      char *newFileName = (char *)malloc(strlen(fileName)+20);
718      strncpy(newFileName, fileName, strlen(fileName)+1);
719      if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
720      return(newFileName);
721    }
722    
723  void  void
724  DataExpanded::dump(const std::string fileName) const  DataExpanded::dump(const std::string fileName) const
725  {  {
    #ifdef PASO_MPI  
    throw DataException("Error - DataExpanded:: dump is not implemented for MPI yet.");  
    #endif  
726     #ifdef USE_NETCDF     #ifdef USE_NETCDF
727     const int ldims=2+DataArrayView::maxRank;     const int ldims=2+DataTypes::maxRank;
728     const NcDim* ncdims[ldims];     const NcDim* ncdims[ldims];
729     NcVar *var, *ids;     NcVar *var, *ids;
730     int rank = getPointDataView().getRank();     int rank = getRank();
731     int type=  getFunctionSpace().getTypeCode();     int type=  getFunctionSpace().getTypeCode();
732     int ndims =0;     int ndims =0;
733     long dims[ldims];     long dims[ldims];
734     const double* d_ptr=&(m_data[0]);     const double* d_ptr=&(m_data[0]);
735     DataArrayView::ShapeType shape = getPointDataView().getShape();     const DataTypes::ShapeType& shape = getShape();
736       int mpi_iam=getFunctionSpace().getDomain().getMPIRank();
737       int mpi_num=getFunctionSpace().getDomain().getMPISize();
738    #ifdef PASO_MPI
739       MPI_Status status;
740    #endif
741    
742     // netCDF error handler     // netCDF error handler
743     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
744     // Create the file.     // Create the file.
745     NcFile dataFile(fileName.c_str(), NcFile::Replace);  #ifdef PASO_MPI
746       if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
747    #endif
748       char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
749       NcFile dataFile(newFileName, NcFile::Replace);
750     // check if writing was successful     // check if writing was successful
751     if (!dataFile.is_valid())     if (!dataFile.is_valid())
752          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 760  DataExpanded::dump(const std::string fil
760     if ( rank >0 ) {     if ( rank >0 ) {
761         dims[0]=shape[0];         dims[0]=shape[0];
762         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
763              throw DataException("Error - DataExpanded:: appending ncdimsion 0 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
764     }     }
765     if ( rank >1 ) {     if ( rank >1 ) {
766         dims[1]=shape[1];         dims[1]=shape[1];
767         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
768              throw DataException("Error - DataExpanded:: appending ncdimsion 1 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
769     }     }
770     if ( rank >2 ) {     if ( rank >2 ) {
771         dims[2]=shape[2];         dims[2]=shape[2];
772         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
773              throw DataException("Error - DataExpanded:: appending ncdimsion 2 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
774     }     }
775     if ( rank >3 ) {     if ( rank >3 ) {
776         dims[3]=shape[3];         dims[3]=shape[3];
777         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
778              throw DataException("Error - DataExpanded:: appending ncdimsion 3 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
779     }     }
780     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
781     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 705  DataExpanded::dump(const std::string fil Line 784  DataExpanded::dump(const std::string fil
784     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
785              throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
786    
787     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )     if (getFunctionSpace().getNumSamples()>0)
788       {
789    
790         if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
791          throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");          throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
792     const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();       const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
793     if (! (ids->put(ids_p,dims[rank+1])) )       if (! (ids->put(ids_p,dims[rank+1])) )
794          throw DataException("Error - DataExpanded:: copy reference id  to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy reference id  to netCDF buffer failed.");
795         if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
    if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )  
796          throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");          throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
797     if (! (var->put(d_ptr,dims)) )       if (! (var->put(d_ptr,dims)) )
798          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
799       }
800    #ifdef PASO_MPI
801       if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
802    #endif
803     #else     #else
804     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.");
805     #endif     #endif
806  }  }
807    
808  void  void  
809  DataExpanded::setTaggedValue(int tagKey,  DataExpanded::setTaggedValue(int tagKey,
810                               const DataArrayView& value)             const DataTypes::ShapeType& pointshape,
811                   const DataTypes::ValueType& value,
812               int dataOffset)
813  {  {
814    int numSamples = getNumSamples();    int numSamples = getNumSamples();
815    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
816    int sampleNo,dataPointNo, i;    int sampleNo,dataPointNo, i;
817    DataArrayView& thisView=getPointDataView();    DataTypes::ValueType::size_type n = getNoValues();
818    DataArrayView::ValueType::size_type n = thisView.noValues();    double* p;
819    double* p,*in=&(value.getData()[0]);    const double* in=&value[0+dataOffset];
820        
821    if (value.noValues() != n) {    if (value.size() != n) {
822      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.");
823    }    }
824    
# Line 740  DataExpanded::setTaggedValue(int tagKey, Line 827  DataExpanded::setTaggedValue(int tagKey,
827      if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {      if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
828          for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {          for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
829              p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);              p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
830              for (int i=0; i<n ;++i) p[i]=in[i];              for (i=0; i<n ;++i) p[i]=in[i];
831          }          }
832      }      }
833    }    }
834  }  }
835    
836    
837    void
838    DataExpanded::reorderByReferenceIDs(int *reference_ids)
839    {
840      int numSamples = getNumSamples();
841      DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
842      int sampleNo, sampleNo2,i;
843      double* p,*p2;
844      register double rtmp;
845      FunctionSpace fs=getFunctionSpace();
846    
847      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
848         const int id_in=reference_ids[sampleNo];
849         const int id=fs.getReferenceIDOfSample(sampleNo);
850         if (id!=id_in) {
851             bool matched=false;
852             for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
853                  if (id == reference_ids[sampleNo2]) {
854                     p=&(m_data[getPointOffset(sampleNo,0)]);
855                     p2=&(m_data[getPointOffset(sampleNo2,0)]);
856                     for (i=0; i<n ;i++) {
857                             rtmp=p[i];
858                             p[i]=p2[i];
859                             p2[i]=rtmp;
860                     }
861                     reference_ids[sampleNo]=id;
862                     reference_ids[sampleNo2]=id_in;
863                     matched=true;
864                     break;
865                  }
866             }
867             if (! matched) {
868                throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
869             }
870         }
871       }
872    }
873    
874    DataTypes::ValueType&
875    DataExpanded::getVector()
876    {
877        return m_data.getData();
878    }
879    
880    const DataTypes::ValueType&
881    DataExpanded::getVector() const
882    {
883        return m_data.getData();
884    }
885    
886  }  // end of namespace  }  // end of namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26