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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26