/[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 97 by jgs, Tue Dec 14 05:39:33 2004 UTC revision 153 by jgs, Tue Oct 25 01:51:20 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>
25  #include <functional>  #include <functional>
26  #include <math.h>  #include <math.h>
27    
28    #include <boost/python/dict.hpp>
29  #include <boost/python/str.hpp>  #include <boost/python/str.hpp>
30  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
31  #include <boost/python/long.hpp>  #include <boost/python/long.hpp>
32    
33  #include "escript/Data/DataException.h"  #include "escript/Data/DataException.h"
   
34  #include "escript/Data/DataExpanded.h"  #include "escript/Data/DataExpanded.h"
35  #include "escript/Data/DataConstant.h"  #include "escript/Data/DataConstant.h"
36  #include "escript/Data/DataTagged.h"  #include "escript/Data/DataTagged.h"
37  #include "escript/Data/DataEmpty.h"  #include "escript/Data/DataEmpty.h"
38  #include "escript/Data/DataArray.h"  #include "escript/Data/DataArray.h"
39  #include "escript/Data/DataAlgorithm.h"  #include "escript/Data/DataProf.h"
40  #include "escript/Data/FunctionSpaceFactory.h"  #include "escript/Data/FunctionSpaceFactory.h"
41  #include "escript/Data/AbstractContinuousDomain.h"  #include "escript/Data/AbstractContinuousDomain.h"
42  #include "escript/Data/UnaryFuncs.h"  #include "escript/Data/UnaryFuncs.h"
# Line 45  using namespace boost::python; Line 46  using namespace boost::python;
46  using namespace boost;  using namespace boost;
47  using namespace escript;  using namespace escript;
48    
49    #if defined DOPROF
50    //
51    // global table of profiling data for all Data objects
52    DataProf dataProfTable;
53    #endif
54    
55  Data::Data()  Data::Data()
56  {  {
57    //    //
# Line 52  Data::Data() Line 59  Data::Data()
59    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
60    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
61    m_data=temp_data;    m_data=temp_data;
62    #if defined DOPROF
63      // create entry in global profiling table for this object
64      profData = dataProfTable.newData();
65    #endif
66  }  }
67    
68  Data::Data(double value,  Data::Data(double value,
# Line 65  Data::Data(double value, Line 76  Data::Data(double value,
76    }    }
77    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
78    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
79    #if defined DOPROF
80      // create entry in global profiling table for this object
81      profData = dataProfTable.newData();
82    #endif
83  }  }
84    
85  Data::Data(double value,  Data::Data(double value,
# Line 75  Data::Data(double value, Line 90  Data::Data(double value,
90    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
91    pair<int,int> dataShape=what.getDataShape();    pair<int,int> dataShape=what.getDataShape();
92    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
93    #if defined DOPROF
94      // create entry in global profiling table for this object
95      profData = dataProfTable.newData();
96    #endif
97  }  }
98    
99  Data::Data(const Data& inData)  Data::Data(const Data& inData)
100  {  {
101    m_data=inData.m_data;    m_data=inData.m_data;
102    #if defined DOPROF
103      // create entry in global profiling table for this object
104      profData = dataProfTable.newData();
105    #endif
106  }  }
107    
108  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 90  Data::Data(const Data& inData, Line 113  Data::Data(const Data& inData,
113    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
114    shared_ptr<DataAbstract> temp_data(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
115    m_data=temp_data;    m_data=temp_data;
116    #if defined DOPROF
117      // create entry in global profiling table for this object
118      profData = dataProfTable.newData();
119    #endif
120  }  }
121    
122  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 99  Data::Data(const Data& inData, Line 126  Data::Data(const Data& inData,
126      m_data=inData.m_data;      m_data=inData.m_data;
127    } else {    } else {
128      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
129      // Note for Lutz, Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
130      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
131      // be able to create an instance of AbstractDomain but that was done      // be able to create an instance of AbstractDomain but that was done
132      // as a boost python work around which may no longer be required.      // as a boost:python work around which may no longer be required.
133      const AbstractDomain& inDataDomain=inData.getDomain();      const AbstractDomain& inDataDomain=inData.getDomain();
134      if  (inDataDomain==functionspace.getDomain()) {      if  (inDataDomain==functionspace.getDomain()) {
135        inDataDomain.interpolateOnDomain(tmp,inData);        inDataDomain.interpolateOnDomain(tmp,inData);
# Line 111  Data::Data(const Data& inData, Line 138  Data::Data(const Data& inData,
138      }      }
139      m_data=tmp.m_data;      m_data=tmp.m_data;
140    }    }
141    #if defined DOPROF
142      // create entry in global profiling table for this object
143      profData = dataProfTable.newData();
144    #endif
145  }  }
146    
147  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(const DataTagged::TagListType& tagKeys,
# Line 125  Data::Data(const DataTagged::TagListType Line 156  Data::Data(const DataTagged::TagListType
156    if (expanded) {    if (expanded) {
157      expand();      expand();
158    }    }
159    #if defined DOPROF
160      // create entry in global profiling table for this object
161      profData = dataProfTable.newData();
162    #endif
163  }  }
164    
165  Data::Data(const numeric::array& value,  Data::Data(const numeric::array& value,
# Line 132  Data::Data(const numeric::array& value, Line 167  Data::Data(const numeric::array& value,
167             bool expanded)             bool expanded)
168  {  {
169    initialise(value,what,expanded);    initialise(value,what,expanded);
170    #if defined DOPROF
171      // create entry in global profiling table for this object
172      profData = dataProfTable.newData();
173    #endif
174  }  }
175    
176  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
# Line 139  Data::Data(const DataArrayView& value, Line 178  Data::Data(const DataArrayView& value,
178             bool expanded)             bool expanded)
179  {  {
180    initialise(value,what,expanded);    initialise(value,what,expanded);
181    #if defined DOPROF
182      // create entry in global profiling table for this object
183      profData = dataProfTable.newData();
184    #endif
185  }  }
186    
187  Data::Data(const object& value,  Data::Data(const object& value,
# Line 147  Data::Data(const object& value, Line 190  Data::Data(const object& value,
190  {  {
191    numeric::array asNumArray(value);    numeric::array asNumArray(value);
192    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
193    #if defined DOPROF
194      // create entry in global profiling table for this object
195      profData = dataProfTable.newData();
196    #endif
197  }  }
198    
199  Data::Data(const object& value,  Data::Data(const object& value,
# Line 167  Data::Data(const object& value, Line 214  Data::Data(const object& value,
214      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
215      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
216    }    }
217    #if defined DOPROF
218      // create entry in global profiling table for this object
219      profData = dataProfTable.newData();
220    #endif
221    }
222    
223    Data::~Data()
224    {
225    
226  }  }
227    
228  escriptDataC  escriptDataC
# Line 185  Data::getDataC() const Line 241  Data::getDataC() const
241    return temp;    return temp;
242  }  }
243    
244  tuple  const boost::python::tuple
245  Data::getShapeTuple() const  Data::getShapeTuple() const
246  {  {
247    const DataArrayView::ShapeType& shape=getDataPointShape();    const DataArrayView::ShapeType& shape=getDataPointShape();
# Line 353  Data::reshapeDataPoint(const DataArrayVi Line 409  Data::reshapeDataPoint(const DataArrayVi
409  Data  Data
410  Data::wherePositive() const  Data::wherePositive() const
411  {  {
412    #if defined DOPROF
413      profData->where++;
414    #endif
415    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
416  }  }
417    
418  Data  Data
419  Data::whereNegative() const  Data::whereNegative() const
420  {  {
421    #if defined DOPROF
422      profData->where++;
423    #endif
424    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
425  }  }
426    
427  Data  Data
428  Data::whereNonNegative() const  Data::whereNonNegative() const
429  {  {
430    #if defined DOPROF
431      profData->where++;
432    #endif
433    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
434  }  }
435    
436  Data  Data
437  Data::whereNonPositive() const  Data::whereNonPositive() const
438  {  {
439    #if defined DOPROF
440      profData->where++;
441    #endif
442    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
443  }  }
444    
445  Data  Data
446  Data::whereZero() const  Data::whereZero() const
447  {  {
448    #if defined DOPROF
449      profData->where++;
450    #endif
451    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));
452  }  }
453    
454  Data  Data
455  Data::whereNonZero() const  Data::whereNonZero() const
456  {  {
457    #if defined DOPROF
458      profData->where++;
459    #endif
460    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));
461  }  }
462    
463  Data  Data
464  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
465  {  {
466    #if defined DOPROF
467      profData->interpolate++;
468    #endif
469    return Data(*this,functionspace);    return Data(*this,functionspace);
470  }  }
471    
# Line 410  Data::probeInterpolation(const FunctionS Line 487  Data::probeInterpolation(const FunctionS
487  Data  Data
488  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
489  {  {
490    #if defined DOPROF
491      profData->grad++;
492    #endif
493    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
494      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
495    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
# Line 443  Data::getDataPointShape() const Line 523  Data::getDataPointShape() const
523    return getPointDataView().getShape();    return getPointDataView().getShape();
524  }  }
525    
526    void
527    Data::fillFromNumArray(const boost::python::numeric::array num_array)
528    {
529      //
530      // check rank
531      if (num_array.getrank()<getDataPointRank())
532          throw DataException("Rank of numarray does not match Data object rank");
533    
534      //
535      // check shape of num_array
536      for (int i=0; i<getDataPointRank(); i++) {
537        if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])
538           throw DataException("Shape of numarray does not match Data object rank");
539      }
540    
541      //
542      // make sure data is expanded:
543      if (!isExpanded()) {
544        expand();
545      }
546    
547      //
548      // and copy over
549      m_data->copyAll(num_array);
550    }
551    
552    const
553    boost::python::numeric::array
554    Data::convertToNumArray()
555    {
556      //
557      // determine the total number of data points
558      int numSamples = getNumSamples();
559      int numDataPointsPerSample = getNumDataPointsPerSample();
560      int numDataPoints = numSamples * numDataPointsPerSample;
561    
562      //
563      // determine the rank and shape of each data point
564      int dataPointRank = getDataPointRank();
565      DataArrayView::ShapeType dataPointShape = getDataPointShape();
566    
567      //
568      // create the numeric array to be returned
569      boost::python::numeric::array numArray(0.0);
570    
571      //
572      // the rank of the returned numeric array will be the rank of
573      // the data points, plus one. Where the rank of the array is n,
574      // the last n-1 dimensions will be equal to the shape of the
575      // data points, whilst the first dimension will be equal to the
576      // total number of data points. Thus the array will consist of
577      // a serial vector of the data points.
578      int arrayRank = dataPointRank + 1;
579      DataArrayView::ShapeType arrayShape;
580      arrayShape.push_back(numDataPoints);
581      for (int d=0; d<dataPointRank; d++) {
582         arrayShape.push_back(dataPointShape[d]);
583      }
584    
585      //
586      // resize the numeric array to the shape just calculated
587      if (arrayRank==1) {
588        numArray.resize(arrayShape[0]);
589      }
590      if (arrayRank==2) {
591        numArray.resize(arrayShape[0],arrayShape[1]);
592      }
593      if (arrayRank==3) {
594        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
595      }
596      if (arrayRank==4) {
597        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
598      }
599      if (arrayRank==5) {
600        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
601      }
602    
603      //
604      // loop through each data point in turn, loading the values for that data point
605      // into the numeric array.
606      int dataPoint = 0;
607      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
608        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
609          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
610          if (dataPointRank==0) {
611            numArray[dataPoint]=dataPointView();
612          }
613          if (dataPointRank==1) {
614            for (int i=0; i<dataPointShape[0]; i++) {
615              numArray[dataPoint][i]=dataPointView(i);
616            }
617          }
618          if (dataPointRank==2) {
619            for (int i=0; i<dataPointShape[0]; i++) {
620              for (int j=0; j<dataPointShape[1]; j++) {
621                numArray[dataPoint][i][j] = dataPointView(i,j);
622              }
623            }
624          }
625          if (dataPointRank==3) {
626            for (int i=0; i<dataPointShape[0]; i++) {
627              for (int j=0; j<dataPointShape[1]; j++) {
628                for (int k=0; k<dataPointShape[2]; k++) {
629                  numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
630                }
631              }
632            }
633          }
634          if (dataPointRank==4) {
635            for (int i=0; i<dataPointShape[0]; i++) {
636              for (int j=0; j<dataPointShape[1]; j++) {
637                for (int k=0; k<dataPointShape[2]; k++) {
638                  for (int l=0; l<dataPointShape[3]; l++) {
639                    numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
640                  }
641                }
642              }
643            }
644          }
645          dataPoint++;
646        }
647      }
648    
649      //
650      // return the loaded array
651      return numArray;
652    }
653    
654    const
655    boost::python::numeric::array
656    Data::convertToNumArrayFromSampleNo(int sampleNo)
657    {
658      //
659      // Check a valid sample number has been supplied
660      if (sampleNo >= getNumSamples()) {
661        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
662      }
663    
664      //
665      // determine the number of data points per sample
666      int numDataPointsPerSample = getNumDataPointsPerSample();
667    
668      //
669      // determine the rank and shape of each data point
670      int dataPointRank = getDataPointRank();
671      DataArrayView::ShapeType dataPointShape = getDataPointShape();
672    
673      //
674      // create the numeric array to be returned
675      boost::python::numeric::array numArray(0.0);
676    
677      //
678      // the rank of the returned numeric array will be the rank of
679      // the data points, plus one. Where the rank of the array is n,
680      // the last n-1 dimensions will be equal to the shape of the
681      // data points, whilst the first dimension will be equal to the
682      // total number of data points. Thus the array will consist of
683      // a serial vector of the data points.
684      int arrayRank = dataPointRank + 1;
685      DataArrayView::ShapeType arrayShape;
686      arrayShape.push_back(numDataPointsPerSample);
687      for (int d=0; d<dataPointRank; d++) {
688         arrayShape.push_back(dataPointShape[d]);
689      }
690    
691      //
692      // resize the numeric array to the shape just calculated
693      if (arrayRank==1) {
694        numArray.resize(arrayShape[0]);
695      }
696      if (arrayRank==2) {
697        numArray.resize(arrayShape[0],arrayShape[1]);
698      }
699      if (arrayRank==3) {
700        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
701      }
702      if (arrayRank==4) {
703        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
704      }
705      if (arrayRank==5) {
706        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
707      }
708    
709      //
710      // loop through each data point in turn, loading the values for that data point
711      // into the numeric array.
712      for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {
713        DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);
714        if (dataPointRank==0) {
715          numArray[dataPoint]=dataPointView();
716        }
717        if (dataPointRank==1) {
718          for (int i=0; i<dataPointShape[0]; i++) {
719            numArray[dataPoint][i]=dataPointView(i);
720          }
721        }
722        if (dataPointRank==2) {
723          for (int i=0; i<dataPointShape[0]; i++) {
724            for (int j=0; j<dataPointShape[1]; j++) {
725              numArray[dataPoint][i][j] = dataPointView(i,j);
726            }
727          }
728        }
729        if (dataPointRank==3) {
730          for (int i=0; i<dataPointShape[0]; i++) {
731            for (int j=0; j<dataPointShape[1]; j++) {
732              for (int k=0; k<dataPointShape[2]; k++) {
733                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
734              }
735            }
736          }
737        }
738        if (dataPointRank==4) {
739          for (int i=0; i<dataPointShape[0]; i++) {
740            for (int j=0; j<dataPointShape[1]; j++) {
741              for (int k=0; k<dataPointShape[2]; k++) {
742                for (int l=0; l<dataPointShape[3]; l++) {
743                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
744                }
745              }
746            }
747          }
748        }
749      }
750    
751      //
752      // return the loaded array
753      return numArray;
754    }
755    
756    const
757    boost::python::numeric::array
758    Data::convertToNumArrayFromDPNo(int sampleNo,
759                                    int dataPointNo)
760    {
761      //
762      // Check a valid sample number has been supplied
763      if (sampleNo >= getNumSamples()) {
764        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
765      }
766    
767      //
768      // Check a valid data point number has been supplied
769      if (dataPointNo >= getNumDataPointsPerSample()) {
770        throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");
771      }
772    
773      //
774      // determine the rank and shape of each data point
775      int dataPointRank = getDataPointRank();
776      DataArrayView::ShapeType dataPointShape = getDataPointShape();
777    
778      //
779      // create the numeric array to be returned
780      boost::python::numeric::array numArray(0.0);
781    
782      //
783      // the shape of the returned numeric array will be the same
784      // as that of the data point
785      int arrayRank = dataPointRank;
786      DataArrayView::ShapeType arrayShape = dataPointShape;
787    
788      //
789      // resize the numeric array to the shape just calculated
790      if (arrayRank==0) {
791        numArray.resize(1);
792      }
793      if (arrayRank==1) {
794        numArray.resize(arrayShape[0]);
795      }
796      if (arrayRank==2) {
797        numArray.resize(arrayShape[0],arrayShape[1]);
798      }
799      if (arrayRank==3) {
800        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
801      }
802      if (arrayRank==4) {
803        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
804      }
805    
806      //
807      // load the values for the data point into the numeric array.
808      DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
809      if (dataPointRank==0) {
810        numArray[0]=dataPointView();
811      }
812      if (dataPointRank==1) {
813        for (int i=0; i<dataPointShape[0]; i++) {
814          numArray[i]=dataPointView(i);
815        }
816      }
817      if (dataPointRank==2) {
818        for (int i=0; i<dataPointShape[0]; i++) {
819          for (int j=0; j<dataPointShape[1]; j++) {
820            numArray[i][j] = dataPointView(i,j);
821          }
822        }
823      }
824      if (dataPointRank==3) {
825        for (int i=0; i<dataPointShape[0]; i++) {
826          for (int j=0; j<dataPointShape[1]; j++) {
827            for (int k=0; k<dataPointShape[2]; k++) {
828              numArray[i][j][k]=dataPointView(i,j,k);
829            }
830          }
831        }
832      }
833      if (dataPointRank==4) {
834        for (int i=0; i<dataPointShape[0]; i++) {
835          for (int j=0; j<dataPointShape[1]; j++) {
836            for (int k=0; k<dataPointShape[2]; k++) {
837              for (int l=0; l<dataPointShape[3]; l++) {
838                numArray[i][j][k][l]=dataPointView(i,j,k,l);
839              }
840            }
841          }
842        }
843      }
844    
845      //
846      // return the loaded array
847      return numArray;
848    }
849    
850  boost::python::numeric::array  boost::python::numeric::array
851  Data::integrate() const  Data::integrate() const
852  {  {
# Line 450  Data::integrate() const Line 854  Data::integrate() const
854    int rank = getDataPointRank();    int rank = getDataPointRank();
855    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
856    
857    #if defined DOPROF
858      profData->integrate++;
859    #endif
860    
861    //    //
862    // calculate the integral values    // calculate the integral values
863    vector<double> integrals(getDataPointSize());    vector<double> integrals(getDataPointSize());
# Line 460  Data::integrate() const Line 868  Data::integrate() const
868    // and load the array with the integral values    // and load the array with the integral values
869    boost::python::numeric::array bp_array(1.0);    boost::python::numeric::array bp_array(1.0);
870    if (rank==0) {    if (rank==0) {
871        bp_array.resize(1);
872      index = 0;      index = 0;
873      bp_array[0] = integrals[index];      bp_array[0] = integrals[index];
874    }    }
# Line 512  Data::integrate() const Line 921  Data::integrate() const
921  Data  Data
922  Data::sin() const  Data::sin() const
923  {  {
924    #if defined DOPROF
925      profData->unary++;
926    #endif
927    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
928  }  }
929    
930  Data  Data
931  Data::cos() const  Data::cos() const
932  {  {
933    #if defined DOPROF
934      profData->unary++;
935    #endif
936    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
937  }  }
938    
939  Data  Data
940  Data::tan() const  Data::tan() const
941  {  {
942    #if defined DOPROF
943      profData->unary++;
944    #endif
945    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
946  }  }
947    
948  Data  Data
949    Data::asin() const
950    {
951    #if defined DOPROF
952      profData->unary++;
953    #endif
954      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
955    }
956    
957    Data
958    Data::acos() const
959    {
960    #if defined DOPROF
961      profData->unary++;
962    #endif
963      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
964    }
965    
966    Data
967    Data::atan() const
968    {
969    #if defined DOPROF
970      profData->unary++;
971    #endif
972      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
973    }
974    
975    Data
976    Data::sinh() const
977    {
978    #if defined DOPROF
979      profData->unary++;
980    #endif
981      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
982    }
983    
984    Data
985    Data::cosh() const
986    {
987    #if defined DOPROF
988      profData->unary++;
989    #endif
990      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
991    }
992    
993    Data
994    Data::tanh() const
995    {
996    #if defined DOPROF
997      profData->unary++;
998    #endif
999      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
1000    }
1001    
1002    Data
1003    Data::asinh() const
1004    {
1005    #if defined DOPROF
1006      profData->unary++;
1007    #endif
1008      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
1009    }
1010    
1011    Data
1012    Data::acosh() const
1013    {
1014    #if defined DOPROF
1015      profData->unary++;
1016    #endif
1017      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
1018    }
1019    
1020    Data
1021    Data::atanh() const
1022    {
1023    #if defined DOPROF
1024      profData->unary++;
1025    #endif
1026      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
1027    }
1028    
1029    Data
1030  Data::log() const  Data::log() const
1031  {  {
1032    #if defined DOPROF
1033      profData->unary++;
1034    #endif
1035    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
1036  }  }
1037    
1038  Data  Data
1039  Data::ln() const  Data::ln() const
1040  {  {
1041    #if defined DOPROF
1042      profData->unary++;
1043    #endif
1044    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
1045  }  }
1046    
1047    Data
1048    Data::sign() const
1049    {
1050    #if defined DOPROF
1051      profData->unary++;
1052    #endif
1053      return escript::unaryOp(*this,escript::fsign);
1054    }
1055    
1056    Data
1057    Data::abs() const
1058    {
1059    #if defined DOPROF
1060      profData->unary++;
1061    #endif
1062      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
1063    }
1064    
1065    Data
1066    Data::neg() const
1067    {
1068    #if defined DOPROF
1069      profData->unary++;
1070    #endif
1071      return escript::unaryOp(*this,negate<double>());
1072    }
1073    
1074    Data
1075    Data::pos() const
1076    {
1077    #if defined DOPROF
1078      profData->unary++;
1079    #endif
1080      Data result;
1081      // perform a deep copy
1082      result.copy(*this);
1083      return result;
1084    }
1085    
1086    Data
1087    Data::exp() const
1088    {
1089    #if defined DOPROF
1090      profData->unary++;
1091    #endif
1092      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1093    }
1094    
1095    Data
1096    Data::sqrt() const
1097    {
1098    #if defined DOPROF
1099      profData->unary++;
1100    #endif
1101      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1102    }
1103    
1104  double  double
1105  Data::Lsup() const  Data::Lsup() const
1106  {  {
1107    #if defined DOPROF
1108      profData->reduction1++;
1109    #endif
1110    //    //
1111    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
1112    return algorithm(DataAlgorithmAdapter<AbsMax>(0));    AbsMax abs_max_func;
1113      return algorithm(abs_max_func,0);
1114    }
1115    
1116    double
1117    Data::Linf() const
1118    {
1119    #if defined DOPROF
1120      profData->reduction1++;
1121    #endif
1122      //
1123      // set the initial absolute minimum value to max double
1124      AbsMin abs_min_func;
1125      return algorithm(abs_min_func,numeric_limits<double>::max());
1126  }  }
1127    
1128  double  double
1129  Data::sup() const  Data::sup() const
1130  {  {
1131    #if defined DOPROF
1132      profData->reduction1++;
1133    #endif
1134    //    //
1135    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1136    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::min()));    FMax fmax_func;
1137      return algorithm(fmax_func,numeric_limits<double>::max()*-1);
1138  }  }
1139    
1140  double  double
1141  Data::inf() const  Data::inf() const
1142  {  {
1143    #if defined DOPROF
1144      profData->reduction1++;
1145    #endif
1146    //    //
1147    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1148    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1149      return algorithm(fmin_func,numeric_limits<double>::max());
1150  }  }
1151    
1152  Data  Data
1153  Data::maxval() const  Data::maxval() const
1154  {  {
1155    // not implemented - will use dp_algorithm  #if defined DOPROF
1156    return (*this);    profData->reduction2++;
1157    #endif
1158      //
1159      // set the initial maximum value to min possible double
1160      FMax fmax_func;
1161      return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1162  }  }
1163    
1164  Data  Data
1165  Data::minval() const  Data::minval() const
1166  {  {
1167    // not implemented - will use dp_algorithm  #if defined DOPROF
1168    return (*this);    profData->reduction2++;
1169    #endif
1170      //
1171      // set the initial minimum value to max possible double
1172      FMin fmin_func;
1173      return dp_algorithm(fmin_func,numeric_limits<double>::max());
1174  }  }
1175    
1176  Data  Data
1177  Data::length() const  Data::length() const
1178  {  {
1179    // not implemented - will use dp_algorithm  #if defined DOPROF
1180    return (*this);    profData->reduction2++;
1181    #endif
1182      Length len_func;
1183      return dp_algorithm(len_func,0);
1184  }  }
1185    
1186  Data  Data
1187  Data::trace() const  Data::trace() const
1188  {  {
1189    // not implemented - will use dp_algorithm  #if defined DOPROF
1190    return (*this);    profData->reduction2++;
1191    #endif
1192      Trace trace_func;
1193      return dp_algorithm(trace_func,0);
1194  }  }
1195    
1196  Data  Data
1197  Data::transpose(int axis) const  Data::transpose(int axis) const
1198  {  {
1199    #if defined DOPROF
1200      profData->reduction2++;
1201    #endif
1202    // not implemented    // not implemented
1203    return (*this);    throw DataException("Error - Data::transpose not implemented yet.");
1204      return Data();
1205  }  }
1206    
1207  Data  const boost::python::tuple
1208  Data::sign() const  Data::mindp() const
1209  {  {
1210    return escript::unaryOp(*this,escript::fsign);    // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1211  }    // abort (for unknown reasons) if there are openmp directives with it in the
1212      // surrounding function
1213    
1214  Data    int SampleNo;
1215  Data::abs() const    int DataPointNo;
 {  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);  
 }  
1216    
1217  Data    calc_mindp(SampleNo,DataPointNo);
1218  Data::neg() const  
1219  {    return make_tuple(SampleNo,DataPointNo);
   return escript::unaryOp(*this,negate<double>());  
1220  }  }
1221    
1222  Data  void
1223  Data::pos() const  Data::calc_mindp(int& SampleNo,
1224                     int& DataPointNo) const
1225  {  {
1226    return (*this);    int i,j;
1227      int lowi=0,lowj=0;
1228      double min=numeric_limits<double>::max();
1229    
1230      Data temp=minval();
1231    
1232      int numSamples=temp.getNumSamples();
1233      int numDPPSample=temp.getNumDataPointsPerSample();
1234    
1235      double next,local_min;
1236      int local_lowi,local_lowj;
1237    
1238      #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1239      {
1240        local_min=min;
1241        #pragma omp for private(i,j) schedule(static)
1242        for (i=0; i<numSamples; i++) {
1243          for (j=0; j<numDPPSample; j++) {
1244            next=temp.getDataPoint(i,j)();
1245            if (next<local_min) {
1246              local_min=next;
1247              local_lowi=i;
1248              local_lowj=j;
1249            }
1250          }
1251        }
1252        #pragma omp critical
1253        if (local_min<min) {
1254          min=local_min;
1255          lowi=local_lowi;
1256          lowj=local_lowj;
1257        }
1258      }
1259    
1260      SampleNo = lowi;
1261      DataPointNo = lowj;
1262  }  }
1263    
1264  Data  void
1265  Data::exp() const  Data::saveDX(std::string fileName) const
1266  {  {
1267    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    boost::python::dict args;
1268      args["data"]=boost::python::object(this);
1269      getDomain().saveDX(fileName,args);
1270      return;
1271  }  }
1272    
1273  Data  void
1274  Data::sqrt() const  Data::saveVTK(std::string fileName) const
1275  {  {
1276    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    boost::python::dict args;
1277      args["data"]=boost::python::object(this);
1278      getDomain().saveVTK(fileName,args);
1279      return;
1280  }  }
1281    
1282  Data&  Data&
1283  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1284  {  {
1285    #if defined DOPROF
1286      profData->binary++;
1287    #endif
1288    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1289    return (*this);    return (*this);
1290  }  }
# Line 644  Data::operator+=(const Data& right) Line 1292  Data::operator+=(const Data& right)
1292  Data&  Data&
1293  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1294  {  {
1295    #if defined DOPROF
1296      profData->binary++;
1297    #endif
1298    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1299    return (*this);    return (*this);
1300  }  }
# Line 651  Data::operator+=(const boost::python::ob Line 1302  Data::operator+=(const boost::python::ob
1302  Data&  Data&
1303  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1304  {  {
1305    #if defined DOPROF
1306      profData->binary++;
1307    #endif
1308    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1309    return (*this);    return (*this);
1310  }  }
# Line 658  Data::operator-=(const Data& right) Line 1312  Data::operator-=(const Data& right)
1312  Data&  Data&
1313  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1314  {  {
1315    #if defined DOPROF
1316      profData->binary++;
1317    #endif
1318    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1319    return (*this);    return (*this);
1320  }  }
# Line 665  Data::operator-=(const boost::python::ob Line 1322  Data::operator-=(const boost::python::ob
1322  Data&  Data&
1323  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1324  {  {
1325    #if defined DOPROF
1326      profData->binary++;
1327    #endif
1328    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1329    return (*this);    return (*this);
1330  }  }
# Line 672  Data::operator*=(const Data& right) Line 1332  Data::operator*=(const Data& right)
1332  Data&  Data&
1333  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1334  {  {
1335    #if defined DOPROF
1336      profData->binary++;
1337    #endif
1338    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1339    return (*this);    return (*this);
1340  }  }
# Line 679  Data::operator*=(const boost::python::ob Line 1342  Data::operator*=(const boost::python::ob
1342  Data&  Data&
1343  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1344  {  {
1345    #if defined DOPROF
1346      profData->binary++;
1347    #endif
1348    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1349    return (*this);    return (*this);
1350  }  }
# Line 686  Data::operator/=(const Data& right) Line 1352  Data::operator/=(const Data& right)
1352  Data&  Data&
1353  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1354  {  {
1355    #if defined DOPROF
1356      profData->binary++;
1357    #endif
1358    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1359    return (*this);    return (*this);
1360  }  }
# Line 693  Data::operator/=(const boost::python::ob Line 1362  Data::operator/=(const boost::python::ob
1362  Data  Data
1363  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1364  {  {
1365    #if defined DOPROF
1366      profData->binary++;
1367    #endif
1368    Data result;    Data result;
1369    result.copy(*this);    result.copy(*this);
1370    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 702  Data::powO(const boost::python::object& Line 1374  Data::powO(const boost::python::object&
1374  Data  Data
1375  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1376  {  {
1377    #if defined DOPROF
1378      profData->binary++;
1379    #endif
1380    Data result;    Data result;
1381    result.copy(*this);    result.copy(*this);
1382    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 709  Data::powD(const Data& right) const Line 1384  Data::powD(const Data& right) const
1384  }  }
1385    
1386  //  //
1387  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1388  Data  Data
1389  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1390  {  {
# Line 722  escript::operator+(const Data& left, con Line 1397  escript::operator+(const Data& left, con
1397  }  }
1398    
1399  //  //
1400  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1401  Data  Data
1402  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1403  {  {
# Line 735  escript::operator-(const Data& left, con Line 1410  escript::operator-(const Data& left, con
1410  }  }
1411    
1412  //  //
1413  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1414  Data  Data
1415  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
1416  {  {
# Line 748  escript::operator*(const Data& left, con Line 1423  escript::operator*(const Data& left, con
1423  }  }
1424    
1425  //  //
1426  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1427  Data  Data
1428  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
1429  {  {
# Line 761  escript::operator/(const Data& left, con Line 1436  escript::operator/(const Data& left, con
1436  }  }
1437    
1438  //  //
1439  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1440  Data  Data
1441  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1442  {  {
# Line 777  escript::operator+(const Data& left, con Line 1452  escript::operator+(const Data& left, con
1452  }  }
1453    
1454  //  //
1455  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1456  Data  Data
1457  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1458  {  {
# Line 793  escript::operator-(const Data& left, con Line 1468  escript::operator-(const Data& left, con
1468  }  }
1469    
1470  //  //
1471  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1472  Data  Data
1473  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1474  {  {
# Line 809  escript::operator*(const Data& left, con Line 1484  escript::operator*(const Data& left, con
1484  }  }
1485    
1486  //  //
1487  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1488  Data  Data
1489  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1490  {  {
# Line 825  escript::operator/(const Data& left, con Line 1500  escript::operator/(const Data& left, con
1500  }  }
1501    
1502  //  //
1503  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1504  Data  Data
1505  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1506  {  {
# Line 838  escript::operator+(const boost::python:: Line 1513  escript::operator+(const boost::python::
1513  }  }
1514    
1515  //  //
1516  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1517  Data  Data
1518  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1519  {  {
# Line 851  escript::operator-(const boost::python:: Line 1526  escript::operator-(const boost::python::
1526  }  }
1527    
1528  //  //
1529  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1530  Data  Data
1531  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1532  {  {
# Line 864  escript::operator*(const boost::python:: Line 1539  escript::operator*(const boost::python::
1539  }  }
1540    
1541  //  //
1542  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1543  Data  Data
1544  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1545  {  {
# Line 877  escript::operator/(const boost::python:: Line 1552  escript::operator/(const boost::python::
1552  }  }
1553    
1554  //  //
 // NOTE: It is essential to specify the namepsace this operator belongs to  
1555  //bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1556  //{  //{
1557  //  /*  //  /*
# Line 939  Data::getItem(const boost::python::objec Line 1613  Data::getItem(const boost::python::objec
1613  Data  Data
1614  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1615  {  {
1616    #if defined DOPROF
1617      profData->slicing++;
1618    #endif
1619    return Data(*this,region);    return Data(*this,region);
1620  }  }
1621    
# Line 955  Data::setItemD(const boost::python::obje Line 1632  Data::setItemD(const boost::python::obje
1632                 const Data& value)                 const Data& value)
1633  {  {
1634    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1635    
1636    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1637    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1638      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
1639    }    }
1640    setSlice(value,slice_region);    if (getFunctionSpace()!=value.getFunctionSpace()) {
1641         setSlice(Data(value,getFunctionSpace()),slice_region);
1642      } else {
1643         setSlice(value,slice_region);
1644      }
1645  }  }
1646    
1647  void  void
1648  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1649                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
1650  {  {
1651    #if defined DOPROF
1652      profData->slicing++;
1653    #endif
1654    Data tempValue(value);    Data tempValue(value);
1655    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1656    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1021  Data::setTaggedValue(int tagKey, Line 1706  Data::setTaggedValue(int tagKey,
1706    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
1707  }  }
1708    
 /*  
 Note: this version removed for now. Not needed, and breaks escript.cpp  
1709  void  void
1710  Data::setTaggedValue(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
1711                       const DataArrayView& value)                              const DataArrayView& value)
1712  {  {
1713    //    //
1714    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
# Line 1039  Data::setTaggedValue(int tagKey, Line 1722  Data::setTaggedValue(int tagKey,
1722    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
1723    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
1724  }  }
1725  */  
1726    int
1727    Data::getTagNumber(int dpno)
1728    {
1729      return m_data->getTagNumber(dpno);
1730    }
1731    
1732    void
1733    Data::setRefValue(int ref,
1734                      const boost::python::numeric::array& value)
1735    {
1736      //
1737      // Construct DataArray from boost::python::object input value
1738      DataArray valueDataArray(value);
1739    
1740      //
1741      // Call DataAbstract::setRefValue
1742      m_data->setRefValue(ref,valueDataArray);
1743    }
1744    
1745    void
1746    Data::getRefValue(int ref,
1747                      boost::python::numeric::array& value)
1748    {
1749      //
1750      // Construct DataArray for boost::python::object return value
1751      DataArray valueDataArray(value);
1752    
1753      //
1754      // Load DataArray with values from data-points specified by ref
1755      m_data->getRefValue(ref,valueDataArray);
1756    
1757      //
1758      // Load values from valueDataArray into return numarray
1759    
1760      // extract the shape of the numarray
1761      int rank = value.getrank();
1762      DataArrayView::ShapeType shape;
1763      for (int i=0; i < rank; i++) {
1764        shape.push_back(extract<int>(value.getshape()[i]));
1765      }
1766    
1767      // and load the numarray with the data from the DataArray
1768      DataArrayView valueView = valueDataArray.getView();
1769    
1770      if (rank==0) {
1771          boost::python::numeric::array temp_numArray(valueView());
1772          value = temp_numArray;
1773      }
1774      if (rank==1) {
1775        for (int i=0; i < shape[0]; i++) {
1776          value[i] = valueView(i);
1777        }
1778      }
1779      if (rank==2) {
1780        for (int i=0; i < shape[0]; i++) {
1781          for (int j=0; j < shape[1]; j++) {
1782            value[i][j] = valueView(i,j);
1783          }
1784        }
1785      }
1786      if (rank==3) {
1787        for (int i=0; i < shape[0]; i++) {
1788          for (int j=0; j < shape[1]; j++) {
1789            for (int k=0; k < shape[2]; k++) {
1790              value[i][j][k] = valueView(i,j,k);
1791            }
1792          }
1793        }
1794      }
1795      if (rank==4) {
1796        for (int i=0; i < shape[0]; i++) {
1797          for (int j=0; j < shape[1]; j++) {
1798            for (int k=0; k < shape[2]; k++) {
1799              for (int l=0; l < shape[3]; l++) {
1800                value[i][j][k][l] = valueView(i,j,k,l);
1801              }
1802            }
1803          }
1804        }
1805      }
1806    
1807    }
1808    
1809    void
1810    Data::archiveData(const std::string fileName)
1811    {
1812      cout << "Archiving Data object to: " << fileName << endl;
1813    
1814      //
1815      // Determine type of this Data object
1816      int dataType = -1;
1817    
1818      if (isEmpty()) {
1819        dataType = 0;
1820        cout << "\tdataType: DataEmpty" << endl;
1821      }
1822      if (isConstant()) {
1823        dataType = 1;
1824        cout << "\tdataType: DataConstant" << endl;
1825      }
1826      if (isTagged()) {
1827        dataType = 2;
1828        cout << "\tdataType: DataTagged" << endl;
1829      }
1830      if (isExpanded()) {
1831        dataType = 3;
1832        cout << "\tdataType: DataExpanded" << endl;
1833      }
1834    
1835      if (dataType == -1) {
1836        throw DataException("archiveData Error: undefined dataType");
1837      }
1838    
1839      //
1840      // Collect data items common to all Data types
1841      int noSamples = getNumSamples();
1842      int noDPPSample = getNumDataPointsPerSample();
1843      int functionSpaceType = getFunctionSpace().getTypeCode();
1844      int dataPointRank = getDataPointRank();
1845      int dataPointSize = getDataPointSize();
1846      int dataLength = getLength();
1847      DataArrayView::ShapeType dataPointShape = getDataPointShape();
1848      int referenceNumbers[noSamples];
1849      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1850        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1851      }
1852      int tagNumbers[noSamples];
1853      if (isTagged()) {
1854        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1855          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1856        }
1857      }
1858    
1859      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1860      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1861      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1862    
1863      //
1864      // Flatten Shape to an array of integers suitable for writing to file
1865      int flatShape[4] = {0,0,0,0};
1866      cout << "\tshape: < ";
1867      for (int dim=0; dim<dataPointRank; dim++) {
1868        flatShape[dim] = dataPointShape[dim];
1869        cout << dataPointShape[dim] << " ";
1870      }
1871      cout << ">" << endl;
1872    
1873      //
1874      // Open archive file
1875      ofstream archiveFile;
1876      archiveFile.open(fileName.data(), ios::out);
1877    
1878      if (!archiveFile.good()) {
1879        throw DataException("archiveData Error: problem opening archive file");
1880      }
1881    
1882      //
1883      // Write common data items to archive file
1884      archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1885      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1886      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1887      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1888      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1889      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1890      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1891      for (int dim = 0; dim < 4; dim++) {
1892        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1893      }
1894      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1895        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1896      }
1897      if (isTagged()) {
1898        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1899          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1900        }
1901      }
1902    
1903      if (!archiveFile.good()) {
1904        throw DataException("archiveData Error: problem writing to archive file");
1905      }
1906    
1907      //
1908      // Archive underlying data values for each Data type
1909      int noValues;
1910      switch (dataType) {
1911        case 0:
1912          // DataEmpty
1913          noValues = 0;
1914          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1915          cout << "\tnoValues: " << noValues << endl;
1916          break;
1917        case 1:
1918          // DataConstant
1919          noValues = m_data->getLength();
1920          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1921          cout << "\tnoValues: " << noValues << endl;
1922          if (m_data->archiveData(archiveFile,noValues)) {
1923            throw DataException("archiveData Error: problem writing data to archive file");
1924          }
1925          break;
1926        case 2:
1927          // DataTagged
1928          noValues = m_data->getLength();
1929          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1930          cout << "\tnoValues: " << noValues << endl;
1931          if (m_data->archiveData(archiveFile,noValues)) {
1932            throw DataException("archiveData Error: problem writing data to archive file");
1933          }
1934          break;
1935        case 3:
1936          // DataExpanded
1937          noValues = m_data->getLength();
1938          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1939          cout << "\tnoValues: " << noValues << endl;
1940          if (m_data->archiveData(archiveFile,noValues)) {
1941            throw DataException("archiveData Error: problem writing data to archive file");
1942          }
1943          break;
1944      }
1945    
1946      if (!archiveFile.good()) {
1947        throw DataException("archiveData Error: problem writing data to archive file");
1948      }
1949    
1950      //
1951      // Close archive file
1952      archiveFile.close();
1953    
1954      if (!archiveFile.good()) {
1955        throw DataException("archiveData Error: problem closing archive file");
1956      }
1957    
1958    }
1959    
1960    void
1961    Data::extractData(const std::string fileName,
1962                      const FunctionSpace& fspace)
1963    {
1964      //
1965      // Can only extract Data to an object which is initially DataEmpty
1966      if (!isEmpty()) {
1967        throw DataException("extractData Error: can only extract to DataEmpty object");
1968      }
1969    
1970      cout << "Extracting Data object from: " << fileName << endl;
1971    
1972      int dataType;
1973      int noSamples;
1974      int noDPPSample;
1975      int functionSpaceType;
1976      int dataPointRank;
1977      int dataPointSize;
1978      int dataLength;
1979      DataArrayView::ShapeType dataPointShape;
1980      int flatShape[4];
1981    
1982      //
1983      // Open the archive file
1984      ifstream archiveFile;
1985      archiveFile.open(fileName.data(), ios::in);
1986    
1987      if (!archiveFile.good()) {
1988        throw DataException("extractData Error: problem opening archive file");
1989      }
1990    
1991      //
1992      // Read common data items from archive file
1993      archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1994      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1995      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1996      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1997      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1998      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1999      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2000      for (int dim = 0; dim < 4; dim++) {
2001        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2002        if (flatShape[dim]>0) {
2003          dataPointShape.push_back(flatShape[dim]);
2004        }
2005      }
2006      int referenceNumbers[noSamples];
2007      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2008        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2009      }
2010      int tagNumbers[noSamples];
2011      if (dataType==2) {
2012        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2013          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2014        }
2015      }
2016    
2017      if (!archiveFile.good()) {
2018        throw DataException("extractData Error: problem reading from archive file");
2019      }
2020    
2021      //
2022      // Verify the values just read from the archive file
2023      switch (dataType) {
2024        case 0:
2025          cout << "\tdataType: DataEmpty" << endl;
2026          break;
2027        case 1:
2028          cout << "\tdataType: DataConstant" << endl;
2029          break;
2030        case 2:
2031          cout << "\tdataType: DataTagged" << endl;
2032          break;
2033        case 3:
2034          cout << "\tdataType: DataExpanded" << endl;
2035          break;
2036        default:
2037          throw DataException("extractData Error: undefined dataType read from archive file");
2038          break;
2039      }
2040    
2041      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2042      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2043      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2044      cout << "\tshape: < ";
2045      for (int dim = 0; dim < dataPointRank; dim++) {
2046        cout << dataPointShape[dim] << " ";
2047      }
2048      cout << ">" << endl;
2049    
2050      //
2051      // Verify that supplied FunctionSpace object is compatible with this Data object.
2052      if ( (fspace.getTypeCode()!=functionSpaceType) ||
2053           (fspace.getNumSamples()!=noSamples) ||
2054           (fspace.getNumDPPSample()!=noDPPSample)
2055         ) {
2056        throw DataException("extractData Error: incompatible FunctionSpace");
2057      }
2058      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2059        if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2060          throw DataException("extractData Error: incompatible FunctionSpace");
2061        }
2062      }
2063      if (dataType==2) {
2064        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2065          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2066            throw DataException("extractData Error: incompatible FunctionSpace");
2067          }
2068        }
2069      }
2070    
2071      //
2072      // Construct a DataVector to hold underlying data values
2073      DataVector dataVec(dataLength);
2074    
2075      //
2076      // Load this DataVector with the appropriate values
2077      int noValues;
2078      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2079      cout << "\tnoValues: " << noValues << endl;
2080      switch (dataType) {
2081        case 0:
2082          // DataEmpty
2083          if (noValues != 0) {
2084            throw DataException("extractData Error: problem reading data from archive file");
2085          }
2086          break;
2087        case 1:
2088          // DataConstant
2089          if (dataVec.extractData(archiveFile,noValues)) {
2090            throw DataException("extractData Error: problem reading data from archive file");
2091          }
2092          break;
2093        case 2:
2094          // DataTagged
2095          if (dataVec.extractData(archiveFile,noValues)) {
2096            throw DataException("extractData Error: problem reading data from archive file");
2097          }
2098          break;
2099        case 3:
2100          // DataExpanded
2101          if (dataVec.extractData(archiveFile,noValues)) {
2102            throw DataException("extractData Error: problem reading data from archive file");
2103          }
2104          break;
2105      }
2106    
2107      if (!archiveFile.good()) {
2108        throw DataException("extractData Error: problem reading from archive file");
2109      }
2110    
2111      //
2112      // Close archive file
2113      archiveFile.close();
2114    
2115      if (!archiveFile.good()) {
2116        throw DataException("extractData Error: problem closing archive file");
2117      }
2118    
2119      //
2120      // Construct an appropriate Data object
2121      DataAbstract* tempData;
2122      switch (dataType) {
2123        case 0:
2124          // DataEmpty
2125          tempData=new DataEmpty();
2126          break;
2127        case 1:
2128          // DataConstant
2129          tempData=new DataConstant(fspace,dataPointShape,dataVec);
2130          break;
2131        case 2:
2132          // DataTagged
2133          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2134          break;
2135        case 3:
2136          // DataExpanded
2137          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2138          break;
2139      }
2140      shared_ptr<DataAbstract> temp_data(tempData);
2141      m_data=temp_data;
2142    }
2143    
2144  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
2145  {  {

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

  ViewVC Help
Powered by ViewVC 1.1.26