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

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

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

revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC revision 126 by jgs, Fri Jul 22 03:53:08 2005 UTC
# Line 18  Line 18 
18  #include "escript/Data/Data.h"  #include "escript/Data/Data.h"
19    
20  #include <iostream>  #include <iostream>
21    #include <fstream>
22  #include <algorithm>  #include <algorithm>
23  #include <vector>  #include <vector>
24  #include <exception>  #include <exception>
# Line 29  Line 30 
30  #include <boost/python/long.hpp>  #include <boost/python/long.hpp>
31    
32  #include "escript/Data/DataException.h"  #include "escript/Data/DataException.h"
   
33  #include "escript/Data/DataExpanded.h"  #include "escript/Data/DataExpanded.h"
34  #include "escript/Data/DataConstant.h"  #include "escript/Data/DataConstant.h"
35  #include "escript/Data/DataTagged.h"  #include "escript/Data/DataTagged.h"
36  #include "escript/Data/DataEmpty.h"  #include "escript/Data/DataEmpty.h"
37  #include "escript/Data/DataArray.h"  #include "escript/Data/DataArray.h"
38  #include "escript/Data/DataAlgorithm.h"  #include "escript/Data/DataProf.h"
39  #include "escript/Data/FunctionSpaceFactory.h"  #include "escript/Data/FunctionSpaceFactory.h"
40  #include "escript/Data/AbstractContinuousDomain.h"  #include "escript/Data/AbstractContinuousDomain.h"
41    #include "escript/Data/UnaryFuncs.h"
42    
43  using namespace std;  using namespace std;
44  using namespace boost::python;  using namespace boost::python;
45  using namespace boost;  using namespace boost;
46  using namespace escript;  using namespace escript;
47    
48    //
49    // global table of profiling data for all Data objects
50    DataProf dataProfTable;
51    
52  Data::Data()  Data::Data()
53  {  {
54    //    //
55    // Default data is type DataEmpty    // Default data is type DataEmpty
56    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
57    m_data=shared_ptr<DataAbstract>(temp);    shared_ptr<DataAbstract> temp_data(temp);
58      m_data=temp_data;
59      // create entry in global profiling table for this object
60      profData = dataProfTable.newData();
61  }  }
62    
63  Data::Data(double value,  Data::Data(double value,
# Line 63  Data::Data(double value, Line 71  Data::Data(double value,
71    }    }
72    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
73    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
74      // create entry in global profiling table for this object
75      profData = dataProfTable.newData();
76  }  }
77    
78  Data::Data(double value,  Data::Data(double value,
# Line 73  Data::Data(double value, Line 83  Data::Data(double value,
83    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
84    pair<int,int> dataShape=what.getDataShape();    pair<int,int> dataShape=what.getDataShape();
85    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
86      // create entry in global profiling table for this object
87      profData = dataProfTable.newData();
88  }  }
89    
90  Data::Data(const Data& inData):  Data::Data(const Data& inData)
   m_data(inData.m_data)  
91  {  {
92      m_data=inData.m_data;
93      // create entry in global profiling table for this object
94      profData = dataProfTable.newData();
95  }  }
96    
97  Data::Data(const Data& inData,  Data::Data(const Data& inData,
98             const DataArrayView::RegionType& region)             const DataArrayView::RegionType& region)
99  {  {
100    //    //
101    // Create data which is a subset(slice) of another Data    // Create Data which is a slice of another Data
102    DataAbstract* tmp=inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
103    m_data=shared_ptr<DataAbstract>(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
104      m_data=temp_data;
105      // create entry in global profiling table for this object
106      profData = dataProfTable.newData();
107  }  }
108    
109  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 96  Data::Data(const Data& inData, Line 113  Data::Data(const Data& inData,
113      m_data=inData.m_data;      m_data=inData.m_data;
114    } else {    } else {
115      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
116      // Note for Lutz, Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
117      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
118      // be able to create an instance of AbstractDomain but that was done      // be able to create an instance of AbstractDomain but that was done
119      // as a boost python work around which may no longer be required.      // as a boost:python work around which may no longer be required.
120      const AbstractDomain& inDataDomain=inData.getDomain();      const AbstractDomain& inDataDomain=inData.getDomain();
121      if  (inDataDomain==functionspace.getDomain()) {      if  (inDataDomain==functionspace.getDomain()) {
122        inDataDomain.interpolateOnDomain(tmp,inData);        inDataDomain.interpolateOnDomain(tmp,inData);
# Line 108  Data::Data(const Data& inData, Line 125  Data::Data(const Data& inData,
125      }      }
126      m_data=tmp.m_data;      m_data=tmp.m_data;
127    }    }
128      // create entry in global profiling table for this object
129      profData = dataProfTable.newData();
130  }  }
131    
132  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(const DataTagged::TagListType& tagKeys,
# Line 117  Data::Data(const DataTagged::TagListType Line 136  Data::Data(const DataTagged::TagListType
136             bool expanded)             bool expanded)
137  {  {
138    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
139    m_data=shared_ptr<DataAbstract>(temp);    shared_ptr<DataAbstract> temp_data(temp);
140      m_data=temp_data;
141    if (expanded) {    if (expanded) {
142      expand();      expand();
143    }    }
144      // create entry in global profiling table for this object
145      profData = dataProfTable.newData();
146  }  }
147    
148  Data::Data(const numeric::array& value,  Data::Data(const numeric::array& value,
# Line 128  Data::Data(const numeric::array& value, Line 150  Data::Data(const numeric::array& value,
150             bool expanded)             bool expanded)
151  {  {
152    initialise(value,what,expanded);    initialise(value,what,expanded);
153      // create entry in global profiling table for this object
154      profData = dataProfTable.newData();
155  }  }
156    
157  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
# Line 135  Data::Data(const DataArrayView& value, Line 159  Data::Data(const DataArrayView& value,
159             bool expanded)             bool expanded)
160  {  {
161    initialise(value,what,expanded);    initialise(value,what,expanded);
162      // create entry in global profiling table for this object
163      profData = dataProfTable.newData();
164  }  }
165    
166  Data::Data(const object& value,  Data::Data(const object& value,
# Line 143  Data::Data(const object& value, Line 169  Data::Data(const object& value,
169  {  {
170    numeric::array asNumArray(value);    numeric::array asNumArray(value);
171    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
172      // create entry in global profiling table for this object
173      profData = dataProfTable.newData();
174  }  }
175    
176  Data::Data(const object& value,  Data::Data(const object& value,
# Line 163  Data::Data(const object& value, Line 191  Data::Data(const object& value,
191      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
192      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
193    }    }
194      // create entry in global profiling table for this object
195      profData = dataProfTable.newData();
196  }  }
197    
198  escriptDataC  escriptDataC
# Line 181  Data::getDataC() const Line 211  Data::getDataC() const
211    return temp;    return temp;
212  }  }
213    
214  tuple  const boost::python::tuple
215  Data::getShapeTuple() const  Data::getShapeTuple() const
216  {  {
217    const DataArrayView::ShapeType& shape=getDataPointShape();    const DataArrayView::ShapeType& shape=getDataPointShape();
# Line 212  Data::copy(const Data& other) Line 242  Data::copy(const Data& other)
242        //        //
243        // Construct a DataExpanded copy        // Construct a DataExpanded copy
244        DataAbstract* newData=new DataExpanded(*temp);        DataAbstract* newData=new DataExpanded(*temp);
245        m_data=shared_ptr<DataAbstract>(newData);        shared_ptr<DataAbstract> temp_data(newData);
246          m_data=temp_data;
247        return;        return;
248      }      }
249    }    }
# Line 220  Data::copy(const Data& other) Line 251  Data::copy(const Data& other)
251      DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());      DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());
252      if (temp!=0) {      if (temp!=0) {
253        //        //
254        // Construct a DataTaggeded copy        // Construct a DataTagged copy
255        DataAbstract* newData=new DataTagged(*temp);        DataAbstract* newData=new DataTagged(*temp);
256        m_data=shared_ptr<DataAbstract>(newData);        shared_ptr<DataAbstract> temp_data(newData);
257          m_data=temp_data;
258        return;        return;
259      }      }
260    }    }
# Line 232  Data::copy(const Data& other) Line 264  Data::copy(const Data& other)
264        //        //
265        // Construct a DataConstant copy        // Construct a DataConstant copy
266        DataAbstract* newData=new DataConstant(*temp);        DataAbstract* newData=new DataConstant(*temp);
267        m_data=shared_ptr<DataAbstract>(newData);        shared_ptr<DataAbstract> temp_data(newData);
268          m_data=temp_data;
269          return;
270        }
271      }
272      {
273        DataEmpty* temp=dynamic_cast<DataEmpty*>(other.m_data.get());
274        if (temp!=0) {
275          //
276          // Construct a DataEmpty copy
277          DataAbstract* newData=new DataEmpty();
278          shared_ptr<DataAbstract> temp_data(newData);
279          m_data=temp_data;
280        return;        return;
281      }      }
282    }    }
# Line 284  Data::isConstant() const Line 328  Data::isConstant() const
328    return (temp!=0);    return (temp!=0);
329  }  }
330    
 Data  
 Data::getSlice(const DataArrayView::RegionType& region) const  
 {  
   return Data(*this,region);  
 }  
   
 void  
 Data::setSlice(const Data& value,  
                const DataArrayView::RegionType& region)  
 {  
   m_data->setSlice(value.m_data.get(), region);  
 }  
   
331  void  void
332  Data::expand()  Data::expand()
333  {  {
334    if (isConstant()) {    if (isConstant()) {
335      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
336      DataAbstract* temp=new DataExpanded(*tempDataConst);      DataAbstract* temp=new DataExpanded(*tempDataConst);
337      m_data=shared_ptr<DataAbstract>(temp);      shared_ptr<DataAbstract> temp_data(temp);
338        m_data=temp_data;
339    } else if (isTagged()) {    } else if (isTagged()) {
340      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
341      DataAbstract* temp=new DataExpanded(*tempDataTag);      DataAbstract* temp=new DataExpanded(*tempDataTag);
342      m_data=shared_ptr<DataAbstract>(temp);      shared_ptr<DataAbstract> temp_data(temp);
343        m_data=temp_data;
344    } else if (isExpanded()) {    } else if (isExpanded()) {
345      //      //
346      // do nothing      // do nothing
# Line 319  Data::expand() Line 352  Data::expand()
352  }  }
353    
354  void  void
 Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  
 {  
   m_data->reshapeDataPoint(shape);  
 }  
   
 void  
355  Data::tag()  Data::tag()
356  {  {
357    if (isConstant()) {    if (isConstant()) {
358      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
359      DataAbstract* temp=new DataTagged(*tempDataConst);      DataAbstract* temp=new DataTagged(*tempDataConst);
360      m_data=shared_ptr<DataAbstract>(temp);      shared_ptr<DataAbstract> temp_data(temp);
361        m_data=temp_data;
362    } else if (isTagged()) {    } else if (isTagged()) {
363      // do nothing      // do nothing
364    } else if (isExpanded()) {    } else if (isExpanded()) {
# Line 342  Data::tag() Line 370  Data::tag()
370    }    }
371  }  }
372    
373    void
374    Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)
375    {
376      m_data->reshapeDataPoint(shape);
377    }
378    
379  Data  Data
380  Data::wherePositive() const  Data::wherePositive() const
381  {  {
382      profData->where++;
383    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
384  }  }
385    
386  Data  Data
387    Data::whereNegative() const
388    {
389      profData->where++;
390      return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
391    }
392    
393    Data
394  Data::whereNonNegative() const  Data::whereNonNegative() const
395  {  {
396      profData->where++;
397    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
398  }  }
399    
400  Data  Data
401  Data::whereNegative() const  Data::whereNonPositive() const
402  {  {
403    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));    profData->where++;
404      return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
405  }  }
406    
407  Data  Data
408  Data::whereZero() const  Data::whereZero() const
409  {  {
410      profData->where++;
411    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));
412  }  }
413    
414  Data  Data
415    Data::whereNonZero() const
416    {
417      profData->where++;
418      return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));
419    }
420    
421    Data
422  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
423  {  {
424      profData->interpolate++;
425    return Data(*this,functionspace);    return Data(*this,functionspace);
426  }  }
427    
# Line 390  Data::probeInterpolation(const FunctionS Line 443  Data::probeInterpolation(const FunctionS
443  Data  Data
444  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
445  {  {
446      profData->grad++;
447    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
448      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
449    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
# Line 423  Data::getDataPointShape() const Line 477  Data::getDataPointShape() const
477    return getPointDataView().getShape();    return getPointDataView().getShape();
478  }  }
479    
480    void
481    Data::fillFromNumArray(const boost::python::numeric::array num_array)
482    {
483      //
484      // check rank:
485      if (num_array.getrank()<getDataPointRank())
486          throw DataException("Rank of numarray does not match Data object rank");
487      //
488      // check rank of num_array
489      for (int i=0; i<getDataPointRank(); i++) {
490        if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])
491           throw DataException("Shape of numarray does not match Data object rank");
492      }
493      //
494      // make sure data is expanded:
495      if (! isExpanded()) expand();
496      //
497      // and copy over:
498      m_data->copyAll(num_array);
499      //
500      // the rank of the returned numeric array will be the rank of
501      // the data points, plus one. Where the rank of the array is n,
502      // the last n-1 dimensions will be equal to the shape of the
503      // data points, whilst the first dimension will be equal to the
504      // total number of data points. Thus the array will consist of
505      // a serial vector of the data points.
506      // int arrayRank = dataPointRank + 1;
507      // DataArrayView::ShapeType arrayShape;
508      // arrayShape.push_back(numDataPoints);
509      // for (int d=0; d<dataPointRank; d++) {
510      //    arrayShape.push_back(dataPointShape[d]);
511      // }
512    
513      //
514      // resize the numeric array to the shape just calculated
515      // if (arrayRank==1) {
516      //   numArray.resize(arrayShape[0]);
517      // }
518      // if (arrayRank==2) {
519      //   numArray.resize(arrayShape[0],arrayShape[1]);
520      // }
521      // if (arrayRank==3) {
522      //   numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
523      // }
524      // if (arrayRank==4) {
525      //   numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
526      // }
527      // if (arrayRank==5) {
528      //   numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
529      // }
530    
531      //
532      // loop through each data point in turn, loading the values for that data point
533      // into the numeric array.
534      // int dataPoint = 0;
535      // for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
536      //   for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
537      //     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
538      //     if (dataPointRank==0) {
539      //       dataPointView()=numArray[dataPoint];
540      //     }
541      //     if (dataPointRank==1) {
542      //       for (int i=0; i<dataPointShape[0]; i++) {
543      //         dataPointView(i)=numArray[dataPoint][i];
544      //       }
545      //     }
546      //     if (dataPointRank==2) {
547      //       for (int i=0; i<dataPointShape[0]; i++) {
548      //         for (int j=0; j<dataPointShape[1]; j++) {
549      //           numArray[dataPoint][i][j] = dataPointView(i,j);
550      //         }
551      //       }
552      //     }
553      //     if (dataPointRank==3) {
554      //       for (int i=0; i<dataPointShape[0]; i++) {
555      //         for (int j=0; j<dataPointShape[1]; j++) {
556      //           for (int k=0; k<dataPointShape[2]; k++) {
557      //             numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
558      //           }
559      //         }
560      //       }
561      //     }
562      //     if (dataPointRank==4) {
563      //       for (int i=0; i<dataPointShape[0]; i++) {
564      //         for (int j=0; j<dataPointShape[1]; j++) {
565      //           for (int k=0; k<dataPointShape[2]; k++) {
566      //             for (int l=0; l<dataPointShape[3]; l++) {
567      //               numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
568      //             }
569      //           }
570      //         }
571      //       }
572      //     }
573      //     dataPoint++;
574      //   }
575      // }
576    
577      //
578      // return the loaded array
579      // return numArray;
580    
581    }
582    
583    const
584    boost::python::numeric::array
585    Data::convertToNumArray()
586    {
587      //
588      // determine the total number of data points
589      int numSamples = getNumSamples();
590      int numDataPointsPerSample = getNumDataPointsPerSample();
591      int numDataPoints = numSamples * numDataPointsPerSample;
592    
593      //
594      // determine the rank and shape of each data point
595      int dataPointRank = getDataPointRank();
596      DataArrayView::ShapeType dataPointShape = getDataPointShape();
597    
598      //
599      // create the numeric array to be returned
600      boost::python::numeric::array numArray(0.0);
601    
602      //
603      // the rank of the returned numeric array will be the rank of
604      // the data points, plus one. Where the rank of the array is n,
605      // the last n-1 dimensions will be equal to the shape of the
606      // data points, whilst the first dimension will be equal to the
607      // total number of data points. Thus the array will consist of
608      // a serial vector of the data points.
609      int arrayRank = dataPointRank + 1;
610      DataArrayView::ShapeType arrayShape;
611      arrayShape.push_back(numDataPoints);
612      for (int d=0; d<dataPointRank; d++) {
613         arrayShape.push_back(dataPointShape[d]);
614      }
615    
616      //
617      // resize the numeric array to the shape just calculated
618      if (arrayRank==1) {
619        numArray.resize(arrayShape[0]);
620      }
621      if (arrayRank==2) {
622        numArray.resize(arrayShape[0],arrayShape[1]);
623      }
624      if (arrayRank==3) {
625        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
626      }
627      if (arrayRank==4) {
628        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
629      }
630      if (arrayRank==5) {
631        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
632      }
633    
634      //
635      // loop through each data point in turn, loading the values for that data point
636      // into the numeric array.
637      int dataPoint = 0;
638      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
639        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
640          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
641          if (dataPointRank==0) {
642            numArray[dataPoint]=dataPointView();
643          }
644          if (dataPointRank==1) {
645            for (int i=0; i<dataPointShape[0]; i++) {
646              numArray[dataPoint][i]=dataPointView(i);
647            }
648          }
649          if (dataPointRank==2) {
650            for (int i=0; i<dataPointShape[0]; i++) {
651              for (int j=0; j<dataPointShape[1]; j++) {
652                numArray[dataPoint][i][j] = dataPointView(i,j);
653              }
654            }
655          }
656          if (dataPointRank==3) {
657            for (int i=0; i<dataPointShape[0]; i++) {
658              for (int j=0; j<dataPointShape[1]; j++) {
659                for (int k=0; k<dataPointShape[2]; k++) {
660                  numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
661                }
662              }
663            }
664          }
665          if (dataPointRank==4) {
666            for (int i=0; i<dataPointShape[0]; i++) {
667              for (int j=0; j<dataPointShape[1]; j++) {
668                for (int k=0; k<dataPointShape[2]; k++) {
669                  for (int l=0; l<dataPointShape[3]; l++) {
670                    numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
671                  }
672                }
673              }
674            }
675          }
676          dataPoint++;
677        }
678      }
679    
680      //
681      // return the loaded array
682      return numArray;
683    }
684    
685    const
686    boost::python::numeric::array
687    Data::convertToNumArrayFromSampleNo(int sampleNo)
688    {
689      //
690      // Check a valid sample number has been supplied
691      if (sampleNo >= getNumSamples()) {
692        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
693      }
694    
695      //
696      // determine the number of data points per sample
697      int numDataPointsPerSample = getNumDataPointsPerSample();
698    
699      //
700      // determine the rank and shape of each data point
701      int dataPointRank = getDataPointRank();
702      DataArrayView::ShapeType dataPointShape = getDataPointShape();
703    
704      //
705      // create the numeric array to be returned
706      boost::python::numeric::array numArray(0.0);
707    
708      //
709      // the rank of the returned numeric array will be the rank of
710      // the data points, plus one. Where the rank of the array is n,
711      // the last n-1 dimensions will be equal to the shape of the
712      // data points, whilst the first dimension will be equal to the
713      // total number of data points. Thus the array will consist of
714      // a serial vector of the data points.
715      int arrayRank = dataPointRank + 1;
716      DataArrayView::ShapeType arrayShape;
717      arrayShape.push_back(numDataPointsPerSample);
718      for (int d=0; d<dataPointRank; d++) {
719         arrayShape.push_back(dataPointShape[d]);
720      }
721    
722      //
723      // resize the numeric array to the shape just calculated
724      if (arrayRank==1) {
725        numArray.resize(arrayShape[0]);
726      }
727      if (arrayRank==2) {
728        numArray.resize(arrayShape[0],arrayShape[1]);
729      }
730      if (arrayRank==3) {
731        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
732      }
733      if (arrayRank==4) {
734        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
735      }
736      if (arrayRank==5) {
737        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
738      }
739    
740      //
741      // loop through each data point in turn, loading the values for that data point
742      // into the numeric array.
743      for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {
744        DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);
745        if (dataPointRank==0) {
746          numArray[dataPoint]=dataPointView();
747        }
748        if (dataPointRank==1) {
749          for (int i=0; i<dataPointShape[0]; i++) {
750            numArray[dataPoint][i]=dataPointView(i);
751          }
752        }
753        if (dataPointRank==2) {
754          for (int i=0; i<dataPointShape[0]; i++) {
755            for (int j=0; j<dataPointShape[1]; j++) {
756              numArray[dataPoint][i][j] = dataPointView(i,j);
757            }
758          }
759        }
760        if (dataPointRank==3) {
761          for (int i=0; i<dataPointShape[0]; i++) {
762            for (int j=0; j<dataPointShape[1]; j++) {
763              for (int k=0; k<dataPointShape[2]; k++) {
764                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
765              }
766            }
767          }
768        }
769        if (dataPointRank==4) {
770          for (int i=0; i<dataPointShape[0]; i++) {
771            for (int j=0; j<dataPointShape[1]; j++) {
772              for (int k=0; k<dataPointShape[2]; k++) {
773                for (int l=0; l<dataPointShape[3]; l++) {
774                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
775                }
776              }
777            }
778          }
779        }
780      }
781    
782      //
783      // return the loaded array
784      return numArray;
785    }
786    
787    const
788    boost::python::numeric::array
789    Data::convertToNumArrayFromDPNo(int sampleNo,
790                                    int dataPointNo)
791    {
792      //
793      // Check a valid sample number has been supplied
794      if (sampleNo >= getNumSamples()) {
795        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
796      }
797    
798      //
799      // Check a valid data point number has been supplied
800      if (dataPointNo >= getNumDataPointsPerSample()) {
801        throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");
802      }
803    
804      //
805      // determine the rank and shape of each data point
806      int dataPointRank = getDataPointRank();
807      DataArrayView::ShapeType dataPointShape = getDataPointShape();
808    
809      //
810      // create the numeric array to be returned
811      boost::python::numeric::array numArray(0.0);
812    
813      //
814      // the shape of the returned numeric array will be the same
815      // as that of the data point
816      int arrayRank = dataPointRank;
817      DataArrayView::ShapeType arrayShape = dataPointShape;
818    
819      //
820      // resize the numeric array to the shape just calculated
821      if (arrayRank==0) {
822        numArray.resize(1);
823      }
824      if (arrayRank==1) {
825        numArray.resize(arrayShape[0]);
826      }
827      if (arrayRank==2) {
828        numArray.resize(arrayShape[0],arrayShape[1]);
829      }
830      if (arrayRank==3) {
831        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
832      }
833      if (arrayRank==4) {
834        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
835      }
836    
837      //
838      // load the values for the data point into the numeric array.
839      DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
840      if (dataPointRank==0) {
841        numArray[0]=dataPointView();
842      }
843      if (dataPointRank==1) {
844        for (int i=0; i<dataPointShape[0]; i++) {
845          numArray[i]=dataPointView(i);
846        }
847      }
848      if (dataPointRank==2) {
849        for (int i=0; i<dataPointShape[0]; i++) {
850          for (int j=0; j<dataPointShape[1]; j++) {
851            numArray[i][j] = dataPointView(i,j);
852          }
853        }
854      }
855      if (dataPointRank==3) {
856        for (int i=0; i<dataPointShape[0]; i++) {
857          for (int j=0; j<dataPointShape[1]; j++) {
858            for (int k=0; k<dataPointShape[2]; k++) {
859              numArray[i][j][k]=dataPointView(i,j,k);
860            }
861          }
862        }
863      }
864      if (dataPointRank==4) {
865        for (int i=0; i<dataPointShape[0]; i++) {
866          for (int j=0; j<dataPointShape[1]; j++) {
867            for (int k=0; k<dataPointShape[2]; k++) {
868              for (int l=0; l<dataPointShape[3]; l++) {
869                numArray[i][j][k][l]=dataPointView(i,j,k,l);
870              }
871            }
872          }
873        }
874      }
875    
876      //
877      // return the loaded array
878      return numArray;
879    }
880    
881  boost::python::numeric::array  boost::python::numeric::array
882  Data::integrate() const  Data::integrate() const
883  {  {
# Line 430  Data::integrate() const Line 885  Data::integrate() const
885    int rank = getDataPointRank();    int rank = getDataPointRank();
886    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
887    
888      profData->integrate++;
889    
890    //    //
891    // calculate the integral values    // calculate the integral values
892    vector<double> integrals(getDataPointSize());    vector<double> integrals(getDataPointSize());
# Line 440  Data::integrate() const Line 897  Data::integrate() const
897    // and load the array with the integral values    // and load the array with the integral values
898    boost::python::numeric::array bp_array(1.0);    boost::python::numeric::array bp_array(1.0);
899    if (rank==0) {    if (rank==0) {
900        bp_array.resize(1);
901      index = 0;      index = 0;
902      bp_array[0] = integrals[index];      bp_array[0] = integrals[index];
903    }    }
# Line 492  Data::integrate() const Line 950  Data::integrate() const
950  Data  Data
951  Data::sin() const  Data::sin() const
952  {  {
953      profData->unary++;
954    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
955  }  }
956    
957  Data  Data
958  Data::cos() const  Data::cos() const
959  {  {
960      profData->unary++;
961    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
962  }  }
963    
964  Data  Data
965  Data::tan() const  Data::tan() const
966  {  {
967      profData->unary++;
968    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
969  }  }
970    
971  Data  Data
972  Data::log() const  Data::log() const
973  {  {
974      profData->unary++;
975    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
976  }  }
977    
978  Data  Data
979  Data::ln() const  Data::ln() const
980  {  {
981      profData->unary++;
982    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
983  }  }
984    
985    Data
986    Data::sign() const
987    {
988      profData->unary++;
989      return escript::unaryOp(*this,escript::fsign);
990    }
991    
992    Data
993    Data::abs() const
994    {
995      profData->unary++;
996      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
997    }
998    
999    Data
1000    Data::neg() const
1001    {
1002      profData->unary++;
1003      return escript::unaryOp(*this,negate<double>());
1004    }
1005    
1006    Data
1007    Data::pos() const
1008    {
1009      profData->unary++;
1010      return (*this);
1011    }
1012    
1013    Data
1014    Data::exp() const
1015    {
1016      profData->unary++;
1017      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1018    }
1019    
1020    Data
1021    Data::sqrt() const
1022    {
1023      profData->unary++;
1024      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1025    }
1026    
1027  double  double
1028  Data::Lsup() const  Data::Lsup() const
1029  {  {
1030      profData->reduction1++;
1031      //
1032      // set the initial absolute maximum value to zero
1033      return algorithm(DataAlgorithmAdapter<AbsMax>(0));
1034    }
1035    
1036    double
1037    Data::Linf() const
1038    {
1039      profData->reduction1++;
1040    //    //
1041    // set the initial absolute maximum value to min possible    // set the initial absolute minimum value to max double
1042    return algorithm(DataAlgorithmAdapter<AbsMax>(numeric_limits<double>::min()));    return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));
1043  }  }
1044    
1045  double  double
1046  Data::sup() const  Data::sup() const
1047  {  {
1048      profData->reduction1++;
1049    //    //
1050    // set the initial maximum value to min possible    // set the initial maximum value to min possible double
1051    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::min()));    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
1052  }  }
1053    
1054  double  double
1055  Data::inf() const  Data::inf() const
1056  {  {
1057      profData->reduction1++;
1058    //    //
1059    // set the initial minimum value to max possible    // set the initial minimum value to max possible double
1060    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
1061  }  }
1062    
1063  Data& Data::operator+=(const Data& right)  Data
1064    Data::maxval() const
1065    {
1066      profData->reduction2++;
1067      //
1068      // set the initial maximum value to min possible double
1069      return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
1070    }
1071    
1072    Data
1073    Data::minval() const
1074    {
1075      profData->reduction2++;
1076      //
1077      // set the initial minimum value to max possible double
1078      return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
1079    }
1080    
1081    Data
1082    Data::length() const
1083    {
1084      profData->reduction2++;
1085      return dp_algorithm(DataAlgorithmAdapter<Length>(0));
1086    }
1087    
1088    Data
1089    Data::trace() const
1090    {
1091      profData->reduction2++;
1092      return dp_algorithm(DataAlgorithmAdapter<Trace>(0));
1093    }
1094    
1095    Data
1096    Data::transpose(int axis) const
1097  {  {
1098      profData->reduction2++;
1099      // not implemented
1100      throw DataException("Error - Data::transpose not implemented yet.");
1101      return Data();
1102    }
1103    
1104    const boost::python::tuple
1105    Data::mindp() const
1106    {
1107      Data temp=minval();
1108    
1109      int numSamples=temp.getNumSamples();
1110      int numDPPSample=temp.getNumDataPointsPerSample();
1111    
1112      int i,j,lowi=0,lowj=0;
1113      double min=numeric_limits<double>::max();
1114    
1115      for (i=0; i<numSamples; i++) {
1116        for (j=0; j<numDPPSample; j++) {
1117          double next=temp.getDataPoint(i,j)();
1118          if (next<min) {
1119            min=next;
1120            lowi=i;
1121            lowj=j;
1122          }
1123        }
1124      }
1125    
1126      return make_tuple(lowi,lowj);
1127    }
1128    
1129    void
1130    Data::saveDX(std::string fileName) const
1131    {
1132      getDomain().saveDX(fileName,*this);
1133      return;
1134    }
1135    
1136    void
1137    Data::saveVTK(std::string fileName) const
1138    {
1139      getDomain().saveVTK(fileName,*this);
1140      return;
1141    }
1142    
1143    Data&
1144    Data::operator+=(const Data& right)
1145    {
1146      profData->binary++;
1147    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1148    return (*this);    return (*this);
1149  }  }
1150    
1151  Data& Data::operator+=(const boost::python::object& right)  Data&
1152    Data::operator+=(const boost::python::object& right)
1153  {  {
1154      profData->binary++;
1155    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1156    return (*this);    return (*this);
1157  }  }
1158    
1159  Data& Data::operator-=(const Data& right)  Data&
1160    Data::operator-=(const Data& right)
1161  {  {
1162      profData->binary++;
1163    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1164    return (*this);    return (*this);
1165  }  }
1166    
1167  Data& Data::operator-=(const boost::python::object& right)  Data&
1168    Data::operator-=(const boost::python::object& right)
1169  {  {
1170      profData->binary++;
1171    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1172    return (*this);    return (*this);
1173  }  }
1174    
1175  Data& Data::operator*=(const Data& right)  Data&
1176    Data::operator*=(const Data& right)
1177  {  {
1178      profData->binary++;
1179    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1180    return (*this);    return (*this);
1181  }  }
1182    
1183  Data& Data::operator*=(const boost::python::object& right)  Data&
1184    Data::operator*=(const boost::python::object& right)
1185  {  {
1186      profData->binary++;
1187    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1188    return (*this);    return (*this);
1189  }  }
1190    
1191  Data& Data::operator/=(const Data& right)  Data&
1192    Data::operator/=(const Data& right)
1193  {  {
1194      profData->binary++;
1195    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1196    return (*this);    return (*this);
1197  }  }
1198    
1199  Data& Data::operator/=(const boost::python::object& right)  Data&
1200    Data::operator/=(const boost::python::object& right)
1201  {  {
1202      profData->binary++;
1203    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1204    return (*this);    return (*this);
1205  }  }
1206    
1207  Data Data::powO(const boost::python::object& right) const  Data
1208    Data::powO(const boost::python::object& right) const
1209  {  {
1210      profData->binary++;
1211    Data result;    Data result;
1212    result.copy(*this);    result.copy(*this);
1213    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1214    return result;    return result;
1215  }  }
1216    
1217  Data Data::powD(const Data& right) const  Data
1218    Data::powD(const Data& right) const
1219  {  {
1220      profData->binary++;
1221    Data result;    Data result;
1222    result.copy(*this);    result.copy(*this);
1223    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 608  Data Data::powD(const Data& right) const Line 1225  Data Data::powD(const Data& right) const
1225  }  }
1226    
1227  //  //
1228  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1229  Data escript::operator+(const Data& left, const Data& right)  Data
1230    escript::operator+(const Data& left, const Data& right)
1231  {  {
1232    Data result;    Data result;
1233    //    //
# Line 620  Data escript::operator+(const Data& left Line 1238  Data escript::operator+(const Data& left
1238  }  }
1239    
1240  //  //
1241  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1242  Data escript::operator-(const Data& left, const Data& right)  Data
1243    escript::operator-(const Data& left, const Data& right)
1244  {  {
1245    Data result;    Data result;
1246    //    //
# Line 632  Data escript::operator-(const Data& left Line 1251  Data escript::operator-(const Data& left
1251  }  }
1252    
1253  //  //
1254  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1255  Data escript::operator*(const Data& left, const Data& right)  Data
1256    escript::operator*(const Data& left, const Data& right)
1257  {  {
1258    Data result;    Data result;
1259    //    //
# Line 644  Data escript::operator*(const Data& left Line 1264  Data escript::operator*(const Data& left
1264  }  }
1265    
1266  //  //
1267  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1268  Data escript::operator/(const Data& left, const Data& right)  Data
1269    escript::operator/(const Data& left, const Data& right)
1270  {  {
1271    Data result;    Data result;
1272    //    //
# Line 656  Data escript::operator/(const Data& left Line 1277  Data escript::operator/(const Data& left
1277  }  }
1278    
1279  //  //
1280  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1281  Data escript::operator+(const Data& left, const boost::python::object& right)  Data
1282    escript::operator+(const Data& left, const boost::python::object& right)
1283  {  {
1284    //    //
1285    // Convert to DataArray format if possible    // Convert to DataArray format if possible
# Line 671  Data escript::operator+(const Data& left Line 1293  Data escript::operator+(const Data& left
1293  }  }
1294    
1295  //  //
1296  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1297  Data escript::operator-(const Data& left, const boost::python::object& right)  Data
1298    escript::operator-(const Data& left, const boost::python::object& right)
1299  {  {
1300    //    //
1301    // Convert to DataArray format if possible    // Convert to DataArray format if possible
# Line 686  Data escript::operator-(const Data& left Line 1309  Data escript::operator-(const Data& left
1309  }  }
1310    
1311  //  //
1312  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1313  Data escript::operator*(const Data& left, const boost::python::object& right)  Data
1314    escript::operator*(const Data& left, const boost::python::object& right)
1315  {  {
1316    //    //
1317    // Convert to DataArray format if possible    // Convert to DataArray format if possible
# Line 701  Data escript::operator*(const Data& left Line 1325  Data escript::operator*(const Data& left
1325  }  }
1326    
1327  //  //
1328  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1329  Data escript::operator/(const Data& left, const boost::python::object& right)  Data
1330    escript::operator/(const Data& left, const boost::python::object& right)
1331  {  {
1332    //    //
1333    // Convert to DataArray format if possible    // Convert to DataArray format if possible
# Line 716  Data escript::operator/(const Data& left Line 1341  Data escript::operator/(const Data& left
1341  }  }
1342    
1343  //  //
1344  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1345  Data escript::operator+(const boost::python::object& left, const Data& right)  Data
1346    escript::operator+(const boost::python::object& left, const Data& right)
1347  {  {
1348    //    //
1349    // Construct the result using the given value and the other parameters    // Construct the result using the given value and the other parameters
# Line 728  Data escript::operator+(const boost::pyt Line 1354  Data escript::operator+(const boost::pyt
1354  }  }
1355    
1356  //  //
1357  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1358  Data escript::operator-(const boost::python::object& left, const Data& right)  Data
1359    escript::operator-(const boost::python::object& left, const Data& right)
1360  {  {
1361    //    //
1362    // Construct the result using the given value and the other parameters    // Construct the result using the given value and the other parameters
# Line 740  Data escript::operator-(const boost::pyt Line 1367  Data escript::operator-(const boost::pyt
1367  }  }
1368    
1369  //  //
1370  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1371  Data escript::operator*(const boost::python::object& left, const Data& right)  Data
1372    escript::operator*(const boost::python::object& left, const Data& right)
1373  {  {
1374    //    //
1375    // Construct the result using the given value and the other parameters    // Construct the result using the given value and the other parameters
# Line 752  Data escript::operator*(const boost::pyt Line 1380  Data escript::operator*(const boost::pyt
1380  }  }
1381    
1382  //  //
1383  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1384  Data escript::operator/(const boost::python::object& left, const Data& right)  Data
1385    escript::operator/(const boost::python::object& left, const Data& right)
1386  {  {
1387    //    //
1388    // Construct the result using the given value and the other parameters    // Construct the result using the given value and the other parameters
# Line 764  Data escript::operator/(const boost::pyt Line 1393  Data escript::operator/(const boost::pyt
1393  }  }
1394    
1395  //  //
1396  // NOTE: It is essential to specify the namepsace this operator belongs to  //bool escript::operator==(const Data& left, const Data& right)
1397  bool escript::operator==(const Data& left, const Data& right)  //{
1398  {  //  /*
1399    /*  //  NB: this operator does very little at this point, and isn't to
1400    NB: this operator does very little at this point, and isn't to  //  be relied on. Requires further implementation.
1401    be relied on. Requires further implementation.  //  */
1402    */  //
1403    //  bool ret;
1404    bool ret;  //
1405    //  if (left.isEmpty()) {
1406    if (left.isEmpty()) {  //    if(!right.isEmpty()) {
1407      if(!right.isEmpty()) {  //      ret = false;
1408        ret = false;  //    } else {
1409      } else {  //      ret = true;
1410        ret = true;  //    }
1411      }  //  }
1412    }  //
1413    //  if (left.isConstant()) {
1414    //    if(!right.isConstant()) {
1415    //      ret = false;
1416    //    } else {
1417    //      ret = true;
1418    //    }
1419    // }
1420    //
1421    //  if (left.isTagged()) {
1422    //   if(!right.isTagged()) {
1423    //      ret = false;
1424    //    } else {
1425    //      ret = true;
1426    //    }
1427    //  }
1428    //
1429    //  if (left.isExpanded()) {
1430    //    if(!right.isExpanded()) {
1431    //      ret = false;
1432    //    } else {
1433    //      ret = true;
1434    //    }
1435    //  }
1436    //
1437    //  return ret;
1438    //}
1439    
1440    if (left.isConstant()) {  Data
1441      if(!right.isConstant()) {  Data::getItem(const boost::python::object& key) const
1442        ret = false;  {
1443      } else {    const DataArrayView& view=getPointDataView();
       ret = true;  
     }  
   }  
1444    
1445    if (left.isTagged()) {    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
     if(!right.isTagged()) {  
       ret = false;  
     } else {  
       ret = true;  
     }  
   }  
1446    
1447    if (left.isExpanded()) {    if (slice_region.size()!=view.getRank()) {
1448      if(!right.isExpanded()) {      throw DataException("Error - slice size does not match Data rank.");
       ret = false;  
     } else {  
       ret = true;  
     }  
1449    }    }
1450    
1451    return ret;    return getSlice(slice_region);
1452  }  }
1453    
1454  Data  Data
1455  Data::getItem(const boost::python::object& key) const  Data::getSlice(const DataArrayView::RegionType& region) const
1456  {  {
1457     const DataArrayView& view=getPointDataView();    profData->slicing++;
1458     DataArrayView::RegionType slice_region=view.getSliceRegion(key);    return Data(*this,region);
    if (slice_region.size()!=view.getRank()) {  
      throw DataException("Error - slice size does not match Data rank.");  
      return Data();  
    }  
    //  
    // Create a new Data which is a slice of this one  
    return getSlice(slice_region);  
1459  }  }
1460    
1461  void  void
1462  Data::setItem(const boost::python::object& key,  Data::setItemO(const boost::python::object& key,
1463                const Data& value)                 const boost::python::object& value)
1464    {
1465      Data tempData(value,getFunctionSpace());
1466      setItemD(key,tempData);
1467    }
1468    
1469    void
1470    Data::setItemD(const boost::python::object& key,
1471                   const Data& value)
1472  {  {
1473    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1474    
1475    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1476    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1477      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
1478    }    }
1479    typeMatch(value);    if (getFunctionSpace()!=value.getFunctionSpace()) {
1480    setSlice(value,slice_region);       setSlice(Data(value,getFunctionSpace()),slice_region);
1481      } else {
1482         setSlice(value,slice_region);
1483      }
1484  }  }
1485    
1486  void  void
1487  Data::typeMatch(const Data& right)  Data::setSlice(const Data& value,
1488                   const DataArrayView::RegionType& region)
1489    {
1490      profData->slicing++;
1491      Data tempValue(value);
1492      typeMatchLeft(tempValue);
1493      typeMatchRight(tempValue);
1494      m_data->setSlice(tempValue.m_data.get(),region);
1495    }
1496    
1497    void
1498    Data::typeMatchLeft(Data& right) const
1499    {
1500      if (isExpanded()){
1501        right.expand();
1502      } else if (isTagged()) {
1503        if (right.isConstant()) {
1504          right.tag();
1505        }
1506      }
1507    }
1508    
1509    void
1510    Data::typeMatchRight(const Data& right)
1511  {  {
   //  
   // match the type of this to the RHS  
1512    if (isTagged()) {    if (isTagged()) {
1513      if (right.isExpanded()) {      if (right.isExpanded()) {
       //  
       // if the right hand side is expanded so must this  
1514        expand();        expand();
1515      }      }
1516    } else if (isConstant()) {    } else if (isConstant()) {
1517      if (right.isExpanded()) {      if (right.isExpanded()) {
       //  
       // if the right hand side is expanded so must this  
1518        expand();        expand();
1519      } else if (right.isTagged()) {      } else if (right.isTagged()) {
       //  
       // if the right hand side is tagged so must this  
1520        tag();        tag();
1521      }      }
1522    }    }
# Line 881  Data::setTaggedValue(int tagKey, Line 1543  Data::setTaggedValue(int tagKey,
1543    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
1544  }  }
1545    
 /*  
 Note: this version removed for now. Not needed, and breaks escript.cpp  
1546  void  void
1547  Data::setTaggedValue(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
1548                       const DataArrayView& value)                              const DataArrayView& value)
1549  {  {
1550    //    //
1551    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
# Line 899  Data::setTaggedValue(int tagKey, Line 1559  Data::setTaggedValue(int tagKey,
1559    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
1560    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
1561  }  }
1562  */  
1563    void
1564    Data::setRefValue(int ref,
1565                      const boost::python::numeric::array& value)
1566    {
1567      //
1568      // Construct DataArray from boost::python::object input value
1569      DataArray valueDataArray(value);
1570    
1571      //
1572      // Call DataAbstract::setRefValue
1573      m_data->setRefValue(ref,valueDataArray);
1574    }
1575    
1576    void
1577    Data::getRefValue(int ref,
1578                      boost::python::numeric::array& value)
1579    {
1580      //
1581      // Construct DataArray for boost::python::object return value
1582      DataArray valueDataArray(value);
1583    
1584      //
1585      // Load DataArray with values from data-points specified by ref
1586      m_data->getRefValue(ref,valueDataArray);
1587    
1588      //
1589      // Load values from valueDataArray into return numarray
1590    
1591      // extract the shape of the numarray
1592      int rank = value.getrank();
1593      DataArrayView::ShapeType shape;
1594      for (int i=0; i < rank; i++) {
1595        shape.push_back(extract<int>(value.getshape()[i]));
1596      }
1597    
1598      // and load the numarray with the data from the DataArray
1599      DataArrayView valueView = valueDataArray.getView();
1600    
1601      if (rank==0) {
1602          boost::python::numeric::array temp_numArray(valueView());
1603          value = temp_numArray;
1604      }
1605      if (rank==1) {
1606        for (int i=0; i < shape[0]; i++) {
1607          value[i] = valueView(i);
1608        }
1609      }
1610      if (rank==2) {
1611        for (int i=0; i < shape[0]; i++) {
1612          for (int j=0; j < shape[1]; j++) {
1613            value[i][j] = valueView(i,j);
1614          }
1615        }
1616      }
1617      if (rank==3) {
1618        for (int i=0; i < shape[0]; i++) {
1619          for (int j=0; j < shape[1]; j++) {
1620            for (int k=0; k < shape[2]; k++) {
1621              value[i][j][k] = valueView(i,j,k);
1622            }
1623          }
1624        }
1625      }
1626      if (rank==4) {
1627        for (int i=0; i < shape[0]; i++) {
1628          for (int j=0; j < shape[1]; j++) {
1629            for (int k=0; k < shape[2]; k++) {
1630              for (int l=0; l < shape[3]; l++) {
1631                value[i][j][k][l] = valueView(i,j,k,l);
1632              }
1633            }
1634          }
1635        }
1636      }
1637    
1638    }
1639    
1640    void
1641    Data::archiveData(const std::string fileName)
1642    {
1643      cout << "Archiving Data object to: " << fileName << endl;
1644    
1645      //
1646      // Determine type of this Data object
1647      int dataType = -1;
1648    
1649      if (isEmpty()) {
1650        dataType = 0;
1651        cout << "\tdataType: DataEmpty" << endl;
1652      }
1653      if (isConstant()) {
1654        dataType = 1;
1655        cout << "\tdataType: DataConstant" << endl;
1656      }
1657      if (isTagged()) {
1658        dataType = 2;
1659        cout << "\tdataType: DataTagged" << endl;
1660      }
1661      if (isExpanded()) {
1662        dataType = 3;
1663        cout << "\tdataType: DataExpanded" << endl;
1664      }
1665    
1666      if (dataType == -1) {
1667        throw DataException("archiveData Error: undefined dataType");
1668      }
1669    
1670      //
1671      // Collect data items common to all Data types
1672      int noSamples = getNumSamples();
1673      int noDPPSample = getNumDataPointsPerSample();
1674      int functionSpaceType = getFunctionSpace().getTypeCode();
1675      int dataPointRank = getDataPointRank();
1676      int dataPointSize = getDataPointSize();
1677      int dataLength = getLength();
1678      DataArrayView::ShapeType dataPointShape = getDataPointShape();
1679      int referenceNumbers[noSamples];
1680      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1681        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1682      }
1683      int tagNumbers[noSamples];
1684      if (isTagged()) {
1685        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1686          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1687        }
1688      }
1689    
1690      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1691      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1692      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1693    
1694      //
1695      // Flatten Shape to an array of integers suitable for writing to file
1696      int flatShape[4] = {0,0,0,0};
1697      cout << "\tshape: < ";
1698      for (int dim=0; dim<dataPointRank; dim++) {
1699        flatShape[dim] = dataPointShape[dim];
1700        cout << dataPointShape[dim] << " ";
1701      }
1702      cout << ">" << endl;
1703    
1704      //
1705      // Open archive file
1706      ofstream archiveFile;
1707      archiveFile.open(fileName.data(), ios::out);
1708    
1709      if (!archiveFile.good()) {
1710        throw DataException("archiveData Error: problem opening archive file");
1711      }
1712    
1713      //
1714      // Write common data items to archive file
1715      archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1716      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1717      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1718      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1719      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1720      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1721      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1722      for (int dim = 0; dim < 4; dim++) {
1723        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1724      }
1725      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1726        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1727      }
1728      if (isTagged()) {
1729        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1730          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1731        }
1732      }
1733    
1734      if (!archiveFile.good()) {
1735        throw DataException("archiveData Error: problem writing to archive file");
1736      }
1737    
1738      //
1739      // Archive underlying data values for each Data type
1740      int noValues;
1741      switch (dataType) {
1742        case 0:
1743          // DataEmpty
1744          noValues = 0;
1745          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1746          cout << "\tnoValues: " << noValues << endl;
1747          break;
1748        case 1:
1749          // DataConstant
1750          noValues = m_data->getLength();
1751          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1752          cout << "\tnoValues: " << noValues << endl;
1753          if (m_data->archiveData(archiveFile,noValues)) {
1754            throw DataException("archiveData Error: problem writing data to archive file");
1755          }
1756          break;
1757        case 2:
1758          // DataTagged
1759          noValues = m_data->getLength();
1760          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1761          cout << "\tnoValues: " << noValues << endl;
1762          if (m_data->archiveData(archiveFile,noValues)) {
1763            throw DataException("archiveData Error: problem writing data to archive file");
1764          }
1765          break;
1766        case 3:
1767          // DataExpanded
1768          noValues = m_data->getLength();
1769          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1770          cout << "\tnoValues: " << noValues << endl;
1771          if (m_data->archiveData(archiveFile,noValues)) {
1772            throw DataException("archiveData Error: problem writing data to archive file");
1773          }
1774          break;
1775      }
1776    
1777      if (!archiveFile.good()) {
1778        throw DataException("archiveData Error: problem writing data to archive file");
1779      }
1780    
1781      //
1782      // Close archive file
1783      archiveFile.close();
1784    
1785      if (!archiveFile.good()) {
1786        throw DataException("archiveData Error: problem closing archive file");
1787      }
1788    
1789    }
1790    
1791    void
1792    Data::extractData(const std::string fileName,
1793                      const FunctionSpace& fspace)
1794    {
1795      //
1796      // Can only extract Data to an object which is initially DataEmpty
1797      if (!isEmpty()) {
1798        throw DataException("extractData Error: can only extract to DataEmpty object");
1799      }
1800    
1801      cout << "Extracting Data object from: " << fileName << endl;
1802    
1803      int dataType;
1804      int noSamples;
1805      int noDPPSample;
1806      int functionSpaceType;
1807      int dataPointRank;
1808      int dataPointSize;
1809      int dataLength;
1810      DataArrayView::ShapeType dataPointShape;
1811      int flatShape[4];
1812    
1813      //
1814      // Open the archive file
1815      ifstream archiveFile;
1816      archiveFile.open(fileName.data(), ios::in);
1817    
1818      if (!archiveFile.good()) {
1819        throw DataException("extractData Error: problem opening archive file");
1820      }
1821    
1822      //
1823      // Read common data items from archive file
1824      archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1825      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1826      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1827      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1828      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1829      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1830      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1831      for (int dim = 0; dim < 4; dim++) {
1832        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1833        if (flatShape[dim]>0) {
1834          dataPointShape.push_back(flatShape[dim]);
1835        }
1836      }
1837      int referenceNumbers[noSamples];
1838      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1839        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1840      }
1841      int tagNumbers[noSamples];
1842      if (dataType==2) {
1843        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1844          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1845        }
1846      }
1847    
1848      if (!archiveFile.good()) {
1849        throw DataException("extractData Error: problem reading from archive file");
1850      }
1851    
1852      //
1853      // Verify the values just read from the archive file
1854      switch (dataType) {
1855        case 0:
1856          cout << "\tdataType: DataEmpty" << endl;
1857          break;
1858        case 1:
1859          cout << "\tdataType: DataConstant" << endl;
1860          break;
1861        case 2:
1862          cout << "\tdataType: DataTagged" << endl;
1863          break;
1864        case 3:
1865          cout << "\tdataType: DataExpanded" << endl;
1866          break;
1867        default:
1868          throw DataException("extractData Error: undefined dataType read from archive file");
1869          break;
1870      }
1871    
1872      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1873      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1874      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1875      cout << "\tshape: < ";
1876      for (int dim = 0; dim < dataPointRank; dim++) {
1877        cout << dataPointShape[dim] << " ";
1878      }
1879      cout << ">" << endl;
1880    
1881      //
1882      // Verify that supplied FunctionSpace object is compatible with this Data object.
1883      if ( (fspace.getTypeCode()!=functionSpaceType) ||
1884           (fspace.getNumSamples()!=noSamples) ||
1885           (fspace.getNumDPPSample()!=noDPPSample)
1886         ) {
1887        throw DataException("extractData Error: incompatible FunctionSpace");
1888      }
1889      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1890        if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
1891          throw DataException("extractData Error: incompatible FunctionSpace");
1892        }
1893      }
1894      if (dataType==2) {
1895        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1896          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
1897            throw DataException("extractData Error: incompatible FunctionSpace");
1898          }
1899        }
1900      }
1901    
1902      //
1903      // Construct a DataVector to hold underlying data values
1904      DataVector dataVec(dataLength);
1905    
1906      //
1907      // Load this DataVector with the appropriate values
1908      int noValues;
1909      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
1910      cout << "\tnoValues: " << noValues << endl;
1911      switch (dataType) {
1912        case 0:
1913          // DataEmpty
1914          if (noValues != 0) {
1915            throw DataException("extractData Error: problem reading data from archive file");
1916          }
1917          break;
1918        case 1:
1919          // DataConstant
1920          if (dataVec.extractData(archiveFile,noValues)) {
1921            throw DataException("extractData Error: problem reading data from archive file");
1922          }
1923          break;
1924        case 2:
1925          // DataTagged
1926          if (dataVec.extractData(archiveFile,noValues)) {
1927            throw DataException("extractData Error: problem reading data from archive file");
1928          }
1929          break;
1930        case 3:
1931          // DataExpanded
1932          if (dataVec.extractData(archiveFile,noValues)) {
1933            throw DataException("extractData Error: problem reading data from archive file");
1934          }
1935          break;
1936      }
1937    
1938      if (!archiveFile.good()) {
1939        throw DataException("extractData Error: problem reading from archive file");
1940      }
1941    
1942      //
1943      // Close archive file
1944      archiveFile.close();
1945    
1946      if (!archiveFile.good()) {
1947        throw DataException("extractData Error: problem closing archive file");
1948      }
1949    
1950      //
1951      // Construct an appropriate Data object
1952      DataAbstract* tempData;
1953      switch (dataType) {
1954        case 0:
1955          // DataEmpty
1956          tempData=new DataEmpty();
1957          break;
1958        case 1:
1959          // DataConstant
1960          tempData=new DataConstant(fspace,dataPointShape,dataVec);
1961          break;
1962        case 2:
1963          // DataTagged
1964          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
1965          break;
1966        case 3:
1967          // DataExpanded
1968          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
1969          break;
1970      }
1971      shared_ptr<DataAbstract> temp_data(tempData);
1972      m_data=temp_data;
1973    }
1974    
1975  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
1976  {  {

Legend:
Removed from v.100  
changed lines
  Added in v.126

  ViewVC Help
Powered by ViewVC 1.1.26