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

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

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

trunk/esys2/escript/src/Data/Data.cpp revision 117 by jgs, Fri Apr 1 05:48:57 2005 UTC trunk/escript/src/Data/Data.cpp revision 286 by gross, Thu Dec 1 06:40: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  boost::python::numeric::array
554  Data::convertToNumArray()  Data::convertToNumArray()
555  {  {
556    //    //
557    // determine the current number of data points    // determine the total number of data points
558    int numSamples = getNumSamples();    int numSamples = getNumSamples();
559    int numDataPointsPerSample = getNumDataPointsPerSample();    int numDataPointsPerSample = getNumDataPointsPerSample();
560    int numDataPoints = numSamples * numDataPointsPerSample;    int numDataPoints = numSamples * numDataPointsPerSample;
# Line 499  Data::convertToNumArray() Line 606  Data::convertToNumArray()
606    int dataPoint = 0;    int dataPoint = 0;
607    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
608      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
   
609        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
   
610        if (dataPointRank==0) {        if (dataPointRank==0) {
611          numArray[dataPoint]=dataPointView();          numArray[dataPoint]=dataPointView();
612        }        }
   
613        if (dataPointRank==1) {        if (dataPointRank==1) {
614          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
615            numArray[dataPoint][i]=dataPointView(i);            numArray[dataPoint][i]=dataPointView(i);
616          }          }
617        }        }
   
618        if (dataPointRank==2) {        if (dataPointRank==2) {
619          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
620            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
# Line 519  Data::convertToNumArray() Line 622  Data::convertToNumArray()
622            }            }
623          }          }
624        }        }
   
625        if (dataPointRank==3) {        if (dataPointRank==3) {
626          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
627            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
# Line 529  Data::convertToNumArray() Line 631  Data::convertToNumArray()
631            }            }
632          }          }
633        }        }
   
634        if (dataPointRank==4) {        if (dataPointRank==4) {
635          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
636            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
# Line 541  Data::convertToNumArray() Line 642  Data::convertToNumArray()
642            }            }
643          }          }
644        }        }
   
645        dataPoint++;        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    
# Line 559  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 622  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::log() const  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::log10() 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::log() 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  Data
1048  Data::sign() const  Data::sign() const
1049  {  {
1050    #if defined DOPROF
1051      profData->unary++;
1052    #endif
1053    return escript::unaryOp(*this,escript::fsign);    return escript::unaryOp(*this,escript::fsign);
1054  }  }
1055    
1056  Data  Data
1057  Data::abs() const  Data::abs() const
1058  {  {
1059    #if defined DOPROF
1060      profData->unary++;
1061    #endif
1062    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
1063  }  }
1064    
1065  Data  Data
1066  Data::neg() const  Data::neg() const
1067  {  {
1068    #if defined DOPROF
1069      profData->unary++;
1070    #endif
1071    return escript::unaryOp(*this,negate<double>());    return escript::unaryOp(*this,negate<double>());
1072  }  }
1073    
1074  Data  Data
1075  Data::pos() const  Data::pos() const
1076  {  {
1077    return (*this);  #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  Data
1087  Data::exp() const  Data::exp() const
1088  {  {
1089    #if defined DOPROF
1090      profData->unary++;
1091    #endif
1092    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1093  }  }
1094    
1095  Data  Data
1096  Data::sqrt() const  Data::sqrt() const
1097  {  {
1098    #if defined DOPROF
1099      profData->unary++;
1100    #endif
1101    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    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  double
1117  Data::Linf() const  Data::Linf() const
1118  {  {
1119    #if defined DOPROF
1120      profData->reduction1++;
1121    #endif
1122    //    //
1123    // set the initial absolute minimum value to max double    // set the initial absolute minimum value to max double
1124    return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));    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>::max()*-1));    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    #if defined DOPROF
1156      profData->reduction2++;
1157    #endif
1158    //    //
1159    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1160    return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    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    #if defined DOPROF
1168      profData->reduction2++;
1169    #endif
1170    //    //
1171    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1172    return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1173  }    return dp_algorithm(fmin_func,numeric_limits<double>::max());
   
 Data  
 Data::length() const  
 {  
   return dp_algorithm(DataAlgorithmAdapter<Length>(0));  
1174  }  }
1175    
1176  Data  Data
1177  Data::trace() const  Data::trace() const
1178  {  {
1179    return dp_algorithm(DataAlgorithmAdapter<Trace>(0));  #if defined DOPROF
1180      profData->reduction2++;
1181    #endif
1182      Trace trace_func;
1183      return dp_algorithm(trace_func,0);
1184  }  }
1185    
1186  Data  Data
1187  Data::transpose(int axis) const  Data::transpose(int axis) const
1188  {  {
1189    #if defined DOPROF
1190      profData->reduction2++;
1191    #endif
1192    // not implemented    // not implemented
1193    throw DataException("Error - Data::transpose not implemented yet.");    throw DataException("Error - Data::transpose not implemented yet.");
1194    return Data();    return Data();
1195  }  }
1196    
1197    const boost::python::tuple
1198    Data::mindp() const
1199    {
1200      // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1201      // abort (for unknown reasons) if there are openmp directives with it in the
1202      // surrounding function
1203    
1204      int SampleNo;
1205      int DataPointNo;
1206    
1207      calc_mindp(SampleNo,DataPointNo);
1208    
1209      return make_tuple(SampleNo,DataPointNo);
1210    }
1211    
1212    void
1213    Data::calc_mindp(int& SampleNo,
1214                     int& DataPointNo) const
1215    {
1216      int i,j;
1217      int lowi=0,lowj=0;
1218      double min=numeric_limits<double>::max();
1219    
1220      Data temp=minval();
1221    
1222      int numSamples=temp.getNumSamples();
1223      int numDPPSample=temp.getNumDataPointsPerSample();
1224    
1225      double next,local_min;
1226      int local_lowi,local_lowj;
1227    
1228      #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1229      {
1230        local_min=min;
1231        #pragma omp for private(i,j) schedule(static)
1232        for (i=0; i<numSamples; i++) {
1233          for (j=0; j<numDPPSample; j++) {
1234            next=temp.getDataPoint(i,j)();
1235            if (next<local_min) {
1236              local_min=next;
1237              local_lowi=i;
1238              local_lowj=j;
1239            }
1240          }
1241        }
1242        #pragma omp critical
1243        if (local_min<min) {
1244          min=local_min;
1245          lowi=local_lowi;
1246          lowj=local_lowj;
1247        }
1248      }
1249    
1250      SampleNo = lowi;
1251      DataPointNo = lowj;
1252    }
1253    
1254  void  void
1255  Data::saveDX(std::string fileName) const  Data::saveDX(std::string fileName) const
1256  {  {
1257    getDomain().saveDX(fileName,*this);    boost::python::dict args;
1258      args["data"]=boost::python::object(this);
1259      getDomain().saveDX(fileName,args);
1260    return;    return;
1261  }  }
1262    
1263  void  void
1264  Data::saveVTK(std::string fileName) const  Data::saveVTK(std::string fileName) const
1265  {  {
1266    getDomain().saveVTK(fileName,*this);    boost::python::dict args;
1267      args["data"]=boost::python::object(this);
1268      getDomain().saveVTK(fileName,args);
1269    return;    return;
1270  }  }
1271    
1272  Data&  Data&
1273  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1274  {  {
1275    #if defined DOPROF
1276      profData->binary++;
1277    #endif
1278    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1279    return (*this);    return (*this);
1280  }  }
# Line 777  Data::operator+=(const Data& right) Line 1282  Data::operator+=(const Data& right)
1282  Data&  Data&
1283  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& 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 784  Data::operator+=(const boost::python::ob Line 1292  Data::operator+=(const boost::python::ob
1292  Data&  Data&
1293  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1294  {  {
1295    #if defined DOPROF
1296      profData->binary++;
1297    #endif
1298    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1299    return (*this);    return (*this);
1300  }  }
# Line 791  Data::operator-=(const Data& right) Line 1302  Data::operator-=(const Data& right)
1302  Data&  Data&
1303  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& 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 798  Data::operator-=(const boost::python::ob Line 1312  Data::operator-=(const boost::python::ob
1312  Data&  Data&
1313  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1314  {  {
1315    #if defined DOPROF
1316      profData->binary++;
1317    #endif
1318    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1319    return (*this);    return (*this);
1320  }  }
# Line 805  Data::operator*=(const Data& right) Line 1322  Data::operator*=(const Data& right)
1322  Data&  Data&
1323  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& 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 812  Data::operator*=(const boost::python::ob Line 1332  Data::operator*=(const boost::python::ob
1332  Data&  Data&
1333  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1334  {  {
1335    #if defined DOPROF
1336      profData->binary++;
1337    #endif
1338    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1339    return (*this);    return (*this);
1340  }  }
# Line 819  Data::operator/=(const Data& right) Line 1342  Data::operator/=(const Data& right)
1342  Data&  Data&
1343  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& 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 826  Data::operator/=(const boost::python::ob Line 1352  Data::operator/=(const boost::python::ob
1352  Data  Data
1353  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1354  {  {
1355    #if defined DOPROF
1356      profData->binary++;
1357    #endif
1358    Data result;    Data result;
1359    result.copy(*this);    result.copy(*this);
1360    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 835  Data::powO(const boost::python::object& Line 1364  Data::powO(const boost::python::object&
1364  Data  Data
1365  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1366  {  {
1367    #if defined DOPROF
1368      profData->binary++;
1369    #endif
1370    Data result;    Data result;
1371    result.copy(*this);    result.copy(*this);
1372    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 842  Data::powD(const Data& right) const Line 1374  Data::powD(const Data& right) const
1374  }  }
1375    
1376  //  //
1377  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1378  Data  Data
1379  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1380  {  {
# Line 855  escript::operator+(const Data& left, con Line 1387  escript::operator+(const Data& left, con
1387  }  }
1388    
1389  //  //
1390  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1391  Data  Data
1392  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1393  {  {
# Line 868  escript::operator-(const Data& left, con Line 1400  escript::operator-(const Data& left, con
1400  }  }
1401    
1402  //  //
1403  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1404  Data  Data
1405  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
1406  {  {
# Line 881  escript::operator*(const Data& left, con Line 1413  escript::operator*(const Data& left, con
1413  }  }
1414    
1415  //  //
1416  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1417  Data  Data
1418  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
1419  {  {
# Line 894  escript::operator/(const Data& left, con Line 1426  escript::operator/(const Data& left, con
1426  }  }
1427    
1428  //  //
1429  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1430  Data  Data
1431  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1432  {  {
# Line 910  escript::operator+(const Data& left, con Line 1442  escript::operator+(const Data& left, con
1442  }  }
1443    
1444  //  //
1445  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1446  Data  Data
1447  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1448  {  {
# Line 926  escript::operator-(const Data& left, con Line 1458  escript::operator-(const Data& left, con
1458  }  }
1459    
1460  //  //
1461  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1462  Data  Data
1463  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1464  {  {
# Line 942  escript::operator*(const Data& left, con Line 1474  escript::operator*(const Data& left, con
1474  }  }
1475    
1476  //  //
1477  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1478  Data  Data
1479  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1480  {  {
# Line 958  escript::operator/(const Data& left, con Line 1490  escript::operator/(const Data& left, con
1490  }  }
1491    
1492  //  //
1493  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1494  Data  Data
1495  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1496  {  {
# Line 971  escript::operator+(const boost::python:: Line 1503  escript::operator+(const boost::python::
1503  }  }
1504    
1505  //  //
1506  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1507  Data  Data
1508  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1509  {  {
# Line 984  escript::operator-(const boost::python:: Line 1516  escript::operator-(const boost::python::
1516  }  }
1517    
1518  //  //
1519  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1520  Data  Data
1521  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1522  {  {
# Line 997  escript::operator*(const boost::python:: Line 1529  escript::operator*(const boost::python::
1529  }  }
1530    
1531  //  //
1532  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1533  Data  Data
1534  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1535  {  {
# Line 1010  escript::operator/(const boost::python:: Line 1542  escript::operator/(const boost::python::
1542  }  }
1543    
1544  //  //
 // NOTE: It is essential to specify the namepsace this operator belongs to  
1545  //bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1546  //{  //{
1547  //  /*  //  /*
# Line 1072  Data::getItem(const boost::python::objec Line 1603  Data::getItem(const boost::python::objec
1603  Data  Data
1604  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1605  {  {
1606    #if defined DOPROF
1607      profData->slicing++;
1608    #endif
1609    return Data(*this,region);    return Data(*this,region);
1610  }  }
1611    
# Line 1088  Data::setItemD(const boost::python::obje Line 1622  Data::setItemD(const boost::python::obje
1622                 const Data& value)                 const Data& value)
1623  {  {
1624    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1625    
1626    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1627    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1628      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
# Line 1103  void Line 1638  void
1638  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1639                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
1640  {  {
1641    #if defined DOPROF
1642      profData->slicing++;
1643    #endif
1644    Data tempValue(value);    Data tempValue(value);
1645    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1646    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1159  Data::setTaggedValue(int tagKey, Line 1697  Data::setTaggedValue(int tagKey,
1697  }  }
1698    
1699  void  void
1700    Data::setTaggedValueFromCPP(int tagKey,
1701                                const DataArrayView& value)
1702    {
1703      //
1704      // Ensure underlying data object is of type DataTagged
1705      tag();
1706    
1707      if (!isTagged()) {
1708        throw DataException("Error - DataTagged conversion failed!!");
1709      }
1710                                                                                                                  
1711      //
1712      // Call DataAbstract::setTaggedValue
1713      m_data->setTaggedValue(tagKey,value);
1714    }
1715    
1716    int
1717    Data::getTagNumber(int dpno)
1718    {
1719      return m_data->getTagNumber(dpno);
1720    }
1721    
1722    void
1723  Data::setRefValue(int ref,  Data::setRefValue(int ref,
1724                    const boost::python::numeric::array& value)                    const boost::python::numeric::array& value)
1725  {  {
# Line 1197  Data::getRefValue(int ref, Line 1758  Data::getRefValue(int ref,
1758    DataArrayView valueView = valueDataArray.getView();    DataArrayView valueView = valueDataArray.getView();
1759    
1760    if (rank==0) {    if (rank==0) {
1761      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");        boost::python::numeric::array temp_numArray(valueView());
1762          value = temp_numArray;
1763    }    }
1764    if (rank==1) {    if (rank==1) {
1765      for (int i=0; i < shape[0]; i++) {      for (int i=0; i < shape[0]; i++) {
# Line 1205  Data::getRefValue(int ref, Line 1767  Data::getRefValue(int ref,
1767      }      }
1768    }    }
1769    if (rank==2) {    if (rank==2) {
1770      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
1771          for (int j=0; j < shape[1]; j++) {
1772            value[i][j] = valueView(i,j);
1773          }
1774        }
1775    }    }
1776    if (rank==3) {    if (rank==3) {
1777      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
1778          for (int j=0; j < shape[1]; j++) {
1779            for (int k=0; k < shape[2]; k++) {
1780              value[i][j][k] = valueView(i,j,k);
1781            }
1782          }
1783        }
1784    }    }
1785    if (rank==4) {    if (rank==4) {
1786      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
1787          for (int j=0; j < shape[1]; j++) {
1788            for (int k=0; k < shape[2]; k++) {
1789              for (int l=0; l < shape[3]; l++) {
1790                value[i][j][k][l] = valueView(i,j,k,l);
1791              }
1792            }
1793          }
1794        }
1795    }    }
1796    
1797  }  }
1798    
 /*  
 Note: this version removed for now. Not needed, and breaks escript.cpp  
1799  void  void
1800  Data::setTaggedValue(int tagKey,  Data::archiveData(const std::string fileName)
                      const DataArrayView& value)  
1801  {  {
1802      cout << "Archiving Data object to: " << fileName << endl;
1803    
1804    //    //
1805    // Ensure underlying data object is of type DataTagged    // Determine type of this Data object
1806    tag();    int dataType = -1;
1807    
1808    if (!isTagged()) {    if (isEmpty()) {
1809      throw DataException("Error - DataTagged conversion failed!!");      dataType = 0;
1810        cout << "\tdataType: DataEmpty" << endl;
1811    }    }
1812                                                                                                                    if (isConstant()) {
1813        dataType = 1;
1814        cout << "\tdataType: DataConstant" << endl;
1815      }
1816      if (isTagged()) {
1817        dataType = 2;
1818        cout << "\tdataType: DataTagged" << endl;
1819      }
1820      if (isExpanded()) {
1821        dataType = 3;
1822        cout << "\tdataType: DataExpanded" << endl;
1823      }
1824    
1825      if (dataType == -1) {
1826        throw DataException("archiveData Error: undefined dataType");
1827      }
1828    
1829    //    //
1830    // Call DataAbstract::setTaggedValue    // Collect data items common to all Data types
1831    m_data->setTaggedValue(tagKey,value);    int noSamples = getNumSamples();
1832      int noDPPSample = getNumDataPointsPerSample();
1833      int functionSpaceType = getFunctionSpace().getTypeCode();
1834      int dataPointRank = getDataPointRank();
1835      int dataPointSize = getDataPointSize();
1836      int dataLength = getLength();
1837      DataArrayView::ShapeType dataPointShape = getDataPointShape();
1838      int referenceNumbers[noSamples];
1839      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1840        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1841      }
1842      int tagNumbers[noSamples];
1843      if (isTagged()) {
1844        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1845          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1846        }
1847      }
1848    
1849      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1850      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1851      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1852    
1853      //
1854      // Flatten Shape to an array of integers suitable for writing to file
1855      int flatShape[4] = {0,0,0,0};
1856      cout << "\tshape: < ";
1857      for (int dim=0; dim<dataPointRank; dim++) {
1858        flatShape[dim] = dataPointShape[dim];
1859        cout << dataPointShape[dim] << " ";
1860      }
1861      cout << ">" << endl;
1862    
1863      //
1864      // Open archive file
1865      ofstream archiveFile;
1866      archiveFile.open(fileName.data(), ios::out);
1867    
1868      if (!archiveFile.good()) {
1869        throw DataException("archiveData Error: problem opening archive file");
1870      }
1871    
1872      //
1873      // Write common data items to archive file
1874      archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1875      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1876      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1877      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1878      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1879      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1880      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1881      for (int dim = 0; dim < 4; dim++) {
1882        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1883      }
1884      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1885        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1886      }
1887      if (isTagged()) {
1888        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1889          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1890        }
1891      }
1892    
1893      if (!archiveFile.good()) {
1894        throw DataException("archiveData Error: problem writing to archive file");
1895      }
1896    
1897      //
1898      // Archive underlying data values for each Data type
1899      int noValues;
1900      switch (dataType) {
1901        case 0:
1902          // DataEmpty
1903          noValues = 0;
1904          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1905          cout << "\tnoValues: " << noValues << endl;
1906          break;
1907        case 1:
1908          // DataConstant
1909          noValues = m_data->getLength();
1910          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1911          cout << "\tnoValues: " << noValues << endl;
1912          if (m_data->archiveData(archiveFile,noValues)) {
1913            throw DataException("archiveData Error: problem writing data to archive file");
1914          }
1915          break;
1916        case 2:
1917          // DataTagged
1918          noValues = m_data->getLength();
1919          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1920          cout << "\tnoValues: " << noValues << endl;
1921          if (m_data->archiveData(archiveFile,noValues)) {
1922            throw DataException("archiveData Error: problem writing data to archive file");
1923          }
1924          break;
1925        case 3:
1926          // DataExpanded
1927          noValues = m_data->getLength();
1928          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1929          cout << "\tnoValues: " << noValues << endl;
1930          if (m_data->archiveData(archiveFile,noValues)) {
1931            throw DataException("archiveData Error: problem writing data to archive file");
1932          }
1933          break;
1934      }
1935    
1936      if (!archiveFile.good()) {
1937        throw DataException("archiveData Error: problem writing data to archive file");
1938      }
1939    
1940      //
1941      // Close archive file
1942      archiveFile.close();
1943    
1944      if (!archiveFile.good()) {
1945        throw DataException("archiveData Error: problem closing archive file");
1946      }
1947    
1948    }
1949    
1950    void
1951    Data::extractData(const std::string fileName,
1952                      const FunctionSpace& fspace)
1953    {
1954      //
1955      // Can only extract Data to an object which is initially DataEmpty
1956      if (!isEmpty()) {
1957        throw DataException("extractData Error: can only extract to DataEmpty object");
1958      }
1959    
1960      cout << "Extracting Data object from: " << fileName << endl;
1961    
1962      int dataType;
1963      int noSamples;
1964      int noDPPSample;
1965      int functionSpaceType;
1966      int dataPointRank;
1967      int dataPointSize;
1968      int dataLength;
1969      DataArrayView::ShapeType dataPointShape;
1970      int flatShape[4];
1971    
1972      //
1973      // Open the archive file
1974      ifstream archiveFile;
1975      archiveFile.open(fileName.data(), ios::in);
1976    
1977      if (!archiveFile.good()) {
1978        throw DataException("extractData Error: problem opening archive file");
1979      }
1980    
1981      //
1982      // Read common data items from archive file
1983      archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1984      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1985      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1986      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1987      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1988      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1989      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1990      for (int dim = 0; dim < 4; dim++) {
1991        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1992        if (flatShape[dim]>0) {
1993          dataPointShape.push_back(flatShape[dim]);
1994        }
1995      }
1996      int referenceNumbers[noSamples];
1997      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1998        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1999      }
2000      int tagNumbers[noSamples];
2001      if (dataType==2) {
2002        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2003          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2004        }
2005      }
2006    
2007      if (!archiveFile.good()) {
2008        throw DataException("extractData Error: problem reading from archive file");
2009      }
2010    
2011      //
2012      // Verify the values just read from the archive file
2013      switch (dataType) {
2014        case 0:
2015          cout << "\tdataType: DataEmpty" << endl;
2016          break;
2017        case 1:
2018          cout << "\tdataType: DataConstant" << endl;
2019          break;
2020        case 2:
2021          cout << "\tdataType: DataTagged" << endl;
2022          break;
2023        case 3:
2024          cout << "\tdataType: DataExpanded" << endl;
2025          break;
2026        default:
2027          throw DataException("extractData Error: undefined dataType read from archive file");
2028          break;
2029      }
2030    
2031      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2032      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2033      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2034      cout << "\tshape: < ";
2035      for (int dim = 0; dim < dataPointRank; dim++) {
2036        cout << dataPointShape[dim] << " ";
2037      }
2038      cout << ">" << endl;
2039    
2040      //
2041      // Verify that supplied FunctionSpace object is compatible with this Data object.
2042      if ( (fspace.getTypeCode()!=functionSpaceType) ||
2043           (fspace.getNumSamples()!=noSamples) ||
2044           (fspace.getNumDPPSample()!=noDPPSample)
2045         ) {
2046        throw DataException("extractData Error: incompatible FunctionSpace");
2047      }
2048      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2049        if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2050          throw DataException("extractData Error: incompatible FunctionSpace");
2051        }
2052      }
2053      if (dataType==2) {
2054        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2055          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2056            throw DataException("extractData Error: incompatible FunctionSpace");
2057          }
2058        }
2059      }
2060    
2061      //
2062      // Construct a DataVector to hold underlying data values
2063      DataVector dataVec(dataLength);
2064    
2065      //
2066      // Load this DataVector with the appropriate values
2067      int noValues;
2068      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2069      cout << "\tnoValues: " << noValues << endl;
2070      switch (dataType) {
2071        case 0:
2072          // DataEmpty
2073          if (noValues != 0) {
2074            throw DataException("extractData Error: problem reading data from archive file");
2075          }
2076          break;
2077        case 1:
2078          // DataConstant
2079          if (dataVec.extractData(archiveFile,noValues)) {
2080            throw DataException("extractData Error: problem reading data from archive file");
2081          }
2082          break;
2083        case 2:
2084          // DataTagged
2085          if (dataVec.extractData(archiveFile,noValues)) {
2086            throw DataException("extractData Error: problem reading data from archive file");
2087          }
2088          break;
2089        case 3:
2090          // DataExpanded
2091          if (dataVec.extractData(archiveFile,noValues)) {
2092            throw DataException("extractData Error: problem reading data from archive file");
2093          }
2094          break;
2095      }
2096    
2097      if (!archiveFile.good()) {
2098        throw DataException("extractData Error: problem reading from archive file");
2099      }
2100    
2101      //
2102      // Close archive file
2103      archiveFile.close();
2104    
2105      if (!archiveFile.good()) {
2106        throw DataException("extractData Error: problem closing archive file");
2107      }
2108    
2109      //
2110      // Construct an appropriate Data object
2111      DataAbstract* tempData;
2112      switch (dataType) {
2113        case 0:
2114          // DataEmpty
2115          tempData=new DataEmpty();
2116          break;
2117        case 1:
2118          // DataConstant
2119          tempData=new DataConstant(fspace,dataPointShape,dataVec);
2120          break;
2121        case 2:
2122          // DataTagged
2123          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2124          break;
2125        case 3:
2126          // DataExpanded
2127          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2128          break;
2129      }
2130      shared_ptr<DataAbstract> temp_data(tempData);
2131      m_data=temp_data;
2132  }  }
 */  
2133    
2134  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
2135  {  {

Legend:
Removed from v.117  
changed lines
  Added in v.286

  ViewVC Help
Powered by ViewVC 1.1.26