/[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

trunk/esys2/escript/src/Data/Data.cpp revision 119 by jgs, Tue Apr 12 04:45:05 2005 UTC trunk/escript/src/Data.cpp revision 790 by bcumming, Wed Jul 26 23:12:34 2006 UTC
# Line 1  Line 1 
1  // $Id$  // $Id$
 /*=============================================================================  
2    
3   ******************************************************************************  /*
4   *                                                                            *   ************************************************************
5   *       COPYRIGHT ACcESS 2004 -  All Rights Reserved                         *   *          Copyright 2006 by ACcESS MNRF                   *
6   *                                                                            *   *                                                          *
7   * This software is the property of ACcESS.  No part of this code             *   *              http://www.access.edu.au                    *
8   * may be copied in any form or by any means without the expressed written    *   *       Primary Business: Queensland, Australia            *
9   * consent of ACcESS.  Copying, use or modification of this software          *   *  Licensed under the Open Software License version 3.0    *
10   * by any unauthorised person is illegal unless that                          *   *     http://www.opensource.org/licenses/osl-3.0.php       *
11   * person has a software license agreement with ACcESS.                       *   *                                                          *
12   *                                                                            *   ************************************************************
13   ******************************************************************************  */
14    #include "Data.h"
 ******************************************************************************/  
15    
16  #include "escript/Data/Data.h"  #include "DataExpanded.h"
17    #include "DataConstant.h"
18    #include "DataTagged.h"
19    #include "DataEmpty.h"
20    #include "DataArray.h"
21    #include "DataArrayView.h"
22    #include "DataProf.h"
23    #include "FunctionSpaceFactory.h"
24    #include "AbstractContinuousDomain.h"
25    #include "UnaryFuncs.h"
26    
 #include <iostream>  
27  #include <fstream>  #include <fstream>
28  #include <algorithm>  #include <algorithm>
29  #include <vector>  #include <vector>
 #include <exception>  
30  #include <functional>  #include <functional>
 #include <math.h>  
31    
32  #include <boost/python/str.hpp>  #include <boost/python/dict.hpp>
33  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
34  #include <boost/python/long.hpp>  #include <boost/python/long.hpp>
35    
 #include "escript/Data/DataException.h"  
 #include "escript/Data/DataExpanded.h"  
 #include "escript/Data/DataConstant.h"  
 #include "escript/Data/DataTagged.h"  
 #include "escript/Data/DataEmpty.h"  
 #include "escript/Data/DataArray.h"  
 #include "escript/Data/DataAlgorithm.h"  
 #include "escript/Data/FunctionSpaceFactory.h"  
 #include "escript/Data/AbstractContinuousDomain.h"  
 #include "escript/Data/UnaryFuncs.h"  
   
36  using namespace std;  using namespace std;
37  using namespace boost::python;  using namespace boost::python;
38  using namespace boost;  using namespace boost;
39  using namespace escript;  using namespace escript;
40    
41    #if defined DOPROF
42    //
43    // global table of profiling data for all Data objects
44    DataProf dataProfTable;
45    #endif
46    
47  Data::Data()  Data::Data()
48  {  {
49    //    //
# Line 52  Data::Data() Line 51  Data::Data()
51    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
52    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
53    m_data=temp_data;    m_data=temp_data;
54      m_protected=false;
55    #if defined DOPROF
56      // create entry in global profiling table for this object
57      profData = dataProfTable.newData();
58    #endif
59  }  }
60    
61  Data::Data(double value,  Data::Data(double value,
# Line 65  Data::Data(double value, Line 69  Data::Data(double value,
69    }    }
70    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
71    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
72      m_protected=false;
73    #if defined DOPROF
74      // create entry in global profiling table for this object
75      profData = dataProfTable.newData();
76    #endif
77  }  }
78    
79  Data::Data(double value,  Data::Data(double value,
# Line 75  Data::Data(double value, Line 84  Data::Data(double value,
84    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
85    pair<int,int> dataShape=what.getDataShape();    pair<int,int> dataShape=what.getDataShape();
86    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
87      m_protected=false;
88    #if defined DOPROF
89      // create entry in global profiling table for this object
90      profData = dataProfTable.newData();
91    #endif
92  }  }
93    
94  Data::Data(const Data& inData)  Data::Data(const Data& inData)
95  {  {
96    m_data=inData.m_data;    m_data=inData.m_data;
97      m_protected=inData.isProtected();
98    #if defined DOPROF
99      // create entry in global profiling table for this object
100      profData = dataProfTable.newData();
101    #endif
102  }  }
103    
104  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 90  Data::Data(const Data& inData, Line 109  Data::Data(const Data& inData,
109    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
110    shared_ptr<DataAbstract> temp_data(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
111    m_data=temp_data;    m_data=temp_data;
112      m_protected=false;
113    #if defined DOPROF
114      // create entry in global profiling table for this object
115      profData = dataProfTable.newData();
116    #endif
117  }  }
118    
119  Data::Data(const Data& inData,  Data::Data(const Data& inData,
120             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
121  {  {
122    #if defined DOPROF
123      // create entry in global profiling table for this object
124      profData = dataProfTable.newData();
125    #endif
126    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
127      m_data=inData.m_data;      m_data=inData.m_data;
128    } else {    } else {
129        #if defined DOPROF
130        profData->interpolate++;
131        #endif
132      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
133      // Note for Lutz, Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
134      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
135      // be able to create an instance of AbstractDomain but that was done      // be able to create an instance of AbstractDomain but that was done
136      // as a boost python work around which may no longer be required.      // as a boost:python work around which may no longer be required.
137      const AbstractDomain& inDataDomain=inData.getDomain();      const AbstractDomain& inDataDomain=inData.getDomain();
138      if  (inDataDomain==functionspace.getDomain()) {      if  (inDataDomain==functionspace.getDomain()) {
139        inDataDomain.interpolateOnDomain(tmp,inData);        inDataDomain.interpolateOnDomain(tmp,inData);
# Line 111  Data::Data(const Data& inData, Line 142  Data::Data(const Data& inData,
142      }      }
143      m_data=tmp.m_data;      m_data=tmp.m_data;
144    }    }
145      m_protected=false;
146  }  }
147    
148  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(const DataTagged::TagListType& tagKeys,
# Line 122  Data::Data(const DataTagged::TagListType Line 154  Data::Data(const DataTagged::TagListType
154    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
155    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
156    m_data=temp_data;    m_data=temp_data;
157      m_protected=false;
158    if (expanded) {    if (expanded) {
159      expand();      expand();
160    }    }
161    #if defined DOPROF
162      // create entry in global profiling table for this object
163      profData = dataProfTable.newData();
164    #endif
165  }  }
166    
167  Data::Data(const numeric::array& value,  Data::Data(const numeric::array& value,
# Line 132  Data::Data(const numeric::array& value, Line 169  Data::Data(const numeric::array& value,
169             bool expanded)             bool expanded)
170  {  {
171    initialise(value,what,expanded);    initialise(value,what,expanded);
172      m_protected=false;
173    #if defined DOPROF
174      // create entry in global profiling table for this object
175      profData = dataProfTable.newData();
176    #endif
177  }  }
178    
179  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
# Line 139  Data::Data(const DataArrayView& value, Line 181  Data::Data(const DataArrayView& value,
181             bool expanded)             bool expanded)
182  {  {
183    initialise(value,what,expanded);    initialise(value,what,expanded);
184      m_protected=false;
185    #if defined DOPROF
186      // create entry in global profiling table for this object
187      profData = dataProfTable.newData();
188    #endif
189  }  }
190    
191  Data::Data(const object& value,  Data::Data(const object& value,
# Line 147  Data::Data(const object& value, Line 194  Data::Data(const object& value,
194  {  {
195    numeric::array asNumArray(value);    numeric::array asNumArray(value);
196    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
197      m_protected=false;
198    #if defined DOPROF
199      // create entry in global profiling table for this object
200      profData = dataProfTable.newData();
201    #endif
202  }  }
203    
204  Data::Data(const object& value,  Data::Data(const object& value,
# Line 167  Data::Data(const object& value, Line 219  Data::Data(const object& value,
219      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
220      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
221    }    }
222      m_protected=false;
223    #if defined DOPROF
224      // create entry in global profiling table for this object
225      profData = dataProfTable.newData();
226    #endif
227    }
228    
229    Data::~Data()
230    {
231    
232  }  }
233    
234  escriptDataC  escriptDataC
# Line 185  Data::getDataC() const Line 247  Data::getDataC() const
247    return temp;    return temp;
248  }  }
249    
250  tuple  const boost::python::tuple
251  Data::getShapeTuple() const  Data::getShapeTuple() const
252  {  {
253    const DataArrayView::ShapeType& shape=getDataPointShape();    const DataArrayView::ShapeType& shape=getDataPointShape();
# Line 288  Data::isTagged() const Line 350  Data::isTagged() const
350    return (temp!=0);    return (temp!=0);
351  }  }
352    
353    /* TODO */
354    /* global reduction -- the local data being empty does not imply that it is empty on other processers*/
355  bool  bool
356  Data::isEmpty() const  Data::isEmpty() const
357  {  {
# Line 303  Data::isConstant() const Line 367  Data::isConstant() const
367  }  }
368    
369  void  void
370    Data::setProtection()
371    {
372       m_protected=true;
373    }
374    
375    bool
376    Data::isProtected() const
377    {
378       return m_protected;
379    }
380    
381    
382    
383    void
384  Data::expand()  Data::expand()
385  {  {
386    if (isConstant()) {    if (isConstant()) {
# Line 353  Data::reshapeDataPoint(const DataArrayVi Line 431  Data::reshapeDataPoint(const DataArrayVi
431  Data  Data
432  Data::wherePositive() const  Data::wherePositive() const
433  {  {
434    #if defined DOPROF
435      profData->where++;
436    #endif
437    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
438  }  }
439    
440  Data  Data
441  Data::whereNegative() const  Data::whereNegative() const
442  {  {
443    #if defined DOPROF
444      profData->where++;
445    #endif
446    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
447  }  }
448    
449  Data  Data
450  Data::whereNonNegative() const  Data::whereNonNegative() const
451  {  {
452    #if defined DOPROF
453      profData->where++;
454    #endif
455    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
456  }  }
457    
458  Data  Data
459  Data::whereNonPositive() const  Data::whereNonPositive() const
460  {  {
461    #if defined DOPROF
462      profData->where++;
463    #endif
464    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
465  }  }
466    
467  Data  Data
468  Data::whereZero() const  Data::whereZero(double tol) const
469  {  {
470    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));  #if defined DOPROF
471      profData->where++;
472    #endif
473      Data dataAbs=abs();
474      return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));
475  }  }
476    
477  Data  Data
478  Data::whereNonZero() const  Data::whereNonZero(double tol) const
479  {  {
480    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));  #if defined DOPROF
481      profData->where++;
482    #endif
483      Data dataAbs=abs();
484      return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));
485  }  }
486    
487  Data  Data
488  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
489  {  {
490    #if defined DOPROF
491      profData->interpolate++;
492    #endif
493    return Data(*this,functionspace);    return Data(*this,functionspace);
494  }  }
495    
# Line 410  Data::probeInterpolation(const FunctionS Line 511  Data::probeInterpolation(const FunctionS
511  Data  Data
512  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
513  {  {
514    #if defined DOPROF
515      profData->grad++;
516    #endif
517    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
518      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
519    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
# Line 443  Data::getDataPointShape() const Line 547  Data::getDataPointShape() const
547    return getPointDataView().getShape();    return getPointDataView().getShape();
548  }  }
549    
550    void
551    Data::fillFromNumArray(const boost::python::numeric::array num_array)
552    {
553      if (isProtected()) {
554            throw DataException("Error - attempt to update protected Data object.");
555      }
556      //
557      // check rank
558      if (num_array.getrank()<getDataPointRank())
559          throw DataException("Rank of numarray does not match Data object rank");
560    
561      //
562      // check shape of num_array
563      for (int i=0; i<getDataPointRank(); i++) {
564        if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])
565           throw DataException("Shape of numarray does not match Data object rank");
566      }
567    
568      //
569      // make sure data is expanded:
570      if (!isExpanded()) {
571        expand();
572      }
573    
574      //
575      // and copy over
576      m_data->copyAll(num_array);
577    }
578    
579    const
580  boost::python::numeric::array  boost::python::numeric::array
581  Data::convertToNumArray()  Data::convertToNumArray()
582  {  {
583    //    //
584    // determine the current number of data points    // determine the total number of data points
585    int numSamples = getNumSamples();    int numSamples = getNumSamples();
586    int numDataPointsPerSample = getNumDataPointsPerSample();    int numDataPointsPerSample = getNumDataPointsPerSample();
587    int numDataPoints = numSamples * numDataPointsPerSample;    int numDataPoints = numSamples * numDataPointsPerSample;
# Line 499  Data::convertToNumArray() Line 633  Data::convertToNumArray()
633    int dataPoint = 0;    int dataPoint = 0;
634    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
635      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
   
636        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
   
637        if (dataPointRank==0) {        if (dataPointRank==0) {
638          numArray[dataPoint]=dataPointView();          numArray[dataPoint]=dataPointView();
639        }        }
   
640        if (dataPointRank==1) {        if (dataPointRank==1) {
641          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
642            numArray[dataPoint][i]=dataPointView(i);            numArray[dataPoint][i]=dataPointView(i);
643          }          }
644        }        }
   
645        if (dataPointRank==2) {        if (dataPointRank==2) {
646          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
647            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
# Line 519  Data::convertToNumArray() Line 649  Data::convertToNumArray()
649            }            }
650          }          }
651        }        }
   
652        if (dataPointRank==3) {        if (dataPointRank==3) {
653          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
654            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
# Line 529  Data::convertToNumArray() Line 658  Data::convertToNumArray()
658            }            }
659          }          }
660        }        }
   
661        if (dataPointRank==4) {        if (dataPointRank==4) {
662          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
663            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
# Line 541  Data::convertToNumArray() Line 669  Data::convertToNumArray()
669            }            }
670          }          }
671        }        }
   
672        dataPoint++;        dataPoint++;
673        }
674      }
675    
676      //
677      // return the loaded array
678      return numArray;
679    }
680    
681    const
682    boost::python::numeric::array
683    Data::convertToNumArrayFromSampleNo(int sampleNo)
684    {
685      //
686      // Check a valid sample number has been supplied
687      if (sampleNo >= getNumSamples()) {
688        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
689      }
690    
691      //
692      // determine the number of data points per sample
693      int numDataPointsPerSample = getNumDataPointsPerSample();
694    
695      //
696      // determine the rank and shape of each data point
697      int dataPointRank = getDataPointRank();
698      DataArrayView::ShapeType dataPointShape = getDataPointShape();
699    
700      //
701      // create the numeric array to be returned
702      boost::python::numeric::array numArray(0.0);
703    
704      //
705      // the rank of the returned numeric array will be the rank of
706      // the data points, plus one. Where the rank of the array is n,
707      // the last n-1 dimensions will be equal to the shape of the
708      // data points, whilst the first dimension will be equal to the
709      // total number of data points. Thus the array will consist of
710      // a serial vector of the data points.
711      int arrayRank = dataPointRank + 1;
712      DataArrayView::ShapeType arrayShape;
713      arrayShape.push_back(numDataPointsPerSample);
714      for (int d=0; d<dataPointRank; d++) {
715         arrayShape.push_back(dataPointShape[d]);
716      }
717    
718      //
719      // resize the numeric array to the shape just calculated
720      if (arrayRank==1) {
721        numArray.resize(arrayShape[0]);
722      }
723      if (arrayRank==2) {
724        numArray.resize(arrayShape[0],arrayShape[1]);
725      }
726      if (arrayRank==3) {
727        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
728      }
729      if (arrayRank==4) {
730        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
731      }
732      if (arrayRank==5) {
733        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
734      }
735    
736      //
737      // loop through each data point in turn, loading the values for that data point
738      // into the numeric array.
739      for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {
740        DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);
741        if (dataPointRank==0) {
742          numArray[dataPoint]=dataPointView();
743        }
744        if (dataPointRank==1) {
745          for (int i=0; i<dataPointShape[0]; i++) {
746            numArray[dataPoint][i]=dataPointView(i);
747          }
748        }
749        if (dataPointRank==2) {
750          for (int i=0; i<dataPointShape[0]; i++) {
751            for (int j=0; j<dataPointShape[1]; j++) {
752              numArray[dataPoint][i][j] = dataPointView(i,j);
753            }
754          }
755        }
756        if (dataPointRank==3) {
757          for (int i=0; i<dataPointShape[0]; i++) {
758            for (int j=0; j<dataPointShape[1]; j++) {
759              for (int k=0; k<dataPointShape[2]; k++) {
760                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
761              }
762            }
763          }
764        }
765        if (dataPointRank==4) {
766          for (int i=0; i<dataPointShape[0]; i++) {
767            for (int j=0; j<dataPointShape[1]; j++) {
768              for (int k=0; k<dataPointShape[2]; k++) {
769                for (int l=0; l<dataPointShape[3]; l++) {
770                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
771                }
772              }
773            }
774          }
775      }      }
776    }    }
777    
# Line 552  Data::convertToNumArray() Line 780  Data::convertToNumArray()
780    return numArray;    return numArray;
781  }  }
782    
783    const
784    boost::python::numeric::array
785    Data::convertToNumArrayFromDPNo(int procNo,
786                                    int sampleNo,
787                                                                    int dataPointNo)
788    
789    {
790        size_t length=0;
791        int i, j, k, l, pos;
792    
793      //
794      // Check a valid sample number has been supplied
795      if (sampleNo >= getNumSamples()) {
796        throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
797      }
798    
799      //
800      // Check a valid data point number has been supplied
801      if (dataPointNo >= getNumDataPointsPerSample()) {
802        throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");
803      }
804    
805      //
806      // determine the rank and shape of each data point
807      int dataPointRank = getDataPointRank();
808      DataArrayView::ShapeType dataPointShape = getDataPointShape();
809    
810      //
811      // create the numeric array to be returned
812      boost::python::numeric::array numArray(0.0);
813    
814      //
815      // the shape of the returned numeric array will be the same
816      // as that of the data point
817      int arrayRank = dataPointRank;
818      DataArrayView::ShapeType arrayShape = dataPointShape;
819    
820      //
821      // resize the numeric array to the shape just calculated
822      if (arrayRank==0) {
823        numArray.resize(1);
824      }
825      if (arrayRank==1) {
826        numArray.resize(arrayShape[0]);
827      }
828      if (arrayRank==2) {
829        numArray.resize(arrayShape[0],arrayShape[1]);
830      }
831      if (arrayRank==3) {
832        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
833      }
834      if (arrayRank==4) {
835        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
836      }
837    
838        // added for the MPI communication
839        length=1;
840        for( i=0; i<arrayRank; i++ )
841            length *= arrayShape[i];
842        double *tmpData = new double[length];
843    
844      //
845      // load the values for the data point into the numeric array.
846    
847        // updated for the MPI case
848        if( get_MPIRank()==procNo ){
849            // create a view of the data if it is stored locally
850            DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
851            
852            // pack the data from the view into tmpData for MPI communication
853            pos=0;
854            switch( dataPointRank ){
855                case 0 :
856                    tmpData[0] = dataPointView();
857                    break;
858                case 1 :        
859                    for( i=0; i<dataPointShape[0]; i++ )
860                        tmpData[i]=dataPointView(i);
861                    break;
862                case 2 :        
863                    for( i=0; i<dataPointShape[0]; i++ )
864                        for( j=0; j<dataPointShape[1]; j++, pos++ )
865                            tmpData[pos]=dataPointView(i,j);
866                    break;
867                case 3 :        
868                    for( i=0; i<dataPointShape[0]; i++ )
869                        for( j=0; j<dataPointShape[1]; j++ )
870                            for( k=0; k<dataPointShape[2]; k++, pos++ )
871                                tmpData[pos]=dataPointView(i,j,k);
872                    break;
873                case 4 :
874                    for( i=0; i<dataPointShape[0]; i++ )
875                        for( j=0; j<dataPointShape[1]; j++ )
876                            for( k=0; k<dataPointShape[2]; k++ )
877                                for( l=0; l<dataPointShape[3]; l++, pos++ )
878                                    tmpData[pos]=dataPointView(i,j,k,l);
879                    break;
880            }
881        }
882    #ifdef PASO_MPI
883            // broadcast the data to all other processes
884            MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
885    #endif
886    
887        // unpack the data
888        switch( dataPointRank ){
889            case 0 :
890                numArray[i]=tmpData[0];
891                break;
892            case 1 :        
893                for( i=0; i<dataPointShape[0]; i++ )
894                    numArray[i]=tmpData[i];
895                break;
896            case 2 :        
897                for( i=0; i<dataPointShape[0]; i++ )
898                    for( j=0; j<dataPointShape[1]; j++ )
899                        tmpData[i+j*dataPointShape[0]];
900                break;
901            case 3 :        
902                for( i=0; i<dataPointShape[0]; i++ )
903                    for( j=0; j<dataPointShape[1]; j++ )
904                        for( k=0; k<dataPointShape[2]; k++ )
905                            tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
906                break;
907            case 4 :
908                for( i=0; i<dataPointShape[0]; i++ )
909                    for( j=0; j<dataPointShape[1]; j++ )
910                        for( k=0; k<dataPointShape[2]; k++ )
911                            for( l=0; l<dataPointShape[3]; l++ )
912                                tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
913                break;
914        }
915    
916        delete [] tmpData;  
917    /*
918      if (dataPointRank==0) {
919        numArray[0]=dataPointView();
920      }
921      if (dataPointRank==1) {
922        for (int i=0; i<dataPointShape[0]; i++) {
923          numArray[i]=dataPointView(i);
924        }
925      }
926      if (dataPointRank==2) {
927        for (int i=0; i<dataPointShape[0]; i++) {
928          for (int j=0; j<dataPointShape[1]; j++) {
929            numArray[i][j] = dataPointView(i,j);
930          }
931        }
932      }
933      if (dataPointRank==3) {
934        for (int i=0; i<dataPointShape[0]; i++) {
935          for (int j=0; j<dataPointShape[1]; j++) {
936            for (int k=0; k<dataPointShape[2]; k++) {
937              numArray[i][j][k]=dataPointView(i,j,k);
938            }
939          }
940        }
941      }
942      if (dataPointRank==4) {
943        for (int i=0; i<dataPointShape[0]; i++) {
944          for (int j=0; j<dataPointShape[1]; j++) {
945            for (int k=0; k<dataPointShape[2]; k++) {
946              for (int l=0; l<dataPointShape[3]; l++) {
947                numArray[i][j][k][l]=dataPointView(i,j,k,l);
948              }
949            }
950          }
951        }
952      }
953    */
954    
955      //
956      // return the loaded array
957      return numArray;
958    }
959    
960  boost::python::numeric::array  boost::python::numeric::array
961  Data::integrate() const  Data::integrate() const
962  {  {
# Line 559  Data::integrate() const Line 964  Data::integrate() const
964    int rank = getDataPointRank();    int rank = getDataPointRank();
965    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
966    
967    #if defined DOPROF
968      profData->integrate++;
969    #endif
970    
971    //    //
972    // calculate the integral values    // calculate the integral values
973    vector<double> integrals(getDataPointSize());    vector<double> integrals(getDataPointSize());
# Line 581  Data::integrate() const Line 990  Data::integrate() const
990      }      }
991    }    }
992    if (rank==2) {    if (rank==2) {
993      bp_array.resize(shape[0],shape[1]);         bp_array.resize(shape[0],shape[1]);
994      for (int i=0; i<shape[0]; i++) {         for (int i=0; i<shape[0]; i++) {
995        for (int j=0; j<shape[1]; j++) {           for (int j=0; j<shape[1]; j++) {
996          index = i + shape[0] * j;             index = i + shape[0] * j;
997          bp_array[i,j] = integrals[index];             bp_array[make_tuple(i,j)] = integrals[index];
998        }           }
999      }         }
1000    }    }
1001    if (rank==3) {    if (rank==3) {
1002      bp_array.resize(shape[0],shape[1],shape[2]);      bp_array.resize(shape[0],shape[1],shape[2]);
# Line 595  Data::integrate() const Line 1004  Data::integrate() const
1004        for (int j=0; j<shape[1]; j++) {        for (int j=0; j<shape[1]; j++) {
1005          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
1006            index = i + shape[0] * ( j + shape[1] * k );            index = i + shape[0] * ( j + shape[1] * k );
1007            bp_array[i,j,k] = integrals[index];            bp_array[make_tuple(i,j,k)] = integrals[index];
1008          }          }
1009        }        }
1010      }      }
# Line 607  Data::integrate() const Line 1016  Data::integrate() const
1016          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
1017            for (int l=0; l<shape[3]; l++) {            for (int l=0; l<shape[3]; l++) {
1018              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
1019              bp_array[i,j,k,l] = integrals[index];              bp_array[make_tuple(i,j,k,l)] = integrals[index];
1020            }            }
1021          }          }
1022        }        }
# Line 622  Data::integrate() const Line 1031  Data::integrate() const
1031  Data  Data
1032  Data::sin() const  Data::sin() const
1033  {  {
1034    #if defined DOPROF
1035      profData->unary++;
1036    #endif
1037    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
1038  }  }
1039    
1040  Data  Data
1041  Data::cos() const  Data::cos() const
1042  {  {
1043    #if defined DOPROF
1044      profData->unary++;
1045    #endif
1046    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
1047  }  }
1048    
1049  Data  Data
1050  Data::tan() const  Data::tan() const
1051  {  {
1052    #if defined DOPROF
1053      profData->unary++;
1054    #endif
1055    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
1056  }  }
1057    
1058  Data  Data
1059  Data::log() const  Data::asin() const
1060    {
1061    #if defined DOPROF
1062      profData->unary++;
1063    #endif
1064      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
1065    }
1066    
1067    Data
1068    Data::acos() const
1069    {
1070    #if defined DOPROF
1071      profData->unary++;
1072    #endif
1073      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
1074    }
1075    
1076    Data
1077    Data::atan() const
1078    {
1079    #if defined DOPROF
1080      profData->unary++;
1081    #endif
1082      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
1083    }
1084    
1085    Data
1086    Data::sinh() const
1087    {
1088    #if defined DOPROF
1089      profData->unary++;
1090    #endif
1091      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
1092    }
1093    
1094    Data
1095    Data::cosh() const
1096    {
1097    #if defined DOPROF
1098      profData->unary++;
1099    #endif
1100      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
1101    }
1102    
1103    Data
1104    Data::tanh() const
1105    {
1106    #if defined DOPROF
1107      profData->unary++;
1108    #endif
1109      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
1110    }
1111    
1112    Data
1113    Data::asinh() const
1114    {
1115    #if defined DOPROF
1116      profData->unary++;
1117    #endif
1118      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
1119    }
1120    
1121    Data
1122    Data::acosh() const
1123    {
1124    #if defined DOPROF
1125      profData->unary++;
1126    #endif
1127      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
1128    }
1129    
1130    Data
1131    Data::atanh() const
1132  {  {
1133    #if defined DOPROF
1134      profData->unary++;
1135    #endif
1136      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
1137    }
1138    
1139    Data
1140    Data::log10() const
1141    {
1142    #if defined DOPROF
1143      profData->unary++;
1144    #endif
1145    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
1146  }  }
1147    
1148  Data  Data
1149  Data::ln() const  Data::log() const
1150  {  {
1151    #if defined DOPROF
1152      profData->unary++;
1153    #endif
1154    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
1155  }  }
1156    
1157  Data  Data
1158  Data::sign() const  Data::sign() const
1159  {  {
1160    #if defined DOPROF
1161      profData->unary++;
1162    #endif
1163    return escript::unaryOp(*this,escript::fsign);    return escript::unaryOp(*this,escript::fsign);
1164  }  }
1165    
1166  Data  Data
1167  Data::abs() const  Data::abs() const
1168  {  {
1169    #if defined DOPROF
1170      profData->unary++;
1171    #endif
1172    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
1173  }  }
1174    
1175  Data  Data
1176  Data::neg() const  Data::neg() const
1177  {  {
1178    #if defined DOPROF
1179      profData->unary++;
1180    #endif
1181    return escript::unaryOp(*this,negate<double>());    return escript::unaryOp(*this,negate<double>());
1182  }  }
1183    
1184  Data  Data
1185  Data::pos() const  Data::pos() const
1186  {  {
1187    return (*this);  #if defined DOPROF
1188      profData->unary++;
1189    #endif
1190      Data result;
1191      // perform a deep copy
1192      result.copy(*this);
1193      return result;
1194  }  }
1195    
1196  Data  Data
1197  Data::exp() const  Data::exp() const
1198  {  {
1199    #if defined DOPROF
1200      profData->unary++;
1201    #endif
1202    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1203  }  }
1204    
1205  Data  Data
1206  Data::sqrt() const  Data::sqrt() const
1207  {  {
1208    #if defined DOPROF
1209      profData->unary++;
1210    #endif
1211    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1212  }  }
1213    
1214  double  double
1215  Data::Lsup() const  Data::Lsup() const
1216  {  {
1217      double localValue, globalValue;
1218    #if defined DOPROF
1219      profData->reduction1++;
1220    #endif
1221    //    //
1222    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
1223    return algorithm(DataAlgorithmAdapter<AbsMax>(0));  
1224      AbsMax abs_max_func;
1225      localValue = algorithm(abs_max_func,0);
1226    #ifdef PASO_MPI
1227      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1228      return globalValue;
1229    #else
1230      return localValue;
1231    #endif
1232  }  }
1233    
1234  double  double
1235  Data::Linf() const  Data::Linf() const
1236  {  {
1237      double localValue, globalValue;
1238    #if defined DOPROF
1239      profData->reduction1++;
1240    #endif
1241    //    //
1242    // set the initial absolute minimum value to max double    // set the initial absolute minimum value to max double
1243    return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));    AbsMin abs_min_func;
1244      localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1245    
1246    #ifdef PASO_MPI
1247      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1248      return globalValue;
1249    #else
1250      return localValue;
1251    #endif
1252  }  }
1253    
1254  double  double
1255  Data::sup() const  Data::sup() const
1256  {  {
1257      double localValue, globalValue;
1258    #if defined DOPROF
1259      profData->reduction1++;
1260    #endif
1261    //    //
1262    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1263    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1264      localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1265    #ifdef PASO_MPI
1266      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1267      return globalValue;
1268    #else
1269      return localValue;
1270    #endif
1271  }  }
1272    
1273  double  double
1274  Data::inf() const  Data::inf() const
1275  {  {
1276      double localValue, globalValue;
1277    #if defined DOPROF
1278      profData->reduction1++;
1279    #endif
1280    //    //
1281    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1282    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1283      localValue = algorithm(fmin_func,numeric_limits<double>::max());
1284    #ifdef PASO_MPI
1285      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1286      return globalValue;
1287    #else
1288      return localValue;
1289    #endif
1290  }  }
1291    
1292    /* TODO */
1293    /* global reduction */
1294  Data  Data
1295  Data::maxval() const  Data::maxval() const
1296  {  {
1297    #if defined DOPROF
1298      profData->reduction2++;
1299    #endif
1300    //    //
1301    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1302    return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1303      return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1304  }  }
1305    
1306  Data  Data
1307  Data::minval() const  Data::minval() const
1308  {  {
1309    #if defined DOPROF
1310      profData->reduction2++;
1311    #endif
1312    //    //
1313    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1314    return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1315  }    return dp_algorithm(fmin_func,numeric_limits<double>::max());
   
 Data  
 Data::length() const  
 {  
   return dp_algorithm(DataAlgorithmAdapter<Length>(0));  
1316  }  }
1317    
1318  Data  Data
1319  Data::trace() const  Data::trace() const
1320  {  {
1321    return dp_algorithm(DataAlgorithmAdapter<Trace>(0));  #if defined DOPROF
1322      profData->reduction2++;
1323    #endif
1324      Trace trace_func;
1325      return dp_algorithm(trace_func,0);
1326    }
1327    
1328    Data
1329    Data::symmetric() const
1330    {
1331         #if defined DOPROF
1332            profData->unary++;
1333         #endif
1334         // check input
1335         DataArrayView::ShapeType s=getDataPointShape();
1336         if (getDataPointRank()==2) {
1337            if(s[0] != s[1])
1338               throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1339         }
1340         else if (getDataPointRank()==4) {
1341            if(!(s[0] == s[2] && s[1] == s[3]))
1342               throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1343         }
1344         else {
1345            throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1346         }
1347         Data ev(0.,getDataPointShape(),getFunctionSpace());
1348         ev.typeMatchRight(*this);
1349         m_data->symmetric(ev.m_data.get());
1350         return ev;
1351    }
1352    
1353    Data
1354    Data::nonsymmetric() const
1355    {
1356         #if defined DOPROF
1357            profData->unary++;
1358         #endif
1359         // check input
1360         DataArrayView::ShapeType s=getDataPointShape();
1361         if (getDataPointRank()==2) {
1362            if(s[0] != s[1])
1363               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1364            DataArrayView::ShapeType ev_shape;
1365            ev_shape.push_back(s[0]);
1366            ev_shape.push_back(s[1]);
1367            Data ev(0.,ev_shape,getFunctionSpace());
1368            ev.typeMatchRight(*this);
1369            m_data->nonsymmetric(ev.m_data.get());
1370            return ev;
1371         }
1372         else if (getDataPointRank()==4) {
1373            if(!(s[0] == s[2] && s[1] == s[3]))
1374               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1375            DataArrayView::ShapeType ev_shape;
1376            ev_shape.push_back(s[0]);
1377            ev_shape.push_back(s[1]);
1378            ev_shape.push_back(s[2]);
1379            ev_shape.push_back(s[3]);
1380            Data ev(0.,ev_shape,getFunctionSpace());
1381            ev.typeMatchRight(*this);
1382            m_data->nonsymmetric(ev.m_data.get());
1383            return ev;
1384         }
1385         else {
1386            throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1387         }
1388    }
1389    
1390    Data
1391    Data::matrixtrace(int axis_offset) const
1392    {
1393         #if defined DOPROF
1394            profData->unary++;
1395         #endif
1396         DataArrayView::ShapeType s=getDataPointShape();
1397         if (getDataPointRank()==2) {
1398            DataArrayView::ShapeType ev_shape;
1399            Data ev(0.,ev_shape,getFunctionSpace());
1400            ev.typeMatchRight(*this);
1401            m_data->matrixtrace(ev.m_data.get(), axis_offset);
1402            return ev;
1403         }
1404         if (getDataPointRank()==3) {
1405            DataArrayView::ShapeType ev_shape;
1406            if (axis_offset==0) {
1407              int s2=s[2];
1408              ev_shape.push_back(s2);
1409            }
1410            else if (axis_offset==1) {
1411              int s0=s[0];
1412              ev_shape.push_back(s0);
1413            }
1414            Data ev(0.,ev_shape,getFunctionSpace());
1415            ev.typeMatchRight(*this);
1416            m_data->matrixtrace(ev.m_data.get(), axis_offset);
1417            return ev;
1418         }
1419         if (getDataPointRank()==4) {
1420            DataArrayView::ShapeType ev_shape;
1421            if (axis_offset==0) {
1422              ev_shape.push_back(s[2]);
1423              ev_shape.push_back(s[3]);
1424            }
1425            else if (axis_offset==1) {
1426              ev_shape.push_back(s[0]);
1427              ev_shape.push_back(s[3]);
1428            }
1429        else if (axis_offset==2) {
1430          ev_shape.push_back(s[0]);
1431          ev_shape.push_back(s[1]);
1432        }
1433            Data ev(0.,ev_shape,getFunctionSpace());
1434            ev.typeMatchRight(*this);
1435        m_data->matrixtrace(ev.m_data.get(), axis_offset);
1436            return ev;
1437         }
1438         else {
1439            throw DataException("Error - Data::matrixtrace can only be calculated for rank 2, 3 or 4 object.");
1440         }
1441    }
1442    
1443    Data
1444    Data::transpose(int axis_offset) const
1445    {
1446    #if defined DOPROF
1447         profData->reduction2++;
1448    #endif
1449         DataArrayView::ShapeType s=getDataPointShape();
1450         DataArrayView::ShapeType ev_shape;
1451         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1452         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1453         int rank=getDataPointRank();
1454         if (axis_offset<0 || axis_offset>rank) {
1455            throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1456         }
1457         for (int i=0; i<rank; i++) {
1458           int index = (axis_offset+i)%rank;
1459           ev_shape.push_back(s[index]); // Append to new shape
1460         }
1461         Data ev(0.,ev_shape,getFunctionSpace());
1462         ev.typeMatchRight(*this);
1463         m_data->transpose(ev.m_data.get(), axis_offset);
1464         return ev;
1465    }
1466    
1467    Data
1468    Data::eigenvalues() const
1469    {
1470         #if defined DOPROF
1471            profData->unary++;
1472         #endif
1473         // check input
1474         DataArrayView::ShapeType s=getDataPointShape();
1475         if (getDataPointRank()!=2)
1476            throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1477         if(s[0] != s[1])
1478            throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1479         // create return
1480         DataArrayView::ShapeType ev_shape(1,s[0]);
1481         Data ev(0.,ev_shape,getFunctionSpace());
1482         ev.typeMatchRight(*this);
1483         m_data->eigenvalues(ev.m_data.get());
1484         return ev;
1485    }
1486    
1487    const boost::python::tuple
1488    Data::eigenvalues_and_eigenvectors(const double tol) const
1489    {
1490         #if defined DOPROF
1491            profData->unary++;
1492         #endif
1493         DataArrayView::ShapeType s=getDataPointShape();
1494         if (getDataPointRank()!=2)
1495            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1496         if(s[0] != s[1])
1497            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1498         // create return
1499         DataArrayView::ShapeType ev_shape(1,s[0]);
1500         Data ev(0.,ev_shape,getFunctionSpace());
1501         ev.typeMatchRight(*this);
1502         DataArrayView::ShapeType V_shape(2,s[0]);
1503         Data V(0.,V_shape,getFunctionSpace());
1504         V.typeMatchRight(*this);
1505         m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1506         return make_tuple(boost::python::object(ev),boost::python::object(V));
1507    }
1508    
1509    const boost::python::tuple
1510    Data::mindp() const
1511    {
1512      // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an
1513      // abort (for unknown reasons) if there are openmp directives with it in the
1514      // surrounding function
1515    
1516      int SampleNo;
1517      int DataPointNo;
1518        int ProcNo;
1519      calc_mindp(ProcNo,SampleNo,DataPointNo);
1520      return make_tuple(ProcNo,SampleNo,DataPointNo);
1521  }  }
1522    
1523  Data  void
1524  Data::transpose(int axis) const  Data::calc_mindp(   int& ProcNo,
1525                    int& SampleNo,
1526            int& DataPointNo) const
1527  {  {
1528    // not implemented    int i,j;
1529    throw DataException("Error - Data::transpose not implemented yet.");    int lowi=0,lowj=0;
1530    return Data();    double min=numeric_limits<double>::max();
1531    
1532      Data temp=minval();
1533    
1534      int numSamples=temp.getNumSamples();
1535      int numDPPSample=temp.getNumDataPointsPerSample();
1536    
1537      double next,local_min;
1538      int local_lowi,local_lowj;
1539    
1540      #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1541      {
1542        local_min=min;
1543        #pragma omp for private(i,j) schedule(static)
1544        for (i=0; i<numSamples; i++) {
1545          for (j=0; j<numDPPSample; j++) {
1546            next=temp.getDataPoint(i,j)();
1547            if (next<local_min) {
1548              local_min=next;
1549              local_lowi=i;
1550              local_lowj=j;
1551            }
1552          }
1553        }
1554        #pragma omp critical
1555        if (local_min<min) {
1556          min=local_min;
1557          lowi=local_lowi;
1558          lowj=local_lowj;
1559        }
1560      }
1561    
1562    #ifdef PASO_MPI
1563        // determine the processor on which the minimum occurs
1564        next = temp.getDataPoint(lowi,lowj)();
1565        int lowProc = 0;
1566        double *globalMins = new double[get_MPISize()+1];
1567        int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1568        
1569        if( get_MPIRank()==0 ){
1570            next = globalMins[lowProc];
1571            for( i=1; i<get_MPISize(); i++ )
1572                if( next>globalMins[i] ){
1573                    lowProc = i;
1574                    next = globalMins[i];
1575                }
1576        }
1577        MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1578    
1579        delete [] globalMins;
1580        ProcNo = lowProc;
1581    #else
1582        ProcNo = 0;
1583    #endif
1584      SampleNo = lowi;
1585      DataPointNo = lowj;
1586  }  }
1587    
1588  void  void
1589  Data::saveDX(std::string fileName) const  Data::saveDX(std::string fileName) const
1590  {  {
1591    getDomain().saveDX(fileName,*this);    boost::python::dict args;
1592      args["data"]=boost::python::object(this);
1593      getDomain().saveDX(fileName,args);
1594    return;    return;
1595  }  }
1596    
1597  void  void
1598  Data::saveVTK(std::string fileName) const  Data::saveVTK(std::string fileName) const
1599  {  {
1600    getDomain().saveVTK(fileName,*this);    boost::python::dict args;
1601      args["data"]=boost::python::object(this);
1602      getDomain().saveVTK(fileName,args);
1603    return;    return;
1604  }  }
1605    
1606  Data&  Data&
1607  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1608  {  {
1609      if (isProtected()) {
1610            throw DataException("Error - attempt to update protected Data object.");
1611      }
1612    #if defined DOPROF
1613      profData->binary++;
1614    #endif
1615    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1616    return (*this);    return (*this);
1617  }  }
# Line 777  Data::operator+=(const Data& right) Line 1619  Data::operator+=(const Data& right)
1619  Data&  Data&
1620  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1621  {  {
1622      if (isProtected()) {
1623            throw DataException("Error - attempt to update protected Data object.");
1624      }
1625    #if defined DOPROF
1626      profData->binary++;
1627    #endif
1628    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1629    return (*this);    return (*this);
1630  }  }
# Line 784  Data::operator+=(const boost::python::ob Line 1632  Data::operator+=(const boost::python::ob
1632  Data&  Data&
1633  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1634  {  {
1635      if (isProtected()) {
1636            throw DataException("Error - attempt to update protected Data object.");
1637      }
1638    #if defined DOPROF
1639      profData->binary++;
1640    #endif
1641    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1642    return (*this);    return (*this);
1643  }  }
# Line 791  Data::operator-=(const Data& right) Line 1645  Data::operator-=(const Data& right)
1645  Data&  Data&
1646  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1647  {  {
1648      if (isProtected()) {
1649            throw DataException("Error - attempt to update protected Data object.");
1650      }
1651    #if defined DOPROF
1652      profData->binary++;
1653    #endif
1654    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1655    return (*this);    return (*this);
1656  }  }
# Line 798  Data::operator-=(const boost::python::ob Line 1658  Data::operator-=(const boost::python::ob
1658  Data&  Data&
1659  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1660  {  {
1661      if (isProtected()) {
1662            throw DataException("Error - attempt to update protected Data object.");
1663      }
1664    #if defined DOPROF
1665      profData->binary++;
1666    #endif
1667    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1668    return (*this);    return (*this);
1669  }  }
# Line 805  Data::operator*=(const Data& right) Line 1671  Data::operator*=(const Data& right)
1671  Data&  Data&
1672  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1673  {  {
1674      if (isProtected()) {
1675            throw DataException("Error - attempt to update protected Data object.");
1676      }
1677    #if defined DOPROF
1678      profData->binary++;
1679    #endif
1680    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1681    return (*this);    return (*this);
1682  }  }
# Line 812  Data::operator*=(const boost::python::ob Line 1684  Data::operator*=(const boost::python::ob
1684  Data&  Data&
1685  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1686  {  {
1687      if (isProtected()) {
1688            throw DataException("Error - attempt to update protected Data object.");
1689      }
1690    #if defined DOPROF
1691      profData->binary++;
1692    #endif
1693    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1694    return (*this);    return (*this);
1695  }  }
# Line 819  Data::operator/=(const Data& right) Line 1697  Data::operator/=(const Data& right)
1697  Data&  Data&
1698  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1699  {  {
1700      if (isProtected()) {
1701            throw DataException("Error - attempt to update protected Data object.");
1702      }
1703    #if defined DOPROF
1704      profData->binary++;
1705    #endif
1706    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1707    return (*this);    return (*this);
1708  }  }
1709    
1710  Data  Data
1711    Data::rpowO(const boost::python::object& left) const
1712    {
1713      if (isProtected()) {
1714            throw DataException("Error - attempt to update protected Data object.");
1715      }
1716    #if defined DOPROF
1717      profData->binary++;
1718    #endif
1719      Data left_d(left,*this);
1720      return left_d.powD(*this);
1721    }
1722    
1723    Data
1724  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1725  {  {
1726    #if defined DOPROF
1727      profData->binary++;
1728    #endif
1729    Data result;    Data result;
1730    result.copy(*this);    result.copy(*this);
1731    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
# Line 835  Data::powO(const boost::python::object& Line 1735  Data::powO(const boost::python::object&
1735  Data  Data
1736  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1737  {  {
1738    #if defined DOPROF
1739      profData->binary++;
1740    #endif
1741    Data result;    Data result;
1742    result.copy(*this);    result.copy(*this);
1743    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1744    return result;    return result;
1745  }  }
1746    
1747    
1748  //  //
1749  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1750  Data  Data
1751  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1752  {  {
# Line 855  escript::operator+(const Data& left, con Line 1759  escript::operator+(const Data& left, con
1759  }  }
1760    
1761  //  //
1762  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1763  Data  Data
1764  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1765  {  {
# Line 868  escript::operator-(const Data& left, con Line 1772  escript::operator-(const Data& left, con
1772  }  }
1773    
1774  //  //
1775  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1776  Data  Data
1777  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
1778  {  {
# Line 881  escript::operator*(const Data& left, con Line 1785  escript::operator*(const Data& left, con
1785  }  }
1786    
1787  //  //
1788  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1789  Data  Data
1790  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
1791  {  {
# Line 894  escript::operator/(const Data& left, con Line 1798  escript::operator/(const Data& left, con
1798  }  }
1799    
1800  //  //
1801  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1802  Data  Data
1803  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1804  {  {
# Line 910  escript::operator+(const Data& left, con Line 1814  escript::operator+(const Data& left, con
1814  }  }
1815    
1816  //  //
1817  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1818  Data  Data
1819  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1820  {  {
# Line 926  escript::operator-(const Data& left, con Line 1830  escript::operator-(const Data& left, con
1830  }  }
1831    
1832  //  //
1833  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1834  Data  Data
1835  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1836  {  {
# Line 942  escript::operator*(const Data& left, con Line 1846  escript::operator*(const Data& left, con
1846  }  }
1847    
1848  //  //
1849  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1850  Data  Data
1851  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1852  {  {
# Line 958  escript::operator/(const Data& left, con Line 1862  escript::operator/(const Data& left, con
1862  }  }
1863    
1864  //  //
1865  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1866  Data  Data
1867  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1868  {  {
# Line 971  escript::operator+(const boost::python:: Line 1875  escript::operator+(const boost::python::
1875  }  }
1876    
1877  //  //
1878  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1879  Data  Data
1880  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1881  {  {
# Line 984  escript::operator-(const boost::python:: Line 1888  escript::operator-(const boost::python::
1888  }  }
1889    
1890  //  //
1891  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1892  Data  Data
1893  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1894  {  {
# Line 997  escript::operator*(const boost::python:: Line 1901  escript::operator*(const boost::python::
1901  }  }
1902    
1903  //  //
1904  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1905  Data  Data
1906  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1907  {  {
# Line 1010  escript::operator/(const boost::python:: Line 1914  escript::operator/(const boost::python::
1914  }  }
1915    
1916  //  //
 // NOTE: It is essential to specify the namepsace this operator belongs to  
1917  //bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1918  //{  //{
1919  //  /*  //  /*
# Line 1055  escript::operator/(const boost::python:: Line 1958  escript::operator/(const boost::python::
1958  //  return ret;  //  return ret;
1959  //}  //}
1960    
1961    /* TODO */
1962    /* global reduction */
1963  Data  Data
1964  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
1965  {  {
# Line 1069  Data::getItem(const boost::python::objec Line 1974  Data::getItem(const boost::python::objec
1974    return getSlice(slice_region);    return getSlice(slice_region);
1975  }  }
1976    
1977    /* TODO */
1978    /* global reduction */
1979  Data  Data
1980  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1981  {  {
1982    #if defined DOPROF
1983      profData->slicing++;
1984    #endif
1985    return Data(*this,region);    return Data(*this,region);
1986  }  }
1987    
1988    /* TODO */
1989    /* global reduction */
1990  void  void
1991  Data::setItemO(const boost::python::object& key,  Data::setItemO(const boost::python::object& key,
1992                 const boost::python::object& value)                 const boost::python::object& value)
# Line 1083  Data::setItemO(const boost::python::obje Line 1995  Data::setItemO(const boost::python::obje
1995    setItemD(key,tempData);    setItemD(key,tempData);
1996  }  }
1997    
1998    /* TODO */
1999    /* global reduction */
2000  void  void
2001  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
2002                 const Data& value)                 const Data& value)
2003  {  {
2004    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
2005    
2006    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
2007    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
2008      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
# Line 1099  Data::setItemD(const boost::python::obje Line 2014  Data::setItemD(const boost::python::obje
2014    }    }
2015  }  }
2016    
2017    /* TODO */
2018    /* global reduction */
2019  void  void
2020  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
2021                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
2022  {  {
2023      if (isProtected()) {
2024            throw DataException("Error - attempt to update protected Data object.");
2025      }
2026    #if defined DOPROF
2027      profData->slicing++;
2028    #endif
2029    Data tempValue(value);    Data tempValue(value);
2030    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
2031    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1137  Data::typeMatchRight(const Data& right) Line 2060  Data::typeMatchRight(const Data& right)
2060    }    }
2061  }  }
2062    
2063    /* TODO */
2064    /* global reduction */
2065  void  void
2066  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2067                       const boost::python::object& value)                       const boost::python::object& value)
2068  {  {
2069      if (isProtected()) {
2070            throw DataException("Error - attempt to update protected Data object.");
2071      }
2072    //    //
2073    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2074    tag();    tag();
# Line 1158  Data::setTaggedValue(int tagKey, Line 2086  Data::setTaggedValue(int tagKey,
2086    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
2087  }  }
2088    
2089    /* TODO */
2090    /* global reduction */
2091    void
2092    Data::setTaggedValueFromCPP(int tagKey,
2093                                const DataArrayView& value)
2094    {
2095      if (isProtected()) {
2096            throw DataException("Error - attempt to update protected Data object.");
2097      }
2098      //
2099      // Ensure underlying data object is of type DataTagged
2100      tag();
2101    
2102      if (!isTagged()) {
2103        throw DataException("Error - DataTagged conversion failed!!");
2104      }
2105                                                                                                                  
2106      //
2107      // Call DataAbstract::setTaggedValue
2108      m_data->setTaggedValue(tagKey,value);
2109    }
2110    
2111    /* TODO */
2112    /* global reduction */
2113    int
2114    Data::getTagNumber(int dpno)
2115    {
2116      return m_data->getTagNumber(dpno);
2117    }
2118    
2119    /* TODO */
2120    /* global reduction */
2121  void  void
2122  Data::setRefValue(int ref,  Data::setRefValue(int ref,
2123                    const boost::python::numeric::array& value)                    const boost::python::numeric::array& value)
2124  {  {
2125      if (isProtected()) {
2126            throw DataException("Error - attempt to update protected Data object.");
2127      }
2128    //    //
2129    // Construct DataArray from boost::python::object input value    // Construct DataArray from boost::python::object input value
2130    DataArray valueDataArray(value);    DataArray valueDataArray(value);
# Line 1171  Data::setRefValue(int ref, Line 2134  Data::setRefValue(int ref,
2134    m_data->setRefValue(ref,valueDataArray);    m_data->setRefValue(ref,valueDataArray);
2135  }  }
2136    
2137    /* TODO */
2138    /* global reduction */
2139  void  void
2140  Data::getRefValue(int ref,  Data::getRefValue(int ref,
2141                    boost::python::numeric::array& value)                    boost::python::numeric::array& value)
# Line 1197  Data::getRefValue(int ref, Line 2162  Data::getRefValue(int ref,
2162    DataArrayView valueView = valueDataArray.getView();    DataArrayView valueView = valueDataArray.getView();
2163    
2164    if (rank==0) {    if (rank==0) {
2165      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");        boost::python::numeric::array temp_numArray(valueView());
2166          value = temp_numArray;
2167    }    }
2168    if (rank==1) {    if (rank==1) {
2169      for (int i=0; i < shape[0]; i++) {      for (int i=0; i < shape[0]; i++) {
# Line 1205  Data::getRefValue(int ref, Line 2171  Data::getRefValue(int ref,
2171      }      }
2172    }    }
2173    if (rank==2) {    if (rank==2) {
2174      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
2175          for (int j=0; j < shape[1]; j++) {
2176            value[i][j] = valueView(i,j);
2177          }
2178        }
2179    }    }
2180    if (rank==3) {    if (rank==3) {
2181      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
2182          for (int j=0; j < shape[1]; j++) {
2183            for (int k=0; k < shape[2]; k++) {
2184              value[i][j][k] = valueView(i,j,k);
2185            }
2186          }
2187        }
2188    }    }
2189    if (rank==4) {    if (rank==4) {
2190      throw DataException("Data::getRefValue error: only rank 1 data handled for now.");      for (int i=0; i < shape[0]; i++) {
2191          for (int j=0; j < shape[1]; j++) {
2192            for (int k=0; k < shape[2]; k++) {
2193              for (int l=0; l < shape[3]; l++) {
2194                value[i][j][k][l] = valueView(i,j,k,l);
2195              }
2196            }
2197          }
2198        }
2199    }    }
2200    
2201  }  }
2202    
 /*  
 Note: this version removed for now. Not needed, and breaks escript.cpp  
 void  
 Data::setTaggedValue(int tagKey,  
                      const DataArrayView& value)  
 {  
   //  
   // Ensure underlying data object is of type DataTagged  
   tag();  
   
   if (!isTagged()) {  
     throw DataException("Error - DataTagged conversion failed!!");  
   }  
                                                                                                                 
   //  
   // Call DataAbstract::setTaggedValue  
   m_data->setTaggedValue(tagKey,value);  
 }  
 */  
   
2203  void  void
2204  Data::archiveData(const std::string fileName)  Data::archiveData(const std::string fileName)
2205  {  {
# Line 1261  Data::archiveData(const std::string file Line 2225  Data::archiveData(const std::string file
2225      dataType = 3;      dataType = 3;
2226      cout << "\tdataType: DataExpanded" << endl;      cout << "\tdataType: DataExpanded" << endl;
2227    }    }
2228    
2229    if (dataType == -1) {    if (dataType == -1) {
2230      throw DataException("archiveData Error: undefined dataType");      throw DataException("archiveData Error: undefined dataType");
2231    }    }
# Line 1274  Data::archiveData(const std::string file Line 2239  Data::archiveData(const std::string file
2239    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
2240    int dataLength = getLength();    int dataLength = getLength();
2241    DataArrayView::ShapeType dataPointShape = getDataPointShape();    DataArrayView::ShapeType dataPointShape = getDataPointShape();
2242    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2243    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2244      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2245    }    }
2246    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2247    if (isTagged()) {    if (isTagged()) {
2248      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2249        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
# Line 1300  Data::archiveData(const std::string file Line 2265  Data::archiveData(const std::string file
2265    cout << ">" << endl;    cout << ">" << endl;
2266    
2267    //    //
2268    // Write common data items to archive file    // Open archive file
2269    ofstream archiveFile;    ofstream archiveFile;
2270    archiveFile.open(fileName.data(), ios::out);    archiveFile.open(fileName.data(), ios::out);
2271    
# Line 1308  Data::archiveData(const std::string file Line 2273  Data::archiveData(const std::string file
2273      throw DataException("archiveData Error: problem opening archive file");      throw DataException("archiveData Error: problem opening archive file");
2274    }    }
2275    
2276      //
2277      // Write common data items to archive file
2278    archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));    archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
2279    archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));    archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
2280    archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));    archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
# Line 1331  Data::archiveData(const std::string file Line 2298  Data::archiveData(const std::string file
2298      throw DataException("archiveData Error: problem writing to archive file");      throw DataException("archiveData Error: problem writing to archive file");
2299    }    }
2300    
   archiveFile.close();  
   
   if (!archiveFile.good()) {  
     throw DataException("archiveData Error: problem closing archive file");  
   }  
   
2301    //    //
2302    // Collect and archive underlying data values for each Data type    // Archive underlying data values for each Data type
2303      int noValues;
2304    switch (dataType) {    switch (dataType) {
2305      case 0:      case 0:
2306        // DataEmpty        // DataEmpty
2307          noValues = 0;
2308          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2309          cout << "\tnoValues: " << noValues << endl;
2310        break;        break;
2311      case 1:      case 1:
2312        // DataConstant        // DataConstant
2313          noValues = m_data->getLength();
2314          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2315          cout << "\tnoValues: " << noValues << endl;
2316          if (m_data->archiveData(archiveFile,noValues)) {
2317            throw DataException("archiveData Error: problem writing data to archive file");
2318          }
2319        break;        break;
2320      case 2:      case 2:
2321        // DataTagged        // DataTagged
2322          noValues = m_data->getLength();
2323          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2324          cout << "\tnoValues: " << noValues << endl;
2325          if (m_data->archiveData(archiveFile,noValues)) {
2326            throw DataException("archiveData Error: problem writing data to archive file");
2327          }
2328        break;        break;
2329      case 3:      case 3:
2330        // DataExpanded        // DataExpanded
2331          noValues = m_data->getLength();
2332          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2333          cout << "\tnoValues: " << noValues << endl;
2334          if (m_data->archiveData(archiveFile,noValues)) {
2335            throw DataException("archiveData Error: problem writing data to archive file");
2336          }
2337        break;        break;
2338    }    }
2339    
2340      if (!archiveFile.good()) {
2341        throw DataException("archiveData Error: problem writing data to archive file");
2342      }
2343    
2344      //
2345      // Close archive file
2346      archiveFile.close();
2347    
2348      if (!archiveFile.good()) {
2349        throw DataException("archiveData Error: problem closing archive file");
2350      }
2351    
2352  }  }
2353    
2354  void  void
# Line 1379  Data::extractData(const std::string file Line 2374  Data::extractData(const std::string file
2374    int flatShape[4];    int flatShape[4];
2375    
2376    //    //
2377    // Open the archive file and read common data items    // Open the archive file
2378    ifstream archiveFile;    ifstream archiveFile;
2379    archiveFile.open(fileName.data(), ios::in);    archiveFile.open(fileName.data(), ios::in);
2380    
# Line 1387  Data::extractData(const std::string file Line 2382  Data::extractData(const std::string file
2382      throw DataException("extractData Error: problem opening archive file");      throw DataException("extractData Error: problem opening archive file");
2383    }    }
2384    
2385      //
2386      // Read common data items from archive file
2387    archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));    archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2388    archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));    archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2389    archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));    archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
# Line 1400  Data::extractData(const std::string file Line 2397  Data::extractData(const std::string file
2397        dataPointShape.push_back(flatShape[dim]);        dataPointShape.push_back(flatShape[dim]);
2398      }      }
2399    }    }
2400    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2401    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2402      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2403    }    }
2404    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2405    if (dataType==2) {    if (dataType==2) {
2406      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2407        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
# Line 1415  Data::extractData(const std::string file Line 2412  Data::extractData(const std::string file
2412      throw DataException("extractData Error: problem reading from archive file");      throw DataException("extractData Error: problem reading from archive file");
2413    }    }
2414    
2415    archiveFile.close();    //
2416      // Verify the values just read from the archive file
   if (!archiveFile.good()) {  
     throw DataException("extractData Error: problem closing archive file");  
   }  
   
2417    switch (dataType) {    switch (dataType) {
2418      case 0:      case 0:
2419        cout << "\tdataType: DataEmpty" << endl;        cout << "\tdataType: DataEmpty" << endl;
# Line 1475  Data::extractData(const std::string file Line 2468  Data::extractData(const std::string file
2468    
2469    //    //
2470    // Load this DataVector with the appropriate values    // Load this DataVector with the appropriate values
2471      int noValues;
2472      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2473      cout << "\tnoValues: " << noValues << endl;
2474    switch (dataType) {    switch (dataType) {
2475      case 0:      case 0:
2476        // DataEmpty        // DataEmpty
2477          if (noValues != 0) {
2478            throw DataException("extractData Error: problem reading data from archive file");
2479          }
2480        break;        break;
2481      case 1:      case 1:
2482        // DataConstant        // DataConstant
2483          if (dataVec.extractData(archiveFile,noValues)) {
2484            throw DataException("extractData Error: problem reading data from archive file");
2485          }
2486        break;        break;
2487      case 2:      case 2:
2488        // DataTagged        // DataTagged
2489          if (dataVec.extractData(archiveFile,noValues)) {
2490            throw DataException("extractData Error: problem reading data from archive file");
2491          }
2492        break;        break;
2493      case 3:      case 3:
2494        // DataExpanded        // DataExpanded
2495          if (dataVec.extractData(archiveFile,noValues)) {
2496            throw DataException("extractData Error: problem reading data from archive file");
2497          }
2498        break;        break;
2499    }    }
2500    
2501      if (!archiveFile.good()) {
2502        throw DataException("extractData Error: problem reading from archive file");
2503      }
2504    
2505      //
2506      // Close archive file
2507      archiveFile.close();
2508    
2509      if (!archiveFile.good()) {
2510        throw DataException("extractData Error: problem closing archive file");
2511      }
2512    
2513    //    //
2514    // Construct an appropriate Data object    // Construct an appropriate Data object
2515    DataAbstract* tempData;    DataAbstract* tempData;
# Line 1513  Data::extractData(const std::string file Line 2533  Data::extractData(const std::string file
2533    }    }
2534    shared_ptr<DataAbstract> temp_data(tempData);    shared_ptr<DataAbstract> temp_data(tempData);
2535    m_data=temp_data;    m_data=temp_data;
   
2536  }  }
2537    
2538  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
# Line 1521  ostream& escript::operator<<(ostream& o, Line 2540  ostream& escript::operator<<(ostream& o,
2540    o << data.toString();    o << data.toString();
2541    return o;    return o;
2542  }  }
2543    
2544    /* Member functions specific to the MPI implementation */
2545    
2546    void
2547    Data::print()
2548    {
2549      int i,j;
2550      
2551      printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2552      for( i=0; i<getNumSamples(); i++ )
2553      {
2554        printf( "[%6d]", i );
2555        for( j=0; j<getNumDataPointsPerSample(); j++ )
2556          printf( "\t%10.7g", (getSampleData(i))[j] );
2557        printf( "\n" );
2558      }
2559    }
2560    
2561    int
2562    Data::get_MPISize() const
2563    {
2564        int error, size;
2565    #ifdef PASO_MPI
2566        error = MPI_Comm_size( get_MPIComm(), &size );
2567    #else
2568        size = 1;
2569    #endif
2570        return size;
2571    }
2572    
2573    int
2574    Data::get_MPIRank() const
2575    {
2576        int error, rank;
2577    #ifdef PASO_MPI
2578        error = MPI_Comm_rank( get_MPIComm(), &rank );
2579    #else
2580        rank = 0;
2581    #endif
2582        return rank;
2583    }
2584    
2585    MPI_Comm
2586    Data::get_MPIComm() const
2587    {
2588    #ifdef PASO_MPI
2589        return MPI_COMM_WORLD;
2590    #else
2591        return -1;
2592    #endif
2593    }
2594    

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

  ViewVC Help
Powered by ViewVC 1.1.26