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

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

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

revision 94 by jgs, Wed Oct 27 00:45:54 2004 UTC revision 151 by jgs, Thu Sep 22 01:55:00 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"
42    
43  using namespace std;  using namespace std;
44  using namespace boost::python;  using namespace boost::python;
45  using namespace boost;  using namespace boost;
46  using namespace escript;  using namespace escript;
47    
48    #if defined DOPROF
49    //
50    // global table of profiling data for all Data objects
51    DataProf dataProfTable;
52    #endif
53    
54  Data::Data()  Data::Data()
55  {  {
56    //    //
57    // Default data is type DataEmpty    // Default data is type DataEmpty
58    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
59    m_data=shared_ptr<DataAbstract>(temp);    shared_ptr<DataAbstract> temp_data(temp);
60      m_data=temp_data;
61    #if defined DOPROF
62      // create entry in global profiling table for this object
63      profData = dataProfTable.newData();
64    #endif
65  }  }
66    
67  Data::Data(double value,  Data::Data(double value,
# Line 63  Data::Data(double value, Line 75  Data::Data(double value,
75    }    }
76    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
77    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
78    #if defined DOPROF
79      // create entry in global profiling table for this object
80      profData = dataProfTable.newData();
81    #endif
82  }  }
83    
84  Data::Data(double value,  Data::Data(double value,
# Line 73  Data::Data(double value, Line 89  Data::Data(double value,
89    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
90    pair<int,int> dataShape=what.getDataShape();    pair<int,int> dataShape=what.getDataShape();
91    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
92    #if defined DOPROF
93      // create entry in global profiling table for this object
94      profData = dataProfTable.newData();
95    #endif
96  }  }
97    
98  Data::Data(const Data& inData):  Data::Data(const Data& inData)
99    m_data(inData.m_data)  {
100  {    m_data=inData.m_data;
101    #if defined DOPROF
102      // create entry in global profiling table for this object
103      profData = dataProfTable.newData();
104    #endif
105  }  }
106    
107  Data::Data(const Data& inData,  Data::Data(const Data& inData,
108             const DataArrayView::RegionType& region)             const DataArrayView::RegionType& region)
109  {  {
110    //    //
111    // Create data which is a subset(slice) of another Data    // Create Data which is a slice of another Data
112    DataAbstract* tmp=inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
113    m_data=shared_ptr<DataAbstract>(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
114      m_data=temp_data;
115    #if defined DOPROF
116      // create entry in global profiling table for this object
117      profData = dataProfTable.newData();
118    #endif
119  }  }
120    
121  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 96  Data::Data(const Data& inData, Line 125  Data::Data(const Data& inData,
125      m_data=inData.m_data;      m_data=inData.m_data;
126    } else {    } else {
127      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
128      // Note for Lutz, Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
129      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
130      // be able to create an instance of AbstractDomain but that was done      // be able to create an instance of AbstractDomain but that was done
131      // as a boost python work around which may no longer be required.      // as a boost:python work around which may no longer be required.
132      const AbstractDomain& inDataDomain=inData.getDomain();      const AbstractDomain& inDataDomain=inData.getDomain();
133      if  (inDataDomain==functionspace.getDomain()) {      if  (inDataDomain==functionspace.getDomain()) {
134        inDataDomain.interpolateOnDomain(tmp,inData);        inDataDomain.interpolateOnDomain(tmp,inData);
# Line 108  Data::Data(const Data& inData, Line 137  Data::Data(const Data& inData,
137      }      }
138      m_data=tmp.m_data;      m_data=tmp.m_data;
139    }    }
140    #if defined DOPROF
141      // create entry in global profiling table for this object
142      profData = dataProfTable.newData();
143    #endif
144  }  }
145    
146  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(const DataTagged::TagListType& tagKeys,
# Line 117  Data::Data(const DataTagged::TagListType Line 150  Data::Data(const DataTagged::TagListType
150             bool expanded)             bool expanded)
151  {  {
152    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
153    m_data=shared_ptr<DataAbstract>(temp);    shared_ptr<DataAbstract> temp_data(temp);
154      m_data=temp_data;
155    if (expanded) {    if (expanded) {
156      expand();      expand();
157    }    }
158    #if defined DOPROF
159      // create entry in global profiling table for this object
160      profData = dataProfTable.newData();
161    #endif
162  }  }
163    
164  Data::Data(const numeric::array& value,  Data::Data(const numeric::array& value,
# Line 128  Data::Data(const numeric::array& value, Line 166  Data::Data(const numeric::array& value,
166             bool expanded)             bool expanded)
167  {  {
168    initialise(value,what,expanded);    initialise(value,what,expanded);
169    #if defined DOPROF
170      // create entry in global profiling table for this object
171      profData = dataProfTable.newData();
172    #endif
173  }  }
174    
175  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
# Line 135  Data::Data(const DataArrayView& value, Line 177  Data::Data(const DataArrayView& value,
177             bool expanded)             bool expanded)
178  {  {
179    initialise(value,what,expanded);    initialise(value,what,expanded);
180    #if defined DOPROF
181      // create entry in global profiling table for this object
182      profData = dataProfTable.newData();
183    #endif
184  }  }
185    
186  Data::Data(const object& value,  Data::Data(const object& value,
# Line 143  Data::Data(const object& value, Line 189  Data::Data(const object& value,
189  {  {
190    numeric::array asNumArray(value);    numeric::array asNumArray(value);
191    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
192    #if defined DOPROF
193      // create entry in global profiling table for this object
194      profData = dataProfTable.newData();
195    #endif
196  }  }
197    
198  Data::Data(const object& value,  Data::Data(const object& value,
# Line 163  Data::Data(const object& value, Line 213  Data::Data(const object& value,
213      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
214      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
215    }    }
216    #if defined DOPROF
217      // create entry in global profiling table for this object
218      profData = dataProfTable.newData();
219    #endif
220    }
221    
222    Data::~Data()
223    {
224    
225  }  }
226    
227  escriptDataC  escriptDataC
# Line 181  Data::getDataC() const Line 240  Data::getDataC() const
240    return temp;    return temp;
241  }  }
242    
243  tuple  const boost::python::tuple
244  Data::getShapeTuple() const  Data::getShapeTuple() const
245  {  {
246    const DataArrayView::ShapeType& shape=getDataPointShape();    const DataArrayView::ShapeType& shape=getDataPointShape();
# Line 212  Data::copy(const Data& other) Line 271  Data::copy(const Data& other)
271        //        //
272        // Construct a DataExpanded copy        // Construct a DataExpanded copy
273        DataAbstract* newData=new DataExpanded(*temp);        DataAbstract* newData=new DataExpanded(*temp);
274        m_data=shared_ptr<DataAbstract>(newData);        shared_ptr<DataAbstract> temp_data(newData);
275          m_data=temp_data;
276        return;        return;
277      }      }
278    }    }
# Line 220  Data::copy(const Data& other) Line 280  Data::copy(const Data& other)
280      DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());      DataTagged* temp=dynamic_cast<DataTagged*>(other.m_data.get());
281      if (temp!=0) {      if (temp!=0) {
282        //        //
283        // Construct a DataTaggeded copy        // Construct a DataTagged copy
284        DataAbstract* newData=new DataTagged(*temp);        DataAbstract* newData=new DataTagged(*temp);
285        m_data=shared_ptr<DataAbstract>(newData);        shared_ptr<DataAbstract> temp_data(newData);
286          m_data=temp_data;
287        return;        return;
288      }      }
289    }    }
# Line 232  Data::copy(const Data& other) Line 293  Data::copy(const Data& other)
293        //        //
294        // Construct a DataConstant copy        // Construct a DataConstant copy
295        DataAbstract* newData=new DataConstant(*temp);        DataAbstract* newData=new DataConstant(*temp);
296        m_data=shared_ptr<DataAbstract>(newData);        shared_ptr<DataAbstract> temp_data(newData);
297          m_data=temp_data;
298          return;
299        }
300      }
301      {
302        DataEmpty* temp=dynamic_cast<DataEmpty*>(other.m_data.get());
303        if (temp!=0) {
304          //
305          // Construct a DataEmpty copy
306          DataAbstract* newData=new DataEmpty();
307          shared_ptr<DataAbstract> temp_data(newData);
308          m_data=temp_data;
309        return;        return;
310      }      }
311    }    }
# Line 284  Data::isConstant() const Line 357  Data::isConstant() const
357    return (temp!=0);    return (temp!=0);
358  }  }
359    
 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);  
 }  
   
360  void  void
361  Data::expand()  Data::expand()
362  {  {
363    if (isConstant()) {    if (isConstant()) {
364      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
365      DataAbstract* temp=new DataExpanded(*tempDataConst);      DataAbstract* temp=new DataExpanded(*tempDataConst);
366      m_data=shared_ptr<DataAbstract>(temp);      shared_ptr<DataAbstract> temp_data(temp);
367        m_data=temp_data;
368    } else if (isTagged()) {    } else if (isTagged()) {
369      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());      DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
370      DataAbstract* temp=new DataExpanded(*tempDataTag);      DataAbstract* temp=new DataExpanded(*tempDataTag);
371      m_data=shared_ptr<DataAbstract>(temp);      shared_ptr<DataAbstract> temp_data(temp);
372        m_data=temp_data;
373    } else if (isExpanded()) {    } else if (isExpanded()) {
374      //      //
375      // do nothing      // do nothing
# Line 319  Data::expand() Line 381  Data::expand()
381  }  }
382    
383  void  void
 Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  
 {  
   m_data->reshapeDataPoint(shape);  
 }  
   
 void  
384  Data::tag()  Data::tag()
385  {  {
386    if (isConstant()) {    if (isConstant()) {
387      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());      DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
388      DataAbstract* temp=new DataTagged(*tempDataConst);      DataAbstract* temp=new DataTagged(*tempDataConst);
389      m_data=shared_ptr<DataAbstract>(temp);      shared_ptr<DataAbstract> temp_data(temp);
390        m_data=temp_data;
391    } else if (isTagged()) {    } else if (isTagged()) {
392      // do nothing      // do nothing
393    } else if (isExpanded()) {    } else if (isExpanded()) {
# Line 342  Data::tag() Line 399  Data::tag()
399    }    }
400  }  }
401    
402    void
403    Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)
404    {
405      m_data->reshapeDataPoint(shape);
406    }
407    
408  Data  Data
409  Data::wherePositive() const  Data::wherePositive() const
410  {  {
411    #if defined DOPROF
412      profData->where++;
413    #endif
414    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
415  }  }
416    
417  Data  Data
418    Data::whereNegative() const
419    {
420    #if defined DOPROF
421      profData->where++;
422    #endif
423      return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
424    }
425    
426    Data
427  Data::whereNonNegative() const  Data::whereNonNegative() const
428  {  {
429    #if defined DOPROF
430      profData->where++;
431    #endif
432    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
433  }  }
434    
435  Data  Data
436  Data::whereNegative() const  Data::whereNonPositive() const
437  {  {
438    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));  #if defined DOPROF
439      profData->where++;
440    #endif
441      return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
442  }  }
443    
444  Data  Data
445  Data::whereZero() const  Data::whereZero() const
446  {  {
447    #if defined DOPROF
448      profData->where++;
449    #endif
450    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));
451  }  }
452    
453  Data  Data
454    Data::whereNonZero() const
455    {
456    #if defined DOPROF
457      profData->where++;
458    #endif
459      return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));
460    }
461    
462    Data
463  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
464  {  {
465    #if defined DOPROF
466      profData->interpolate++;
467    #endif
468    return Data(*this,functionspace);    return Data(*this,functionspace);
469  }  }
470    
# Line 390  Data::probeInterpolation(const FunctionS Line 486  Data::probeInterpolation(const FunctionS
486  Data  Data
487  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
488  {  {
489    #if defined DOPROF
490      profData->grad++;
491    #endif
492    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
493      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
494    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
# Line 423  Data::getDataPointShape() const Line 522  Data::getDataPointShape() const
522    return getPointDataView().getShape();    return getPointDataView().getShape();
523  }  }
524    
525    void
526    Data::fillFromNumArray(const boost::python::numeric::array num_array)
527    {
528      //
529      // check rank
530      if (num_array.getrank()<getDataPointRank())
531          throw DataException("Rank of numarray does not match Data object rank");
532    
533      //
534      // check shape of num_array
535      for (int i=0; i<getDataPointRank(); i++) {
536        if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])
537           throw DataException("Shape of numarray does not match Data object rank");
538      }
539    
540      //
541      // make sure data is expanded:
542      if (!isExpanded()) {
543        expand();
544      }
545    
546      //
547      // and copy over
548      m_data->copyAll(num_array);
549    }
550    
551    const
552    boost::python::numeric::array
553    Data::convertToNumArray()
554    {
555      //
556      // determine the total number of data points
557      int numSamples = getNumSamples();
558      int numDataPointsPerSample = getNumDataPointsPerSample();
559      int numDataPoints = numSamples * numDataPointsPerSample;
560    
561      //
562      // determine the rank and shape of each data point
563      int dataPointRank = getDataPointRank();
564      DataArrayView::ShapeType dataPointShape = getDataPointShape();
565    
566      //
567      // create the numeric array to be returned
568      boost::python::numeric::array numArray(0.0);
569    
570      //
571      // the rank of the returned numeric array will be the rank of
572      // the data points, plus one. Where the rank of the array is n,
573      // the last n-1 dimensions will be equal to the shape of the
574      // data points, whilst the first dimension will be equal to the
575      // total number of data points. Thus the array will consist of
576      // a serial vector of the data points.
577      int arrayRank = dataPointRank + 1;
578      DataArrayView::ShapeType arrayShape;
579      arrayShape.push_back(numDataPoints);
580      for (int d=0; d<dataPointRank; d++) {
581         arrayShape.push_back(dataPointShape[d]);
582      }
583    
584      //
585      // resize the numeric array to the shape just calculated
586      if (arrayRank==1) {
587        numArray.resize(arrayShape[0]);
588      }
589      if (arrayRank==2) {
590        numArray.resize(arrayShape[0],arrayShape[1]);
591      }
592      if (arrayRank==3) {
593        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
594      }
595      if (arrayRank==4) {
596        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
597      }
598      if (arrayRank==5) {
599        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
600      }
601    
602      //
603      // loop through each data point in turn, loading the values for that data point
604      // into the numeric array.
605      int dataPoint = 0;
606      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
607        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
608          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
609          if (dataPointRank==0) {
610            numArray[dataPoint]=dataPointView();
611          }
612          if (dataPointRank==1) {
613            for (int i=0; i<dataPointShape[0]; i++) {
614              numArray[dataPoint][i]=dataPointView(i);
615            }
616          }
617          if (dataPointRank==2) {
618            for (int i=0; i<dataPointShape[0]; i++) {
619              for (int j=0; j<dataPointShape[1]; j++) {
620                numArray[dataPoint][i][j] = dataPointView(i,j);
621              }
622            }
623          }
624          if (dataPointRank==3) {
625            for (int i=0; i<dataPointShape[0]; i++) {
626              for (int j=0; j<dataPointShape[1]; j++) {
627                for (int k=0; k<dataPointShape[2]; k++) {
628                  numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
629                }
630              }
631            }
632          }
633          if (dataPointRank==4) {
634            for (int i=0; i<dataPointShape[0]; i++) {
635              for (int j=0; j<dataPointShape[1]; j++) {
636                for (int k=0; k<dataPointShape[2]; k++) {
637                  for (int l=0; l<dataPointShape[3]; l++) {
638                    numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
639                  }
640                }
641              }
642            }
643          }
644          dataPoint++;
645        }
646      }
647    
648      //
649      // return the loaded array
650      return numArray;
651    }
652    
653    const
654    boost::python::numeric::array
655    Data::convertToNumArrayFromSampleNo(int sampleNo)
656    {
657      //
658      // Check a valid sample number has been supplied
659      if (sampleNo >= getNumSamples()) {
660        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
661      }
662    
663      //
664      // determine the number of data points per sample
665      int numDataPointsPerSample = getNumDataPointsPerSample();
666    
667      //
668      // determine the rank and shape of each data point
669      int dataPointRank = getDataPointRank();
670      DataArrayView::ShapeType dataPointShape = getDataPointShape();
671    
672      //
673      // create the numeric array to be returned
674      boost::python::numeric::array numArray(0.0);
675    
676      //
677      // the rank of the returned numeric array will be the rank of
678      // the data points, plus one. Where the rank of the array is n,
679      // the last n-1 dimensions will be equal to the shape of the
680      // data points, whilst the first dimension will be equal to the
681      // total number of data points. Thus the array will consist of
682      // a serial vector of the data points.
683      int arrayRank = dataPointRank + 1;
684      DataArrayView::ShapeType arrayShape;
685      arrayShape.push_back(numDataPointsPerSample);
686      for (int d=0; d<dataPointRank; d++) {
687         arrayShape.push_back(dataPointShape[d]);
688      }
689    
690      //
691      // resize the numeric array to the shape just calculated
692      if (arrayRank==1) {
693        numArray.resize(arrayShape[0]);
694      }
695      if (arrayRank==2) {
696        numArray.resize(arrayShape[0],arrayShape[1]);
697      }
698      if (arrayRank==3) {
699        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
700      }
701      if (arrayRank==4) {
702        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
703      }
704      if (arrayRank==5) {
705        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
706      }
707    
708      //
709      // loop through each data point in turn, loading the values for that data point
710      // into the numeric array.
711      for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {
712        DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);
713        if (dataPointRank==0) {
714          numArray[dataPoint]=dataPointView();
715        }
716        if (dataPointRank==1) {
717          for (int i=0; i<dataPointShape[0]; i++) {
718            numArray[dataPoint][i]=dataPointView(i);
719          }
720        }
721        if (dataPointRank==2) {
722          for (int i=0; i<dataPointShape[0]; i++) {
723            for (int j=0; j<dataPointShape[1]; j++) {
724              numArray[dataPoint][i][j] = dataPointView(i,j);
725            }
726          }
727        }
728        if (dataPointRank==3) {
729          for (int i=0; i<dataPointShape[0]; i++) {
730            for (int j=0; j<dataPointShape[1]; j++) {
731              for (int k=0; k<dataPointShape[2]; k++) {
732                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
733              }
734            }
735          }
736        }
737        if (dataPointRank==4) {
738          for (int i=0; i<dataPointShape[0]; i++) {
739            for (int j=0; j<dataPointShape[1]; j++) {
740              for (int k=0; k<dataPointShape[2]; k++) {
741                for (int l=0; l<dataPointShape[3]; l++) {
742                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
743                }
744              }
745            }
746          }
747        }
748      }
749    
750      //
751      // return the loaded array
752      return numArray;
753    }
754    
755    const
756    boost::python::numeric::array
757    Data::convertToNumArrayFromDPNo(int sampleNo,
758                                    int dataPointNo)
759    {
760      //
761      // Check a valid sample number has been supplied
762      if (sampleNo >= getNumSamples()) {
763        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
764      }
765    
766      //
767      // Check a valid data point number has been supplied
768      if (dataPointNo >= getNumDataPointsPerSample()) {
769        throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");
770      }
771    
772      //
773      // determine the rank and shape of each data point
774      int dataPointRank = getDataPointRank();
775      DataArrayView::ShapeType dataPointShape = getDataPointShape();
776    
777      //
778      // create the numeric array to be returned
779      boost::python::numeric::array numArray(0.0);
780    
781      //
782      // the shape of the returned numeric array will be the same
783      // as that of the data point
784      int arrayRank = dataPointRank;
785      DataArrayView::ShapeType arrayShape = dataPointShape;
786    
787      //
788      // resize the numeric array to the shape just calculated
789      if (arrayRank==0) {
790        numArray.resize(1);
791      }
792      if (arrayRank==1) {
793        numArray.resize(arrayShape[0]);
794      }
795      if (arrayRank==2) {
796        numArray.resize(arrayShape[0],arrayShape[1]);
797      }
798      if (arrayRank==3) {
799        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
800      }
801      if (arrayRank==4) {
802        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
803      }
804    
805      //
806      // load the values for the data point into the numeric array.
807      DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
808      if (dataPointRank==0) {
809        numArray[0]=dataPointView();
810      }
811      if (dataPointRank==1) {
812        for (int i=0; i<dataPointShape[0]; i++) {
813          numArray[i]=dataPointView(i);
814        }
815      }
816      if (dataPointRank==2) {
817        for (int i=0; i<dataPointShape[0]; i++) {
818          for (int j=0; j<dataPointShape[1]; j++) {
819            numArray[i][j] = dataPointView(i,j);
820          }
821        }
822      }
823      if (dataPointRank==3) {
824        for (int i=0; i<dataPointShape[0]; i++) {
825          for (int j=0; j<dataPointShape[1]; j++) {
826            for (int k=0; k<dataPointShape[2]; k++) {
827              numArray[i][j][k]=dataPointView(i,j,k);
828            }
829          }
830        }
831      }
832      if (dataPointRank==4) {
833        for (int i=0; i<dataPointShape[0]; i++) {
834          for (int j=0; j<dataPointShape[1]; j++) {
835            for (int k=0; k<dataPointShape[2]; k++) {
836              for (int l=0; l<dataPointShape[3]; l++) {
837                numArray[i][j][k][l]=dataPointView(i,j,k,l);
838              }
839            }
840          }
841        }
842      }
843    
844      //
845      // return the loaded array
846      return numArray;
847    }
848    
849  boost::python::numeric::array  boost::python::numeric::array
850  Data::integrate() const  Data::integrate() const
851  {  {
# Line 430  Data::integrate() const Line 853  Data::integrate() const
853    int rank = getDataPointRank();    int rank = getDataPointRank();
854    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
855    
856    #if defined DOPROF
857      profData->integrate++;
858    #endif
859    
860    //    //
861    // calculate the integral values    // calculate the integral values
862    vector<double> integrals(getDataPointSize());    vector<double> integrals(getDataPointSize());
# Line 440  Data::integrate() const Line 867  Data::integrate() const
867    // and load the array with the integral values    // and load the array with the integral values
868    boost::python::numeric::array bp_array(1.0);    boost::python::numeric::array bp_array(1.0);
869    if (rank==0) {    if (rank==0) {
870        bp_array.resize(1);
871      index = 0;      index = 0;
872      bp_array[0] = integrals[index];      bp_array[0] = integrals[index];
873    }    }
# Line 492  Data::integrate() const Line 920  Data::integrate() const
920  Data  Data
921  Data::sin() const  Data::sin() const
922  {  {
923    #if defined DOPROF
924      profData->unary++;
925    #endif
926    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
927  }  }
928    
929  Data  Data
930  Data::cos() const  Data::cos() const
931  {  {
932    #if defined DOPROF
933      profData->unary++;
934    #endif
935    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
936  }  }
937    
938  Data  Data
939  Data::tan() const  Data::tan() const
940  {  {
941    #if defined DOPROF
942      profData->unary++;
943    #endif
944    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
945  }  }
946    
947  Data  Data
948    Data::asin() const
949    {
950    #if defined DOPROF
951      profData->unary++;
952    #endif
953      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
954    }
955    
956    Data
957    Data::acos() const
958    {
959    #if defined DOPROF
960      profData->unary++;
961    #endif
962      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
963    }
964    
965    Data
966    Data::atan() const
967    {
968    #if defined DOPROF
969      profData->unary++;
970    #endif
971      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
972    }
973    
974    Data
975    Data::sinh() const
976    {
977    #if defined DOPROF
978      profData->unary++;
979    #endif
980      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
981    }
982    
983    Data
984    Data::cosh() const
985    {
986    #if defined DOPROF
987      profData->unary++;
988    #endif
989      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
990    }
991    
992    Data
993    Data::tanh() const
994    {
995    #if defined DOPROF
996      profData->unary++;
997    #endif
998      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
999    }
1000    
1001    Data
1002    Data::asinh() const
1003    {
1004    #if defined DOPROF
1005      profData->unary++;
1006    #endif
1007      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
1008    }
1009    
1010    Data
1011    Data::acosh() const
1012    {
1013    #if defined DOPROF
1014      profData->unary++;
1015    #endif
1016      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
1017    }
1018    
1019    Data
1020    Data::atanh() const
1021    {
1022    #if defined DOPROF
1023      profData->unary++;
1024    #endif
1025      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
1026    }
1027    
1028    Data
1029  Data::log() const  Data::log() const
1030  {  {
1031    #if defined DOPROF
1032      profData->unary++;
1033    #endif
1034    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
1035  }  }
1036    
1037  Data  Data
1038  Data::ln() const  Data::ln() const
1039  {  {
1040    #if defined DOPROF
1041      profData->unary++;
1042    #endif
1043    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
1044  }  }
1045    
1046    Data
1047    Data::sign() const
1048    {
1049    #if defined DOPROF
1050      profData->unary++;
1051    #endif
1052      return escript::unaryOp(*this,escript::fsign);
1053    }
1054    
1055    Data
1056    Data::abs() const
1057    {
1058    #if defined DOPROF
1059      profData->unary++;
1060    #endif
1061      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
1062    }
1063    
1064    Data
1065    Data::neg() const
1066    {
1067    #if defined DOPROF
1068      profData->unary++;
1069    #endif
1070      return escript::unaryOp(*this,negate<double>());
1071    }
1072    
1073    Data
1074    Data::pos() const
1075    {
1076    #if defined DOPROF
1077      profData->unary++;
1078    #endif
1079      Data result;
1080      // perform a deep copy
1081      result.copy(*this);
1082      return result;
1083    }
1084    
1085    Data
1086    Data::exp() const
1087    {
1088    #if defined DOPROF
1089      profData->unary++;
1090    #endif
1091      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1092    }
1093    
1094    Data
1095    Data::sqrt() const
1096    {
1097    #if defined DOPROF
1098      profData->unary++;
1099    #endif
1100      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1101    }
1102    
1103  double  double
1104  Data::Lsup() const  Data::Lsup() const
1105  {  {
1106    //  #if defined DOPROF
1107    // set the initial absolute maximum value to min possible    profData->reduction1++;
1108    return algorithm(DataAlgorithmAdapter<AbsMax>(numeric_limits<double>::min()));  #endif
1109      //
1110      // set the initial absolute maximum value to zero
1111      AbsMax abs_max_func;
1112      return algorithm(abs_max_func,0);
1113    }
1114    
1115    double
1116    Data::Linf() const
1117    {
1118    #if defined DOPROF
1119      profData->reduction1++;
1120    #endif
1121      //
1122      // set the initial absolute minimum value to max double
1123      AbsMin abs_min_func;
1124      return algorithm(abs_min_func,numeric_limits<double>::max());
1125  }  }
1126    
1127  double  double
1128  Data::sup() const  Data::sup() const
1129  {  {
1130    //  #if defined DOPROF
1131    // set the initial maximum value to min possible    profData->reduction1++;
1132    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::min()));  #endif
1133      //
1134      // set the initial maximum value to min possible double
1135      FMax fmax_func;
1136      return algorithm(fmax_func,numeric_limits<double>::max()*-1);
1137  }  }
1138    
1139  double  double
1140  Data::inf() const  Data::inf() const
1141  {  {
1142    //  #if defined DOPROF
1143    // set the initial minimum value to max possible    profData->reduction1++;
1144    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));  #endif
1145      //
1146      // set the initial minimum value to max possible double
1147      FMin fmin_func;
1148      return algorithm(fmin_func,numeric_limits<double>::max());
1149  }  }
1150    
1151  Data& Data::operator+=(const Data& right)  Data
1152    Data::maxval() const
1153    {
1154    #if defined DOPROF
1155      profData->reduction2++;
1156    #endif
1157      //
1158      // set the initial maximum value to min possible double
1159      FMax fmax_func;
1160      return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1161    }
1162    
1163    Data
1164    Data::minval() const
1165    {
1166    #if defined DOPROF
1167      profData->reduction2++;
1168    #endif
1169      //
1170      // set the initial minimum value to max possible double
1171      FMin fmin_func;
1172      return dp_algorithm(fmin_func,numeric_limits<double>::max());
1173    }
1174    
1175    Data
1176    Data::length() const
1177    {
1178    #if defined DOPROF
1179      profData->reduction2++;
1180    #endif
1181      Length len_func;
1182      return dp_algorithm(len_func,0);
1183    }
1184    
1185    Data
1186    Data::trace() const
1187    {
1188    #if defined DOPROF
1189      profData->reduction2++;
1190    #endif
1191      Trace trace_func;
1192      return dp_algorithm(trace_func,0);
1193    }
1194    
1195    Data
1196    Data::transpose(int axis) const
1197  {  {
1198    #if defined DOPROF
1199      profData->reduction2++;
1200    #endif
1201      // not implemented
1202      throw DataException("Error - Data::transpose not implemented yet.");
1203      return Data();
1204    }
1205    
1206    const boost::python::tuple
1207    Data::mindp() const
1208    {
1209      // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1210      // abort (for unknown reasons) if there are openmp directives with it in the
1211      // surrounding function
1212    
1213      int SampleNo;
1214      int DataPointNo;
1215    
1216      calc_mindp(SampleNo,DataPointNo);
1217    
1218      return make_tuple(SampleNo,DataPointNo);
1219    }
1220    
1221    void
1222    Data::calc_mindp(int& SampleNo,
1223                     int& DataPointNo) const
1224    {
1225      int i,j;
1226      int lowi=0,lowj=0;
1227      double min=numeric_limits<double>::max();
1228    
1229      Data temp=minval();
1230    
1231      int numSamples=temp.getNumSamples();
1232      int numDPPSample=temp.getNumDataPointsPerSample();
1233    
1234      double next,local_min;
1235      int local_lowi,local_lowj;
1236    
1237      #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1238      {
1239        local_min=min;
1240        #pragma omp for private(i,j) schedule(static)
1241        for (i=0; i<numSamples; i++) {
1242          for (j=0; j<numDPPSample; j++) {
1243            next=temp.getDataPoint(i,j)();
1244            if (next<local_min) {
1245              local_min=next;
1246              local_lowi=i;
1247              local_lowj=j;
1248            }
1249          }
1250        }
1251        #pragma omp critical
1252        if (local_min<min) {
1253          min=local_min;
1254          lowi=local_lowi;
1255          lowj=local_lowj;
1256        }
1257      }
1258    
1259      SampleNo = lowi;
1260      DataPointNo = lowj;
1261    }
1262    
1263    void
1264    Data::saveDX(std::string fileName) const
1265    {
1266      getDomain().saveDX(fileName,*this);
1267      return;
1268    }
1269    
1270    void
1271    Data::saveVTK(std::string fileName) const
1272    {
1273      getDomain().saveVTK(fileName,*this);
1274      return;
1275    }
1276    
1277    Data&
1278    Data::operator+=(const Data& right)
1279    {
1280    #if defined DOPROF
1281      profData->binary++;
1282    #endif
1283    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1284    return (*this);    return (*this);
1285  }  }
1286    
1287  Data& Data::operator+=(const boost::python::object& right)  Data&
1288    Data::operator+=(const boost::python::object& right)
1289  {  {
1290    #if defined DOPROF
1291      profData->binary++;
1292    #endif
1293    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1294    return (*this);    return (*this);
1295  }  }
1296    
1297  Data& Data::operator-=(const Data& right)  Data&
1298    Data::operator-=(const Data& right)
1299  {  {
1300    #if defined DOPROF
1301      profData->binary++;
1302    #endif
1303    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1304    return (*this);    return (*this);
1305  }  }
1306    
1307  Data& Data::operator-=(const boost::python::object& right)  Data&
1308    Data::operator-=(const boost::python::object& right)
1309  {  {
1310    #if defined DOPROF
1311      profData->binary++;
1312    #endif
1313    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1314    return (*this);    return (*this);
1315  }  }
1316    
1317  Data& Data::operator*=(const Data& right)  Data&
1318    Data::operator*=(const Data& right)
1319  {  {
1320    #if defined DOPROF
1321      profData->binary++;
1322    #endif
1323    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1324    return (*this);    return (*this);
1325  }  }
1326    
1327  Data& Data::operator*=(const boost::python::object& right)  Data&
1328    Data::operator*=(const boost::python::object& right)
1329  {  {
1330    #if defined DOPROF
1331      profData->binary++;
1332    #endif
1333    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1334    return (*this);    return (*this);
1335  }  }
1336    
1337  Data& Data::operator/=(const Data& right)  Data&
1338    Data::operator/=(const Data& right)
1339  {  {
1340    #if defined DOPROF
1341      profData->binary++;
1342    #endif
1343    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1344    return (*this);    return (*this);
1345  }  }
1346    
1347  Data& Data::operator/=(const boost::python::object& right)  Data&
1348    Data::operator/=(const boost::python::object& right)
1349  {  {
1350    #if defined DOPROF
1351      profData->binary++;
1352    #endif
1353    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1354    return (*this);    return (*this);
1355  }  }
1356    
1357  Data Data::powO(const boost::python::object& right) const  Data
1358    Data::powO(const boost::python::object& right) const
1359  {  {
1360    #if defined DOPROF
1361      profData->binary++;
1362    #endif
1363    Data result;    Data result;
1364    result.copy(*this);    result.copy(*this);
1365    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1366    return result;    return result;
1367  }  }
1368    
1369  Data Data::powD(const Data& right) const  Data
1370    Data::powD(const Data& right) const
1371  {  {
1372    #if defined DOPROF
1373      profData->binary++;
1374    #endif
1375    Data result;    Data result;
1376    result.copy(*this);    result.copy(*this);
1377    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 608  Data Data::powD(const Data& right) const Line 1379  Data Data::powD(const Data& right) const
1379  }  }
1380    
1381  //  //
1382  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1383  Data escript::operator+(const Data& left, const Data& right)  Data
1384    escript::operator+(const Data& left, const Data& right)
1385  {  {
1386    Data result;    Data result;
1387    //    //
# Line 620  Data escript::operator+(const Data& left Line 1392  Data escript::operator+(const Data& left
1392  }  }
1393    
1394  //  //
1395  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1396  Data escript::operator-(const Data& left, const Data& right)  Data
1397    escript::operator-(const Data& left, const Data& right)
1398  {  {
1399    Data result;    Data result;
1400    //    //
# Line 632  Data escript::operator-(const Data& left Line 1405  Data escript::operator-(const Data& left
1405  }  }
1406    
1407  //  //
1408  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1409  Data escript::operator*(const Data& left, const Data& right)  Data
1410    escript::operator*(const Data& left, const Data& right)
1411  {  {
1412    Data result;    Data result;
1413    //    //
# Line 644  Data escript::operator*(const Data& left Line 1418  Data escript::operator*(const Data& left
1418  }  }
1419    
1420  //  //
1421  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1422  Data escript::operator/(const Data& left, const Data& right)  Data
1423    escript::operator/(const Data& left, const Data& right)
1424  {  {
1425    Data result;    Data result;
1426    //    //
# Line 656  Data escript::operator/(const Data& left Line 1431  Data escript::operator/(const Data& left
1431  }  }
1432    
1433  //  //
1434  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1435  Data escript::operator+(const Data& left, const boost::python::object& right)  Data
1436    escript::operator+(const Data& left, const boost::python::object& right)
1437  {  {
1438    //    //
1439    // Convert to DataArray format if possible    // Convert to DataArray format if possible
# Line 671  Data escript::operator+(const Data& left Line 1447  Data escript::operator+(const Data& left
1447  }  }
1448    
1449  //  //
1450  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1451  Data escript::operator-(const Data& left, const boost::python::object& right)  Data
1452    escript::operator-(const Data& left, const boost::python::object& right)
1453  {  {
1454    //    //
1455    // Convert to DataArray format if possible    // Convert to DataArray format if possible
# Line 686  Data escript::operator-(const Data& left Line 1463  Data escript::operator-(const Data& left
1463  }  }
1464    
1465  //  //
1466  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1467  Data escript::operator*(const Data& left, const boost::python::object& right)  Data
1468    escript::operator*(const Data& left, const boost::python::object& right)
1469  {  {
1470    //    //
1471    // Convert to DataArray format if possible    // Convert to DataArray format if possible
# Line 701  Data escript::operator*(const Data& left Line 1479  Data escript::operator*(const Data& left
1479  }  }
1480    
1481  //  //
1482  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1483  Data escript::operator/(const Data& left, const boost::python::object& right)  Data
1484    escript::operator/(const Data& left, const boost::python::object& right)
1485  {  {
1486    //    //
1487    // Convert to DataArray format if possible    // Convert to DataArray format if possible
# Line 716  Data escript::operator/(const Data& left Line 1495  Data escript::operator/(const Data& left
1495  }  }
1496    
1497  //  //
1498  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1499  Data escript::operator+(const boost::python::object& left, const Data& right)  Data
1500    escript::operator+(const boost::python::object& left, const Data& right)
1501  {  {
1502    //    //
1503    // Construct the result using the given value and the other parameters    // Construct the result using the given value and the other parameters
# Line 728  Data escript::operator+(const boost::pyt Line 1508  Data escript::operator+(const boost::pyt
1508  }  }
1509    
1510  //  //
1511  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1512  Data escript::operator-(const boost::python::object& left, const Data& right)  Data
1513    escript::operator-(const boost::python::object& left, const Data& right)
1514  {  {
1515    //    //
1516    // Construct the result using the given value and the other parameters    // Construct the result using the given value and the other parameters
# Line 740  Data escript::operator-(const boost::pyt Line 1521  Data escript::operator-(const boost::pyt
1521  }  }
1522    
1523  //  //
1524  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1525  Data escript::operator*(const boost::python::object& left, const Data& right)  Data
1526    escript::operator*(const boost::python::object& left, const Data& right)
1527  {  {
1528    //    //
1529    // Construct the result using the given value and the other parameters    // Construct the result using the given value and the other parameters
# Line 752  Data escript::operator*(const boost::pyt Line 1534  Data escript::operator*(const boost::pyt
1534  }  }
1535    
1536  //  //
1537  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1538  Data escript::operator/(const boost::python::object& left, const Data& right)  Data
1539    escript::operator/(const boost::python::object& left, const Data& right)
1540  {  {
1541    //    //
1542    // Construct the result using the given value and the other parameters    // Construct the result using the given value and the other parameters
# Line 764  Data escript::operator/(const boost::pyt Line 1547  Data escript::operator/(const boost::pyt
1547  }  }
1548    
1549  //  //
1550  // NOTE: It is essential to specify the namepsace this operator belongs to  //bool escript::operator==(const Data& left, const Data& right)
1551  bool escript::operator==(const Data& left, const Data& right)  //{
1552  {  //  /*
1553    /*  //  NB: this operator does very little at this point, and isn't to
1554    NB: this operator does very little at this point, and isn't to  //  be relied on. Requires further implementation.
1555    be relied on. Requires further implementation.  //  */
1556    */  //
1557    //  bool ret;
1558    bool ret;  //
1559    //  if (left.isEmpty()) {
1560    if (left.isEmpty()) {  //    if(!right.isEmpty()) {
1561      if(!right.isEmpty()) {  //      ret = false;
1562        ret = false;  //    } else {
1563      } else {  //      ret = true;
1564        ret = true;  //    }
1565      }  //  }
1566    }  //
1567    //  if (left.isConstant()) {
1568    //    if(!right.isConstant()) {
1569    //      ret = false;
1570    //    } else {
1571    //      ret = true;
1572    //    }
1573    // }
1574    //
1575    //  if (left.isTagged()) {
1576    //   if(!right.isTagged()) {
1577    //      ret = false;
1578    //    } else {
1579    //      ret = true;
1580    //    }
1581    //  }
1582    //
1583    //  if (left.isExpanded()) {
1584    //    if(!right.isExpanded()) {
1585    //      ret = false;
1586    //    } else {
1587    //      ret = true;
1588    //    }
1589    //  }
1590    //
1591    //  return ret;
1592    //}
1593    
1594    if (left.isConstant()) {  Data
1595      if(!right.isConstant()) {  Data::getItem(const boost::python::object& key) const
1596        ret = false;  {
1597      } else {    const DataArrayView& view=getPointDataView();
       ret = true;  
     }  
   }  
1598    
1599    if (left.isTagged()) {    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
     if(!right.isTagged()) {  
       ret = false;  
     } else {  
       ret = true;  
     }  
   }  
1600    
1601    if (left.isExpanded()) {    if (slice_region.size()!=view.getRank()) {
1602      if(!right.isExpanded()) {      throw DataException("Error - slice size does not match Data rank.");
       ret = false;  
     } else {  
       ret = true;  
     }  
1603    }    }
1604    
1605    return ret;    return getSlice(slice_region);
1606  }  }
1607    
1608  Data  Data
1609  Data::getItem(const boost::python::object& key) const  Data::getSlice(const DataArrayView::RegionType& region) const
1610  {  {
1611     const DataArrayView& view=getPointDataView();  #if defined DOPROF
1612     DataArrayView::RegionType slice_region=view.getSliceRegion(key);    profData->slicing++;
1613     if (slice_region.size()!=view.getRank()) {  #endif
1614       throw DataException("Error - slice size does not match Data rank.");    return Data(*this,region);
1615       return Data();  }
1616     }  
1617     //  void
1618     // Create a new Data which is a slice of this one  Data::setItemO(const boost::python::object& key,
1619     return getSlice(slice_region);                 const boost::python::object& value)
1620    {
1621      Data tempData(value,getFunctionSpace());
1622      setItemD(key,tempData);
1623  }  }
1624    
1625  void  void
1626  Data::setItem(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
1627                const Data& value)                 const Data& value)
1628  {  {
1629    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1630    
1631    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1632    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1633      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
1634    }    }
1635    typeMatch(value);    if (getFunctionSpace()!=value.getFunctionSpace()) {
1636    setSlice(value,slice_region);       setSlice(Data(value,getFunctionSpace()),slice_region);
1637      } else {
1638         setSlice(value,slice_region);
1639      }
1640  }  }
1641    
1642  void  void
1643  Data::typeMatch(const Data& right)  Data::setSlice(const Data& value,
1644                   const DataArrayView::RegionType& region)
1645    {
1646    #if defined DOPROF
1647      profData->slicing++;
1648    #endif
1649      Data tempValue(value);
1650      typeMatchLeft(tempValue);
1651      typeMatchRight(tempValue);
1652      m_data->setSlice(tempValue.m_data.get(),region);
1653    }
1654    
1655    void
1656    Data::typeMatchLeft(Data& right) const
1657    {
1658      if (isExpanded()){
1659        right.expand();
1660      } else if (isTagged()) {
1661        if (right.isConstant()) {
1662          right.tag();
1663        }
1664      }
1665    }
1666    
1667    void
1668    Data::typeMatchRight(const Data& right)
1669  {  {
   //  
   // match the type of this to the RHS  
1670    if (isTagged()) {    if (isTagged()) {
1671      if (right.isExpanded()) {      if (right.isExpanded()) {
       //  
       // if the right hand side is expanded so must this  
1672        expand();        expand();
1673      }      }
1674    } else if (isConstant()) {    } else if (isConstant()) {
1675      if (right.isExpanded()) {      if (right.isExpanded()) {
       //  
       // if the right hand side is expanded so must this  
1676        expand();        expand();
1677      } else if (right.isTagged()) {      } else if (right.isTagged()) {
       //  
       // if the right hand side is tagged so must this  
1678        tag();        tag();
1679      }      }
1680    }    }
# Line 881  Data::setTaggedValue(int tagKey, Line 1701  Data::setTaggedValue(int tagKey,
1701    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
1702  }  }
1703    
 /*  
 Note: this version removed for now. Not needed, and breaks escript.cpp  
1704  void  void
1705  Data::setTaggedValue(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
1706                       const DataArrayView& value)                              const DataArrayView& value)
1707  {  {
1708    //    //
1709    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
# Line 899  Data::setTaggedValue(int tagKey, Line 1717  Data::setTaggedValue(int tagKey,
1717    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
1718    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
1719  }  }
1720  */  
1721    int
1722    Data::getTagNumber(int dpno)
1723    {
1724      return m_data->getTagNumber(dpno);
1725    }
1726    
1727    void
1728    Data::setRefValue(int ref,
1729                      const boost::python::numeric::array& value)
1730    {
1731      //
1732      // Construct DataArray from boost::python::object input value
1733      DataArray valueDataArray(value);
1734    
1735      //
1736      // Call DataAbstract::setRefValue
1737      m_data->setRefValue(ref,valueDataArray);
1738    }
1739    
1740    void
1741    Data::getRefValue(int ref,
1742                      boost::python::numeric::array& value)
1743    {
1744      //
1745      // Construct DataArray for boost::python::object return value
1746      DataArray valueDataArray(value);
1747    
1748      //
1749      // Load DataArray with values from data-points specified by ref
1750      m_data->getRefValue(ref,valueDataArray);
1751    
1752      //
1753      // Load values from valueDataArray into return numarray
1754    
1755      // extract the shape of the numarray
1756      int rank = value.getrank();
1757      DataArrayView::ShapeType shape;
1758      for (int i=0; i < rank; i++) {
1759        shape.push_back(extract<int>(value.getshape()[i]));
1760      }
1761    
1762      // and load the numarray with the data from the DataArray
1763      DataArrayView valueView = valueDataArray.getView();
1764    
1765      if (rank==0) {
1766          boost::python::numeric::array temp_numArray(valueView());
1767          value = temp_numArray;
1768      }
1769      if (rank==1) {
1770        for (int i=0; i < shape[0]; i++) {
1771          value[i] = valueView(i);
1772        }
1773      }
1774      if (rank==2) {
1775        for (int i=0; i < shape[0]; i++) {
1776          for (int j=0; j < shape[1]; j++) {
1777            value[i][j] = valueView(i,j);
1778          }
1779        }
1780      }
1781      if (rank==3) {
1782        for (int i=0; i < shape[0]; i++) {
1783          for (int j=0; j < shape[1]; j++) {
1784            for (int k=0; k < shape[2]; k++) {
1785              value[i][j][k] = valueView(i,j,k);
1786            }
1787          }
1788        }
1789      }
1790      if (rank==4) {
1791        for (int i=0; i < shape[0]; i++) {
1792          for (int j=0; j < shape[1]; j++) {
1793            for (int k=0; k < shape[2]; k++) {
1794              for (int l=0; l < shape[3]; l++) {
1795                value[i][j][k][l] = valueView(i,j,k,l);
1796              }
1797            }
1798          }
1799        }
1800      }
1801    
1802    }
1803    
1804    void
1805    Data::archiveData(const std::string fileName)
1806    {
1807      cout << "Archiving Data object to: " << fileName << endl;
1808    
1809      //
1810      // Determine type of this Data object
1811      int dataType = -1;
1812    
1813      if (isEmpty()) {
1814        dataType = 0;
1815        cout << "\tdataType: DataEmpty" << endl;
1816      }
1817      if (isConstant()) {
1818        dataType = 1;
1819        cout << "\tdataType: DataConstant" << endl;
1820      }
1821      if (isTagged()) {
1822        dataType = 2;
1823        cout << "\tdataType: DataTagged" << endl;
1824      }
1825      if (isExpanded()) {
1826        dataType = 3;
1827        cout << "\tdataType: DataExpanded" << endl;
1828      }
1829    
1830      if (dataType == -1) {
1831        throw DataException("archiveData Error: undefined dataType");
1832      }
1833    
1834      //
1835      // Collect data items common to all Data types
1836      int noSamples = getNumSamples();
1837      int noDPPSample = getNumDataPointsPerSample();
1838      int functionSpaceType = getFunctionSpace().getTypeCode();
1839      int dataPointRank = getDataPointRank();
1840      int dataPointSize = getDataPointSize();
1841      int dataLength = getLength();
1842      DataArrayView::ShapeType dataPointShape = getDataPointShape();
1843      int referenceNumbers[noSamples];
1844      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1845        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
1846      }
1847      int tagNumbers[noSamples];
1848      if (isTagged()) {
1849        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1850          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
1851        }
1852      }
1853    
1854      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
1855      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
1856      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
1857    
1858      //
1859      // Flatten Shape to an array of integers suitable for writing to file
1860      int flatShape[4] = {0,0,0,0};
1861      cout << "\tshape: < ";
1862      for (int dim=0; dim<dataPointRank; dim++) {
1863        flatShape[dim] = dataPointShape[dim];
1864        cout << dataPointShape[dim] << " ";
1865      }
1866      cout << ">" << endl;
1867    
1868      //
1869      // Open archive file
1870      ofstream archiveFile;
1871      archiveFile.open(fileName.data(), ios::out);
1872    
1873      if (!archiveFile.good()) {
1874        throw DataException("archiveData Error: problem opening archive file");
1875      }
1876    
1877      //
1878      // Write common data items to archive file
1879      archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1880      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1881      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1882      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1883      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1884      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1885      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
1886      for (int dim = 0; dim < 4; dim++) {
1887        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1888      }
1889      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1890        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
1891      }
1892      if (isTagged()) {
1893        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1894          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
1895        }
1896      }
1897    
1898      if (!archiveFile.good()) {
1899        throw DataException("archiveData Error: problem writing to archive file");
1900      }
1901    
1902      //
1903      // Archive underlying data values for each Data type
1904      int noValues;
1905      switch (dataType) {
1906        case 0:
1907          // DataEmpty
1908          noValues = 0;
1909          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1910          cout << "\tnoValues: " << noValues << endl;
1911          break;
1912        case 1:
1913          // DataConstant
1914          noValues = m_data->getLength();
1915          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1916          cout << "\tnoValues: " << noValues << endl;
1917          if (m_data->archiveData(archiveFile,noValues)) {
1918            throw DataException("archiveData Error: problem writing data to archive file");
1919          }
1920          break;
1921        case 2:
1922          // DataTagged
1923          noValues = m_data->getLength();
1924          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1925          cout << "\tnoValues: " << noValues << endl;
1926          if (m_data->archiveData(archiveFile,noValues)) {
1927            throw DataException("archiveData Error: problem writing data to archive file");
1928          }
1929          break;
1930        case 3:
1931          // DataExpanded
1932          noValues = m_data->getLength();
1933          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1934          cout << "\tnoValues: " << noValues << endl;
1935          if (m_data->archiveData(archiveFile,noValues)) {
1936            throw DataException("archiveData Error: problem writing data to archive file");
1937          }
1938          break;
1939      }
1940    
1941      if (!archiveFile.good()) {
1942        throw DataException("archiveData Error: problem writing data to archive file");
1943      }
1944    
1945      //
1946      // Close archive file
1947      archiveFile.close();
1948    
1949      if (!archiveFile.good()) {
1950        throw DataException("archiveData Error: problem closing archive file");
1951      }
1952    
1953    }
1954    
1955    void
1956    Data::extractData(const std::string fileName,
1957                      const FunctionSpace& fspace)
1958    {
1959      //
1960      // Can only extract Data to an object which is initially DataEmpty
1961      if (!isEmpty()) {
1962        throw DataException("extractData Error: can only extract to DataEmpty object");
1963      }
1964    
1965      cout << "Extracting Data object from: " << fileName << endl;
1966    
1967      int dataType;
1968      int noSamples;
1969      int noDPPSample;
1970      int functionSpaceType;
1971      int dataPointRank;
1972      int dataPointSize;
1973      int dataLength;
1974      DataArrayView::ShapeType dataPointShape;
1975      int flatShape[4];
1976    
1977      //
1978      // Open the archive file
1979      ifstream archiveFile;
1980      archiveFile.open(fileName.data(), ios::in);
1981    
1982      if (!archiveFile.good()) {
1983        throw DataException("extractData Error: problem opening archive file");
1984      }
1985    
1986      //
1987      // Read common data items from archive file
1988      archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1989      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1990      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
1991      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
1992      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
1993      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
1994      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
1995      for (int dim = 0; dim < 4; dim++) {
1996        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
1997        if (flatShape[dim]>0) {
1998          dataPointShape.push_back(flatShape[dim]);
1999        }
2000      }
2001      int referenceNumbers[noSamples];
2002      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2003        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2004      }
2005      int tagNumbers[noSamples];
2006      if (dataType==2) {
2007        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2008          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2009        }
2010      }
2011    
2012      if (!archiveFile.good()) {
2013        throw DataException("extractData Error: problem reading from archive file");
2014      }
2015    
2016      //
2017      // Verify the values just read from the archive file
2018      switch (dataType) {
2019        case 0:
2020          cout << "\tdataType: DataEmpty" << endl;
2021          break;
2022        case 1:
2023          cout << "\tdataType: DataConstant" << endl;
2024          break;
2025        case 2:
2026          cout << "\tdataType: DataTagged" << endl;
2027          break;
2028        case 3:
2029          cout << "\tdataType: DataExpanded" << endl;
2030          break;
2031        default:
2032          throw DataException("extractData Error: undefined dataType read from archive file");
2033          break;
2034      }
2035    
2036      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2037      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2038      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2039      cout << "\tshape: < ";
2040      for (int dim = 0; dim < dataPointRank; dim++) {
2041        cout << dataPointShape[dim] << " ";
2042      }
2043      cout << ">" << endl;
2044    
2045      //
2046      // Verify that supplied FunctionSpace object is compatible with this Data object.
2047      if ( (fspace.getTypeCode()!=functionSpaceType) ||
2048           (fspace.getNumSamples()!=noSamples) ||
2049           (fspace.getNumDPPSample()!=noDPPSample)
2050         ) {
2051        throw DataException("extractData Error: incompatible FunctionSpace");
2052      }
2053      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2054        if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2055          throw DataException("extractData Error: incompatible FunctionSpace");
2056        }
2057      }
2058      if (dataType==2) {
2059        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2060          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2061            throw DataException("extractData Error: incompatible FunctionSpace");
2062          }
2063        }
2064      }
2065    
2066      //
2067      // Construct a DataVector to hold underlying data values
2068      DataVector dataVec(dataLength);
2069    
2070      //
2071      // Load this DataVector with the appropriate values
2072      int noValues;
2073      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2074      cout << "\tnoValues: " << noValues << endl;
2075      switch (dataType) {
2076        case 0:
2077          // DataEmpty
2078          if (noValues != 0) {
2079            throw DataException("extractData Error: problem reading data from archive file");
2080          }
2081          break;
2082        case 1:
2083          // DataConstant
2084          if (dataVec.extractData(archiveFile,noValues)) {
2085            throw DataException("extractData Error: problem reading data from archive file");
2086          }
2087          break;
2088        case 2:
2089          // DataTagged
2090          if (dataVec.extractData(archiveFile,noValues)) {
2091            throw DataException("extractData Error: problem reading data from archive file");
2092          }
2093          break;
2094        case 3:
2095          // DataExpanded
2096          if (dataVec.extractData(archiveFile,noValues)) {
2097            throw DataException("extractData Error: problem reading data from archive file");
2098          }
2099          break;
2100      }
2101    
2102      if (!archiveFile.good()) {
2103        throw DataException("extractData Error: problem reading from archive file");
2104      }
2105    
2106      //
2107      // Close archive file
2108      archiveFile.close();
2109    
2110      if (!archiveFile.good()) {
2111        throw DataException("extractData Error: problem closing archive file");
2112      }
2113    
2114      //
2115      // Construct an appropriate Data object
2116      DataAbstract* tempData;
2117      switch (dataType) {
2118        case 0:
2119          // DataEmpty
2120          tempData=new DataEmpty();
2121          break;
2122        case 1:
2123          // DataConstant
2124          tempData=new DataConstant(fspace,dataPointShape,dataVec);
2125          break;
2126        case 2:
2127          // DataTagged
2128          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2129          break;
2130        case 3:
2131          // DataExpanded
2132          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2133          break;
2134      }
2135      shared_ptr<DataAbstract> temp_data(tempData);
2136      m_data=temp_data;
2137    }
2138    
2139  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
2140  {  {

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

  ViewVC Help
Powered by ViewVC 1.1.26