/[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 150 by jgs, Thu Sep 15 03:44:45 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::asin() const
896    {
897      profData->unary++;
898      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
899    }
900    
901    Data
902    Data::acos() const
903    {
904      profData->unary++;
905      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
906    }
907    
908    Data
909    Data::atan() const
910    {
911      profData->unary++;
912      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
913    }
914    
915    Data
916    Data::sinh() const
917    {
918      profData->unary++;
919      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
920    }
921    
922    Data
923    Data::cosh() const
924    {
925      profData->unary++;
926      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
927    }
928    
929    Data
930    Data::tanh() const
931    {
932      profData->unary++;
933      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
934    }
935    
936    Data
937    Data::asinh() const
938    {
939      profData->unary++;
940      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
941    }
942    
943    Data
944    Data::acosh() const
945    {
946      profData->unary++;
947      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
948    }
949    
950    Data
951    Data::atanh() const
952    {
953      profData->unary++;
954      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
955    }
956    
957    Data
958  Data::log() const  Data::log() const
959  {  {
960      profData->unary++;
961    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
962  }  }
963    
964  Data  Data
965  Data::ln() const  Data::ln() const
966  {  {
967      profData->unary++;
968    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
969  }  }
970    
971  Data  Data
972  Data::sign() const  Data::sign() const
973  {  {
974      profData->unary++;
975    return escript::unaryOp(*this,escript::fsign);    return escript::unaryOp(*this,escript::fsign);
976  }  }
977    
978  Data  Data
979  Data::abs() const  Data::abs() const
980  {  {
981      profData->unary++;
982    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
983  }  }
984    
985  Data  Data
986  Data::neg() const  Data::neg() const
987  {  {
988      profData->unary++;
989    return escript::unaryOp(*this,negate<double>());    return escript::unaryOp(*this,negate<double>());
990  }  }
991    
992  Data  Data
993  Data::pos() const  Data::pos() const
994  {  {
995    return (*this);    profData->unary++;
996      Data result;
997      // perform a deep copy
998      result.copy(*this);
999      return result;
1000  }  }
1001    
1002  Data  Data
1003  Data::exp() const  Data::exp() const
1004  {  {
1005      profData->unary++;
1006    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1007  }  }
1008    
1009  Data  Data
1010  Data::sqrt() const  Data::sqrt() const
1011  {  {
1012      profData->unary++;
1013    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1014  }  }
1015    
1016  double  double
1017  Data::Lsup() const  Data::Lsup() const
1018  {  {
1019      profData->reduction1++;
1020    //    //
1021    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
1022    return algorithm(DataAlgorithmAdapter<AbsMax>(0));    AbsMax abs_max_func;
1023      return algorithm(abs_max_func,0);
1024  }  }
1025    
1026  double  double
1027  Data::Linf() const  Data::Linf() const
1028  {  {
1029      profData->reduction1++;
1030    //    //
1031    // set the initial absolute minimum value to max double    // set the initial absolute minimum value to max double
1032    return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));    AbsMin abs_min_func;
1033      return algorithm(abs_min_func,numeric_limits<double>::max());
1034  }  }
1035    
1036  double  double
1037  Data::sup() const  Data::sup() const
1038  {  {
1039      profData->reduction1++;
1040    //    //
1041    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1042    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1043      return algorithm(fmax_func,numeric_limits<double>::max()*-1);
1044  }  }
1045    
1046  double  double
1047  Data::inf() const  Data::inf() const
1048  {  {
1049      profData->reduction1++;
1050    //    //
1051    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1052    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1053      return algorithm(fmin_func,numeric_limits<double>::max());
1054  }  }
1055    
1056  Data  Data
1057  Data::maxval() const  Data::maxval() const
1058  {  {
1059      profData->reduction2++;
1060    //    //
1061    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1062    return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1063      return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1064  }  }
1065    
1066  Data  Data
1067  Data::minval() const  Data::minval() const
1068  {  {
1069      profData->reduction2++;
1070    //    //
1071    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1072    return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1073      return dp_algorithm(fmin_func,numeric_limits<double>::max());
1074  }  }
1075    
1076  Data  Data
1077  Data::length() const  Data::length() const
1078  {  {
1079    return dp_algorithm(DataAlgorithmAdapter<Length>(0));    profData->reduction2++;
1080      Length len_func;
1081      return dp_algorithm(len_func,0);
1082  }  }
1083    
1084  Data  Data
1085  Data::trace() const  Data::trace() const
1086  {  {
1087    return dp_algorithm(DataAlgorithmAdapter<Trace>(0));    profData->reduction2++;
1088      Trace trace_func;
1089      return dp_algorithm(trace_func,0);
1090  }  }
1091    
1092  Data  Data
1093  Data::transpose(int axis) const  Data::transpose(int axis) const
1094  {  {
1095      profData->reduction2++;
1096    // not implemented    // not implemented
1097    throw DataException("Error - Data::transpose not implemented yet.");    throw DataException("Error - Data::transpose not implemented yet.");
1098    return Data();    return Data();
1099  }  }
1100    
1101    const boost::python::tuple
1102    Data::mindp() const
1103    {
1104      // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1105      // abort (for unknown reasons) if there are openmp directives with it in the
1106      // surrounding function
1107    
1108      int SampleNo;
1109      int DataPointNo;
1110    
1111      calc_mindp(SampleNo,DataPointNo);
1112    
1113      return make_tuple(SampleNo,DataPointNo);
1114    }
1115    
1116    void
1117    Data::calc_mindp(int& SampleNo,
1118                     int& DataPointNo) const
1119    {
1120      int i,j;
1121      int lowi=0,lowj=0;
1122      double min=numeric_limits<double>::max();
1123    
1124      Data temp=minval();
1125    
1126      int numSamples=temp.getNumSamples();
1127      int numDPPSample=temp.getNumDataPointsPerSample();
1128    
1129      double next,local_min;
1130      int local_lowi,local_lowj;
1131    
1132      #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1133      {
1134        local_min=min;
1135        #pragma omp for private(i,j) schedule(static)
1136        for (i=0; i<numSamples; i++) {
1137          for (j=0; j<numDPPSample; j++) {
1138            next=temp.getDataPoint(i,j)();
1139            if (next<local_min) {
1140              local_min=next;
1141              local_lowi=i;
1142              local_lowj=j;
1143            }
1144          }
1145        }
1146        #pragma omp critical
1147        if (local_min<min) {
1148          min=local_min;
1149          lowi=local_lowi;
1150          lowj=local_lowj;
1151        }
1152      }
1153    
1154      SampleNo = lowi;
1155      DataPointNo = lowj;
1156    }
1157    
1158  void  void
1159  Data::saveDX(std::string fileName) const  Data::saveDX(std::string fileName) const
1160  {  {
# Line 770  Data::saveVTK(std::string fileName) cons Line 1172  Data::saveVTK(std::string fileName) cons
1172  Data&  Data&
1173  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1174  {  {
1175      profData->binary++;
1176    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1177    return (*this);    return (*this);
1178  }  }
# Line 777  Data::operator+=(const Data& right) Line 1180  Data::operator+=(const Data& right)
1180  Data&  Data&
1181  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1182  {  {
1183      profData->binary++;
1184    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1185    return (*this);    return (*this);
1186  }  }
# Line 784  Data::operator+=(const boost::python::ob Line 1188  Data::operator+=(const boost::python::ob
1188  Data&  Data&
1189  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1190  {  {
1191      profData->binary++;
1192    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1193    return (*this);    return (*this);
1194  }  }
# Line 791  Data::operator-=(const Data& right) Line 1196  Data::operator-=(const Data& right)
1196  Data&  Data&
1197  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1198  {  {
1199      profData->binary++;
1200    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1201    return (*this);    return (*this);
1202  }  }
# Line 798  Data::operator-=(const boost::python::ob Line 1204  Data::operator-=(const boost::python::ob
1204  Data&  Data&
1205  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1206  {  {
1207      profData->binary++;
1208    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1209    return (*this);    return (*this);
1210  }  }
# Line 805  Data::operator*=(const Data& right) Line 1212  Data::operator*=(const Data& right)
1212  Data&  Data&
1213  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1214  {  {
1215      profData->binary++;
1216    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1217    return (*this);    return (*this);
1218  }  }
# Line 812  Data::operator*=(const boost::python::ob Line 1220  Data::operator*=(const boost::python::ob
1220  Data&  Data&
1221  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1222  {  {
1223      profData->binary++;
1224    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1225    return (*this);    return (*this);
1226  }  }
# Line 819  Data::operator/=(const Data& right) Line 1228  Data::operator/=(const Data& right)
1228  Data&  Data&
1229  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1230  {  {
1231      profData->binary++;
1232    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1233    return (*this);    return (*this);
1234  }  }
# Line 826  Data::operator/=(const boost::python::ob Line 1236  Data::operator/=(const boost::python::ob
1236  Data  Data
1237  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1238  {  {
1239      profData->binary++;
1240    Data result;    Data result;
1241    result.copy(*this);    result.copy(*this);
1242    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 835  Data::powO(const boost::python::object& Line 1246  Data::powO(const boost::python::object&
1246  Data  Data
1247  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1248  {  {
1249      profData->binary++;
1250    Data result;    Data result;
1251    result.copy(*this);    result.copy(*this);
1252    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 842  Data::powD(const Data& right) const Line 1254  Data::powD(const Data& right) const
1254  }  }
1255    
1256  //  //
1257  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1258  Data  Data
1259  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1260  {  {
# Line 855  escript::operator+(const Data& left, con Line 1267  escript::operator+(const Data& left, con
1267  }  }
1268    
1269  //  //
1270  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1271  Data  Data
1272  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1273  {  {
# Line 868  escript::operator-(const Data& left, con Line 1280  escript::operator-(const Data& left, con
1280  }  }
1281    
1282  //  //
1283  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1284  Data  Data
1285  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
1286  {  {
# Line 881  escript::operator*(const Data& left, con Line 1293  escript::operator*(const Data& left, con
1293  }  }
1294    
1295  //  //
1296  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1297  Data  Data
1298  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
1299  {  {
# Line 894  escript::operator/(const Data& left, con Line 1306  escript::operator/(const Data& left, con
1306  }  }
1307    
1308  //  //
1309  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1310  Data  Data
1311  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1312  {  {
# Line 910  escript::operator+(const Data& left, con Line 1322  escript::operator+(const Data& left, con
1322  }  }
1323    
1324  //  //
1325  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1326  Data  Data
1327  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1328  {  {
# Line 926  escript::operator-(const Data& left, con Line 1338  escript::operator-(const Data& left, con
1338  }  }
1339    
1340  //  //
1341  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1342  Data  Data
1343  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1344  {  {
# Line 942  escript::operator*(const Data& left, con Line 1354  escript::operator*(const Data& left, con
1354  }  }
1355    
1356  //  //
1357  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1358  Data  Data
1359  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1360  {  {
# Line 958  escript::operator/(const Data& left, con Line 1370  escript::operator/(const Data& left, con
1370  }  }
1371    
1372  //  //
1373  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1374  Data  Data
1375  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1376  {  {
# Line 971  escript::operator+(const boost::python:: Line 1383  escript::operator+(const boost::python::
1383  }  }
1384    
1385  //  //
1386  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1387  Data  Data
1388  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1389  {  {
# Line 984  escript::operator-(const boost::python:: Line 1396  escript::operator-(const boost::python::
1396  }  }
1397    
1398  //  //
1399  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1400  Data  Data
1401  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1402  {  {
# Line 997  escript::operator*(const boost::python:: Line 1409  escript::operator*(const boost::python::
1409  }  }
1410    
1411  //  //
1412  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1413  Data  Data
1414  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1415  {  {
# Line 1010  escript::operator/(const boost::python:: Line 1422  escript::operator/(const boost::python::
1422  }  }
1423    
1424  //  //
 // NOTE: It is essential to specify the namepsace this operator belongs to  
1425  //bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1426  //{  //{
1427  //  /*  //  /*
# Line 1072  Data::getItem(const boost::python::objec Line 1483  Data::getItem(const boost::python::objec
1483  Data  Data
1484  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1485  {  {
1486      profData->slicing++;
1487    return Data(*this,region);    return Data(*this,region);
1488  }  }
1489    
# Line 1088  Data::setItemD(const boost::python::obje Line 1500  Data::setItemD(const boost::python::obje
1500                 const Data& value)                 const Data& value)
1501  {  {
1502    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1503    
1504    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1505    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1506      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
# Line 1103  void Line 1516  void
1516  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1517                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
1518  {  {
1519      profData->slicing++;
1520    Data tempValue(value);    Data tempValue(value);
1521    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1522    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1159  Data::setTaggedValue(int tagKey, Line 1573  Data::setTaggedValue(int tagKey,
1573  }  }
1574    
1575  void  void
1576    Data::setTaggedValueFromCPP(int tagKey,
1577                                const DataArrayView& value)
1578    {
1579      //
1580      // Ensure underlying data object is of type DataTagged
1581      tag();
1582    
1583      if (!isTagged()) {
1584        throw DataException("Error - DataTagged conversion failed!!");
1585      }
1586                                                                                                                  
1587      //
1588      // Call DataAbstract::setTaggedValue
1589      m_data->setTaggedValue(tagKey,value);
1590    }
1591    
1592    int
1593    Data::getTagNumber(int dpno)
1594    {
1595      return m_data->getTagNumber(dpno);
1596    }
1597    
1598    void
1599  Data::setRefValue(int ref,  Data::setRefValue(int ref,
1600                    const boost::python::numeric::array& value)                    const boost::python::numeric::array& value)
1601  {  {
# Line 1197  Data::getRefValue(int ref, Line 1634  Data::getRefValue(int ref,
1634    DataArrayView valueView = valueDataArray.getView();    DataArrayView valueView = valueDataArray.getView();
1635    
1636    if (rank==0) {    if (rank==0) {
1637      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");        boost::python::numeric::array temp_numArray(valueView());
1638          value = temp_numArray;
1639    }    }
1640    if (rank==1) {    if (rank==1) {
1641      for (int i=0; i < shape[0]; i++) {      for (int i=0; i < shape[0]; i++) {
# Line 1205  Data::getRefValue(int ref, Line 1643  Data::getRefValue(int ref,
1643      }      }
1644    }    }
1645    if (rank==2) {    if (rank==2) {
1646      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
1647          for (int j=0; j < shape[1]; j++) {
1648            value[i][j] = valueView(i,j);
1649          }
1650        }
1651    }    }
1652    if (rank==3) {    if (rank==3) {
1653      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
1654          for (int j=0; j < shape[1]; j++) {
1655            for (int k=0; k < shape[2]; k++) {
1656              value[i][j][k] = valueView(i,j,k);
1657            }
1658          }
1659        }
1660    }    }
1661    if (rank==4) {    if (rank==4) {
1662      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
1663          for (int j=0; j < shape[1]; j++) {
1664            for (int k=0; k < shape[2]; k++) {
1665              for (int l=0; l < shape[3]; l++) {
1666                value[i][j][k][l] = valueView(i,j,k,l);
1667              }
1668            }
1669          }
1670        }
1671    }    }
1672    
1673  }  }
1674    
 /*  
 Note: this version removed for now. Not needed, and breaks escript.cpp  
1675  void  void
1676  Data::setTaggedValue(int tagKey,  Data::archiveData(const std::string fileName)
                      const DataArrayView& value)  
1677  {  {
1678      cout << "Archiving Data object to: " << fileName << endl;
1679    
1680    //    //
1681    // Ensure underlying data object is of type DataTagged    // Determine type of this Data object
1682    tag();    int dataType = -1;
1683    
1684    if (!isTagged()) {    if (isEmpty()) {
1685      throw DataException("Error - DataTagged conversion failed!!");      dataType = 0;
1686        cout << "\tdataType: DataEmpty" << endl;
1687    }    }
1688                                                                                                                    if (isConstant()) {
1689        dataType = 1;
1690        cout << "\tdataType: DataConstant" << endl;
1691      }
1692      if (isTagged()) {
1693        dataType = 2;
1694        cout << "\tdataType: DataTagged" << endl;
1695      }
1696      if (isExpanded()) {
1697        dataType = 3;
1698        cout << "\tdataType: DataExpanded" << endl;
1699      }
1700    
1701      if (dataType == -1) {
1702        throw DataException("archiveData Error: undefined dataType");
1703      }
1704    
1705    //    //
1706    // Call DataAbstract::setTaggedValue    // Collect data items common to all Data types
1707    m_data->setTaggedValue(tagKey,value);    int noSamples = getNumSamples();
1708      int noDPPSample = getNumDataPointsPerSample();
1709      int functionSpaceType = getFunctionSpace().getTypeCode();
1710      int dataPointRank = getDataPointRank();
1711      int dataPointSize = getDataPointSize();
1712      int dataLength = getLength();
1713      DataArrayView::ShapeType dataPointShape = getDataPointShape();
1714      int referenceNumbers[noSamples];
1715      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1716        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1717      }
1718      int tagNumbers[noSamples];
1719      if (isTagged()) {
1720        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1721          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1722        }
1723      }
1724    
1725      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1726      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1727      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1728    
1729      //
1730      // Flatten Shape to an array of integers suitable for writing to file
1731      int flatShape[4] = {0,0,0,0};
1732      cout << "\tshape: < ";
1733      for (int dim=0; dim<dataPointRank; dim++) {
1734        flatShape[dim] = dataPointShape[dim];
1735        cout << dataPointShape[dim] << " ";
1736      }
1737      cout << ">" << endl;
1738    
1739      //
1740      // Open archive file
1741      ofstream archiveFile;
1742      archiveFile.open(fileName.data(), ios::out);
1743    
1744      if (!archiveFile.good()) {
1745        throw DataException("archiveData Error: problem opening archive file");
1746      }
1747    
1748      //
1749      // Write common data items to archive file
1750      archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1751      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1752      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1753      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1754      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1755      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1756      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1757      for (int dim = 0; dim < 4; dim++) {
1758        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1759      }
1760      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1761        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1762      }
1763      if (isTagged()) {
1764        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1765          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1766        }
1767      }
1768    
1769      if (!archiveFile.good()) {
1770        throw DataException("archiveData Error: problem writing to archive file");
1771      }
1772    
1773      //
1774      // Archive underlying data values for each Data type
1775      int noValues;
1776      switch (dataType) {
1777        case 0:
1778          // DataEmpty
1779          noValues = 0;
1780          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1781          cout << "\tnoValues: " << noValues << endl;
1782          break;
1783        case 1:
1784          // DataConstant
1785          noValues = m_data->getLength();
1786          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1787          cout << "\tnoValues: " << noValues << endl;
1788          if (m_data->archiveData(archiveFile,noValues)) {
1789            throw DataException("archiveData Error: problem writing data to archive file");
1790          }
1791          break;
1792        case 2:
1793          // DataTagged
1794          noValues = m_data->getLength();
1795          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1796          cout << "\tnoValues: " << noValues << endl;
1797          if (m_data->archiveData(archiveFile,noValues)) {
1798            throw DataException("archiveData Error: problem writing data to archive file");
1799          }
1800          break;
1801        case 3:
1802          // DataExpanded
1803          noValues = m_data->getLength();
1804          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1805          cout << "\tnoValues: " << noValues << endl;
1806          if (m_data->archiveData(archiveFile,noValues)) {
1807            throw DataException("archiveData Error: problem writing data to archive file");
1808          }
1809          break;
1810      }
1811    
1812      if (!archiveFile.good()) {
1813        throw DataException("archiveData Error: problem writing data to archive file");
1814      }
1815    
1816      //
1817      // Close archive file
1818      archiveFile.close();
1819    
1820      if (!archiveFile.good()) {
1821        throw DataException("archiveData Error: problem closing archive file");
1822      }
1823    
1824    }
1825    
1826    void
1827    Data::extractData(const std::string fileName,
1828                      const FunctionSpace& fspace)
1829    {
1830      //
1831      // Can only extract Data to an object which is initially DataEmpty
1832      if (!isEmpty()) {
1833        throw DataException("extractData Error: can only extract to DataEmpty object");
1834      }
1835    
1836      cout << "Extracting Data object from: " << fileName << endl;
1837    
1838      int dataType;
1839      int noSamples;
1840      int noDPPSample;
1841      int functionSpaceType;
1842      int dataPointRank;
1843      int dataPointSize;
1844      int dataLength;
1845      DataArrayView::ShapeType dataPointShape;
1846      int flatShape[4];
1847    
1848      //
1849      // Open the archive file
1850      ifstream archiveFile;
1851      archiveFile.open(fileName.data(), ios::in);
1852    
1853      if (!archiveFile.good()) {
1854        throw DataException("extractData Error: problem opening archive file");
1855      }
1856    
1857      //
1858      // Read common data items from archive file
1859      archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1860      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1861      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1862      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1863      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1864      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1865      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1866      for (int dim = 0; dim < 4; dim++) {
1867        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1868        if (flatShape[dim]>0) {
1869          dataPointShape.push_back(flatShape[dim]);
1870        }
1871      }
1872      int referenceNumbers[noSamples];
1873      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1874        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1875      }
1876      int tagNumbers[noSamples];
1877      if (dataType==2) {
1878        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1879          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1880        }
1881      }
1882    
1883      if (!archiveFile.good()) {
1884        throw DataException("extractData Error: problem reading from archive file");
1885      }
1886    
1887      //
1888      // Verify the values just read from the archive file
1889      switch (dataType) {
1890        case 0:
1891          cout << "\tdataType: DataEmpty" << endl;
1892          break;
1893        case 1:
1894          cout << "\tdataType: DataConstant" << endl;
1895          break;
1896        case 2:
1897          cout << "\tdataType: DataTagged" << endl;
1898          break;
1899        case 3:
1900          cout << "\tdataType: DataExpanded" << endl;
1901          break;
1902        default:
1903          throw DataException("extractData Error: undefined dataType read from archive file");
1904          break;
1905      }
1906    
1907      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1908      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1909      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1910      cout << "\tshape: < ";
1911      for (int dim = 0; dim < dataPointRank; dim++) {
1912        cout << dataPointShape[dim] << " ";
1913      }
1914      cout << ">" << endl;
1915    
1916      //
1917      // Verify that supplied FunctionSpace object is compatible with this Data object.
1918      if ( (fspace.getTypeCode()!=functionSpaceType) ||
1919           (fspace.getNumSamples()!=noSamples) ||
1920           (fspace.getNumDPPSample()!=noDPPSample)
1921         ) {
1922        throw DataException("extractData Error: incompatible FunctionSpace");
1923      }
1924      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1925        if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
1926          throw DataException("extractData Error: incompatible FunctionSpace");
1927        }
1928      }
1929      if (dataType==2) {
1930        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1931          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
1932            throw DataException("extractData Error: incompatible FunctionSpace");
1933          }
1934        }
1935      }
1936    
1937      //
1938      // Construct a DataVector to hold underlying data values
1939      DataVector dataVec(dataLength);
1940    
1941      //
1942      // Load this DataVector with the appropriate values
1943      int noValues;
1944      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
1945      cout << "\tnoValues: " << noValues << endl;
1946      switch (dataType) {
1947        case 0:
1948          // DataEmpty
1949          if (noValues != 0) {
1950            throw DataException("extractData Error: problem reading data from archive file");
1951          }
1952          break;
1953        case 1:
1954          // DataConstant
1955          if (dataVec.extractData(archiveFile,noValues)) {
1956            throw DataException("extractData Error: problem reading data from archive file");
1957          }
1958          break;
1959        case 2:
1960          // DataTagged
1961          if (dataVec.extractData(archiveFile,noValues)) {
1962            throw DataException("extractData Error: problem reading data from archive file");
1963          }
1964          break;
1965        case 3:
1966          // DataExpanded
1967          if (dataVec.extractData(archiveFile,noValues)) {
1968            throw DataException("extractData Error: problem reading data from archive file");
1969          }
1970          break;
1971      }
1972    
1973      if (!archiveFile.good()) {
1974        throw DataException("extractData Error: problem reading from archive file");
1975      }
1976    
1977      //
1978      // Close archive file
1979      archiveFile.close();
1980    
1981      if (!archiveFile.good()) {
1982        throw DataException("extractData Error: problem closing archive file");
1983      }
1984    
1985      //
1986      // Construct an appropriate Data object
1987      DataAbstract* tempData;
1988      switch (dataType) {
1989        case 0:
1990          // DataEmpty
1991          tempData=new DataEmpty();
1992          break;
1993        case 1:
1994          // DataConstant
1995          tempData=new DataConstant(fspace,dataPointShape,dataVec);
1996          break;
1997        case 2:
1998          // DataTagged
1999          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2000          break;
2001        case 3:
2002          // DataExpanded
2003          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2004          break;
2005      }
2006      shared_ptr<DataAbstract> temp_data(tempData);
2007      m_data=temp_data;
2008  }  }
 */  
2009    
2010  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
2011  {  {

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

  ViewVC Help
Powered by ViewVC 1.1.26