/[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 110 by jgs, Mon Feb 14 04:14:42 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    
27  #include <iostream>  #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
581    Data::convertToNumArray()
582    {
583      //
584      // determine the total number of data points
585      int numSamples = getNumSamples();
586      int numDataPointsPerSample = getNumDataPointsPerSample();
587      int numDataPoints = numSamples * numDataPointsPerSample;
588    
589      //
590      // determine the rank and shape of each data point
591      int dataPointRank = getDataPointRank();
592      DataArrayView::ShapeType dataPointShape = getDataPointShape();
593    
594      //
595      // create the numeric array to be returned
596      boost::python::numeric::array numArray(0.0);
597    
598      //
599      // the rank of the returned numeric array will be the rank of
600      // the data points, plus one. Where the rank of the array is n,
601      // the last n-1 dimensions will be equal to the shape of the
602      // data points, whilst the first dimension will be equal to the
603      // total number of data points. Thus the array will consist of
604      // a serial vector of the data points.
605      int arrayRank = dataPointRank + 1;
606      DataArrayView::ShapeType arrayShape;
607      arrayShape.push_back(numDataPoints);
608      for (int d=0; d<dataPointRank; d++) {
609         arrayShape.push_back(dataPointShape[d]);
610      }
611    
612      //
613      // resize the numeric array to the shape just calculated
614      if (arrayRank==1) {
615        numArray.resize(arrayShape[0]);
616      }
617      if (arrayRank==2) {
618        numArray.resize(arrayShape[0],arrayShape[1]);
619      }
620      if (arrayRank==3) {
621        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
622      }
623      if (arrayRank==4) {
624        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
625      }
626      if (arrayRank==5) {
627        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
628      }
629    
630      //
631      // loop through each data point in turn, loading the values for that data point
632      // into the numeric array.
633      int dataPoint = 0;
634      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
635        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
636          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
637          if (dataPointRank==0) {
638            numArray[dataPoint]=dataPointView();
639          }
640          if (dataPointRank==1) {
641            for (int i=0; i<dataPointShape[0]; i++) {
642              numArray[dataPoint][i]=dataPointView(i);
643            }
644          }
645          if (dataPointRank==2) {
646            for (int i=0; i<dataPointShape[0]; i++) {
647              for (int j=0; j<dataPointShape[1]; j++) {
648                numArray[dataPoint][i][j] = dataPointView(i,j);
649              }
650            }
651          }
652          if (dataPointRank==3) {
653            for (int i=0; i<dataPointShape[0]; i++) {
654              for (int j=0; j<dataPointShape[1]; j++) {
655                for (int k=0; k<dataPointShape[2]; k++) {
656                  numArray[dataPoint][i][j][k]=dataPointView(i,j,k);
657                }
658              }
659            }
660          }
661          if (dataPointRank==4) {
662            for (int i=0; i<dataPointShape[0]; i++) {
663              for (int j=0; j<dataPointShape[1]; j++) {
664                for (int k=0; k<dataPointShape[2]; k++) {
665                  for (int l=0; l<dataPointShape[3]; l++) {
666                    numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);
667                  }
668                }
669              }
670            }
671          }
672          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    
778      //
779      // return the loaded array
780      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 450  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 472  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 486  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 498  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 513  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
1235    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
1243      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    return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));  #if defined DOPROF
1298      profData->reduction2++;
1299    #endif
1300      //
1301      // set the initial maximum value to min possible double
1302      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    return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));  #if defined DOPROF
1310  }    profData->reduction2++;
1311    #endif
1312  Data    //
1313  Data::length() const    // set the initial minimum value to max possible double
1314  {    FMin fmin_func;
1315    return dp_algorithm(DataAlgorithmAdapter<Length>(0));    return dp_algorithm(fmin_func,numeric_limits<double>::max());
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 656  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 663  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 670  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 677  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 684  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 691  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 698  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 714  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 734  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 747  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 760  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 773  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 789  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 805  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 821  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 837  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 850  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 863  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 876  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 889  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 934  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 948  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 962  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 978  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 1016  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 1037  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 1050  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 1076  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 1084  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  
2203  void  void
2204  Data::setTaggedValue(int tagKey,  Data::archiveData(const std::string fileName)
                      const DataArrayView& value)  
2205  {  {
2206      cout << "Archiving Data object to: " << fileName << endl;
2207    
2208    //    //
2209    // Ensure underlying data object is of type DataTagged    // Determine type of this Data object
2210    tag();    int dataType = -1;
2211    
2212    if (!isTagged()) {    if (isEmpty()) {
2213      throw DataException("Error - DataTagged conversion failed!!");      dataType = 0;
2214        cout << "\tdataType: DataEmpty" << endl;
2215    }    }
2216                                                                                                                    if (isConstant()) {
2217        dataType = 1;
2218        cout << "\tdataType: DataConstant" << endl;
2219      }
2220      if (isTagged()) {
2221        dataType = 2;
2222        cout << "\tdataType: DataTagged" << endl;
2223      }
2224      if (isExpanded()) {
2225        dataType = 3;
2226        cout << "\tdataType: DataExpanded" << endl;
2227      }
2228    
2229      if (dataType == -1) {
2230        throw DataException("archiveData Error: undefined dataType");
2231      }
2232    
2233    //    //
2234    // Call DataAbstract::setTaggedValue    // Collect data items common to all Data types
2235    m_data->setTaggedValue(tagKey,value);    int noSamples = getNumSamples();
2236      int noDPPSample = getNumDataPointsPerSample();
2237      int functionSpaceType = getFunctionSpace().getTypeCode();
2238      int dataPointRank = getDataPointRank();
2239      int dataPointSize = getDataPointSize();
2240      int dataLength = getLength();
2241      DataArrayView::ShapeType dataPointShape = getDataPointShape();
2242      vector<int> referenceNumbers(noSamples);
2243      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2244        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2245      }
2246      vector<int> tagNumbers(noSamples);
2247      if (isTagged()) {
2248        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2249          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
2250        }
2251      }
2252    
2253      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2254      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2255      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2256    
2257      //
2258      // Flatten Shape to an array of integers suitable for writing to file
2259      int flatShape[4] = {0,0,0,0};
2260      cout << "\tshape: < ";
2261      for (int dim=0; dim<dataPointRank; dim++) {
2262        flatShape[dim] = dataPointShape[dim];
2263        cout << dataPointShape[dim] << " ";
2264      }
2265      cout << ">" << endl;
2266    
2267      //
2268      // Open archive file
2269      ofstream archiveFile;
2270      archiveFile.open(fileName.data(), ios::out);
2271    
2272      if (!archiveFile.good()) {
2273        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));
2279      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
2280      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2281      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2282      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2283      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2284      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
2285      for (int dim = 0; dim < 4; dim++) {
2286        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2287      }
2288      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2289        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2290      }
2291      if (isTagged()) {
2292        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2293          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2294        }
2295      }
2296    
2297      if (!archiveFile.good()) {
2298        throw DataException("archiveData Error: problem writing to archive file");
2299      }
2300    
2301      //
2302      // Archive underlying data values for each Data type
2303      int noValues;
2304      switch (dataType) {
2305        case 0:
2306          // DataEmpty
2307          noValues = 0;
2308          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2309          cout << "\tnoValues: " << noValues << endl;
2310          break;
2311        case 1:
2312          // 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;
2320        case 2:
2321          // 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;
2329        case 3:
2330          // 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;
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
2355    Data::extractData(const std::string fileName,
2356                      const FunctionSpace& fspace)
2357    {
2358      //
2359      // Can only extract Data to an object which is initially DataEmpty
2360      if (!isEmpty()) {
2361        throw DataException("extractData Error: can only extract to DataEmpty object");
2362      }
2363    
2364      cout << "Extracting Data object from: " << fileName << endl;
2365    
2366      int dataType;
2367      int noSamples;
2368      int noDPPSample;
2369      int functionSpaceType;
2370      int dataPointRank;
2371      int dataPointSize;
2372      int dataLength;
2373      DataArrayView::ShapeType dataPointShape;
2374      int flatShape[4];
2375    
2376      //
2377      // Open the archive file
2378      ifstream archiveFile;
2379      archiveFile.open(fileName.data(), ios::in);
2380    
2381      if (!archiveFile.good()) {
2382        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));
2388      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2389      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2390      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2391      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2392      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2393      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2394      for (int dim = 0; dim < 4; dim++) {
2395        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2396        if (flatShape[dim]>0) {
2397          dataPointShape.push_back(flatShape[dim]);
2398        }
2399      }
2400      vector<int> referenceNumbers(noSamples);
2401      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2402        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2403      }
2404      vector<int> tagNumbers(noSamples);
2405      if (dataType==2) {
2406        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2407          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2408        }
2409      }
2410    
2411      if (!archiveFile.good()) {
2412        throw DataException("extractData Error: problem reading from archive file");
2413      }
2414    
2415      //
2416      // Verify the values just read from the archive file
2417      switch (dataType) {
2418        case 0:
2419          cout << "\tdataType: DataEmpty" << endl;
2420          break;
2421        case 1:
2422          cout << "\tdataType: DataConstant" << endl;
2423          break;
2424        case 2:
2425          cout << "\tdataType: DataTagged" << endl;
2426          break;
2427        case 3:
2428          cout << "\tdataType: DataExpanded" << endl;
2429          break;
2430        default:
2431          throw DataException("extractData Error: undefined dataType read from archive file");
2432          break;
2433      }
2434    
2435      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2436      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2437      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2438      cout << "\tshape: < ";
2439      for (int dim = 0; dim < dataPointRank; dim++) {
2440        cout << dataPointShape[dim] << " ";
2441      }
2442      cout << ">" << endl;
2443    
2444      //
2445      // Verify that supplied FunctionSpace object is compatible with this Data object.
2446      if ( (fspace.getTypeCode()!=functionSpaceType) ||
2447           (fspace.getNumSamples()!=noSamples) ||
2448           (fspace.getNumDPPSample()!=noDPPSample)
2449         ) {
2450        throw DataException("extractData Error: incompatible FunctionSpace");
2451      }
2452      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2453        if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2454          throw DataException("extractData Error: incompatible FunctionSpace");
2455        }
2456      }
2457      if (dataType==2) {
2458        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2459          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2460            throw DataException("extractData Error: incompatible FunctionSpace");
2461          }
2462        }
2463      }
2464    
2465      //
2466      // Construct a DataVector to hold underlying data values
2467      DataVector dataVec(dataLength);
2468    
2469      //
2470      // 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) {
2475        case 0:
2476          // DataEmpty
2477          if (noValues != 0) {
2478            throw DataException("extractData Error: problem reading data from archive file");
2479          }
2480          break;
2481        case 1:
2482          // DataConstant
2483          if (dataVec.extractData(archiveFile,noValues)) {
2484            throw DataException("extractData Error: problem reading data from archive file");
2485          }
2486          break;
2487        case 2:
2488          // DataTagged
2489          if (dataVec.extractData(archiveFile,noValues)) {
2490            throw DataException("extractData Error: problem reading data from archive file");
2491          }
2492          break;
2493        case 3:
2494          // DataExpanded
2495          if (dataVec.extractData(archiveFile,noValues)) {
2496            throw DataException("extractData Error: problem reading data from archive file");
2497          }
2498          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
2515      DataAbstract* tempData;
2516      switch (dataType) {
2517        case 0:
2518          // DataEmpty
2519          tempData=new DataEmpty();
2520          break;
2521        case 1:
2522          // DataConstant
2523          tempData=new DataConstant(fspace,dataPointShape,dataVec);
2524          break;
2525        case 2:
2526          // DataTagged
2527          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2528          break;
2529        case 3:
2530          // DataExpanded
2531          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2532          break;
2533      }
2534      shared_ptr<DataAbstract> temp_data(tempData);
2535      m_data=temp_data;
2536  }  }
 */  
2537    
2538  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
2539  {  {
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.110  
changed lines
  Added in v.790

  ViewVC Help
Powered by ViewVC 1.1.26