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

temp/escript/src/DataExpanded.cpp revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC trunk/escript/src/DataExpanded.cpp revision 2005 by jfenwick, Mon Nov 10 01:21:39 2008 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2008 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15    #include "Data.h"
16  #include "DataExpanded.h"  #include "DataExpanded.h"
17  #include "DataException.h"  #include "DataException.h"
18  #include "DataConstant.h"  #include "DataConstant.h"
# 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)    : parent(what,DataTypes::shapeFromNumArray(value))
40  {  {
   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]));  
   }  
41    //    //
42    // initialise the data array for this object    // initialise the data array for this object
43    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    initialise(what.getNumSamples(),what.getNumDPPSample());
44    //    //
45    // copy the given value to every data point    // copy the given value to every data point
46    copy(value);    copy(value);
47  }  }
48    
49  DataExpanded::DataExpanded(const DataExpanded& other)  DataExpanded::DataExpanded(const DataExpanded& other)
50    : DataAbstract(other.getFunctionSpace()),    : parent(other.getFunctionSpace(), other.getShape()),
51    m_data(other.m_data)    m_data(other.m_data)
52  {  {
   //  
   // create the view for the data  
   DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());  
   setPointDataView(temp);  
53  }  }
54    
55  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
56    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(), other.getShape())
57  {  {
58    //    //
59    // initialise the data array for this object    // initialise the data array for this object
60    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
61    //    //
62    // DataConstant only has one value, copy this to every data point    // DataConstant only has one value, copy this to every data point
63    copy(other.getPointDataView());    copy(other);
64  }  }
65    
66  DataExpanded::DataExpanded(const DataTagged& other)  DataExpanded::DataExpanded(const DataTagged& other)
67    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(), other.getShape())
68  {  {
69    //    //
70    // initialise the data array for this object    // initialise the data array for this object
71    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
72    //    //
73    // 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
74    // value from the given DataTagged object    // value from the given DataTagged object
75    int i,j;    int i,j;
76    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
77    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
78    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
79    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
80      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
81        try {        try {
82          getPointDataView().copy(getPointOffset(i,j),             DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(),
83                                  other.getPointDataView(),                                  other.getVector(),
84                                  other.getPointOffset(i,j));                                  other.getPointOffset(i,j));
85        }        }
86        catch (std::exception& e) {        catch (std::exception& e) {
# Line 96  DataExpanded::DataExpanded(const DataTag Line 91  DataExpanded::DataExpanded(const DataTag
91  }  }
92    
93  DataExpanded::DataExpanded(const DataExpanded& other,  DataExpanded::DataExpanded(const DataExpanded& other,
94                             const DataArrayView::RegionType& region)                             const DataTypes::RegionType& region)
95    : DataAbstract(other.getFunctionSpace())    : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
96  {  {
97    //    //
98    // get the shape of the slice    // get the shape of the slice
99    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));  //   DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
100    //    //
101    // initialise this Data object to the shape of the slice    // initialise this Data object to the shape of the slice
102    initialise(shape,other.getNumSamples(),other.getNumDPPSample());    initialise(other.getNumSamples(),other.getNumDPPSample());
103    //    //
104    // copy the data    // copy the data
105    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
106    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
107    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
108    int i,j;    int i,j;
109    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
110    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
111      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
112        try {        try {
113          getPointDataView().copySlice(getPointOffset(i,j),  //         getPointDataView().copySlice(getPointOffset(i,j),
114                                       other.getPointDataView(),  //                                      other.getPointDataView(),
115    //                                      other.getPointOffset(i,j),
116    //                                      region_loop_range);
117            DataTypes::copySlice(getVector(),getShape(),getPointOffset(i,j),
118                                         other.getVector(),
119                         other.getShape(),
120                                       other.getPointOffset(i,j),                                       other.getPointOffset(i,j),
121                                       region_loop_range);                                       region_loop_range);
122        }        }
# Line 127  DataExpanded::DataExpanded(const DataExp Line 127  DataExpanded::DataExpanded(const DataExp
127    }    }
128  }  }
129    
130  DataExpanded::DataExpanded(const DataArrayView& value,  // DataExpanded::DataExpanded(const DataArrayView& value,
131                             const FunctionSpace& what)  //                            const FunctionSpace& what)
132    : DataAbstract(what)  //   : DataAbstract(what)
133  {  // {
134    //  //   //
135    // get the shape of the given data value  //   // get the shape of the given data value
136    DataArrayView::ShapeType tempShape=value.getShape();  //   DataTypes::ShapeType tempShape=value.getShape();
137    //  //   //
138    // initialise this Data object to the shape of the given data value  //   // initialise this Data object to the shape of the given data value
139    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());  //   initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
140    //  //   //
141    // copy the given value to every data point  //   // copy the given value to every data point
142    copy(value);  //   copy(value);
143  }  // }
144    
145  DataExpanded::DataExpanded(const FunctionSpace& what,  DataExpanded::DataExpanded(const FunctionSpace& what,
146                             const DataArrayView::ShapeType &shape,                             const DataTypes::ShapeType &shape,
147                             const DataArrayView::ValueType &data)                             const DataTypes::ValueType &data)
148    : DataAbstract(what)    : parent(what,shape)
149  {  {
150    //    EsysAssert(data.size()%getNoValues()==0,
151    // create the view of the data                   "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
152    initialise(shape,what.getNumSamples(),what.getNumDPPSample());  
153    //    if (data.size()==getNoValues())
154    // copy the data in the correct format    {
155    m_data.getData()=data;       ValueType& vec=m_data.getData();
156         //
157         // create the view of the data
158         initialise(what.getNumSamples(),what.getNumDPPSample());
159         // now we copy this value to all elements
160         for (int i=0;i<getLength();)
161         {
162        for (int j=0;j<getNoValues();++j,++i)
163        {
164            vec[i]=data[j];
165        }
166         }
167      }
168      else
169      {
170         //
171         // copy the data in the correct format
172         m_data.getData()=data;
173      }
174    
175    
176  }  }
177    
178  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
# Line 160  DataExpanded::~DataExpanded() Line 180  DataExpanded::~DataExpanded()
180  }  }
181    
182  DataAbstract*  DataAbstract*
183  DataExpanded::getSlice(const DataArrayView::RegionType& region) const  DataExpanded::deepCopy()
184    {
185      return new DataExpanded(*this);
186    }
187    
188    
189    DataAbstract*
190    DataExpanded::getSlice(const DataTypes::RegionType& region) const
191  {  {
192    return new DataExpanded(*this,region);    return new DataExpanded(*this,region);
193  }  }
194    
195  void  void
196  DataExpanded::setSlice(const DataAbstract* value,  DataExpanded::setSlice(const DataAbstract* value,
197                         const DataArrayView::RegionType& region)                         const DataTypes::RegionType& region)
198  {  {
199    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
200    if (tempDataExp==0) {    if (tempDataExp==0) {
# Line 175  DataExpanded::setSlice(const DataAbstrac Line 202  DataExpanded::setSlice(const DataAbstrac
202    }    }
203    //    //
204    // get shape of slice    // get shape of slice
205    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));    DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
206    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
207    //    //
208    // check shape    // check shape
209    if (getPointDataView().getRank()!=region.size()) {    if (getRank()!=region.size()) {
210      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
211    }    }
212    if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {    if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
213      throw DataException (value->getPointDataView().createShapeErrorMessage(      throw DataException (DataTypes::createShapeErrorMessage(
214          "Error - Couldn't copy slice due to shape mismatch.",shape));          "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
215    }    }
216    //    //
217    // copy the data from the slice into this object    // copy the data from the slice into this object
218    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataTypes::ValueType::size_type numRows=m_data.getNumRows();
219    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataTypes::ValueType::size_type numCols=m_data.getNumCols();
220    int i, j;    int i, j;
221      ValueType& vec=getVector();
222      const ShapeType& mshape=getShape();
223      const ValueType& tVec=tempDataExp->getVector();
224      const ShapeType& tShape=tempDataExp->getShape();
225    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
226    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
227      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
228        getPointDataView().copySliceFrom(getPointOffset(i,j),  /*      getPointDataView().copySliceFrom(getPointOffset(i,j),
229                                         tempDataExp->getPointDataView(),                                         tempDataExp->getPointDataView(),
230                                         tempDataExp->getPointOffset(i,j),                                         tempDataExp->getPointOffset(i,j),
231                                           region_loop_range);*/
232            DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
233                                           tVec,
234                           tShape,
235                                           tempDataExp->getPointOffset(i,j),
236                                         region_loop_range);                                         region_loop_range);
237    
238      }      }
239    }    }
240  }  }
241    
242  void  void
243  DataExpanded::copy(const DataArrayView& value)  DataExpanded::copy(const DataConstant& value)
244  {  {
245      EsysAssert((checkShape(getShape(), value.getShape())),
246                     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
247    
248    //    //
249    // copy a single value to every data point in this object    // copy a single value to every data point in this object
250    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
# Line 217  DataExpanded::copy(const DataArrayView& Line 257  DataExpanded::copy(const DataArrayView&
257        // DOASSERT is on which of course will play        // DOASSERT is on which of course will play
258        // havoc with the omp threads. Run single threaded        // havoc with the omp threads. Run single threaded
259        // if using DOASSERT.        // if using DOASSERT.
260        getPointDataView().copy(getPointOffset(i,j),value);        //getPointDataView().copy(getPointOffset(i,j),value);
261          DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(), value.getVector(), 0);
262      }      }
263    }    }
264  }  }
265    
266    
267    // void
268    // DataExpanded::copy(const DataArrayView& value)
269    // {
270    //   //
271    //   // copy a single value to every data point in this object
272    //   int nRows=m_data.getNumRows();
273    //   int nCols=m_data.getNumCols();
274    //   int i,j;
275    //   #pragma omp parallel for private(i,j) schedule(static)
276    //   for (i=0;i<nRows;i++) {
277    //     for (j=0;j<nCols;j++) {
278    //       // NOTE: An exception may be thown from this call if
279    //       // DOASSERT is on which of course will play
280    //       // havoc with the omp threads. Run single threaded
281    //       // if using DOASSERT.
282    //       getPointDataView().copy(getPointOffset(i,j),value);
283    //     }
284    //   }
285    // }
286    
287  void  void
288  DataExpanded::copy(const boost::python::numeric::array& value)  DataExpanded::copy(const boost::python::numeric::array& value)
289  {  {
290    
291    // extract the shape of the numarray    // extract the shape of the numarray
292    DataArrayView::ShapeType tempShape;    DataTypes::ShapeType tempShape;
293    for (int i=0; i < value.getrank(); i++) {    for (int i=0; i < value.getrank(); i++) {
294      tempShape.push_back(extract<int>(value.getshape()[i]));      tempShape.push_back(extract<int>(value.getshape()[i]));
295    }    }
296    
297    // get the space for the data vector    // get the space for the data vector
298    int len = DataArrayView::noValues(tempShape);  //   int len = DataTypes::noValues(tempShape);
299    DataVector temp_data(len, 0.0, len);  //   DataVector temp_data(len, 0.0, len);
300    DataArrayView temp_dataView(temp_data, tempShape);  //   DataArrayView temp_dataView(temp_data, tempShape);
301    temp_dataView.copy(value);  //   temp_dataView.copy(value);
302    
303    //    //
304    // check the input shape matches this shape    // check the input shape matches this shape
305    if (!getPointDataView().checkShape(temp_dataView.getShape())) {    if (!DataTypes::checkShape(getShape(),tempShape)) {
306      throw DataException(getPointDataView().createShapeErrorMessage(      throw DataException(DataTypes::createShapeErrorMessage(
307                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
308                          temp_dataView.getShape()));                          tempShape,getShape()));
309    }    }
310    //    //
311    // now copy over the data    // now copy over the data
312    copy(temp_dataView);    //copy(temp_dataView);
313      getVector().copyFromNumArray(value);
314  }  }
315    
316    
317  void  void
318  DataExpanded::initialise(const DataArrayView::ShapeType& shape,  DataExpanded::initialise(int noSamples,
                          int noSamples,  
319                           int noDataPointsPerSample)                           int noDataPointsPerSample)
320  {  {
321      if (noSamples==0)     //retain the default empty object
322      {
323         return;
324      }
325    //    //
326    // resize data array to the required size    // resize data array to the required size
327    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);  
328  }  }
329    
330  string  string
# Line 269  DataExpanded::toString() const Line 332  DataExpanded::toString() const
332  {  {
333    stringstream temp;    stringstream temp;
334    FunctionSpace fs=getFunctionSpace();    FunctionSpace fs=getFunctionSpace();
335    //  
336    // create a temporary view as the offset will be changed    int offset=0;
   DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());  
337    for (int i=0;i<m_data.getNumRows();i++) {    for (int i=0;i<m_data.getNumRows();i++) {
338      for (int j=0;j<m_data.getNumCols();j++) {      for (int j=0;j<m_data.getNumCols();j++) {
339        tempView.setOffset(getPointOffset(i,j));        offset=getPointOffset(i,j);
340        stringstream suffix;        stringstream suffix;
341        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
342        temp << tempView.toString(suffix.str());        temp << DataTypes::pointToString(getVector(),getShape(),offset,suffix.str());
343        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {        if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
344          temp << endl;          temp << endl;
345        }        }
346      }      }
347    }    }
348      string result=temp.str();
349      if (result.empty())
350      {
351        return "(data contains no samples)\n";
352      }
353    return temp.str();    return temp.str();
354  }  }
355    
356  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
357  DataExpanded::getPointOffset(int sampleNo,  DataExpanded::getPointOffset(int sampleNo,
358                               int dataPointNo) const                               int dataPointNo) const
359  {  {
360    return m_data.index(sampleNo,dataPointNo);    return m_data.index(sampleNo,dataPointNo);
361  }  }
362    
363  DataArrayView  DataTypes::ValueType::size_type
364  DataExpanded::getDataPoint(int sampleNo,  DataExpanded::getPointOffset(int sampleNo,
365                             int dataPointNo)                               int dataPointNo)
366  {  {
367    DataArrayView temp(m_data.getData(),getPointDataView().getShape(),getPointOffset(sampleNo,dataPointNo));    return m_data.index(sampleNo,dataPointNo);
   return temp;  
368  }  }
369    
370  DataArrayView::ValueType::size_type  DataTypes::ValueType::size_type
371  DataExpanded::getLength() const  DataExpanded::getLength() const
372  {  {
373    return m_data.size();    return m_data.size();
374  }  }
375    
 int  
 DataExpanded::archiveData(ofstream& archiveFile,  
                           const DataArrayView::ValueType::size_type noValues) const  
 {  
   return(m_data.archiveData(archiveFile, noValues));  
 }  
376    
 int  
 DataExpanded::extractData(ifstream& archiveFile,  
                           const DataArrayView::ValueType::size_type noValues)  
 {  
   return(m_data.extractData(archiveFile, noValues));  
 }  
377    
378  void  void
379  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 381  DataExpanded::copyToDataPoint(const int
381    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
382    int numSamples = getNumSamples();    int numSamples = getNumSamples();
383    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
384    int dataPointRank = getPointDataView().getRank();    int dataPointRank = getRank();
385    ShapeType dataPointShape = getPointDataView().getShape();    ShapeType dataPointShape = getShape();
386    if (numSamples*numDataPointsPerSample > 0) {    if (numSamples*numDataPointsPerSample > 0) {
387       //TODO: global error handling       //TODO: global error handling
388       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {       if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
# Line 337  DataExpanded::copyToDataPoint(const int Line 391  DataExpanded::copyToDataPoint(const int
391       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
392             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
393       }       }
394       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
395         ValueType& vec=getVector();
396       if (dataPointRank==0) {       if (dataPointRank==0) {
397           dataPointView()=value;           vec[offset]=value;
398       } else if (dataPointRank==1) {       } else if (dataPointRank==1) {
399          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
400              dataPointView(i)=value;              vec[offset+i]=value;
401          }          }
402       } else if (dataPointRank==2) {       } else if (dataPointRank==2) {
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                dataPointView(i,j)=value;                vec[offset+getRelIndex(dataPointShape,i,j)]=value;
406             }             }
407          }          }
408       } else if (dataPointRank==3) {       } else if (dataPointRank==3) {
409          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
410             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
411                for (int k=0; k<dataPointShape[2]; k++) {                for (int k=0; k<dataPointShape[2]; k++) {
412                   dataPointView(i,j,k)=value;                   vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
413                }                }
414             }             }
415          }          }
# Line 363  DataExpanded::copyToDataPoint(const int Line 418  DataExpanded::copyToDataPoint(const int
418             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
419               for (int k=0; k<dataPointShape[2]; k++) {               for (int k=0; k<dataPointShape[2]; k++) {
420                 for (int l=0; l<dataPointShape[3]; l++) {                 for (int l=0; l<dataPointShape[3]; l++) {
421                    dataPointView(i,j,k,l)=value;                    vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
422                 }                 }
423               }               }
424             }             }
# Line 377  DataExpanded::copyToDataPoint(const int Line 432  DataExpanded::copyToDataPoint(const int
432    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
433    int numSamples = getNumSamples();    int numSamples = getNumSamples();
434    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
435    int dataPointRank = getPointDataView().getRank();    int dataPointRank = getRank();
436    ShapeType dataPointShape = getPointDataView().getShape();    const ShapeType& shape = getShape();
437    //    //
438    // check rank:    // check rank:
439    if (value.getrank()!=dataPointRank)    if (value.getrank()!=dataPointRank)
# Line 391  DataExpanded::copyToDataPoint(const int Line 446  DataExpanded::copyToDataPoint(const int
446       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {       if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
447             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");             throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
448       }       }
449       DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);       ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
450         ValueType& vec=getVector();
451       if (dataPointRank==0) {       if (dataPointRank==0) {
452           dataPointView()=extract<double>(value[0]);           vec[offset]=extract<double>(value[0]);
453       } else if (dataPointRank==1) {       } else if (dataPointRank==1) {
454          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<shape[0]; i++) {
455              dataPointView(i)=extract<double>(value[i]);              vec[offset+i]=extract<double>(value[i]);
456          }          }
457       } else if (dataPointRank==2) {       } else if (dataPointRank==2) {
458          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<shape[0]; i++) {
459             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<shape[1]; j++) {
460                dataPointView(i,j)=extract<double>(value[i][j]);                vec[offset+getRelIndex(shape,i,j)]=extract<double>(value[i][j]);
461             }             }
462          }          }
463       } else if (dataPointRank==3) {       } else if (dataPointRank==3) {
464          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<shape[0]; i++) {
465             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<shape[1]; j++) {
466                for (int k=0; k<dataPointShape[2]; k++) {                for (int k=0; k<shape[2]; k++) {
467                   dataPointView(i,j,k)=extract<double>(value[i][j][k]);                   vec[offset+getRelIndex(shape,i,j,k)]=extract<double>(value[i][j][k]);
468                }                }
469             }             }
470          }          }
471       } else if (dataPointRank==4) {       } else if (dataPointRank==4) {
472           for (int i=0; i<dataPointShape[0]; i++) {           for (int i=0; i<shape[0]; i++) {
473             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<shape[1]; j++) {
474               for (int k=0; k<dataPointShape[2]; k++) {               for (int k=0; k<shape[2]; k++) {
475                 for (int l=0; l<dataPointShape[3]; l++) {                 for (int l=0; l<shape[3]; l++) {
476                    dataPointView(i,j,k,l)=extract<double>(value[i][j][k][l]);                    vec[offset+getRelIndex(shape,i,j,k,l)]=extract<double>(value[i][j][k][l]);
477                 }                 }
478               }               }
479             }             }
# Line 431  DataExpanded::copyAll(const boost::pytho Line 487  DataExpanded::copyAll(const boost::pytho
487    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
488    int numSamples = getNumSamples();    int numSamples = getNumSamples();
489    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
490    int dataPointRank = getPointDataView().getRank();    int dataPointRank = getRank();
491    ShapeType dataPointShape = getPointDataView().getShape();    const ShapeType& dataPointShape = getShape();
492    //    //
493    // check rank:    // check rank:
494    if (value.getrank()!=dataPointRank+1)    if (value.getrank()!=dataPointRank+1)
# Line 440  DataExpanded::copyAll(const boost::pytho Line 496  DataExpanded::copyAll(const boost::pytho
496    if (value.getshape()[0]!=numSamples*numDataPointsPerSample)    if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
497         throw DataException("leading dimension of numarray is too small");         throw DataException("leading dimension of numarray is too small");
498    //    //
499      ValueType& vec=getVector();
500    int dataPoint = 0;    int dataPoint = 0;
501    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
502      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
503        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);        ValueType::size_type offset=getPointOffset(sampleNo, dataPointNo);
504        if (dataPointRank==0) {        if (dataPointRank==0) {
505           dataPointView()=extract<double>(value[dataPoint]);           vec[offset]=extract<double>(value[dataPoint]);
506        } else if (dataPointRank==1) {        } else if (dataPointRank==1) {
507           for (int i=0; i<dataPointShape[0]; i++) {           for (int i=0; i<dataPointShape[0]; i++) {
508              dataPointView(i)=extract<double>(value[dataPoint][i]);              vec[offset+i]=extract<double>(value[dataPoint][i]);
509           }           }
510        } else if (dataPointRank==2) {        } else if (dataPointRank==2) {
511           for (int i=0; i<dataPointShape[0]; i++) {           for (int i=0; i<dataPointShape[0]; i++) {
512             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
513               dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);           vec[offset+getRelIndex(dataPointShape,i,j)]=extract<double>(value[dataPoint][i][j]);
514             }             }
515           }           }
516         } else if (dataPointRank==3) {         } else if (dataPointRank==3) {
517           for (int i=0; i<dataPointShape[0]; i++) {           for (int i=0; i<dataPointShape[0]; i++) {
518             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
519               for (int k=0; k<dataPointShape[2]; k++) {               for (int k=0; k<dataPointShape[2]; k++) {
520                   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]);
521               }               }
522             }             }
523           }           }
# Line 469  DataExpanded::copyAll(const boost::pytho Line 526  DataExpanded::copyAll(const boost::pytho
526             for (int j=0; j<dataPointShape[1]; j++) {             for (int j=0; j<dataPointShape[1]; j++) {
527               for (int k=0; k<dataPointShape[2]; k++) {               for (int k=0; k<dataPointShape[2]; k++) {
528                 for (int l=0; l<dataPointShape[3]; l++) {                 for (int l=0; l<dataPointShape[3]; l++) {
529                   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]);
530                 }                 }
531               }               }
532             }             }
# Line 489  DataExpanded::symmetric(DataAbstract* ev Line 546  DataExpanded::symmetric(DataAbstract* ev
546    if (temp_ev==0) {    if (temp_ev==0) {
547      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).");
548    }    }
549    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
550    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
551      ValueType& evVec=temp_ev->getVector();
552      const ShapeType& evShape=temp_ev->getShape();
553    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
554    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
555      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
556           DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
557                                      evView,ev->getPointOffset(sampleNo,dataPointNo));                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
558      }      }
559    }    }
560  }  }
# Line 509  DataExpanded::nonsymmetric(DataAbstract* Line 568  DataExpanded::nonsymmetric(DataAbstract*
568    if (temp_ev==0) {    if (temp_ev==0) {
569      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).");
570    }    }
571    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
572    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
573      ValueType& evVec=temp_ev->getVector();
574      const ShapeType& evShape=temp_ev->getShape();
575    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
576    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
577      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
578           DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
579                                      evView,ev->getPointOffset(sampleNo,dataPointNo));                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
580      }      }
581    }    }
582  }  }
# Line 529  DataExpanded::trace(DataAbstract* ev, in Line 590  DataExpanded::trace(DataAbstract* ev, in
590    if (temp_ev==0) {    if (temp_ev==0) {
591      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).");
592    }    }
593    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
594    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
595      ValueType& evVec=temp_ev->getVector();
596      const ShapeType& evShape=temp_ev->getShape();
597    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
598    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
599      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
600           DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
601                                      evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
602      }      }
603    }    }
604  }  }
# Line 550  DataExpanded::transpose(DataAbstract* ev Line 613  DataExpanded::transpose(DataAbstract* ev
613    if (temp_ev==0) {    if (temp_ev==0) {
614      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).");
615    }    }
616    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
617    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
618      ValueType& evVec=temp_ev->getVector();
619      const ShapeType& evShape=temp_ev->getShape();
620    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
621    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
622      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
623           DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
624                                      evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
625      }      }
626    }    }
627  }  }
# Line 571  DataExpanded::swapaxes(DataAbstract* ev, Line 636  DataExpanded::swapaxes(DataAbstract* ev,
636    if (temp_ev==0) {    if (temp_ev==0) {
637      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).");
638    }    }
639    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
640    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
641      ValueType& evVec=temp_ev->getVector();
642      const ShapeType& evShape=temp_ev->getShape();
643    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
644    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
645      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
646           DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
647                                      evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
648      }      }
649    }    }
650  }  }
# Line 591  DataExpanded::eigenvalues(DataAbstract* Line 658  DataExpanded::eigenvalues(DataAbstract*
658    if (temp_ev==0) {    if (temp_ev==0) {
659      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).");
660    }    }
661    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
662    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
663      ValueType& evVec=temp_ev->getVector();
664      const ShapeType& evShape=temp_ev->getShape();
665    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
666    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
667      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
668           DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
669                                      evView,ev->getPointOffset(sampleNo,dataPointNo));                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
670      }      }
671    }    }
672  }  }
# Line 615  DataExpanded::eigenvalues_and_eigenvecto Line 684  DataExpanded::eigenvalues_and_eigenvecto
684    if (temp_V==0) {    if (temp_V==0) {
685      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).");
686    }    }
687    DataArrayView& thisView=getPointDataView();    ValueType& vec=getVector();
688    DataArrayView& evView=ev->getPointDataView();    const ShapeType& shape=getShape();
689    DataArrayView& VView=V->getPointDataView();    ValueType& evVec=temp_ev->getVector();
690      const ShapeType& evShape=temp_ev->getShape();
691      ValueType& VVec=temp_V->getVector();
692      const ShapeType& VShape=temp_V->getShape();
693    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
694    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
695      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
696           DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),           DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
697                                                       evView,ev->getPointOffset(sampleNo,dataPointNo),                                      evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
698                                                       VView,V->getPointOffset(sampleNo,dataPointNo),                                      VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
                                                      tol);  
699      }      }
700    }    }
701  }  }
702    
703  void  void
704  DataExpanded::setToZero(){  DataExpanded::setToZero(){
705    // TODO: Surely there is a more efficient way to do this????
706    // Why is there no memset here? Parallel issues?
707    int numSamples = getNumSamples();    int numSamples = getNumSamples();
708    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
709    DataArrayView& thisView=getPointDataView();    DataTypes::ValueType::size_type n = getNoValues();
   DataArrayView::ValueType::size_type n = thisView.noValues();  
710    double* p;    double* p;
711    int  sampleNo,dataPointNo, i;    int  sampleNo,dataPointNo, i;
712    #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)    #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
713    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
714      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
715          p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);          p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
716          for (int i=0; i<n ;++i) p[i]=0.;          for (i=0; i<n ;++i) p[i]=0.;
717      }      }
718    }    }
719  }  }
720    
721    /* Append MPI rank to file name if multiple MPI processes */
722    char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
723      /* Make plenty of room for the mpi_rank number and terminating '\0' */
724      char *newFileName = (char *)malloc(strlen(fileName)+20);
725      strncpy(newFileName, fileName, strlen(fileName)+1);
726      if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
727      return(newFileName);
728    }
729    
730  void  void
731  DataExpanded::dump(const std::string fileName) const  DataExpanded::dump(const std::string fileName) const
732  {  {
    #ifdef PASO_MPI  
    throw DataException("Error - DataExpanded:: dump is not implemented for MPI yet.");  
    #endif  
733     #ifdef USE_NETCDF     #ifdef USE_NETCDF
734     const int ldims=2+DataArrayView::maxRank;     const int ldims=2+DataTypes::maxRank;
735     const NcDim* ncdims[ldims];     const NcDim* ncdims[ldims];
736     NcVar *var, *ids;     NcVar *var, *ids;
737     int rank = getPointDataView().getRank();     int rank = getRank();
738     int type=  getFunctionSpace().getTypeCode();     int type=  getFunctionSpace().getTypeCode();
739     int ndims =0;     int ndims =0;
740     long dims[ldims];     long dims[ldims];
741     const double* d_ptr=&(m_data[0]);     const double* d_ptr=&(m_data[0]);
742     DataArrayView::ShapeType shape = getPointDataView().getShape();     const DataTypes::ShapeType& shape = getShape();
743       int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
744       int mpi_num=getFunctionSpace().getDomain()->getMPISize();
745    #ifdef PASO_MPI
746       MPI_Status status;
747    #endif
748    
749    #ifdef PASO_MPI
750       /* Serialize NetCDF I/O */
751       if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
752    #endif
753    
754     // netCDF error handler     // netCDF error handler
755     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
756     // Create the file.     // Create the file.
757     NcFile dataFile(fileName.c_str(), NcFile::Replace);     char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
758       NcFile dataFile(newFileName, NcFile::Replace);
759     // check if writing was successful     // check if writing was successful
760     if (!dataFile.is_valid())     if (!dataFile.is_valid())
761          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 769  DataExpanded::dump(const std::string fil
769     if ( rank >0 ) {     if ( rank >0 ) {
770         dims[0]=shape[0];         dims[0]=shape[0];
771         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )         if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
772              throw DataException("Error - DataExpanded:: appending ncdimsion 0 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
773     }     }
774     if ( rank >1 ) {     if ( rank >1 ) {
775         dims[1]=shape[1];         dims[1]=shape[1];
776         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )         if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
777              throw DataException("Error - DataExpanded:: appending ncdimsion 1 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
778     }     }
779     if ( rank >2 ) {     if ( rank >2 ) {
780         dims[2]=shape[2];         dims[2]=shape[2];
781         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )         if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
782              throw DataException("Error - DataExpanded:: appending ncdimsion 2 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
783     }     }
784     if ( rank >3 ) {     if ( rank >3 ) {
785         dims[3]=shape[3];         dims[3]=shape[3];
786         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )         if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
787              throw DataException("Error - DataExpanded:: appending ncdimsion 3 to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
788     }     }
789     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
790     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 793  DataExpanded::dump(const std::string fil
793     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
794              throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");              throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
795    
796     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )     if (getFunctionSpace().getNumSamples()>0)
797       {
798    
799         if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
800          throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");          throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
801     const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();       const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
802     if (! (ids->put(ids_p,dims[rank+1])) )       if (! (ids->put(ids_p,dims[rank+1])) )
803          throw DataException("Error - DataExpanded:: copy reference id  to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy reference id  to netCDF buffer failed.");
804         if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
    if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )  
805          throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");          throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
806     if (! (var->put(d_ptr,dims)) )       if (! (var->put(d_ptr,dims)) )
807          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
808       }
809    #ifdef PASO_MPI
810       if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
811    #endif
812     #else     #else
813     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.");
814     #endif     #endif
815  }  }
816    
817  void  void  
818  DataExpanded::setTaggedValue(int tagKey,  DataExpanded::setTaggedValue(int tagKey,
819                               const DataArrayView& value)             const DataTypes::ShapeType& pointshape,
820                   const DataTypes::ValueType& value,
821               int dataOffset)
822  {  {
823    int numSamples = getNumSamples();    int numSamples = getNumSamples();
824    int numDataPointsPerSample = getNumDPPSample();    int numDataPointsPerSample = getNumDPPSample();
825    int sampleNo,dataPointNo, i;    int sampleNo,dataPointNo, i;
826    DataArrayView& thisView=getPointDataView();    DataTypes::ValueType::size_type n = getNoValues();
827    DataArrayView::ValueType::size_type n = thisView.noValues();    double* p;
828    double* p,*in=&(value.getData()[0]);    const double* in=&value[0+dataOffset];
829        
830    if (value.noValues() != n) {    if (value.size() != n) {
831      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.");
832    }    }
833    
# Line 740  DataExpanded::setTaggedValue(int tagKey, Line 836  DataExpanded::setTaggedValue(int tagKey,
836      if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {      if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
837          for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {          for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
838              p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);              p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
839              for (int i=0; i<n ;++i) p[i]=in[i];              for (i=0; i<n ;++i) p[i]=in[i];
840          }          }
841      }      }
842    }    }
843  }  }
844    
845    
846    void
847    DataExpanded::reorderByReferenceIDs(int *reference_ids)
848    {
849      int numSamples = getNumSamples();
850      DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
851      int sampleNo, sampleNo2,i;
852      double* p,*p2;
853      register double rtmp;
854      FunctionSpace fs=getFunctionSpace();
855    
856      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
857         const int id_in=reference_ids[sampleNo];
858         const int id=fs.getReferenceIDOfSample(sampleNo);
859         if (id!=id_in) {
860             bool matched=false;
861             for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
862                  if (id == reference_ids[sampleNo2]) {
863                     p=&(m_data[getPointOffset(sampleNo,0)]);
864                     p2=&(m_data[getPointOffset(sampleNo2,0)]);
865                     for (i=0; i<n ;i++) {
866                             rtmp=p[i];
867                             p[i]=p2[i];
868                             p2[i]=rtmp;
869                     }
870                     reference_ids[sampleNo]=id;
871                     reference_ids[sampleNo2]=id_in;
872                     matched=true;
873                     break;
874                  }
875             }
876             if (! matched) {
877                throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
878             }
879         }
880       }
881    }
882    
883    DataTypes::ValueType&
884    DataExpanded::getVector()
885    {
886        return m_data.getData();
887    }
888    
889    const DataTypes::ValueType&
890    DataExpanded::getVector() const
891    {
892        return m_data.getData();
893    }
894    
895  }  // end of namespace  }  // end of namespace
Legend:
Removed from v.1387  
changed lines
  Added in v.2005
  ViewVC Help
Powered by ViewVC 1.1.26