/[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 149 by jgs, Thu Sep 1 03:31:39 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    //
49    // global table of profiling data for all Data objects
50    DataProf dataProfTable;
51    
52  Data::Data()  Data::Data()
53  {  {
54    //    //
# Line 52  Data::Data() Line 56  Data::Data()
56    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
57    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
58    m_data=temp_data;    m_data=temp_data;
59      // create entry in global profiling table for this object
60      profData = dataProfTable.newData();
61  }  }
62    
63  Data::Data(double value,  Data::Data(double value,
# Line 65  Data::Data(double value, Line 71  Data::Data(double value,
71    }    }
72    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
73    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
74      // create entry in global profiling table for this object
75      profData = dataProfTable.newData();
76  }  }
77    
78  Data::Data(double value,  Data::Data(double value,
# Line 75  Data::Data(double value, Line 83  Data::Data(double value,
83    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
84    pair<int,int> dataShape=what.getDataShape();    pair<int,int> dataShape=what.getDataShape();
85    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
86      // create entry in global profiling table for this object
87      profData = dataProfTable.newData();
88  }  }
89    
90  Data::Data(const Data& inData)  Data::Data(const Data& inData)
91  {  {
92    m_data=inData.m_data;    m_data=inData.m_data;
93      // create entry in global profiling table for this object
94      profData = dataProfTable.newData();
95  }  }
96    
97  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 90  Data::Data(const Data& inData, Line 102  Data::Data(const Data& inData,
102    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
103    shared_ptr<DataAbstract> temp_data(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
104    m_data=temp_data;    m_data=temp_data;
105      // create entry in global profiling table for this object
106      profData = dataProfTable.newData();
107  }  }
108    
109  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 99  Data::Data(const Data& inData, Line 113  Data::Data(const Data& inData,
113      m_data=inData.m_data;      m_data=inData.m_data;
114    } else {    } else {
115      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
116      // Note for Lutz, Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
117      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
118      // be able to create an instance of AbstractDomain but that was done      // be able to create an instance of AbstractDomain but that was done
119      // as a boost python work around which may no longer be required.      // as a boost:python work around which may no longer be required.
120      const AbstractDomain& inDataDomain=inData.getDomain();      const AbstractDomain& inDataDomain=inData.getDomain();
121      if  (inDataDomain==functionspace.getDomain()) {      if  (inDataDomain==functionspace.getDomain()) {
122        inDataDomain.interpolateOnDomain(tmp,inData);        inDataDomain.interpolateOnDomain(tmp,inData);
# Line 111  Data::Data(const Data& inData, Line 125  Data::Data(const Data& inData,
125      }      }
126      m_data=tmp.m_data;      m_data=tmp.m_data;
127    }    }
128      // create entry in global profiling table for this object
129      profData = dataProfTable.newData();
130  }  }
131    
132  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(const DataTagged::TagListType& tagKeys,
# Line 125  Data::Data(const DataTagged::TagListType Line 141  Data::Data(const DataTagged::TagListType
141    if (expanded) {    if (expanded) {
142      expand();      expand();
143    }    }
144      // create entry in global profiling table for this object
145      profData = dataProfTable.newData();
146  }  }
147    
148  Data::Data(const numeric::array& value,  Data::Data(const numeric::array& value,
# Line 132  Data::Data(const numeric::array& value, Line 150  Data::Data(const numeric::array& value,
150             bool expanded)             bool expanded)
151  {  {
152    initialise(value,what,expanded);    initialise(value,what,expanded);
153      // create entry in global profiling table for this object
154      profData = dataProfTable.newData();
155  }  }
156    
157  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
# Line 139  Data::Data(const DataArrayView& value, Line 159  Data::Data(const DataArrayView& value,
159             bool expanded)             bool expanded)
160  {  {
161    initialise(value,what,expanded);    initialise(value,what,expanded);
162      // create entry in global profiling table for this object
163      profData = dataProfTable.newData();
164  }  }
165    
166  Data::Data(const object& value,  Data::Data(const object& value,
# Line 147  Data::Data(const object& value, Line 169  Data::Data(const object& value,
169  {  {
170    numeric::array asNumArray(value);    numeric::array asNumArray(value);
171    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
172      // create entry in global profiling table for this object
173      profData = dataProfTable.newData();
174  }  }
175    
176  Data::Data(const object& value,  Data::Data(const object& value,
# Line 167  Data::Data(const object& value, Line 191  Data::Data(const object& value,
191      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
192      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
193    }    }
194      // create entry in global profiling table for this object
195      profData = dataProfTable.newData();
196  }  }
197    
198  escriptDataC  escriptDataC
# Line 185  Data::getDataC() const Line 211  Data::getDataC() const
211    return temp;    return temp;
212  }  }
213    
214  tuple  const boost::python::tuple
215  Data::getShapeTuple() const  Data::getShapeTuple() const
216  {  {
217    const DataArrayView::ShapeType& shape=getDataPointShape();    const DataArrayView::ShapeType& shape=getDataPointShape();
# Line 353  Data::reshapeDataPoint(const DataArrayVi Line 379  Data::reshapeDataPoint(const DataArrayVi
379  Data  Data
380  Data::wherePositive() const  Data::wherePositive() const
381  {  {
382      profData->where++;
383    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
384  }  }
385    
386  Data  Data
387  Data::whereNegative() const  Data::whereNegative() const
388  {  {
389      profData->where++;
390    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
391  }  }
392    
393  Data  Data
394  Data::whereNonNegative() const  Data::whereNonNegative() const
395  {  {
396      profData->where++;
397    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
398  }  }
399    
400  Data  Data
401  Data::whereNonPositive() const  Data::whereNonPositive() const
402  {  {
403      profData->where++;
404    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
405  }  }
406    
407  Data  Data
408  Data::whereZero() const  Data::whereZero() const
409  {  {
410      profData->where++;
411    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));
412  }  }
413    
414  Data  Data
415  Data::whereNonZero() const  Data::whereNonZero() const
416  {  {
417      profData->where++;
418    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));
419  }  }
420    
421  Data  Data
422  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
423  {  {
424      profData->interpolate++;
425    return Data(*this,functionspace);    return Data(*this,functionspace);
426  }  }
427    
# Line 410  Data::probeInterpolation(const FunctionS Line 443  Data::probeInterpolation(const FunctionS
443  Data  Data
444  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
445  {  {
446      profData->grad++;
447    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
448      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
449    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
# Line 443  Data::getDataPointShape() const Line 477  Data::getDataPointShape() const
477    return getPointDataView().getShape();    return getPointDataView().getShape();
478  }  }
479    
480    void
481    Data::fillFromNumArray(const boost::python::numeric::array num_array)
482    {
483      //
484      // check rank
485      if (num_array.getrank()<getDataPointRank())
486          throw DataException("Rank of numarray does not match Data object rank");
487    
488      //
489      // check shape of num_array
490      for (int i=0; i<getDataPointRank(); i++) {
491        if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])
492           throw DataException("Shape of numarray does not match Data object rank");
493      }
494    
495      //
496      // make sure data is expanded:
497      if (!isExpanded()) {
498        expand();
499      }
500    
501      //
502      // and copy over
503      m_data->copyAll(num_array);
504    }
505    
506    const
507  boost::python::numeric::array  boost::python::numeric::array
508  Data::convertToNumArray()  Data::convertToNumArray()
509  {  {
510    //    //
511    // determine the current number of data points    // determine the total number of data points
512    int numSamples = getNumSamples();    int numSamples = getNumSamples();
513    int numDataPointsPerSample = getNumDataPointsPerSample();    int numDataPointsPerSample = getNumDataPointsPerSample();
514    int numDataPoints = numSamples * numDataPointsPerSample;    int numDataPoints = numSamples * numDataPointsPerSample;
# Line 499  Data::convertToNumArray() Line 560  Data::convertToNumArray()
560    int dataPoint = 0;    int dataPoint = 0;
561    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
562      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
   
563        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
   
564        if (dataPointRank==0) {        if (dataPointRank==0) {
565          numArray[dataPoint]=dataPointView();          numArray[dataPoint]=dataPointView();
566        }        }
   
567        if (dataPointRank==1) {        if (dataPointRank==1) {
568          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
569            numArray[dataPoint][i]=dataPointView(i);            numArray[dataPoint][i]=dataPointView(i);
570          }          }
571        }        }
   
572        if (dataPointRank==2) {        if (dataPointRank==2) {
573          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
574            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
# Line 519  Data::convertToNumArray() Line 576  Data::convertToNumArray()
576            }            }
577          }          }
578        }        }
   
579        if (dataPointRank==3) {        if (dataPointRank==3) {
580          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
581            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
# Line 529  Data::convertToNumArray() Line 585  Data::convertToNumArray()
585            }            }
586          }          }
587        }        }
   
588        if (dataPointRank==4) {        if (dataPointRank==4) {
589          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
590            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
# Line 541  Data::convertToNumArray() Line 596  Data::convertToNumArray()
596            }            }
597          }          }
598        }        }
   
599        dataPoint++;        dataPoint++;
600        }
601      }
602    
603      //
604      // return the loaded array
605      return numArray;
606    }
607    
608    const
609    boost::python::numeric::array
610    Data::convertToNumArrayFromSampleNo(int sampleNo)
611    {
612      //
613      // Check a valid sample number has been supplied
614      if (sampleNo >= getNumSamples()) {
615        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
616      }
617    
618      //
619      // determine the number of data points per sample
620      int numDataPointsPerSample = getNumDataPointsPerSample();
621    
622      //
623      // determine the rank and shape of each data point
624      int dataPointRank = getDataPointRank();
625      DataArrayView::ShapeType dataPointShape = getDataPointShape();
626    
627      //
628      // create the numeric array to be returned
629      boost::python::numeric::array numArray(0.0);
630    
631      //
632      // the rank of the returned numeric array will be the rank of
633      // the data points, plus one. Where the rank of the array is n,
634      // the last n-1 dimensions will be equal to the shape of the
635      // data points, whilst the first dimension will be equal to the
636      // total number of data points. Thus the array will consist of
637      // a serial vector of the data points.
638      int arrayRank = dataPointRank + 1;
639      DataArrayView::ShapeType arrayShape;
640      arrayShape.push_back(numDataPointsPerSample);
641      for (int d=0; d<dataPointRank; d++) {
642         arrayShape.push_back(dataPointShape[d]);
643      }
644    
645      //
646      // resize the numeric array to the shape just calculated
647      if (arrayRank==1) {
648        numArray.resize(arrayShape[0]);
649      }
650      if (arrayRank==2) {
651        numArray.resize(arrayShape[0],arrayShape[1]);
652      }
653      if (arrayRank==3) {
654        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
655      }
656      if (arrayRank==4) {
657        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
658      }
659      if (arrayRank==5) {
660        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
661      }
662    
663      //
664      // loop through each data point in turn, loading the values for that data point
665      // into the numeric array.
666      for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {
667        DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);
668        if (dataPointRank==0) {
669          numArray[dataPoint]=dataPointView();
670        }
671        if (dataPointRank==1) {
672          for (int i=0; i<dataPointShape[0]; i++) {
673            numArray[dataPoint][i]=dataPointView(i);
674          }
675        }
676        if (dataPointRank==2) {
677          for (int i=0; i<dataPointShape[0]; i++) {
678            for (int j=0; j<dataPointShape[1]; j++) {
679              numArray[dataPoint][i][j] = dataPointView(i,j);
680            }
681          }
682        }
683        if (dataPointRank==3) {
684          for (int i=0; i<dataPointShape[0]; i++) {
685            for (int j=0; j<dataPointShape[1]; j++) {
686              for (int k=0; k<dataPointShape[2]; k++) {
687                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
688              }
689            }
690          }
691        }
692        if (dataPointRank==4) {
693          for (int i=0; i<dataPointShape[0]; i++) {
694            for (int j=0; j<dataPointShape[1]; j++) {
695              for (int k=0; k<dataPointShape[2]; k++) {
696                for (int l=0; l<dataPointShape[3]; l++) {
697                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
698                }
699              }
700            }
701          }
702        }
703      }
704    
705      //
706      // return the loaded array
707      return numArray;
708    }
709    
710    const
711    boost::python::numeric::array
712    Data::convertToNumArrayFromDPNo(int sampleNo,
713                                    int dataPointNo)
714    {
715      //
716      // Check a valid sample number has been supplied
717      if (sampleNo >= getNumSamples()) {
718        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
719      }
720    
721      //
722      // Check a valid data point number has been supplied
723      if (dataPointNo >= getNumDataPointsPerSample()) {
724        throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");
725      }
726    
727      //
728      // determine the rank and shape of each data point
729      int dataPointRank = getDataPointRank();
730      DataArrayView::ShapeType dataPointShape = getDataPointShape();
731    
732      //
733      // create the numeric array to be returned
734      boost::python::numeric::array numArray(0.0);
735    
736      //
737      // the shape of the returned numeric array will be the same
738      // as that of the data point
739      int arrayRank = dataPointRank;
740      DataArrayView::ShapeType arrayShape = dataPointShape;
741    
742      //
743      // resize the numeric array to the shape just calculated
744      if (arrayRank==0) {
745        numArray.resize(1);
746      }
747      if (arrayRank==1) {
748        numArray.resize(arrayShape[0]);
749      }
750      if (arrayRank==2) {
751        numArray.resize(arrayShape[0],arrayShape[1]);
752      }
753      if (arrayRank==3) {
754        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
755      }
756      if (arrayRank==4) {
757        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
758      }
759    
760      //
761      // load the values for the data point into the numeric array.
762      DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
763      if (dataPointRank==0) {
764        numArray[0]=dataPointView();
765      }
766      if (dataPointRank==1) {
767        for (int i=0; i<dataPointShape[0]; i++) {
768          numArray[i]=dataPointView(i);
769        }
770      }
771      if (dataPointRank==2) {
772        for (int i=0; i<dataPointShape[0]; i++) {
773          for (int j=0; j<dataPointShape[1]; j++) {
774            numArray[i][j] = dataPointView(i,j);
775          }
776        }
777      }
778      if (dataPointRank==3) {
779        for (int i=0; i<dataPointShape[0]; i++) {
780          for (int j=0; j<dataPointShape[1]; j++) {
781            for (int k=0; k<dataPointShape[2]; k++) {
782              numArray[i][j][k]=dataPointView(i,j,k);
783            }
784          }
785        }
786      }
787      if (dataPointRank==4) {
788        for (int i=0; i<dataPointShape[0]; i++) {
789          for (int j=0; j<dataPointShape[1]; j++) {
790            for (int k=0; k<dataPointShape[2]; k++) {
791              for (int l=0; l<dataPointShape[3]; l++) {
792                numArray[i][j][k][l]=dataPointView(i,j,k,l);
793              }
794            }
795          }
796      }      }
797    }    }
798    
# Line 559  Data::integrate() const Line 808  Data::integrate() const
808    int rank = getDataPointRank();    int rank = getDataPointRank();
809    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
810    
811      profData->integrate++;
812    
813    //    //
814    // calculate the integral values    // calculate the integral values
815    vector<double> integrals(getDataPointSize());    vector<double> integrals(getDataPointSize());
# Line 622  Data::integrate() const Line 873  Data::integrate() const
873  Data  Data
874  Data::sin() const  Data::sin() const
875  {  {
876      profData->unary++;
877    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
878  }  }
879    
880  Data  Data
881  Data::cos() const  Data::cos() const
882  {  {
883      profData->unary++;
884    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
885  }  }
886    
887  Data  Data
888  Data::tan() const  Data::tan() const
889  {  {
890      profData->unary++;
891    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
892  }  }
893    
894  Data  Data
895  Data::log() const  Data::log() const
896  {  {
897      profData->unary++;
898    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
899  }  }
900    
901  Data  Data
902  Data::ln() const  Data::ln() const
903  {  {
904      profData->unary++;
905    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
906  }  }
907    
908  Data  Data
909  Data::sign() const  Data::sign() const
910  {  {
911      profData->unary++;
912    return escript::unaryOp(*this,escript::fsign);    return escript::unaryOp(*this,escript::fsign);
913  }  }
914    
915  Data  Data
916  Data::abs() const  Data::abs() const
917  {  {
918      profData->unary++;
919    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
920  }  }
921    
922  Data  Data
923  Data::neg() const  Data::neg() const
924  {  {
925      profData->unary++;
926    return escript::unaryOp(*this,negate<double>());    return escript::unaryOp(*this,negate<double>());
927  }  }
928    
929  Data  Data
930  Data::pos() const  Data::pos() const
931  {  {
932    return (*this);    profData->unary++;
933      Data result;
934      // perform a deep copy
935      result.copy(*this);
936      return result;
937  }  }
938    
939  Data  Data
940  Data::exp() const  Data::exp() const
941  {  {
942      profData->unary++;
943    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
944  }  }
945    
946  Data  Data
947  Data::sqrt() const  Data::sqrt() const
948  {  {
949      profData->unary++;
950    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
951  }  }
952    
953  double  double
954  Data::Lsup() const  Data::Lsup() const
955  {  {
956      profData->reduction1++;
957    //    //
958    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
959    return algorithm(DataAlgorithmAdapter<AbsMax>(0));    AbsMax abs_max_func;
960      return algorithm(abs_max_func,0);
961  }  }
962    
963  double  double
964  Data::Linf() const  Data::Linf() const
965  {  {
966      profData->reduction1++;
967    //    //
968    // set the initial absolute minimum value to max double    // set the initial absolute minimum value to max double
969    return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));    AbsMin abs_min_func;
970      return algorithm(abs_min_func,numeric_limits<double>::max());
971  }  }
972    
973  double  double
974  Data::sup() const  Data::sup() const
975  {  {
976      profData->reduction1++;
977    //    //
978    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
979    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
980      return algorithm(fmax_func,numeric_limits<double>::max()*-1);
981  }  }
982    
983  double  double
984  Data::inf() const  Data::inf() const
985  {  {
986      profData->reduction1++;
987    //    //
988    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
989    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
990      return algorithm(fmin_func,numeric_limits<double>::max());
991  }  }
992    
993  Data  Data
994  Data::maxval() const  Data::maxval() const
995  {  {
996      profData->reduction2++;
997    //    //
998    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
999    return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1000      return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1001  }  }
1002    
1003  Data  Data
1004  Data::minval() const  Data::minval() const
1005  {  {
1006      profData->reduction2++;
1007    //    //
1008    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1009    return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1010      return dp_algorithm(fmin_func,numeric_limits<double>::max());
1011  }  }
1012    
1013  Data  Data
1014  Data::length() const  Data::length() const
1015  {  {
1016    return dp_algorithm(DataAlgorithmAdapter<Length>(0));    profData->reduction2++;
1017      Length len_func;
1018      return dp_algorithm(len_func,0);
1019  }  }
1020    
1021  Data  Data
1022  Data::trace() const  Data::trace() const
1023  {  {
1024    return dp_algorithm(DataAlgorithmAdapter<Trace>(0));    profData->reduction2++;
1025      Trace trace_func;
1026      return dp_algorithm(trace_func,0);
1027  }  }
1028    
1029  Data  Data
1030  Data::transpose(int axis) const  Data::transpose(int axis) const
1031  {  {
1032      profData->reduction2++;
1033    // not implemented    // not implemented
1034    throw DataException("Error - Data::transpose not implemented yet.");    throw DataException("Error - Data::transpose not implemented yet.");
1035    return Data();    return Data();
1036  }  }
1037    
1038    const boost::python::tuple
1039    Data::mindp() const
1040    {
1041      // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1042      // abort (for unknown reasons) if there are openmp directives with it in the
1043      // surrounding function
1044    
1045      int SampleNo;
1046      int DataPointNo;
1047    
1048      calc_mindp(SampleNo,DataPointNo);
1049    
1050      return make_tuple(SampleNo,DataPointNo);
1051    }
1052    
1053    void
1054    Data::calc_mindp(int& SampleNo,
1055                     int& DataPointNo) const
1056    {
1057      int i,j;
1058      int lowi=0,lowj=0;
1059      double min=numeric_limits<double>::max();
1060    
1061      Data temp=minval();
1062    
1063      int numSamples=temp.getNumSamples();
1064      int numDPPSample=temp.getNumDataPointsPerSample();
1065    
1066      double next,local_min;
1067      int local_lowi,local_lowj;
1068    
1069      #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1070      {
1071        local_min=min;
1072        #pragma omp for private(i,j) schedule(static)
1073        for (i=0; i<numSamples; i++) {
1074          for (j=0; j<numDPPSample; j++) {
1075            next=temp.getDataPoint(i,j)();
1076            if (next<local_min) {
1077              local_min=next;
1078              local_lowi=i;
1079              local_lowj=j;
1080            }
1081          }
1082        }
1083        #pragma omp critical
1084        if (local_min<min) {
1085          min=local_min;
1086          lowi=local_lowi;
1087          lowj=local_lowj;
1088        }
1089      }
1090    
1091      SampleNo = lowi;
1092      DataPointNo = lowj;
1093    }
1094    
1095  void  void
1096  Data::saveDX(std::string fileName) const  Data::saveDX(std::string fileName) const
1097  {  {
# Line 770  Data::saveVTK(std::string fileName) cons Line 1109  Data::saveVTK(std::string fileName) cons
1109  Data&  Data&
1110  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1111  {  {
1112      profData->binary++;
1113    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1114    return (*this);    return (*this);
1115  }  }
# Line 777  Data::operator+=(const Data& right) Line 1117  Data::operator+=(const Data& right)
1117  Data&  Data&
1118  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1119  {  {
1120      profData->binary++;
1121    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1122    return (*this);    return (*this);
1123  }  }
# Line 784  Data::operator+=(const boost::python::ob Line 1125  Data::operator+=(const boost::python::ob
1125  Data&  Data&
1126  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1127  {  {
1128      profData->binary++;
1129    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1130    return (*this);    return (*this);
1131  }  }
# Line 791  Data::operator-=(const Data& right) Line 1133  Data::operator-=(const Data& right)
1133  Data&  Data&
1134  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1135  {  {
1136      profData->binary++;
1137    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1138    return (*this);    return (*this);
1139  }  }
# Line 798  Data::operator-=(const boost::python::ob Line 1141  Data::operator-=(const boost::python::ob
1141  Data&  Data&
1142  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1143  {  {
1144      profData->binary++;
1145    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1146    return (*this);    return (*this);
1147  }  }
# Line 805  Data::operator*=(const Data& right) Line 1149  Data::operator*=(const Data& right)
1149  Data&  Data&
1150  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1151  {  {
1152      profData->binary++;
1153    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1154    return (*this);    return (*this);
1155  }  }
# Line 812  Data::operator*=(const boost::python::ob Line 1157  Data::operator*=(const boost::python::ob
1157  Data&  Data&
1158  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1159  {  {
1160      profData->binary++;
1161    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1162    return (*this);    return (*this);
1163  }  }
# Line 819  Data::operator/=(const Data& right) Line 1165  Data::operator/=(const Data& right)
1165  Data&  Data&
1166  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1167  {  {
1168      profData->binary++;
1169    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1170    return (*this);    return (*this);
1171  }  }
# Line 826  Data::operator/=(const boost::python::ob Line 1173  Data::operator/=(const boost::python::ob
1173  Data  Data
1174  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1175  {  {
1176      profData->binary++;
1177    Data result;    Data result;
1178    result.copy(*this);    result.copy(*this);
1179    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 835  Data::powO(const boost::python::object& Line 1183  Data::powO(const boost::python::object&
1183  Data  Data
1184  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1185  {  {
1186      profData->binary++;
1187    Data result;    Data result;
1188    result.copy(*this);    result.copy(*this);
1189    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 842  Data::powD(const Data& right) const Line 1191  Data::powD(const Data& right) const
1191  }  }
1192    
1193  //  //
1194  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1195  Data  Data
1196  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1197  {  {
# Line 855  escript::operator+(const Data& left, con Line 1204  escript::operator+(const Data& left, con
1204  }  }
1205    
1206  //  //
1207  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1208  Data  Data
1209  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1210  {  {
# Line 868  escript::operator-(const Data& left, con Line 1217  escript::operator-(const Data& left, con
1217  }  }
1218    
1219  //  //
1220  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1221  Data  Data
1222  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
1223  {  {
# Line 881  escript::operator*(const Data& left, con Line 1230  escript::operator*(const Data& left, con
1230  }  }
1231    
1232  //  //
1233  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1234  Data  Data
1235  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
1236  {  {
# Line 894  escript::operator/(const Data& left, con Line 1243  escript::operator/(const Data& left, con
1243  }  }
1244    
1245  //  //
1246  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1247  Data  Data
1248  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1249  {  {
# Line 910  escript::operator+(const Data& left, con Line 1259  escript::operator+(const Data& left, con
1259  }  }
1260    
1261  //  //
1262  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1263  Data  Data
1264  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1265  {  {
# Line 926  escript::operator-(const Data& left, con Line 1275  escript::operator-(const Data& left, con
1275  }  }
1276    
1277  //  //
1278  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1279  Data  Data
1280  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1281  {  {
# Line 942  escript::operator*(const Data& left, con Line 1291  escript::operator*(const Data& left, con
1291  }  }
1292    
1293  //  //
1294  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1295  Data  Data
1296  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1297  {  {
# Line 958  escript::operator/(const Data& left, con Line 1307  escript::operator/(const Data& left, con
1307  }  }
1308    
1309  //  //
1310  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1311  Data  Data
1312  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1313  {  {
# Line 971  escript::operator+(const boost::python:: Line 1320  escript::operator+(const boost::python::
1320  }  }
1321    
1322  //  //
1323  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1324  Data  Data
1325  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1326  {  {
# Line 984  escript::operator-(const boost::python:: Line 1333  escript::operator-(const boost::python::
1333  }  }
1334    
1335  //  //
1336  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1337  Data  Data
1338  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1339  {  {
# Line 997  escript::operator*(const boost::python:: Line 1346  escript::operator*(const boost::python::
1346  }  }
1347    
1348  //  //
1349  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1350  Data  Data
1351  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1352  {  {
# Line 1010  escript::operator/(const boost::python:: Line 1359  escript::operator/(const boost::python::
1359  }  }
1360    
1361  //  //
 // NOTE: It is essential to specify the namepsace this operator belongs to  
1362  //bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1363  //{  //{
1364  //  /*  //  /*
# Line 1072  Data::getItem(const boost::python::objec Line 1420  Data::getItem(const boost::python::objec
1420  Data  Data
1421  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1422  {  {
1423      profData->slicing++;
1424    return Data(*this,region);    return Data(*this,region);
1425  }  }
1426    
# Line 1088  Data::setItemD(const boost::python::obje Line 1437  Data::setItemD(const boost::python::obje
1437                 const Data& value)                 const Data& value)
1438  {  {
1439    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1440    
1441    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1442    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1443      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
# Line 1103  void Line 1453  void
1453  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1454                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
1455  {  {
1456      profData->slicing++;
1457    Data tempValue(value);    Data tempValue(value);
1458    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1459    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1159  Data::setTaggedValue(int tagKey, Line 1510  Data::setTaggedValue(int tagKey,
1510  }  }
1511    
1512  void  void
1513    Data::setTaggedValueFromCPP(int tagKey,
1514                                const DataArrayView& value)
1515    {
1516      //
1517      // Ensure underlying data object is of type DataTagged
1518      tag();
1519    
1520      if (!isTagged()) {
1521        throw DataException("Error - DataTagged conversion failed!!");
1522      }
1523                                                                                                                  
1524      //
1525      // Call DataAbstract::setTaggedValue
1526      m_data->setTaggedValue(tagKey,value);
1527    }
1528    
1529    int
1530    Data::getTagNumber(int dpno)
1531    {
1532      return m_data->getTagNumber(dpno);
1533    }
1534    
1535    void
1536  Data::setRefValue(int ref,  Data::setRefValue(int ref,
1537                    const boost::python::numeric::array& value)                    const boost::python::numeric::array& value)
1538  {  {
# Line 1197  Data::getRefValue(int ref, Line 1571  Data::getRefValue(int ref,
1571    DataArrayView valueView = valueDataArray.getView();    DataArrayView valueView = valueDataArray.getView();
1572    
1573    if (rank==0) {    if (rank==0) {
1574      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");        boost::python::numeric::array temp_numArray(valueView());
1575          value = temp_numArray;
1576    }    }
1577    if (rank==1) {    if (rank==1) {
1578      for (int i=0; i < shape[0]; i++) {      for (int i=0; i < shape[0]; i++) {
# Line 1205  Data::getRefValue(int ref, Line 1580  Data::getRefValue(int ref,
1580      }      }
1581    }    }
1582    if (rank==2) {    if (rank==2) {
1583      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
1584          for (int j=0; j < shape[1]; j++) {
1585            value[i][j] = valueView(i,j);
1586          }
1587        }
1588    }    }
1589    if (rank==3) {    if (rank==3) {
1590      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
1591          for (int j=0; j < shape[1]; j++) {
1592            for (int k=0; k < shape[2]; k++) {
1593              value[i][j][k] = valueView(i,j,k);
1594            }
1595          }
1596        }
1597    }    }
1598    if (rank==4) {    if (rank==4) {
1599      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
1600          for (int j=0; j < shape[1]; j++) {
1601            for (int k=0; k < shape[2]; k++) {
1602              for (int l=0; l < shape[3]; l++) {
1603                value[i][j][k][l] = valueView(i,j,k,l);
1604              }
1605            }
1606          }
1607        }
1608    }    }
1609    
1610  }  }
1611    
 /*  
 Note: this version removed for now. Not needed, and breaks escript.cpp  
1612  void  void
1613  Data::setTaggedValue(int tagKey,  Data::archiveData(const std::string fileName)
                      const DataArrayView& value)  
1614  {  {
1615      cout << "Archiving Data object to: " << fileName << endl;
1616    
1617    //    //
1618    // Ensure underlying data object is of type DataTagged    // Determine type of this Data object
1619    tag();    int dataType = -1;
1620    
1621    if (!isTagged()) {    if (isEmpty()) {
1622      throw DataException("Error - DataTagged conversion failed!!");      dataType = 0;
1623        cout << "\tdataType: DataEmpty" << endl;
1624    }    }
1625                                                                                                                    if (isConstant()) {
1626        dataType = 1;
1627        cout << "\tdataType: DataConstant" << endl;
1628      }
1629      if (isTagged()) {
1630        dataType = 2;
1631        cout << "\tdataType: DataTagged" << endl;
1632      }
1633      if (isExpanded()) {
1634        dataType = 3;
1635        cout << "\tdataType: DataExpanded" << endl;
1636      }
1637    
1638      if (dataType == -1) {
1639        throw DataException("archiveData Error: undefined dataType");
1640      }
1641    
1642    //    //
1643    // Call DataAbstract::setTaggedValue    // Collect data items common to all Data types
1644    m_data->setTaggedValue(tagKey,value);    int noSamples = getNumSamples();
1645      int noDPPSample = getNumDataPointsPerSample();
1646      int functionSpaceType = getFunctionSpace().getTypeCode();
1647      int dataPointRank = getDataPointRank();
1648      int dataPointSize = getDataPointSize();
1649      int dataLength = getLength();
1650      DataArrayView::ShapeType dataPointShape = getDataPointShape();
1651      int referenceNumbers[noSamples];
1652      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1653        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1654      }
1655      int tagNumbers[noSamples];
1656      if (isTagged()) {
1657        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1658          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1659        }
1660      }
1661    
1662      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1663      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1664      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1665    
1666      //
1667      // Flatten Shape to an array of integers suitable for writing to file
1668      int flatShape[4] = {0,0,0,0};
1669      cout << "\tshape: < ";
1670      for (int dim=0; dim<dataPointRank; dim++) {
1671        flatShape[dim] = dataPointShape[dim];
1672        cout << dataPointShape[dim] << " ";
1673      }
1674      cout << ">" << endl;
1675    
1676      //
1677      // Open archive file
1678      ofstream archiveFile;
1679      archiveFile.open(fileName.data(), ios::out);
1680    
1681      if (!archiveFile.good()) {
1682        throw DataException("archiveData Error: problem opening archive file");
1683      }
1684    
1685      //
1686      // Write common data items to archive file
1687      archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1688      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1689      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1690      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1691      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1692      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1693      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1694      for (int dim = 0; dim < 4; dim++) {
1695        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1696      }
1697      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1698        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1699      }
1700      if (isTagged()) {
1701        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1702          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1703        }
1704      }
1705    
1706      if (!archiveFile.good()) {
1707        throw DataException("archiveData Error: problem writing to archive file");
1708      }
1709    
1710      //
1711      // Archive underlying data values for each Data type
1712      int noValues;
1713      switch (dataType) {
1714        case 0:
1715          // DataEmpty
1716          noValues = 0;
1717          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1718          cout << "\tnoValues: " << noValues << endl;
1719          break;
1720        case 1:
1721          // DataConstant
1722          noValues = m_data->getLength();
1723          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1724          cout << "\tnoValues: " << noValues << endl;
1725          if (m_data->archiveData(archiveFile,noValues)) {
1726            throw DataException("archiveData Error: problem writing data to archive file");
1727          }
1728          break;
1729        case 2:
1730          // DataTagged
1731          noValues = m_data->getLength();
1732          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1733          cout << "\tnoValues: " << noValues << endl;
1734          if (m_data->archiveData(archiveFile,noValues)) {
1735            throw DataException("archiveData Error: problem writing data to archive file");
1736          }
1737          break;
1738        case 3:
1739          // DataExpanded
1740          noValues = m_data->getLength();
1741          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1742          cout << "\tnoValues: " << noValues << endl;
1743          if (m_data->archiveData(archiveFile,noValues)) {
1744            throw DataException("archiveData Error: problem writing data to archive file");
1745          }
1746          break;
1747      }
1748    
1749      if (!archiveFile.good()) {
1750        throw DataException("archiveData Error: problem writing data to archive file");
1751      }
1752    
1753      //
1754      // Close archive file
1755      archiveFile.close();
1756    
1757      if (!archiveFile.good()) {
1758        throw DataException("archiveData Error: problem closing archive file");
1759      }
1760    
1761    }
1762    
1763    void
1764    Data::extractData(const std::string fileName,
1765                      const FunctionSpace& fspace)
1766    {
1767      //
1768      // Can only extract Data to an object which is initially DataEmpty
1769      if (!isEmpty()) {
1770        throw DataException("extractData Error: can only extract to DataEmpty object");
1771      }
1772    
1773      cout << "Extracting Data object from: " << fileName << endl;
1774    
1775      int dataType;
1776      int noSamples;
1777      int noDPPSample;
1778      int functionSpaceType;
1779      int dataPointRank;
1780      int dataPointSize;
1781      int dataLength;
1782      DataArrayView::ShapeType dataPointShape;
1783      int flatShape[4];
1784    
1785      //
1786      // Open the archive file
1787      ifstream archiveFile;
1788      archiveFile.open(fileName.data(), ios::in);
1789    
1790      if (!archiveFile.good()) {
1791        throw DataException("extractData Error: problem opening archive file");
1792      }
1793    
1794      //
1795      // Read common data items from archive file
1796      archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1797      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1798      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1799      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1800      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1801      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1802      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1803      for (int dim = 0; dim < 4; dim++) {
1804        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1805        if (flatShape[dim]>0) {
1806          dataPointShape.push_back(flatShape[dim]);
1807        }
1808      }
1809      int referenceNumbers[noSamples];
1810      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1811        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1812      }
1813      int tagNumbers[noSamples];
1814      if (dataType==2) {
1815        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1816          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1817        }
1818      }
1819    
1820      if (!archiveFile.good()) {
1821        throw DataException("extractData Error: problem reading from archive file");
1822      }
1823    
1824      //
1825      // Verify the values just read from the archive file
1826      switch (dataType) {
1827        case 0:
1828          cout << "\tdataType: DataEmpty" << endl;
1829          break;
1830        case 1:
1831          cout << "\tdataType: DataConstant" << endl;
1832          break;
1833        case 2:
1834          cout << "\tdataType: DataTagged" << endl;
1835          break;
1836        case 3:
1837          cout << "\tdataType: DataExpanded" << endl;
1838          break;
1839        default:
1840          throw DataException("extractData Error: undefined dataType read from archive file");
1841          break;
1842      }
1843    
1844      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1845      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1846      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1847      cout << "\tshape: < ";
1848      for (int dim = 0; dim < dataPointRank; dim++) {
1849        cout << dataPointShape[dim] << " ";
1850      }
1851      cout << ">" << endl;
1852    
1853      //
1854      // Verify that supplied FunctionSpace object is compatible with this Data object.
1855      if ( (fspace.getTypeCode()!=functionSpaceType) ||
1856           (fspace.getNumSamples()!=noSamples) ||
1857           (fspace.getNumDPPSample()!=noDPPSample)
1858         ) {
1859        throw DataException("extractData Error: incompatible FunctionSpace");
1860      }
1861      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1862        if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
1863          throw DataException("extractData Error: incompatible FunctionSpace");
1864        }
1865      }
1866      if (dataType==2) {
1867        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1868          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
1869            throw DataException("extractData Error: incompatible FunctionSpace");
1870          }
1871        }
1872      }
1873    
1874      //
1875      // Construct a DataVector to hold underlying data values
1876      DataVector dataVec(dataLength);
1877    
1878      //
1879      // Load this DataVector with the appropriate values
1880      int noValues;
1881      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
1882      cout << "\tnoValues: " << noValues << endl;
1883      switch (dataType) {
1884        case 0:
1885          // DataEmpty
1886          if (noValues != 0) {
1887            throw DataException("extractData Error: problem reading data from archive file");
1888          }
1889          break;
1890        case 1:
1891          // DataConstant
1892          if (dataVec.extractData(archiveFile,noValues)) {
1893            throw DataException("extractData Error: problem reading data from archive file");
1894          }
1895          break;
1896        case 2:
1897          // DataTagged
1898          if (dataVec.extractData(archiveFile,noValues)) {
1899            throw DataException("extractData Error: problem reading data from archive file");
1900          }
1901          break;
1902        case 3:
1903          // DataExpanded
1904          if (dataVec.extractData(archiveFile,noValues)) {
1905            throw DataException("extractData Error: problem reading data from archive file");
1906          }
1907          break;
1908      }
1909    
1910      if (!archiveFile.good()) {
1911        throw DataException("extractData Error: problem reading from archive file");
1912      }
1913    
1914      //
1915      // Close archive file
1916      archiveFile.close();
1917    
1918      if (!archiveFile.good()) {
1919        throw DataException("extractData Error: problem closing archive file");
1920      }
1921    
1922      //
1923      // Construct an appropriate Data object
1924      DataAbstract* tempData;
1925      switch (dataType) {
1926        case 0:
1927          // DataEmpty
1928          tempData=new DataEmpty();
1929          break;
1930        case 1:
1931          // DataConstant
1932          tempData=new DataConstant(fspace,dataPointShape,dataVec);
1933          break;
1934        case 2:
1935          // DataTagged
1936          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
1937          break;
1938        case 3:
1939          // DataExpanded
1940          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
1941          break;
1942      }
1943      shared_ptr<DataAbstract> temp_data(tempData);
1944      m_data=temp_data;
1945  }  }
 */  
1946    
1947  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
1948  {  {

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

  ViewVC Help
Powered by ViewVC 1.1.26