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

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

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

revision 97 by jgs, Tue Dec 14 05:39:33 2004 UTC revision 123 by jgs, Fri Jul 8 04:08:13 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    const
481    boost::python::numeric::array
482    Data::convertToNumArray()
483    {
484      //
485      // determine the total number of data points
486      int numSamples = getNumSamples();
487      int numDataPointsPerSample = getNumDataPointsPerSample();
488      int numDataPoints = numSamples * numDataPointsPerSample;
489    
490      //
491      // determine the rank and shape of each data point
492      int dataPointRank = getDataPointRank();
493      DataArrayView::ShapeType dataPointShape = getDataPointShape();
494    
495      //
496      // create the numeric array to be returned
497      boost::python::numeric::array numArray(0.0);
498    
499      //
500      // the rank of the returned numeric array will be the rank of
501      // the data points, plus one. Where the rank of the array is n,
502      // the last n-1 dimensions will be equal to the shape of the
503      // data points, whilst the first dimension will be equal to the
504      // total number of data points. Thus the array will consist of
505      // a serial vector of the data points.
506      int arrayRank = dataPointRank + 1;
507      DataArrayView::ShapeType arrayShape;
508      arrayShape.push_back(numDataPoints);
509      for (int d=0; d<dataPointRank; d++) {
510         arrayShape.push_back(dataPointShape[d]);
511      }
512    
513      //
514      // resize the numeric array to the shape just calculated
515      if (arrayRank==1) {
516        numArray.resize(arrayShape[0]);
517      }
518      if (arrayRank==2) {
519        numArray.resize(arrayShape[0],arrayShape[1]);
520      }
521      if (arrayRank==3) {
522        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
523      }
524      if (arrayRank==4) {
525        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
526      }
527      if (arrayRank==5) {
528        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
529      }
530    
531      //
532      // loop through each data point in turn, loading the values for that data point
533      // into the numeric array.
534      int dataPoint = 0;
535      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
536        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
537          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
538          if (dataPointRank==0) {
539            numArray[dataPoint]=dataPointView();
540          }
541          if (dataPointRank==1) {
542            for (int i=0; i<dataPointShape[0]; i++) {
543              numArray[dataPoint][i]=dataPointView(i);
544            }
545          }
546          if (dataPointRank==2) {
547            for (int i=0; i<dataPointShape[0]; i++) {
548              for (int j=0; j<dataPointShape[1]; j++) {
549                numArray[dataPoint][i][j] = dataPointView(i,j);
550              }
551            }
552          }
553          if (dataPointRank==3) {
554            for (int i=0; i<dataPointShape[0]; i++) {
555              for (int j=0; j<dataPointShape[1]; j++) {
556                for (int k=0; k<dataPointShape[2]; k++) {
557                  numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
558                }
559              }
560            }
561          }
562          if (dataPointRank==4) {
563            for (int i=0; i<dataPointShape[0]; i++) {
564              for (int j=0; j<dataPointShape[1]; j++) {
565                for (int k=0; k<dataPointShape[2]; k++) {
566                  for (int l=0; l<dataPointShape[3]; l++) {
567                    numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
568                  }
569                }
570              }
571            }
572          }
573          dataPoint++;
574        }
575      }
576    
577      //
578      // return the loaded array
579      return numArray;
580    }
581    
582    const
583    boost::python::numeric::array
584    Data::convertToNumArrayFromSampleNo(int sampleNo)
585    {
586      //
587      // Check a valid sample number has been supplied
588      if (sampleNo >= getNumSamples()) {
589        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
590      }
591    
592      //
593      // determine the number of data points per sample
594      int numDataPointsPerSample = getNumDataPointsPerSample();
595    
596      //
597      // determine the rank and shape of each data point
598      int dataPointRank = getDataPointRank();
599      DataArrayView::ShapeType dataPointShape = getDataPointShape();
600    
601      //
602      // create the numeric array to be returned
603      boost::python::numeric::array numArray(0.0);
604    
605      //
606      // the rank of the returned numeric array will be the rank of
607      // the data points, plus one. Where the rank of the array is n,
608      // the last n-1 dimensions will be equal to the shape of the
609      // data points, whilst the first dimension will be equal to the
610      // total number of data points. Thus the array will consist of
611      // a serial vector of the data points.
612      int arrayRank = dataPointRank + 1;
613      DataArrayView::ShapeType arrayShape;
614      arrayShape.push_back(numDataPointsPerSample);
615      for (int d=0; d<dataPointRank; d++) {
616         arrayShape.push_back(dataPointShape[d]);
617      }
618    
619      //
620      // resize the numeric array to the shape just calculated
621      if (arrayRank==1) {
622        numArray.resize(arrayShape[0]);
623      }
624      if (arrayRank==2) {
625        numArray.resize(arrayShape[0],arrayShape[1]);
626      }
627      if (arrayRank==3) {
628        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
629      }
630      if (arrayRank==4) {
631        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
632      }
633      if (arrayRank==5) {
634        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
635      }
636    
637      //
638      // loop through each data point in turn, loading the values for that data point
639      // into the numeric array.
640      for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {
641        DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);
642        if (dataPointRank==0) {
643          numArray[dataPoint]=dataPointView();
644        }
645        if (dataPointRank==1) {
646          for (int i=0; i<dataPointShape[0]; i++) {
647            numArray[dataPoint][i]=dataPointView(i);
648          }
649        }
650        if (dataPointRank==2) {
651          for (int i=0; i<dataPointShape[0]; i++) {
652            for (int j=0; j<dataPointShape[1]; j++) {
653              numArray[dataPoint][i][j] = dataPointView(i,j);
654            }
655          }
656        }
657        if (dataPointRank==3) {
658          for (int i=0; i<dataPointShape[0]; i++) {
659            for (int j=0; j<dataPointShape[1]; j++) {
660              for (int k=0; k<dataPointShape[2]; k++) {
661                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
662              }
663            }
664          }
665        }
666        if (dataPointRank==4) {
667          for (int i=0; i<dataPointShape[0]; i++) {
668            for (int j=0; j<dataPointShape[1]; j++) {
669              for (int k=0; k<dataPointShape[2]; k++) {
670                for (int l=0; l<dataPointShape[3]; l++) {
671                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
672                }
673              }
674            }
675          }
676        }
677      }
678    
679      //
680      // return the loaded array
681      return numArray;
682    }
683    
684    const
685    boost::python::numeric::array
686    Data::convertToNumArrayFromDPNo(int sampleNo,
687                                    int dataPointNo)
688    {
689      //
690      // Check a valid sample number has been supplied
691      if (sampleNo >= getNumSamples()) {
692        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
693      }
694    
695      //
696      // Check a valid data point number has been supplied
697      if (dataPointNo >= getNumDataPointsPerSample()) {
698        throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");
699      }
700    
701      //
702      // determine the rank and shape of each data point
703      int dataPointRank = getDataPointRank();
704      DataArrayView::ShapeType dataPointShape = getDataPointShape();
705    
706      //
707      // create the numeric array to be returned
708      boost::python::numeric::array numArray(0.0);
709    
710      //
711      // the shape of the returned numeric array will be the same
712      // as that of the data point
713      int arrayRank = dataPointRank;
714      DataArrayView::ShapeType arrayShape = dataPointShape;
715    
716      //
717      // resize the numeric array to the shape just calculated
718      if (arrayRank==0) {
719        numArray.resize(1);
720      }
721      if (arrayRank==1) {
722        numArray.resize(arrayShape[0]);
723      }
724      if (arrayRank==2) {
725        numArray.resize(arrayShape[0],arrayShape[1]);
726      }
727      if (arrayRank==3) {
728        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
729      }
730      if (arrayRank==4) {
731        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
732      }
733    
734      //
735      // load the values for the data point into the numeric array.
736      DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
737      if (dataPointRank==0) {
738        numArray[0]=dataPointView();
739      }
740      if (dataPointRank==1) {
741        for (int i=0; i<dataPointShape[0]; i++) {
742          numArray[i]=dataPointView(i);
743        }
744      }
745      if (dataPointRank==2) {
746        for (int i=0; i<dataPointShape[0]; i++) {
747          for (int j=0; j<dataPointShape[1]; j++) {
748            numArray[i][j] = dataPointView(i,j);
749          }
750        }
751      }
752      if (dataPointRank==3) {
753        for (int i=0; i<dataPointShape[0]; i++) {
754          for (int j=0; j<dataPointShape[1]; j++) {
755            for (int k=0; k<dataPointShape[2]; k++) {
756              numArray[i][j][k]=dataPointView(i,j,k);
757            }
758          }
759        }
760      }
761      if (dataPointRank==4) {
762        for (int i=0; i<dataPointShape[0]; i++) {
763          for (int j=0; j<dataPointShape[1]; j++) {
764            for (int k=0; k<dataPointShape[2]; k++) {
765              for (int l=0; l<dataPointShape[3]; l++) {
766                numArray[i][j][k][l]=dataPointView(i,j,k,l);
767              }
768            }
769          }
770        }
771      }
772    
773      //
774      // return the loaded array
775      return numArray;
776    }
777    
778  boost::python::numeric::array  boost::python::numeric::array
779  Data::integrate() const  Data::integrate() const
780  {  {
# Line 450  Data::integrate() const Line 782  Data::integrate() const
782    int rank = getDataPointRank();    int rank = getDataPointRank();
783    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
784    
785      profData->integrate++;
786    
787    //    //
788    // calculate the integral values    // calculate the integral values
789    vector<double> integrals(getDataPointSize());    vector<double> integrals(getDataPointSize());
# Line 460  Data::integrate() const Line 794  Data::integrate() const
794    // and load the array with the integral values    // and load the array with the integral values
795    boost::python::numeric::array bp_array(1.0);    boost::python::numeric::array bp_array(1.0);
796    if (rank==0) {    if (rank==0) {
797        bp_array.resize(1);
798      index = 0;      index = 0;
799      bp_array[0] = integrals[index];      bp_array[0] = integrals[index];
800    }    }
# Line 512  Data::integrate() const Line 847  Data::integrate() const
847  Data  Data
848  Data::sin() const  Data::sin() const
849  {  {
850      profData->unary++;
851    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
852  }  }
853    
854  Data  Data
855  Data::cos() const  Data::cos() const
856  {  {
857      profData->unary++;
858    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
859  }  }
860    
861  Data  Data
862  Data::tan() const  Data::tan() const
863  {  {
864      profData->unary++;
865    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
866  }  }
867    
868  Data  Data
869  Data::log() const  Data::log() const
870  {  {
871      profData->unary++;
872    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
873  }  }
874    
875  Data  Data
876  Data::ln() const  Data::ln() const
877  {  {
878      profData->unary++;
879    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
880  }  }
881    
882    Data
883    Data::sign() const
884    {
885      profData->unary++;
886      return escript::unaryOp(*this,escript::fsign);
887    }
888    
889    Data
890    Data::abs() const
891    {
892      profData->unary++;
893      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
894    }
895    
896    Data
897    Data::neg() const
898    {
899      profData->unary++;
900      return escript::unaryOp(*this,negate<double>());
901    }
902    
903    Data
904    Data::pos() const
905    {
906      profData->unary++;
907      return (*this);
908    }
909    
910    Data
911    Data::exp() const
912    {
913      profData->unary++;
914      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
915    }
916    
917    Data
918    Data::sqrt() const
919    {
920      profData->unary++;
921      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
922    }
923    
924  double  double
925  Data::Lsup() const  Data::Lsup() const
926  {  {
927      profData->reduction1++;
928    //    //
929    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
930    return algorithm(DataAlgorithmAdapter<AbsMax>(0));    return algorithm(DataAlgorithmAdapter<AbsMax>(0));
931  }  }
932    
933  double  double
934    Data::Linf() const
935    {
936      profData->reduction1++;
937      //
938      // set the initial absolute minimum value to max double
939      return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));
940    }
941    
942    double
943  Data::sup() const  Data::sup() const
944  {  {
945      profData->reduction1++;
946    //    //
947    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
948    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::min()));    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
949  }  }
950    
951  double  double
952  Data::inf() const  Data::inf() const
953  {  {
954      profData->reduction1++;
955    //    //
956    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
957    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
# Line 566  Data::inf() const Line 960  Data::inf() const
960  Data  Data
961  Data::maxval() const  Data::maxval() const
962  {  {
963    // not implemented - will use dp_algorithm    profData->reduction2++;
964    return (*this);    //
965      // set the initial maximum value to min possible double
966      return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
967  }  }
968    
969  Data  Data
970  Data::minval() const  Data::minval() const
971  {  {
972    // not implemented - will use dp_algorithm    profData->reduction2++;
973    return (*this);    //
974      // set the initial minimum value to max possible double
975      return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
976  }  }
977    
978  Data  Data
979  Data::length() const  Data::length() const
980  {  {
981    // not implemented - will use dp_algorithm    profData->reduction2++;
982    return (*this);    return dp_algorithm(DataAlgorithmAdapter<Length>(0));
983  }  }
984    
985  Data  Data
986  Data::trace() const  Data::trace() const
987  {  {
988    // not implemented - will use dp_algorithm    profData->reduction2++;
989    return (*this);    return dp_algorithm(DataAlgorithmAdapter<Trace>(0));
990  }  }
991    
992  Data  Data
993  Data::transpose(int axis) const  Data::transpose(int axis) const
994  {  {
995      profData->reduction2++;
996    // not implemented    // not implemented
997    return (*this);    throw DataException("Error - Data::transpose not implemented yet.");
998      return Data();
999  }  }
1000    
1001  Data  const boost::python::tuple
1002  Data::sign() const  Data::mindp() const
1003  {  {
1004    return escript::unaryOp(*this,escript::fsign);    Data temp=minval();
 }  
1005    
1006  Data    int numSamples=temp.getNumSamples();
1007  Data::abs() const    int numDPPSample=temp.getNumDataPointsPerSample();
 {  
   return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);  
 }  
1008    
1009  Data    int i,j,lowi=0,lowj=0;
1010  Data::neg() const    double min=numeric_limits<double>::max();
 {  
   return escript::unaryOp(*this,negate<double>());  
 }  
1011    
1012  Data    for (i=0; i<numSamples; i++) {
1013  Data::pos() const      for (j=0; j<numDPPSample; j++) {
1014  {        double next=temp.getDataPoint(i,j)();
1015    return (*this);        if (next<min) {
1016            min=next;
1017            lowi=i;
1018            lowj=j;
1019          }
1020        }
1021      }
1022    
1023      return make_tuple(lowi,lowj);
1024  }  }
1025    
1026  Data  void
1027  Data::exp() const  Data::saveDX(std::string fileName) const
1028  {  {
1029    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    getDomain().saveDX(fileName,*this);
1030      return;
1031  }  }
1032    
1033  Data  void
1034  Data::sqrt() const  Data::saveVTK(std::string fileName) const
1035  {  {
1036    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    getDomain().saveVTK(fileName,*this);
1037      return;
1038  }  }
1039    
1040  Data&  Data&
1041  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1042  {  {
1043      profData->binary++;
1044    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1045    return (*this);    return (*this);
1046  }  }
# Line 644  Data::operator+=(const Data& right) Line 1048  Data::operator+=(const Data& right)
1048  Data&  Data&
1049  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1050  {  {
1051      profData->binary++;
1052    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1053    return (*this);    return (*this);
1054  }  }
# Line 651  Data::operator+=(const boost::python::ob Line 1056  Data::operator+=(const boost::python::ob
1056  Data&  Data&
1057  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1058  {  {
1059      profData->binary++;
1060    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1061    return (*this);    return (*this);
1062  }  }
# Line 658  Data::operator-=(const Data& right) Line 1064  Data::operator-=(const Data& right)
1064  Data&  Data&
1065  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1066  {  {
1067      profData->binary++;
1068    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1069    return (*this);    return (*this);
1070  }  }
# Line 665  Data::operator-=(const boost::python::ob Line 1072  Data::operator-=(const boost::python::ob
1072  Data&  Data&
1073  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1074  {  {
1075      profData->binary++;
1076    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1077    return (*this);    return (*this);
1078  }  }
# Line 672  Data::operator*=(const Data& right) Line 1080  Data::operator*=(const Data& right)
1080  Data&  Data&
1081  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1082  {  {
1083      profData->binary++;
1084    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1085    return (*this);    return (*this);
1086  }  }
# Line 679  Data::operator*=(const boost::python::ob Line 1088  Data::operator*=(const boost::python::ob
1088  Data&  Data&
1089  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1090  {  {
1091      profData->binary++;
1092    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1093    return (*this);    return (*this);
1094  }  }
# Line 686  Data::operator/=(const Data& right) Line 1096  Data::operator/=(const Data& right)
1096  Data&  Data&
1097  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1098  {  {
1099      profData->binary++;
1100    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1101    return (*this);    return (*this);
1102  }  }
# Line 693  Data::operator/=(const boost::python::ob Line 1104  Data::operator/=(const boost::python::ob
1104  Data  Data
1105  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1106  {  {
1107      profData->binary++;
1108    Data result;    Data result;
1109    result.copy(*this);    result.copy(*this);
1110    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 702  Data::powO(const boost::python::object& Line 1114  Data::powO(const boost::python::object&
1114  Data  Data
1115  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1116  {  {
1117      profData->binary++;
1118    Data result;    Data result;
1119    result.copy(*this);    result.copy(*this);
1120    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 709  Data::powD(const Data& right) const Line 1122  Data::powD(const Data& right) const
1122  }  }
1123    
1124  //  //
1125  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1126  Data  Data
1127  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1128  {  {
# Line 722  escript::operator+(const Data& left, con Line 1135  escript::operator+(const Data& left, con
1135  }  }
1136    
1137  //  //
1138  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1139  Data  Data
1140  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1141  {  {
# Line 735  escript::operator-(const Data& left, con Line 1148  escript::operator-(const Data& left, con
1148  }  }
1149    
1150  //  //
1151  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1152  Data  Data
1153  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
1154  {  {
# Line 748  escript::operator*(const Data& left, con Line 1161  escript::operator*(const Data& left, con
1161  }  }
1162    
1163  //  //
1164  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1165  Data  Data
1166  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
1167  {  {
# Line 761  escript::operator/(const Data& left, con Line 1174  escript::operator/(const Data& left, con
1174  }  }
1175    
1176  //  //
1177  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1178  Data  Data
1179  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1180  {  {
# Line 777  escript::operator+(const Data& left, con Line 1190  escript::operator+(const Data& left, con
1190  }  }
1191    
1192  //  //
1193  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1194  Data  Data
1195  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1196  {  {
# Line 793  escript::operator-(const Data& left, con Line 1206  escript::operator-(const Data& left, con
1206  }  }
1207    
1208  //  //
1209  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1210  Data  Data
1211  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1212  {  {
# Line 809  escript::operator*(const Data& left, con Line 1222  escript::operator*(const Data& left, con
1222  }  }
1223    
1224  //  //
1225  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1226  Data  Data
1227  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1228  {  {
# Line 825  escript::operator/(const Data& left, con Line 1238  escript::operator/(const Data& left, con
1238  }  }
1239    
1240  //  //
1241  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1242  Data  Data
1243  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1244  {  {
# Line 838  escript::operator+(const boost::python:: Line 1251  escript::operator+(const boost::python::
1251  }  }
1252    
1253  //  //
1254  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1255  Data  Data
1256  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1257  {  {
# Line 851  escript::operator-(const boost::python:: Line 1264  escript::operator-(const boost::python::
1264  }  }
1265    
1266  //  //
1267  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1268  Data  Data
1269  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1270  {  {
# Line 864  escript::operator*(const boost::python:: Line 1277  escript::operator*(const boost::python::
1277  }  }
1278    
1279  //  //
1280  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1281  Data  Data
1282  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1283  {  {
# Line 877  escript::operator/(const boost::python:: Line 1290  escript::operator/(const boost::python::
1290  }  }
1291    
1292  //  //
 // NOTE: It is essential to specify the namepsace this operator belongs to  
1293  //bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1294  //{  //{
1295  //  /*  //  /*
# Line 939  Data::getItem(const boost::python::objec Line 1351  Data::getItem(const boost::python::objec
1351  Data  Data
1352  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1353  {  {
1354      profData->slicing++;
1355    return Data(*this,region);    return Data(*this,region);
1356  }  }
1357    
# Line 955  Data::setItemD(const boost::python::obje Line 1368  Data::setItemD(const boost::python::obje
1368                 const Data& value)                 const Data& value)
1369  {  {
1370    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1371    
1372    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1373    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1374      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
1375    }    }
1376    setSlice(value,slice_region);    if (getFunctionSpace()!=value.getFunctionSpace()) {
1377         setSlice(Data(value,getFunctionSpace()),slice_region);
1378      } else {
1379         setSlice(value,slice_region);
1380      }
1381  }  }
1382    
1383  void  void
1384  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1385                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
1386  {  {
1387      profData->slicing++;
1388    Data tempValue(value);    Data tempValue(value);
1389    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1390    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1021  Data::setTaggedValue(int tagKey, Line 1440  Data::setTaggedValue(int tagKey,
1440    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
1441  }  }
1442    
 /*  
 Note: this version removed for now. Not needed, and breaks escript.cpp  
1443  void  void
1444  Data::setTaggedValue(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
1445                       const DataArrayView& value)                              const DataArrayView& value)
1446  {  {
1447    //    //
1448    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
# Line 1039  Data::setTaggedValue(int tagKey, Line 1456  Data::setTaggedValue(int tagKey,
1456    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
1457    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
1458  }  }
1459  */  
1460    void
1461    Data::setRefValue(int ref,
1462                      const boost::python::numeric::array& value)
1463    {
1464      //
1465      // Construct DataArray from boost::python::object input value
1466      DataArray valueDataArray(value);
1467    
1468      //
1469      // Call DataAbstract::setRefValue
1470      m_data->setRefValue(ref,valueDataArray);
1471    }
1472    
1473    void
1474    Data::getRefValue(int ref,
1475                      boost::python::numeric::array& value)
1476    {
1477      //
1478      // Construct DataArray for boost::python::object return value
1479      DataArray valueDataArray(value);
1480    
1481      //
1482      // Load DataArray with values from data-points specified by ref
1483      m_data->getRefValue(ref,valueDataArray);
1484    
1485      //
1486      // Load values from valueDataArray into return numarray
1487    
1488      // extract the shape of the numarray
1489      int rank = value.getrank();
1490      DataArrayView::ShapeType shape;
1491      for (int i=0; i < rank; i++) {
1492        shape.push_back(extract<int>(value.getshape()[i]));
1493      }
1494    
1495      // and load the numarray with the data from the DataArray
1496      DataArrayView valueView = valueDataArray.getView();
1497    
1498      if (rank==0) {
1499        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1500      }
1501      if (rank==1) {
1502        for (int i=0; i < shape[0]; i++) {
1503          value[i] = valueView(i);
1504        }
1505      }
1506      if (rank==2) {
1507        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1508      }
1509      if (rank==3) {
1510        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1511      }
1512      if (rank==4) {
1513        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1514      }
1515    
1516    }
1517    
1518    void
1519    Data::archiveData(const std::string fileName)
1520    {
1521      cout << "Archiving Data object to: " << fileName << endl;
1522    
1523      //
1524      // Determine type of this Data object
1525      int dataType = -1;
1526    
1527      if (isEmpty()) {
1528        dataType = 0;
1529        cout << "\tdataType: DataEmpty" << endl;
1530      }
1531      if (isConstant()) {
1532        dataType = 1;
1533        cout << "\tdataType: DataConstant" << endl;
1534      }
1535      if (isTagged()) {
1536        dataType = 2;
1537        cout << "\tdataType: DataTagged" << endl;
1538      }
1539      if (isExpanded()) {
1540        dataType = 3;
1541        cout << "\tdataType: DataExpanded" << endl;
1542      }
1543    
1544      if (dataType == -1) {
1545        throw DataException("archiveData Error: undefined dataType");
1546      }
1547    
1548      //
1549      // Collect data items common to all Data types
1550      int noSamples = getNumSamples();
1551      int noDPPSample = getNumDataPointsPerSample();
1552      int functionSpaceType = getFunctionSpace().getTypeCode();
1553      int dataPointRank = getDataPointRank();
1554      int dataPointSize = getDataPointSize();
1555      int dataLength = getLength();
1556      DataArrayView::ShapeType dataPointShape = getDataPointShape();
1557      int referenceNumbers[noSamples];
1558      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1559        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1560      }
1561      int tagNumbers[noSamples];
1562      if (isTagged()) {
1563        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1564          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1565        }
1566      }
1567    
1568      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1569      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1570      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1571    
1572      //
1573      // Flatten Shape to an array of integers suitable for writing to file
1574      int flatShape[4] = {0,0,0,0};
1575      cout << "\tshape: < ";
1576      for (int dim=0; dim<dataPointRank; dim++) {
1577        flatShape[dim] = dataPointShape[dim];
1578        cout << dataPointShape[dim] << " ";
1579      }
1580      cout << ">" << endl;
1581    
1582      //
1583      // Open archive file
1584      ofstream archiveFile;
1585      archiveFile.open(fileName.data(), ios::out);
1586    
1587      if (!archiveFile.good()) {
1588        throw DataException("archiveData Error: problem opening archive file");
1589      }
1590    
1591      //
1592      // Write common data items to archive file
1593      archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1594      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1595      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1596      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1597      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1598      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1599      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1600      for (int dim = 0; dim < 4; dim++) {
1601        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1602      }
1603      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1604        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1605      }
1606      if (isTagged()) {
1607        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1608          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1609        }
1610      }
1611    
1612      if (!archiveFile.good()) {
1613        throw DataException("archiveData Error: problem writing to archive file");
1614      }
1615    
1616      //
1617      // Archive underlying data values for each Data type
1618      int noValues;
1619      switch (dataType) {
1620        case 0:
1621          // DataEmpty
1622          noValues = 0;
1623          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1624          cout << "\tnoValues: " << noValues << endl;
1625          break;
1626        case 1:
1627          // DataConstant
1628          noValues = m_data->getLength();
1629          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1630          cout << "\tnoValues: " << noValues << endl;
1631          if (m_data->archiveData(archiveFile,noValues)) {
1632            throw DataException("archiveData Error: problem writing data to archive file");
1633          }
1634          break;
1635        case 2:
1636          // DataTagged
1637          noValues = m_data->getLength();
1638          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1639          cout << "\tnoValues: " << noValues << endl;
1640          if (m_data->archiveData(archiveFile,noValues)) {
1641            throw DataException("archiveData Error: problem writing data to archive file");
1642          }
1643          break;
1644        case 3:
1645          // DataExpanded
1646          noValues = m_data->getLength();
1647          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1648          cout << "\tnoValues: " << noValues << endl;
1649          if (m_data->archiveData(archiveFile,noValues)) {
1650            throw DataException("archiveData Error: problem writing data to archive file");
1651          }
1652          break;
1653      }
1654    
1655      if (!archiveFile.good()) {
1656        throw DataException("archiveData Error: problem writing data to archive file");
1657      }
1658    
1659      //
1660      // Close archive file
1661      archiveFile.close();
1662    
1663      if (!archiveFile.good()) {
1664        throw DataException("archiveData Error: problem closing archive file");
1665      }
1666    
1667    }
1668    
1669    void
1670    Data::extractData(const std::string fileName,
1671                      const FunctionSpace& fspace)
1672    {
1673      //
1674      // Can only extract Data to an object which is initially DataEmpty
1675      if (!isEmpty()) {
1676        throw DataException("extractData Error: can only extract to DataEmpty object");
1677      }
1678    
1679      cout << "Extracting Data object from: " << fileName << endl;
1680    
1681      int dataType;
1682      int noSamples;
1683      int noDPPSample;
1684      int functionSpaceType;
1685      int dataPointRank;
1686      int dataPointSize;
1687      int dataLength;
1688      DataArrayView::ShapeType dataPointShape;
1689      int flatShape[4];
1690    
1691      //
1692      // Open the archive file
1693      ifstream archiveFile;
1694      archiveFile.open(fileName.data(), ios::in);
1695    
1696      if (!archiveFile.good()) {
1697        throw DataException("extractData Error: problem opening archive file");
1698      }
1699    
1700      //
1701      // Read common data items from archive file
1702      archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1703      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1704      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1705      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1706      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1707      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1708      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1709      for (int dim = 0; dim < 4; dim++) {
1710        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1711        if (flatShape[dim]>0) {
1712          dataPointShape.push_back(flatShape[dim]);
1713        }
1714      }
1715      int referenceNumbers[noSamples];
1716      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1717        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1718      }
1719      int tagNumbers[noSamples];
1720      if (dataType==2) {
1721        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1722          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1723        }
1724      }
1725    
1726      if (!archiveFile.good()) {
1727        throw DataException("extractData Error: problem reading from archive file");
1728      }
1729    
1730      //
1731      // Verify the values just read from the archive file
1732      switch (dataType) {
1733        case 0:
1734          cout << "\tdataType: DataEmpty" << endl;
1735          break;
1736        case 1:
1737          cout << "\tdataType: DataConstant" << endl;
1738          break;
1739        case 2:
1740          cout << "\tdataType: DataTagged" << endl;
1741          break;
1742        case 3:
1743          cout << "\tdataType: DataExpanded" << endl;
1744          break;
1745        default:
1746          throw DataException("extractData Error: undefined dataType read from archive file");
1747          break;
1748      }
1749    
1750      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1751      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1752      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1753      cout << "\tshape: < ";
1754      for (int dim = 0; dim < dataPointRank; dim++) {
1755        cout << dataPointShape[dim] << " ";
1756      }
1757      cout << ">" << endl;
1758    
1759      //
1760      // Verify that supplied FunctionSpace object is compatible with this Data object.
1761      if ( (fspace.getTypeCode()!=functionSpaceType) ||
1762           (fspace.getNumSamples()!=noSamples) ||
1763           (fspace.getNumDPPSample()!=noDPPSample)
1764         ) {
1765        throw DataException("extractData Error: incompatible FunctionSpace");
1766      }
1767      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1768        if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
1769          throw DataException("extractData Error: incompatible FunctionSpace");
1770        }
1771      }
1772      if (dataType==2) {
1773        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1774          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
1775            throw DataException("extractData Error: incompatible FunctionSpace");
1776          }
1777        }
1778      }
1779    
1780      //
1781      // Construct a DataVector to hold underlying data values
1782      DataVector dataVec(dataLength);
1783    
1784      //
1785      // Load this DataVector with the appropriate values
1786      int noValues;
1787      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
1788      cout << "\tnoValues: " << noValues << endl;
1789      switch (dataType) {
1790        case 0:
1791          // DataEmpty
1792          if (noValues != 0) {
1793            throw DataException("extractData Error: problem reading data from archive file");
1794          }
1795          break;
1796        case 1:
1797          // DataConstant
1798          if (dataVec.extractData(archiveFile,noValues)) {
1799            throw DataException("extractData Error: problem reading data from archive file");
1800          }
1801          break;
1802        case 2:
1803          // DataTagged
1804          if (dataVec.extractData(archiveFile,noValues)) {
1805            throw DataException("extractData Error: problem reading data from archive file");
1806          }
1807          break;
1808        case 3:
1809          // DataExpanded
1810          if (dataVec.extractData(archiveFile,noValues)) {
1811            throw DataException("extractData Error: problem reading data from archive file");
1812          }
1813          break;
1814      }
1815    
1816      if (!archiveFile.good()) {
1817        throw DataException("extractData Error: problem reading from archive file");
1818      }
1819    
1820      //
1821      // Close archive file
1822      archiveFile.close();
1823    
1824      if (!archiveFile.good()) {
1825        throw DataException("extractData Error: problem closing archive file");
1826      }
1827    
1828      //
1829      // Construct an appropriate Data object
1830      DataAbstract* tempData;
1831      switch (dataType) {
1832        case 0:
1833          // DataEmpty
1834          tempData=new DataEmpty();
1835          break;
1836        case 1:
1837          // DataConstant
1838          tempData=new DataConstant(fspace,dataPointShape,dataVec);
1839          break;
1840        case 2:
1841          // DataTagged
1842          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
1843          break;
1844        case 3:
1845          // DataExpanded
1846          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
1847          break;
1848      }
1849      shared_ptr<DataAbstract> temp_data(tempData);
1850      m_data=temp_data;
1851    
1852    }
1853    
1854  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
1855  {  {

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

  ViewVC Help
Powered by ViewVC 1.1.26