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

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

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

revision 117 by jgs, Fri Apr 1 05:48:57 2005 UTC revision 148 by jgs, Tue Aug 23 01:24:31 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      // check rank of num_array
489      for (int i=0; i<getDataPointRank(); i++) {
490        if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])
491           throw DataException("Shape of numarray does not match Data object rank");
492      }
493      //
494      // make sure data is expanded:
495      if (! isExpanded()) expand();
496      //
497      // and copy over:
498      m_data->copyAll(num_array);
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      //       dataPointView()=numArray[dataPoint];
540      //     }
541      //     if (dataPointRank==1) {
542      //       for (int i=0; i<dataPointShape[0]; i++) {
543      //         dataPointView(i)=numArray[dataPoint][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    
583    const
584  boost::python::numeric::array  boost::python::numeric::array
585  Data::convertToNumArray()  Data::convertToNumArray()
586  {  {
587    //    //
588    // determine the current number of data points    // determine the total number of data points
589    int numSamples = getNumSamples();    int numSamples = getNumSamples();
590    int numDataPointsPerSample = getNumDataPointsPerSample();    int numDataPointsPerSample = getNumDataPointsPerSample();
591    int numDataPoints = numSamples * numDataPointsPerSample;    int numDataPoints = numSamples * numDataPointsPerSample;
# Line 499  Data::convertToNumArray() Line 637  Data::convertToNumArray()
637    int dataPoint = 0;    int dataPoint = 0;
638    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
639      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
   
640        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
   
641        if (dataPointRank==0) {        if (dataPointRank==0) {
642          numArray[dataPoint]=dataPointView();          numArray[dataPoint]=dataPointView();
643        }        }
   
644        if (dataPointRank==1) {        if (dataPointRank==1) {
645          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
646            numArray[dataPoint][i]=dataPointView(i);            numArray[dataPoint][i]=dataPointView(i);
647          }          }
648        }        }
   
649        if (dataPointRank==2) {        if (dataPointRank==2) {
650          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
651            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
# Line 519  Data::convertToNumArray() Line 653  Data::convertToNumArray()
653            }            }
654          }          }
655        }        }
   
656        if (dataPointRank==3) {        if (dataPointRank==3) {
657          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
658            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
# Line 529  Data::convertToNumArray() Line 662  Data::convertToNumArray()
662            }            }
663          }          }
664        }        }
   
665        if (dataPointRank==4) {        if (dataPointRank==4) {
666          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
667            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
# Line 541  Data::convertToNumArray() Line 673  Data::convertToNumArray()
673            }            }
674          }          }
675        }        }
   
676        dataPoint++;        dataPoint++;
677        }
678      }
679    
680      //
681      // return the loaded array
682      return numArray;
683    }
684    
685    const
686    boost::python::numeric::array
687    Data::convertToNumArrayFromSampleNo(int sampleNo)
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      // determine the number of data points per sample
697      int numDataPointsPerSample = getNumDataPointsPerSample();
698    
699      //
700      // determine the rank and shape of each data point
701      int dataPointRank = getDataPointRank();
702      DataArrayView::ShapeType dataPointShape = getDataPointShape();
703    
704      //
705      // create the numeric array to be returned
706      boost::python::numeric::array numArray(0.0);
707    
708      //
709      // the rank of the returned numeric array will be the rank of
710      // the data points, plus one. Where the rank of the array is n,
711      // the last n-1 dimensions will be equal to the shape of the
712      // data points, whilst the first dimension will be equal to the
713      // total number of data points. Thus the array will consist of
714      // a serial vector of the data points.
715      int arrayRank = dataPointRank + 1;
716      DataArrayView::ShapeType arrayShape;
717      arrayShape.push_back(numDataPointsPerSample);
718      for (int d=0; d<dataPointRank; d++) {
719         arrayShape.push_back(dataPointShape[d]);
720      }
721    
722      //
723      // resize the numeric array to the shape just calculated
724      if (arrayRank==1) {
725        numArray.resize(arrayShape[0]);
726      }
727      if (arrayRank==2) {
728        numArray.resize(arrayShape[0],arrayShape[1]);
729      }
730      if (arrayRank==3) {
731        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
732      }
733      if (arrayRank==4) {
734        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
735      }
736      if (arrayRank==5) {
737        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
738      }
739    
740      //
741      // loop through each data point in turn, loading the values for that data point
742      // into the numeric array.
743      for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {
744        DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);
745        if (dataPointRank==0) {
746          numArray[dataPoint]=dataPointView();
747        }
748        if (dataPointRank==1) {
749          for (int i=0; i<dataPointShape[0]; i++) {
750            numArray[dataPoint][i]=dataPointView(i);
751          }
752        }
753        if (dataPointRank==2) {
754          for (int i=0; i<dataPointShape[0]; i++) {
755            for (int j=0; j<dataPointShape[1]; j++) {
756              numArray[dataPoint][i][j] = dataPointView(i,j);
757            }
758          }
759        }
760        if (dataPointRank==3) {
761          for (int i=0; i<dataPointShape[0]; i++) {
762            for (int j=0; j<dataPointShape[1]; j++) {
763              for (int k=0; k<dataPointShape[2]; k++) {
764                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
765              }
766            }
767          }
768        }
769        if (dataPointRank==4) {
770          for (int i=0; i<dataPointShape[0]; i++) {
771            for (int j=0; j<dataPointShape[1]; j++) {
772              for (int k=0; k<dataPointShape[2]; k++) {
773                for (int l=0; l<dataPointShape[3]; l++) {
774                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
775                }
776              }
777            }
778          }
779        }
780      }
781    
782      //
783      // return the loaded array
784      return numArray;
785    }
786    
787    const
788    boost::python::numeric::array
789    Data::convertToNumArrayFromDPNo(int sampleNo,
790                                    int dataPointNo)
791    {
792      //
793      // Check a valid sample number has been supplied
794      if (sampleNo >= getNumSamples()) {
795        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
796      }
797    
798      //
799      // Check a valid data point number has been supplied
800      if (dataPointNo >= getNumDataPointsPerSample()) {
801        throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");
802      }
803    
804      //
805      // determine the rank and shape of each data point
806      int dataPointRank = getDataPointRank();
807      DataArrayView::ShapeType dataPointShape = getDataPointShape();
808    
809      //
810      // create the numeric array to be returned
811      boost::python::numeric::array numArray(0.0);
812    
813      //
814      // the shape of the returned numeric array will be the same
815      // as that of the data point
816      int arrayRank = dataPointRank;
817      DataArrayView::ShapeType arrayShape = dataPointShape;
818    
819      //
820      // resize the numeric array to the shape just calculated
821      if (arrayRank==0) {
822        numArray.resize(1);
823      }
824      if (arrayRank==1) {
825        numArray.resize(arrayShape[0]);
826      }
827      if (arrayRank==2) {
828        numArray.resize(arrayShape[0],arrayShape[1]);
829      }
830      if (arrayRank==3) {
831        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
832      }
833      if (arrayRank==4) {
834        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
835      }
836    
837      //
838      // load the values for the data point into the numeric array.
839      DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
840      if (dataPointRank==0) {
841        numArray[0]=dataPointView();
842      }
843      if (dataPointRank==1) {
844        for (int i=0; i<dataPointShape[0]; i++) {
845          numArray[i]=dataPointView(i);
846        }
847      }
848      if (dataPointRank==2) {
849        for (int i=0; i<dataPointShape[0]; i++) {
850          for (int j=0; j<dataPointShape[1]; j++) {
851            numArray[i][j] = dataPointView(i,j);
852          }
853        }
854      }
855      if (dataPointRank==3) {
856        for (int i=0; i<dataPointShape[0]; i++) {
857          for (int j=0; j<dataPointShape[1]; j++) {
858            for (int k=0; k<dataPointShape[2]; k++) {
859              numArray[i][j][k]=dataPointView(i,j,k);
860            }
861          }
862        }
863      }
864      if (dataPointRank==4) {
865        for (int i=0; i<dataPointShape[0]; i++) {
866          for (int j=0; j<dataPointShape[1]; j++) {
867            for (int k=0; k<dataPointShape[2]; k++) {
868              for (int l=0; l<dataPointShape[3]; l++) {
869                numArray[i][j][k][l]=dataPointView(i,j,k,l);
870              }
871            }
872          }
873      }      }
874    }    }
875    
# Line 559  Data::integrate() const Line 885  Data::integrate() const
885    int rank = getDataPointRank();    int rank = getDataPointRank();
886    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
887    
888      profData->integrate++;
889    
890    //    //
891    // calculate the integral values    // calculate the integral values
892    vector<double> integrals(getDataPointSize());    vector<double> integrals(getDataPointSize());
# Line 622  Data::integrate() const Line 950  Data::integrate() const
950  Data  Data
951  Data::sin() const  Data::sin() const
952  {  {
953      profData->unary++;
954    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
955  }  }
956    
957  Data  Data
958  Data::cos() const  Data::cos() const
959  {  {
960      profData->unary++;
961    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
962  }  }
963    
964  Data  Data
965  Data::tan() const  Data::tan() const
966  {  {
967      profData->unary++;
968    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
969  }  }
970    
971  Data  Data
972  Data::log() const  Data::log() const
973  {  {
974      profData->unary++;
975    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
976  }  }
977    
978  Data  Data
979  Data::ln() const  Data::ln() const
980  {  {
981      profData->unary++;
982    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
983  }  }
984    
985  Data  Data
986  Data::sign() const  Data::sign() const
987  {  {
988      profData->unary++;
989    return escript::unaryOp(*this,escript::fsign);    return escript::unaryOp(*this,escript::fsign);
990  }  }
991    
992  Data  Data
993  Data::abs() const  Data::abs() const
994  {  {
995      profData->unary++;
996    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
997  }  }
998    
999  Data  Data
1000  Data::neg() const  Data::neg() const
1001  {  {
1002      profData->unary++;
1003    return escript::unaryOp(*this,negate<double>());    return escript::unaryOp(*this,negate<double>());
1004  }  }
1005    
1006  Data  Data
1007  Data::pos() const  Data::pos() const
1008  {  {
1009    return (*this);    profData->unary++;
1010      Data result;
1011      // perform a deep copy
1012      result.copy(*this);
1013      return result;
1014  }  }
1015    
1016  Data  Data
1017  Data::exp() const  Data::exp() const
1018  {  {
1019      profData->unary++;
1020    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1021  }  }
1022    
1023  Data  Data
1024  Data::sqrt() const  Data::sqrt() const
1025  {  {
1026      profData->unary++;
1027    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1028  }  }
1029    
1030  double  double
1031  Data::Lsup() const  Data::Lsup() const
1032  {  {
1033      profData->reduction1++;
1034    //    //
1035    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
1036    return algorithm(DataAlgorithmAdapter<AbsMax>(0));    AbsMax abs_max_func;
1037      return algorithm(abs_max_func,0);
1038  }  }
1039    
1040  double  double
1041  Data::Linf() const  Data::Linf() const
1042  {  {
1043      profData->reduction1++;
1044    //    //
1045    // set the initial absolute minimum value to max double    // set the initial absolute minimum value to max double
1046    return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));    AbsMin abs_min_func;
1047      return algorithm(abs_min_func,numeric_limits<double>::max());
1048  }  }
1049    
1050  double  double
1051  Data::sup() const  Data::sup() const
1052  {  {
1053      profData->reduction1++;
1054    //    //
1055    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1056    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1057      return algorithm(fmax_func,numeric_limits<double>::max()*-1);
1058  }  }
1059    
1060  double  double
1061  Data::inf() const  Data::inf() const
1062  {  {
1063      profData->reduction1++;
1064    //    //
1065    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1066    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1067      return algorithm(fmin_func,numeric_limits<double>::max());
1068  }  }
1069    
1070  Data  Data
1071  Data::maxval() const  Data::maxval() const
1072  {  {
1073      profData->reduction2++;
1074    //    //
1075    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1076    return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1077      return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1078  }  }
1079    
1080  Data  Data
1081  Data::minval() const  Data::minval() const
1082  {  {
1083      profData->reduction2++;
1084    //    //
1085    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1086    return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1087      return dp_algorithm(fmin_func,numeric_limits<double>::max());
1088  }  }
1089    
1090  Data  Data
1091  Data::length() const  Data::length() const
1092  {  {
1093    return dp_algorithm(DataAlgorithmAdapter<Length>(0));    profData->reduction2++;
1094      Length len_func;
1095      return dp_algorithm(len_func,0);
1096  }  }
1097    
1098  Data  Data
1099  Data::trace() const  Data::trace() const
1100  {  {
1101    return dp_algorithm(DataAlgorithmAdapter<Trace>(0));    profData->reduction2++;
1102      Trace trace_func;
1103      return dp_algorithm(trace_func,0);
1104  }  }
1105    
1106  Data  Data
1107  Data::transpose(int axis) const  Data::transpose(int axis) const
1108  {  {
1109      profData->reduction2++;
1110    // not implemented    // not implemented
1111    throw DataException("Error - Data::transpose not implemented yet.");    throw DataException("Error - Data::transpose not implemented yet.");
1112    return Data();    return Data();
1113  }  }
1114    
1115    const boost::python::tuple
1116    Data::mindp() const
1117    {
1118      // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1119      // abort (for unknown reasons) if there are openmp directives with it in the
1120      // surrounding function
1121    
1122      int SampleNo;
1123      int DataPointNo;
1124    
1125      calc_mindp(SampleNo,DataPointNo);
1126    
1127      return make_tuple(SampleNo,DataPointNo);
1128    }
1129    
1130    void
1131    Data::calc_mindp(int& SampleNo,
1132                     int& DataPointNo) const
1133    {
1134      int i,j;
1135      int lowi=0,lowj=0;
1136      double min=numeric_limits<double>::max();
1137    
1138      Data temp=minval();
1139    
1140      int numSamples=temp.getNumSamples();
1141      int numDPPSample=temp.getNumDataPointsPerSample();
1142    
1143      double next,local_min;
1144      int local_lowi,local_lowj;
1145    
1146      #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1147      {
1148        local_min=min;
1149        #pragma omp for private(i,j) schedule(static)
1150        for (i=0; i<numSamples; i++) {
1151          for (j=0; j<numDPPSample; j++) {
1152            next=temp.getDataPoint(i,j)();
1153            if (next<local_min) {
1154              local_min=next;
1155              local_lowi=i;
1156              local_lowj=j;
1157            }
1158          }
1159        }
1160        #pragma omp critical
1161        if (local_min<min) {
1162          min=local_min;
1163          lowi=local_lowi;
1164          lowj=local_lowj;
1165        }
1166      }
1167    
1168      SampleNo = lowi;
1169      DataPointNo = lowj;
1170    }
1171    
1172  void  void
1173  Data::saveDX(std::string fileName) const  Data::saveDX(std::string fileName) const
1174  {  {
# Line 770  Data::saveVTK(std::string fileName) cons Line 1186  Data::saveVTK(std::string fileName) cons
1186  Data&  Data&
1187  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1188  {  {
1189      profData->binary++;
1190    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1191    return (*this);    return (*this);
1192  }  }
# Line 777  Data::operator+=(const Data& right) Line 1194  Data::operator+=(const Data& right)
1194  Data&  Data&
1195  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1196  {  {
1197      profData->binary++;
1198    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1199    return (*this);    return (*this);
1200  }  }
# Line 784  Data::operator+=(const boost::python::ob Line 1202  Data::operator+=(const boost::python::ob
1202  Data&  Data&
1203  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1204  {  {
1205      profData->binary++;
1206    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1207    return (*this);    return (*this);
1208  }  }
# Line 791  Data::operator-=(const Data& right) Line 1210  Data::operator-=(const Data& right)
1210  Data&  Data&
1211  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1212  {  {
1213      profData->binary++;
1214    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1215    return (*this);    return (*this);
1216  }  }
# Line 798  Data::operator-=(const boost::python::ob Line 1218  Data::operator-=(const boost::python::ob
1218  Data&  Data&
1219  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1220  {  {
1221      profData->binary++;
1222    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1223    return (*this);    return (*this);
1224  }  }
# Line 805  Data::operator*=(const Data& right) Line 1226  Data::operator*=(const Data& right)
1226  Data&  Data&
1227  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1228  {  {
1229      profData->binary++;
1230    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1231    return (*this);    return (*this);
1232  }  }
# Line 812  Data::operator*=(const boost::python::ob Line 1234  Data::operator*=(const boost::python::ob
1234  Data&  Data&
1235  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1236  {  {
1237      profData->binary++;
1238    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1239    return (*this);    return (*this);
1240  }  }
# Line 819  Data::operator/=(const Data& right) Line 1242  Data::operator/=(const Data& right)
1242  Data&  Data&
1243  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1244  {  {
1245      profData->binary++;
1246    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1247    return (*this);    return (*this);
1248  }  }
# Line 826  Data::operator/=(const boost::python::ob Line 1250  Data::operator/=(const boost::python::ob
1250  Data  Data
1251  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1252  {  {
1253      profData->binary++;
1254    Data result;    Data result;
1255    result.copy(*this);    result.copy(*this);
1256    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 835  Data::powO(const boost::python::object& Line 1260  Data::powO(const boost::python::object&
1260  Data  Data
1261  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1262  {  {
1263      profData->binary++;
1264    Data result;    Data result;
1265    result.copy(*this);    result.copy(*this);
1266    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 842  Data::powD(const Data& right) const Line 1268  Data::powD(const Data& right) const
1268  }  }
1269    
1270  //  //
1271  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1272  Data  Data
1273  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1274  {  {
# Line 855  escript::operator+(const Data& left, con Line 1281  escript::operator+(const Data& left, con
1281  }  }
1282    
1283  //  //
1284  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1285  Data  Data
1286  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1287  {  {
# Line 868  escript::operator-(const Data& left, con Line 1294  escript::operator-(const Data& left, con
1294  }  }
1295    
1296  //  //
1297  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1298  Data  Data
1299  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
1300  {  {
# Line 881  escript::operator*(const Data& left, con Line 1307  escript::operator*(const Data& left, con
1307  }  }
1308    
1309  //  //
1310  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1311  Data  Data
1312  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
1313  {  {
# Line 894  escript::operator/(const Data& left, con Line 1320  escript::operator/(const Data& left, con
1320  }  }
1321    
1322  //  //
1323  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1324  Data  Data
1325  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1326  {  {
# Line 910  escript::operator+(const Data& left, con Line 1336  escript::operator+(const Data& left, con
1336  }  }
1337    
1338  //  //
1339  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1340  Data  Data
1341  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1342  {  {
# Line 926  escript::operator-(const Data& left, con Line 1352  escript::operator-(const Data& left, con
1352  }  }
1353    
1354  //  //
1355  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1356  Data  Data
1357  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1358  {  {
# Line 942  escript::operator*(const Data& left, con Line 1368  escript::operator*(const Data& left, con
1368  }  }
1369    
1370  //  //
1371  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1372  Data  Data
1373  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1374  {  {
# Line 958  escript::operator/(const Data& left, con Line 1384  escript::operator/(const Data& left, con
1384  }  }
1385    
1386  //  //
1387  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1388  Data  Data
1389  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1390  {  {
# Line 971  escript::operator+(const boost::python:: Line 1397  escript::operator+(const boost::python::
1397  }  }
1398    
1399  //  //
1400  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1401  Data  Data
1402  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1403  {  {
# Line 984  escript::operator-(const boost::python:: Line 1410  escript::operator-(const boost::python::
1410  }  }
1411    
1412  //  //
1413  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1414  Data  Data
1415  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1416  {  {
# Line 997  escript::operator*(const boost::python:: Line 1423  escript::operator*(const boost::python::
1423  }  }
1424    
1425  //  //
1426  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1427  Data  Data
1428  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1429  {  {
# Line 1010  escript::operator/(const boost::python:: Line 1436  escript::operator/(const boost::python::
1436  }  }
1437    
1438  //  //
 // NOTE: It is essential to specify the namepsace this operator belongs to  
1439  //bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1440  //{  //{
1441  //  /*  //  /*
# Line 1072  Data::getItem(const boost::python::objec Line 1497  Data::getItem(const boost::python::objec
1497  Data  Data
1498  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1499  {  {
1500      profData->slicing++;
1501    return Data(*this,region);    return Data(*this,region);
1502  }  }
1503    
# Line 1088  Data::setItemD(const boost::python::obje Line 1514  Data::setItemD(const boost::python::obje
1514                 const Data& value)                 const Data& value)
1515  {  {
1516    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1517    
1518    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1519    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1520      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
# Line 1103  void Line 1530  void
1530  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1531                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
1532  {  {
1533      profData->slicing++;
1534    Data tempValue(value);    Data tempValue(value);
1535    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1536    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1159  Data::setTaggedValue(int tagKey, Line 1587  Data::setTaggedValue(int tagKey,
1587  }  }
1588    
1589  void  void
1590    Data::setTaggedValueFromCPP(int tagKey,
1591                                const DataArrayView& value)
1592    {
1593      //
1594      // Ensure underlying data object is of type DataTagged
1595      tag();
1596    
1597      if (!isTagged()) {
1598        throw DataException("Error - DataTagged conversion failed!!");
1599      }
1600                                                                                                                  
1601      //
1602      // Call DataAbstract::setTaggedValue
1603      m_data->setTaggedValue(tagKey,value);
1604    }
1605    
1606    void
1607  Data::setRefValue(int ref,  Data::setRefValue(int ref,
1608                    const boost::python::numeric::array& value)                    const boost::python::numeric::array& value)
1609  {  {
# Line 1197  Data::getRefValue(int ref, Line 1642  Data::getRefValue(int ref,
1642    DataArrayView valueView = valueDataArray.getView();    DataArrayView valueView = valueDataArray.getView();
1643    
1644    if (rank==0) {    if (rank==0) {
1645      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");        boost::python::numeric::array temp_numArray(valueView());
1646          value = temp_numArray;
1647    }    }
1648    if (rank==1) {    if (rank==1) {
1649      for (int i=0; i < shape[0]; i++) {      for (int i=0; i < shape[0]; i++) {
# Line 1205  Data::getRefValue(int ref, Line 1651  Data::getRefValue(int ref,
1651      }      }
1652    }    }
1653    if (rank==2) {    if (rank==2) {
1654      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
1655          for (int j=0; j < shape[1]; j++) {
1656            value[i][j] = valueView(i,j);
1657          }
1658        }
1659    }    }
1660    if (rank==3) {    if (rank==3) {
1661      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
1662          for (int j=0; j < shape[1]; j++) {
1663            for (int k=0; k < shape[2]; k++) {
1664              value[i][j][k] = valueView(i,j,k);
1665            }
1666          }
1667        }
1668    }    }
1669    if (rank==4) {    if (rank==4) {
1670      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
1671          for (int j=0; j < shape[1]; j++) {
1672            for (int k=0; k < shape[2]; k++) {
1673              for (int l=0; l < shape[3]; l++) {
1674                value[i][j][k][l] = valueView(i,j,k,l);
1675              }
1676            }
1677          }
1678        }
1679    }    }
1680    
1681  }  }
1682    
 /*  
 Note: this version removed for now. Not needed, and breaks escript.cpp  
1683  void  void
1684  Data::setTaggedValue(int tagKey,  Data::archiveData(const std::string fileName)
                      const DataArrayView& value)  
1685  {  {
1686      cout << "Archiving Data object to: " << fileName << endl;
1687    
1688    //    //
1689    // Ensure underlying data object is of type DataTagged    // Determine type of this Data object
1690    tag();    int dataType = -1;
1691    
1692    if (!isTagged()) {    if (isEmpty()) {
1693      throw DataException("Error - DataTagged conversion failed!!");      dataType = 0;
1694        cout << "\tdataType: DataEmpty" << endl;
1695    }    }
1696                                                                                                                    if (isConstant()) {
1697        dataType = 1;
1698        cout << "\tdataType: DataConstant" << endl;
1699      }
1700      if (isTagged()) {
1701        dataType = 2;
1702        cout << "\tdataType: DataTagged" << endl;
1703      }
1704      if (isExpanded()) {
1705        dataType = 3;
1706        cout << "\tdataType: DataExpanded" << endl;
1707      }
1708    
1709      if (dataType == -1) {
1710        throw DataException("archiveData Error: undefined dataType");
1711      }
1712    
1713    //    //
1714    // Call DataAbstract::setTaggedValue    // Collect data items common to all Data types
1715    m_data->setTaggedValue(tagKey,value);    int noSamples = getNumSamples();
1716      int noDPPSample = getNumDataPointsPerSample();
1717      int functionSpaceType = getFunctionSpace().getTypeCode();
1718      int dataPointRank = getDataPointRank();
1719      int dataPointSize = getDataPointSize();
1720      int dataLength = getLength();
1721      DataArrayView::ShapeType dataPointShape = getDataPointShape();
1722      int referenceNumbers[noSamples];
1723      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1724        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1725      }
1726      int tagNumbers[noSamples];
1727      if (isTagged()) {
1728        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1729          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1730        }
1731      }
1732    
1733      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1734      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1735      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1736    
1737      //
1738      // Flatten Shape to an array of integers suitable for writing to file
1739      int flatShape[4] = {0,0,0,0};
1740      cout << "\tshape: < ";
1741      for (int dim=0; dim<dataPointRank; dim++) {
1742        flatShape[dim] = dataPointShape[dim];
1743        cout << dataPointShape[dim] << " ";
1744      }
1745      cout << ">" << endl;
1746    
1747      //
1748      // Open archive file
1749      ofstream archiveFile;
1750      archiveFile.open(fileName.data(), ios::out);
1751    
1752      if (!archiveFile.good()) {
1753        throw DataException("archiveData Error: problem opening archive file");
1754      }
1755    
1756      //
1757      // Write common data items to archive file
1758      archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1759      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1760      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1761      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1762      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1763      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1764      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1765      for (int dim = 0; dim < 4; dim++) {
1766        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1767      }
1768      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1769        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1770      }
1771      if (isTagged()) {
1772        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1773          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1774        }
1775      }
1776    
1777      if (!archiveFile.good()) {
1778        throw DataException("archiveData Error: problem writing to archive file");
1779      }
1780    
1781      //
1782      // Archive underlying data values for each Data type
1783      int noValues;
1784      switch (dataType) {
1785        case 0:
1786          // DataEmpty
1787          noValues = 0;
1788          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1789          cout << "\tnoValues: " << noValues << endl;
1790          break;
1791        case 1:
1792          // DataConstant
1793          noValues = m_data->getLength();
1794          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1795          cout << "\tnoValues: " << noValues << endl;
1796          if (m_data->archiveData(archiveFile,noValues)) {
1797            throw DataException("archiveData Error: problem writing data to archive file");
1798          }
1799          break;
1800        case 2:
1801          // DataTagged
1802          noValues = m_data->getLength();
1803          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1804          cout << "\tnoValues: " << noValues << endl;
1805          if (m_data->archiveData(archiveFile,noValues)) {
1806            throw DataException("archiveData Error: problem writing data to archive file");
1807          }
1808          break;
1809        case 3:
1810          // DataExpanded
1811          noValues = m_data->getLength();
1812          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1813          cout << "\tnoValues: " << noValues << endl;
1814          if (m_data->archiveData(archiveFile,noValues)) {
1815            throw DataException("archiveData Error: problem writing data to archive file");
1816          }
1817          break;
1818      }
1819    
1820      if (!archiveFile.good()) {
1821        throw DataException("archiveData Error: problem writing data to archive file");
1822      }
1823    
1824      //
1825      // Close archive file
1826      archiveFile.close();
1827    
1828      if (!archiveFile.good()) {
1829        throw DataException("archiveData Error: problem closing archive file");
1830      }
1831    
1832    }
1833    
1834    void
1835    Data::extractData(const std::string fileName,
1836                      const FunctionSpace& fspace)
1837    {
1838      //
1839      // Can only extract Data to an object which is initially DataEmpty
1840      if (!isEmpty()) {
1841        throw DataException("extractData Error: can only extract to DataEmpty object");
1842      }
1843    
1844      cout << "Extracting Data object from: " << fileName << endl;
1845    
1846      int dataType;
1847      int noSamples;
1848      int noDPPSample;
1849      int functionSpaceType;
1850      int dataPointRank;
1851      int dataPointSize;
1852      int dataLength;
1853      DataArrayView::ShapeType dataPointShape;
1854      int flatShape[4];
1855    
1856      //
1857      // Open the archive file
1858      ifstream archiveFile;
1859      archiveFile.open(fileName.data(), ios::in);
1860    
1861      if (!archiveFile.good()) {
1862        throw DataException("extractData Error: problem opening archive file");
1863      }
1864    
1865      //
1866      // Read common data items from archive file
1867      archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1868      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1869      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1870      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1871      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1872      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1873      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1874      for (int dim = 0; dim < 4; dim++) {
1875        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1876        if (flatShape[dim]>0) {
1877          dataPointShape.push_back(flatShape[dim]);
1878        }
1879      }
1880      int referenceNumbers[noSamples];
1881      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1882        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1883      }
1884      int tagNumbers[noSamples];
1885      if (dataType==2) {
1886        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1887          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1888        }
1889      }
1890    
1891      if (!archiveFile.good()) {
1892        throw DataException("extractData Error: problem reading from archive file");
1893      }
1894    
1895      //
1896      // Verify the values just read from the archive file
1897      switch (dataType) {
1898        case 0:
1899          cout << "\tdataType: DataEmpty" << endl;
1900          break;
1901        case 1:
1902          cout << "\tdataType: DataConstant" << endl;
1903          break;
1904        case 2:
1905          cout << "\tdataType: DataTagged" << endl;
1906          break;
1907        case 3:
1908          cout << "\tdataType: DataExpanded" << endl;
1909          break;
1910        default:
1911          throw DataException("extractData Error: undefined dataType read from archive file");
1912          break;
1913      }
1914    
1915      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1916      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1917      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1918      cout << "\tshape: < ";
1919      for (int dim = 0; dim < dataPointRank; dim++) {
1920        cout << dataPointShape[dim] << " ";
1921      }
1922      cout << ">" << endl;
1923    
1924      //
1925      // Verify that supplied FunctionSpace object is compatible with this Data object.
1926      if ( (fspace.getTypeCode()!=functionSpaceType) ||
1927           (fspace.getNumSamples()!=noSamples) ||
1928           (fspace.getNumDPPSample()!=noDPPSample)
1929         ) {
1930        throw DataException("extractData Error: incompatible FunctionSpace");
1931      }
1932      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1933        if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
1934          throw DataException("extractData Error: incompatible FunctionSpace");
1935        }
1936      }
1937      if (dataType==2) {
1938        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1939          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
1940            throw DataException("extractData Error: incompatible FunctionSpace");
1941          }
1942        }
1943      }
1944    
1945      //
1946      // Construct a DataVector to hold underlying data values
1947      DataVector dataVec(dataLength);
1948    
1949      //
1950      // Load this DataVector with the appropriate values
1951      int noValues;
1952      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
1953      cout << "\tnoValues: " << noValues << endl;
1954      switch (dataType) {
1955        case 0:
1956          // DataEmpty
1957          if (noValues != 0) {
1958            throw DataException("extractData Error: problem reading data from archive file");
1959          }
1960          break;
1961        case 1:
1962          // DataConstant
1963          if (dataVec.extractData(archiveFile,noValues)) {
1964            throw DataException("extractData Error: problem reading data from archive file");
1965          }
1966          break;
1967        case 2:
1968          // DataTagged
1969          if (dataVec.extractData(archiveFile,noValues)) {
1970            throw DataException("extractData Error: problem reading data from archive file");
1971          }
1972          break;
1973        case 3:
1974          // DataExpanded
1975          if (dataVec.extractData(archiveFile,noValues)) {
1976            throw DataException("extractData Error: problem reading data from archive file");
1977          }
1978          break;
1979      }
1980    
1981      if (!archiveFile.good()) {
1982        throw DataException("extractData Error: problem reading from archive file");
1983      }
1984    
1985      //
1986      // Close archive file
1987      archiveFile.close();
1988    
1989      if (!archiveFile.good()) {
1990        throw DataException("extractData Error: problem closing archive file");
1991      }
1992    
1993      //
1994      // Construct an appropriate Data object
1995      DataAbstract* tempData;
1996      switch (dataType) {
1997        case 0:
1998          // DataEmpty
1999          tempData=new DataEmpty();
2000          break;
2001        case 1:
2002          // DataConstant
2003          tempData=new DataConstant(fspace,dataPointShape,dataVec);
2004          break;
2005        case 2:
2006          // DataTagged
2007          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2008          break;
2009        case 3:
2010          // DataExpanded
2011          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2012          break;
2013      }
2014      shared_ptr<DataAbstract> temp_data(tempData);
2015      m_data=temp_data;
2016  }  }
 */  
2017    
2018  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
2019  {  {

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

  ViewVC Help
Powered by ViewVC 1.1.26