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

trunk/esys2/escript/src/Data/DataExpanded.cpp revision 97 by jgs, Tue Dec 14 05:39:33 2004 UTC trunk/escript/src/DataExpanded.cpp revision 580 by gross, Wed Mar 8 05:45:51 2006 UTC
# Line 13  Line 13 
13   ******************************************************************************   ******************************************************************************
14  */  */
15    
16  #include "escript/Data/DataException.h"  #include "DataExpanded.h"
17  #include "escript/Data/DataExpanded.h"  #include "DataException.h"
18  #include "escript/Data/DataConstant.h"  #include "DataConstant.h"
19  #include "escript/Data/DataTagged.h"  #include "DataTagged.h"
 #include "escript/Data/DataArrayView.h"  
20    
21  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
22    
 #include <iostream>  
   
23  using namespace std;  using namespace std;
24  using namespace boost::python;  using namespace boost::python;
25  using namespace boost;  using namespace boost;
# Line 35  DataExpanded::DataExpanded(const boost:: Line 32  DataExpanded::DataExpanded(const boost::
32  {  {
33    DataArrayView::ShapeType tempShape;    DataArrayView::ShapeType tempShape;
34    //    //
35    // extract the shape from the python array    // extract the shape of the python numarray
36    for (int i=0; i<value.getrank(); ++i) {    for (int i=0; i<value.getrank(); i++) {
     //cout << extract<int>(value.getshape()[i]) << endl;  
37      tempShape.push_back(extract<int>(value.getshape()[i]));      tempShape.push_back(extract<int>(value.getshape()[i]));
38    }    }
39      //
40      // initialise the data array for this object
41    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
42    //    //
43    // copy the given value to every data point    // copy the given value to every data point
# Line 49  DataExpanded::DataExpanded(const boost:: Line 47  DataExpanded::DataExpanded(const boost::
47  DataExpanded::DataExpanded(const DataExpanded& other)  DataExpanded::DataExpanded(const DataExpanded& other)
48    : DataAbstract(other.getFunctionSpace()),    : DataAbstract(other.getFunctionSpace()),
49    m_data(other.m_data)    m_data(other.m_data)
50  {      {
51    //    //
52    // create the view for the data    // create the view for the data
53    DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());    DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());
# Line 58  DataExpanded::DataExpanded(const DataExp Line 56  DataExpanded::DataExpanded(const DataExp
56    
57  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
58    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace())
59  {      {
60      //
61      // initialise the data array for this object
62    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
63    //    //
64    // DataConstant only has one value    // DataConstant only has one value, copy this to every data point
65    copy(other.getPointDataView());    copy(other.getPointDataView());
66  }  }
67    
68  DataExpanded::DataExpanded(const DataTagged& other)  DataExpanded::DataExpanded(const DataTagged& other)
69    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace())
70  {      {
71      //
72      // initialise the data array for this object
73    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
74      //
75      // for each data point in this object, extract and copy the corresponding data
76      // value from the given DataTagged object
77    int i,j;    int i,j;
78    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
79    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
80  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
81    for (i=0;i<numRows;++i) {    for (i=0;i<numRows;i++) {
82      for (j=0;j<numCols;++j) {      for (j=0;j<numCols;j++) {
83        try {        try {
84          getPointDataView().copy(getPointOffset(i,j), other.getPointDataView(), other.getPointOffset(i,j));          getPointDataView().copy(getPointOffset(i,j),
85                                    other.getPointDataView(),
86                                    other.getPointOffset(i,j));
87        }        }
88        catch (std::exception& e) {        catch (std::exception& e) {
89          cout << e.what() << endl;          cout << e.what() << endl;
# Line 97  DataExpanded::DataExpanded(const DataExp Line 104  DataExpanded::DataExpanded(const DataExp
104    initialise(shape,other.getNumSamples(),other.getNumDPPSample());    initialise(shape,other.getNumSamples(),other.getNumDPPSample());
105    //    //
106    // copy the data    // copy the data
107      DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
108    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
109    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
110    int i,j;    int i,j;
111  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
112    for (i=0;i<numRows;++i) {    for (i=0;i<numRows;i++) {
113      for (j=0;j<numCols;++j) {      for (j=0;j<numCols;j++) {
114        try {        try {
115          getPointDataView().copySlice(getPointOffset(i,j),other.getPointDataView(),other.getPointOffset(i,j),region);          getPointDataView().copySlice(getPointOffset(i,j),
116                                         other.getPointDataView(),
117                                         other.getPointOffset(i,j),
118                                         region_loop_range);
119        }        }
120        catch (std::exception& e) {        catch (std::exception& e) {
121          cout << e.what() << endl;          cout << e.what() << endl;
# Line 117  DataExpanded::DataExpanded(const DataArr Line 128  DataExpanded::DataExpanded(const DataArr
128                             const FunctionSpace& what)                             const FunctionSpace& what)
129    : DataAbstract(what)    : DataAbstract(what)
130  {  {
131      //
132      // get the shape of the given data value
133    DataArrayView::ShapeType tempShape=value.getShape();    DataArrayView::ShapeType tempShape=value.getShape();
134      //
135      // initialise this Data object to the shape of the given data value
136    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
137      //
138      // copy the given value to every data point
139    copy(value);    copy(value);
140  }  }
141    
142    DataExpanded::DataExpanded(const FunctionSpace& what,
143                               const DataArrayView::ShapeType &shape,
144                               const DataArrayView::ValueType &data)
145      : DataAbstract(what)
146    {
147      //
148      // create the view of the data
149      initialise(shape,what.getNumSamples(),what.getNumDPPSample());
150      //
151      // copy the data in the correct format
152      m_data.getData()=data;
153    }
154    
155  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
156  {  {
157  }  }
# Line 129  DataExpanded::~DataExpanded() Line 159  DataExpanded::~DataExpanded()
159  void  void
160  DataExpanded::reshapeDataPoint(const DataArrayView::ShapeType& shape)  DataExpanded::reshapeDataPoint(const DataArrayView::ShapeType& shape)
161  {  {
   //  
   // reshape a rank zero data point  
162    if (getPointDataView().getRank()!=0) {    if (getPointDataView().getRank()!=0) {
163      stringstream temp;      stringstream temp;
164      temp << "Error - Can only reshape Data with data points of rank 0. "      temp << "Error - Can only reshape Data with data points of rank 0. "
# Line 138  DataExpanded::reshapeDataPoint(const Dat Line 166  DataExpanded::reshapeDataPoint(const Dat
166           << getPointDataView().getRank();           << getPointDataView().getRank();
167      throw DataException(temp.str());      throw DataException(temp.str());
168    }    }
169    DataBlocks2D newData(getNumSamples(),getNumDPPSample(), DataArrayView::noValues(shape));    //
170      // create the new DataBlocks2D data array, and a corresponding DataArrayView
171      DataBlocks2D newData(getNumSamples(),getNumDPPSample(),DataArrayView::noValues(shape));
172    DataArrayView newView(newData.getData(),shape);    DataArrayView newView(newData.getData(),shape);
173    //    //
174    // Copy the original data to every value for the new shape    // Copy the original data to every value for the new shape
175    int i,j;    int i,j;
176    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
177    int nCols=m_data.getNumCols();    int nCols=m_data.getNumCols();
178  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
179    for (i=0;i<nRows;++i) {    for (i=0;i<nRows;i++) {
180      for (j=0;j<nCols;++j) {      for (j=0;j<nCols;j++) {
       //  
       // Copy the data into the specified offset  
181        // NOTE: An exception may be thown from this call if        // NOTE: An exception may be thown from this call if
182        // DOASSERT is on which of course will play        // DOASSERT is on which of course will play
183        // havoc with the omp threads. Run single threaded        // havoc with the omp threads. Run single threaded
184        // if using DOASSERT.        // if using DOASSERT.
185        newView.copy(newData.index(i,j),m_data.getData()[m_data.index(i,j)]);        newView.copy(newData.index(i,j),m_data(i,j));
186      }      }
187    }    }
188      // swap the new data array for the original
189    m_data.Swap(newData);    m_data.Swap(newData);
190      // set the corresponding DataArrayView
191    DataArrayView temp(m_data.getData(),shape);    DataArrayView temp(m_data.getData(),shape);
192    setPointDataView(temp);    setPointDataView(temp);
193  }  }
# Line 176  DataExpanded::setSlice(const DataAbstrac Line 206  DataExpanded::setSlice(const DataAbstrac
206    if (tempDataExp==0) {    if (tempDataExp==0) {
207      throw DataException("Programming error - casting to DataExpanded.");      throw DataException("Programming error - casting to DataExpanded.");
208    }    }
209      //
210      // get shape of slice
211      DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
212      DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
213      //
214      // check shape
215    if (getPointDataView().getRank()!=region.size()) {    if (getPointDataView().getRank()!=region.size()) {
216      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
217    }    }
218    if (!value->getPointDataView().checkShape(DataArrayView::getResultSliceShape(region))) {    if (tempDataExp->getPointDataView().getRank()>0 and !value->getPointDataView().checkShape(shape)) {
219      throw DataException (value->getPointDataView().createShapeErrorMessage(      throw DataException (value->getPointDataView().createShapeErrorMessage(
220          "Error - Couldn't copy slice due to shape mismatch.",DataArrayView::getResultSliceShape(region)));          "Error - Couldn't copy slice due to shape mismatch.",shape));
221    }    }
222    //    //
223    // copy the data    // copy the data from the slice into this object
224    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
225    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
226    int i, j;    int i, j;
227  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
228    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
229      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
230        getPointDataView().copySliceFrom(getPointOffset(i,j),tempDataExp->getPointDataView(),tempDataExp->getPointOffset(i,j),region);        getPointDataView().copySliceFrom(getPointOffset(i,j),
231                                           tempDataExp->getPointDataView(),
232                                           tempDataExp->getPointOffset(i,j),
233                                           region_loop_range);
234      }      }
235    }    }
236  }  }
# Line 200  void Line 239  void
239  DataExpanded::copy(const DataArrayView& value)  DataExpanded::copy(const DataArrayView& value)
240  {  {
241    //    //
242    // Copy a single value to every data point    // copy a single value to every data point in this object
243    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
244    int nCols=m_data.getNumCols();    int nCols=m_data.getNumCols();
245    int i, j;    int i,j;
246  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
247    for (i=0;i<nRows;++i) {    for (i=0;i<nRows;i++) {
248      for (j=0;j<nCols;++j) {      for (j=0;j<nCols;j++) {
       //  
       // Copy the data into the specified offset  
249        // NOTE: An exception may be thown from this call if        // NOTE: An exception may be thown from this call if
250        // DOASSERT is on which of course will play        // DOASSERT is on which of course will play
251        // havoc with the omp threads. Run single threaded        // havoc with the omp threads. Run single threaded
# Line 222  void Line 259  void
259  DataExpanded::copy(const boost::python::numeric::array& value)  DataExpanded::copy(const boost::python::numeric::array& value)
260  {  {
261    //    //
262    // first convert the numarray into a DataArrayView format    // first convert the numarray into a DataArray object
263    DataArray temp(value);    DataArray temp(value);
264    //    //
265    // check the input shape matches this shape, this will throw an exception    // check the input shape matches this shape
266    if (!getPointDataView().checkShape(temp.getView().getShape())) {    if (!getPointDataView().checkShape(temp.getView().getShape())) {
267      throw DataException(getPointDataView().createShapeErrorMessage(      throw DataException(getPointDataView().createShapeErrorMessage(
268                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
269                          temp.getView().getShape()));                          temp.getView().getShape()));
270    }    }
271    //    //
272    // now copy over the entire data structure    // now copy over the data
273    copy(temp.getView());    copy(temp.getView());
274  }  }
275    
# Line 242  DataExpanded::initialise(const DataArray Line 279  DataExpanded::initialise(const DataArray
279                           int noDataPointsPerSample)                           int noDataPointsPerSample)
280  {  {
281    //    //
282    // resize data to the required size    // resize data array to the required size
283    m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));    m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));
284    //    //
285    // create a point data viewer of the data    // create the data view of the data array
286    DataArrayView temp(m_data.getData(),shape);    DataArrayView temp(m_data.getData(),shape);
287    setPointDataView(temp);    setPointDataView(temp);
288  }  }
# Line 257  DataExpanded::toString() const Line 294  DataExpanded::toString() const
294    //    //
295    // create a temporary view as the offset will be changed    // create a temporary view as the offset will be changed
296    DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());    DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
297    for (int i=0;i<m_data.getNumRows();++i) {    for (int i=0;i<m_data.getNumRows();i++) {
298      for (int j=0;j<m_data.getNumCols();++j) {      for (int j=0;j<m_data.getNumCols();j++) {
299        tempView.setOffset(m_data.index(i,j));        tempView.setOffset(m_data.index(i,j));
300        stringstream suffix;        stringstream suffix;
301        suffix << "(" << i << "," << j << ")";        suffix << "(" << i << "," << j << ")";
# Line 289  DataExpanded::getDataPoint(int sampleNo, Line 326  DataExpanded::getDataPoint(int sampleNo,
326  DataArrayView::ValueType::size_type  DataArrayView::ValueType::size_type
327  DataExpanded::getLength() const  DataExpanded::getLength() const
328  {  {
329    return m_data.getData().size();    return m_data.size();
330    }
331    
332    void
333    DataExpanded::setRefValue(int ref,
334                              const DataArray& value)
335    {
336      //
337      // Get the number of samples and data-points per sample.
338      int numSamples = getNumSamples();
339      int numDPPSample = getNumDPPSample();
340    
341      //
342      // Determine the sample number which corresponds to this reference number.
343      int sampleNo = -1;
344      int tempRef = -1;
345      for (int n=0; n<numSamples; n++) {
346        tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);
347        if (tempRef == ref) {
348          sampleNo = n;
349          break;
350        }
351      }
352      if (sampleNo == -1) {
353        throw DataException("DataExpanded::setRefValue error: invalid ref number supplied.");
354      }
355    
356      for (int n=0; n<numDPPSample; n++) {
357        //
358        // Get *each* data-point in the sample in turn.
359        DataArrayView pointView = getDataPoint(sampleNo, n);
360        //
361        // Assign the values in the DataArray to this data-point.
362        pointView.copy(value.getView());
363      }
364    }
365    
366    void
367    DataExpanded::getRefValue(int ref,
368                              DataArray& value)
369    {
370      //
371      // Get the number of samples and data-points per sample.
372      int numSamples = getNumSamples();
373      int numDPPSample = getNumDPPSample();
374    
375      //
376      // Determine the sample number which corresponds to this reference number
377      int sampleNo = -1;
378      int tempRef = -1;
379      for (int n=0; n<numSamples; n++) {
380        tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);
381        if (tempRef == ref) {
382          sampleNo = n;
383          break;
384        }
385      }
386      if (sampleNo == -1) {
387        throw DataException("DataExpanded::getRefValue error: invalid ref number supplied.");
388      }
389    
390      //
391      // Get the *first* data-point associated with this sample number.
392      DataArrayView pointView = getDataPoint(sampleNo, 0);
393    
394      //
395      // Load the values from this data-point into the DataArray
396      value.getView().copy(pointView);
397    }
398    
399    int
400    DataExpanded::archiveData(ofstream& archiveFile,
401                              const DataArrayView::ValueType::size_type noValues) const
402    {
403      return(m_data.archiveData(archiveFile, noValues));
404    }
405    
406    int
407    DataExpanded::extractData(ifstream& archiveFile,
408                              const DataArrayView::ValueType::size_type noValues)
409    {
410      return(m_data.extractData(archiveFile, noValues));
411    }
412    
413    void
414    DataExpanded::copyAll(const boost::python::numeric::array& value) {
415      //
416      // Get the number of samples and data-points per sample.
417      int numSamples = getNumSamples();
418      int numDataPointsPerSample = getNumDPPSample();
419      int dataPointRank = getPointDataView().getRank();
420      ShapeType dataPointShape = getPointDataView().getShape();
421      //
422      // check rank:
423      if (value.getrank()!=dataPointRank+1)
424           throw DataException("Rank of numarray does not match Data object rank");
425      if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
426           throw DataException("leading dimension of numarray is too small");
427      //
428      int dataPoint = 0;
429      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
430        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
431          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
432          if (dataPointRank==0) {
433             dataPointView()=extract<double>(value[dataPoint]);
434          } else if (dataPointRank==1) {
435             for (int i=0; i<dataPointShape[0]; i++) {
436                dataPointView(i)=extract<double>(value[dataPoint][i]);
437             }
438          } else if (dataPointRank==2) {
439             for (int i=0; i<dataPointShape[0]; i++) {
440               for (int j=0; j<dataPointShape[1]; j++) {
441                 dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
442               }
443             }
444           } else if (dataPointRank==3) {
445             for (int i=0; i<dataPointShape[0]; i++) {
446               for (int j=0; j<dataPointShape[1]; j++) {
447                 for (int k=0; k<dataPointShape[2]; k++) {
448                     dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);
449                 }
450               }
451             }
452           } else if (dataPointRank==4) {
453             for (int i=0; i<dataPointShape[0]; i++) {
454               for (int j=0; j<dataPointShape[1]; j++) {
455                 for (int k=0; k<dataPointShape[2]; k++) {
456                   for (int l=0; l<dataPointShape[3]; l++) {
457                     dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);
458                   }
459                 }
460               }
461             }
462          }
463          dataPoint++;
464        }
465      }
466    }
467    void
468    DataExpanded::eigenvalues(DataAbstract* ev)
469    {
470      int numSamples = getNumSamples();
471      int numDataPointsPerSample = getNumDPPSample();
472      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
473      if (temp_ev==0) {
474        throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
475      }
476      DataArrayView& thisView=getPointDataView();
477      DataArrayView& evView=ev->getPointDataView();
478      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
479      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
480        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
481             DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
482                                        evView,ev->getPointOffset(sampleNo,dataPointNo));
483        }
484      }
485    }
486    void
487    DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
488    {
489      int numSamples = getNumSamples();
490      int numDataPointsPerSample = getNumDPPSample();
491      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
492      if (temp_ev==0) {
493        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
494      }
495      DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
496      if (temp_V==0) {
497        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
498      }
499      DataArrayView& thisView=getPointDataView();
500      DataArrayView& evView=ev->getPointDataView();
501      DataArrayView& VView=V->getPointDataView();
502      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
503      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
504        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
505             DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
506                                                         evView,ev->getPointOffset(sampleNo,dataPointNo),
507                                                         VView,V->getPointOffset(sampleNo,dataPointNo),
508                                                         tol);
509        }
510      }
511  }  }
512    
513  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.97  
changed lines
  Added in v.580

  ViewVC Help
Powered by ViewVC 1.1.26