/[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 94 by jgs, Wed Oct 27 00:45:54 2004 UTC revision 119 by jgs, Tue Apr 12 04:45:05 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"
# Line 38  Line 38 
38  #include "escript/Data/DataAlgorithm.h"  #include "escript/Data/DataAlgorithm.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"
42    
43  using namespace std;  using namespace std;
44  using namespace boost::python;  using namespace boost::python;
# Line 49  Data::Data() Line 50  Data::Data()
50    //    //
51    // Default data is type DataEmpty    // Default data is type DataEmpty
52    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
53    m_data=shared_ptr<DataAbstract>(temp);    shared_ptr<DataAbstract> temp_data(temp);
54      m_data=temp_data;
55  }  }
56    
57  Data::Data(double value,  Data::Data(double value,
# Line 75  Data::Data(double value, Line 77  Data::Data(double value,
77    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
78  }  }
79    
80  Data::Data(const Data& inData):  Data::Data(const Data& inData)
   m_data(inData.m_data)  
81  {  {
82      m_data=inData.m_data;
83  }  }
84    
85  Data::Data(const Data& inData,  Data::Data(const Data& inData,
86             const DataArrayView::RegionType& region)             const DataArrayView::RegionType& region)
87  {  {
88    //    //
89    // Create data which is a subset(slice) of another Data    // Create Data which is a slice of another Data
90    DataAbstract* tmp=inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
91    m_data=shared_ptr<DataAbstract>(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
92      m_data=temp_data;
93  }  }
94    
95  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 117  Data::Data(const DataTagged::TagListType Line 120  Data::Data(const DataTagged::TagListType
120             bool expanded)             bool expanded)
121  {  {
122    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
123    m_data=shared_ptr<DataAbstract>(temp);    shared_ptr<DataAbstract> temp_data(temp);
124      m_data=temp_data;
125    if (expanded) {    if (expanded) {
126      expand();      expand();
127    }    }
# Line 212  Data::copy(const Data& other) Line 216  Data::copy(const Data& other)
216        //        //
217        // Construct a DataExpanded copy        // Construct a DataExpanded copy
218        DataAbstract* newData=new DataExpanded(*temp);        DataAbstract* newData=new DataExpanded(*temp);
219        m_data=shared_ptr<DataAbstract>(newData);        shared_ptr<DataAbstract> temp_data(newData);
220          m_data=temp_data;
221        return;        return;
222      }      }
223    }    }
# Line 220  Data::copy(const Data& other) Line 225  Data::copy(const Data& other)
225      DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());      DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());
226      if (temp!=0) {      if (temp!=0) {
227        //        //
228        // Construct a DataTaggeded copy        // Construct a DataTagged copy
229        DataAbstract* newData=new DataTagged(*temp);        DataAbstract* newData=new DataTagged(*temp);
230        m_data=shared_ptr<DataAbstract>(newData);        shared_ptr<DataAbstract> temp_data(newData);
231          m_data=temp_data;
232        return;        return;
233      }      }
234    }    }
# Line 232  Data::copy(const Data& other) Line 238  Data::copy(const Data& other)
238        //        //
239        // Construct a DataConstant copy        // Construct a DataConstant copy
240        DataAbstract* newData=new DataConstant(*temp);        DataAbstract* newData=new DataConstant(*temp);
241        m_data=shared_ptr<DataAbstract>(newData);        shared_ptr<DataAbstract> temp_data(newData);
242          m_data=temp_data;
243          return;
244        }
245      }
246      {
247        DataEmpty* temp=dynamic_cast<DataEmpty*>(other.m_data.get());
248        if (temp!=0) {
249          //
250          // Construct a DataEmpty copy
251          DataAbstract* newData=new DataEmpty();
252          shared_ptr<DataAbstract> temp_data(newData);
253          m_data=temp_data;
254        return;        return;
255      }      }
256    }    }
# Line 284  Data::isConstant() const Line 302  Data::isConstant() const
302    return (temp!=0);    return (temp!=0);
303  }  }
304    
 Data  
 Data::getSlice(const DataArrayView::RegionType& region) const  
 {  
   return Data(*this,region);  
 }  
   
 void  
 Data::setSlice(const Data& value,  
                const DataArrayView::RegionType& region)  
 {  
   m_data->setSlice(value.m_data.get(), region);  
 }  
   
305  void  void
306  Data::expand()  Data::expand()
307  {  {
308    if (isConstant()) {    if (isConstant()) {
309      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
310      DataAbstract* temp=new DataExpanded(*tempDataConst);      DataAbstract* temp=new DataExpanded(*tempDataConst);
311      m_data=shared_ptr<DataAbstract>(temp);      shared_ptr<DataAbstract> temp_data(temp);
312        m_data=temp_data;
313    } else if (isTagged()) {    } else if (isTagged()) {
314      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
315      DataAbstract* temp=new DataExpanded(*tempDataTag);      DataAbstract* temp=new DataExpanded(*tempDataTag);
316      m_data=shared_ptr<DataAbstract>(temp);      shared_ptr<DataAbstract> temp_data(temp);
317        m_data=temp_data;
318    } else if (isExpanded()) {    } else if (isExpanded()) {
319      //      //
320      // do nothing      // do nothing
# Line 319  Data::expand() Line 326  Data::expand()
326  }  }
327    
328  void  void
 Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  
 {  
   m_data->reshapeDataPoint(shape);  
 }  
   
 void  
329  Data::tag()  Data::tag()
330  {  {
331    if (isConstant()) {    if (isConstant()) {
332      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
333      DataAbstract* temp=new DataTagged(*tempDataConst);      DataAbstract* temp=new DataTagged(*tempDataConst);
334      m_data=shared_ptr<DataAbstract>(temp);      shared_ptr<DataAbstract> temp_data(temp);
335        m_data=temp_data;
336    } else if (isTagged()) {    } else if (isTagged()) {
337      // do nothing      // do nothing
338    } else if (isExpanded()) {    } else if (isExpanded()) {
# Line 342  Data::tag() Line 344  Data::tag()
344    }    }
345  }  }
346    
347    void
348    Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)
349    {
350      m_data->reshapeDataPoint(shape);
351    }
352    
353  Data  Data
354  Data::wherePositive() const  Data::wherePositive() const
355  {  {
# Line 349  Data::wherePositive() const Line 357  Data::wherePositive() const
357  }  }
358    
359  Data  Data
360    Data::whereNegative() const
361    {
362      return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
363    }
364    
365    Data
366  Data::whereNonNegative() const  Data::whereNonNegative() const
367  {  {
368    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
369  }  }
370    
371  Data  Data
372  Data::whereNegative() const  Data::whereNonPositive() const
373  {  {
374    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
375  }  }
376    
377  Data  Data
# Line 367  Data::whereZero() const Line 381  Data::whereZero() const
381  }  }
382    
383  Data  Data
384    Data::whereNonZero() const
385    {
386      return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));
387    }
388    
389    Data
390  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
391  {  {
392    return Data(*this,functionspace);    return Data(*this,functionspace);
# Line 424  Data::getDataPointShape() const Line 444  Data::getDataPointShape() const
444  }  }
445    
446  boost::python::numeric::array  boost::python::numeric::array
447    Data::convertToNumArray()
448    {
449      //
450      // determine the current number of data points
451      int numSamples = getNumSamples();
452      int numDataPointsPerSample = getNumDataPointsPerSample();
453      int numDataPoints = numSamples * numDataPointsPerSample;
454    
455      //
456      // determine the rank and shape of each data point
457      int dataPointRank = getDataPointRank();
458      DataArrayView::ShapeType dataPointShape = getDataPointShape();
459    
460      //
461      // create the numeric array to be returned
462      boost::python::numeric::array numArray(0.0);
463    
464      //
465      // the rank of the returned numeric array will be the rank of
466      // the data points, plus one. Where the rank of the array is n,
467      // the last n-1 dimensions will be equal to the shape of the
468      // data points, whilst the first dimension will be equal to the
469      // total number of data points. Thus the array will consist of
470      // a serial vector of the data points.
471      int arrayRank = dataPointRank + 1;
472      DataArrayView::ShapeType arrayShape;
473      arrayShape.push_back(numDataPoints);
474      for (int d=0; d<dataPointRank; d++) {
475         arrayShape.push_back(dataPointShape[d]);
476      }
477    
478      //
479      // resize the numeric array to the shape just calculated
480      if (arrayRank==1) {
481        numArray.resize(arrayShape[0]);
482      }
483      if (arrayRank==2) {
484        numArray.resize(arrayShape[0],arrayShape[1]);
485      }
486      if (arrayRank==3) {
487        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
488      }
489      if (arrayRank==4) {
490        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
491      }
492      if (arrayRank==5) {
493        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
494      }
495    
496      //
497      // loop through each data point in turn, loading the values for that data point
498      // into the numeric array.
499      int dataPoint = 0;
500      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
501        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
502    
503          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
504    
505          if (dataPointRank==0) {
506            numArray[dataPoint]=dataPointView();
507          }
508    
509          if (dataPointRank==1) {
510            for (int i=0; i<dataPointShape[0]; i++) {
511              numArray[dataPoint][i]=dataPointView(i);
512            }
513          }
514    
515          if (dataPointRank==2) {
516            for (int i=0; i<dataPointShape[0]; i++) {
517              for (int j=0; j<dataPointShape[1]; j++) {
518                numArray[dataPoint][i][j] = dataPointView(i,j);
519              }
520            }
521          }
522    
523          if (dataPointRank==3) {
524            for (int i=0; i<dataPointShape[0]; i++) {
525              for (int j=0; j<dataPointShape[1]; j++) {
526                for (int k=0; k<dataPointShape[2]; k++) {
527                  numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
528                }
529              }
530            }
531          }
532    
533          if (dataPointRank==4) {
534            for (int i=0; i<dataPointShape[0]; i++) {
535              for (int j=0; j<dataPointShape[1]; j++) {
536                for (int k=0; k<dataPointShape[2]; k++) {
537                  for (int l=0; l<dataPointShape[3]; l++) {
538                    numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
539                  }
540                }
541              }
542            }
543          }
544    
545          dataPoint++;
546    
547        }
548      }
549    
550      //
551      // return the loaded array
552      return numArray;
553    }
554    
555    boost::python::numeric::array
556  Data::integrate() const  Data::integrate() const
557  {  {
558    int index;    int index;
# Line 440  Data::integrate() const Line 569  Data::integrate() const
569    // and load the array with the integral values    // and load the array with the integral values
570    boost::python::numeric::array bp_array(1.0);    boost::python::numeric::array bp_array(1.0);
571    if (rank==0) {    if (rank==0) {
572        bp_array.resize(1);
573      index = 0;      index = 0;
574      bp_array[0] = integrals[index];      bp_array[0] = integrals[index];
575    }    }
# Line 519  Data::ln() const Line 649  Data::ln() const
649    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
650  }  }
651    
652    Data
653    Data::sign() const
654    {
655      return escript::unaryOp(*this,escript::fsign);
656    }
657    
658    Data
659    Data::abs() const
660    {
661      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
662    }
663    
664    Data
665    Data::neg() const
666    {
667      return escript::unaryOp(*this,negate<double>());
668    }
669    
670    Data
671    Data::pos() const
672    {
673      return (*this);
674    }
675    
676    Data
677    Data::exp() const
678    {
679      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
680    }
681    
682    Data
683    Data::sqrt() const
684    {
685      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
686    }
687    
688  double  double
689  Data::Lsup() const  Data::Lsup() const
690  {  {
691    //    //
692    // set the initial absolute maximum value to min possible    // set the initial absolute maximum value to zero
693    return algorithm(DataAlgorithmAdapter<AbsMax>(numeric_limits<double>::min()));    return algorithm(DataAlgorithmAdapter<AbsMax>(0));
694    }
695    
696    double
697    Data::Linf() const
698    {
699      //
700      // set the initial absolute minimum value to max double
701      return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));
702  }  }
703    
704  double  double
705  Data::sup() const  Data::sup() const
706  {  {
707    //    //
708    // set the initial maximum value to min possible    // set the initial maximum value to min possible double
709    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::min()));    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
710  }  }
711    
712  double  double
713  Data::inf() const  Data::inf() const
714  {  {
715    //    //
716    // set the initial minimum value to max possible    // set the initial minimum value to max possible double
717    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
718  }  }
719    
720  Data& Data::operator+=(const Data& right)  Data
721    Data::maxval() const
722    {
723      //
724      // set the initial maximum value to min possible double
725      return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));
726    }
727    
728    Data
729    Data::minval() const
730    {
731      //
732      // set the initial minimum value to max possible double
733      return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));
734    }
735    
736    Data
737    Data::length() const
738    {
739      return dp_algorithm(DataAlgorithmAdapter<Length>(0));
740    }
741    
742    Data
743    Data::trace() const
744    {
745      return dp_algorithm(DataAlgorithmAdapter<Trace>(0));
746    }
747    
748    Data
749    Data::transpose(int axis) const
750    {
751      // not implemented
752      throw DataException("Error - Data::transpose not implemented yet.");
753      return Data();
754    }
755    
756    void
757    Data::saveDX(std::string fileName) const
758    {
759      getDomain().saveDX(fileName,*this);
760      return;
761    }
762    
763    void
764    Data::saveVTK(std::string fileName) const
765    {
766      getDomain().saveVTK(fileName,*this);
767      return;
768    }
769    
770    Data&
771    Data::operator+=(const Data& right)
772  {  {
773    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
774    return (*this);    return (*this);
775  }  }
776    
777  Data& Data::operator+=(const boost::python::object& right)  Data&
778    Data::operator+=(const boost::python::object& right)
779  {  {
780    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
781    return (*this);    return (*this);
782  }  }
783    
784  Data& Data::operator-=(const Data& right)  Data&
785    Data::operator-=(const Data& right)
786  {  {
787    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
788    return (*this);    return (*this);
789  }  }
790    
791  Data& Data::operator-=(const boost::python::object& right)  Data&
792    Data::operator-=(const boost::python::object& right)
793  {  {
794    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
795    return (*this);    return (*this);
796  }  }
797    
798  Data& Data::operator*=(const Data& right)  Data&
799    Data::operator*=(const Data& right)
800  {  {
801    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
802    return (*this);    return (*this);
803  }  }
804    
805  Data& Data::operator*=(const boost::python::object& right)  Data&
806    Data::operator*=(const boost::python::object& right)
807  {  {
808    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
809    return (*this);    return (*this);
810  }  }
811    
812  Data& Data::operator/=(const Data& right)  Data&
813    Data::operator/=(const Data& right)
814  {  {
815    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
816    return (*this);    return (*this);
817  }  }
818    
819  Data& Data::operator/=(const boost::python::object& right)  Data&
820    Data::operator/=(const boost::python::object& right)
821  {  {
822    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
823    return (*this);    return (*this);
824  }  }
825    
826  Data Data::powO(const boost::python::object& right) const  Data
827    Data::powO(const boost::python::object& right) const
828  {  {
829    Data result;    Data result;
830    result.copy(*this);    result.copy(*this);
# Line 599  Data Data::powO(const boost::python::obj Line 832  Data Data::powO(const boost::python::obj
832    return result;    return result;
833  }  }
834    
835  Data Data::powD(const Data& right) const  Data
836    Data::powD(const Data& right) const
837  {  {
838    Data result;    Data result;
839    result.copy(*this);    result.copy(*this);
# Line 609  Data Data::powD(const Data& right) const Line 843  Data Data::powD(const Data& right) const
843    
844  //  //
845  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namepsace this operator belongs to
846  Data escript::operator+(const Data& left, const Data& right)  Data
847    escript::operator+(const Data& left, const Data& right)
848  {  {
849    Data result;    Data result;
850    //    //
# Line 621  Data escript::operator+(const Data& left Line 856  Data escript::operator+(const Data& left
856    
857  //  //
858  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namepsace this operator belongs to
859  Data escript::operator-(const Data& left, const Data& right)  Data
860    escript::operator-(const Data& left, const Data& right)
861  {  {
862    Data result;    Data result;
863    //    //
# Line 633  Data escript::operator-(const Data& left Line 869  Data escript::operator-(const Data& left
869    
870  //  //
871  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namepsace this operator belongs to
872  Data escript::operator*(const Data& left, const Data& right)  Data
873    escript::operator*(const Data& left, const Data& right)
874  {  {
875    Data result;    Data result;
876    //    //
# Line 645  Data escript::operator*(const Data& left Line 882  Data escript::operator*(const Data& left
882    
883  //  //
884  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namepsace this operator belongs to
885  Data escript::operator/(const Data& left, const Data& right)  Data
886    escript::operator/(const Data& left, const Data& right)
887  {  {
888    Data result;    Data result;
889    //    //
# Line 657  Data escript::operator/(const Data& left Line 895  Data escript::operator/(const Data& left
895    
896  //  //
897  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namepsace this operator belongs to
898  Data escript::operator+(const Data& left, const boost::python::object& right)  Data
899    escript::operator+(const Data& left, const boost::python::object& right)
900  {  {
901    //    //
902    // Convert to DataArray format if possible    // Convert to DataArray format if possible
# Line 672  Data escript::operator+(const Data& left Line 911  Data escript::operator+(const Data& left
911    
912  //  //
913  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namepsace this operator belongs to
914  Data escript::operator-(const Data& left, const boost::python::object& right)  Data
915    escript::operator-(const Data& left, const boost::python::object& right)
916  {  {
917    //    //
918    // Convert to DataArray format if possible    // Convert to DataArray format if possible
# Line 687  Data escript::operator-(const Data& left Line 927  Data escript::operator-(const Data& left
927    
928  //  //
929  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namepsace this operator belongs to
930  Data escript::operator*(const Data& left, const boost::python::object& right)  Data
931    escript::operator*(const Data& left, const boost::python::object& right)
932  {  {
933    //    //
934    // Convert to DataArray format if possible    // Convert to DataArray format if possible
# Line 702  Data escript::operator*(const Data& left Line 943  Data escript::operator*(const Data& left
943    
944  //  //
945  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namepsace this operator belongs to
946  Data escript::operator/(const Data& left, const boost::python::object& right)  Data
947    escript::operator/(const Data& left, const boost::python::object& right)
948  {  {
949    //    //
950    // Convert to DataArray format if possible    // Convert to DataArray format if possible
# Line 717  Data escript::operator/(const Data& left Line 959  Data escript::operator/(const Data& left
959    
960  //  //
961  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namepsace this operator belongs to
962  Data escript::operator+(const boost::python::object& left, const Data& right)  Data
963    escript::operator+(const boost::python::object& left, const Data& right)
964  {  {
965    //    //
966    // Construct the result using the given value and the other parameters    // Construct the result using the given value and the other parameters
# Line 729  Data escript::operator+(const boost::pyt Line 972  Data escript::operator+(const boost::pyt
972    
973  //  //
974  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namepsace this operator belongs to
975  Data escript::operator-(const boost::python::object& left, const Data& right)  Data
976    escript::operator-(const boost::python::object& left, const Data& right)
977  {  {
978    //    //
979    // Construct the result using the given value and the other parameters    // Construct the result using the given value and the other parameters
# Line 741  Data escript::operator-(const boost::pyt Line 985  Data escript::operator-(const boost::pyt
985    
986  //  //
987  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namepsace this operator belongs to
988  Data escript::operator*(const boost::python::object& left, const Data& right)  Data
989    escript::operator*(const boost::python::object& left, const Data& right)
990  {  {
991    //    //
992    // Construct the result using the given value and the other parameters    // Construct the result using the given value and the other parameters
# Line 753  Data escript::operator*(const boost::pyt Line 998  Data escript::operator*(const boost::pyt
998    
999  //  //
1000  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namepsace this operator belongs to
1001  Data escript::operator/(const boost::python::object& left, const Data& right)  Data
1002    escript::operator/(const boost::python::object& left, const Data& right)
1003  {  {
1004    //    //
1005    // Construct the result using the given value and the other parameters    // Construct the result using the given value and the other parameters
# Line 765  Data escript::operator/(const boost::pyt Line 1011  Data escript::operator/(const boost::pyt
1011    
1012  //  //
1013  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namepsace this operator belongs to
1014  bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1015  {  //{
1016    /*  //  /*
1017    NB: this operator does very little at this point, and isn't to  //  NB: this operator does very little at this point, and isn't to
1018    be relied on. Requires further implementation.  //  be relied on. Requires further implementation.
1019    */  //  */
1020    //
1021    bool ret;  //  bool ret;
1022    //
1023    if (left.isEmpty()) {  //  if (left.isEmpty()) {
1024      if(!right.isEmpty()) {  //    if(!right.isEmpty()) {
1025        ret = false;  //      ret = false;
1026      } else {  //    } else {
1027        ret = true;  //      ret = true;
1028      }  //    }
1029    }  //  }
1030    //
1031    //  if (left.isConstant()) {
1032    //    if(!right.isConstant()) {
1033    //      ret = false;
1034    //    } else {
1035    //      ret = true;
1036    //    }
1037    // }
1038    //
1039    //  if (left.isTagged()) {
1040    //   if(!right.isTagged()) {
1041    //      ret = false;
1042    //    } else {
1043    //      ret = true;
1044    //    }
1045    //  }
1046    //
1047    //  if (left.isExpanded()) {
1048    //    if(!right.isExpanded()) {
1049    //      ret = false;
1050    //    } else {
1051    //      ret = true;
1052    //    }
1053    //  }
1054    //
1055    //  return ret;
1056    //}
1057    
1058    if (left.isConstant()) {  Data
1059      if(!right.isConstant()) {  Data::getItem(const boost::python::object& key) const
1060        ret = false;  {
1061      } else {    const DataArrayView& view=getPointDataView();
       ret = true;  
     }  
   }  
1062    
1063    if (left.isTagged()) {    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
     if(!right.isTagged()) {  
       ret = false;  
     } else {  
       ret = true;  
     }  
   }  
1064    
1065    if (left.isExpanded()) {    if (slice_region.size()!=view.getRank()) {
1066      if(!right.isExpanded()) {      throw DataException("Error - slice size does not match Data rank.");
       ret = false;  
     } else {  
       ret = true;  
     }  
1067    }    }
1068    
1069    return ret;    return getSlice(slice_region);
1070  }  }
1071    
1072  Data  Data
1073  Data::getItem(const boost::python::object& key) const  Data::getSlice(const DataArrayView::RegionType& region) const
1074    {
1075      return Data(*this,region);
1076    }
1077    
1078    void
1079    Data::setItemO(const boost::python::object& key,
1080                   const boost::python::object& value)
1081  {  {
1082     const DataArrayView& view=getPointDataView();    Data tempData(value,getFunctionSpace());
1083     DataArrayView::RegionType slice_region=view.getSliceRegion(key);    setItemD(key,tempData);
    if (slice_region.size()!=view.getRank()) {  
      throw DataException("Error - slice size does not match Data rank.");  
      return Data();  
    }  
    //  
    // Create a new Data which is a slice of this one  
    return getSlice(slice_region);  
1084  }  }
1085    
1086  void  void
1087  Data::setItem(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
1088                const Data& value)                 const Data& value)
1089  {  {
1090    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1091    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1092    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1093      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
1094    }    }
1095    typeMatch(value);    if (getFunctionSpace()!=value.getFunctionSpace()) {
1096    setSlice(value,slice_region);       setSlice(Data(value,getFunctionSpace()),slice_region);
1097      } else {
1098         setSlice(value,slice_region);
1099      }
1100    }
1101    
1102    void
1103    Data::setSlice(const Data& value,
1104                   const DataArrayView::RegionType& region)
1105    {
1106      Data tempValue(value);
1107      typeMatchLeft(tempValue);
1108      typeMatchRight(tempValue);
1109      m_data->setSlice(tempValue.m_data.get(),region);
1110  }  }
1111    
1112  void  void
1113  Data::typeMatch(const Data& right)  Data::typeMatchLeft(Data& right) const
1114    {
1115      if (isExpanded()){
1116        right.expand();
1117      } else if (isTagged()) {
1118        if (right.isConstant()) {
1119          right.tag();
1120        }
1121      }
1122    }
1123    
1124    void
1125    Data::typeMatchRight(const Data& right)
1126  {  {
   //  
   // match the type of this to the RHS  
1127    if (isTagged()) {    if (isTagged()) {
1128      if (right.isExpanded()) {      if (right.isExpanded()) {
       //  
       // if the right hand side is expanded so must this  
1129        expand();        expand();
1130      }      }
1131    } else if (isConstant()) {    } else if (isConstant()) {
1132      if (right.isExpanded()) {      if (right.isExpanded()) {
       //  
       // if the right hand side is expanded so must this  
1133        expand();        expand();
1134      } else if (right.isTagged()) {      } else if (right.isTagged()) {
       //  
       // if the right hand side is tagged so must this  
1135        tag();        tag();
1136      }      }
1137    }    }
# Line 881  Data::setTaggedValue(int tagKey, Line 1158  Data::setTaggedValue(int tagKey,
1158    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
1159  }  }
1160    
1161    void
1162    Data::setRefValue(int ref,
1163                      const boost::python::numeric::array& value)
1164    {
1165      //
1166      // Construct DataArray from boost::python::object input value
1167      DataArray valueDataArray(value);
1168    
1169      //
1170      // Call DataAbstract::setRefValue
1171      m_data->setRefValue(ref,valueDataArray);
1172    }
1173    
1174    void
1175    Data::getRefValue(int ref,
1176                      boost::python::numeric::array& value)
1177    {
1178      //
1179      // Construct DataArray for boost::python::object return value
1180      DataArray valueDataArray(value);
1181    
1182      //
1183      // Load DataArray with values from data-points specified by ref
1184      m_data->getRefValue(ref,valueDataArray);
1185    
1186      //
1187      // Load values from valueDataArray into return numarray
1188    
1189      // extract the shape of the numarray
1190      int rank = value.getrank();
1191      DataArrayView::ShapeType shape;
1192      for (int i=0; i < rank; i++) {
1193        shape.push_back(extract<int>(value.getshape()[i]));
1194      }
1195    
1196      // and load the numarray with the data from the DataArray
1197      DataArrayView valueView = valueDataArray.getView();
1198    
1199      if (rank==0) {
1200        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1201      }
1202      if (rank==1) {
1203        for (int i=0; i < shape[0]; i++) {
1204          value[i] = valueView(i);
1205        }
1206      }
1207      if (rank==2) {
1208        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1209      }
1210      if (rank==3) {
1211        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1212      }
1213      if (rank==4) {
1214        throw DataException("Data::getRefValue error: only rank 1 data handled for now.");
1215      }
1216    
1217    }
1218    
1219  /*  /*
1220  Note: this version removed for now. Not needed, and breaks escript.cpp  Note: this version removed for now. Not needed, and breaks escript.cpp
1221  void  void
# Line 901  Data::setTaggedValue(int tagKey, Line 1236  Data::setTaggedValue(int tagKey,
1236  }  }
1237  */  */
1238    
1239    void
1240    Data::archiveData(const std::string fileName)
1241    {
1242      cout << "Archiving Data object to: " << fileName << endl;
1243    
1244      //
1245      // Determine type of this Data object
1246      int dataType = -1;
1247    
1248      if (isEmpty()) {
1249        dataType = 0;
1250        cout << "\tdataType: DataEmpty" << endl;
1251      }
1252      if (isConstant()) {
1253        dataType = 1;
1254        cout << "\tdataType: DataConstant" << endl;
1255      }
1256      if (isTagged()) {
1257        dataType = 2;
1258        cout << "\tdataType: DataTagged" << endl;
1259      }
1260      if (isExpanded()) {
1261        dataType = 3;
1262        cout << "\tdataType: DataExpanded" << endl;
1263      }
1264      if (dataType == -1) {
1265        throw DataException("archiveData Error: undefined dataType");
1266      }
1267    
1268      //
1269      // Collect data items common to all Data types
1270      int noSamples = getNumSamples();
1271      int noDPPSample = getNumDataPointsPerSample();
1272      int functionSpaceType = getFunctionSpace().getTypeCode();
1273      int dataPointRank = getDataPointRank();
1274      int dataPointSize = getDataPointSize();
1275      int dataLength = getLength();
1276      DataArrayView::ShapeType dataPointShape = getDataPointShape();
1277      int referenceNumbers[noSamples];
1278      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1279        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1280      }
1281      int tagNumbers[noSamples];
1282      if (isTagged()) {
1283        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1284          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1285        }
1286      }
1287    
1288      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1289      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1290      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1291    
1292      //
1293      // Flatten Shape to an array of integers suitable for writing to file
1294      int flatShape[4] = {0,0,0,0};
1295      cout << "\tshape: < ";
1296      for (int dim=0; dim<dataPointRank; dim++) {
1297        flatShape[dim] = dataPointShape[dim];
1298        cout << dataPointShape[dim] << " ";
1299      }
1300      cout << ">" << endl;
1301    
1302      //
1303      // Write common data items to archive file
1304      ofstream archiveFile;
1305      archiveFile.open(fileName.data(), ios::out);
1306    
1307      if (!archiveFile.good()) {
1308        throw DataException("archiveData Error: problem opening archive file");
1309      }
1310    
1311      archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1312      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1313      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1314      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1315      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1316      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1317      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1318      for (int dim = 0; dim < 4; dim++) {
1319        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1320      }
1321      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1322        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1323      }
1324      if (isTagged()) {
1325        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1326          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1327        }
1328      }
1329    
1330      if (!archiveFile.good()) {
1331        throw DataException("archiveData Error: problem writing to archive file");
1332      }
1333    
1334      archiveFile.close();
1335    
1336      if (!archiveFile.good()) {
1337        throw DataException("archiveData Error: problem closing archive file");
1338      }
1339    
1340      //
1341      // Collect and archive underlying data values for each Data type
1342      switch (dataType) {
1343        case 0:
1344          // DataEmpty
1345          break;
1346        case 1:
1347          // DataConstant
1348          break;
1349        case 2:
1350          // DataTagged
1351          break;
1352        case 3:
1353          // DataExpanded
1354          break;
1355      }
1356    
1357    }
1358    
1359    void
1360    Data::extractData(const std::string fileName,
1361                      const FunctionSpace& fspace)
1362    {
1363      //
1364      // Can only extract Data to an object which is initially DataEmpty
1365      if (!isEmpty()) {
1366        throw DataException("extractData Error: can only extract to DataEmpty object");
1367      }
1368    
1369      cout << "Extracting Data object from: " << fileName << endl;
1370    
1371      int dataType;
1372      int noSamples;
1373      int noDPPSample;
1374      int functionSpaceType;
1375      int dataPointRank;
1376      int dataPointSize;
1377      int dataLength;
1378      DataArrayView::ShapeType dataPointShape;
1379      int flatShape[4];
1380    
1381      //
1382      // Open the archive file and read common data items
1383      ifstream archiveFile;
1384      archiveFile.open(fileName.data(), ios::in);
1385    
1386      if (!archiveFile.good()) {
1387        throw DataException("extractData Error: problem opening archive file");
1388      }
1389    
1390      archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1391      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1392      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1393      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1394      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1395      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1396      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1397      for (int dim = 0; dim < 4; dim++) {
1398        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1399        if (flatShape[dim]>0) {
1400          dataPointShape.push_back(flatShape[dim]);
1401        }
1402      }
1403      int referenceNumbers[noSamples];
1404      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1405        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1406      }
1407      int tagNumbers[noSamples];
1408      if (dataType==2) {
1409        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1410          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1411        }
1412      }
1413    
1414      if (!archiveFile.good()) {
1415        throw DataException("extractData Error: problem reading from archive file");
1416      }
1417    
1418      archiveFile.close();
1419    
1420      if (!archiveFile.good()) {
1421        throw DataException("extractData Error: problem closing archive file");
1422      }
1423    
1424      switch (dataType) {
1425        case 0:
1426          cout << "\tdataType: DataEmpty" << endl;
1427          break;
1428        case 1:
1429          cout << "\tdataType: DataConstant" << endl;
1430          break;
1431        case 2:
1432          cout << "\tdataType: DataTagged" << endl;
1433          break;
1434        case 3:
1435          cout << "\tdataType: DataExpanded" << endl;
1436          break;
1437        default:
1438          throw DataException("extractData Error: undefined dataType read from archive file");
1439          break;
1440      }
1441    
1442      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1443      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1444      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1445      cout << "\tshape: < ";
1446      for (int dim = 0; dim < dataPointRank; dim++) {
1447        cout << dataPointShape[dim] << " ";
1448      }
1449      cout << ">" << endl;
1450    
1451      //
1452      // Verify that supplied FunctionSpace object is compatible with this Data object.
1453      if ( (fspace.getTypeCode()!=functionSpaceType) ||
1454           (fspace.getNumSamples()!=noSamples) ||
1455           (fspace.getNumDPPSample()!=noDPPSample)
1456         ) {
1457        throw DataException("extractData Error: incompatible FunctionSpace");
1458      }
1459      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1460        if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
1461          throw DataException("extractData Error: incompatible FunctionSpace");
1462        }
1463      }
1464      if (dataType==2) {
1465        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1466          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
1467            throw DataException("extractData Error: incompatible FunctionSpace");
1468          }
1469        }
1470      }
1471    
1472      //
1473      // Construct a DataVector to hold underlying data values
1474      DataVector dataVec(dataLength);
1475    
1476      //
1477      // Load this DataVector with the appropriate values
1478      switch (dataType) {
1479        case 0:
1480          // DataEmpty
1481          break;
1482        case 1:
1483          // DataConstant
1484          break;
1485        case 2:
1486          // DataTagged
1487          break;
1488        case 3:
1489          // DataExpanded
1490          break;
1491      }
1492    
1493      //
1494      // Construct an appropriate Data object
1495      DataAbstract* tempData;
1496      switch (dataType) {
1497        case 0:
1498          // DataEmpty
1499          tempData=new DataEmpty();
1500          break;
1501        case 1:
1502          // DataConstant
1503          tempData=new DataConstant(fspace,dataPointShape,dataVec);
1504          break;
1505        case 2:
1506          // DataTagged
1507          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
1508          break;
1509        case 3:
1510          // DataExpanded
1511          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
1512          break;
1513      }
1514      shared_ptr<DataAbstract> temp_data(tempData);
1515      m_data=temp_data;
1516    
1517    }
1518    
1519  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
1520  {  {
1521    o << data.toString();    o << data.toString();

Legend:
Removed from v.94  
changed lines
  Added in v.119

  ViewVC Help
Powered by ViewVC 1.1.26